Issue 
A&A
Volume 638, June 2020



Article Number  A107  
Number of page(s)  10  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201936895  
Published online  19 June 2020 
Simulations of stellar winds from Xray bursts
Characterization of solutions and observable variables
^{1}
Departament de Física, EEBE, Universitat Politècnica de Catalunya, c/Eduard Maristany 16, 08019 Barcelona, Spain
email: gloria.sala@upc.edu
^{2}
Institut d’Estudis Espacials de Catalunya, c/Gran Capità 24, Ed. Nexus201, 08034 Barcelona, Spain
Received:
11
October
2019
Accepted:
29
April
2020
Context. Photospheric radius expansion during Xray bursts can be used to measure neutron star radii and help constrain the equation of state of neutron star matter. Understanding the stellar wind dynamics is important for interpreting observations, and therefore stellar wind models, though studied in past decades, have regained interest and need to be revisited with updated data and methods.
Aims. Here, we study the radiative wind model in the context of Xray bursts with modern techniques and physics input. We focus on characterization of the solutions and the study of observable magnitudes as a function of free model parameters.
Methods. We implemented a spherically symmetric nonrelativistic wind model in a stationary regime, with updated opacity tables and modern numerical techniques. Total mass and energy outflows (Ṁ, Ė) were treated as free parameters.
Results. A highresolution parameterspace exploration was performed to allow better characterization of observable magnitudes. High correlation was found between different photospheric magnitudes and free parameters. For instance, the photospheric ratio of gravitational energy outflow to radiative luminosity is directly proportional to the photospheric wind velocity.
Conclusions. The correlations found here could help determine the physical conditions of the inner layers, where nuclear reactions take place, by means of observable photospheric values. Further studies are needed to determine the range of physical conditions in which the correlations are valid.
Key words: Xrays: bursts / stars: winds, outflows / stars: neutron / methods: numerical
© ESO 2020
1. Introduction
Type I normal (short) Xray bursts (XRBs) are highly energetic and recurrent thermonuclear events occurring on the envelope of accreting neutron stars in binary systems where the secondary star is usually a main sequence star or red giant. Most observed XRBs have short orbital periods in the range 0.2−15 h^{1}. As a result, the secondary star overfills its Roche lobe and masstransfer ensues through the inner Lagrangian point (L1) of the system. The material stripped from the secondary has angular momentum such that it forms an accretion disk around the neutron star. Viscous forces then progressively remove angular momentum from the disk forcing the material to spiral in and pile up on top of the neutron star. The accreted material accumulates under mildly degenerate conditions, driving a temperature increase and the onset of nuclear reactions. As a result, a thermonuclear runaway occurs, generating a massive luminosity increase as well as nucleosynthesis of heavier elements, mostly around A = 64 (see, e.g., José et al. 2010; Fisker et al. 2008; Woosley et al. 2004).
Type I XRBs are observationally characterized by a quick and sharp luminosity rise in about 1 − 10 s, a burst duration ranging from ∼10 to 100 s, a total energy output of about 10^{39} erg, and recurrence periods from hours to days. The luminosity can sometimes rise up to one thousand times the persistent (accretion) luminosity, although the typical increase is in the hundreds of times this latter, reaching values of the order of 10^{38} erg s^{−1}. Another observable based on the light curve is the ratio between timeintegrated persistent and burst fluxes, α, typically in the range ∼40 − 100, corresponding to the ratio between the gravitational potential energy released by matter falling onto the neutron star during the accretion stage (GM_{ns}/R_{ns} ∼ 200 MeV per nucleon) and the nuclear energy generated in the burst (about 5 MeV per nucleon, for a solar mixture burned all the way up to the Fegroup nuclei). The presence of heavy elements can, theoretically, be detected (see Weinberg et al. 2006; Chang et al. 2005, 2006; Bildsten et al. 2003) in the form of absorption features in the spectrum, which mostly lies in the Xray range. For further information on XRBs, see Strohmayer & Bildsten (2006), Keek & in’t Zand (2008), Galloway et al. (2008), José (2016).
The mechanism powering XRBs bears a clear resemblance to that for classical novae, but unlike these latter, the high surface gravity of the neutron star prevents, in principle, the explosive ejection of material. However, the luminosity can approach or even exceed the Eddington limit, which may lead to the ejection of some material through a radiationdriven wind.
Stellar winds have been studied in different scenarios throughout most of twentieth century, and in a variety of forms (see Parker 1965; Żytkow 1972; Castor et al. 1975). The simplest models assume spherical symmetry and stationary wind and can be broadly classified according to the main driving mechanism, namely gas pressure or radiation, although magnetic fields can also play an important role. For a wind to be radiatively driven, high luminosity and high opacity must be present.
In the framework of neutron stars, several studies of radiationdriven winds have been performed since the early 1980s with varying hypotheses, approximations, and calculation techniques. Ebisuzaki et al. (1983), Kato (1983), Quinn & Paczynski (1985) and Joss & Melia (1987) all used nonrelativistic models with an approximated formula for opacity as a function of temperature. The models adopted different boundary conditions both at the photosphere and at the base of the wind envelope, as well as different treatments of the sonic point singularity (see Sect. 2). Studies based on general relativistic models were performed by Turolla et al. (1986), Paczynski & Proszynski (1986). In a more recent work, Yu & Weinberg (2018) used MESA code (see Paxton et al. 2011) to perform a timedependent hydrodynamic simulation of the wind envelope following a hydrostatic burst rise.
Advances in computational power and numerical techniques have allowed several studies since the 2000s to perform hydrodynamic simulations of XRBs with extended nuclear reaction networks (see Fisker et al. 2008; Woosley et al. 2004; José et al. 2010). These studies have shown that XRBs synthesize a large variety of protonrich nuclei. The potential impact of XRBs on galactic abundances is still a matter of debate and relies on the high overproduction of some particular isotopes. It has been suggested that if a tiny fraction of the accreted envelope is ejected through radiationdriven winds, XRBs may potentially be the source of some light pnuclei, such as ^{92, 94}Mo and ^{96, 98}Ru (see Schatz et al. 1998, 2001), which are systematically underproduced in all canonical scenarios proposed to date for the synthesis of such species. However, it is still not clear whether XRBs contribute to the galactic abundances, because the strong gravitational pull of the neutron star prevents the direct ejection of matter, only potentially viable through a radiationdriven wind. Additionally, the study of XRBs wind can lead to a more accurate determination of neutron stars radii, and help constrain the equation of state of neutron star matter (see Özel et al. 2010; Steiner et al. 2010). Thus, the interest in stellar winds in the context of XRBs has been renewed. For this reason, we proceeded to reanalyze these wind mechanisms, in an attempt to improve (with respect to some of the previously mentioned earlier studies from the 1980s) on some aspects of the input physics (e.g., updated opacity) and the numerical techniques (critical point treatment), and to extend the analysis of solutions, explore the parameter space in more detail, and characterize the observable magnitudes.
The paper is organized as follows. The basic wind model, input physics, and simulation setup are described in Sect. 2 and in Appendix A. A full account of the results obtained from our simulations of radiativedriven winds is presented in Sect. 3. Finally, the significance of our results and our main conclusions are summarized in Sect. 4.
2. Model, input physics, and initial setup
Without any prior assumption on the temperature profile throughout the envelope, the spherically symmetric stationary wind equations constitute a set of nonlinear coupled differential equations. The boundary conditions are often given at more than one point and with implicit expressions. Another common feature is the appearance of a singular point where the wind becomes supersonic (see Sect. 2.1). Different approaches have been implemented to deal with the numerical difficulties involved in solving the equations close to this singular point (see Żytkow 1972; Joss & Melia 1987; Paczynski & Proszynski 1986).
The simulations reported in this work rely on a stationary, nonrelativistic wind model with spherical symmetry. Neutron stars have a typical radius of a few times their Schwarzschild radius, and so while general relativity corrections on gravity could be of importance close to its surface, they diminish rapidly in a wind envelope that extends hundreds of kilometers. Paczynski & Proszynski (1986) showed that the main effects of general relativity on the wind structure can be relevant for solutions with lower mass outflows, resulting in a more extended (up to a factor x4) and cooler envelope. However, they rapidly become less significant for higher mass outflows (> 5 × 10^{17} g s^{−1}), giving corrections of only a few percent. We therefore leave the effects of general relativity for a followup work and study the Newtonian approach here as a stepping stone. The radiationdriven wind, treated as a fully ionized perfect gas, is assumed to be optically thick and in local thermal equilibrium (LTE) with radiation. For a more detailed treatment of the optically thin regions and their impact on observables, see for example Joss & Melia (1987).
In this framework, the basic hydrodynamic equations for mass (1), energy (2), momentum conservation (3), and radiative energy transport (4) become:
where the variables r, T, v, ρ, L_{R} stand for the radial coordinate, temperature, wind velocity, gas density, and radiative luminosity, respectively. The total pressure P, and specific enthalpy h, include the contributions of gas and radiation, in the form:
The opacity κ(T, ρ) must be provided externally either from a radiation transport theoretical model, phenomenological relations, or experimental values. Opacity can have different sources depending on the microscopic processes involved; in general, these depend on photon frequency. The tables calculated by OPAL (Iglesias & Rogers 1996) take these aspects into account to give a Rosseland mean opacity as a function of temperature, density, and composition of the gas only. These tables have become the standard source for stellar opacity in recent decades and are used in the present work^{2}.
The total mass and energy outflows (Ṁ, Ė) arise as integration constants from the mass and energy conservation laws and can be determined by imposing conditions at the base of the wind. Therefore, Ṁ, Ė are model parameters. During an Xray burst these will vary in time, but are assumed to be approximately constant across the wind for a fixed time. The mass of the neutron star core M is considered fixed and constitutes the only relevant source of gravity (the contribution of the envelope mass can be neglected). Lastly the mean molecular mass μ is a function of the mass fractions X_{i} of the different species present in the envelope. The rest of the symbols have their usual meaning: G for gravitational constant, c for speed of light, k for Boltzmann constant, for radiation energy density constant (σ is the StefanBoltzmann constant) and m_{A} for atomic mass unit.
By differentiating the mass conservation Eq. (1), for a constant Ṁ, one can obtain:
One can also calculate, from Eq. (5), the pressure differential appearing as the last term in the momentum Eq. (3):
and by using Eq. (6) to replace the density differential:
After some regrouping, both differential Eqs. (3) and (4) can be put in the form:
where we introduced the luminosity ratio (L_{Edd} is the local Eddington luminosity) and the density and radiative luminosity can be obtained using the mass and energy conservation Eqs. (1) and (2):
These constitute a system of four equations, two algebraic (11) and (12) and two firstorder ordinary differential Eqs. (9) and (10) that can be solved for v, T, ρ, L_{R} as functions of the independent variable r.
2.1. Boundary conditions
Two boundary conditions are required to solve the system of two firstorder differential equations. The first condition is given by the definition of the photosphere, where the gas temperature is equal to the effective temperature T_{eff}, and can be written as:
which shall be referred to as “temperature condition”. Hereafter, subindex “ph” denotes evaluation at the photosphere.
The location of the photosphere in radius is not known a priori so this adds an extra unknown parameter. This is usually given by the point at which the optical depth is about unity, which would require the integration of the atmosphere above the photospheric surface. In the case of an extended atmosphere it is usual to take a locally defined effective optical depth τ^{*} = κρr that takes a minimum at the photosphere with an approximate value of (see Kovetz 1998). Kato & Hachisu (1994) showed that this minimum usually corresponds to the point where this effective optical depth takes a value closer to . Accordingly, we adopt the following “optical condition”:
These conditions, (13) and (14), when replaced in the radiation diffusion Eq. (4) give rise to a condition for the temperature gradient at the photosphere:
which can be used alternatively instead of either the optical or temperature conditions.
The second boundary condition arises from analyzing the topology of solutions to the momentum equation. From Eq. (9), the wind velocity gradient can be expressed as:
from where it is clear that the momentum equation presents a singularity when:
that is, when the wind velocity is equal to the local isothermal sound speed (hereafter, singularity condition). If the conditions were such that this singular value were met, in order to have a nonsingular physical value for the velocity gradient we must simultaneously require the following “regularity condition”:
so that the quotient is finite. A solution may or may not pass through this “critical point” a priori, but the only physically acceptable solution with a wind velocity that approaches zero towards the bottom of the envelope and always increases outwards is the one that indeed does pass through this point. Furthermore, the wind becomes supersonic above this critical point (also known as the sonic point). It is not known a priori where this critical point will lie, but together conditions (17) and (18) provide both the extra condition needed to locate it and the second boundary condition for the integration constant. We use the subscript “cr” to refer to variables evaluated at the critical sonic point.
We therefore have four unknown parameters v_{cr}, r_{cr}, T_{ph}, r_{ph} with four conditions to solve them: two at the photosphere (temperature and optical conditions) and two at the critical point (singularity and regularity conditions). The system constitutes a twopoint boundaryvalue problem with nonlinear firstorder differential equations. Standard numerical methods for dealing with them are the relaxation method and the shooting method.
However, the resulting solutions may not always reach physical values compatible with those of a neutron star at the wind base. For instance, a high temperature may be found at radii much larger than those of typical neutron stars; on the other hand, velocity may not be realistically low enough (or temperature/density high enough) when reaching a typical neutron star radius. Previous works relied on simplified assumptions for the wind base conditions, such as for example a fixed value of density or temperature, or an energy generation rate given by integration of simple nuclear models. We find that these choices for a wind base may not always be compatible with some of the recent and more realistic hydrodynamic simulations for Xray bursts (see José et al. 2010).
In order to accurately determine the conditions at the wind base, hydrodynamic simulations of the nuclear burning layers are likely required, but this is clearly out of the scope of this paper. Here, we search for a condition for the wind base that is realistic for a neutron star undergoing an Xray burst, while at the same time leave some margin for the natural variability of the physical conditions. Setting a range of values for all the physical variables seems like a good approach at first, but the choice of several restriction values seems arbitrary and turned out to not always be compatible with selfconsistent solutions. For instance, a particular choice of maximum velocity and minimum temperature at the wind base sometimes ended up giving solutions that were not deemed stationary (according to the characteristic time criterion explained in Sect. 2.2). Instead, we chose a criterion that is physically based, gives similar wind structure for all solutions, and reduces the amount of variables to be restricted, thus reducing arbitrariness, while at the same time still being compatible with the expected conditions for Xray bursts. We set the wind base at the innermost point where the radiation pressure gradient becomes larger than the gas pressure gradient, that is where:
Above this point, radiation pressure becomes more important than gas pressure as the driving mechanism. Furthermore, in all solutions explored, the velocity gradient term in the momentum conservation Eq. (3) is still negligible where the equality in Eq. (19) is met, and so the physical conditions can be considered to approximately match those of a static envelope. It can also be shown that at this point .
Solutions that meet the wind base condition in Eq. (19) at radii compatible with possible neutron stars (7 − 20 km) are then selected and considered as possible candidates for further study. In Sects. 3.1 and 3.2 we discuss the validity of these solutions for describing XRBs.
2.2. Numerical procedure
The integration method used to solve the differential equations is the adaptivestep RungeKutta (RK45). In order to deal with the twopoint boundaryvalue problem (photosphere and critical point) we use the shooting method as follows. Integration starts at a critical point with radius r_{cr} and velocity v_{cr} satisfying both singularity and regularity conditions for a temporarily assumed T_{cr}. Integration proceeds outwards until the optical photospheric condition is met (κρr ≃ 8/3). There we parametrize a distance to the photospheric temperature condition with a value ϕ:
and then use a suitable rootfinding method for ϕ in terms of the chosen starting value for T_{cr}. The numerical difficulties of starting the integration from the critical point were solved through a change of variables that cancels the singular denominator in the velocity gradient (see Appendix A).
The integration from the critical point to the photosphere corresponds to the region where the wind is supersonic. Once a supersonic profile is found, we go back to the critical point and proceed integrating inwards. The integration is stopped either when the temperature exceeds a given limit (T > 10^{10} K) or the Schwarzschild radius is reached.
The solutions are checked for some hypothesis consistency and classified accordingly, in a similar fashion to that in Quinn & Paczynski (1985), as follows. In order for the wind to be considered thick, most of the acceleration must happen at high optical depth. It is enough to require that the sonic point lies at τ ≫ 1. Solutions for which the effective optical depth of the critical point is too small () are deemed as optically “thin”. The stationary hypothesis is also checked by comparing the dynamical characteristic times Δt above and below the critical point, where
which is an estimate of the travel time of a fluid element between arbitrary radii r_{1} and r_{2}. For the solutions to be considered stationary, Δt_{above} ≪ Δt_{below}, which by use of mass conservation translates into the Δm_{above} ≪ Δm_{below} requirement for the envelope mass. Solutions with Δm_{above} > 0.1 Δm_{below} were considered to be nonstationary. We note that the value of Δm_{below} depends on the choice of the inner boundary for the wind base, and the classification of a particular solution may change with a different cutoff for the integration inwards or when imposing the nuclear burning shell conditions. In every solution found, the photospheric velocity never exceeded a few percent of the speed of light (v_{ph} ≲ 0.03c), and so the nonrelativistic hypothesis is safe. The hypothesis consistency criteria mentioned above are fairly common and can be found in several works (see e.g. Joss & Melia 1987; Kato & Hachisu 1994).
The procedure was repeated for different values of the wind parameters (Ṁ, Ė). In all the simulations reported in this work, a neutron star mass of M = 1.4 M_{⊙} has been adopted, together with an envelope composition given by X = 0, Y = 0.9, Z = 0.1. Although recent simulations of XRBs (José et al. 2010; Fisker et al. 2008; Woosley et al. 2004) show higher metallicities, the aforementioned values were the highest available on the opacity tables used. More realistic values will be implemented in future work as opacity data become available.
3. Results
Simulations were run for a wide range of values in the wind parameter space (Ṁ, Ė), exploring over 1000 points. We normalized the energy output Ė in terms of a constant Eddington luminosity
where κ_{o} = 0.2 cm^{2} g^{−1} is the electron scattering opacity, and X_{H} is the hydrogen abundance. For a neutron star with M = 1.4 M_{⊙} and no hydrogen in the wind envelope, L_{o} ≃ 3.51 × 10^{38} erg s^{−1}.
3.1. Wind profiles
A suite of different wind profiles is displayed in Fig. 1. Each plot presents a selection of solutions that have a wind base compatible with a neutron star radius of R_{ns} = 13 km. The first two plots in Fig. 1 are solutions for velocity and temperature, respectively. The third plot is a profile of characteristic time Δt, which is the time it takes for a fluid element to reach the photosphere from a given radius. The bottom plot is the luminosity ratio Γ = L_{R}/L_{Edd} (where L_{Edd} is the local Eddington luminosity).
Fig. 1. Wind profiles obtained that are compatible with a neutron star radius of 13 km (vertical dashed line), with different values of parameters (Ṁ, Ė). Values of mass outflow Ṁ are indicated by the line color, and values of energy outflow Ė, in units of L_{o} (see text), are indicated with labels next to a black circle (•), which marks the location of the critical sonic point in each curve. Top to bottom: velocity, temperature, characteristic time, and luminosity ratio Γ, all presented as a function of radius. All curves end at the photosphere. 

Open with DEXTER 
The general pattern observed is that solutions with smaller energy outflows are characterized by larger radii, but also a higher mass output and longer characteristic timescales. This is due to the fact that a fluid element requires longer distances to accelerate to a sufficiently high velocity when its mass is larger or the available energy is smaller. Profiles with higher mass outflow (lower energy outflow) may not be suitable for describing short XRBs, since the total characteristic time (from wind base to photosphere) would be larger than the typical burst duration.
Another significant behavior is that the radiative luminosity starts very low in comparison with the local Eddington limit at low radii, suggesting that the wind is driven by gas pressure rather than radiation in the inner regions. The luminosity ratio Γ = L_{R}/L_{Edd} then rises sharply, before reaching the critical point and flattens very close to unity after that, remaining flat for the rest of the supersonic profile, except for some minor bumps (about 3% maximum excess) in some models close to the photosphere.
3.2. Parameter space exploration
Figures 2–4 show photospheric, critical point and wind base values, respectively, as a function of model parameters (Ṁ, Ė). Temperature and radius are plotted in all three figures. Luminosity ratio Γ is plotted in photosphere and critical point figures (Γ_{wb} ≃ 0.5 is constant). Wind velocity is only shown for the photosphere, because at the critical point, , and wind base values are negligible with values smaller than a few centimeters per second. Other relevant plots presented are effective optical depth τ^{*} for critical point only, and wind base density and characteristic time. All figures show both the acceptable (selfconsistent) solutions and those discarded for being either optically thin or nonstationary solutions, according to the criteria explained in Sect. 2.2. The subset of acceptable solutions that are compatible with possible neutron stars is also indicated, as well as the solutions plotted in Fig. 1. We note that the criterion for stationary wind depends on the ratio of envelope mass above and below the critical point, and thus the choice of the wind base affects this classification.
Fig. 2. Wind parameter space sweep. Color coded values of photospheric magnitudes for different values of parameters (Ṁ, Ė). From top to bottom: temperature, velocity, radius, and luminosity ratio Γ. Points marked with up or down triangles (△/▿) correspond to “thin” and nonstationary solutions, respectively. Points marked with filled circles (•) are selfconsistent solution candidates. The × symbol denotes no solutions found for the given boundary conditions. The fully colored area marks solutions whose wind base is compatible with possible neutron star radii (7 − 20 km), and white diamonds inside of it (⋄) mark selected solutions plotted in Fig. 1. 

Open with DEXTER 
Fig. 3. Same as Fig. 2, but for values at the critical (sonic) point, with the exception of the velocity plot that has been replaced by an effective optical depth τ^{*} = κρr plot. 

Open with DEXTER 
Fig. 4. Same as Fig. 2, but for values of temperature, density, radius, and characteristic time Δt at the wind base. Selfconsistent solutions outside the area of compatibility with neutron stars are not colored here in order to better resolve values in the area of interest. 

Open with DEXTER 
We observe a clear pattern: as energy output Ė increases and mass output Ṁ decreases, photospheric temperature rises, in the exact opposite fashion to the photospheric radius, which decreases. We study this correlation further below in Sect. 3.3. The photospheric wind velocity gradient with respect to parameters seems to be roughly orthogonal to that of temperature or radius, only achieving a few percent of speed of light, and peaking outside the selfconsistent solutions zone. The photospheric radiative luminosity is very close to local Eddington luminosity, as shown by the Γ plot, except for lowenergy and highmass outflows, where it peaks at up to a 3% excess. This is due to a small bump in opacity tables that lowers the local Eddington luminosity.
As for the critical point (Fig. 3), no obvious relationship between temperature and radius appears. The change in critical radius seems to go mostly inversely with Ė, while the change in critical temperature seems to have a bigger component in the Ṁ direction, especially at high Ė. By definition of critical point, a critical velocity plot would give the same information as the temperature one, instead we show a plot for effective optical depth. The critical point lies at higher optical depth with bigger mass outflow, hinting at a stretching of the envelope as mentioned above. The radiative luminosity is always close to the local Eddington limit, but just slightly lower, by no more than ∼0.1%.
The temperature and density values at the wind base (as defined by Eq. (19)) are shown in Fig. 4. Their values are in a reasonable range for XRBs conditions (see José et al. 2010; Fisker et al. 2008; Woosley et al. 2004) in the area of radial compatibility with neutron stars. However, not all of the models in this area have wind base characteristic times (i.e. the time required for a fluid element to reach the photosphere from the wind base) that are well suited for describing short XRBs, whose duration is of the order of few hundreds of seconds. These models could still represent winds in longduration XRBs. Nevertheless, it is possible that the high values of Δt found are a result of the choice of inner boundary condition, as mentioned in Sect. 2.2. In the profiles shown in Fig. 1, the velocities quickly go to zero at the wind base, and therefore wind material spends nearly all its outflow time just above the wind base. This is probably an artifact of the timeindependent treatment, which may no longer be accurate close to the inner boundary.
3.3. Photospheric correlations
A close correlation was found between several photospheric magnitudes across solutions. First, for most profiles, the photospheric luminosity was found to be almost equal to the local Eddington luminosity, as can be seen in Fig. 5 where the distribution peaks sharply around Γ_{ph} = 1, and the values expressed in terms of mean and standard deviation are: Γ_{ph} = 1.0018 ± 0.0055. This can also be seen in the color maps of the wind parameter space in Fig. 2, where the coloring for photospheric Γ is mostly homogeneous, except for a small rise in the lowenergy and highmass outflows region where a bump in opacity table causes the local Eddington luminosity value to drop a little. The value of radiative luminosity (and therefore also opacity) at the photosphere also shows little variation across models, with the mean and standard deviation being: L_{R, ph}/L_{o} = 1.025 ± 0.022 and κ_{ph}/κ_{o} = 0.981 ± 0.024.
Fig. 5. Distribution of photospheric luminosity ratio Γ across acceptable solutions with varying values of model parameters. 

Open with DEXTER 
As both photospheric κ and L_{R} show very little variation, two correlations arise from the photospheric boundary conditions. Assuming a fairly constant L_{R}, the effective temperature condition gives:
Similarly, for the optical depth condition and κ_{ph} ≃ const.,
Furthermore, by using both photospheric conditions, the equations for conservation of mass and energy, and assuming that the radiative luminosity is close to the Eddington limit, we can arrive at the following expression, which relates observable photospheric magnitudes with wind model parameters.
Such correlations can be derived considering the following. First, using mass conservation to replace Ṁ, and assuming , it is easy to see that:
This relationship is valid whenever the luminosity is equal to the local Eddington limit. Further, particularly at the photosphere, using the optical condition, we have:
This accounts for the first equality of expression (25).
Now let us take the energy conservation equation and rewrite it in the form:
where we have introduced the notation for local sound speed and for escape velocity. For the wind to be optically thick, the critical point must lie at much greater optical depth than the photosphere. Since the critical point is the place where the wind becomes supersonic, this often translates into having v^{2} ≫ s^{2} at the photosphere. Lastly, in almost every solution obtained here, photospheric escape velocity is still around one order of magnitude above wind velocity, which means s^{2} ≪ v^{2} ≪ u^{2}. The last relation is relaxed as Ṁ increases and Ė decreases. So far, neglecting the corresponding terms we have at the photosphere:
If we now use the photospheric boundary conditions rewritten in a convenient form:
by taking the product of both and dividing by 4GM, we get:
and if again L_{R} = L_{Edd}, we get at the photosphere:
By replacing this value in Eq. (29) we arrive finally at the relation
which is the second equality in expression (25).
For solutions for which the approximation v^{2} ≪ u^{2} weakens, one should add the term corresponding to wind kinetic energy :
where we used relation (27) to write the equation in terms of v/c. As v ≪ c, and photospheric luminosity varies weakly across solutions (L_{R, ph} ∼ L_{o}), the correction term will only become relevant for relatively high values of Ṁ. That is, if
4. Conclusions and discussion
In the previous section, we present several results, including a set of wind profiles (solutions for the wind model), a highresolution parameter space exploration characterizing different physical magnitudes at points of interest, and several correlations found for photospheric (observable) values.
The wind profiles found in this work (Sect. 3.1 and Fig. 1) are qualitatively similar to those found in previous works (see, e.g., Quinn & Paczynski 1985). These previous studies frequently relied on an approximated expression for the opacity, in the form:
where κ_{es} = 0.2 cm^{2} g^{−1} is the electronscattering opacity, X is the hydrogen mass fraction, α^{−1} = 4.5 × 10^{8} K, and ξ is a parameter for which different values haven been adopted (ξ = 1, Kato 1983; Ebisuzaki et al. 1983, ξ ≃ 0.86, Quinn & Paczynski 1985). Figure 6 shows a comparison of the opacity profiles obtained for several solutions reported in this work with values obtained through the approximate Eq. (36) in some of the aforementioned previous studies. Equation (36) neither takes into account the expected density dependence nor metallicity effects^{3}. According to more recent studies, a high concentration of heavy elements is expected in the envelope after XRBs (see, for example, José et al. 2010; Woosley et al. 2004; Fisker et al. 2008), and therefore high opacity. The inclusion of updated opacity tables was a necessary improvement in this regard.
Fig. 6. Opacity profiles vs. temperature. Acceptable solutions obtained in this work with wind base at R_{ns} = 13 km and different values of log Ṁ (continuous lines), compared to other prescriptions adopted in previous studies (dashed lines) with different values of ξ (see text). Opacity is normalized to κ_{es} = 0.2 cm^{2} g^{−1}. 

Open with DEXTER 
At sufficiently high temperatures (close to the wind base), opacity drops considerably in all prescriptions, making the contribution of radiation to wind acceleration less important, while favoring the effect of gas pressure. As temperature and density drop outwards, the gas pressure gradient diminishes, while the opacity increases, lowering the local Eddington limit. As a result, the wind is further accelerated and sustained by radiation. This increase in opacity is higher close to the photosphere with opacity tables than those reported in previous studies, even for the trial Z = 0.1 used in this work. The differences may seem relatively small here, namely on the order of a few percent, but they are expected to rise even more for higher metallicity. As we plan to study stellar winds with higher metallicities and different mixtures of heavy elements given by XRBs simulations in future work, we consider the inclusion of opacity tables to be an important step towards that goal.
With regard to the parameter space exploration presented in this work (Sect. 3.2, Figs. 2–4), it served two purposes. The first was to find the regions for which acceptable solutions exist for the wind model used, by testing part of the assumed hypothesis (i.e., optically thick and stationary wind), in a similar way to Quinn & Paczynski (1985). The main improvement in this regard is the number of solutions obtained in the present work with respect to these latter authors. This allows for a higher resolution, which helps, for instance, to better determine these regions. The other purpose of the parameter space exploration was the characterization of physical magnitudes at points of interest, especially at the photosphere (high resolution enabled this too). The dependence of these photospheric values on the wind parameters serves as a link between the physics of the layers close to the base of the envelope and observable magnitudes at the photosphere in the following way: The parameters explored, namely mass and energy outflow (Ṁ, Ė), are fully determined when imposing suitable boundary conditions for all physical variables at the base of the wind envelope. To this end, we would need to rely on hydrodynamic XRB models of the underlying layers, where nuclear reactions take place. We intend to do so in future work. In order for the current wind model to be applied safely and consistently, the boundary conditions provided by a given XRB model should result in mass and energy outflows that lie in the acceptable region. Once the parameters are fixed at the base by the XRB model, the corresponding wind solution gives the values of each observable magnitude at the photosphere. However, each value of a single photospheric magnitude taken separately does not uniquely determine the (Ṁ, Ė) pair.
The correlations found in Sect. 3.3 (expression (25)) relate the photospheric magnitudes (that can either be directly observed or calculated from observable ones) with model parameters (Ṁ, Ė) in a more direct and unique way. For instance, the radiative luminosity L_{R} can be estimated if the distance to the XRB source is known, while wind velocity v could be estimated from the Doppler shift in spectral absorption lines during the XRB. With both velocity and luminosity, one could calculate the total energy outflow Ė by means of Eq. (25). An estimate of photospheric temperature (e.g., by fitting a blackbody frequency distribution for the spectrum), can be used to calculate the photospheric radius by means of the r ∼ T^{−2} photospheric correlation or directly through the definition of effective temperature (Eq. (13)). Once we have r, v, and L_{R} at the photosphere, calculation of the mass outflow Ṁ is possible through Eq. (27), or from the mass conservation (Eq. (1)) with the density given by the ρ ∼ r^{−1} photospheric correlation. Furthermore, since nuclear reactions in XRBs take place at the interface of the envelope with the neutron star core, these correlations could provide a technique for determining the radius of the neutron star indirectly by observation of the photospheric magnitudes alone.
In summary, here we explore the possible solutions of a nonrelativistic radiative wind model in a typical XRB scenario at high resolution in the parameter space (Ṁ, Ė). The wind profiles suggest a transition from gas pressure (in the inner regions) to radiation pressure (as the wind becomes supersonic) as the main driving mechanism. The inclusion of updated microphysics data with higher metallicities and higher opacities, is important here because it reinforces the role of radiation in the wind acceleration, facilitating the abovementioned transition. Radiative luminosity was found to remain very close to the local Eddington limit in this regime, and at the photosphere in particular. A high correlation was found between different observable (photospheric) magnitudes and the model parameters. These correlations can be used to link observational data to the physical conditions of the underlying layers where nuclear reactions take place. However, given some of the simplifying hypothesis of the model employed (nonrelativistic regime, optically thick wind, LTE, etc.), it remains to be determined whether or not these correlations hold once more complex studies are performed. Hydrodynamic models of XRBs with nuclear reactions are required to impose boundary conditions at the base of the wind envelope and will determine the model parameters. We will leave this task for a followup study, where we will explore the inclusion of general relativistic effects that may become relevant in the inner layers, and additional regimes (e.g., different chemical abundance patterns, higher envelope metallicities).
Exceptions include GX 13+1 (592.8 h), Cir X1 (398.4 h), and Cyg X2 (236.2 h), see http://www.sron.nl/~jeanz/bursterlist.html for an updated list of known galactic type I XRB systems.
Opacity tables employed are limited to log T < 8.7. For higher temperatures, opacity was extrapolated using a formula introduced by Paczynski (1983).
See, however, Joss & Melia (1987), Paczynski (1983), for other densitydependent prescriptions which still do not include metallicity effects.
Acknowledgments
The authors would like to thank Arman Aryaeipour and Manuel Linares for fruitful discussions. This work has been partially supported by the Spanish MINECO grant AYA2017–86274–P, by the E.U. FEDER funds, and by the AGAUR/Generalitat de Catalunya grant SGR661/2017. This article benefited from discussions within the “ChETEC” COST Action (CA16117).
References
 Bildsten, L., Chang, P., & Paerels, F. 2003, ApJ, 591, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Chang, P., Bildsten, L., & Wasserman, I. 2005, ApJ, 629, 998 [NASA ADS] [CrossRef] [Google Scholar]
 Chang, P., Morsink, S., Bildsten, L., & Wasserman, I. 2006, ApJ, 636, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Ebisuzaki, T., Hanawa, T., & Sugimoto, D. 1983, PASJ, 35, 17 [NASA ADS] [Google Scholar]
 Fisker, J. L., Schatz, H., & Thielemann, F.K. 2008, ApJS, 174, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 José, J. 2016, Stellar Explosions: Hydrodynamics and Nucleosynthesis (Boca Raton, FL: CRC/Taylor and Francis) [CrossRef] [Google Scholar]
 José, J., Moreno, F., Parikh, A., & Iliadis, C. 2010, ApJS, 189, 204 [NASA ADS] [CrossRef] [Google Scholar]
 Joss, P. C., & Melia, F. 1987, ApJ, 312, 700 [NASA ADS] [CrossRef] [Google Scholar]
 Kato, M. 1983, PASJ, 35, 33 [Google Scholar]
 Kato, M., & Hachisu, I. 1994, ApJ, 437, 802 [NASA ADS] [CrossRef] [Google Scholar]
 Keek, L., & in’t Zand, J. J. M. 2008, Proceedings of the 7th INTEGRAL Workshop, 32 [Google Scholar]
 Kovetz, A. 1998, ApJ, 495, 401 [NASA ADS] [CrossRef] [Google Scholar]
 Özel, F., Baym, G., & Güver, T. 2010, Phys. Rev. D, 82, 101301 [NASA ADS] [CrossRef] [Google Scholar]
 Paczynski, B. 1983, ApJ, 267, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Paczynski, B., & Proszynski, M. 1986, ApJ, 302, 519 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1965, Space Sci. Rev., 4, 666 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Quinn, T., & Paczynski, B. 1985, ApJ, 289, 634 [CrossRef] [Google Scholar]
 Schatz, H., Aprahamian, A., Görres, J., et al. 1998, Phys. Rep., 294, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Schatz, H., Aprahamian, A., Barnard, V., et al. 2001, Phys. Rev. Lett., 86, 3471 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, ApJ, 722, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Strohmayer, T., & Bildsten, L. 2006, New Views of Thermonuclear Bursts (Cambridge: Cambridge Univ. Press), 39, 113 [Google Scholar]
 Turolla, R., Nobili, L., & Calvani, M. 1986, ApJ, 303, 573 [CrossRef] [Google Scholar]
 Weinberg, N. N., Bildsten, L., & Schatz, H. 2006, ApJ, 639, 1018 [NASA ADS] [CrossRef] [Google Scholar]
 Woosley, S. E., Heger, A., Cumming, A., et al. 2004, ApJS, 151, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, H., & Weinberg, N. N. 2018, ApJ, 863, 53 [CrossRef] [Google Scholar]
 Żytkow, A. 1972, Acta Astron., 22, 103 [Google Scholar]
Appendix A: Critical point substitution
In order to avoid the numerical difficulties involved in calculating the velocity gradient at the critical point we used a substitution. The velocity gradient in the wind equations is expressed in the form of a ratio:
with N, D being the numerator and denominator which should both become zero at the critical point in order to have a regular solution. Explicitly, the singularity condition is given by . Hereafter, the local isothermal sound speed will be denoted by s, where .
Consider the following function and its derivative:
In differential form we get:
If we now set , then and it can be easily found that:
The radial gradient of this new variable y then satisfies:
and replacing expression (A.1) for the velocity gradient, the denominator D is canceled, giving:
This derivative no longer presents a singularity and can be integrated instead of the momentum equation for the new variable y.
However, transforming back to v presents one inconvenience. The function y(x) has one global minimum at x = 1 with y_{min} = 1, which corresponds to the critical point (where v = s). From there, y increases monotonically towards both x = 0 and x → ∞ (see Fig. A.1). Its inverse function therefore has two branches:
Fig. A.1. Variable substitution. Change of variables y(x), in red, used to integrate close to the critical point (x, y) = (1, 1), and the two branches of its inverse function (green and blue). 

Open with DEXTER 
with the following properties for all y:
Care must be taken over the choice of which branch to take. It is not known a priori whether at a certain point the wind should be supersonic (x_{+}) or subsonic (x_{−}), unless integration is started from the critical sonic point in a particular direction. When integrating away from it, the value of y must increase, so one must ensure that the radial gradient takes a positive value for r ≳ r_{cr} and negative value for r ≲ r_{cr} and prevent any numerical error near the critical point that may give a different result.
All Figures
Fig. 1. Wind profiles obtained that are compatible with a neutron star radius of 13 km (vertical dashed line), with different values of parameters (Ṁ, Ė). Values of mass outflow Ṁ are indicated by the line color, and values of energy outflow Ė, in units of L_{o} (see text), are indicated with labels next to a black circle (•), which marks the location of the critical sonic point in each curve. Top to bottom: velocity, temperature, characteristic time, and luminosity ratio Γ, all presented as a function of radius. All curves end at the photosphere. 

Open with DEXTER  
In the text 
Fig. 2. Wind parameter space sweep. Color coded values of photospheric magnitudes for different values of parameters (Ṁ, Ė). From top to bottom: temperature, velocity, radius, and luminosity ratio Γ. Points marked with up or down triangles (△/▿) correspond to “thin” and nonstationary solutions, respectively. Points marked with filled circles (•) are selfconsistent solution candidates. The × symbol denotes no solutions found for the given boundary conditions. The fully colored area marks solutions whose wind base is compatible with possible neutron star radii (7 − 20 km), and white diamonds inside of it (⋄) mark selected solutions plotted in Fig. 1. 

Open with DEXTER  
In the text 
Fig. 3. Same as Fig. 2, but for values at the critical (sonic) point, with the exception of the velocity plot that has been replaced by an effective optical depth τ^{*} = κρr plot. 

Open with DEXTER  
In the text 
Fig. 4. Same as Fig. 2, but for values of temperature, density, radius, and characteristic time Δt at the wind base. Selfconsistent solutions outside the area of compatibility with neutron stars are not colored here in order to better resolve values in the area of interest. 

Open with DEXTER  
In the text 
Fig. 5. Distribution of photospheric luminosity ratio Γ across acceptable solutions with varying values of model parameters. 

Open with DEXTER  
In the text 
Fig. 6. Opacity profiles vs. temperature. Acceptable solutions obtained in this work with wind base at R_{ns} = 13 km and different values of log Ṁ (continuous lines), compared to other prescriptions adopted in previous studies (dashed lines) with different values of ξ (see text). Opacity is normalized to κ_{es} = 0.2 cm^{2} g^{−1}. 

Open with DEXTER  
In the text 
Fig. A.1. Variable substitution. Change of variables y(x), in red, used to integrate close to the critical point (x, y) = (1, 1), and the two branches of its inverse function (green and blue). 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.