Issue 
A&A
Volume 569, September 2014



Article Number  A23  
Number of page(s)  7  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201424272  
Published online  11 September 2014 
Timedependent modeling of extended thin decretion disks of critically rotating stars
^{1}
Department of Theoretical Physics and AstrophysicsMasaryk
University,
Kotlářská 2,
611 37
Brno,
Czech Republic
email:
petrk@physics.muni.cz
^{2}
Institut für Physik und Astronomie, Universität
Potsdam, KarlLiebknechtStraße
24/25, 14476
PotsdamGolm,
Germany
Received:
26
May
2014
Accepted:
13
July
2014
Context. During their evolution massive stars can reach the phase of critical rotation when a further increase in rotational speed is no longer possible. Direct centrifugal ejection from a critically or nearcritically rotating surface forms a gaseous equatorial decretion disk. Anomalous viscosity provides the efficient mechanism for transporting the angular momentum outwards. The outer part of the disk can extend up to a very large distance from the parent star.
Aims. We study the evolution of density, radial and azimuthal velocity, and angular momentum loss rate of equatorial decretion disks out to very distant regions. We investigate how the physical characteristics of the disk depend on the distribution of temperature and viscosity.
Methods. We calculated stationary models using the NewtonRaphson method. For timedependent hydrodynamic modeling we developed the numerical code based on an explicit finite difference scheme on an Eulerian grid including full NavierStokes shear viscosity.
Results. The sonic point distance and the maximum angular momentum loss rate strongly depend on the temperature profile and are almost independent of viscosity. The rotational velocity at large radii rapidly drops accordingly to temperature and viscosity distribution. The total amount of disk mass and the disk angular momentum increase with decreasing temperature and viscosity.
Conclusions. The timedependent onedimensional models basically confirm the results obtained in the stationary models as well as the assumptions of the analytical approximations. Including full NavierStokes viscosity we systematically avoid the rotational velocity sign change at large radii. The unphysical drop of the rotational velocity and angular momentum loss at large radii (present in some models) can be avoided in the models with decreasing temperature and viscosity.
Key words: stars: massloss / stars: evolution / stars: rotation / hydrodynamics
© ESO, 2014
1. Introduction
The outflowing disks may be formed around various types of stars, such as Be stars, B[e] stars, and possibly luminous blue variables (e.g., Smith & Townsend 2007) and asymptotic giant branch stars (see, e.g., Matt et al. 2000; de Ruyter et al. 2006). Observational evidence supports the idea that these disks are Keplerian (rotationally supported) gaseous disks (Carciofi & Bjorkman 2008). The solutions for timeindependent viscous decretion disk structures where one assumes the disk to be isothermal and in vertical hydrostatic equilibrium (e.g., Okazaki 2001; Carciofi & Bjorkman 2008) show that, because of highly supersonic rotational velocity, the disks are geometrically very thin, and the disk opening angle is only a few degrees in the region close to the star. These studies naturally support the idea of the viscous decretion disk model leading to the formation of nearKeplerian disks around critically rotating stars. However, since the specific angular momentum in the nearKeplerian disk increases with the radius, decretion disks are unlikely nearKeplerian far from the star. A natural expectation is that a disk which is Keplerian near a star becomes angularmomentum conserving far from the star, although this transitional feature of decretion disks is not very well understood theoretically nor has it been observationally confirmed. Most of the models analyze the inner parts of the disk, whereas the evolution close to the sonic point or even in the supersonic regions up to the possible outer disk edge has not been, to our knowledge, very well studied.
In this paper we study the characteristics of the outflowing disks of critically rotating stars. The mass loss rate is determined by the angular momentum loss rate needed to keep the star at critical rotation (Krtička et al. 2011, see also Kurfürst 2012; Kurfürst & Krtička 2012). The basic scenario follows the model of the viscous decretion disk proposed by Lee et al. (1991, see also Okazaki 2001. In this model Lee et al. (1991) obtained a steady structure of viscous decretion disks around Be stars in thermal and radiative equilibrium. As the physics of accretion disks is quite similar, we follow the main principles (Pringle 1981; Frank et al. 1992) used in accretion disk theory for the description of decretion disks. The main uncertainties are the viscous coupling and the temperature distribution. Although some recent models indicate a constant value of the viscosity throughout such disks (Penna et al. 2012), we investigate the cases when the viscous coupling varies outward as a certain power law. The disk temperature is mainly affected by the irradiation from the central star (Lee et al. 1991). Motivated by the nonlocal thermodynamic equilibrium (NLTE) simulations (Carciofi & Bjorkman 2008) we extrapolate the temperature distribution up to quite distant regions as a power law.
2. Basic equations and parameterization
The disk is described in a cylindrical coordinate system assuming axial symmetry (e.g., Lee et al. 1991; Okazaki 2001; Krtička et al. 2011). The mass conservation (continuity) equation in the geometrically thin case is (we denote the cylindrical radius as R, the spherical radius as r) (1)where Σ is the surface density defined as (2)ρ is the density, and V_{R} is the radial component of velocity. The equation of radial momentum conservation in the thin disk approximation is (3)where V_{φ} is the azimuthal component of velocity, a is the speed of sound, a^{2} = kT/(μm_{u}), μ is the mean molecular weight (μ = 0.62 for the ionized hydrogen gas), m_{u} is the atomic mass unit, and M is stellar mass. The last term on the righthand side comes from the description of the gravitational force in a cylindrical coordinate system in the thin disk approximation, i.e., the disk thickness is negligible with respect to radial distance of any point studied (Matsumoto et al. 1984; Krtička et al. 2011; Kurfürst & Krtička 2012). From the conservation of azimuthal component of momentum we have (4)where f_{visc} means the viscous force per unit volume, exerted by the outer disk segment on the inner disk segment. In the axisymmetric (∂/∂φ = 0) and geometrically thin case the viscous force density is f_{visc} = ∂T_{Rφ}/∂R + 2T_{Rφ}/R (e.g., Mihalas & Mihalas 1984), where T_{Rφ} denotes the Rφ component of the stress tensor. Viscous force density f_{visc} is usually represented by using the first order linear viscosity term (see Eq. (8)) with the adopted ShakuraSunyaev α parameter (Shakura & Sunyaev 1973). We also examine the cases with nonconstant α parameters. Taking into account the turbulent motion of the gas, we can write for the kinematic viscosity ν (Frank et al. 1992) (5)with H denoting the typical vertical scaleheight of the disk, H^{2} = a^{2}R^{3}/ (GM) in Keplerian case.
The NLTE simulations (e.g., Carciofi & Bjorkman 2008) show that the radial temperature distribution in the very inner regions (up to few stellar radii) corresponds to a flat blackbody reprocessing disk due to the optically thick nature of this inner part , where T_{0} is the disk temperature at R = R_{eq}. As the disk becomes vertically optically thin, the disk temperature rises to the optically thin radiative equilibrium temperature with the average of about 60% of T_{eff}. The temperature radial profile at larger radii is nearly isothermal with relatively mild temperature decrease (about 1000 K from 10 to 50 stellar radii). Millar & Marlborough (1998, 1999) also found the radial temperature distribution up to 100 stellar radii to be nearly isothermal. We approximate these dependencies by a radial power law (6)where p is a free parameter (0 ≤ p< 0.5). Using the same power law temperature decline we also extrapolate the radial temperature structure of the outer part of the disk.
The viscosity also influences the temperature, but the contribution of the viscous heating in the disk is very small (practically negligible) compared with the heating that comes from radiative flux from the star (e.g., Smak 1989; Carciofi & Bjorkman 2008). The viscous heating dominates in the inner disk regions (Lee et al. 1991) only in the case of enormous value of mass loss rate (Ṁ ≥ 10^{5} M_{⊙} yr^{1}). This is the reason why the viscosity is parameterized via temperature independent α parameter.
The radial profile of α is not quite certain. The common agreement is that the value of the α parameter of turbulent viscosity is less than 1, since for α> 1 the rapid thermalization due to shocks would lead again to α ≤ 1 (Shakura & Sunyaev 1973). Most authors use the value around 0.1; some of the most recent works focused on the disk viscosity problem (e.g., Penna et al. 2012) find the α viscosity coefficient as constant, α ≈ 0.025. However, here we regard it as an open question if the α should be considered as constant throughout the disk or not. Therefore, we introduce (7)where α_{0} is the viscosity of the inner region of the disk near the stellar surface, n is a free parameter of the radial viscosity dependence, n> 0. The kinematic viscosity ν(R) is related to temperature via Eq. (5).
The first order linear viscosity term in Eq. (4)gives (8)while including the full viscosity term we obtain the righthand side of Eq. (4)in the form (9)The same expression is obtained by employing the angular momentum equation (Pringle 1981; Frank et al. 1992) (10)(used in numerical scheme; see Sect. 4.2), where Ω is the angular velocity (Ω = V_{φ}/R) and G is the viscous torque acting between two neighboring disk segments, (11)Most authors use the concept of some outer disk radius R_{out}, where, because of the radiation pressure for example, the disk matter may be completely driven outward (Lee et al. 1991). Using the reasonable parameterization of T and α, we are not dependent on the choice of R_{out} and can calculate the radial profiles of the characteristics as a direct solution of the hydrodynamic equations.
3. Analytical estimates of the radial thin disk structure
To obtain a general idea about the behavior of the main characteristics of the system we consider the stationary form of the basic hydrodynamic equations (Okazaki 2001; Krtička et al. 2011). Integrating Eq. (4)(using the term from Eq. (9)) and dividing this by the stationary mass conservation equation (RΣV_{R} = const.) we derive for α = const.(12)Inclusion of only the first order linear viscosity term (see Eq. (8)) gives a similar relation without the second term in the bracket. In the innermost part of the disk, V_{R} ≪ a. Consequently, the second term on the lefthand side of Eq. (12)dominates, therefore V_{R} ~ R and Σ ~ R^{2} from the continuity equation (Okazaki 2001). In the inner region of the disk the radial pressure gradient is negligible compared with the gravitational force, hence from the momentum Eq. (3)Keplerian rotation follows, V_{φ} ~ R^{−1/2}. Since Ṁ = const. (see Eq. (1)), we have in this region . Equation (12)becomes (13)where for Keplerian rotational velocity the numerical factor γ = 3/2. Obviously, the second order linear viscosity term in the Keplerian case represents one half of the corresponding first order linear term assuming the same α parameter, and therefore cannot be neglected. In the distant region near the sonic point (V_{R} ≈ a) the radial velocity increases; therefore, the first term on the lefthand side of Eq. (12)fully dominates, hence RV_{φ} = const. and the disk is angular momentum conserving.
In very distant supersonic regions with nearly flat disk temperature distributions (V_{R} ≫ a, V_{R} ≫ V_{φ}, ∂a^{2}/∂R ≈ 0) and negligible gravity the radial momentum Eq. (3)with use of mass conservation Eq. (1)implies logarithmic radial dependence, . The second term on the lefthand side of Eq. (13)rises, consequently the azimuthal velocity has to substantially decrease and may in general become even negative. This is, however, not possible according to the logarithmic term in Eq. (12). The numerical simulations (see Figs. 5–7) prove that the azimuthal velocity even in extremely distant regions does not change its direction. Equation (13)indicates that for steeper viscosity and temperature decrease (lower α and a^{2} in distant regions) the region of azimuthal velocity drop moves outwards.
From the stationary form of the radial momentum and mass conservation Eqs. (3) and (1)with the help of Eq. (6), the sonic point condition (Okazaki 2001; Krtička et al. 2011) is fullfilled at the sonic point radius (14)Substituting , where is the Keplerian velocity (Krtička et al. 2011, see also Figs. 1 and 2) we derive an estimate of the sonic point radius (15)Radius of the maximum of the disk angular momentum loss roughly corresponds to the sonic point radius (see Figs. 1–7). Since , we have (16)From the above equations we can see that adding cooling (higher p) can substantially increase the sonic point distance and consequently the angular momentum loss of the disk for a fixed Ṁ. For example, for p = 0.2 the maximum loss rate of angular momentum increases roughly 2 times and for p = 0.4 it increases roughly 5–6 times in comparison with the isothermal case, according to the type of star.
4. Numerical methods
4.1. Stationary calculation
For the stationary models (Kurfürst & Krtička 2012) we solved the system of Eqs. (1), (3), and (4)omitting the explicitly timedependent terms, using the NewtonRaphson method (e.g., Krtička 2003). For the purpose of this study we selected the star with the following parameters corresponding to the mainsequence spectral Type B0 (Harmanec 1988): T_{eff} = 30 000 K,M = 14.5 M_{⊙},R_{⋆} = 5.8 R_{⊙}. To solve the system of hydrodynamic equations (supplemented by the sonic point condition Eq. (14)) numerically we used the socalled shooting method, based on changing the inner boundary (photospheric) radial velocity in order to find a proper branch of the solution. The azimuthal velocity at the inner disk boundary (stellar equatorial surface) corresponds to Keplerian velocity. The solution of the system of hydrodynamic equations is independent of the scaling of the surface density Σ (cf. Krtička et al. 2011) and the mass loss rate Ṁ is treated as a free parameter in our calculations. In our models the inner boundary radius is the equatorial radius of the critically rotating star, R_{eq} = 3/2 R_{⋆}.
Numerical problems occured when we attempted to involve the second order viscosity term according to Eq. (9). Despite its complete analytic linearization in the Jacobi matrix, the solutions suffered from severe vibrations or perturbations mainly in the proximity of and above the sonic point, and so we used only the first order viscosity term Eq. (8)in the azimuthal momentum Eq. (4). For the numerical calculation we selected a radial grid consisting of 300–1000 grid points according to various initial conditions. We used the numerical package LAPACK (Anderson et al. 1999) to solve the system of linearized equations.
4.2. Timedependent calculation
For the timedependent calculations we write the lefthand sides of hydrodynamic equations in conservative form (see, e.g., Norman & Winkler 1986; Hirsch 1988; Stone & Norman 1992; Feldmeier 1995; LeVeque et al. 1998) (17)where the quantities u=Σ, ΣV, R × ΣV and F(u)=ΣV, ΣV ⊗ V, R × ΣV ⊗ V for the mass, momentum, and angular momentum conservation equations, respectively (with × denoting the vector product and ⊗ denoting the tensor product). Assuming axial symmetry of the thin disk, ∂/∂φ = 0, all functions are only radially and time dependent. The radial component of mass conservation equation is given in Eq. (1)and the radial component of momentum equation (see Eq. (3)) in its conservative form gives (18)where P is the isothermal gas pressure, P = a^{2}Σ. The explicit form of the angular momentum equation (see Eq. (10)) in this case is (19)the term denotes the density of the viscous force in a form derived in Eq. (9). Because we parameterize the disk temperature profile via Eq. (6), we do not employ the energy equation for the calculation in this case.
For the timedependent calculations we extended the onedimensional hydrodynamic code of Feldmeier (1995). Following Norman & Winkler (1986), the angular momentum advection flux acts as the azimuthal component of momentum flow. We nevertheless do not use the consistent advection schema (Norman et al. 1980) as it is described in detail in Norman & Winkler (1986), but employ Eq. (10)in its explicit form. Equations (1), (3), and (10)with use of Eq. (11)are discretized using timeexplicit operatorsplitting and finite difference method on staggered radial grids (see LeVeque et al. 1998). The advection fluxes are calculated on the boundaries of control volumes of these grids (see, e.g., Roache 1982; LeVeque et al. 2002) using van Leer’s monotonic interpolation (van Leer 1977, 1982).
In the source steps regarding the righthand sides of Eqs. (18)and (19)we accelerate the fluid by the action of external forces (gravity) and internal pressure forces on the gas momenta (see Norman & Winkler 1986). Involving the radial artificial viscosity term Q (see Eq. (23)) we solve finitedifference approximations to the following differential equations (see Norman & Winkler 1986), where Π = ΣV_{R} is the radial momentum density, J = ΣRV_{φ} is the angular momentum density, and is the second order viscosity term derived in Eq. (9). We adopt the artificial viscosity Q in the explicit form (Norman & Winkler 1986, see also, e.g., Caramana et al. 1998) (23)where V_{R} is the radial velocity component, a is the sound speed, the lower index i denotes the ith spatial grid step. The second term scaled by a constant C_{2} = 1.0 is the quadratic artificial viscosity (see Caramana et al. 1998) used in compressive zones. The linear viscosity term should be sparingly used for damping oscillations in stagnant regions of the flow (Norman & Winkler 1986). We use this term with C_{1} = 0.5 rarely when some oscillations may occur near the inner boundary region (near stellar equatorial surface), either in cases with low α_{0} viscosity parameter (α_{0}< 0.02) or in cases with steeper temperature decrease (p> 0.2).
For the timedependent modeling we employ two types of stars: main sequence stars of spectral Type B0 (Harmanec 1988) with parameters introduced in Sect. 5.1 and a Pop III star with the following parameters (Marigo et al. 2001; Ekström et al. 2008): T_{eff} = 30 000 K,M = 50 M_{⊙},R_{⋆} = 30 R_{⊙}. The calculations were extended up to a considerable distance from the parent star, although this may in most cases be a purely hypothetical issue due to the low disk density in such regions. According to the analytical prescriptions introduced in Sects. 2 and 3, we solve the set of Eqs. (1)–(4)with the use of Eq. (9).
The inner boundary (stellar equatorial) values were adopted on the equatorial radius of the critically rotating star (R_{eq} = 3/2 R_{⋆}) in following way. The estimation of inner boundary surface density Σ(R_{eq}) is implemented as a fixed boundary value; in the case of a critically rotating B0type star the isothermal disk surface density is Σ(R_{eq}) = 1.6 × 10^{2} g cm^{2}, roughly corresponding to Ṁ ≈ 10^{9} M_{⊙} yr^{1} derived by Granada et al. (2013). In the case of a critically rotating Pop III type star the isothermal disk surface density is Σ(R_{eq}) = 1.6 × 10^{5} g cm^{2}, which roughly corresponds to Ṁ ≈ 10^{6} M_{⊙} yr^{1} (Ekström et al. 2008). Similarly to the stationary calculations described in Sect. 4.1, the timedependent models are independent of the scaling of the surface density Σ. The inner boundary condition for V_{R} is free, i.e., the quantity is extrapolated from mesh interior values as a 0th order extrapolation. As the inner boundary condition for V_{φ} we assume the Keplerian velocity (critically rotating stellar equatorial surface). The outer boundary conditions are considered as outflowing for all quantities.
We set the initial surface density profile to Σ ~ R^{2} (Okazaki 2001). We start the numerical calculation with zero initial gas radial velocity and Keplerian rotational velocity throughout the whole disk. Our numerical test confirmed that the final solution is independent of the initial conditions and does not depend on the actual inner boundary condition for V_{R}.
5. Results of numerical models
5.1. Stationary calculations
In stationary calculations the first order linear viscosity was assumed according to Eq. (8). The model with constant viscosity profile and p = 0 in Fig. 1 shows a rapid decrease in the rotational velocity and the angular momentum loss at large radii. These quantities even drop to negative values in the case of the adopted first order linear viscosity prescription (see Sect. 3). The velocity drop is caused by the increase of the second term in Eq. (13)at large radii. We consider the drop to be unphysical. As a solution to this problem within stationary calculations we introduce the models with power law viscosity decline. Up to a certain value of p parameter in temperature power law profile the models avoid the rapid rotational velocity drop (show constant angular momentum loss rate) in supersonic region (Fig. 2).
5.2. Disk evolution time
In the timedependent models we recognize the wave that converges the initial state of calculated quantities to their final stationary state. Because of the initial density distribution we may regard the wave as physical (not only numerical artefact). We assume that during the disk developing phase a similar transforming wave occurs and its amplitude and velocity depends on physical conditions (namely the density distribution) in the stellar surroundings. There might be a possibility to observe some bow shocks at the boundary between the developing disk and the interstellar medium (even though the disk radial velocity in the distant regions is about one order of magnitude lower than in the case of linedriven stellar winds). The wave may also determine the timescale of the Be star disk growth and dissipation phases (e.g., Guinan & Hayes 1984; Štefl et al. 2003). In the subsonic region the wave establishes nearly hydrostatic equilibrium in the radial direction (Eq. (3)) and the wave speed approximately equals the sound speed. In the supersonic region this wave propagates as a shock wave. Its propagation speed can be approximated as (see Fig. 3), where the subscripts 0 and 1 denote the values in front of and behind the shock front, respectively (Zel’dovich & Raizer 1967). We regard the shock propagation time as the dynamical time t_{dyn} ≈ R/D = 0.3R/a. For example for the distance 10^{4}R_{eq} the isothermal constant viscosity B0type star disk model gives t_{dyn} ≈ 40 yr. The dynamical time is almost independent of the viscosity while it significantly increases with decreasing temperature.
Fig. 1 Dependence of the relative radial and azimuthal velocities and the relative angular momentum loss rate on radius for various temperature profiles calculated by the method described in Sect. 4.1. Constant viscosity α = 0.025 is assumed. Arrows denote the sonic point. 
Fig. 3 Snapshot of the radial velocity and the surface density transforming wave in the case of a B0type star isothermal constant viscosity model (see Fig. 5), the time t_{dyn} ≈ 40 yr. The wave propagation velocity is denoted as D. In the model, the ratio Σ_{1}/ Σ_{0} (surface densities behind and in front of the wavefront) is about one order of magnitude and slightly increases with the distance. 
Fig. 4 Comparison of the density wave propagation time (lower two branches denoted as t_{dyn}) with the disk viscous time (upper two branches denoted as t_{visc}), for B0type star isothermal (p = 0) constant viscosity model and the model with decreasing (p = 0.4) temperature profile (see Fig. 6), in dependence on radius. The disk viscous time is calculated from Eq. (24). Since the rotational velocity V_{φ} is adopted from the models, the graph of t_{visc} is cut off in the region of the rapid rotational velocity drop. The plotted values of t_{dyn} are adopted from the models. 
We associate the disk evolution time with disk viscous time (Okazaki 2001; Maeder 2009) (24)i.e., the timescale on which matter diffuses through the disk under the effect of the viscous torques (Frank et al. 1992). In the isothermal constant viscosity case the same model gives t_{visc} ~ 10^{2} yr for the sonic point radius. The viscous time significantly grows with temperature and viscosity. The comparison of the two times t_{dyn} and t_{visc} is in Fig. 4.
We also investigated the case with arbitrarily low nonzero initial surface density values in the whole computational domain. We assumed the initial (Keplerian) value of V_{φ} up to only a few tens of stellar radii followed by a discontinuous jump down to zero. Even using these initial conditions the disk evolves to a large distance and it converges to the proper final state. The density in the outer disk radius region is lower than the initial density value, forming a rarefaction wave that really extends radially with time (but more work on this point is needed).
5.3. Stationary state reached by timedependent models
From the calculations it follows that the profiles of surface density and radial velocity as well as the sonic point distance (where V_{R}/a = 1) very weakly depend on the viscosity parameter n. The outer limit of Keplerian rotation velocity region (V_{φ} ~ R^{0.5}) is almost independent of the viscosity parameter. The calculations nevertheless show strong dependence of the outer edge of the region where the rotational velocity behaves as angular momentum conserving (V_{φ} ~ R^{1}) (i.e., of the distance where the rotational velocity begins to rapidly drop) on viscosity parameter (for a given temperature profile). For a selected range of viscosity parameter n the distance of this region differs approximately by one order of magnitude (see, e.g., Fig. 5). In the models the location of this rapid rotational velocity drop does not exceed the radius where the disk equatorial density falls to averaged interstellar medium density (its mean density we assume as 10^{23} g cm^{3} e.g., Misiriotis et al. 2006; Maeder 2009). At this distance a kinetic plasma modeling would likely be required; moreover, the interaction of the disk with interstellar medium has to be taken into account.
Within timedependent calculations we examined the differences in the numerical results between the two different prescriptions for viscous torque (Eqs. (8)and (9)). Similarly to the steady disk calculations (see Sect. 5.1), the first order linear viscosity calculations (Eq. (8)) exhibit the rapid rotational velocity drop to negative values and consequently indicate very slow convergence to zero. The rotational velocity profiles calculated involving the second order linear viscosity term (Eq. (9)) confirm the analytical result V_{φ}> 0 throughout the entire disk range (see Sect. 3).
Fig. 5 Dependence of scaled vertically integrated density and radial and azimuthal disk velocities and the scaled angular momentum loss rate on radius in the case of isothermal disk (p = 0) of selected B0type star for various radial viscosity profiles (various n) in a final stationary state of timedependent models (the rapid drop of rotational velocity and the angular momentum loss rate in the outer disk region is a stationary jump, not a shock wave). 
Figure 5 illustrates the isothermal case (p = 0) of a B0type star (Ṁ ≈ 10^{9} M_{⊙} yr^{1}, see Sect. 4.2) with α(R_{eq}) = 0.025 (Penna et al. 2012). The calculated radius of the sonic point R_{s} ≈ 550 R_{eq} roughly corresponds to the analytical prediction from Eq. (15)with R_{s} being approximately 480 R_{eq}. The maximum angular momentum loss rate (see Eq. (16)) is independent of viscosity while it strongly depends on the profile of temperature. Since roughly equals the angular momentum loss rate at the sonic point radius, we assume the total angular momentum contained in the disk to be (25)and we therefore regard R_{s} as the disk outer edge. Analogously it also determines the mass of the disk (26)Comparing for example the total disk angular momentum J_{disk} with the total stellar angular momentum where a nondimensional parameter η = 0.05 (Meynet & Maeder 2006), in this case the ratio J_{disk}/J_{⋆} = 1.2 × 10^{6}.
Fig. 6 As in Fig. 5, however for temperature decreasing with a power law with p = 0.4. Inner boundary viscosity α(R_{eq}) = 0.025 is considered. The characteristic radii (sonic point distance, outer disk radius) are in this case significantly larger. 
Figure 6 shows the case of decreasing temperature profile (p = 0.4) of the B0type star with the same viscosity profiles as in Fig. 5. The sonic point distance is roughly R_{s} ≈ 31 500 R_{eq} for all the viscosity profiles, which is about two orders of magnitude larger than in the isothermal case. In this model the ratio J_{disk}/J_{⋆} = 7.9 × 10^{5}, which is similarly about two orders of magnitude larger than in the isothermal case.
Fig. 7 Comparison of the radial profiles of relative surface density and relative velocities and the radial profiles of angular momentum loss in the case of decreasing viscosity (n = 0.2) for isothermal disk (p = 0) and for outward decreasing temperature profile (p = 0.4) in a final stationary state of timedependent disk models for Pop III star. 
Figure 7 shows the calculated profiles of the pop III star’s disk (Ṁ ≈ 10^{6} M_{⊙} yr^{1}, see Sect. 4.2) with various temperature parameters for fixed decreasing viscosity (n = 0.2). The graph clearly shows a strong dependence of radii of the sonic point and of the rapid rotational velocity drop, as well as of the slopes of surface density and radial velocity on temperature profile. The sonic point radius is located at R_{s} ≈ 360 R_{eq} in the isothermal case and ≈16 500 R_{eq} in case of radially decreasing temperature with selected parameter p = 0.4. The use of the same dimensionless parameter η = 0.05 gives the ratio J_{disk}/J_{⋆} = 1.2 × 10^{3} for the isothermal model and J_{disk}/J_{⋆} = 0.2 for the model with decreasing temperature. In the latter case the disk carries away a significant fraction of stellar angular momentum and the star may not have enough angular momentum to develop the disk fully. In this case the stellar evolution has to be calculated together with the disk evolution. The calculations support the conclusion that the unphysical drop of the rotational velocity can be avoided in the models with radially decreasing viscosity parameter and temperature.
Within the timedependent calculations we also examined subcritically rotating stars modifying the boundary condition for V_{φ}. For the inner boundary value of the azimuthal velocity (where V_{K}(R_{eq}) denotes the Keplerian velocity at the stellar equator), the models precisely converge in the supersonic region. However, in the case when the boundary rotational velocity is only slightly higher than the above limit () there occur (more or less regular) pulsations in the density and radial velocity profiles in the region close to the star. For the density (and consequently the radial velocity) profile is unstable and gradually declines; for lower V_{φ}(R_{eq}) the decrease in density is faster.
6. Conclusions
We calculated axisymmetrical, vertically integrated onedimensional timedependent models of decretion disks of critically rotating stars. For this purpose we developed a numerical code for timedependent hydrodynamical modeling that includes full NavierStokes viscosity. We extrapolate the disk temperature profiles obtained by NLTE simulations (e.g., Carciofi & Bjorkman 2008) by the parameterized profiles. Various temperature profiles give different slopes of integrated density decrease throughout most of the disk. The radial dependence of the disk equatorial density may in various regions differ from the parameterized density profiles generally used in models dealing with the disk thermal structure (e.g., Sigut et al. 2009; McGill et al. 2013). Since the radial profile of the α viscosity parameter is not quite certain, we parameterize it via an independent power law radial dependence.
The timedependent onedimensional models confirm the basic results obtained in the stationary models in respect of the sonic point distance and of the distance of the disk outer edge (i.e., the radius where the rotational velocity begins to rapidly decrease) on parameterized temperature and viscosity profiles. The sonic point is located at larger radii in the models with steeper temperature decrease while its radius very weakly depends on the viscosity profile. The sonic radius strongly depends on both temperature and viscosity profiles and does not exceed the distance where we expect the disk equatorial density may drop to the average density of the interstellar medium. Consequently, the total angular momentum contained in the disk and the mass of the disk increase with the decreasing temperature and viscosity profiles. The analytical relations provided in Sect. 3 give adequate approximations of the numerical models. The unphysical drop of the rotational velocity and angular momentum loss at large radii, which is present in the isothermal models with constant viscosity parameters, can be avoided in the models with decreasing temperature and viscosity parameters.
Acknowledgments
The access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the program “Projects of Large Infrastructure for Research, Development, and Innovations” (LM2010005) is appreciated. This work was supported by the grant GA ČR 1310589S.
References
 Anderson, E., Bai, Z., Bischof, C., et al. 1999, LAPACK Users’ Guide, (Philadelphia: SIAM) [Google Scholar]
 Caramana, E. J., Shashkov, M. J., & Whalen, P. P. 1998, J. Comp. Phys., 144, 70 [Google Scholar]
 Carciofi, A. C., & Bjorkman, J. E. 2008, ApJ, 684, 1374 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Ruyter, S., van Winckel, H., Maas, T., et al. 2006, A&A, 448, 641 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ekström, S., Meynet, G., Maeder, A., & Barblan, F. 2008, A&A, 478, 467 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feldmeier, A. 1995, A&A, 299, 523 [NASA ADS] [Google Scholar]
 Frank, J., King, A., & Raine, D. 1992, Accretion power in astrophysics (Cambridge Astrophysical Series), vol. 21 [Google Scholar]
 Guinan, E. F., & Hayes, D. P. 1984, ApJ, 287, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Granada, A., Ekström, S., Georgy, C., et al. 2013, A&A, 553, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harmanec, P. 1988, Bull. Astron. Inst. Czech., 39, 329 [Google Scholar]
 Hirsch, Ch. 1988, Numerical Computation of Internal and External Flows (John Wiley & Sons Ltd.), 1 [Google Scholar]
 Krtička, J. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner (San Francisco: ASP), 259 [Google Scholar]
 Krtička, J., Owocki, S. P., & Meynet, G. 2011, A&A, 527, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kurfürst, P. 2012, in From Interacting Binaries to Exoplanets: Essential Modeling Tools, eds. M. T. Richards, & I. Hubeny (Cambridge: Cambridge University Press), IAU Symp., 282, 257 [Google Scholar]
 Kurfürst, P., & Krtička, J. 2012, in Circumstellar Dynamics at High Resolution, eds. A. Carciofi, & Th. Rivinius (San Francisco: ASP), 464, 223 [Google Scholar]
 Lee, U., Saio, H., & Osaki, Y. 1991, MNRAS, 250, 432 [NASA ADS] [CrossRef] [Google Scholar]
 LeVeque, R. J. 2002, Finite Volume Methods for Hyperbolic Problems (Cambridge University Press) [Google Scholar]
 LeVeque, R. J., Mihalas, D., Dorfi. E. A., & Müller, E. 1998, Computational Methods for Astrophysical Fluid Flow (Berlin: SpringerVerlag) [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars (Berlin Heidelberg: SpringerVerlag) [Google Scholar]
 Marigo, P., Girardi, L., Chiosi, C., & Wood, P. R. 2001, A&A, 371, 152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matsumoto, R. 1984, PASJ, 36, 71 [NASA ADS] [Google Scholar]
 Matt, S., Balick, B., Winglee, R., et al. 2000, ApJ, 545, 965 [NASA ADS] [CrossRef] [Google Scholar]
 McGill, M. A., Sigut, T. A. A., & Jones, C. E. 2013, ApJS, 204, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Meynet, G., & Maeder, A. 2006, Stars with the B[e] Phenomenon, eds. M. Kraus, & A. S. Miroschnichenko, ASPC Conf. Ser., 355, 27 [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundation of Radiation Hydrodynamics (New York: Oxford University Press) [Google Scholar]
 Millar, C. E., & Marlborough, J. M. 1998, ApJ, 494, 715 [NASA ADS] [CrossRef] [Google Scholar]
 Millar, C. E., & Marlborough, J. M. 1999, ApJ, 516, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Misiriotis, A., Xilouris, E. M., Papamastorakis, J., et al. 2006, A&A, 459, 113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Norman, M. L., & Winkler, K.H. A. 1986, Astrophysical Radiation Hydrodynamics. NATO Advanced Research Workshop on Astrophysical Radiation Hydrodynamics, held in Garching bei Munchen, Germany, August, 1982, eds. K.H. A. Winkler, & M. L. Norman [Google Scholar]
 Norman, M. L., Wilson, J. R., & Barton, R. T. 1980, ApJ, 239, 968 [NASA ADS] [CrossRef] [Google Scholar]
 Okazaki, A. T. 2001, PASJ, 53, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Penna, R. F., Sadowski, A., Kulkarni, A. K., & Narayan, R. 2012, MNRAS, 428, 2255 [Google Scholar]
 Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Roache, P. J. 1982, Computational Fluid Dynamics (New Mexico: Hermosa Publishers) [Google Scholar]
 Sigut, T. A. A., McGill, M. A., & Jones, C. E. 2009, ApJ, 699, 1973 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Smak, J. 1989, Acta Astron., 39, 201 [NASA ADS] [Google Scholar]
 Smith, N., & Townsend, R. H. D. 2007, ApJ, 666, 967 [NASA ADS] [CrossRef] [Google Scholar]
 Štefl, S., Baade, D., Rivinius, T., et al. 2003, A&A, 402, 253 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992, Ap&SS, 80, 753 [Google Scholar]
 van Leer, B. 1977, J. Comput. Phys., 23, 276 [Google Scholar]
 van Leer, B. 1982, Fluxvector splitting for the Euler equations, Proc. Int. Conf. on Numerical Methods in Fluid Dynamics, 8th, Aachen, West Germany, June 28–July 2 (Berlin: SpringerVerlag) [Google Scholar]
 Zel’dovich, Ya. B., & Raizer, Yu. P. 1967, Physics of shock waves and hightemperature hydrodynamic phenomena (New York: Academic Press) [Google Scholar]
All Figures
Fig. 1 Dependence of the relative radial and azimuthal velocities and the relative angular momentum loss rate on radius for various temperature profiles calculated by the method described in Sect. 4.1. Constant viscosity α = 0.025 is assumed. Arrows denote the sonic point. 

In the text 
Fig. 2 As in Fig. 1, but with variable α parameter, α ~ R^{0.2}. 

In the text 
Fig. 3 Snapshot of the radial velocity and the surface density transforming wave in the case of a B0type star isothermal constant viscosity model (see Fig. 5), the time t_{dyn} ≈ 40 yr. The wave propagation velocity is denoted as D. In the model, the ratio Σ_{1}/ Σ_{0} (surface densities behind and in front of the wavefront) is about one order of magnitude and slightly increases with the distance. 

In the text 
Fig. 4 Comparison of the density wave propagation time (lower two branches denoted as t_{dyn}) with the disk viscous time (upper two branches denoted as t_{visc}), for B0type star isothermal (p = 0) constant viscosity model and the model with decreasing (p = 0.4) temperature profile (see Fig. 6), in dependence on radius. The disk viscous time is calculated from Eq. (24). Since the rotational velocity V_{φ} is adopted from the models, the graph of t_{visc} is cut off in the region of the rapid rotational velocity drop. The plotted values of t_{dyn} are adopted from the models. 

In the text 
Fig. 5 Dependence of scaled vertically integrated density and radial and azimuthal disk velocities and the scaled angular momentum loss rate on radius in the case of isothermal disk (p = 0) of selected B0type star for various radial viscosity profiles (various n) in a final stationary state of timedependent models (the rapid drop of rotational velocity and the angular momentum loss rate in the outer disk region is a stationary jump, not a shock wave). 

In the text 
Fig. 6 As in Fig. 5, however for temperature decreasing with a power law with p = 0.4. Inner boundary viscosity α(R_{eq}) = 0.025 is considered. The characteristic radii (sonic point distance, outer disk radius) are in this case significantly larger. 

In the text 
Fig. 7 Comparison of the radial profiles of relative surface density and relative velocities and the radial profiles of angular momentum loss in the case of decreasing viscosity (n = 0.2) for isothermal disk (p = 0) and for outward decreasing temperature profile (p = 0.4) in a final stationary state of timedependent disk models for Pop III star. 

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