© M. Curé et al. 2022

1 Introduction

Classical Be stars (CBes) are fast rotating main-sequence B-type stars that form an equatorial gas rotating disk. The inner part of the disk, usually ≲20 stellar radii, is geometrically thin and optically thick and rotates in a quasi-Keplerian orbit (Quirrenbach et al. 1997; Meilland et al. 2012; Rivinius et al. 2013). These disks, referred to as decretion disks, are built from mass ejected from the equatorial stellar surface that acquires sufficient velocity and angular momentum to orbit the star. Once the material is ejected from the star, it is generally governed by gravity and viscous forces (for the latest review, see Rivinius et al. 2013, and references therein).

Currently, the best theoretical framework describing the evolution of these disks once formed is the viscous decretion disk model (VDD) proposed by Lee et al. (1991). In this model, viscosity acts to shuffle the angular momentum of the circumstellar material, and with the assumption of Keplerian gas motion, a steady and thermally stable structure is obtained. The strength of viscosity is parameterized as α, according to the α-disk model proposed by Shakura & Sunyaev (1973) for accretion disks. In the case of decretion disks, the mass loss rate in charge of feeding the disk has an opposite sign. The α parameter dictates the timescales over which Be star decretion disks evolve and dissipate. Haubois et al. (2012) used the VDD model to study the temporal evolution of the disk density in Be stars and found that a small fraction (~1%) of the mass ejected by the star acquires sufficient angular momentum to move outward, orbiting at increasing radii, while most of the ejected mass falls back onto the star. Rímulo et al. (2018) analyzed 81 outburst events in 54 CBes by modeling light curves with models that solve the radiative transfer problem in the nonlocal thermodynamic equilibrium regime in two or three dimensions. They found that the viscosity parameter was higher during disk build-up (α = 0.63) than during the dissipation phase (α = 0.29). This means that the formation timescale is faster than the dissipation phase. In addition to viscosity, radiative ablation may systematically remove, by means of radiative acceleration, some material from the inner part of CBes disks, thus forming a disk-wind type of structure both above and below the disk (Kee et al. 2016, 2018a,b; Kee & Kuiper 2019).

On the other hand, Okazaki (2001) solved the hydrodynamic equations for a viscous disk through the α-disk model, assuming isothermal conditions and a radiative force given by an ensemble of optically thin lines. His goal was to find the transition between near-Keplerian and angular momentum conserving motion inside the disk. Since the specific angular momentum increases with the distance from the stellar photosphere making it difficult to hold Keplerian motion at larger distances. He found for three different α values (1.0, 0.1, and 0.01) the existence of a transonic solution. The sonic point was located far from the stellar surface at values larger than 100 stellar radii. Inside this region, the outflow was found to be highly subsonic. For a particular set of line-force parameters Okazaki (2001) found that for α ≥ 0.95 the topology of the sonic point is nodal, and for α < 0.9 it is saddle. Finally, he found that in the inner subsonic region the disk is near-Keplerian, while in the outer subsonic and supersonic regions the angular momentum is conserved. We note that in his work the specific angular momentum at the sonic point ℓ_s is set as an eigenvalue, making the angular velocity a fixed quantity.

Krtička et al. (2011) studied the mass-loss rate occurring inside a decretion disk that stems from the angular momentum loss due to the condition of maintaining the star rotating close to and below its critical rotational speed. Then, from a specific angular momentum, the mass-loss rate is related with the outer disk size. In their work they discussed three different physical processes that affect the outer disk: binarity, thermal expansion in supersonic flows, and radiative ablation in the inner disk, and they also presented how to implement these considerations about the mass-loss rates in decretion disks in stellar evolution codes. Another study of the outer disk in Be stars was the work by Kurfürst et al. (2014), who studied the dependency between the physical features of large disks and their temperature and viscosity distribution. In spite of the relevance of the results of these two works, in their calculations they do not include radiation force in the wind equations.

Klement et al. (2017) noted the importance of extending the study of the density structure using VDD in the outer parts of the disk. Since most observational constraints come from optical and IR wavelengths, the outer regions are detectable only at radio wavelengths. In their work, they compiled data from ultraviolet up to radio wavelengths for the spectral energy distribution and found a turn-down in the SED at 10–60 μm in a sample of six Be stars. To reproduce the observations considering a VDD model, they used a truncated disk with radii between 26 and 108 stellar radii. They concluded that tidal forces from a binary companion are the only mechanism (at close distances to the star) that can truncate the disk.

In this paper we revisit the work of Okazaki (2001), but consider another solution schema that does not involve angular variables, such as an eigenvalue, and which will finally lead us to conclusions that can respond to more recent works about observational characteristics in the flux distribution of Be stars, and the effects of rotational speed and the velocity and density wind structure.

This work is organized as follows. Section 2 reviews the VDD model and proposes a new procedure to solve the equation of motion. In Sect. 3 we analyze the topology at the singularity from the equation of motion. Then in Sect. 4 a numerical procedure is proposed to obtain the different viscous transonic decretion disk solutions. In Sect. 5 numerical calculations are performed to understand the influence of the different parameters. Finally, in Sects. 6 and 7 we give a discussion and our conclusions, respectively.

2 Viscous Decretion Disks Hydrodynamic Model

Based on the VDD scenario proposed by Lee et al. (1991), Okazaki (2001) studied the influence of a radiative force due to an ensemble of optically thin lines described by the ad hoc model from Chen & Marlborough (1994). In addition, the Shakura–Sunyaev prescription for the viscous stress was adopted. For the sake of consistency, we follow the same nomenclature and equations from Okazaki (2001, see derivation details therein). The geometrically thin circumstellar disk of a Be star is in steady state, and is symmetric about the rotational axis and the equatorial plane.

The disk equations in cylindrical coordinates (r, ϕ, z) follow mass conservation $$\dot{M}+2\pi r{V}_{\text{r}}\Sigma =0,$$(1)

where Ṁ is the mass loss rate, Σ is the vertical integrated density, and V_r is the radial component of the vertical averaged velocity. The momentum conservation, radial and angular components are, respectively, $$-V_r\frac{\text{d}{V}_r}{\text{d}r}+\frac{V_{\varphi }^2}{r}-\frac{GM}{r^2}-\frac{1}{\Sigma }\frac{\text{d\hspace{0.17em}}W}{\text{d}r}+g_{\text{rad}}+\frac{3W}{2r\Sigma }=0,$$(2) $$V_r\frac{\text{d}{V}_{\varphi }}{\text{d}r}+\frac{V_r{V}_{\varphi }}{r}-\frac{1}{r^2\Sigma }\frac{\text{d}\left(r^2{t}_{r\varphi }\right)}{\text{d}r}=0.$$(3)

Here V_ϕ is the angular component of the vertical (z-axis) averaged velocity, M is the stellar mass, G is the gravitational constant, W is the pressure, and t_rϕ is the r − ϕ component of the viscous stress. The radiative force (vertically averaged) is given by F_rad (see below). The state equation of an ideal gas reads $$W=c_{\text{s}}^2\Sigma ,$$(4)

where $c_{\text{s}}^2$ is the isothermal sound speed. Finally, the Shakura–Sunyaev viscosity prescription is given in terms of the r − ϕ component of the viscous stress $$t_{\text{r}\varphi }=-\alpha W,$$(5)

where α is the viscosity parameter.

2.1 Radiative Acceleration

The radiative acceleration description used in the VDD model was proposed by Chen & Marlborough (1994). The radiative acceleration is produced by an ensemble of optically thin lines $$g_{\text{rad}}=\frac{GM\Gamma }{r^2}+\frac{GM\left(1-\Gamma \right)}{r^2}\eta {\left(\frac{r}{R}\right)}^{ϵ},$$(6)

where Г is the Eddington factor due to electron scattering. The parameters ϵ and η characterize the decay rate and magnitude of the radiative (line) force, respectively, and R is the stellar radius.

2.2 Okazaki’s Solution Procedure

Rearranging Eqs. (1)-(5), a constant of motion is found: $$\ell +\alpha c_{\text{s}}^2\frac{r}{V_{\text{r}}}=C.$$(7)

Here ℓ is the angular momentum per mass unit (specific angular momentum) and C is an unknown constant of motion or eigenvalue. The solution of this problem relies on a suitable form for calculating the constant C. Okazaki (2001) adopted the evaluation of C at the sonic point (r = r_s) $$C={\ell }_{\text{s}}+\alpha c_{\text{s}}^2\frac{r_{\text{s}}}{c_{\text{s}}},$$(8)

where ℓ_s is the specific angular momentum at the sonic point. Now the eigenvalue of this problem is ℓ_s. Then, solving Eq. (7) for ℓ and eliminating t_rϕ, W, and Σ from Eqs. (1)-(3) we obtain Okazaki’s equation of motion (hereafter OEoM) $$\left(1-\frac{c_{\text{s}}^2}{V_{\text{r}}^2}\right)V_{\text{r}}\frac{\text{d}{V}_{\text{r}}}{\text{d}r}=\frac{{\ell }^2}{r^3}+\frac{5}{2}\frac{c_{\text{s}}^2}{r}+g_{\text{eff}},$$(9)

where ℓ is given by $$\ell ={\ell }_{\text{s}}+\alpha c_{\text{s}}^2\left(\frac{r_{\text{s}}}{c_{\text{s}}}-\frac{r}{V_{\text{r}}}\right),$$(10)

and the effective gravity g_eff reads $$\begin{array}{ll}{g}_{\text{eff}}\hfill & =-\frac{GM}{r^2}+g_{\text{rad}}\hfill \\ \hfill & =-\frac{GM\left(1-\Gamma \right)}{r^2}\left(1-\eta {\left(\frac{r}{R}\right)}^{ϵ}\right).\hfill \end{array}$$(11)

His work was advocated to solve the set of Eqs. (9)-(11).

Our main criticism to this approach is the evaluation of C at the sonic point for ℓ_s. This means that now ℓ_s is the eigenvalue of the OEoM and corresponds to a rotational quantity (i.e., it fixes the value of the stellar rotational velocity). Thus, for an individual star with a set of stellar and line-force parameters we obtain a solution of the OEoM only for one specific value for the stellar rotational speed.

At the stellar surface r = R, the azimuthal component of the vertically averaged velocity is equal to the stellar rotational velocity, V_ϕ(r = R) = V_rot ≡ ΩV_crit. Here V_crit is the critical rotational speed: $$V_{\text{crit}}=\sqrt{\frac{GM\left(1-\Gamma \right)}{R}}.$$(12)

Thus, the variable $V_{\varphi }/\sqrt{\left(GM/R\right)}$ at the stellar surface is $$\frac{V_{\varphi }\left(R\right)}{\sqrt{\left(GM/R\right)}}=\Omega \sqrt{\left(1-\Gamma \right)}.$$(13)

For classical Be stars, the Eddington factor Г is a very small value and can be neglected, obtaining $$\frac{V_{\varphi }\left(R\right)}{\sqrt{\left(GM/R\right)}}\simeq \Omega .$$(14)

A close inspection of Fig. 2b from Okazaki (2001), specifically the upper left corner for the variable $V_{\varphi }\left(R\right)/\sqrt{\left(GM/R\right)}$, shows that for different values of the line-force parameters, the rotational speed of the star is also different. In other words, the eigenvalue of the problem ℓ_s, which is directly related with the stellar rotational quantity Ω, depends on the line-force parameters η and ϵ. This solution scheme is completely different from the one used in m-CAK theory (see, e.g., Curé 2004, and references therein), where Ω is one of the stellar parameters and it is not determined from the solution of the OEoM (Eq. (9)).

2.3 A New Solution Procedure

In this section we propose an alternative procedure to solve the equation of motion, where the stellar rotational speed is an input stellar parameter and not an eigenvalue. In addition, in Okazakis’s procedure, the constant C has to be evaluated at the sonic point, at this point the value of the radial velocity is known, but not its location; therefore, there are two unknown quantities, ℓ_s and r_s, that determine the value of C.

2.3.1 Evaluation of the Constant of Motion C

Then evaluating C (Eq. (7)) at the stellar surface r = R, we get $$C=\ell \left(R\right)+\alpha c_{\text{s}}^2\frac{R}{V_{\text{r}}\left(R\right)},$$(15)

where ℓ(R) = ℓ_R = R V_ϕ = R Ω V_crit is a known quantity. Thus, ℓ has now the following expression: $$\ell ={\ell }_{\text{R}}+\alpha c_{\text{s}}^2\left(\frac{R}{V_{\text{r}}}-\frac{r}{V_{\text{r}}}\right).$$(16)

In this new solution schema, V_r = V_r(R) is the eigenvalue.

In order obtain a solution to this problem, we have to solve Eq. (9), together with Eqs. (11) and (16). Here we do not have a typical eigenvalue problem for the unknown V_r, as in the m-CAK theory. In this case, V_r is also the initial condition at the stellar surface of this nonlinear fist-order differential equation.

2.3.2 Dimensionless Variables

This problem can be expressed in a dimensionless form. First, we define¹ u = −R/r and w(u) = V_r(r)/c_s with the constants $${\gamma }_1=\frac{GM}{c_{\text{s}}^{2}R},$$(17) $${\gamma }_2=\frac{1}{c_{\text{s}}^2}{\left(\frac{{\ell }_{\text{R}}}{R}\right)}^2,$$(18)

and $${\gamma }_3=\frac{\alpha c_{s}R}{{\ell }_{\text{R}}}.$$(19)

Then, we substitute Eqs. (11) and (16) in Eq. (9), and we obtain the dimensionless equation of motion (hereafter EoM) $$\begin{array}{ll}\left(1-w^2\right)\frac{\text{d}w}{\text{d}u}=wF\left(u,w;\lambda \right)\hfill & -1<u<0\hfill \end{array},$$(20)

where the function F ≡ F(u, w; λ) is defined as $$F\left(u,w;\lambda \right)=\frac{5}{2}{u}^{-1}+{\gamma }_1[1-\eta \left(-u\right)-\in ]+{\gamma }_2{\left[1+{\gamma }_3\left(\lambda +\frac{1}{uw}\right)\right]}^{2}u.$$(21)

We redefine the eigenvalue of this problem as $$\lambda =\frac{c_s}{V_{\text{r}}}=\frac{1}{w\left(-1\right)}.$$(22)

Thus, the initial condition at the stellar surface for the EoM is $$w\left(-1\right)=1/\lambda .$$(23)

In order to find a transonic wind solution (i.e., a solution that starts with a low speed value at the stellar surface and reaches a speed higher than the sound speed at large distances), the solution of this problem must pass through a singularity from the EoM. In addition, the dependency on λ is nonlinear in Eq. (20) and also in the initial condition. Therefore, in the next section we analyze the topology at the singularity.

3 Topology of the EoM

Following the detailed topological analysis of the CAK (Curé & Rial 2004) and m-CAK (Curé & Rial 2007) cases, in order to obtain a transonic solution the location of the singular point is obtained after the analysis of the singularity condition. In addition, to assure a smooth transition from the subsonic solution branch to the supersonic branch, an extra condition must be imposed, namely the regularity condition.

3.1 Singularity and Regularity Conditions

Rearranging Eq. (20) and defining a function H, we obtain $$H\equiv H\left(u,w,w\text{'};\lambda \right)=\left(1-w^2\right)w\text{'}-wF\left(u,w;\lambda \right),$$(24)

where w′ = dw/du. Then, the mathematical definition of the singularity condition reads $$\frac{\partial }{\partial w\text{'}}H=0,$$(25)

which gives (1 − w²) = 0, or the corresponding decretion solution w = +1.

The regularity condition must be imposed at this singular(or sonic) point. This condition is $$\frac{\text{d}}{\text{d}u}H=\frac{\partial }{\partial u}H+w\text{'}\frac{\partial }{\partial w}H=0.$$(26)

Equations (25) and (26) are only valid at the sonic point: u = u_s, w(u_s) = 1, and w′(u_s) = w′_s. Solving w′_s from Eq. (26), we obtain an analytical expression for the value of the velocity gradient (w′_s) at this point: $$w_s^{\text{'}}=-\frac{1}{4}{F}_w\left(u_s,1;\lambda \right)\pm \frac{1}{4}\sqrt{F_w{\left(u_s,1;\lambda \right)}^2-8F_u\left(u_s,1;\lambda \right)}.$$(27)

Here F_u(u_s, 1; λ) = ∂F/∂u and F_w(u_s, 1; λ) = ∂F/∂w. Hereafter, we use F_u and F_w instead of F_u(u_s, 1; λ) and F_w(u_s, 1; λ), respectively. The partial derivatives of F are as follows: $$\begin{array}{ll}{F}_u=\hfill & -\frac{5}{2u^2}-{\gamma }_1\eta ϵ{\left(-u\right)}^{-\left(1+ϵ\right)}\hfill \\ \hfill & +{\gamma }_2{\left(1+{\gamma }_3\left(\lambda +\frac{1}{u}\right)\right)}^2\hfill \\ \hfill & -2{\gamma }_2{\gamma }_3{u}^{-1}\left(1+{\gamma }_3\left(\lambda +\frac{1}{u}\right)\right)\hfill \end{array}$$(28)

and $$F_w=-2{\gamma }_2{\gamma }_3\left(1+{\gamma }_3\left(\lambda +\frac{1}{u}\right)\right).$$(29)

3.2 Solution Branches

A necessary condition for the existence of w′_s > 0 (decretion) comes from the argument of the square root in Eq. (27): $$\Delta F={\left(F_w\right)}^2-8F_u\ge 0.$$(30)

With this condition (Eq. (30)), from Eqs. (21) and (27) there exist two branches for w′_s, denoted w′_s+ and w′_s−, where w′_s+ ≥ w′_s−. When ∆F = (F_w)² − 8F_u = 0, the two branches are equal: w′_s+ = w′_s−.

3.3 Classification of the Singular Point

The standard classification of singular points, following Amann (1990), is shown in Table 1. To exemplify this classification we select throughout this work the same stellar parameters from Okazaki (2001) for a B0 main-sequence star: M = 17.8 M_⊙, R = 7.41 R_⊙, and T_eff = 28 000 K.

Figure 1 shows contour plots in terms of r/R and λ for the following conditions: the black solid line corresponds to F = 0, F_u = 0 is shown as a short-dashed straight line, and $F_w^2-8F_u=0$ as a long-dashed straight line. The spiral topology corresponds to the segment of the curve F = 0 that is to the left of the long-dashed straight line, the nodal topology corresponds to the segment that lies between the two straight lines, and the saddle (or X-type) topology corresponds to the segment that is to the right of the short-dashed straight line. In the inset of this figure, which shows a zoomed-in image, we see that we can restrict the range of the eigenvalue λ in order to solve numerically the EoM.

Fig. 1 Different segments for the singular point classification (see text for details). The parameters used are Ω = 0.9, α = 1.0, ϵ = 0.1, and η = 0.5.

3.4 Determining the Range of the Eigenvalues

From Fig. 1, defining ${\widehat{r}}_{s1}$ and λ₁ as the values of r_s/R and λ at the intersection between F = 0 and $F_w^2-8$ F_u = 0, we obtain ${\widehat{r}}_{s1}=65.74$ and λ₁ = 129.22. Similarly for the other intersection, we obtain ${\widehat{r}}_{s2}=74.25$ and λ₂ = 130.82. Using this information we can search for the values of the velocity gradient (w′_s±) and the eigenvalue (λ) at the sonic point given by Eq. (27).

Figure 2 shows the value of w′_s as a function of the eigenvalue. In order to obtain this value we need to set the value of λ and then evaluate w′_s± for the values of ${\widehat{r}}_s\left(=1/u_s\right)$ in the range ${\widehat{r}}_{s1}<{\widehat{r}}_s<{\widehat{r}}_{s2}$ for the nodal solution (shown as the black dashed line for w′_s+ and the gray dashed line for w′_s−). On the other hand, for the saddle solution, we set ${\widehat{r}}_s$ in the range ${\widehat{r}}_s>{\widehat{r}}_{s2}$. The solid black line shows the branch w′_s+ and gray solid line shows the branch w′_s−, which is not a physical solution for a transonic decretion disk, because w′_s− < 0.

Fig. 2 Velocity gradient at the singular point as a function of the eigenvalue λ. The solid lines are for the saddle solutions and the dashed lines for the nodal solutions. The parameters used are the same as in Fig. 1 (see text for details).

4 Numerical Procedure

Knowing the topology of the nonlinear EoM, we can now describe the numerical procedure to obtain the different viscous transonic decretion disk solutions. Our proposed procedure is the following. First, we define a grid of eigenvalues (λ_i) in the range (λ₁, λ₂) for nodal solutions and (aλ₁, λ₂) for the saddle solution. A typical value of a is 0.7, according our calculations. Second, we integrate from the stellar surface u = −1, with the initial condition w(−1) = 1/λ_i, up to w(u) = 1 −-ϵ₀. When this value of the velocity is attained, the location of the singular point u = u_s is therefore known. A typical value of ϵ₀ is 10⁻⁶. Third, for the entire grid of λ_; solutions that reached the sonic point, we compare the values of the numerical velocity gradient at the sonic point (w′_s)_i with the values of w′_s− or w′_s+, depending on the solution type, obtained from Eq. (27), calculating the absolute error ϕ_i(λ_i) = ‖(w′_s)_i −-w′_s±‖. It should be noted that when calculating w′_s± from this equation the value of u_s must lie in the segment of the nodal or saddle solutions depending on the type of solution sought (see Fig. 1). Finally, the eigenvalue λ_i corresponds to min(ϕ_i).

Based on this numerical methodology we are now in a condition to perform numerical experiments to understand the behavior of the solutions in terms of the different parameters involved in the EoM.

5 Numerical Calculations

In this section we perform numerical calculations to understand the influence of the different parameters involved in the EoM. All our calculations hereafter are performed for a typical B0 main-sequence star: M = 17.8 M_⊙, R = 7.41 R_⊙, and T_eff = 28 000 K with T_disk = 0.5 T_eff and Ω = 0.9.

5.1 Solutions without Line Force

Following Okazaki (2001), we solve first the EoM without line force, and therefore set η = 0 to describe a pure viscous outflow. Figure 3 shows the normalized velocity profile w(r) as a function of r/R. The eigenvalues are λ = 771.10; 3734.65; 28498.51 for α = 1.0; 0.1; 0.01, respectively. The location of the sonic point is similar for all solutions and lie in the range 425 ≤ r_s/R ≤ 479 (see the inset in Fig. 3).

Fig. 3 Dimensionless velocity profile (w) for the case without line force (η = 0). The black solid line corresponds to α = 1, the black dotted line is α = 0.1, and the black dashed line is α = 0.01. The inset shows a zoomed-in image of the region around the sonic point. The light gray horizontal solid line represents the sound speed value (w = 1). The stellar parameters are given at the beginning of Sect. 5.

Fig. 4 Dimensionless velocity profiles (w) for the two indistinguishable nodal solutions as a function of the r/R coordinate, with α = 1.0 and line-force parameters ϵ = 0.1 and η = 0.5. The inset shows a zoomed-in image of the two solutions around the singular point. The black solid line correspond to the w′s+ branch solution and the dashed line to the w′s− branch solution. The light gray solid line represents the sound speed value (w = 1). The stellar parameters are given at the beginning of Sect. 5.

5.2 Viscosity Parameter α

Analogously to the results found by Okazaki (2001), the main impact of the viscosity parameter α is the topology of the singular point.

5.2.1 Nodal Solutions

High values of α implies that the sonic point has a nodal topology. This topology shows two branches (dashed lines in Fig. 2). Figure 4 show the velocity profiles for both nodal solutions. The viscosity and line-force parameters are α = 1.0, ϵ = 0.1, η = 0.5, and T_disk = 0.5 T_eff. The two nodal solutions are indistinguishable from each other on the scale of Fig. 4, but in the zoomed-in image around the singular point it is possible to distinguish each one separately. The solution with w′_s+ (black solid line) has λ = 223.15 and r_s = 107.78 R; the w′_s− solution (dashed line) has λ = 222.91 and r_s = 106.93 R. For any practical purpose, there is almost no difference in the behavior of the two nodal solutions, and they have almost the same eigenvalue; therefore, we can select either of them.

Fig. 5 Dimensionless velocity profile (w) (solid black line) for the saddle solution with α = 0.2 and line-force parameters ϵ = 0.1 and η = 0.5. The light gray solid line represents the sound speed value (w = 1). The stellar parameters are given at the beginning of Sect. 5.

5.2.2 Saddle Solution

Depending on the value of α, the topology switches between saddle and nodal. Using same parameters as the previous nodal case, but with α = 0.2, the resulting saddle solution is shown in Fig. 5. We clearly see here that the location of the sonic point is almost the same as both nodal solutions; here we have λ = 679.63 and r_s = 109.62 R.

Although the nodal and saddle solutions, shown in Figs. 4 and 5, only differ in the values of the viscosity parameter α and, consequently, in their topology, the location of the sonic point hardly differs and the velocity profiles w(r) are very similar.

In the next subsections we investigate in detail the role of the different parameters in this viscous transonic decretion outflows.

5.3 Type of Solution Depending on the Value of η

The line-force strength and behavior is determined by the η and ϵ parameters, respectively. Here we study the variation of the η parameter with α = 0.2. Figure 6 shows w(r) for different values of the η parameter with ϵ = 0.1. We confirm that the larger η is, the nearer the sonic point is to the stellar surface. The specific locations are r_s/R = 614.8; 109.6; 43.9 for η = 0.1; 0.5; 0.6, respectively. When η = 0.6 the sonic point is located at a distance lower than 50 R, but such a strong force is very unlikely. In Sect. 5.5 we discuss other combinations of line-force parameters and disk temperatures that give similar results about the location of the singular point.

We studied in terms of the α parameter when the solution is nodal or saddle as a function of the η line-force parameter. We find that the switching zone is very narrow, as shown in Fig. 7 and summarized in Table 2.

All the solutions we obtained from the EoM are physical solutions. The characterization of the stability of the steady state can be analyzed by means of the evolution of perturbed waves in the time-dependent equation of motion (see different approaches in Criminale et al. 2018, and references therein), but this type of study is beyond the scope of this work.

5.4 Type of Solution Depending on the Value of Ω

In order to study the influence of the centrifugal force on this viscous disk decretion model, we calculate for different values of Ω the dimensionless velocity profile w as function of r (see Fig. 8). These w(r) profiles show a very unexpected result. They all converge to a unique solution after a very small distance above the stellar surface, r < 1.0004 R, as shown in the inset. Similarly, Fig. 9 shows the behavior of V_ϕ(r)/V_crit as function of r. Again, all solutions converge to a unique solution very near the stellar surface.

A steep gradient in V_ϕ(r) might cause a boundary layer, provoking that the specific angular momentum just outside of this layer is settled into a value, which is smoothly connected to the outer parts of the disk. The results shown in Fig. 9 can be used to constrain the value of Ω (or the specific angular momentum) to a much smaller region (i.e., to obtain a smooth solution, in this case Ω ~ 0.7, where the gradients are not steep).

Fig. 6 Dimensionless velocity profile w as function of r/R with ϵ = 0.1 and α = 0.2. Here η = 0.1 is shown as a black solid line, η = 0.5 as a black dotted line, and η = 0. 6 as a black dashed line. The light gray solid line represents the sound speed value w = 1 . The stellar parameters are given at the beginning of Sect. 5.

Fig. 7 Values of α for nodal and saddle regions as function of η with ϵ = 0.1. The solid black line denotes the lower limit of α for nodal solutions. The light gray line denotes the upper limit of saddle solutions. Hatched zones represents the nodal (upper) and saddle (lower) regions. The stellar parameters are given at the beginning of Sect. 5 (see text for details).

Fig. 8 Dimensionless velocity profile, w as function of r/R for different values of Ω with α = 0.2, ϵ = 0.1, and η = 0.5. The gray solid line is Ω = 0.99, the dot-dashed line is Ω = 0.9, the dotted line is Ω = 0.8, the dashed line is Ω = 0.7, and the solid black line is Ω = 0.6. The light gray horizontal line represents the sound speed value, w = 1 . The stellar parameters are given at the beginning of Sect. 5 (see text for details).

Fig. 9 Vϕ(r)/Vcrit as function of r for the same set of parameters a Fig. 8. All solutions converge to one in a very small region above th stellar surface (see text for details).

5.5 Type of Solution Depending on the Line-Force Parameters and T_disk

We analyze the behavior of the ϵ parameter and the disk temperature T_disk in the VDD model. The parameter ϵ describes the decay of the line force moving outward from the stellar surface.

Table 3 summarizes the location of the sonic point and eigenvalue in terms of the value of ϵ. The dependence of the line force in terms of r, when Г → 0, is given by (see Eq. (6)) $$g_{\text{rad}}\to \frac{GM}{r^2}\eta {\left(\frac{r}{R}\right)}^{\in }.$$(31)

Thus, the larger ϵ is, the slower its decay is in terms of r, and due to a greater line force the location of the sonic point lies closer to the stellar surface.

Table4 summarizes the location of the sonic point and eigen-value in terms of the value of the disk temperature T_disk and different combinations of line-force parameters η and ϵ. The influence of T_disk in the location of the sonic point is not determinant when the value of ϵ is high. On the other hand, for low values of ϵ the influence of T_disk is quite important, reducing the location of r_s from ~240R to ~ 160R when T_disk/T_eff increases from 0.5 to 0.9 when η = ϵ = 0.2.

As we pointed out in Sect. 5.3, it is very unlikely to have such a strong line force with η = 0.6. However, it is indeed possible with plausible values of ϵ, η, and T_disk to obtain transonic disk solutions that might explain the results of Klement et al. (2017).

6 Discussion

In this work we showed the dependence of the velocity field for different parameters involved in the EoM from the VDD model. However, the standard methodology used to obtain an observable, such as Hα, is to use the (volumetric) density (ρ) or the vertical integrated density (Σ) as a function of r/R in a Keplerian orbit as input in radiative transport codes, such as BEDISK (Sigut & Jones 2007) or HDUST (Carciofi & Bjorkman 2006). The standard modeling of ρ is $$\rho \left(r\right)={\rho }_0{\left(\frac{r}{R}\right)}^{-m}{e}^{-{\left(z/H\right)}^2},$$(32)

where ρ(R) = ρ₀ is the density at the stellar surface. The definition of Σ follows after integrating ρ(r) in the z direction: $$\Sigma \left(r\right)={\Sigma }_0{\left(\frac{r}{R}\right)}^{-n}.$$(33)

Here² n = m − 3/2, and Σ(R) = Σ₀ is the vertical integrated density at the base of the wind.

Figure 10 shows different values of Σ(r)/Σ(R) as a function of r for the models shown in Sect. 5.3 (see also Fig. 6). The solid gray line represents the fit of Eq. (33) to the vertical integrated density profile for η = 0.6. We clearly see in this figure that none of the solutions maintains the behavior described by Eq. (33) in the entire range of r, but only in the range R < r ≲ 10 R. The results for the fits are n = 0.745 for η = 0.6, n = 0.649 for η = 0.5, and n = 0.554 for η = 0.1. All these values of n are calculated in the interval R ≤ r ≤ 20 R. Fits for the interval R ≤ r ≤ 100 R, give values of n that differ in the fourth decimal place with the previous fit interval. These results shows that the standard formulae giving by ρ(r) or Σ(r) are fairly good approximations to our results, but only for r ≲ 10 R.

A possible explanation for the different behaviors shown in Fig. 10, with different locations of the sonic point for each η (see Sect. 5.3) is the following. Even if the Shakura–Sunyaev viscosity model is applicable to supersonic regions, the viscosity is just as inefficient there as the angular momentum transfer mechanism because the advection timescale in the supersonic region is much shorter than the viscous timescale, which makes the flow to be angular-momentum conserving.

In addition, Klement et al. (2017) studied the observed SED turn-down by means of the VDD model, concluding that it can be explained with a truncated disk. They argued about two possible explanations of this: tidal forces from a close binary companion or a velocity profile with a transonic transition not too far from the star. They assumed that a binary companion is the most probable scenario, but for most of their six sources binarity remained undetected. Nevertheless, we showed in this work, that for a suitable combination of η ϵ, and/or T_disk it is possible to have singular point locations that explain the SED turn-down without needing a close companion.

In reference to the influence of the stellar rotational speed (in terms of the critical speed) Ω, we found that the viscosity effects collapse all the solutions to (almost) unique w(r) and υ_ϕ(r) profiles in a very small region above the stellar surface, as shown in Figs. 8 and 9. The results from υ_ϕ(r) shown in Fig. 9 can be used to restrict the range of Ω leading to a smooth solution in order to avoid the formation of a boundary layer.

7 Conclusions

We revisited the Okazaki (2001) work of a viscous transonic decretion disk model with the inclusion of a radiative acceleration produced by an ensemble of optically thin lines, described by the Chen & Marlborough (1994) model. We developed a new solution approach, where the eigenvalue of the problem is no longer an angular quantity, but a radial quantity. After a detailed topological analysis of the steady-state equation of motion, three possible physical solutions were found: two nodal solutions, where the viscosity parameter α is larger (nodal region in Fig. 7), and one saddle solution for lower values of α (saddle region in Fig. 7). The two nodal solutions are almost indistinguishable from each other. The value of the viscosity parameter α, given by the Shakura–Sunyaev model, is not determinant for the location of the sonic point. Other parameters, especially those of the line force, η, and ϵ, directly influence the solution of the EoM. Finally, any effect of the stellar rotation is rapidly damped close to the stellar surface due to viscosity. In addition, in order to obtain only smooth solutions the range of Ω should be restricted.

As a future work we will describe the line force using the standard CAK (and its improvements) theory (Castor et al. 1975) instead of the ad hoc model from Chen & Marlborough (1994). With this more detailed model for the line acceleration, it will be possible to incorporate rapid rotational effects such as oblate shape and gravity darkening to better describe a transonic VDD model of Be stars. Finally, to study the observable such as the SED or the line Hα (see Klement et al. 2017) we plan to use the density profile Σ(r) as input in BEDISK and/or HDUST to calculate these observables.

Fig. 10 Σ(r)/Σ(R) as function of r with α = 0.2 and ϵ = 0.1 (all stellar parameters are given at the beginning of Sect. 5). Here η = 0.1 is shown as a black solid line, η = 0.5 as a black dotted line, and η = 0.6 as a black dashed line. The light gray line represents the fit of Eq. (33) for the η = 0.6 case (see text for details).

Acknowledgements

The authors would like to thank the referee, Atsuo Okazaki, for his thoughtful comments and suggestions to improve this work. M.C. and C.A. acknowledge the support from Centro de Astrofísica de Valparaíso. M.C., C.A. and I.A. thank the support from FONDECYT project 1190485. M.C. and C.A. thank to project ANID-FAPESP 2019/13354-1. IA is also grateful for the support from FONDECYT project 11190147. C.A. thanks the support from FONDECYT project 11190945. This project has also received funding from the European Unions Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie grant Agreement N° 823734. This work has been possible thanks to the use of AWS-U.Chile-NLHPC credits. Powered@NLHPC: this research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).

References

All Tables

All Figures

	Fig. 1 Different segments for the singular point classification (see text for details). The parameters used are Ω = 0.9, α = 1.0, ϵ = 0.1, and η = 0.5.
In the text

	Fig. 2 Velocity gradient at the singular point as a function of the eigenvalue λ. The solid lines are for the saddle solutions and the dashed lines for the nodal solutions. The parameters used are the same as in Fig. 1 (see text for details).
In the text

Fig. 3 Dimensionless velocity profile (w) for the case without line force (η = 0). The black solid line corresponds to α = 1, the black dotted line is α = 0.1, and the black dashed line is α = 0.01. The inset shows a zoomed-in image of the region around the sonic point. The light gray horizontal solid line represents the sound speed value (w = 1). The stellar parameters are given at the beginning of Sect. 5.

In the text

Fig. 4 Dimensionless velocity profiles (w) for the two indistinguishable nodal solutions as a function of the r/R coordinate, with α = 1.0 and line-force parameters ϵ = 0.1 and η = 0.5. The inset shows a zoomed-in image of the two solutions around the singular point. The black solid line correspond to the w′s+ branch solution and the dashed line to the w′s− branch solution. The light gray solid line represents the sound speed value (w = 1). The stellar parameters are given at the beginning of Sect. 5.

In the text

	Fig. 5 Dimensionless velocity profile (w) (solid black line) for the saddle solution with α = 0.2 and line-force parameters ϵ = 0.1 and η = 0.5. The light gray solid line represents the sound speed value (w = 1). The stellar parameters are given at the beginning of Sect. 5.
In the text

	Fig. 6 Dimensionless velocity profile w as function of r/R with ϵ = 0.1 and α = 0.2. Here η = 0.1 is shown as a black solid line, η = 0.5 as a black dotted line, and η = 0. 6 as a black dashed line. The light gray solid line represents the sound speed value w = 1 . The stellar parameters are given at the beginning of Sect. 5.
In the text

	Fig. 7 Values of α for nodal and saddle regions as function of η with ϵ = 0.1. The solid black line denotes the lower limit of α for nodal solutions. The light gray line denotes the upper limit of saddle solutions. Hatched zones represents the nodal (upper) and saddle (lower) regions. The stellar parameters are given at the beginning of Sect. 5 (see text for details).
In the text

Fig. 8 Dimensionless velocity profile, w as function of r/R for different values of Ω with α = 0.2, ϵ = 0.1, and η = 0.5. The gray solid line is Ω = 0.99, the dot-dashed line is Ω = 0.9, the dotted line is Ω = 0.8, the dashed line is Ω = 0.7, and the solid black line is Ω = 0.6. The light gray horizontal line represents the sound speed value, w = 1 . The stellar parameters are given at the beginning of Sect. 5 (see text for details).

In the text

	Fig. 9 Vϕ(r)/Vcrit as function of r for the same set of parameters a Fig. 8. All solutions converge to one in a very small region above th stellar surface (see text for details).
In the text

	Fig. 10 Σ(r)/Σ(R) as function of r with α = 0.2 and ϵ = 0.1 (all stellar parameters are given at the beginning of Sect. 5). Here η = 0.1 is shown as a black solid line, η = 0.5 as a black dotted line, and η = 0.6 as a black dashed line. The light gray line represents the fit of Eq. (33) for the η = 0.6 case (see text for details).
In the text