Analytical Solutions for Radiation-Driven Winds in Massive Stars II: The $\delta$-slow Regime

Accurate mass-loss rates and terminal velocities from massive stars winds are essential to obtain synthetic spectra from radiative transfer calculations and to determine the evolutionary path of massive stars. From a theoretical point of view, analytical expressions for the wind parameters and velocity profile would have many advantages over numerical calculations that solve the complex non-linear set of hydrodynamic equations. In a previous work, we obtained an analytical description for the fast wind regime. Now, we propose an approximate expression for the line-force in terms of new parameters and obtain a velocity profile closed-form solution (in terms of the Lambert $W$ function) for the $\delta$-slow regime. Using this analytical velocity profile, we were able to obtain the mass-loss rates based on the m-CAK theory. Moreover, we established a relation between this new set of line-force parameters with the known stellar and m-CAK line-force parameters. To this purpose, we calculated a grid of numerical hydrodynamical models and performed a multivariate multiple regression. The numerical and our descriptions lead to good agreement between their values.


INTRODUCTION
The knowledge of stellar wind properties of massive stars is fundamental for understanding stellar evolution processes, different evolutionary scenarios and enrichment of star's nearby environments.
Accurate wind parameters (mass-loss rate and terminal velocity) are crucial for the study of the wind properties of massive stars. Insights into the physics of stellar winds are attained by studying the effects of wind parameters on the emergent line spectrum and by comparing the latter with observations. From a theoretical point of view, this implies to solve highly non-linear equations in which the radiation field and hydrodynamics are strongly coupled.
Winds of massive stars are driven by the transfer of momentum from the radiation field to the plasma by scattering processes in the spectral lines (Lucy & Solomon 1970). Currently, these winds are best described by the m-CAK theory (Castor et al. 1975;Friend & Abbott 1986;Pauldrach et al. 1986).
Generally, there are many approximations that reduce considerably the complexity of the computation of the hydrodynamic and the NLTE radiative transfer solutions. One example is the extensive use of a simple analytical approximation for the velocity field, the so-★ E-mail: ignacio.araya@umayor.cl called -law, first proposed by Lamers & Rogerson (1978). A value of 0.8 − 1.2, generally agrees very well with the m-CAK numerical hydrodynamic solution (Lamers & Cassinelli 1999). This value of is determined empirically by fitting the observed line profile with a synthetic one. This approximation has been proved to be very effective and efficient to describe the winds of O-and early B-type supergiants. However, in the case of late B-and A-type supergiants there is a clear tendency towards higher values of , even with values larger than 3, leading to inconsistencies with respect to the hydrodynamic theory (Stahl et al. 1991;Verdugo et al. 1999;Crowther et al. 2006;Lefever et al. 2007;Markova & Puls 2008;Searle et al. 2008;Haucke et al. 2018). Therefore, accurate analytical approximations of the m-CAK hydrodynamic equations are indispensable to have a self-consistent coupling between the hydrodynamics and multidimensional radiative transfer problems in moving media.
For the case of the fast regime (standard m-CAK solution), this issue was addressed by Villata (1992), Müller & Vink (2008) and Araya et al. (2014). The aim of this work is to extend the procedure of Araya et al. (2014) to the -slow regime. The -slow solution 1 , found by Curé et al. (2011), is based on the m-CAK theory, that describes the wind velocity profile when the ionization-related lineforce parameter takes higher values than the ones provided by the standard m-CAK solution (see, e.g., Lamers & Cassinelli 1999, and references therein). High values of , even larger than 1/3, which corresponds to a wind with neutral hydrogen as a trace element (Puls et al. 2000), are expected in strong ionization gradients (see also Kudritzki 2002). The -slow solution is characterized by low terminal speeds ( ∞ ) and might explain the obtained values for late-B and A-type supergiants. This solution also seems to fit quite well the observed anomalous correlation between the terminal and escape velocities found in A supergiants, as well as their corresponding wind momentum-luminosity relationship (Curé et al. 2011).
With the purpose to have an approximate solution from the hydrodynamic, Araya et al. (2014) developed an expression in terms of the stellar and m-CAK line-force parameters ( , , and ) and applied it to the fast regime. This expression, based on the works of Müller & Vink (2008) and Villata (1992), describes the line acceleration as function of the radial distance, allowing to solve analytically the hydrodynamic stationary equation of motion. The use of expressions for both radiation force and velocity profile as a function of the lineforce parameters can provide a clear view into how the line-driven mechanism is related with the hydrodynamics.
On the other hand, it is important to obtain a simple representation of the radiation force and the derived slow solutions under such different ionization conditions. Therefore, a significant contribution of this work consists in offering a quick way to generate an analytical expression to estimate mass-loss rates for these alternative wind regimes. There are currently no parametric expressions that can be used for this purpose without the need to fully solve the hydrodynamic equations.
This work is organized as follows: Section 2 presents briefly the hydrodynamic equations for line-driven winds and the dimensionless form of the equation of motion. In Section 3, the basic concepts developed by Müller & Vink (2008) are recapitulated including their line acceleration term as function of the radial distance. Then, this line acceleration term is modified with the purpose to obtain a better agreement with the -slow solution. In Section 4, a recipe to obtain the line acceleration parameters (required by the line acceleration term) is developed, based on a grid of hydrodynamic models and a multivariate multiple regression. Then, an analytical expression for the -slow solution is developed and compared with the numerical models described in Section 5. In Section 6, we give our conclusions. In addition, a recipe to derive the mass-loss rate based on our expression is provided in Appendix A.

THE STANDARD HYDRODYNAMICAL WIND MODEL
The CAK theory for line-driven winds was originally developed by Castor et al. (1975). This theory describes, for a point source, a stationary, one-dimensional, non-rotating, isothermal, outflowing wind with spherical symmetry. Adopting these assumptions, and neglecting the effects of viscosity, heat conduction and magnetic fields, the equations of mass conservation and radial momentum state: Here is the fluid radial velocity, / = is the velocity gradient and line is the line acceleration. All other variables have their standard meaning (for a detailed derivation and definitions of variables, constants and functions, see Curé 2004). The so called m-CAK theory, which include the effects of rotation and a disk-like source, was developed by Friend & Abbott (1986) and Pauldrach et al. (1986), based on a general expression from Abbott (1982) for the line force: where the coefficient (eigenvalue) depends on the mass-loss rate and the line-force parameter (see Eq. A5). ( ) is the dilution factor, 11 is the electron number density in units of 10 −11 cm −3 , and FD is the finite disk correction factor. The m-CAK line-force parameters are: , and .
The momentum equation (Eq. 2) can be expressed in a dimensionless form (see e.g., Müller & Vink 2008;Araya et al. 2014) as:ˆ= withˆ= / * ,ˆ= / andˆc rit = esc / √ 2. Here * is the stellar radius, is the isothermal sound speed,ˆc rit is the dimensionless rotational break-up velocity and esc is the escape velocity. The dimensionless line acceleration reads: Using Eq. 1 together with the equation of state for an ideal gas ( = 2 ), the dimensionless equation of motion is: In general, the calculation of the line acceleration involves the coupling of hydrodynamics with the radiative transport in NLTE. A very successful approach is to calculate the line acceleration using the Sobolev approximation. The pioneering work of Castor et al. (1975) laid the foundations of CAK theory and later improvements (m-CAK). A further description was done by Feldmeier (1998) who extended the CAK approach using a second order Sobolev approximation, i.e, line = line ( , , , ). However, in this work, to obtain an analytical expression of the -slow solution, we will use a radial dependence for the line acceleration following the methodology used by Araya et al. (2014), i.e., line = line ( ). This approach allows to obtain an analytical expression for the velocity field in terms of the Lambert function (see Section 3).

LINE ACCELERATION
In this section we review the basic concepts developed by Müller & Vink (2008, hereafter MV08) to derive, later on, a general analytical expression for the velocity profile in the frame of the -slow radiation-driven wind regime for massive stars. We demonstrate that this expression enables to integrate the equation of motion (Eq. 6) leading to an analytical expression for the -slow wind velocity profile.

The Fast Regime Approximation
In the framework of m-CAK stellar wind theory, MV08 present a mathematical expression for the line acceleration via a parameterized description that depends only on the radial coordinate. Using Monte Carlo multi-line radiative transfer calculations (de Koter et al. 1997;Vink et al. 1999) and a velocity profile from a -law, these authors computed the line acceleration. Then, the numerical line acceleration, which collect all the physically motivated mathematical properties for the radiative line acceleration term, is expressed by by the following function: whereˆ0, 1 ,ˆ0, and are the MV08 line acceleration parameters.
It is important to note that these parameters, lack of any physical meaning, and besides, are not directly related to , and parameters from m-CAK theory. Replacing Eq. 7 in Eq. 6, the dimensionless equation of motion are derived and a fully analytical velocity profile is obtained (see MV08 for details about the methodology used to obtain this solution) by means of the Lambert W-function (Corless et al. 1993(Corless et al. , 1996Cranmer 2004).
The line acceleration expression given by MV08 (Eq. 7) results in a good approximation for the m-CAK line force for ≤ 0.2, but this expression fails for -slow solutions, when 0.25. Overall, this approximation gives a poor agreement with respect to the numerical -slow solution (from m-CAK theory). The numerical solutions are obtained from the stationary hydrodynamic code H (Curé 2004). Araya et al. (2014) developed an analytical solution for the velocity of the fast wind regime in terms of the stellar and m-CAK line-force parameters combining the methodology from MV08 and the line acceleration proposed by Villata (1992). Unfortunately, this expression also fails when the line force parameter is higher than about 0.3, because in this case a term from the proposed line acceleration expression turns complex. From a mathematical point of view, high values of would require high values of in order to obtain an expression with real values, but such kind of values would be totally unphysical.

The New -slow Regime Approximation
In view of the unsatisfactory results obtained when applying the approximate description of the wind velocity for the -slow case, we decided to modify the functional form of the line acceleration given by MV08 in order to obtain a better description of the -slow wind. Thus, our proposed line acceleration is the following: whereˆ0, 1 , 2 , and are the new set of line acceleration parameters.
The new expression follows the same mathematical properties as MV08's but the inclusion of the 2 parameter yields to a better agreement with the numerical line acceleration from the m-CAK model.
Based on this new definition for the radiation force, the new dimensionless equation of motion reads: The same methodology developed by MV08 is employed to solve the new equation of motion and the solution is given through the Lambert function, with being 2 1 the Gauss hypergeometric function. Note that the constant of integration vanishes due to the subtraction between the integrals atˆandˆc. The critical (or sonic) point,ˆc, is obtained numerically making the RHS of Eq. 9 equal zero. Finally, taking into account the numerical solution from H as reference, a good agreement is obtained with our expression (Eq. 10) for the velocity profile.

LINE ACCELERATION PARAMETERS
In Araya et al. (2014) a relationship between the MV08 line-force parameters (ˆ0, 1 ,ˆ0, and ) and the stellar and m-CAK lineforce parameters was given. This relationship is an easy-to-use and versatile method to compute the velocity profile analytically, because both stellar and m-CAK line force parameters are already available for a wide range of spectral types (see, Abbott 1982;Pauldrach et al. 1986;Lamers & Cassinelli 1999;Noebauer & Sim 2015;Gormaz-Matamala et al. 2019;Lattimer & Cranmer 2021).
To derive a similar relationship, now for the -slow regime, we created a grid of m-CAK hydrodynamic models and develop that relationship applying a multivariate multiple regression (MMR Rencher & Christensen 2012;Mardia et al. 1980).

Grid of Hydrodynamic Models
We built a H grid of stellar models for -slow solutions. The grid points were selected to cover the region of the eff -log diagram where the B-and A-type supergiants are located.
For each given pair of stellar parameters ( eff , log ), the stellar radius was calculated from bol by means of the flux-weighted gravity-luminosity relationship (Kudritzki et al. 2003(Kudritzki et al. , 2008, but in addition we added 20 values for stellar radius (from 5 to 100 to 60 without rotation (Ekström et al. 2008), while the black lines correspond to the zero age main-sequence (ZAMS) and the terminal age main-sequence (TAMS). in steps of 5 ). The surface gravities comprise the range of log = 2.7 down to about 90% of the Eddington limit, in steps of 0.15 dex. We considered 22 effective temperature grid points, ranging from 9 000 K to 19 500 K, in steps of 500 K. These eff and log ranges were adopted to describe mainly the wind of intermediate and late B supergiants.
The m-CAK line-force parameters used for each set of ( eff , log ) values are given in Table 1. We considered only high values of in order to obtain -slow solutions.
Then, a huge combinations of parameters were executed in H -, considering the standard boundary condition at the stellar surface, for the optical depth, * = 2/3. In addition, it is worth noting that only some combinations of all parameters used in H converged to a physical stationary solution, i.e., we obtained 141 067 -slow solutions from our initial set (about a 2% of our initial input). In the eff -log plane, see Fig. 1, we show in green dots all converged models, whereas blue dots indicates that no -slow solution was achieved for the given combination of parameters. Furthermore, the number of converged -slow solutions, in the eff -log plane, shows that most of the models are concentrated in the region of log ≥ 1.65, with a peak around eff = 14 kK and log = 2.4. Also, few models converged with values of ≤ 0.28 and ≥ 0.57. This behavior must be considered at the moment to define the limits of our approximation for -slow solutions.
Finally, for each hydrodynamic model we fitted (Least Squares) the m-CAK line acceleration (ˆl ine ) with our proposed line acceleration expression (Eq. 8) in order to obtain the corresponding new line acceleration parameters (ˆ0, 1 , 2 , and ).

Multivariate Multiple Regression
To derive the relationship for the new line acceleration parameters (ˆ0, 1 , 2 and ) as function of stellar ( eff , log g, * / ) and m-CAK line-force parameters ( , , ) a MMR is applied to our grid of models.
A multiple multivariate regression model is: where is a × matrix of data in the dependent variables, is a × (1 + ) matrix of regression: a first column of 1's and in the remaining columns the data of the independent variables, is a (1 + ) × matrix of parameters (the intercept and parameters, one for each of the independent variables), and is a × matrix of measurement error. The model is the same for each dependent variable ( , =1, . . . , ), but with different coefficients ( , =0, . . . , ; =0, . . . , ), i.e., = 0 + 1 eff + 2 log g + 3 * / + (14) where represents the measurement errors. Each row of represents an observation of each of the measured response variable. Additional assumptions in the model are that the expectation of is given by ( ) = or ( ) = 0, and the covariance matrix of the vectors in the rows of is Σ, that is, the columns in can be correlated. Also, there is an assumption of normality about the response variables that allows to perform the hypothesis testing in regression.
For our problem, the dependent variables areˆ0, 1 , 2 and , and the independent variables are eff , log g, * / , , and . The database has = 141 067 records.
A data transformation is necessary to obtain a good fit of the linear model. Thus, a Box-Cox transformation (Seber &   explained by the regression) given in Table 2. Therefore, the regression explains almost all the variability ofˆ0 .27 0 , a large amount of the variability of ( 1 + 1) 5.3 and ( + 1) −3.56 , and a minor proportion of 0.45 2 .
This new relationship for the line acceleration parameters (ˆ0, 1 , 2 and ) as function of stellar and m-CAK line force parameters is valid only for -slow solutions, specifically for values of between 0.29 and 0.35. Therefore, it cannot be compared or used with others parametrizations obtained using an approximation for the velocity profile of fast solution (see e.g. Muĳres et al. 2012).

THE APPROXIMATIVE SOLUTION
Once we know the relationship (estimated model) between the line acceleration parameters as a function of the stellar and m-CAK lineforce parameters, we can use Eq. 10 to obtain the velocity profile of the -slow wind in terms of the Lambert W-function.
We point out that considering the number of converged models for some values of and , we limit our approximation to values of between 0.45 and 0.55, and values of between 0.29 and 0.35. In addition, we could expect a lower precision for values of log lower than 1.65.
In the following of this section, we discuss the accuracy of the terminal velocities and the derivation of mass-loss rates obtained using this analytical treatment.

Terminal Velocity
To measure the goodness of fit of the estimated model, the terminal velocity obtained by H is compared with our formulated solution.
We consider two terminal velocity vectors: H ∞ defined as the terminal velocity calculated with H (hereafter "true terminal velocity") and A ∞ as the terminal velocity obtained from the our solution atˆ= / * = 100, i.e., The relative error of the estimated terminal velocity A ∞ with respect to the true terminal velocity is calculated by: We obtain that the 0.90 quantile of the distribution of the relative error are below 21%, and the 0.95 quantile of them are below 27% ( 0.95 = 27.32).

Mass-loss Rate
Although our solution is developed to obtain a wind velocity profile, we can derive a recipe to obtain a mass-loss rate. This recipe is based on the m-CAK theory, specifically the work of Curé (2004), where the velocity profile is described by our proposed solution (Eq. 10). The full procedure is explained in Appendix A.
Then, similar to the procedure performed for the terminal velocity, we measure the goodness of fit of the estimated mass loss rates of the models by comparing the values (vector) calculated with H , H , and the ones obtained with our solution, A .
The relative error of the estimated mass-loss rate A with respect to the true mass-loss rate H was calculated analogously to the velocity error (Eq. 20). In comparison with the terminal velocities, the mass-loss rates have slightly higher relative errors. We observe that most of the data are below ∼ 63%. The 0.90 quantile of the distribution of the relative error is below 39% and the 0.95 quantile is about 46% ( 0.95 = 46.40).
Finally, the recipe for the calculation of the estimated mass-loss rate as a function of stellar and m-CAK line-force parameters ( eff , log , * / , , and ) is the following: (i) Computeˆ0, 1 , 2 and from Eqs. 15 through 18, calculating their respective inverse functions.
(ii) Calculate ( ) from the analytical expression given in Eq. 10.
(iii) Obtain A using ( ) from (ii) and its gradient in the m-CAK theory (see Appendix A).

DISCUSSION AND CONCLUSIONS
In the frame of the -slow wind regime, we have proposed a new approximate expression for the line force based on the MV08 methodology. This new expression is a pure function of the radial coordinate and depends on the following parameters:ˆ0, 1 , 2 , and . With this line-force we derived an analytical expressions for the velocity profile, terminal velocity and a recipe for mass-loss rate (based on m-CAK theory and our velocity approximation). Furthermore, after generating a grid of hydrodynamic models, we apply a multivariate multiple regression to obtain a relationship among these new line-force parameters with the stellar ( eff , log , and * / ) and m-CAK line-force parameters ( , , and ).
The m-CAK line force parameters should be in principle selfconsistently calculated coupling the hydrodynamics with the contribution to the line-acceleration from hundreds of thousand spectral lines (Lattimer & Cranmer 2021;Gormaz-Matamala et al. 2019;Pauldrach 2003, and references therein). This type of calculations has not been performed for the -slow regime, so far. Nevertheless, based on preliminary line-profile fittings, using the -slow solution, Cidale et al. (2017) found that the value of is in the same range as in the fast regime, while is a factor 2-3 lower.
Notwithstanding we can perform a test of our solution using and parameters from the fast regime. To this purpose, we consider stellar and wind parameters from the work of Curé et al. (2011), where they explore the influence of ionization changes throughout the wind in the velocity profile for theoretical models of A-type supergiant stars. Thus, we select the models that match our grid extension (dismissing the region where ≤ 0.28 and ≥ 0.57) in order to compare it to our expression. In addition, with purpose to test the full range of our work, we also consider the parameters from models of Venero et al. (2016), where they perform a numerical study of hydrodynamic solutions within the -slow domain, based on fundamental parameters of typical B supergiants stars.
The stellar and wind parameters from the mentioned works are listed in Table 3. This table also gives the values of the mass-loss rate and terminal velocity obtained from our analytical solution together with those values calculated from hydrodynamic results (H code). All hydrodynamic models are calculated without stellar rotation.
The accuracy of our approach is reflected in the low relative errors for the mass-loss rate and terminal velocity obtained from our solution and the hydrodynamical code. For the terminal velocity we obtain a relative error mean and median of 15.6% and 15%, respectively. In the case of the mass-loss rate a relative error mean and median of 10.4% and 7.8% are obtained, respectively.
The use of approximate expressions that describe closely the hydrodynamics of stellar winds give the advantage of solving the radiative transfer problem for moving media in an easy way. In particular, this new expression might properly describe the winds of late B-and A-type supergiants, without considering a -law with high values ( 3) that lack of any physical justification in the frame of m-CAK fast solution.
The new expressions for the -slow solutions together with the previously derived expression for the fast solutions (Araya et al. 2014) provide an easy-to-use procedure to calculate m-CAK wind hydrodynamics.
Furthermore, it is important to remark that these expressions that represent the hydrodynamics of the wind can be also applied to stellar evolution codes, where mass loss rates are necessary to estimate the evolutionary phases of a star.
In future we plan to consider the stellar rotation into our expressions and, in addition, compare the synthetic line profiles calculated from wind velocity profiles using a hydrodynamic code and our solutions.  Curé et al. (2011) and Venero et al. (2016), respectively. 10). Now from Eqs. A12 and A13 we can solve the singular point location, = . Note that 0.1 to assure a -slow solution (Curé et al. 2011).
Finally, the mass loss rate is solved from the eigenvalue, , when the singular point is replaced in Eq. A2. This paper has been typeset from a T E X/L A T E X file prepared by the author.