Irradiation of Astrophysical Objects -
SED and Flux Effects on Thermally Driven Winds
We develop a general method for the self consistent calculation of the hydrodynamics of an astrophysical object irradiated by a radiation field with an arbitrary strength and spectral energy distribution (SED). Using the XSTAR photoionization code, we calculate heating and cooling rates as a function of gas photoionization parameter and temperature for several examples of SEDs: bremsstrahlung, blackbody, hard and soft state XRBs, Type 1 and Type 2 AGN. As an application of our method we study the hydrodynamics of 1-dimensional spherical winds heated by a uniform radiation field using the code Athena++. We find that in all cases explored a wind settles into a transonic, steady state. The wind evolves along the radiative heating equilibrium curve until adiabatic cooling effects become important and the flow departs from radiative equilibrium. If the flow is heated very rapidly, for example as in a thermally unstable regime, the corresponding column density of gas is low. Perhaps one of the most intriguing results of our work is the two stage acceleration of the wind that happens when there are two thermally unstable regions and the flux is relatively high. The efficiency with which the radiation field transfers energy to the wind is dependent on the SED of the external source, particularly the relative flux of soft X-rays. These results suggest that detailed photoionization calculations are essential not only to predict spectra but also to properly capture the flow dynamics.
keywords:hydrodynamics - radiation: dynamics - methods: numerical - stars: winds, outflows - galaxies: active - X-rays: binaries
The gas dynamics of many astrophysical systems is affected by radiation from an external source. Planets and moons are irradiated by their host star (e.g., Johnstone, Hollenbach & Bally 1998; Alexander, Clarke & Pringle 2006; Owen et al. 2010) and white dwarves can be evaporated by their companion (e.g., Liu, Meyer & Meyer-Hofmeister 1995). Likewise outflows are radiatively driven by AGN (e.g., Begelman, McKee & Shields 1983 henceforth BMS83)and X-ray binaries (e.g., Basko & Sunyaev 1973; Basko et al. 1977; London, McCray & Auer 1981; Luketic et al. 2010) and the outer layers of accretion discs are irradiated by their inner regions. It is therefore important to develop very general yet robust numerical methods for simulating the heating of hydrodynamic systems by external radiation fields. We developed such a method using the photoionization code XSTAR (Bautista & Kallman 2001) and MHD code Athena++ (Gardiner & Stone 2005, 2008).
One possible mechanism for launching winds observed in AGNs and YSOs is thermal driving (see Lamers & Cassinelli 1999 and references therein). In a thermally driven wind, the thermal energy of the gas comes to dominate the gravitational binding energy of the central object and the wind launches. If the thermal energy at the base is insufficient to launch the wind this scenario can still operate provided an external radiation source, such as the central object, heats the gas to sufficiently high temperatures.
The innermost regions of AGN and XRBs are known to be strong emitters of X-rays. These can heat the gas and provide the thermal energy required for wind launching. An interesting question therefore is how the SED of the external source affects the outflowing wind? This has been previously explored in the context of disc photoevaporation models (e.g., Bally & Scoville 1982; Hollenbach, Johnstone & Shu 1993; Owen, Clarke & Ercolano 2012 and references therein). There are important observational implications because radiative heating can provide a link between the X-ray spectra emitted by the central region and the lines produced further out in the wind.
In the most basic physical picture, gas at the base of the wind is cold and slow and thermal processes dominate over the dynamical processes. The external radiation source heats the wind and the gas evolves along the radiative equilibrium curve. Once the gas is hot enough, dynamical processes such as adiabatic expansion become important and the gas departs from equilibrium. Precisely when this occurs depends on the flux of incident radiation. This is observationally interesting because the line emission and absorption depends on the temperature and ionization state of the wind.
Real systems are complex (e.g., they are non-spherical, could be optically thick and rotate). However, to illustrate our method and to isolate the effects of different SEDs on thermal winds we study a simple problem of 1-dimensional spherical winds in the optically thin limit. Observed SEDs are typically a superposition of more basic spectra, such as a low temperature blackbody and high frequency powerlaw. Therefore we first consider more basic SEDs such as bremsstrahlung and blackbody to elucidate the basic physics at play. We then use observationally motivated SEDs from Type 1 and Type 2 AGN and soft state and hard state XRBs and generate the corresponding heating rates using the XSTAR photoionization code as a function of gas temperature, T and ionization, . These heating rates are used to find spherically symmetric wind solutions using the MHD code Athena++. These are used to determine the effects of radiation flux and driving SED on the observationally interesting wind properties such as mass flux and absorption measure distribution.
The structure of this paper is as follows. In Section 2, we describe our numerical methods including the photoionization code XSTAR, the MHD code Athena++ and our implementation of the heating/cooling module. In Section 3, we apply these methods to find spherically symmetric winds thermally driven by different types of external radiation fields and explore their observational features. In Section 4, we summarize our findings and discuss the limitations and future applications of these methods, in particular their relevance to 2D and 3D problems.
2 Numerical Methods
We are interested in performing hydrodynamic simulations where the radiative heating/cooling is calculated self-consistently using a photoionization code. We implement a heating/cooling module into Athena++ where the energy equation is sourced by a net cooling rate that accounts for the radiative heating and cooling of the gas by an external source. The rate is calculated self-consistently using the XSTAR photoionization code as a function of the temperature and ionization parameter calculated by Athena++. Photoionization calculations are computationaly expensive so we precalculate on a grid . We then interpolate between points on this grid within Athena++ to calculate the net heating. Below we briefly describe our XSTAR simulation, Athena++ simulation and implementation of the heating/cooling term.
2.1 Photoionization - Xstar
To accurately model photoionization processes we used version 2.3 of the code XSTAR (Bautista & Kallman 2001). By assuming a certain SED of incident flux of radiation on a box of gas containing known species of gas, XSTAR calculates the heating/cooling rate as a function of temperature and ionization fraction in the gas. For a review of methods for modelling photoionized plasmas see Kallman (2010).
For ease in comparing results for different SEDs we relate the mean photon energy to the X-ray temperature via
where is the Boltzman constant and we use as the low energy X-ray cutoff. This ensures that at large ionization fractions, when Compton processes dominate, the equilibrium curves will match up for all SEDs. Though we find that the gas never reaches this Compton temperature of in transonic solutions, our goal was to have the greatest possible uniformity across our runs.
Within our XSTAR simulation we assume the heating/cooling in the gas is due to Compton, X-ray, Bremsstrahlung and line processes. Our main model dependent assumptions are 1) Composition of the gas. 2) Available atomic transitions. 3) Source SED.
We use the elemental abundances described in Verner, Barthel & Tytler (1994) where elements with abundances of hydrogen were set to zero. We found that gas composition can play an important role, namely in determining the equilibrium curve. The AGN cases exhibited an unphysical downwards dip at low ionization parameter. For these cases we therefore used the abundances in Lodders, Palme & Gail (2009). We expect most variations due to composition to occur at low ionization/temperature where line processes are important. Because we are mainly interested in the relationship between the equilibrium curve and the wind we chose to leave a detailed study of the gas composition on the equilibrium curve for future work.
Comprehensive atomic data is crucial for determining the thermal properties of gases hotter than . XSTARv2.35 incorporates over 200 000 lines for gas at a given . As photoionization codes improve by incorporating better atomic data and more comprehensive lists on atomic transitions our methods for calculating net heating will only improve.
Our main focus within this study is how the source SED affects the net heating rate and in turn affects the driven wind. We describe our philosophy regarding choice of SEDs in 2.3.2.
2.2 Hydrodynamics - Athena++
2.2.1 Basic Equations
The basic equations for single fluid hydrodynamics are
where is the fluid density, the velocity, a diagonal tensor with components P the gas pressure, is the gravitational potential of the central object and is the energy where is the internal energy. is the net cooling rate described in Section 2.3 and is assumed to be a function of temperature and ionization parameter
where is the X-ray flux which we assume to be uniform, is the mean molecular weight and is the proton mass. We take an equation of state where . The isothermal sound speed is and the adiabatic sound speed . We can compute the temperature from the internal energy via .
2.2.2 Simulation Parameters
We choose parameters for our different runs to be able to consistently compare simulations between different classes of SEDs. A useful measure for thermally driven winds is the hydrodynamic escape parameter , which measures the ratio of gravitational to thermal energy. At the base of the wind, for the gas is gravitationally bound and no thermally driven outflow is expected. However for thermally driven hydrodynamic winds will be produced (Stone & Proga 2009). We fix the central object mass to help in comparing results across our different runs. Our results could be applied to systems with different masses, for instance AGN systems, if the relevant length (see eq 4) and ionizing flux (see eq 5) are scaled appropriately. At the base of the wind, for most runs, we set the initial temperature and set . This sets the inner radius . The exception is for the Blackbody SED, because the equilibrium curve is thermally unstable at this temperature, we set which corresponds to a and is still in the regime where winds can be thermally driven. These temperatures ensure that we are on the flat part of the S curve in the cold phase for all SEDs and also slightly away from any thermally unstable region. By initializing the gas in equilibrium near the star, the ionisation parameter is fixed for each SED.
We explored the effects of varying and values along the base of the radiative equilibrium S curve. For the mBl case with we compared our fiducial case with and cases. We found that at large radii the velocity differed by and and temperature varied by and respectively. We concluded that the outflows at large radii are largely insensitive to position on the radiative equilibrium curve, as one would expect for transonic flows.
A physically relevant length scale for this problem is the Compton radius
which is the radius at which gas at has internal energy equal to its gravitational potential energy. Above we have set . For a stellar mass central object and our fiducial choice of Compton temperature this corresponds to a Compton radius .
Near the Compton temperature, radiative processes will be dominated by Compton heating and cooling. The Compton heating rate can be used to define a critical radiation flux
where . This corresponds to the flux at which X-rays can heat the gas to via Compton processes in the sound crossing time of . This is the analogue of the critical luminosity of BMS83 (eq 2.12) for a uniform radiation field. We explore simulations where the flux ranges from .
We are interested in wind solutions that have been launched due to heating. Therefore we set the radius and the gas is initially at which are favourable conditions for radiative wind launching. By varying the flux we can study the effects of the particular SED on the launched winds.
Assuming a uniform radiation field allows us to define a fiducial density at the base of the wind for each SED via
Because the ionization parameter is fixed, a change in flux is equivalent to a change in density of the gas. A summary of our simulation parameters for all runs are shown in Table 1. In Section 4, we discuss the effects of non-uniform radiation fields, caused by geometric or optical depth effects.
|Modified Blondin (mBl)||None||0.28||3.43||5.24||2.50||5||12.9|
|Soft State XRB (XRB1)||Isobaric||0.28||3.43||24.09||2.50||5||2.81|
|Hard State XRB (XRB2)||Isobaric||0.75||1.29||12.29||2.50||5||9.02|
|Type 1 AGN (AGN1)||Isobaric||1.00||0.96||2.16||2.50||5||59.2|
|Type 2 AGN (AGN2)||Isobaric||1.52||0.63||6.47||2.50||5||24.4|
2.2.3 Grid & Boundary Conditions
The simulation region extends from . By definition of the critical flux, we expect the wind to heat to roughly over these length scales. We use a logarithmically spaced grid of points and a scale factor that defines the grid spacing recursively via .
At the inner boundary we impose outflow boundary conditions on and while keeping the density fixed at in the first active zone. As a result the ionization parameter is fixed in the first active zone. Hence we do not enforce radiative equilibrium in the first active zone except at the initial time because the temperature can evolve away from its initial equilibrium value. At the outer boundary we impose outflow boundary conditions on , and .
Within Athena++ we assume that all the microphysics incorporated into XSTAR is captured by a net cooling rate that depends only on the macroscopic gas properties, temperature and ionization parameter . We have implicitly assumed that the microphysical processes occur on much shorter time scales than the dynamical time governing the hydro simulation which justifies this approach. We give a brief description of our Athena++ heating module in Section 2.3.1 and provide a detailed description along with some numerical testing in Appendix A. We describe our different choices of SED providing this heating in Section 2.3.2.
2.3.1 Athena++ Implementation
First we calculate on a logarithmically spaced grid of size over the range and using XSTAR. This takes roughly one week running on 12 cores for each SED.
At every timestep in the hydro simulation we must calculate the rate of change in energy of the gas from radiative processes, the term appearing in the energy equation (2c). This is related to the heating rate and cooling rate calculated by XSTAR via
where At every cell center location we calculate the ionization parameter using (3). We then bilinearly interpolate over the nearest grid points in from our XSTAR grid. A backwards Euler method was used to solve for the heating rate under the assumption that the change in internal energy is due solely to the net heating and not the dynamics of the gas. Using the backwards Euler method ensures that the gas reaches equilibrium smoothly. This is important in this problem because the base of the wind is in equilibrium and this prevents any transients from entering the solution. A detailed description of this method is provided in Appendix A.
2.3.2 Driving Radiation Field
Our goal is to use observed SEDs to generate net heating rates to generate realistic hydrodynamics models. We assume this radiation field is uniform throughout the simulation region, as might happen for instance if the central object is surrounded by an optically thin gas and illuminated by a large number of uniform point sources. We first use simpler heating and cooling models, derived from idealized SED’s that have equilibrium curves with progresively more complicated features.
We first consider the analytic heating model proposed by Blondin (1994) as a fit to thermal Bremsstrahlung where
where we have used an X-ray temperature and a line temperature . We have chosen the X-ray temperature so that . Above the terms represent the parametrized heating and cooling from Compton processes (), X-rays (), Bremsstrahlung () and lines (). They depend only on the temperature of the gas T and the ionization parameter .
In our Modified Blondin (mBl) model we set , . This corresponds to the parametrization of “Model C” in Higginbottom & Proga (2015) but with an increased X-ray temperature. These parameters were chosen so that the equilibrium curve is free of the thermal instability (TI) (Field 1965). In the Blondin (Bl) model we set which is the original fit proposed by Blondin (1994) but with increased X-ray temperature. The equilibrium curve is unstable to isobaric perturbations. These analytic models were also used to test our interpolative heating scheme (see Appendix A). The Bremsstrahlung (Brem) model uses an SED of the same name as input to XSTAR for generating net heating rates. The Blackbody (BB) model provides a case where the equilibrium curve is thermally unstable to isochoric perturbations (Kallman & McCray 1982; Buff & McCray 1974). It can be viewed as either a true blackbody source from optically thin emission or as a proxy for an SED with a low energy cutoff say from absorption by cold gas near the source.
The remaining models are based on SEDs obtained from observations. We consider an X-ray binary in a soft state (XRB1) and a hard state (XRB2) (Trigo et al. 2013) and Active Galactic Nuclei SEDs from Type 1 (AGN1) and Type 2 (AGN2) AGN (Mehdipour et al. 2015). We note that the shape of the equilibrium curve is highly dependent on the atomic abundances used, in particular at low values of the ionization parameter where the cooling is largely dominated by lines. As an illustration of our method we used AGN spectra for NGC 5548 at two different epochs. However, for other applications like modeling AGN feedback one might want to use more typical AGN spectra (e.g. Sazonov, Ostriker & Sunyaev 2004; Sazonov et al. 2005). We show plots of these SEDs and the corresponding equilibrium curves output by XSTAR in Fig. 1. We indicate the initial conditions at the base of the wind for each SED with a dot.
At the base of the wind the velocity is subsonic and gravitational energy dominates over thermal energy () so basic energy considerations require that energy be added for a transonic wind to launch. This energy is provided by an external radiation field which differs according to each SED.
Irrespective of SED, the flows follow qualitatively similar behaviour - near the base of the wind they are in radiative equilibrium and follow the equilibrium curve until dynamical processes become important and they begin to radiatively heat. We briefly describe the dynamics of the flow in Section 3.1 for completenesss. However, we are primarily interested in the bulk properties of the flow - the mass flux and efficiency of momentum and energy tranfer to the wind. In particular, we compare how the different radiation fields affect these bulk observable properties in Section 3.3. We then discuss the spectroscopic properties of these flows by exploring their signatures on the absorption measure distribution Section 3.4
3.1 Dynamical Variables
For each SED we find stationary wind solutions for a range of subcritical radiation fluxes. The density is monotonically decreasing and velocity almost always monotonically increasing yielding a constant mass flux at all radii as expected for a stationary flow. The exception is at thermally unstable points where there is a drop in outwards pressure, leading to a dip in the velocity profile and a dip in the mass flux. The flows are qualitativly similar whereby they follow the equilibrium curve near the base of the wind before dynamical processes become important and they begin to heat. We summarize this for all SED models in Fig. 2 where we show the phase space evolution of each solution. Each color corresponds to a fixed fraction of the critical radiation flux for each model = (purple), (blue), (green), (orange), (red), (pink) and (black). For larger incident flux, equilibrium is maintained for larger ionizations allowing the flow to reach higher temperatures. Our definition of critical flux (5) appears consistent with these results since simulations with critical flux (purple) asymptotically approaches the Compton temperature.
As the flow accelerates adiabatic cooling becomes an important contributor to the energy balance, namely the gas temperature is not determined by the balance between radiative heating and cooling only but rather radiative heating, radiative cooling and adiabatic cooling. Consequently, the gas temperature is below the radiatve equilibrium temperature. When flux is sufficiently high (e.g. see the purple curves corresponding to the critical flux) the wind can undergo a second phase of such acceleration (Fig. 3). When this occurs there is a large spike in the heating rate with most of the energy going into thermal energy. This thermal energy is used by the flow to first exit the gravitational potential well. The kinetic energy then dominates over the gravitational energy and thermal energy goes to further increase the velocity of the wind. At large distances the flow is well developed with total energy equipartitioned between thermal and kinetic energy. We note in particular that at the sonic point the rate of heating satisfies (see Fig. 4), as expected from general energetic arguments for externally heated transonic winds (Holzer & Axford 1970; Lamers & Cassinelli 1999).
The physics of radiatively heated winds can be understood by considering the energetics of the wind. Generally, there is a competiton between the thermal energy and gravitational energy per unit mass. If gravitational energy dominates we have an accreting Bondi solution, but if thermal energy dominates then we have a thermal wind solution (Parker 1958). Because we are interested in studying the effects of heating, we considered initial conditions where the ratio of gravitational to thermal energy, as measured by . Therefore, in order to avoid an accreting Bondi solution, must come to dominate the energy of the gas.
Suppose we have reached a stationary solution , then we may express (2) as the Bernoulli function (e.g Lamers & Cassinelli 1999)
We may interpret each term on the left hand side as telling us the change in energy per unit mass of the wind in the kinetic, thermal and gravitational energy as we move outwards in the wind. The radiative heating acts as a source of heat . We ensure that the energy budget is satisfied by plotting both sides of (9) which is possible because we have an analytic expression for the heating.
In the lefthand panel of Fig. 4 we plot , , and the total energy for the representative cases and for the mBl SED. In the right hand panel of Fig. 4, we plot , , and . Since we are examining steady state solutions this allows us to see where energy is injected into the flow.
At small radii the thermal energy follows the equilibrium curve. When the thermal energy approximately equals the magnitude of gravitational energy the total energy is approximately zero (indicated by the spike in ). Between this point and the sonic point energy is primarilly injected into the thermal component. This is where the thermal component falls off the radiative equilibrium curve, since the flow velocity is non-negligible and adiabatic cooling dominates over radiative cooling.
3.3 Bulk Properties
We are interested in how the bulk properties of the flow are affected by the radiation field. In particular, we are interested in the flux of mass, particles, momentum and kinetic energy at the outer boundary and the corresponding efficiencies
The efficiencies have been defined by normalizing each of the bulk wind fluxes by the corresponding flux of the radiation field. We interpret as an efficiency because it is a measure of the coupling strength between the radiation field and the gas flow for that particular flux.
We plot each of the above fluxes and efficiencies in Fig. 5. First we notice that for a given SED, kinematic fluxes scale with radiative flux . Averaging over SEDs we find , and . From the heating rate argument used to derive the critical flux we expect . If there is an equipartition of kinetic and thermal energy we would therefore expect . However, we see this dependence is in fact much weaker with Starting with the particle flux, each subsequent quantity carries an additional power of which from the basic scaling argument is expected to add an additional power of . Comparing different SEDs we notice that there are two classes of SED - the main group and the BB SED which is somewhat of an outlier. We note that the BB SED lacks soft photons relative to all the other models, which leads to an overall decrease in the dynamical fluxes. This is a general trend with the other SEDs whereas the AGN SEDs have dynamical fluxes roughly an order of magnitude larger for fixed total flux. As a figure of merit for AGN, for Brem and for BB. The BB fluxes resemble those of the other SEDs if we plot instead of total flux. This reflects the microscopic properties of the heating i.e. in the subsonic part of the wind, where dynamical fluxes are determined, it is the soft X-ray photons which play a dominant role in heating the wind. This can also be seen from the equilibrium curves where the base of the S-curve is determined by line cooling and heating due to photoionization, which are sensitive to photons from UV to soft X-rays.
The efficiencies also have a powerlaw scaling with , and . Unlike the dynamical fluxes, the efficiencies are relatively independent of the radiation flux. In particular, increasing the efficiency of mass, particle and momentum flux is easier to achieve by changing the external SED rather than increasing the radiative flux. The efficiency in driving the wind is more complex than the basic parameters characterizing the radiation field. For instance 5 models have the same X-ray temperature (mBl, Bl, Brem, BB and XRB1). Even ignoring (BB), the outlying SED, which is much weaker in soft photons, the remaining SEDs still show a spread of roughly half an order of magnitude in efficiency.
This suggests that accurately modeling the gross wind properties requires accurate modeling of the radiation field. In particular, simply knowing the X-ray temperature and total flux is insufficient because as we have shown, efficiencies are weakly dependent on flux. For a fixed X-ray temperture, the efficiency can vary by almost an order of magnitude. The dynamical fluxes at large distances are a consequence of those same fluxes at the sonic point. However, because the sonic point occurs far away from the Compton radius, the primary heating mechanism is not Compton heating but rather photoionization. This heating rate is sensitive to the soft X-ray part of the radiation field and therefore a careful modeling of the gas is necessary.
3.4 Absorption Measure Distribution
One may test the validity of theoretical wind solution by comparing wind densities, velocities and photoionization with observed values. However, in addition to finding that solutions have gas in the observed parameter space one should ensure that a sufficiently high column density of gas with those properties is present to in fact be observable. One such measure is the Absorption Measure Distribution (AMD; Holczer, Behar & Kaspi 2007) which measures the absorbing column of gas at a particular ionization state along a line of sight. It is defined as
where is the total hydrogen column density along the line of sight and is the characteristic length for particle number variations and we have used our uniform ionizing flux model (see eq. (3)). This is observationally interesting because particular lines will only be absorbed by a gas in a particular ionization state (see for example Adhikari et al. 2015). However, if the gas density or characteristic length or both are too low any such spectral features will be absent. We emphasize that important comparisons can be made to observations by detecting individual lines at a particular ionization state and that comparison with the entire distribution is unnecessary.
We find that all our models exhibit dips in the AMD where they experience rapid heating. This can occur in two ways - 1) when the equilibrium curve is formally unstable or 2) dynamical processes in the flow become important, triggering heating. The AMD behaves qualitatively differently in these different flows. In the case of a formally stable equilibrium curve (mBl), a dip occurs when dynamical processes become important and the flow falls out of equilibrium. At large radii the AMD follows a powerlaw, largely independent of any dips at smaller radii. The thermal instability leads to sharper dips because of the larger heating rates experienced by the gas. After the thermal instability, the gas can reach a new stable intermediate temperature state (see XRB1 for example). The AMD then reaches a level commensurate to where it was before the unstable zone. If no such phase exists (see Brem for example) then the AMD will remain small at these ionizations. This is because along the equilibrium curve the dynamical time scales are less important so higher densities of the gas can exist. The AMD can therefore probe several interesting features of the equilibrium curve. It can tell us about what parts of the gas are thermally unstable and also indicate that parts of the gas exist along an intermediate stable temperature. These predictions require both accurate modeling of the heating rates (to find the thermally unstable regions) and hydrodynamic flow (to see what parts of the equilibrium curve are probed).
Using the photionization code XSTAR we used both idealized and observed SEDs to generate the heating/cooling rates of a gas as a function of temperature and ionization parameter. We implemented this into a heating module in the hydrodynamics code Athena++ to model spherically symmetric, thermally driven winds. The general behaviour and evolution can be characterized as follows.
1) We find stationary solutions that follow the classic picture of thermally driven winds - energy is added to the base of the wind in the subsonic region to overcome the gravitational potential and to accelerate the wind to supersonic velocities. The mass loss rate is determined by the amount of energy deposited at and below the sonic point, whereas other properties such as the wind velocity and energy flux are determined by the total amount of energy deposited in the supersonic part of the flow.
2) The mass, particle, momentum and kinetic energy fluxes are functions of the radiation field flux and the nature of the radiation field. The efficiency of driving the wind is a weak function of radiation flux but strongly dependent on the details of the SED.
3) The AMD provides a useful tool for determining the behaviour of the gas. When the gas experiences rapid heating there is a steep drop in the AMD. If gas exists at a stable intermediate temperature the AMD will have a bump at that ionization state.
To lowest order, one can assume that heating is dominated by Compton processes and the heating rate is simply parametrized by the X-ray temperature. This is a valid approximation at high temperatures and correctly predicts that the flow heats to the Compton temperature for a critical flux of incident radiation. However, at low temperatures, , Compton processes are subdominant to photoionization heating and cooling due to line and free-free emission. Therefore the X-ray temperature is insufficient to predict the state of the gas. One may then go one step further and compute the equilibrium curve for a given SED. The state of the flow is approximated by traversing the equilibrium curve. Regions where the equilibrium temperature is very sensitive to the photoionization parameter indicate highly efficient net heating that can lead even to catastrophic heating with TI being an example. A high energy deposition rate results in a flow acceleration and that in turn leads to adiabatic cooling. When the flow traverses another such region, as is possible for sufficiently high illuminating flux, it experiences another phase of acceleration. Consequently, the state of the flow is not properly captured by the radiative equilibrium solution. Our simulations show that falling out of radiative equilibrium has several interesting effects on the flow and these have important observational consequences.
Let us first compare how models with the lowest level of complexity behave. The mBl, Bl, Brem and XRB1 models all have the same . They also begin with the same gravitational and thermal energy at the base of the wind. However they have nearly an order of magnitude difference in the various fluxes and corresponding efficiencies as shown in Fig. 5. For all cases, the sonic point occurs at a temperature , hence the relevant temperature scale for determining mass flux at the critical point is not the X-ray temperature. Models with the same X-ray temperature therefore will not necesarily have similar mass flux, as our simulations show. The other wind fluxes are proportional to the mass flux times additional powers of the velocity . Therefore they differ by roughly an order of mangitude as well between the various models.
What about a model where we simply traverse the equilibrium curve? This is equivalent to assuming we provide a large enough flux for the radiative processes to dominate over any hydrodynamics. As a proxy for this model we will use the critical luminosity cases. However, as was shown in Section 3.3 the various wind fluxes have power-law scaling and the efficiencies . Furthermore, it fails to predict when the flow falls out of equilibrium. This means that certain intermediate stable temperature branches which would naively be expected to exist from the equilibrium curve might never be visited because insufficient flux is heating the flow. This is observationally one of the most interesting features of this type of system because it provides a link between the macroscopic gas properties (temperature and ionization parameter) and the underlying microphysics.
In high flux cases, the flow has more opportunities to traverse steep portions of the S curve as it explores a larger portion of it. Though not necessarily formally unstable to the TI, these regions where temperature varies strongly with ionization parameter lead to an acceleration of the flow. High flux cases can therefore experience multiple phases of acceleration, provided they remain in radiative equilibrium and the S curve features multiple steep regions. We found this was the case with the mBl SED where in the highest flux runs had two stages of acceleration. The AMD features two dips, indicating the ionization state of the gas when it accelerates.
We have assumed an optically thin wind, though a posteriori we can estimate the validity of this approximation. Except when , all cases have column densities and therefore are optically thin throughout. The critical case is optically thin in the supersonic part. This indicates that the optically thin heating assumption is not always valid, however it can be used as a starting point for this type of calculation.
Our results suggest that complete spectral information is important when modeling radiative heating and thermally driven flows. This is especially imprtant for systems like AGN. Besides the Type 1 and Type 2 AGN cases, we may also interpret the BB case as an AGN which has nearby cold gas that absorbs the soft photons. The AGN1 and AGN2 cases have a difference of a few in their various wind fluxes and efficiencies. However the BB case is several orders of magnitude smaller, primarily due to highly suppressed soft photons. When modeling observed AGN outflows this is important because any uncertainty in the observed flux will result in an uncertainty in the modeled outflow.
We found numerical artifacts affecting our solution in several ways. This suggests the need for new physics to properly resolve the flow. When the wind becomes non-adiabatic and falls of the equilibrium curve, there is a drop in the length scale associated with the dynamical variables. This becomes more pronounced at higher radiation fluxes and can lead to not finding a stationary numerical solution. Likewise when solutions are thermally unstable there is a decrease in the length scale of dynamical variables - if it falls below the grid scale the TI is not numerically well resolved. This cannot be adressed with additional resolution as the regions are geometrically thin. Future studies of wind solutions should include the effects of thermal conduction because of relatively steep temperature gradients (see for example Rozanska & Czerny 1996).
Our future studies will include 2D and 3D simulations of disc winds driven by the radiation force. In addition to a considerably more complex geometry and dynamics due to a rotating disc for instance, accurate modeling of the radiative heating is also critical. If the disc is irradited by a central source or the inner parts of the disc, the outer disc will heat. This affects the strength of the radiation force and therefore the radiatively driven wind. However it has been shown that shadowing of the source by an inner disc wind is an important effect in this type of setup (e.g., Proga, Stone & Kallman 2000 and Higginbottom et al. 2014). If line driving is important, calculating the force multiplier resulting from the lines responsible for the heating/cooling as well as additional lines that are negligible energetically but important for momentum transfer. By using observed SEDs as an input to XSTAR we will model self-consistently the heating/cooling in the system. We also note that this method is quite general and can be applied to other radiatively heated systems for which an SED has been observed. For example one could consider the net heating of planetary atmospheres by the host star.
This work was supported by NASA under ATP grant NNX14AK44G. We thank Maria Diaz Trigo for providing the XRB SEDs and Missagh Mehdipour for providing the AGN SEDs used in our work. We also thank Tim Kallman for providing us with XSTAR and many valuable discussions.
-  Alexander, R., Clarke, C., Pringle, J., 2006, MNRAS, 422, L82.
-  Adhikari, T.P., Rozanska, A., Sobolewsa, M., Czerny, B., 2015, ApJ
-  Bally, J., Scoville, N.Z., 1982, ApJ, 255, 497.
-  Basko, M.M., Sunyaev, R.A., 1973, Ap&SS, 23:117B.
-  Basko, M.M., Sunyaev, R.A., Hatchett, S., McCray, R., 1977, ApJ, 215:276B.
-  Bautista, M.A., Kallman, T.R., 2001, ApJ, 134:139-149.
-  Begelman, M.C., McKee, C.F., Shields, G.A., 1983, ApJ 271:70-88.
-  Blondin, J.M., 1994, ApJ, 435, 756.
-  Buff, J., McCray, R., 1974, ApJ 188:L37
-  Field, G.B., 1965, ApJ 142:531.
-  Gardiner,T.A., Stone. J.M., 2005, J. Comput. Phys., 205, 509
-  Gardiner,T.A., Stone. J.M., 2008, J. Comput. Phys., 227, 4123
-  Higginbottom, N., Proga, D., Knigge, C., Long, K.S., Matthews, J.H., Sim, S.A., 2014, ApJ, 789:19H.
-  Higginbottom, N., Proga, D., 2015, ApJ, 807:107.
-  Hollenbach, D., Johnstone, D., Shu, F., 1993, ASP Conf. Series, 35:26.
-  Holczer, T., Behar, E., Kaspi, S., 2007, ApJ 663: 799-807.
-  Holzer, T.E., Axford, W.I., 1970, An Rev. Astr. Ap. 8, 31.
-  Johnstone, H., Hollenbach, D., Bally, J., 1998, ApJ, 499, 758.
-  Kallman, T.R., McCray, R., 1982, ApJ 50: 263-317.
-  Lamers, H.J.G.L.M., Cassinelli, J.P., 1999, Introduction to Stellar Winds, Cambridge UP.
-  Liu, F.K., Meyer, F., Meyer-Hofmeister, E., 1995, A&A, 300, 823-833.
-  Lodders, K., Palme, H., Gail, H.P., 2009, Landolt Brnstein (Springer), 44
-  London, R., McCray, R., Auer, L.H., 1981, ApJ, 243:970L.
-  Luketic, S., Proga, D., Kallamn, T., Raymond, J., Miller, J., 2010, ApJ, 719, 515.
-  Mehdipour, M., Kaastra, J.S., Miller, M.C., et al., 2015, A&A 575:A22
-  Owen, J., Clarke, C.J., Ercolano, B., 2012, MNRAS, 422, 1880-1901.
-  Owen, J., Ercolano, B., Clarke, C., Alexander, R., 2010, MNRAS, 401, 1415.
-  Parker, E.N., 1958, ApJ, 128:664.
-  Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., Numerical Recipes in C, Cambridge UP, 1992.
-  Proga, D., Stone, J.M., Kallman, T.R., 2000, ApJ, 543:686-696.
-  Rozanska, A., Czerny, B., 1996., Acta Astron., 46, 233.
-  Sazonov, S.Y., Ostriker, J.P., Ciotti, L., Sunyaev, R.A., 2005, MNRAS, 358: 168-180.
-  Sazonov, S.Y., Ostriker, J.P., Sunyaev, R.A., 2004, MNRAS, 347:144-156.
-  Stone, J.M., Proga, D., 2009, ApJ 694:205-213.
-  Trigo, M.D., Boirin, L., Migliari, S., Miller-Jones, J., Parmar, A., Sidoli, L., 2013, IAUS 290:25T
-  Verner, D.A., Barthel, P.D., Tytler, D., 1994, A&A 108:287-340
Appendix A Heating Function
Our Athena++ simulation calculates radiative heating rates self consistently by interpolating the net heating rates precalculated using the XSTAR photoionization code for a range of temperatures and ionization parameters. We describe this interpolation scheme in Section A.1 and the results of various performance tests in Section A.2.
We calculate the net heating rate using XSTAR for an logarithmically spaced grid in a domain and . We take , , and . We will refer to this as our heating table.
Within Athena++, at every cell and timestep, the net heating must be calculated. This requires two interpolations: 1) In phase space to determine the nearest points in the heating table 2) In T space to implicitly solve for the heating rate consistent with the hydro time step.
a.1.1 Space - Bilinear Interpolation
When we require the heating rate for ionization and temperature we 1) Determine the (i,j) location in the heating table for the four pairs of points , , and nearest to . 2) Bilinearly interpolate between these points to estimate .
We have generated the heating table uniformally in logarithmic space according to
Therefore when we are finding the indices we do not have to traverse the entire table we can simply invert the above expressions. We then use the standard bilinear interpolation scheme on these four points.
a.1.2 Space - Backward Euler Interpolation
We discretize the energy evolution equation (2c) under the assumption that changes in internal energy are due to the heating/cooling and express this in terms of the implicit relation
This is the backward Euler method where the heating rate during the th time step is calculated using the temperature at the th time step. The alternative forward Euler method calculates the heating rate during the th time step using the temperature at the th time step. The advantage of backward Euler method is that as the equilibrium temperature is reached the contribution from the heating will also be decreasing so equilibrium is reached more smoothly. The disadvantage is is now implicitly defined by (14) and a root finder is needed to solve for it.
We implemented a two step bisection root finder. The process of finding a zero of a function i.e solving for such that consists of two steps - 1) bracketing the zero, i.e. finding numbers a,b such that . A zero will then lie in the interval [a,b]. 2) Refining this interval such that so that we can say that we have located the location of the zero to a precision .
Step 1) consists of checking for the existence of a root in progresively larger intervals where and where and T is the current temperture. We then check the condition and if it is not satisfied continue decreasing (increasing) () until it is satisfied. Step 2) uses the zbrent rootfinding algorithm (Press et al. 1992) to solve for the root to an accuracy .
We compared how both methods performed when calculating the heating (cooling) of a box of uniformally ionized gas at temperature below (above) its equilibrium temperature. The backwards Euler method did not oscillate above/below the equilibrium temperature as with forward Euler. This is an important property because we expect wind solutions to be at equilibrium near their base. A comparable level of performance by the forward Euler method requires decreasing the timestep, which in practice entails more computational expense than employing the root finder.
To test our algorithm and implementation we use the analytic heating/cooling rates given by (8) (Blondin 1994). We compare simulation results where this forumula is used explicitly and simulations where our algorithm is used to interpolate over heating/cooling rates calculated on an grid.
We perform two main tests 1) Compute the relative error in the heating/cooling rate at different locations in the () phase space to test the bilinear interpolation scheme. 2) Compute the relative error in the kinematic quantities and heating rates for radiatively heated 1D stationary winds in Athena++.
a.2.1 Bilinear Interpolation
We consider grids of points , logarithmically spaced in the plane. We consider grids of size N = 50, 100 & 200. We then evaluate the heating/cooling rate in a logarithmically spaced grid of size in the plane using the analytic expression (8) and by using our bilinear interpolation scheme on the grid. We ensure this grid is shifted relative to the grid used to calculate heating rates so we are testing the interpolation scheme.
We show results for N = 50, 100 & 200 in Fig. 7. We have plotted the relative error in the numerically interpolated net heating rate (subscript n), normalized to the average analytic heating and cooling (subscript a)
We normalize the error in this way to not exagerate the importance of numerical errors near or at equilibrium where .
We note that the error contours are stratified, particularly at small N. This is because near the calculated points the bilinear approximation will perform best. As we distance ourselves from this point and before we approach the next point on the grid the approximation will worsen.
We see that if we want to achieve roughly 1% accuracy we need roughly a grid. In order to achieve 0.1% accuracy we need a grid of .
a.2.2 Stationary Solution
We consider grids of points , logarithmically spaced in the plane. We consider grids of size N = 10, 50, 100 & 500. We perform Athena++ simulations where the heating/cooling rate is determined by the bilinear interpolation scheme on this grid. After Athena++ has found a stationary solution, where the hydrodynamic time step becomes constant, we compare these to the stationary solutions found when heating/cooling is computed analytically. We plot the relative % error in the late time density, velocity, temperature and heating rate in Fig. (8).
Relative errors are large near the star where the solution is subsonic. However, past the sonic point the relative error shows a clear trend in the accuracy of the simulations. In particular a grid yields roughly 1% accuracy while a grid yields a 0.1% accuracy. Interestingly the accuracy in the final solution is comparable to the accuracy in the plane discussed above.
We note that the solution follows the equilibrium curve near the star before falling off and staying slightly below it in the heating regime. The location where the wind is no longer in equilibrium corresponds to the radius at which the relative % error converges to its large radius value. It makes sense that inside this radius the % error has large variations, because the analytic value should be near zero and therefore the relative errors are amplified.