Modeling of the Wind/Disk Outﬂow from Be Stars

: The objective of this work is to reproduce the formation of the fast polar wind and viscous disk outﬂow from Be stars in a uniﬁed physical picture. Numerical modeling of the plasma outﬂow from fast rotating stars was performed taking into account the acceleration of the plasma due to scattering of the radiation of the star in lines of plasma ions and excitation of the hydrodynamic turbulence in the outﬂow. The fast polar wind naturally arises in this picture with an expected ﬂow rate. For the ﬁrst time, it is shown that a disk-like outﬂow with a relatively high level of turbulence is formed at the equator of fast rotating stars emitting radiation-driven wind. However, the level of turbulent viscosity is well below the level necessary for the formation of a Keplerian disk.


Introduction
Be stars have attracted special interest in recent years due to the discovery of transient radiation from binary systems containing a massive O or B class star and a pulsar or a black hole [1]. Binary system PSR1256-63/LS2883 is one of the representatives of such objects. It consists of a pulsar and massive O-type companion [2][3][4]. Interaction of the pulsar wind with the outflow from the Be star leads to radiation in the range from the radio to hard gamma rays with photon energy up to 1 TEV [5][6][7][8][9][10][11][12][13]. GeV flares 30 days after the passage of the periastron in 2010 [14], 2014 [13] and 2017 [15,16] are one of the most mysterious phenomena in this system. So far, the problem of the flares has not been solved [17,18]. One of the reasons for this could be our poor understanding of the physics of the outflow from Be stars.
Be stars produce an outflow with the rate of mass loss in the range 10 −10 -10 −7 M , yr −1 [19,20]. The outflow consists of two components, fast polar wind and slow equatorial flow. The mass flow in the disk can vary in time in the range 10 −10 to 10 −7 M , yr −1 [21][22][23]. Typically, the rotation speed of Be stars is 70% of the critical speed [24]. Sometimes, it can reach 90% of the critical speed [25,26].
At present, there is a common consent regarding physics of the polar wind. The polar wind from Be stars is accelerated by radiation from the star due to the scattering of photons in lines of ions of the wind. The radiative force provides a high velocity of wind and a high rate of mass loss in agreement with observations.
The situation with the disk outflow is much worse. Theoretical arguments [27][28][29] and observations [30][31][32] are in favor of the model of viscous disk-like outflow at the equatorial plane-decretion disk. In this disk, transport of the angular momentum from the rotating star is provided by turbulent viscosity. This transport is accompanied by matter outflow from the surface of the star [33] to the outer edge of the disc [34]. This leads to the formation of a quasi Keplerian disk [35,36]. However, up to now, there are no convincing solutions of the two most important problems of Be star disks [24,37,38]. One of the largest unknowns is how to supply the disk by matter. Although the star rotates close to 70% of the critical speed, this is not sufficient to provide the rate of the mass flow comparable with the observed one. Our recent simulations took into account only gravitation and rotation. They gave much lower flow rate than the observed one [39].
There are only a few mechanisms to eject gas from the photosphere with a high flow rate. Photospheric pulsations can provide ejection of the gas [40]. However, it is not clear if the pulsations have amplitudes large enough to lift gas off the star. Subphotospheric motion can produce outbursts of material [41]. The magnetic field can also amplify mass loss [42].
The next problem is how to provide transfer of the angular momentum from the star to the disk. The star surface rotates slower than the inner edge of the Keplerian disk. To overcome this problem, the mechanism of the angular momentum transfer by waves [43] or by magnetic field [44,45] was proposed. However, the magnetic field should be very strong to provide an efficient angular momentum transfer. Therefore, a strong magnetic field is observationally ruled out [44].
Without a solution to these problems, it is difficult to consider the idea of the viscous decretion disk as the completed theoretical model. Up to now, the physics of the polar wind [27][28][29] and that of the decretion disk are considered separately [46][47][48] as different physical phenomena. At the same time, the same physical processes take place in both polar wind and decretion wind. They are line acceleration, centrifugal force, gravitational force, gas pressure and turbulence. It is interesting to consider how these processes result in the formation of completely different types of outflow in the polar region and slow quasi Keplerian flow at the equator in one picture. Some ideas in this regard were developed in [49], where the model of the wind compressed disk was proposed, and in [50], exploring the effect of bi-stability of the wind at the threshold of effective temperature of the 25,000 K star. However, these works took into account only radiative force. It is impossible to obtain the viscous disk in this way. The role of turbulence, which certainly plays key role in formation of the decretion disks, still has not been investigated. The first step in this process has been conducted in our recent work [39], where we considered the outflow, neglecting the process of the line accretion of the wind. We found that all the outflow occurs in this case from a small region of the star surface located at the equator. The turbulence results in the formation of the disk-like outflow. However, the flow rate in this model appeared to be too low in comparison with the observed values. In this work, we added the process of the line acceleration in the simplest CAK form [51] with constant force multipliers and neglecting gravity darkening, finite size of the star and azimuthal and polar components of the radiative force to our model. Taking into account that the turbulence has to dramatically modify the dynamics of plasma at the equator, it is unlikely that all the neglected processes can essentially affect this modification.
The paper is organized as follows. In the second section, the model and basic equations determining the dynamics of the wind and turbulence are described, including the test problem used for the verification of the code and model. In the third section, the main results are presented, which are discussed in the last section.

Methods
Calculations are performed in the frameworks of the model developed earlier in [39]. The main new process included in the model is the radiative force due to the scattering of photons radiated by the star in lines of ions of the matter of the wind.

Radiative Acceleration
The radiative force f rad can be separated into two components: f C due to Thompson scattering of photon continuum radiation on electrons and f L due to scattering of photons in lines of ions. f C can be presented in the following form: where ρ is the density of the wind, G is gravitational constant, r is the radius vector from the center of the star to a point, Γ is the continuum Eddington factor defined as which represents the ratio of the star's luminosity L * over the Eddington luminosity. Here, c is the light velocity, M * is the stellar mass, or σ T is the Thompson mass crosssection.
The radiative force f L is determined in the CAK model [51] as a sum over all the lines where k L is the mass absorption coefficient of an individual line, ν 0 is the center line frequency, and v th is the local thermal speed of the ions. Sobolev optical depth τ is defined as follows: Generally speaking, this force is not radial, provided that the wind has azimuthal velocity and the radiation from the star is not uniform due to gravity darkening [28,52]. We neglect these effects. In this work, we use widely accepted simplified point source approximation in which the line force is presented in the following form: where k, α, and δ are force multiplier coefficients. As typical parameters for Be star LS2883, we accept k ∼ 0.15, α ∼ 0.64 and δ ∼ 0.08 [53]. W(r) is the dilution factor defined as [51] In fact, force multiplier coefficients are not constants. Accurate calculations in non-LTE models show that all the varies with coordinates depend on the effective temperature of the star T e f f [29,54,55]. However, for our goal of the investigation how hydrodynamical turbulence affects the dynamics of the wind, original CAK approximation is enough to draw unambiguous conclusions. For the same reason, the deviation of the factor 10 −11 n e W δ from 1 is neglected. However, the oblateness of the star due to rotation was taken into account.

System of Equations
The gravitational attraction, rotation of the star, pressure gradients, turbulence and line and continuum acceleration are taken into account. The problem is solved in the rotating frame system. The system of equations consists of the equation of continuity and Navier-Stocks equation [56] ∂ρv i ∂t where µ T is an eddy dynamical viscosity, and p and ρ are the pressure and density, respectively. ϕ is the total gravitational and centrifugal potential of the following form Here, ω is the angular velocity of the star; R and r are the polar and cylindrical radius, respectively; M * is the mass of the star; and G is the gravitational constant. The turbulence is described in the simplest k-ε model. One equation for intensity of turbulence k is as follows: ∂ρk ∂t where where k is the turbulent viscosity, and C µ , C ε 1 , C ε 2 , σ k , σ E are the model constants [57].
The equation of state of the isothermal gas is p = ρC 2 , where C is constant sound velocity, which is connected to the effective temperature of the star T e f f as follows, where k B is the Boltzmann constant, and m p is the mass of the particles of the wind.

Dimensionless Variables
As in our previous work [39], we consider the problem in dimensionless variables. Velocities are expressed in V K -Kepler velocity at the pole of the star; the distance is measured by the radius of the star at the pole R pole ; density is measured in units of typical density of Be stars at the photosphere ρ ph ; and dimensionless mass flow rateṁ is defined as follows:ṁ =Ṁ/(4πρ ph V K R 2 pole ), (12) whereṀ is the rate of the mass flow. We changed the definition ofṁ in comparison with the definition provided in our previous work [39]. Parameters of Be star in the binary system PSR1256-63/LS2883 are given in [4,58,59] in a relatively wide range. For accuracy, we used the parameters collected in Table 1. They provide for the dimensionless squared sound velocity v 2 s = 5.7 × 10 −4 andṁ = 1.6 × 10 −3 . These values are used below only as the reference point in the space of the dimensionless parameters to understand where we are located in this space.
It is convenient to represent the line radiative force in dimensionless form as well. According to [60], without rotation, radially expanding cold wind driven only by radiation and gravitation has the following relationship between mass flow rateṀ 0 and luminosity of the star L * : Introducing dimensionless mass flow rateṁ 0 =Ṁ 0 /(4πρ ph V K R 2 pole ), the line radiative force in dimensionless variables can be rewritten in the following form: The physical sense ofṁ 0 is evident. This is the dimensionless mass flux, which could create the nonrotating star of a given mass and luminosity.
Thus, the majority of the radiation-driven winds are defined by only 4 dimensionless parameters:ṁ 0 , v 2 s , ω/ω c and α, provided that the radiative force is taken for a point source in the CAK form [51] with constant coefficients. ω c is the critical angular speed of rotation.

Boundary Conditions
The flow is solved in axisymmetric approximation in the 2-dimensional plane (r,z). The computational domain includes upper and lower hemispheres. Rotation leads to the oblateness of the star. According to [38], the shape of the stellar surface is determined by the following equation: To reproduce this shape, nonorthogonal mesh was used, unlike the conventional approach used in [29,49,61], where the orthogonal mesh was used in the spherical system of coordinates and the shape of the surface was reproduced by the following steps.
At the surface of the star, we specify only the dimensionless density equal to 1 and temperature, which is defined by v s . We do not specify the velocity of plasma here. The value of the plasma velocity is defined in the process of solution. Therefore, the density of the mass flux from every point of the stellar surface is determined by the solution rather than specified by hand.
The outer surface is a sphere with a radius that greatly exceeds that of the critical points in the flow. According to [62], the critical points are located rather close to the star at the distance ∼1.5. We considered the radius of the outer surface equal to 30.

Verification
ANSYS CFX package [63] was used for the numerical modeling. The code and the computer model were verified on the problem of the spherically symmetric outflow driven by radiative force from the nonrotating star. For example, we present the dependence of the ratio V R on R in Figure 1. The calculations were performed for the following parameters: and v 2 s = 2 × 10 −2 (dashed line). According to [62], the maximum of these curves provides the location of the critical radius R c of the flow. The conditions at R c define the mass flow rate from the star: in the first case, R c = 1.63 and in the second case R c = 1.8.
Dimensionless mass flux equalsṁ = 2.46 × 10 −2 at v 2 s = 8 × 10 −3 andṁ = 2.56 × 10 −2 at v 2 s = 2 × 10 −2 . The difference in the mass flux is due to difference in the sound velocity. According to [62] (note the misprint in Equation (8.121) of this work), mass flux from the star with first corrections due to nonzero sound velocity equalṡ According to this equation, the mass flux should be 2.453 × 10 −2 at v 2 s = 8 × 10 −3 and 2.546 × 10 −2 at v 2 s = 2 × 10 −2 . In all cases, the difference between the calculated mass flux and analytical predictions is less than 1%. The coincidence of the predictions with the calculated mass fluxes makes us confident in the validity of the code and the model.

Results
We present the results for the following parameters of the star. The angular speed of rotation of the star in dimensionless variables is ω = 0.43546, which corresponds to 80% of the critical speed. v 2 s = 4 × 10 −3 andṁ 0 = 2.4 × 10 −2 . The last parameters are almost one order of magnitude greater than the parameters typical for the Be star in the binary system PSR1256-63/LS2883. There are technical problems in the simulation of the radiation-driven outflow at smaller values ofṁ 0 and v 2 s . The gradients of the velocity and pressure at the star surface becomes so strong that any small deviation of the shape of the star surface from the equipotential one produces remarkable problems in the solution. An additional improvement of the model is needed for simulation of the flow at smallerṁ 0 and v 2 s . Figure 2 presents the structure of the flow in the upper hemisphere. The most impressive observation is that the flow dramatically differs from the flow obtained by us [39] when the radiative force was neglected. The flow lines are practically radial. However, the flow is not spherically symmetric. The density of plasma at the equator exceeds the density of the plasma at the pole. In this regard, our results reproduce the results of the works [29,49,61], where the radiative force was taken in the same form.

Structure of the Flow
Indeed, the excess of the density at the equator is due to compression of the flow to the equator. This conclusion is supported by the structure of the flow lines near the star surface shown in Figure 3. It follows from this figure that close to the star, there is a component of velocity directed to the equator, which results in tohe compression of the flow to the equator.  At our parameters, this compression is not very strong. Figure 4 shows the the distribution of the product ρ · R 2 along a line perpendicular to the equator at 3 different distances from the star: at radius 2,4 and 8 along the equator. It is seen that a relatively dense disk is formed. The ratio between the density at the equator and at the edge of the disk exceeds 6 at r = 8. The density of the mass flux is slower at the equator than at the pole because the radial velocity is reduced at the equator in compared with the velocity at the pole as it follows from Figure 5. This figure shows that the gradients of velocity (and density as well) are huge at the stellar surface, which is a challenge for any numerical code.

Role of Turbulence
The most important result of our study is the fact that the flow at the equator is unstable and, therefore, naturally produces turbulence. The instability of the flow is visible from numerical simulations of the flow without turbulence. However, the demonstration of these results is beyond the scope of the present work. Instability of the flow at the equator was expected. Analytical studies shows that at relatively large angular rotations, there is no steady state solution of the problem crossing the so-called critical point [64]. Incorporation in to the model of the equations of turbulence suppresses the time variability of the flow at the stellar surface, making it steady state. Figure 6 shows the distribution of the so called Shakura-Sunyaev viscosity parameter [65] α SS , defined as follows: where τ rϕ is the azimuthal component of the turbulent stress tensor defined as follows: µ T is taken from the solution of equations for k (10) and ε (11). The color shows α SS in the logarithmic scale. At the stellar surface, α SS achieves a comfortable value of 0.2. However, very quickly, α SS decreases to low values, which already cannot affect the dynamics of the flow. Thus, the turbulence is generated only very closely to the stellar surface. It is advected down stream of the flow with a damp-forming disk-like structure with a full opening angle close to 30 • . It is interesting that there is another source of turbulence in the flow. The interaction of the slow radial flow in the disk and fast radial flow in the polar wind leads to the Kelvin-Helmholtz instability amplifying turbulence at the boundary of the disk.  Figure 7 shows that the turbulence only slightly increases the angular momentum of the plasm very close to the stellar surface. Here, α SS achieves a remarkable value close to 0.2. Downstream, the wind the matter moves with practically constant angular momentum as with ordinary wind. Turbulent tensions created by a pure hydrodynamic turbulence are too small to provide a remarkable increase in the angular momentum.

Discussion
The basic objective of the work was to consider how the gravitation, fast rotation of the star, radiative force and turbulence can result in qualitatively different outflows: fast polar wind at the poles in comparison to the slow Keplerian viscous disk at the equator. For this, we firstly modeled the outflow driven by all these processes simultaneously.
We obtained denser and slower flow at the equator in comparison with the polar region. This agrees with all the works performed in the same approximation [29,49,61]. The formation of denser flow at the equator occurs due to the choice of the radiative force for a point-like source. Azimuthal and polar components of the force are equal to zero. More perfect calculations taking into account gravity darkening and finite size of the star lead to an opposite picture. According to [27,29], the flow is concentrated to the pole without the formation of any signature of the disk .
The most important result concerns turbulence excitation in the flow. We found that the turbulence is excited in the region at the equator forming a disk with a full opening angle 30 • . This was expected because analytical analysis predicts that at the fast rotation of the star, there is no steady state solution at the equator [64]. This is likely to mean that only nonstationary flow can exist near the equator. The level of turbulence at the stellar surface appears comfortably high-Shakura-Sunyeav α SS ∼ 0.2 here. However, the generation of the hydrodynamical turbulence is damped very quickly downstream of the flow. Therefore, we observed only a small increase in the angular momentum of the plasma close to the stellar surface. Almost everywhere downstream of the flow, the angular momentum of the plasma is practically conserved. The viscous tensions produced by the hydrodynamic turbulence is too small to provide a remarkable growth of the angular momentum of the plasma.
In this work, we intentionally studied the case of pure hydrodynamic turbulence as the existing observational data are in favor of very low magnetic field on the surface of Be stars. Therefore, it is reasonable to consider the case of pure hydrodynamic turbulence as the first step in the solution of the problem. We obtained that the hydrodynamic turbulence is excited near the equator. However, the level of turbulence is not high enough to provide an efficient transfer of the angular momentum from the star to the disk. This is also an expected result as the same situation takes place with accretion disks. The magnetic field plays a key role in the generation of turbulent viscosity in the accretion disks.
We already discussed the role of the chaotic magnetic field in the decretion disks [39]. This field should be generated inside the disk rather than advected from the surface of the star. Starting with the works by [65] and by [66], it was evident that a small-scale magnetic field generated inside the disk due to turbulent motion is responsible for the viscosity in the disk. Generation of the chaotic magnetic field due to the Magneto Rotational Instability (MRI) [67,68] is the most promising mechanism for the accretion and apparently for the decretion disks. Computer modeling of the magnetic field generation in the disks due to this mechanism [69,70] confirms this expectation.
It is interesting that at the temperatures corresponding to v 2 s ≈ 3 × 10 12 cm 2 s −2 , the value of the magnetic field which can produce turbulent viscosity with α ∼ 1 is estimated as B ∼ 4πρ ph v 2 s , which gives 20 G at the inner edge of the disk. On one hand, this field is too small for the detection by instruments from Earth. On the other hand, this field is high enough to affect the dynamics of interaction of the pulsar wind and 1 TeV electrons with the disk.