Radiation-Driven Wind Hydrodynamics of Massive Stars: A Review

Mass loss from massive stars plays a determining role in their evolution through the upper Hertzsprung-Russell diagram. The hydrodynamic theory that describes their steady-state winds is the line-driven wind theory (m-CAK). From this theory, the mass loss rate and the velocity profile of the wind can be derived, and estimating these properly will have a profound impact on quantitative spectroscopy analyses from the spectra of these objects. Currently, the so-called beta-law, which is an approximation for the fast solution, is widely used instead of m-CAK hydrodynamics, and when the derived value is beta greater than 1.2, there is no hydrodynamic justification for these values. This review focuses on (1) a detailed topological analysis of the equation of motion (EoM), (2) solving the EoM numerically for all three different (fast and two slow) wind solutions, (3) deriving analytical approximations for the velocity profile via the LambertW function and (4) presenting a discussion of the applicability of the slow solutions.


Introduction
At the beginning of the XX century, Johnson [1,2] and Milne [3] argued that the force on ions in the atmosphere of a luminous star could be responsible for the ejection of these ions from the star. They also argued that the ejected ions should carry with them the corresponding number of electrons, strictly there should be no charge-current, but they did not realize at that time that the collisional coupling between ions and protons would drag the rest of the plasma (mostly fully ionized Hydrogen), with them as well, at least to supersonic velocities, and this theory was laid aside. It was Chandrasekhar [4,5], who in the context of globular cluster dynamics, developed the theory of collisions due to an inverse square law, and Spitzer [6] applied Chandrasekhar's theory for collisions between charged particles. Morton [7] was the first to report far-ultraviolet observations of three OB supergiants from an Aerobee-sounding rocket. After this came Copernicus, the first satellite with a telescope on board, and since then it has been possible to obtain stellar spectra in the ultraviolet (UV) region. Morton [7] found that the resonance lines of C IV, N V and Si IV showed the typical P-Cygni 1 profiles. He found that the displacements in the profiles of C IV λλ1549. 5 and Si IV λλ1402. 8 corresponding to outflow velocities in the range 1500-3000 km/s. Snow and Morton [9] showed through a detailed survey that stars brighter than M bol ∼ −6 have strong P-Cygni profiles in their spectra and therefore lose mass. The same conclusion was arrived at by Abbott [10], who compared the radiative force with the gravitational force and concluded that radiative forces could initialize and maintain the mass-loss process for stars with an initial mass at the zero-age main sequence (ZAMS) of about 15 M or greater. the influence of the finite cone angle correction on the dynamics of the wind [described in the appendix from 18]. They found a much better agreement between the improved or modified CAK theory (hereafter m-CAK) and the observations for the mass loss rate and the terminal velocity in a large domain in the Hertzsprung-Russell diagram.
The equation of motion of the m-CAK theory is a highly non-linear differential equation that has singular points, eigenvalues and solution branches [see 16,[23][24][25][26][27]. Since it is challenging to solve this differential equation numerically, PPK found that the velocity field, v(r) from the m-CAK theory, can be described by a simple approximation, known as the β law approximation (see below). In addition, Kudritzki et al. [28, hereafter KPPA] developed analytical approximations for the localization of the critical point, mass loss rate and terminal velocity with an agreement within 5% for v ∞ and 10% forṀ, when compared to the correct numerical calculations.
Radiation-driven stellar winds are hydrodynamic phenomena involving the flow of the outer layers of the atmospheres of massive stars. This review is focused on describing the investigation of the m-CAK hydrodynamic theory, its topology and its three known physical solutions.
Section 2 presents the theory to calculate the radiation (line) force via an analytical description thanks to the Sobolev approximation. Section 3 introduces the m-CAK hydrodynamic theory, and its topological description is given in section 4. Section 5 shows all three known physical solutions, whilst section 6 present analytical approximate solutions based on the LambertW function. Finally, in section 7, we summarise the main topics of this review and discuss the applicability of slow solutions.

The Radiation Force
The exact calculation for the radiation force requires a knowledge of the radiation field (in all the lines and continua) and of the physical processes (scattering, absorption and emission) that contribute to the exchange of energy and momentum throughout the wind. The radiation field is represented by the monochromatic specific intensity I ν (µ), µ is the cosine of the angle between the incoming beam and the velocity vector of the interacting particles. Thus, the radiation force per unit of volume, at a distance r, exerted on a point particle per unit of time is equal to momentum removed from the incident radiation field (κρ I(µ) µ/c) integrated over all the scattering directions. This force is given by where the absorption coefficient κ ν is given in units of cm 2 g −1 . The net flux density comes from the interaction processes, integrated over the whole spectral range, between the radiation field emitted by the photosphere and the stellar wind of mass density ρ at the distance r. Here, it is assumed that the emissivity (thermal emission and photon scattering) in the expanding atmosphere is isotropic. Therefore, no net momentum change occurs from this process (see [29], Chapter 20). The absorption coefficient κ ν consists of three main contributions: where κ Th represents the Thomson scattering, κ cont ν the contribution of bound-free and free-free transitions and κ line ν the sum of all line absorption coefficients at frequency ν. The radiation force can be calculated by state-of-the-art non-local thermodynamic equilibrium (NLTE) radiative transfer codes such as FASTWIND [30,31], CMFGEN [20,[32][33][34] or POWR [35,36], but these calculations depend on the velocity and density profile used to describe the wind.

Radiative force due to electron scattering
The interaction between photons and free electrons is described by a Compton process (an excellent review of this process, including Monte-Carlo calculations can be found in [37]). If photons with energy hν m e c 2 are scattered by Maxwellian electrons 3 having kT m e c 2 , the frequency shift will be very small, but if the scattering process is repeated many times, the small amounts of energy exchanged between the electrons and photons can build up and give rise to substantial effects.
In the non-relativistic limit without the influence of quantum effects (hν m e c 2 ), and neglecting the possible effects described above, the scattering cross-section is frequencyindependent and called the Thomson cross-section, namely: The value of this cross-section is σ Th = 6.65 × 10 −25 cm 2 and the absorption coefficient is, therefore: κ Th ρ = n e σ Th .
Using this value (κ Th ) in Eq. (2) and integrating Eq. (1), we obtain the contribution of the Thomson scattering to the radiation force, where L is the luminosity of the star. The radiative acceleration on the electrons is then It is useful to define the ratio of the Thomson scattering force and the gravitational force by: here G is the gravitational constant and M * the star's mass. In the standard one-component description of stellar winds, the force over the density of the plasma is given by: being ρ = m p n p + ∑ ions (m i n i ) + n e m e the mass density. The principal contribution of the ions comes from Helium, and neglecting the electrons, n e m e , the density is Here A He is the atomic mass of a Helium atom, Y He is the relative abundance by number of Helium with respect to Hydrogen (the latter being described by the subscript p), and m p is the proton mass. Based on the conservation of charge, it is possible to express the electron number density as n e = n p (1 + q He Y He ), where q He = 0, 1, 2 depending on the Helium ionisation state.
Thus, the ratio n e /ρ is: (10) and the acceleration is: or Quite often, the canonical value of κ Th = 0.34 cm 2 g −1 is adopted, which follows from assuming a fully ionised plasma at solar abundance. In addition, since the continuum of OB stars, near its maximum, is also optically thin in the lines, the contribution of the continuum to the total radiative force is neglected. The next section provides a general description of the line force based on the Sobolev approximation (see, e.g., Lamers and Cassinelli [8] or Hubeny and Mihalas [29]).

Radiative force due to lines.
The contribution to the radiation force due to the spectral lines in the wind of massive stars is provided by the momentum transfer of photons (via absorption and re-emission processes in optically thick lines) mainly from the most dominant ions (i.e., C, O, N, and the Fe-group). The proper calculation of the line force (per unit volume) is given by: where φ is the Gaussian absorption profile. The summation is over all the line transitions (l), assuming non-overlapping lines, for which the wind is optically thick. κ l is the opacity coefficient (in cm 2 g −1 ) of lines formed between levels l (lower) and u (upper) with energy h ν 0 , κ l ρ = π e 2 m e c f l n l 1 − n u g l n l g u .
The number density n l and n u of ions in levels l and u are given in cm −3 , g l and g u are the corresponding statistical weights, and f l is the oscillator strength of the line. The CAK theory allows us to find an analytical expression for the line force in a moving media with large velocity gradients in terms of the macroscopic variables using the Sobolev approximation. However, this expression only applies to radiating flows in the non-relativistic regime.

The Sobolev approximation
In a moving plasma like the stellar wind, the interaction of radiation with matter can be better understood as follows. Let's consider a single spectral line thermally broadened with a rest wavelength λ 1 . A photon emitted from the stellar surface, with wavelength λ * < λ 1 , propagates without interacting with the matter until, due to the Doppler shift, it is scattered at the blue edge of the line in question. Due to the expansion of the wind, the particles viewed from any direction from a certain position always appear to be receding. This means that independent of the scattered direction of the photon (forward or backwards), the distance travelled always causes its comoving wavelength to be red-shifted. After many scatterings, the photon's wavelength has been shifted to the line's red edge, and the interaction of this photon with the line (λ 1 ) ceases. The region in the wind where an incoming photon can interact with the ions is called the interaction zone. It is also well known that the winds of massive stars reach terminal velocities of several times the sound speed, and the point at which the wind velocity is equal to the sound speed (the sonic point) is very near to the photosphere. This means that almost all the region where stellar winds are found is supersonic. This description corresponds to the Sobolev approximation [17], where all the relevant physical quantities, such as the opacity, source function, etc., are considered constant in the interaction zone, i.e., the width of the interaction zone is small compared with a characteristic flow length. Thus, for a generic Doppler-broadened line profile, the Sobolev-length, L s , is defined as: where T eff is the star's effective temperature, v th = 2 k B T eff /m p the thermal speed of the protons, k B the Boltzmann constant. A characteristic length of the flow is Typical values of thermal velocities in OB-types stars are about 7 − 20 km/s while terminal velocities are about 1000-3000 km/s, (see, e.g., Lamers and Cassinelli [8], Puls et al. [12]). More recent measurements of terminal velocities based on observations performed in the frame of the ULLYSES collaboration [38] have been accomplished by Hawcroft et al. [39].

The line-force due to a single-line
Castor [18] analysed in detail the Sobolev approximation in the context of stellar winds and showed that the force produced by the incoming radiation due to a single line can be expressed as 4 : where ∆ν d = v th ν/c corresponds to the Doppler shift, F ν is the flux of the radiation field at frequency ν, k l is the monochromatic line absorption coefficient per unit mass, and is the optical depth. Evaluating the optical depth for a normalized Gaussian profile and using the Sobolev approximation, we find: With this expression, we can interpret the RHS of (17) as: i) (F ν ∆ν d /c) is the rate of momentum emitted by the star per unit area at frequency ν with bandwidth ∆ν d , ii) (τ l /k l ) = ρ v th /(dv/dr) represents the amount of mass that can absorb this momentum, and iii) (1 − e −τ l ) is the probability that such an absorption occurs.
Then, by defining, where σ e = σ Th n e /ρ corresponds to the Thomson scattering absorption coefficient per density. In a moving medium, t represents the optical depth that a line will have if its opacity is equal to its electron scattering opacity. Based on this definition, it is possible to rewrite t as where η l = k l /σ e . The first factor in (21) is related only to line properties, and the second only to dynamic variables of the wind.
2.2.3. The line-force due to a statistical distribution of line strength The total line force, due to the addition of all the single lines of the ions, for a point star approximation and for non-overlapping single lines, is given by: Expressing (22) in terms of ∆ν d = ν v th /c, and the relation, F = L/4πr 2 , we obtain Abbott [10] was the first to compile and publish a list of ca. 250 000 lines for atoms from H to Zn in ionisation stages I to VI. Based on such a line list [22,40,41], it is possible to derive a line strength distribution function [24,42]. This distribution can be described as follows: and represents the number of lines in the line-strength interval (k l , k l + ∆k l ) obtained from the total spectrum and weighted by the flux mean of line strength (L ν ν/L). Notice that in Eq. (24) the distribution in frequency space of the lines is independent from the distribution in line strength. An alternative formulation of the line statistic is given by Gayley [43] [see also 22]. The logarithm of the number of lines can be fitted by a linear function, namely: being N 0 the number of lines (strong and weak) that effectively contribute to the line force. Typical values of the parameter α are 0.45 ≤ α ≤ 0.7 [8,42]. Notice that line force parameters are not free but depend on the transfer problem in each individual star [see 22,41,42,[44][45][46], for a detailed description of the calculation of the line-force parameters]. Extending the sum in Eq. (23) to an integral, we obtain the line force expression: Neglecting the lower limit of the integral, a valid approximation for stars of type OB, replacing it by zero and integrating, the line force becomes: where Γ is the Γ-function. Then, dividing by the total density, we obtain the standard form of the line acceleration, with here, Γ e is the radiative acceleration due to Thomson scattering in terms of the gravitational acceleration, andṀ is the mass-loss rate. Here the continuity equation has been used, and the variables, such as N 0 or Γ(α), have been collected into the constant k. Note that this expression for the line-force (Eq. 28) only takes interactions between ions and radially emitted photons into account [16,18].

The correction factor
CAK (see their appendix) discussed qualitatively the effect on the line force that the proper shape of the star (non-radial incoming photons) would have on the wind kinematics. Later PPK and FA independently investigated the influence of this effect, known as the finite disk correction factor, thereby developing the m-CAK theory.
The expression of the line force for incoming photons from an arbitrary direction, for a radial flow velocity field, comes from the definition of Eq. (20), thus, where Inserting t σ in Eq. (26), instead of t, and integrating, we obtain the following expression for the line force: where CF is the correction factor, defined as the ratio of the force due to the non-radial contributions to that of a point star approximation, namely: here µ c = (1 − R 2 * /r 2 ), where R * is the stellar radius, and v = dv/dr. In appendix A, we summarised some properties of the correction factor.

The Ionization Balance
In his work, Abbott [10], assumed local thermodynamic equilibrium (LTE) and used the modified Saha formula (see Hubeny and Mihalas [29]) to take into account the dilution of the radiation field and the possible difference between the electron kinetic temperature T e and the radiation temperature T r . Due to the changes in the ionisation throughout the wind, Abbott fitted the line force not only in terms of (r 2 v v ) α (see Eq. 28) but also as a function of the ratio n e /W(r), where is the dilution factor. He found that the functional dependence of this quotient in the line force is: where the electron number density, n e , is given in units of 10 11 cm −3 . This proportionality means that the greater the density, the lower the ionisation level. In view of the fact that the lower ionisation levels have more line transitions, usually at the maximum of the radiation field, the line force increases with increasing density. Values of this δ line-force parameters for the fast solution (see below) are in the range 0 < δ 0.2 [8], but for a pure Hydrogen atmosphere, the value is δ = 1/3 as Puls et al. [42] demonstrated.

The m-CAK Hydrodynamic model
The 1-D m-CAK stationary model for line-driven stellar winds considers the following assumptions: an isothermal fluid in spherical symmetry and neglecting the influence of viscosity effects, heat conduction and magnetic fields.
The stationary continuity and momentum conservation equations are: being P the fluid pressure, v φ = v rot R * /r, where v rot is the stellar rotational speed at the equator. In addition, g line (ρ, dv/dr, n e ) corresponds to the acceleration due to an ensemble of lines.
The standard or m-CAK parameterization of the line force [10,23,24] is the following: where the coefficient C is given by Eq. (29). Substituting the density from the mass conservation equation (Eq. 36) into the momentum equation (Eq. 37), we obtain the equation of motion (EoM).
Transforming to dimensionless variables, that is: where a is the isothermal sound speed of an ideal gas, P = a 2 ρ.
Using these new variables, the EoM now reads: where the constants are the following: being Z He the number of free electrons provided by Helium, v esc the escape velocity, and the function g is defined as: In order to find a physical wind solution of the EoM (Eq. 42), i.e., starting from the photosphere with a small velocity and reaching infinite with a supersonic velocity, we first need to understand the topology of this equation.

Topological Analysis
As mentioned previously, the first wind model was developed by Parker [14] for the sun. The m-CAK model has a driving force (line force) that depends not only on the radial coordinate r (or u) but also on the velocity and the velocity gradient. These characteristics complicate the study of the EoM's topology that gives rise to the different solutions.
Mathematically, singular points are located where the singularity condition is satisfied, i.e.: and these locations form the locus of singular points. At these specific points, in order to have a smooth wind solution between solution branches, a regularity condition must be imposed, namely: Using the following coordinate transformation: we can now solve Equations (42), (48) and (49), only valid simultaneously at one singular point, obtaining the following set of Equations: derivation details and definitions of f 1 (u, Z), f 2 (u, Z) and f 3 (u, Z) are summarised in appendix B. Solving for Y and C from the set of Equations (52), (53) and (54), we obtain: and These last two Eqs. are generalisations of the relations found by KPPA (see their Eq. [21] and Eq. [34] for Y and Eq. [20] and Eq. [44] for the eigenvalue), but now including the rotational speed of the star.

The critical point function R
The set of Eqs. (52), (53) and (54) are only valid at the singular point for the unknowns: Y s , C s , Z s and u s 5 . Due to the fact that there are only three equations and four unknowns, it is not possible to solve them. Nevertheless, from this set of Equations, we can derive the function, R(u, Z), defined by: where f 123 (u, Z) has the following definition: The locus of singular points, u s , is given by the points which are solutions of the following equation: It should be noted that no approximation has been made in the derivation of the above topological equations.
To determine which is the location of the singular point in the locus of points that satisfies, R(u, Z) = 0, and therefore determine the values of Y s ,C s , Z s and u s , we need to set a boundary condition at the stellar surface. For this boundary condition, the most used are: Set the density at the stellar surface to a specific value, Usually this base density is in the range 10 −8 g cm −3 to 10 −13 g cm −3 . For some examples, see the works of de Araujo and de Freitas Pacheco [47], Friend and MacGregor [48], Madura et al. [49], Curé [50], Araya et al. [51]. ii) Set the optical depth integral to a specific value, i.e., Employing one of these boundary conditions at the stellar surface plus the regularity condition at the singular point, we can solve from the EoM (Eq. 42) the velocity profile, w = w(u), together with the value of the eigenvalue, C , and therefore the mass loss rate,Ṁ.

Types of solutions
We developed a numerical code that discretizes by finite differences the EoM and, using the Newton-Raphson method, iterates to a numerical solution, this code is called HYDWIND and is described in more detail in Curé [50] [see also 52].
After performing the topological analysis of the EoM, we were able, thanks to HYDWIND, to find the numerical solutions of all three known m-CAK physical solutions: fast, Ω-slow and δ-slow solutions.

Fast solution
From the pioneering work of CAK and its improvements from FA and PPK, the code HYDWIND is able to obtain the standard solution of the m-CAK theory, and we called it hereafter the fast solution.
The function R(u, Z) is shown in Fig. 1 for a non-rotating star (a rot = 0). The plane R(u, Z) = 0 is plotted in light grey (right panel) and the intersection of both functions, which corresponds to the locus of singular points, is plotted with black lines. The locus of singular points for the fast solution is the one that starts at Z ∼ 0 and u = −1. The other locus of singular points will be discussed below. Knowing the topology of the m-CAK model, specifically the locus of singular points, we now solve the EoM for the velocity profile, v(r), or equivalently w = w(u), and the eigenvalue C , which is proportional to the mass loss rate,Ṁ. Then, the wind parameters obtained for this model are: a terminal velocity (v ∞ ) of 3467 km/s and a mass-loss rate (Ṁ) of 2.456 × 10 −6 M /yr. Figure 2 shows the velocity profile of this model as a function of log(r/R * − 1) (left panel) and as a function of u (right panel). The location of the singular point (r s ) is shown with a red dot, and it is located near the stellar surface (r s = 1.029 R * or u s = −0.9719). At this point, the wind velocity is 181.4 km/s, a highly supersonic speed (a = 24.17 km/s).  This steep velocity gradient is due to the rapid increase of the line force just above the stellar surface, as shown in Fig. 3, where the sound speed is reached at r = 1.006 R * and the maximum of g line is reached at r = 1.3 R * .  Figure 3. The radiative acceleration, g line , for a typical O5 V star without rotation as function of r/R * , for r < 10 R * .
As previously mentioned, the wind parameters (v ∞ andṀ) must be calculated within the framework of the radiative transport problem. However, to understand the complex non-linear dependence on the wind parameters from the line-force parameters, in the following figures, we show how the wind parameters depend in terms of each one of the line-force parameters. Figure 4 shows the dependence of the wind parameters, v ∞ (left panel) andṀ (right panel) as a function of the line force parameter α, using the same stellar parameters and keeping the line force parameters k and δ fixed. There is an increase in the values of both wind parameters as α increases. The dependence of the wind parameters as a function of the line force parameter k is shown in Fig. 5. In this case, wind parameters also increase as k increases. It is clearly seen in Fig. 5 that the terminal velocity depends only slightly on the value of k rather than the mass-loss rate, which has a significant impact on the value of k [53]. Finally, the dependence of the wind parameters as a function of the line force parameter δ is shown in Fig. 6. We observe that the terminal velocity has a decreasing behaviour when the parameter δ increase, while the mass-loss rate can have a decreasing or increasing behaviour. This behaviour depends on the parameter k, for low values of k the mass-loss rate decreases while δ increases, but for larger values, the behaviour is reversed. In the next subsections, we will discuss each of the slow solutions. The combined effect of the line-force parameters and the three physical wind solutions is discussed in detail in Venero et al. [53].
Overall, the fast solution from the m-CAK theory has been very successful in explaining the terminal velocities and mass-loss rates of massive stars [see 8,11-13].

Ω-slow solution
The original CAK model considered the star point approximation, i.e., all the photons are radially directed over the wind plasma. In that work, CAK only discussed the effect of the finite disk of the star seen by an observer in the wind. FA and PPK implemented the finite disk correction factor and solved the EoM. In both works, they also studied the influence of rotation in the equatorial plane of a rotating star, but they could not obtain solutions for rapidly rotating stars. The reason was found by Curé [50], for Ω 0.75, where Ω = v rot /v crit . From this value of Ω, the fast solution ceases to exist, and another type of solution is found. This solution called the Ω-slow solution, is characterised by a slower and denser wind in comparison with the fast solution.
It is well known that Be stars are the fastest rotators among stars [54]. Thus, in this section, we will study the topology and the wind solutions for a typical B2.5 V star with the following stellar and line force parameters: T eff = 20000K, log g = 4.11, R/R = 4.0, k = 0.61, α = 0.5, and δ = 0.0. The lower (surface) boundary condition is fixed at ρ * = 8.7 × 10 −13 g/cm 3 . In addition, the distortion of the shape of the star caused by its high rotational speed and gravity-darkening effects are not considered.     only Ω-slow solutions; the location of the singular point is almost independent of Ω. This is a consequence of the shape of the locus curve of singular points (see Fig. 7). This locus is located almost at a constant value of u.
The Ω-slow solutions are only valid in the equatorial plane in this 1D m-CAK model. Notice that this model does not take into account the oblateness and gravity-darkening effects. See Araya et al. [55] for the implementation in the 1D model (equatorial plane) and Cranmer and Owocki [56] for the implementation in the 2D model. In the equatorial plane, the higher Ω, the greater the centrifugal force and, consequently, the lesser the effective gravity. Therefore, the higher Ω is, the higher the rate of mass loss and, through the continuity equation, the higher the wind density, as shown in Fig. 9.  Figure 9. Wind density profiles, ρ(u) (in gr/cm 3 ) versus u, for fast and Ω-slow solutions. The colour scheme is the same as the one used in Fig. 8

δ-slow solution
The δ-slow solution was found numerically by Curé et al. [27]. This solution, based on the m-CAK theory, describes the wind velocity profile when the ionization-related line-force parameter δ takes larger values, δ 0.28. These values are larger than the ones provided by the standard m-CAK solution [see, 8, and references therein]. Nevertheless, Puls et al. [42] calculated the value of δ for a pure Hydrogen atmosphere finding a value of δ ∼ 1/3. These high values of δ are also found in atmospheres and winds with extremely low metallicities [see, 57].
The δ-slow solution, as well as the Ω-slow solution, is characterized by low velocities. This solution could explain the velocities obtained for late-B and A-type supergiant stars and seems to fit well the observed anomalous correlation between the terminal and escape velocities found in A supergiant stars [27]. To present the topological analysis of this type of solution, we adopt the model T19 from Venero et al. [53]. The stellar and line-force parameters are T eff = 19000K, log g = 2.5, R/R = 40, k = 0.32, and α = 0.5. We use τ * = 2/3 as a boundary condition at the stellar surface. In

The β-law approximation
In the work of PPK, after obtaining the numerical solution of the EoM, they assumed a power law approximation to describe the velocity profile only as a function of the radial coordinate r. This approximation is known as the β-law approximation and has the following expression: where v ∞ is the terminal velocity and the value of β determine the shape of the velocity profile.
In the context of stellar wind diagnostics, these parameters are considered fit parameters that must be determined through spectral line fitting. Usually, the range used for the β parameter is 0.7 β 4 [11]. Figure 12 shows the velocity profile of the fast solution for the stellar and line-force parameters given at the beginning of section 5.1, together with six different values of β for a β-law velocity profile. From this figure, we can conclude that the fast solution cannot be described properly by a β-law with β > 1.2.
On the other hand, Figure 13 shows the velocity profile for the δ-slow solution for the stellar and line-force parameters given at the beginning of section 5.3, and δ = 0.32. Also, the same β-law profiles of Fig. 12

Analytical wind solutions
Using an analytical expression to represent the radiation force and solve the equation of motion analytically offers numerous advantages over the numerical integration of the EoM. These formulae can be used in all cases where good first estimates are needed; for example, it gives the advantage of solving the radiative transfer problem for moving media in an easier way.
KPPA were one of the first to develop analytical solutions for radiation-driven winds considering the finite disk correction factor in the line force. Based on these solutions, they provided approximated analytical expressions for the terminal velocity and the mass loss rate in terms of the stellar parameters (L, M * and R * ), the line force parameters (k, α and δ) and the free parameter β from the β-law (they adopted β = 1.0 for their results).
Other authors have tried to simplify the complicated numerical treatment from the theory. Villata [58] (hereafter V92), with the purpose of simplifying the integration of the EoM, derived an approximated expression for the line acceleration term, which depends only on the radial coordinate. Müller and Vink [59] (hereafter MV08) presented an analytical expression for the velocity field using a parameterized description for the line acceleration that (as in V92's) also depends on the radial coordinate. These line-acceleration expressions do not depend on the velocity or the velocity gradient, as the standard m-CAK description does. Araya et al. [60] proposed to achieve a complete analytical description of the 1-D hydrodynamical solution for radiation-driven winds in the fast regime by gathering the advantages of both previous approximations (the use of known parameters and the LambertW function). In addition, Araya et al. [61] developed an analytical solution for the δ-slow regime. To date, no approximation using the LambertW function has been performed for the Ω-slow regime, we expect to do this in the future.
In the next sections, we describe the results and methodology used to solve analytically the equation of motion for fast and δ-slow regimes.

Solution of the dimensionless equation of motion
In this section, we recapitulated the methodology described by MV08 to obtain the dimensionless equation of motion.
In a dimensionless form, the momentum equation can be expressed as follows, wherer is a dimensionless radial coordinater = r/R * , and the dimensionless velocities (in units of the isothermal sound speed a) are: here, v crit is the rotational break-up velocity in the case of a rotating star. It is usually determined by dividing the effective escape velocity, v esc , by a factor of √ 2. Similarly, a dimensionless line acceleration can be written as follows:ĝ Using the continuity equation and the equation of state of an ideal gas, the dimensionless equation of motion reads as follows: Lastly, a 1-D velocity profile is derived analytically in terms of LambertW function [62][63][64]. See MV08 for a detailed description of the methodology used to arrive at this solution. This analytical solution is expressed as follows: where In the last equation appears the parameterr c , which represents the position of the sonic (or critical) point. Also, the LambertW function (W j ) has only two real branches, indicated by the sub-index j, where j = 0, −1. These two branches coincide at the sonic point,r c , i.e., j = 0 for 1 ≤r ≤r c −1 forr >r c , A regularity condition must be imposed, as in the m-CAK case, since the LHS of the equation of motion (Eq. 67) vanishes atv = 1 (singularity condition in CAK formalism). This is equivalent to ensuring that the RHS of Eq. (67) vanishes atr =r c . Therefore, andr c is obtained by solving this last equation. Finally, the velocity profile is derived using the function x(r), Eq. (69), into Eq. (68).

The fast regime
KPPA's analytical study of radiation-driven stellar winds allowed V92 to derive an approximate expression for the line acceleration term. In this case, the line acceleration is only dependent on the radial coordinate, and it reads as follows: with According to Kudritzki et al. [65], the exponent β can be calculated based on the force multiplier parameters and the escape velocity, v esc : with v esc in km/s. Then, using Eq. (72) in its dimensionless form (Eq. 66) and inserting it into the dimensionless equation of motion (Eq. 67), it yields: Based on V92's approximation of the line acceleration, this differential equation can be viewed as a solar-like differential equation of motion. Hence, the singular point is the sonic point. Additionally, V92's equation of motion does not have eigenvalues, which means it doesn't depend explicitly on the star's mass loss rate.
Using a standard numerical integration method, V92 solved the equation of motion and obtained terminal velocities that were within 3-4 % of those computed by PPK and KPP.
A parametrized description of line acceleration was presented years later by MV08 that is dependent on the radial coordinate (like V92's). The line acceleration in MV08 was determined using Monte-Carlo multi-line radiative transfer calculations [66,67] and a β law. Following this, the line acceleration was fitted using the following formula: whereĝ 0 , δ 1 ,r 0 and γ are the acceleration line parameters. Then, the solution of the equation of motion, based on their methodology and line acceleration, is:v with As a result of the approximations described above, the velocity profile can be represented analytically, greatly simplifying the solution of the equation of motion.
Also, it is relevant to note that each of the mentioned approximations has its own advantages and disadvantages. Even though Villata's approximation of the radiation force is general and can directly be applied to describe any massive star's wind, the momentum equation still needs to be solved numerically. With MV08's approximation, the equation of motion can be analytically solved, based onĝ 0 , δ 1 ,r 0 , and γ parameters of the star. Nevertheless, it is still necessary to perform Monte Carlo multi-line radiative transfer calculations in order to determine these parameters.
This methodology to solve the equation of motion was revisited by Araya et al. [60], and in order to derive a fully analytical expression combining V92's expression of the equation of motion, with the methodology developed by MV08. This analytical solution is,v with where As was mentioned in the previous section,r c can be obtained numerically, making the RHS of Eq. (75) equal zero. In order to obtain the terminal velocity in a simpler way, we can use the average value ofr c (r c = 1.0026) obtained by Araya et al. [60]. Note that this value can be used only in the supersonic region. Equation (79) has the advantage that it is based not only on the LambertW function but also on stellar parameters and the line force parameters. For a wide range of spectral types, stellar and force multiplier parameters are given [see, e.g., 8,10,22,24,41,44].
By comparing the analytical solution to the 1-D hydrodynamic code HYDWIND, the accuracy of the analytical solution can be tested. Figure 15 compares the results obtained with our analytical approximation to those obtained with the hydrodynamics for four stars taken from Araya et al. [60]. Both solutions have similar behaviours. However, as shown by Araya et al. [60], the analytical approximation close to the stellar surface (subsonic region) is not good enough. In the same way, Figure 16 compares the numerical and analytical velocity profiles near to the stellar surface for Ori.  There is a limitation to this analytical expression when the line force parameter δ exceeds about 0.3. This is due to the complexity of a term in the proposed line acceleration expression.
To obtain an expression with real values, high values of δ would require high values of α. However, such kind of α values would be totally unphysical (α > 1). As an illustration of the dependence of this expression on the parameters α and δ, Fig. 17 shows the domain of the complex and real regions when this expression is evaluated to the given line acceleration term.

The δ-slow regime
Considering the results obtained when using an approximate description of the wind velocity for the δ-slow case, Araya et al. [61] modified the function of the line acceleration given by MV08 to better describe of the δ-slow wind.
As a result, the proposed line acceleration is: whereĝ 0 , δ 1 , δ 2 , and γ are the new set of line acceleration parameters.
It is notable that the δ 2 parameter, which has been incorporated into this new expression, provides a much better agreement with the numerical line acceleration obtained from the m-CAK model in the δ-slow regime compared with the one from MV08.
Based on this new definition of the radiation force, the new dimensionless equation of motion reads: The LambertW function is used to solve the equation of motion, Eq. (83) following the same methodology developed by MV08,v with being 2 F 1 the Gauss Hypergeometric function. The critical (or sonic) point,r c , is obtained numerically, making the RHS of Eq. (83) equal zero.
Ultimately, this expression for the velocity profile is in quite satisfactory agreement with the numerical solution from HYDWIND.
As described in Araya et al. [60], a relationship was established between the MV08 lineforce parameters (ĝ 0 , δ 1 ,r 0 , and γ) and the stellar and m-CAK line-force parameters. In addition to being easy to use, this relationship provides a straightforward and versatile method of calculating velocity profiles analytically for a wide range of spectral types since both stellar and m-CAK line force parameters are available [see, 8,10,22,24,41,44].
A similar relationship can be derived for the δ-slow regime using m-CAK hydrodynamic models, that is, creating a grid of HYDWIND models for δ-slow solutions. These models are then analysed using a multivariate multiple regression analysis [ MMR 68,69].
To develop this hydrodynamic grid, the stellar radius is calculated from M bol for each pair of stellar parameters (T eff , log g) by using the flux weighted gravity-luminosity relationship [70,71]. Additionally, a total of 20 stellar radius values were added (ranging from 5 R to 100 R in steps of 5 R ). Surface gravities are in the range of log g = 2.7 down to about 90% of Eddington's limit in steps of 0.15 dex. Effective temperatures are between 9 000 K to 19 500 K, in steps of 500 K. The range of this grid has been chosen to cover the region of the T eff -log g diagram that contains B-and A-type supergiants. In Table 1, the m-CAK line-force parameters for each set of (T eff , log g) values are listed. For the purpose of obtaining δ-slow solutions, only high values of δ are considered. For the T eff -log g plane (see Fig. 18), we show in blue dots all converged models. log g Figure 18. Hydrodynamic models in the T eff -log g plane. Blue dots represent the converged solutions. Grey solid lines are the evolutionary tracks for stars of 7M to 60M without rotation [72], and black lines represent the zero-age main sequence (ZAMS) and the terminal age main-sequence (TAMS). In order to obtain the new line acceleration parameters (ĝ 0 , δ 1 , δ 2 , and γ) for each model, the m-CAK line acceleration was fitted, using Least Squares, with the proposed line acceleration expression (Eq. 82). Then, a MMR is applied to the grid of models in order to derive the relationship between new line acceleration parameters (ĝ 0 , δ 1 , δ 2 and γ) and stellar (T eff , log g, R * /R ) and m-CAK line-force parameters (k, α, δ). The estimated parameters are:  After the estimated values for each dependent variable,ĝ 0.27 0 , (δ 1 + 1) 5.3 , δ 0.45 2 , (γ + 1) −3.56 , are obtained they are transformed intoĝ 0 , δ 1 , δ 2 , and γ through their respective inverse functions. Finally, we can use these parameters in Eq. (84) to calculate the velocity profile.
The velocity profiles obtained via HYDWIND code and the analytical solution are shown in Fig. 19 for one model with δ-slow solution. The model is taken from Curé et al. [27]. Remember, however, that this relationship holds only for δ-slow solutions, especially for values of δ between 0.29 and 0.35. In addition, considering the number of converged models in the grid, the authors recommend using this expression for values of α between 0.45 and 0.55.
Finally, it is important to remark that both analytical solutions for the velocity profile, fast and δ-slow, do depend only on the stellar (T eff , log g, R * /R ) and m-CAK line-force parameters (k, α, δ). Regarding to the mass-loss rate, V92 proposed an expression to obtain it (in terms of the stellar and m-CAK parameters) following the approximations given by Kudritzki et al. [28]. Also, Araya et al. [61], in the appendix, proposed a method to obtain the mass-loss based on Curé [50].

Summary and discussion
Observations over the past decades have shown that the basic wind parameters behave as predicted by the theory. This fundamental agreement between observations and theory provides strong evidence that the winds from massive stars are driven by radiation pressure. This has given the m-CAK theory a well-established status in the massive star community. However, several issues are contentious and still unclear, such as the calibration of the Wind-momentum Luminosity Relationship (WLR) [73], disks of Be stars, wind parameters determination, the applicability of the called slow wind solutions, among others. All these issues are the focus of massive stars research.
This review is focused on the theoretical and numerical research of wind hydrodynamics of massive stars based on the m-CAK theory, with particular emphasis on its topology and hydrodynamic solutions.
We presented a topological analysis of the one-dimensional m-CAK hydrodynamic model and its three known hydrodynamic solutions, the fast, Ω-slow and δ-slow solutions. From a topological point of view, slow solutions are obtained from a new branch of solutions with a locus of singular points far from the stellar surface, unlike fast solutions with a family of singular points near the stellar surface.
We continued analyzing the dependence of the line force parameters (k, α, and δ) with the wind parameters (mass loss rate and terminal velocity) in order to understand the complex non-linear dependence between these parameters. In the case of α, there is an increase in both wind parameters as this parameter increases. This behaviour is similar to the k parameter, but the dependence is very slight for terminal velocity, while the mass loss rate has a significant impact. For the δ parameter, the terminal velocity has a decreasing behaviour when this parameter increases, while the mass-loss rate can have a decreasing or increasing behaviour, which depends on the parameter k. When k is low, mass-loss rates decrease, while δ increases, whereas when k is large, the opposite occurs.
In addition, we compared the β-law with the hydrodynamic solutions. We concluded that the fast solution could not be adequately described by a β-law with β > 1.2, while the δ-slow solution cannot be described by any β-law.
Furthermore, we presented two analytical expressions for the solution of line-driven winds in terms of the stellar and line-force parameters. The expressions are addressed to obtain the fast and δ-slow velocity regimes in a simple way. Both solutions are based on the LambertW function and an approximative expression for the wind line acceleration as a function of the radial distance. The importance of an analytical solution lies in its simplicity in studying the properties of the wind instead of solving complex hydrodynamic equations. In addition, these analytical expressions can be used in radiative transfer or stellar evolution codes [see, e.g., 74].
Concerning the applicability of the slows solutions, in the case of the Ω-slow solution, their behaviour suggest that it can play a paramount role in the ejection of material to the equatorial circumstellar environment of Be and B[e] stars. Be stars are a unique set of massive stars whose main distinguishing characteristics are rapid rotation and the presence of a dense, gaseous circumstellar disk orbiting in a quasi-Keplerian fashion. There is a long-standing problem in understanding the formation of disks in Be stars; this is one of the major areas of ongoing research in Be stars. The gaseous disks are not remnants of the objects' protostellar environments, nor are they formed through the accretion of material [54,75]. On the contrary, the equatorial gas consists of a decretion disk formed from a material originating from the central star.
As was stated above, attempts to solve this problem have been made without much success, for example, the link between the line-driven winds and these discs, called the Wind Compressed Disk [76], but the work of Owocki et al. [77] was the first to show this is not a viable mechanism for rapidly rotating stars due to the non-radial line force components. The most accepted model to successfully reproduce many Be star observables is the viscous decretion disk (VDD) model, developed by Lee et al. [78] and examined by Okazaki [79], Bjorkman and Carciofi [25], Krtička et al. [80] and Curé et al. [81]. Currently, how that material is ejected into the equatorial plane and how sufficient angular momentum is transferred to the disk to maintain quasi-Keplerian rotation are among the primary unresolved questions currently driving classical Be star research.
Araya et al. [55] studied the Ω-slow wind solution and its relation with the disks of Be stars. Overall, this work is an extension to the study done by Silaj et al. [82], where they precisely investigated if the density distribution provided by the Ω-slow wind solution could adequately describe the physical conditions to form a dense disc in Keplerian rotation via angular momentum transfer. They considered a two-component wind model, i.e., a fast, thin wind in the polar latitudes and a Ω-slow, dense wind in the equatorial regions. Based on the equatorial density distributions, Hα line profiles are generated and compared with an ad-hoc emission profile, which agrees with the observations. In addition, their calculations assumed three different scenarios related to the shape (oblate correction factor) and the star's brightness (gravity darkening). Finally, they found that under certain conditions (related to the line-force parameter of the wind), a significant Hα line profile could be produced, and maybe the line-driven winds through the Ω-slow solution can have an essential role in the disc formation of Be stars.
In addition, Araya et al. [51] studied the zone where the classical m-CAK fast solution ceases to exist, and the Ω-slow solution emerges at rapid rotational speeds. This study used two hydrodynamic codes with time-dependent and stationary approaches. They found that both solutions can co-exist in this transition region, which depends exclusively on the initial conditions. In addition, they performed base-density wind perturbations to test the stability of the solution within the co-existence region. A switch of the solution was found under certain perturbation conditions. The results are addressed to a possible role in the ejection of extra material into the equatorial plane via pulsation modes, where the Ω-slow solution can play an important role.
A current weakness of this m-CAK model is that it does not consider the role of viscosity and its influence on angular momentum transport. This mechanism might explain the formation of a Keplerian disc.
On the other hand, the δ-slow solution is promising to explain the discrepancies of the wind parameters between observations and theory in late B-and A-type supergiant stars. According to the findings of Venero et al. [53], these suggest that the terminal velocity of early and mid-B supergiants agrees with the results seen from fast outflowing winds. By contrast, the results obtained for late B supergiants and, mainly, those obtained for early A supergiants agree with the results achieved for δ-slow stationary outflowing wind regimes. Then, the δ-slow solution enables us to describe the global features of the wind quite well, such as mass-loss rates and terminal velocities of moderately and slowly rotating B supergiants.
Conversely, Venero et al. [53] stated that the δ-slow solution seems not to be present in stars with T eff > 21 000 K. This restricts the possibility of a switch between fast and slow regimes at such temperatures. Consequently, this would be a physical explanation for why an empirical bi-stability jump can be observed around 21 000 K in B supergiants [83]. From a theoretical perspective, a velocity jump has also been found using Monte Carlo modelling and the co-moving frame method [see, e.g., 13,84,85] In addition, it is generally accepted that most O and early B-type stars can be adequately modelled with a β velocity law with 0.7 β 1.2. However, supergiants A and B exhibit β values that tend towards higher values, often β ≥ 2 [see, e.g., [86][87][88][89][90]. Venero et al. 2023 (in preparation) propose that δ-slow solutions might explain these winds. They show that the δ-slow regime could adequately fit the Hα line profile of B supergiants with the same accuracy as that obtained using a β-law with β ≥ 2, but now with a hydrodynamic explanation of the velocity profile used.
The investigation carried out in the latest works inspired us to go deeper into the possible role of slow wind solutions with respect to the unresolved questions related to massive stars. In view of our results, we are encouraged to develop this line of research further. In the case of the Ω-slow solution and its link to Be stars, or possibly to B[e] stars, it would require 2D/3D models for a better understanding to take into account non-radial forces, the effects of stellar distortion and gravity darkening. These considerations could change, in turn, the nature of the Ω-slow solution or the behaviour regarding the co-existence of solutions and a switch between them. The δ-slow solution could play an essential role in understanding the winds of B-and Atype supergiants. Moreover, this solution is expected to solve the disagreement between the observations and theory for these stars and, in this way, calibrate the wind-momentum luminosity relationship.
As we mentioned previously. In the standard procedure for finding stellar and wind parameters, the β-law (β, v ∞ ) and a mass-loss rate (Ṁ) are three 'free' input parameters in radiative transfer codes, comparing synthetic spectra with the observed spectra of a star. The β law comes from an approximation of the fast wind solutions, and the values of β should be in a restricted interval. To improve the hydrodynamic approximation used in this standard procedure, we have developed two hydrodynamic procedures to derive stellar and wind parameters: • The self-consistent CAK procedure [44], based on the m-CAK model. Here we iteratively calculate the line-force parameters using the atomic line database from CMFGEN, coupled with the m-CAK hydrodynamic until convergence. We obtain the line-force parameters and, therefore, the velocity profile and the mass-loss rate. Thus, none of the input parameters are 'free', but self-consistently calculated. • The LambertW procedure [45]. In this procedure, we start using a β-law and a value foṙ M in CMFGEN. After convergence, we calculated the line acceleration as a function of r, and using the LambertW function we obtain a new velocity profile. This is inserted in CMFGEN and the cycle is repeated until convergence. In this LambertW procedure, the only input parameter is the mass-loss rate.
We expect that these two alternatives, which reduce the number of input parameters, will in the future, have a significant impact on the standard procedures for obtaining stellar and wind parameters of massive stars.