Wave mediated angular momentum transport in astrophysical boundary layers
AAstronomy & Astrophysics manuscript no. aa26005-15 c (cid:13)
ESO 2018November 14, 2018
Wave mediated angular momentum transport in astrophysicalboundary layers
Marius Hertfelder (cid:63) and Wilhelm Kley Institut für Astronomie und Astrophysik, Abt. Computational Physics, Auf der Morgenstelle 10, 72076 Tübingen, GermanyReceived xx / Accepted xx
ABSTRACT
Context.
Disk accretion onto weakly magnetized stars leads to the formation of a boundary layer (BL) where the gas loses its excesskinetic energy and settles onto the star. There are still many open questions concerning the BL, for instance the transport of angularmomentum (AM) or the vertical structure.
Aims.
It is the aim of this work to investigate the AM transport in the BL where the magneto-rotational instability (MRI) is notoperating owing to the increasing angular velocity Ω ( r ) with radius. We will therefore search for an appropriate mechanism andexamine its e ffi ciency and implications. Methods.
We perform 2D numerical hydrodynamical simulations in a cylindrical coordinate system ( r , ϕ ) for a thin, vertically inte-grated accretion disk around a young star. We employ a realistic equation of state and include both cooling from the disk surfaces andradiation transport in radial and azimuthal direction. The viscosity in the disk is treated by the α -model; in the BL there is no viscosityterm included. Results.
We find that our setup is unstable to the sonic instability which sets in shortly after the simulations have been started. Acousticwaves are generated and traverse the domain, developing weak shocks in the vicinity of the BL. Furthermore, the system undergoesrecurrent outbursts where the activity in the disk increases strongly. The instability and the waves do not die out for over 2000 orbits.
Conclusions.
There is indeed a purely hydrodynamical mechanism that enables AM transport in the BL. It is e ffi cient and wavemediated; however, this renders it a non-local transport method, which means that models of a e ff ective local viscosity like the α -viscosity are probably not applicable in the BL. A variety of further implications of the non-local AM transport are discussed. Key words. accretion, accretion disks – hydrodynamics – instabilities – waves – methods: numerical – stars: protostars
1. Introduction
One ubiquitous phenomenon in astrophysics is the accretion ofmatter on a central object via an accretion disk. This process canbe observed for a variety of central objects, such as young stars,compact objects like white dwarfs or neutron stars, and activegalactic nuclei. The physics of the accretion disk itself is reason-ably well understood. However, the tiny region where the accre-tion disk connects to the star is still one of the major problems. Inthis region, called the boundary layer (BL), a smooth connectionis established between the disk, which rotates nearly at Keple-rian frequency, and the star, which in general rotates much moreslowly. It is of great dynamical and thermal importance since thegas undergoes a supersonic velocity drop in a very confined re-gion of a few percent of the stellar radius and loses up to onehalf of the total accretion energy during that process (Klu´zniak1987).The BL has been studied extensively for over 40 years, bothin an analytical (e.g. Bertout & Regev 1992) and a numericalapproach (e.g. Kley & Hensler 1987). Simulations involving theBL around a young star have been performed by Kley & Lin(1996) and Kley & Lin (1999), including full radiation transport.Most of the work that has been done on the BL (see e.g. Hert-felder et al. 2013, for a brief review of the history) assumes thatthe gas first slows down in the equatorial plane of the disk andthen spreads around the star. This is the classical picture of theBL. However, there is an alternative theory called the spread- (cid:63) e-mail: [email protected] ing layer (Inogamov & Sunyaev 1999, 2010), which was ini-tially developed for neutron stars. Within this concept the gasfirst spreads around the star because of the ram pressure in theBL and then slows down on the whole surface of the star.The majority of the simulations of the BL start from thepremise that the observed angular momentum (hereafter AM)transport in the BL is driven by local turbulent stresses andthe authors consequently adopt a α -viscosity model (Shakura &Sunyaev 1973) in their simulations. Sometimes the classical α -model is modified slightly in order to prevent supersonic infallvelocities (Kley & Papaloizou 1997) or, for instance, to take intoaccount the radial pressure scale height in the BL (Popham &Narayan 1995). The utilization of a local viscosity model waslater justified for the disk by the discovery of the magnetorota-tional instability (MRI) which creates turbulence that acts like agenuine viscosity on macroscopic scales (Velikhov 1959; Chan-drasekhar 1960; Balbus & Hawley 1991, 1998; Balbus 2003;Balbus & Lesa ff re 2008). However, as pointed out by Godon(1995) and Abramowicz et al. (1996) and recently shown by Pes-sah & Chan (2012), if the angular velocity increases with radius,d Ω / d r >
0, the MRI is e ff ectively damped out and the asso-ciated AM transport oscillates about zero. Since this situationclearly applies for the BL, we do not expect to obtain su ffi cientAM transport through the MRI. There have been various alter-native transport mechanisms proposed, among them the Kelvin-Helmholtz instability (Kippenhahn & Thomas 1978), the baro-clinic instability (Fujimoto 1993), and the Tayler-Spruit dynamo(Tayler 1973; Spruit 2002; Piro & Bildsten 2004), none of which Article number, page 1 of 18 a r X i v : . [ a s t r o - ph . S R ] M a y & A proofs: manuscript no. aa26005-15 have yet been proven to e ffi ciently transport mass and AM in theBL.In this work, we will focus on the AM transport in the BLthrough non-axisymmetric instabilities and therefore carry out2D simulations in the midplane of the star-disk system usingcylindrical coordinates. Recently, a promising candidate for thetransport has been proposed and investigated in a series of pa-pers. According to this theory, the steep velocity drop in theBL is prone to the sonic instability (Glatzel 1988; Belyaev &Rafikov 2012), which is an instability of a supersonic shearlayer, much like the Papaloizou-Pringle instability (Papaloizou& Pringle 1984; Narayan et al. 1987). Acoustic waves are ex-cited in the BL as a consequence of the sonic instability andAM can be transported by these modes in an e ffi cient way. Thishas been demonstrated for 2D flows (Belyaev et al. 2012), 3Dflows in cylindrical coordinates (Belyaev et al. 2013a), and evenfor 3D magnetohydrodynamical flows (Belyaev et al. 2013b).The wave mediated AM transport implies that this process isintrinsically non-local since the waves can potentially travel along way before they dissipate and release the AM to the fluid.Therefore, it is problematic to describe the AM transport in theBL by means of a local viscosity like the α -model. Although theabove-mentioned simulations are already quite sophisticated, theauthors make some simplifications that constrain the validity oftheir models. This is the point where we step in with this paperand relax three of the simplifications made. First, we make useof a realistic equation of state instead of an isothermal one andpropagate the temperature of the system, as well. We use full ra-diative di ff usion in the disk plane and an approximation for thevertical flux, and hence employ a quasi-3D radiation transport.Second, we use realistic, state-of-the-art initial models like thosealready published in Hertfelder et al. (2013) as a starting pointfor the simulations presented here. Third, mass constantly entersthe simulation domain from the outer boundary and thus there isa net mass flow through the disk, as is the case for a real accre-tion disk. Because of the high computational costs the treatmentof the radiation brings with it, we perform 2D simulations anddo not yet include magnetic fields. Thus, the next steps will be toextend the simulations presented here to 3D and also to includemagnetic fields.The paper is organized as follows. In Sect. 2, the basic phys-ical model is described. We give an overview of the equationsused for our simulations and present the assumptions we havemade. Furthermore, we briefly review how the AM transport dueto stresses in the fluid can be measured and quantized. Section 3is devoted to the numerical method that we employ in our simu-lations. We discuss the numerical code, the initial and boundaryconditions, as well as the model parameters. In Sect. 4 we startthe analysis of the results. We first discuss the sonic instabil-ity and the acoustic modes that are excited in the BL. Then, inSect. 5, we investigate how these instabilities a ff ect the physi-cal picture of the BL and look deeply into the AM transport andlong-term behavior of our models. We conclude with Sect. 5.
2. Physics
In this section, we present the physical foundations for the sim-ulations we have performed and that are described later in thispublication. In order to investigate non-axisymmetric behaviorin the disk plane, we utilized the vertically integrated Navier-Stokes equations. Since the focus of our current work lies oninstability developing in the BL, we describe the flow of a per-turbed, compressible fluid and also the derivation of the resultingReynolds stresses.
In order to obtain a set of 2D Navier-Stokes equations includ-ing radiation transport, we apply a cylindrical coordinate systemwith coordinates ( r , ϕ, z ) and integrate the 3D equations over thevertical coordinate z . The star lies in the center of the coordinatesystem and the midplane of the disk coincides with z =
0. If weassume a Gaussian profile for the mass density ρ , we can definethe 2D surface density by Σ = (cid:90) ∞−∞ ρ d z = √ π ρ ( r , ϕ, H , (1)where we introduce the e ff ective pressure scale height H , whichis a measure of the vertical extent of the disk. Under the assump-tion of hydrostatic balance and an isothermal equation of state inthe z -direction, the pressure scale height reads H = c s Ω K , (2)with c s and Ω K = (cid:112) GM ∗ / r being the midplane soundspeedand the Keplerian angular velocity, respectively. The inverse of Ω K ( R ∗ ) is proportional to the period of a Keplerian orbit at thesurface of the star, P ∗ = π/ Ω K ( R ∗ ). Derivatives with respect to z are neglected, as is the z -component of the velocity, u z . Clearly,this is not a reasonable assumption if one wants to investigate theoverall structure of the BL. However, it is still justified since werestrict ourselves to the study of non-axisymmetric flow structurein this work.The vertically integrated continuity equation reads ∂ Σ ∂ t + r ∂ ( r Σ u r ) ∂ r + r ∂ ( Σ u ϕ ) ∂ϕ = . (3)We express the momenta equations in terms of the conservativevariables and introduce the radial momentum density s = Σ u r and the angular momentum density h = Σ ru ϕ ≡ Σ r Ω . Thesevariables reflect the physically conserved quantities, as opposedto the primitive variables u r and u ϕ . The conservation of momen-tum in radial direction is then given by ∂ s ∂ t + r ∂ ( rsu r ) ∂ r + r ∂ ( su ϕ ) ∂ϕ − Σ u ϕ r = − ∂ p ∂ r + r ∂ ( r σ rr ) ∂ r + r ∂ ( σ r ϕ ) ∂ϕ − r σ ϕϕ − Σ ∂ Φ ∂ r , (4)and the equation in azimuthal direction reads ∂ h ∂ t + r ∂ ( rhu r ) ∂ r + r ∂ ( hu ϕ ) ∂ϕ = − ∂ p ∂ϕ + r ∂ ( r σ r ϕ ) ∂ r + ∂ ( σ ϕϕ ) ∂ϕ − Σ ∂ Φ ∂ϕ . (5)Here p denotes the vertically integrated pressure p = Σ R G T /µ ,where R G = k B / m H with Boltzmann’s constant k B and the massof hydrogen m H , T is the temperature, and µ is the mean molec-ular weight; Φ = − GM ∗ / r is the gravitational potential ( G and M ∗ are the gravitational constant and stellar mass). The compo-nents of the vertically integrated viscous stress tensor, σ i j , can befound in e.g. Masset (2002). The equation for the conservation Article number, page 2 of 18. Hertfelder and W. Kley: Wave mediated AM transport in astrophysical BLs of energy is given by ∂ e ∂ t + r ∂ ( reu r ) ∂ r + r ∂ ( eu ϕ ) ∂ϕ = − p (cid:34) r ∂ ru r ∂ r + r ∂ u ϕ ∂ϕ (cid:35) + σ rr ∂ u r ∂ r + σ r ϕ (cid:34) r ∂ u r ∂ϕ + r ∂ Ω ∂ r (cid:35) + σ ϕϕ (cid:34) ∂ Ω ∂ϕ + u r r (cid:35) − σ SB T ff − Hr (cid:34) ∂ ( rF r ) ∂ r + ∂ F ϕ ∂ϕ (cid:35) , (6)where σ SB and T e ff are the Stefan-Boltzmann constant and thee ff ective temperature, and the factor H follows from the verticalintegration.The last line of Eq. (6) incorporates the radiative flux in theflux-limited di ff usion approximation, F = − K R ∇ T . (7)Utilizing this approach, we assume that the BL is optically thickboth in the radial and vertical directions, where the inclusionof a flux-limiter also allows for the treatment of optically thinregions. Equation (7) has exactly the same mathematical formas molecular heat conduction in a gas. The e ff ective radiativeconductivity K R is K R = λ acT κ R ρ , (8)where κ R , a , and c are the Rosseland mean opacity, the radiationconstant and the speed of light, and λ is the flux limiter (Lever-more & Pomraning 1981; Levermore 1984). We adopt the for-mulation given by Kley (1989) for λ . We assume that the disklocally radiates like a blackbody of temperature T e ff and asso-ciate the vertical flux with the blackbody luminosity divided bythe radiating area, F z = σ SB T ff . The factor 2 in Eq. (6) stemsfrom the fact that the disk has two sides. The vertical flux pro-vides the cooling of the disk, while the other parts of the energyequation generate heat through viscous dissipation (second linein Eq. 6) or redistribute it in the disk plane (the terms with F r and F ϕ ). The e ff ective temperature T e ff is evaluated according toHubeny (1990), who approximates the optical depth in the ver-tical direction of the disk (see Sect. 5.3). By approximating theflux in the vertical direction as well, we employ a quasi-3D radi-ation transport.We also adopt a piecewise prescription for the Rosselandmean opacity, where in each temperature regime κ R is given bya power-law approximation, κ R = κ ρ a T b , as described in Lin &Papaloizou (1985) and Bell & Lin (1994). Each of the regimesis characterized by a prevailing mechanism that dominates theopacity in this region. The parameters κ , a , and b for the variousopacity regimes are listed in e.g. Müller & Kley (2012). In the simulations presented here, we strictly distinguish be-tween the BL and the disk with regard to viscosity. As has al-ready been mentioned in Sect. 1, the AM transport in the diskis mainly driven by turbulent stresses due to the MRI and cantherefore be parametrized with the classic α -viscosity approachby Shakura & Sunyaev (1973). The e ff ective kinematic viscosityis written as ν = α c s H (9)in this ansatz, where α is a dimensionless parameter whose valueis derived from MRI simulations. The cause of the AM in the BL and its investigation is, however, the subject of this work. Con-sequently, we do not make use of an additional viscosity term inthe BL region. The entire AM transport in the BL is thereforedue solely to the Reynolds stresses exerted by the gas.From a numerical point of view, the distinction between theBL and the disk is made by searching for the global maximumof the radial angular velocity profile, Ω ( r max ) = Ω max . This hasto be done for every angle ϕ in order to prevent a suppressionof non-axisymmetric perturbations. We can then computation-ally realize the viscosity as pointed out above by utilizing thefollowing conditions: – If r < r max , we set α = BL region ). – If r > r max , α is set to the default value, α = .
01, and theviscous stress tensor does not vanish ( disk region ). Gas flows in general can be characterized by a laminar bulk flowwith superimposed perturbations. One typical approach is to de-compose the flow variables f into an averaged or filtered part ¯ f and a fluctuating or subfiltered, unresolved part f (cid:48) : f ( r , t ) = ¯ f ( r , t ) + f (cid:48) ( r , t ) . (10)Within the classical Reynolds approach, the averaging is per-formed over an ensemble of systems, thus ¯ f is an ensemblemean. In the case of a statistically steady or stationary flow field,or when it is homogeneous, ergodicity can be assumed. Thismeans that both a su ffi ciently long time average over many timescales and a spatial average over many length scales correspondto the ensemble mean. Often, practical flows are inhomogeneousand the ensemble mean is replaced by a time average:¯ f τ ( r , t ; τ ) = τ (cid:90) tt − τ f ( r , t (cid:48) ) d t (cid:48) . (11)If, on the other hand, there is a symmetry in at least one direction,as there is in our case, a spatial average is also a viable solution.The actual average process in the Reynolds decomposition isa simple arithmetic mean over space or time (see Eq. 11). A morecommon approach in the case of compressible flows, however, isto use a density-weighted Reynolds average for the filtered partof the variables (Gatski & Bonnet 2013; Peyret 1996; Balbus& Papaloizou 1999). These density-weighted, or Favre variablesare given by a decomposition, which reads f = ˜ f + f (cid:48)(cid:48) with ˜ f = ρ f ¯ ρ , (12)where the bar denotes a Reynolds average and the tilde a Favreaverage.After performing a decomposition (of any kind) of the flowfield, one can then derive the Navier-Stokes equations for themean or resolved field. They are called the Reynolds-averagedNavier-Stokes (RANS) equations, or, in the case of a Favremean, the Favre-averaged Navier-Stokes (FANS) equations (e.g.Gatski & Bonnet 2013; Peyret 1996). Those mean field equa-tions contain additional terms, however, that are dependent onthe fluctuating part of the flow variables. In many engineeringapplications where simulations of the mean field are performed,those parts are unknown and have to be dealt with by using ap-propriate closure relations or models. Here, on the other hand,we are running direct numerical simulations and solve for the Article number, page 3 of 18 & A proofs: manuscript no. aa26005-15 actual values of the flow field variables. From these simulations,we can compute the additional terms, which appear in the aver-aged equations.From the FANS equations, we obtain the Reynolds stressesexerted by the fluctuating motion of the gas. For accretion disks,the dominant contribution is given by R r ϕ = Σ u r u ϕ − Σ ˜ u r ˜ u ϕ , (13)where we have already performed a vertical integration. Thequantity R r ϕ appears in the angular momentum equation (5),where it adds to σ r ϕ . We note that for the second part in Eq. (13),it is incorrect to use either the mean values of the laminar flow(i.e. the initial values, when starting from a relaxed model) orassume ˜ u r ≈
0. Our simulations showed that the instability-dominated state has other mean values than the laminar flowstate, and that both terms in Eq. (13) are approximately of thesame order of magnitude. Our simulations also showed that forthe quasi-stationary state both temporal and azimuthal averagingprocedures yield the same mean values since the evolution of theBL is significantly slower than an orbital period. Therefore, wemostly used azimuthal averaging for the mean values since it iscomputationally far less expensive.The vertically integrated Reynolds stress, which has the di-mensions of a 2D pressure, can be expressed by a dimensionlessparameter α Re via α Re = R r ϕ Σ c . (14)Comparing the strength of the Reynolds stresses and the viscousstresses, we can identify R r ϕ = − σ r ϕ and calculate an e ff ectiveviscous α ν from this equation. This relation is given by α ν = − R r ϕ Σ c s (cid:96) r ∂ Ω ∂ r , (15)where (cid:96) is the turbulent scale height in the BL, which is approx-imated to be the smaller of the disk scale height H and the radialpressure scale height (see Papaloizou & Stanley 1986; Popham& Narayan 1995). The mass and AM transport in the disk are two important quan-tities when we are interested in the e ff ectivity of the Reynoldsstresses in the BL. Both result from the 1D stationary hydrody-namic equations, where they play the role of eigenvalues. Themass flux is obtained from the continuity equation and reads˙ M = π r Σ u r . (16)The angular momentum flux follows from the momentum equa-tion in azimuthal direction and is given by˙ J = π r Σ u r u ϕ − π r σ r ϕ + π r R r ϕ . (17)Most commonly, the AM flux ˙ J is expressed as a dimensionlessparameter j , which is the AM flux normalized to the advectiveAM flux at the surface of the star: j = ˙ J ˙ J ∗ , ˙ J ∗ = ˙ MR ∗ Ω K ( R ∗ ) . (18)Either flux is defined such that if mass or AM are traveling out-ward, ˙ M and ˙ J are positive. The three terms in Eq. (17) can be identified as the advective AM transport, the viscous transport,and the transport due to Reynolds stresses. Both ˙ M and j areconstant with respect to space and time only if the system hasreached a stationary state, i.e. ∂/∂ t =
0. Then, ˙ M equals theimposed mass flux that is set as a boundary condition. For non-stationary states, ˙ M and j represent the local mass and AM fluxat a certain time.
3. Numerics
We used the code
FARGO (Masset 2000) with modifications ofBaruteau (2008) for the simulations presented in this paper.
FARGO is a specialized hydrodynamics code for simulations ofan accretion disk in cylindrical coordinates ( r , ϕ ) that exploitsthe advantages of the FARGO algorithm (Masset 2000). TheFARGO algorithm is a special method for the azimuthal trans-port in di ff erentially rotating disks that is able to speed up thesimulations by up to one order of magnitude and exhibits alower numerical viscosity than the usual transport algorithm.The speed-up works best when the velocity perturbations aresmaller than the bulk flow in azimuthal direction. Since the aver-age azimuthal velocity u ϕ ( r ) changes significantly in the courseof the simulations, the background rotational profile for theFARGO algorithm is updated at each time step in order to gainfull benefit of the method. The algorithm has also been shown toproduce valid results when shocks are created in the disk (Masset2000). The radiative di ff usion in the flux-limited approximation(Eq. 7) is solved implicitly using a successive over-relaxation(SOR) method (see Müller & Kley 2013). The code has beentested and used extensively for a variety of applications. In orderto verify its suitability for simulations of the BL, we comparedit to the 1D BL code described in Hertfelder et al. (2013) andfound a perfect agreement for the 1D models. To model the BL with a net mass flow through the disk wehave implemented specific boundary conditions in
FARGO . Atthe outer boundary we impose a mass flux ˙ M which goes intothe disk and whose absolute value can be chosen as a parameter.We require the azimuthal velocity to be Keplerian at the outerboundary, u ϕ ( r out ) = r out Ω K ( r out ). Modeling the inner bound-ary, i.e. the beginning of the star, is more di ffi cult. Since we aretreating the problem in cylindrical geometry in the disk plane,we cannot simultaneously simulate the outer parts of the star.Therefore, we start the simulation domain of all our simulationsat r in = R ∗ , where R ∗ is the radius of the chosen star, and setthe azimuthal velocity to the velocity of the rigidly rotating star, u ϕ ( r in ) = R ∗ Ω ∗ , i.e. we impose a no-slip boundary condition in ϕ -direction. The radial velocity is set to a certain fraction F ofthe Keplerian velocity u r ( r in ) = −F r in Ω K ( r in ) . (19)This is obviously an outflow BC for the radial velocity at theinner boundary and it is necessary since the matter which is in-serted in the domain from the outer boundary must be allowedto leave again at the inner boundary. The factor F determineshow deep the inner boundary r in is located inside the star: thesmaller F is, the deeper the domain ranges into the star. The ac-tual choice of F is a trade-o ff between containing a su ffi cientpart of the stellar envelope in the domain and avoiding densities Article number, page 4 of 18. Hertfelder and W. Kley: Wave mediated AM transport in astrophysical BLs and temperatures that are too high at the inner boundary to pre-vent the simulations from being slowed down too much. A valueof F = − is used in our simulations. We have tested di ff er-ent values of F and found that it does not have any influenceon the initial models apart from shifting the radial profiles by asmall amount in r . In particular, it does not have an e ff ect on themass flux through the inner boundary once the initial model is inan equilibrium state. Instead, the density adjusts such that ˙ M atthe inner boundary equals the value at the outer boundary. Apartfrom the special boundary treatment mentioned above, we im-pose zero-gradient or Neuman-type boundary conditions for theremaining variables, which means that the normal derivative atthe boundary vanishes.The 2D simulations were started from fully relaxed 1D ax-isymmetric BL models which were generated using the samemethods as in Hertfelder et al. (2013). Those simulations usea classical α -parametrization for the viscosity with a uniformvalue of α = .
01 both in the disk and the BL. We employedthese models to have realistic profiles for the density and tem-peratures as initial conditions. After the 1D models reached astationary state we extended them to 2D input data for
FARGO andran them again with a small number of azimuthal cells until theywere fully relaxed. Following this burn-in procedure we startedthe actual simulations from these axisymmetric initial conditionswith full resolution in r and ϕ . Prior to the start, the radial veloc-ity was randomly perturbed by 1% of the Keplerian velocity atthe surface of the star. The initial profiles are shown in Fig. 10below (dashed lines). In this paper we focus on the non-axisymmetric instabilities andthe AM transport in the BL. We chose the physical system of aBL around a young star, because the ratio of the mass and theradius of the star, M ∗ / R ∗ , is much smaller than, for instance, fora white dwarf, which we considered in a previous paper aboutthe BL. This ratio determines the radial scale height in the BLand the smaller H is, the higher the resolution must be in orderto resolve the BL correctly. Therefore, the computational costsincrease enormously when one goes to more compact objects.However, we expect what we find here to be the generic case.The young star we consider in this paper has one solar mass, M ∗ = M (cid:12) , and is three times the radius of the sun, R ∗ = R (cid:12) .It has an e ff ective temperature of T ∗ = L (cid:12) . The star is not rotating, i.e. Ω ∗ =
0, which isslightly below the observations (Cohen & Kuhi 1979), but waschosen for the sake of a large velocity drop. We impose a massaccretion rate of ˙ M = − M (cid:12) / yr, which is a typical value for TTauri stars such as HL Tau (Hayashi et al. 1993; Lin et al. 1994),and is two orders of magnitude smaller than the accretion rateFU Orionis stars can reach in outburst phases (Hartmann et al.1993). We took α = .
01 for the viscosity parameter in the diskand α = – Simulation A: It has a resolution of 512 ×
512 and the radialdomain ranges from R ∗ to 2 R ∗ . This simulation features thelowest resolution of all runs performed. – Simulation B: This is the reference simulation, and unlessstated otherwise, we always refer to it in the text. It has asimulation domain of [ R ∗ , R ∗ ] and a resolution of 1024 × – Simulation C: Here, the number of grid cells is doubled ineach direction and hence the resolution is given by 2048 × – Simulation D: In order to ensure that the radial domain insimulations A-C is not too small, we increased it to [ R ∗ , R ∗ ].While the first three simulations utilized an arithmetic grid,we chose a logarithmic spacing in radius for this run alongwith a resolution of 1024 ×
4. Boundary layer instability
In this section, we will summarize the key features of the sonicinstability and its associated acoustic waves in the BL (Glatzel1988; Belyaev & Rafikov 2012; Belyaev et al. 2012, 2013a) anddemonstrate that our models are also prone to the sonic instabil-ity.
The sonic instability arises in supersonic shear flows and is anon-local instability similar to the Papaloizou-Pringle instability(Papaloizou & Pringle 1984; Narayan et al. 1987; Glatzel 1988).It is the analogue of the classical KH instability in the regime of asupersonic drop in the velocity of a compressible, stratified fluid.The sonic instability has recently been studied in the context ofthe BL extensively by Belyaev & Rafikov (2012).There are two ways by which the sonic instability operates:The first is the overreflection of sound waves from a critical layerthat corresponds to the corotation resonance in a rotating fluid,i.e. the radius where the pattern frequency of the mode matchesthe rotation frequency of the bulk flow. This mechanism resultsin a dispersion relation that features several sharp peaks andhence discrete modes are clearly favored, especially if the den-sity contrast across the shear interface is high (see e.g. panel (n)of Fig. 3 in Belyaev & Rafikov 2012). Each mode can be associ-ated with a pseudo-energy, or conserved action (angular momen-tum), which can leak through the critical layer when a mode isreflected at the corotation resonance (Narayan et al. 1987). Thetunneled wave has the opposite sign of the action and, since it is aconserved quantity, the reflected wave must undergo an increasein its action in order to globally compensate for the additionalnegative action. Thus, the reflected wave is amplified. Further-more, if a mode is trapped such that it performs multiple reflec-tions at the corotation resonance, this amplification mechanism,called the corotation amplifier (e.g. Mark 1976), eventually leadsto instability.The other destabilizing mechanism is the radiation of en-ergy away from the BL. If a mode inside the BL has a negativeaction density, it becomes even more negative through the radi-ation mechanism and, consequently, is amplified, which resultsin further instability. The dispersion relation of this destabilizingmechanism is a smooth function of the wave vector. Therefore,a continuum of modes are associated with it. This mechanismdominates the overreflection for a small density contrast on ei-ther side of the interface (see again Fig. 3 of Belyaev & Rafikov2012).Each mechanism, the overreflection and the radiation, ismost e ffi cient in the regime of large density contrast or equaldensity, respectively, and they operate together. The growth ratesof the sonic instability are approximately 1 / P ∗ , and the depen-dence of the growth rate on the density contrast is small, de-spite the di ff erent destabilizing mechanisms. This further dis-tinguishes the sonic instability from the KH instability. In the Article number, page 5 of 18 & A proofs: manuscript no. aa26005-15 r a d i a l k i n e t i c E n e r g y [ e r g ] simulation Asimulation Bsimulation Crate: 0.605rate: 1.087rate: 1.090 Fig. 1.
The total radial kinetic energy as a function of the very firstorbits for three simulations mentioned in Sect. 3.3, which di ff er only intheir spatial resolution. At around orbit 5, the sonic instability starts tooperate in all three cases, and as a result, the radial kinetic energy risesrapidly. The growth rates were obtained by applying an exponential fitto the first exponential increase. numerical simulations of the BL, a very high numerical res-olution is necessary to correctly reproduce the growth rate ofthe sonic instability. Such a high resolution is computationallyvery demanding and costly, especially for investigations of thelong-term evolution of the BL. However, as has been shown byBelyaev & Rafikov (2012) and revisited in Belyaev et al. (2012),a lower resolution merely leads to a slower growth of the insta-bility.We confirm this point by comparing the growth rates of threesimulations with di ff erent spatial resolution. Figure 1 shows theevolution of the radial kinetic energy when the sonic instabilitystarts to set in for the three resolutions 512 × × × ∼ . and 2048 cells both yieldan almost identical growth rate of 1 .
09 per orbit. The shear layerwith positive gradient of Ω is initially resolved with approxi-mately 150 and 300 cells in radial direction, respectively. In sim-ulation A, we only have about 75 cells across the velocity drop.We also verify by Fig. 1 that the growth rates of the sonic insta-bility, which are of the order of 1 per orbit, are actually very high.We will now focus again on our reference simulation, althoughall simulations behave basically the same. After about ten orbits,the rise of the radial kinetic energy is interrupted briefly and thencontinues with a slightly di ff erent growth rate. The paused in-crease in the instability and the change in growth rate is due tothe change in density and velocity in the shearing region withtime. Since the mass inside the BL is delivered from the disk,the density in the vicinity of the velocity drop rises. The veloc-ity itself changes as a result of the initiated angular momentumtransport. As has already been mentioned and is elucidated inBelyaev & Rafikov (2012), the growth rate of the sonic instabil-ity shows a weak dependence on the density ratio and the modesdepend on the slope of the velocity drop, as well. Thus, a changein the instability mechanism can be observed at this time.We now have a closer look at how the sonic instability man-ifests itself in the radial velocity and consider Fig. 2 for this pur-pose. The upper panel is a snapshot taken after five orbits and the lower panel shows u r at t =
10 orbits. In both panels, the veloc-ity is normalized to the Keplerian velocity at the surface of thestar. Figure 2a) features a distinct boundary at r crit ≈ .
19, whichcoincides with the inflection point in Ω ( r ) and we can clearly seesound waves on either side, propagating away from the bound-ary. The pattern speed equals the rotation rate of the unperturbedflow at the location of the boundary and thus we associate it withthe critical layer, or corotation resonance, which we have alreadyintroduced. It is also clearly noticeable that the modes undergoa phase shift across the critical layer. Both features are a clearsignature of the sonic instability and we can compare Fig. 2a) toFig. 2 of Belyaev et al. (2012) and Fig. 5a) of Belyaev & Rafikov(2012) to verify that our model is indeed unstable to the sonic in-stability.Figure 2b) depicts the radial velocity five orbits later and theappearance has already noticeably change. The amplitude of thesound waves in the lower panel is at least a factor of 2 largerthan that in the upper panel and amounts to a large percentageof the azimuthal velocity. The picture shows a distinct pattern ofwave front crossings and reflections for r < r crit and the criticallayer at r crit = .
19 has already begun to smear out at this pointin time. When the color scale is saturated we recognize threeregions where the amplitude of the waves is dropping owing tothe wave fronts crossing each other. They lie at r ≈ . , . .
07. Inside the corotation resonance, the waves are reflectednear the inner boundary and at the critical layer, and are thustrapped between the two radii. We have traced three wave frontswith the black lines to illustrate the pattern that is exactly theoverreflection mechanism that we have already elucidated. Thepattern number of the sonic instability (15 for the simulationshown in Fig. 2) does not depend on details of the inner bound-ary condition. This has been verified by running the same simu-lation with identical perturbations with modified boundary con-ditions, where all variables are kept at their initial values for alltime (“do-nothing” BCs). The pattern number is, instead, set bythe initial random perturbation. For a di ff erent seed of the ini-tial perturbations the pattern number will be di ff erent. When thetrapped waves undergo a reflection at the critical layer, part of thewave action leaks through the critical layer and has the oppositesign compared to the reflected wave, whose amplitude grows asa consequence of the global conservation of wave action. Thus,the whole process is self-sustaining and even grows with time.The leakage of wave action is also clearly visible in Figure 2b)for r > r crit . Now the reason for the change in the growth rate ofthe sonic instability after about ten orbits (see Fig. 1) becomesapparent: It is simply the overreflection mechanism setting in,gaining in strength, and growing dominant. In Section 4.1 we described the linear growth regime of the sonicinstability and analyzed its manifestation in our models. Afteronly 20 −
30 orbits the sonic instability saturates and the sys-tem remains in the excited state with a radial kinetic energy thatis approximately three orders of magnitude larger than beforethe sonic instability set in. There are still considerable variationsin the radial kinetic energy (see Fig. 6); however, we suggestthat they a result of a secondary KH like instability that sets inperiodically and dies out shortly after. This theory is examinedthoroughly in Section 5.2.After the sonic instability has saturated we observe wavesthat are acoustic in nature propagate through the largest part ofthe domain. Those modes have been excited by the sonic insta-bility and do not die out for the whole simulation time, which
Article number, page 6 of 18. Hertfelder and W. Kley: Wave mediated AM transport in astrophysical BLs
Fig. 2. r - ϕ plot of the radial velocity in units of u K ( R ∗ ) at t = t =
10 (b) orbits. The velocity in the upper panel, a), has been multipliedby a factor of 2 in order to use the same colormap as below. Both snapshots show the linear stage of the sonic instability and the middle branch inaction. The sound waves are propagating toward the inner and the outer edge of the domain in an oblique fashion. The three black lines in panelb) illustrate the wavefronts. The arrows denote the direction of the unperturbed azimuthal flow. amounts to over 2000 orbits. Considerable e ff ort has been ex-pended by Belyaev & Rafikov (2012) and Belyaev et al. (2012,2013a) in order to analyze these acoustic modes by means oflinear stability analysis and derive proper dispersion relations.Since this undertaking is quite involved we limit ourselves to re-viewing those elements that are crucial for the understanding ofour results.One finds that the acoustic modes of the BL are related to thewaves in a plane parallel vortex sheet in the supersonic regime(Belyaev et al. 2013a). For a discontinuous drop of velocity inthe isothermal vortex sheet in absence of stratification or rota-tional e ff ects it can be shown that the linearized Euler equationsyield three di ff erent expressions for the dispersion relation ω ,which are denoted the lower, middle , and upper branch (Belyaev& Rafikov 2012). The three branches di ff er in the direction ofpropagation of the modes on either side of the interface (ve-locity drop) and their general appearance is depicted in Fig. 4of Belyaev et al. (2013a). The lower and upper branch are verysimilar to each other. They both show a region where the wavespropagate parallel to the interface and a region where the wavevector forms an angle between 0 and π/ ϕ -direction in the BL and both in r and in ϕ in the disk, where k r / k ϕ (cid:29)
1, however. For the upper branch, the opposite ap-plies. Modes of the middle branch propagate in an oblique man- ner both in the disk and in the BL with the same wave vector.Furthermore, the waves undergo a phase shift of π at the inter-face. This holds true for all branches. The dispersion relations ofthe three branches of the vortex sheet can then be extended byincluding radial stratification and the e ff ects of rotation (Corio-lis force) and one can find approximate dispersion relations fora more realistic case of a BL. This was done in Belyaev et al.(2013a) for the isothermal BL. It would be interesting to see howa realistic equation of state and the radiation transport a ff ect thedispersion relations. However, this is beyond the scope of thispaper.Figure 3 is a r - ϕ plot of the radial velocity after 45 orbits.Three white dashed contour lines of the azimuthal velocity u ϕ have been overdrawn. All velocities are, as usual, given in unitsof the Keplerian velocity at the surface of the star. The imageshows a standing wave in azimuthal direction for radii (cid:46) . r - and ϕ -directions. Additionally, the amplitudeof the wave changes sign at r ≈ .
17, i.e. the phase is shifted by π/
2. This pattern is a clear signature of the lower branch of theacoustic modes, which we can observe in action here. The con-tour lines of u ϕ reveal that the interface between the disk and theBL, or in other words the region where the velocity starts to de-viate substantially from Keplerian, is strongly deformed. The de-formation of the interface is evidently due to the large scale vor-tices that reside at the base of the BL (see Fig. 3). At r = . Article number, page 7 of 18 & A proofs: manuscript no. aa26005-15
Fig. 3.
Two-dimensional plot of the radial velocity u r at t =
45 orbits.The direction of the unperturbed flow is from bottom to top, i.e. in thedirection of growing ϕ . The dashed white lines are contour lines of theazimuthal velocity at the levels 0 . , .
3, and 0 . u K ( R ∗ ). The radial velocity developsa pattern with vortices in the BL and shocks at the top of the BL. Thispattern most likely resembles the lower branch of the acoustic modesdiscussed in Belyaev et al. (2013a). sound speed amounts to ∼
8% of the Keplerian velocity, whereasthe gas rotates with 70% u K ( R ∗ ). Therefore, the azimuthal mo-tion is still nearly ten times super sonic and information aboutthe interface deformation cannot propagate upstream. Hence theoblique shocks, which are clearly visible in the radial velocity,have developed. The whole system described above rotates ina prograde way with a pattern speed of Ω p = ω/ m , where ω is the oscillation frequency in the inertial frame and m the az-imuthal wavenumber. For the episode described here (orbit 45), Ω p ≈ . kr / m (cid:29)
1, where k and m are the radial and az-imuthal wave numbers, respectively. Therefore, one can applythe WKB approximation and assume perturbations of the form δ f ∝ exp (cid:2) i ( kr + m ϕ − ω t ) (cid:3) . (20)A linear stability analysis for a fluid disk with a polytropic equa-tion of state then yields the dispersion relation( ω − m Ω ) = κ + k c , (21)where κ = (cid:16) Ω r dd r ( r Ω ) (cid:17) . is the epicyclic frequency and ω/ m =Ω p (Binney & Tremaine 2008). Of course, Eq. (21) is only anapproximation to the dispersion relation of the waves in ourfully radiative disk. However, it will mainly be used to explainthe wave dynamics in the disk. Furthermore, the WKB approx-imation is no longer valid in the BL where we have seen that kr ≈
0. In the disk we can assume a Keplerian rotation pro-file, Ω ∝ r − . , and consequently the epicyclic frequency isreal. Therefore, there is a region around the corotation resonance Ω = Ω p , flanked by two Lindblad resonances
Ω = Ω p ± κ/ m ,where k is imaginary and the waves are evanescent. Formally,according to Eq. (21) the wave number k becomes zero at the Lindblad resonances. However, the WKB approximation breaksdown at these points since the assumption kr / m (cid:29) < Ω p <
1, there are twocorotation resonances, one in the BL and one in the disk, ei-ther of them flanked by two Lindblad resonances. Assumingthe disk rotates with Keplerian frequency, we can easily derivethe corotation and Lindblad resonance, which are then locatedat r cor = Ω − / and r Lind = ( Ω p m / ( m ± − / . The waves orthe weak shocks we discussed earlier launched at the BL nowpropagate into the disk until they reach the Lindblad resonances,where they are reflected and propagate back toward the BL. Inthe BL the waves might be reflected by Lindblad and corotationresonances as well; however, Eq. (21) is not applicable there.Through this process of consecutive reflections the waves canbecome trapped between the BL and the disk.This idea of the trapped modes has been proposed byBelyaev et al. (2012) and was indeed observed in their simu-lations. In Fig. 3 we do not find such clear evidence for reflectedshocks. Reflection of the wave depends, however, on whetherthe modes are absorbed before they reach the evanescent regionin the disk. We observe that the modes in our simulation canonly travel until r ≈ . α -viscosityapplied in the disk (see Sect. 2.2). We have run simulations withlarger domains in order to investigate a possible influence of theouter boundary on the reflection of the waves. These test runsconfirm the early damping of the waves observed in the simula-tions presented in this work. The e ff ective viscosity in the disk islikely to play a major role in this damping process. The actual in-teraction between the waves and the MRI turbulence in the diskcan, however, only be investigated in real MHD simulations.In Fig. 4 we plot the velocities u r and u ϕ as a function of ( r , ϕ )at t =
173 orbits, a more energetic state than the one depicted inFig. 3 (see Fig. 6 below). It is immediately apparent that the stateat t =
173 orbits is far more violent than the situation shown inFig. 3. The amplitude of the radial velocity has increased by afactor of 2 and the shocks have grown more intense and vigor-ous. While the general pattern of the radial velocity remains un-changed and still constitutes the lower branch with the vorticesat the base of the BL and the outgoing waves at the top of the BL,the wavefronts are now heavily deformed by the shocks. Theseshocks are a consequence of the strongly distorted interface ascan be seen in Fig. 4b). The lines of constant u ϕ in the r - ϕ space(the contours) have a wave-like shape with many irregularitiesand the transition from zero rotation rate to Keplerian rotationis spread over ∼ . R ∗ and has broadened considerably. Thoughthe radial velocity is more chaotic in Fig. 4b) than in Fig. 3, thereis still a dominant mode visible and the azimuthal wave numbersare m =
14 and 24, respectively. As we shall see later, the an-gular momentum transport by the waves is very e ffi cient in thisstage.For completeness, in Fig. 5 we show a r - ϕ -plot of the radialvelocity that spans the whole simulation domain. The snapshotwas taken after 200 orbits while the system was in a quiet state.The pattern of the lower branch with its vortices at the base of theBL is clearly visible. We can also see that the waves emergingfrom the top of the BL are indeed tightly wound and that thespirals are leading. Figure 5 shows another azimuthal wave with m = m = π simulations. Thewavenumber of the high-frequency pattern at the base of the BLis given by m = Article number, page 8 of 18. Hertfelder and W. Kley: Wave mediated AM transport in astrophysical BLs
Fig. 4. r - ϕ plot of the radial velocity (a) and the azimuthal velocity (b) at t =
173 orbits. Both quantities are given in units of u K ( R ∗ ). The dashedwhite lines in both plots trace the contour of u ϕ = .
15 and 0 . u r and the heavy deformation of the interface (see b) are clearly visible. Fig. 5.
Full-disk plot of u r at t =
200 orbits. The direction of the unper-turbed flow is counterclockwise. The radial velocity is given in units of u K ( R ∗ ). The middle branch is a transient mode and dominates duringthe linear growth of the sonic instability. It can be observed inthe early stages of the simulation and the oblique wavefronts oneither side of the interface are clearly visible in Fig. 2. Both panela) and b) show features of the middle branch, yet they are moredistinct in b). The upper branch, on the other hand, might beobserved during the late times of the simulations. We defer thediscussion of this issue to Sect. 5.5.
5. The effects of the instabilities
In the previous section we discuss the hydrodynamical instabil-ities in the BL in detail. Now we show the implications of the r a d i a l k i n e t i c E n e r g y [ e r g ] Fig. 6.
Temporal evolution of the radial kinetic energy of our referencesimulation. The time is given in units of the Keplerian orbit at r = R ∗ , P ∗ . The switching between quiet and outburst states described in thetext is clearly noticeable. We have marked the exact times for some ofthe outbursts for easier comparison with the text. BL instabilities and the physical results of our simulations of theradiative BL around a young star.
Figure 6 visualizes the total radial kinetic energy, E kin,rad = (cid:82) disk 12 Σ u r d A , for which we summed the radial kinetic energy ofevery cell in the simulation domain. For a given time, the radialkinetic energy reflects the strength of the instabilities. There-fore, we can monitor the time-dependence of the instabilitiesby tracking E kin,rad . The reference simulation shows a stronglytime-dependent behavior that can be characterized as a recurringseries of states of outburst and quiescence. When the kinetic en- Article number, page 9 of 18 & A proofs: manuscript no. aa26005-15 ergy reaches a maximum, the instabilities are vigorous and theinterface between the BL and the disk is heavily deformed (seeFig. 4). In the quiescence states, on the other hand, the amplitudeof the velocity perturbations drops considerably, often by a fac-tor of more than 10, and the interface deformation is less stronglypronounced. The strength of the outbursts does not follow a clearpattern. In the beginning, the outbursts get stronger, with the sec-ond outburst at t ≈
173 orbits being the most powerful. Duringthis outburst the amplitude of the perturbations in the radial ve-locity reaches over 20% of the Keplerian velocity at the surfaceof the star. Thus, the amplitude of the radial velocity becomescomparable to that in azimuthal direction, or even exceeds it for r (cid:46) .
1. The third outburst at t ≈ M increases by more than twoorders of magnitude. In the high activity phases, ˙ M distinctlyrises in the outer parts of the simulation domain as well. The an-gular momentum transport, j , shows a similar trend. At times ofstrong instability, the transport of angular momentum, which ismainly directed to the center of the disk, also grows remarkably,quite often by a factor of about 100.Unlike in a stationary state where mass and AM are trans-ported to the center of the disk (see Eq. 17), we observe thatthe transport of mass does not have a preferred direction in cer-tain regions of the simulation domain. We consider for this pur-pose the upper panel of Fig. 7. There is a clear boundary at r ≈ . − . R ∗ , which is approximately given by the locationof the maximum of Ω and divides the rather organized, inwarddirected mass flux in the disk from the alternating transport inthe BL, where mass is sloshing about back and forth in radiusbecause of the passage of waves which have alternating signsof mass transport. This is clearly visible in the form of the di-agonal features in the upper panel of Fig. 7. Thus, the direc-tion of the mass transport in the BL is both a function of timeand radius, and the variability with time is on the order of oneorbit. After roughly 600 orbits, the strength of the mass trans-port ceases slightly, only gaining strength again during the sub-sequent outburst periods. This development matches that of theradial kinetic energy (Fig. 6) where we can also detect a changein activity starting after 600 orbits. After 1200 orbits another ma-jor change in the mass transport pattern takes place. The wholesystem seems to calm down and the sloshing behavior with itsfrequent changes in direction ceases.The AM flux does not feature a variability as intense as themass flux; however, the diagonal features that have also beendiscussed in the context of ˙ M are still present. An analysis ofthe individual components of j confirms that these features arisefrom the advective part of the AM flux, which shows a sloshingbehavior similar to the mass flux. Apart from these features, j assumes mostly negative values (see the lower panel in Fig. 7),meaning that AM is transported inward. This evidence agreeswith our expectations that in the disk j is negative for the rea-sons pointed out earlier. One very distinct feature in the graph r [ R ∗ ] Ω [ Ω K ( R ∗ ) ] Orbit 0Orbit 10Orbit 50Orbit 200Orbit 600Orbit 800Orbit 1000
Fig. 8.
The velocity in azimuthal direction in units of the Keplerian ve-locity at the stellar surface for seven di ff erent orbits. The curves demon-strate the spin-up of the gas for r (cid:46) . r (cid:38) . of the AM flux is the very pronounced boundary that starts atthe beginning of the simulation at r close to 1 . R ∗ and movesinward up to r ≈ . R ∗ during the first 800 orbits. This linecan be characterized by an enhanced AM flux in its near vicin-ity and a shift in direction from inward ( r < r line ) to outward( r > r line ). Figure 8 visualizes the destination of the AM whichis accumulated at this interface as a result of transport from bothsides. It is spent to spin up the gas at r (cid:46) .
2. Hence, the BL isbroadened considerably during the course of the simulation. TheAM needed for the spin-up of the inner layers is extracted fromthe faster rotating part at r (cid:38) .
2, which is, consequently, sloweddown. By comparing Fig. 7 and Fig. 8, we can deduce that theAM transport is the strongest where the gradient of Ω is highest,i.e. where the shearing is strongest.It can be recognized from Fig. 7 that AM is mostly trans-ported inward although the mass flux is also directed outward.Especially when ˙ M > M over the wholeazimuthal domain for a given radius can therefore be positive(meaning that mass is transported outward). However, accord-ing to Eq. (17), the total advective AM flux can easily remainnegative if u ϕ is smaller at azimuthal locations of positive ˙ M than where it is negative. This is indeed very often the case andwe have picked an illustrative example to stress this point (seeFig. 9). It is a snapshot of the mass accretion rate, the radial ve-locity, and the angular velocity at orbit 145, where we plottedthe dependence on azimuth for r = . R ∗ . The figure clearlyshows that whenever ˙ M ( ϕ ) is positive and large, u ϕ ( ϕ ) is small.The comparison between the radial and the azimuthal velocityreveals that the anti-correlation between u r and u ϕ is such that Article number, page 10 of 18. Hertfelder and W. Kley: Wave mediated AM transport in astrophysical BLs
Fig. 7.
The upper panel shows the temporal evolution of the mass accretion rate in the radial direction, ˙ M ( r , t ), in units of solar masses per year. Itis obtained from the mass accretion rate per cell (see Eq. 16) by integrating over the azimuthal direction. The lower panel likewise visualizes theAM flux normalized to the advective AM flux at the stellar surface (see Eq. 18), which is also integrated over azimuth. We note that we applieda so-called symlog scaling for the plot, which means that we have a logarithmic scaling for both the positive and the negative values. Only in thevicinity of 0 is the scale linear. A negative value means that mass or AM is transported to the center of the disk. The phases of strong activity inthe radial velocity (see Fig. 6, marked with the black dashed lines here) correspond to phases of increased mass and AM transport. The diagonalfeatures in the upper panel are a clear signature of waves, which have alternating signs of mass transport. This sloshing behavior spreads to theAM flux (lower panel) via the advective component of j . an extreme of the azimuthal velocity coincides with a root of theradial one. Therefore, the two patterns that are similar, althoughnot perfectly equal, are shifted by one quarter of the azimuthalwavelength of the m = We now take a closer look at the cause of the periodic outburstsand the behavior of the primitive variables in the disk duringsuch an outburst and we consider Fig. 10 for this purpose. Thegraphic shows the velocity in azimuthal direction, the surfacedensity, the pressure, and the temperature in the midplane of thedisk for five di ff erent times that are centered around the outburstat t =
173 orbits and are not equidistant. These time steps coverthe second outburst in Fig. 6. This part of the radial kinetic en-ergy evolution is shown as a zoom-in in the lower right panel ofFig. 10.At the beginning of the outburst, the surface density has in-creased quite strongly, especially in the vicinity of r ≈ . ffi cient enoughto guide it all through. This situation is universal for all ob-served outbursts, which are preceded by a phase of quiescenceduring which the activity and the mass transport in the BL arenot as strong as in outburst phases. Therefore, the whole pro-cess described here using the example of the first outburst can beequally applied to the other outbursts. During the accumulationof the mass just outside the BL, the azimuthal velocity profile issmooth and quasi-stationary. The high amount of mass leads toa high temperature and hence a high pressure, which means that u ϕ is remarkably sub-Keplerian in the first half of the domainsince the gas is stabilized by pressure.After a su ffi cient amount of mass has accumulated, a run-away process sets in during which u ϕ and the other quantitieschange very rapidly. According to Figs. 7 and 10 a great amountof mass is guided through the BL onto the star over such an out-burst. The process is initiated outside the BL at about r ≈ . ff ect of the deceleration is much stronger,mass moves inward. As a consequence of the decreasing density, Article number, page 11 of 18 & A proofs: manuscript no. aa26005-15 u ϕ [ u K ( R ∗ ) ] Orbit 160Orbit 170Orbit 173Orbit 180Orbit 190 Σ [ g c m − ] r [ R ∗ ] P [ d y n ec m ] r [ R ∗ ] T c [ K ]
120 140 160 180 20001234567 radial kin. Energy
Fig. 10.
Plots of the azimuthal velocity, the surface density, the pressure, and the midplane temperature for five time steps centered around thesecond outburst. The dashed lines denote the initial conditions of the simulations. The dash-dotted line in the upper left panel gives the azimuthalvelocity for pure Keplerian rotation. In the lower right panel we have included a zoom-in of Fig. 6 for better identification of the exact location ofthe plot times. the pressure decreases and the gas is even less stabilized and canaccelerate its inward movement. After the density is satisfacto-rily depleted, the Reynolds stresses decrease and the violent ac-cretion process comes to an end. Then the density grows againbecause of the lower stresses and the azimuthal velocity is slowlyregains its pre-outburst profile.The variation of the midplane temperature during the out-burst phase is consistent with what we have already discussedfor the other physical quantities. At the beginning, despite thedecreasing density, the temperature rises for r (cid:46) . R ∗ and evendevelops a small bump around r = . R ∗ . The increase in tem-perature is due to the dissipation of energy through the growingReynolds stresses and indicates the region where the energy re-lease by the waves is largest. This radius is in agreement withwhat we expected from the plot of the azimuthal velocity if wealso take into account that the energy is redistributed by radia-tive di ff usion. Subsequently, the surface density drops more andmore and the Reynolds stresses decrease as well. Accordingly,the temperature also drops, and even drops below the level ofthe initial state. After enough mass has again piled up around theBL, the temperature will rise to its pre-outburst level.The outburst sequence shown here for the second outburst ofthe reference simulation is exemplary for the general behaviorof outbursts. What we have just discussed might be less clearlyvisible or even more distinct, depending on the strength of theoutburst. However, the mechanism is the same throughout. From simulation C, which has a resolution of 2048 × r < . R ∗ andthe depletion of mass is more e ff ective. Thus, it is possible thatthe instabilities become stronger with increasing resolution. Theincreasing strength of the outbursts might instead be due to theother modes becoming excited through the initial perturbations,as is typical for classical KH instability.We propose the following scenario for the physical cause ofan outburst such as the one described above. The radial veloc-ity has developed weak shocks at the top of the BL. Every timea fluid element passes through that shock it is slowed down bythese oblique shocks. The cumulative e ff ect of all these shockpassages leads to a decrease in the azimuthal velocity of the gasand creates a wide plateau in the velocity profile (see e.g. orbit160 in Fig. 10). Such a flat region with a weak velocity gradi-ent is prone to the classical Kelvin-Helmholtz (KH) instability,which sets in and induces the rapid changes in the BL. In ab-sence of rotation, it follows from Rayleigh’s stability equation that a necessary condition for instability is the occurrence of aninflection point in the velocity. A stronger form of this condi-tion was given by Fjørtoft (1950), who showed that the velocityprofile not only needs to have an inflection point, but also has tosatisfy the condition u (cid:48)(cid:48) ( u − u s ) < u s is the velocity at the inflection point. However, this require- Article number, page 12 of 18. Hertfelder and W. Kley: Wave mediated AM transport in astrophysical BLs ˙ M [ − M s o l y r − ] u r [ u K ( R ∗ ) ] π/ π/ π/ π/ π/ π π/ π/ π/ π/ π/ π Azimuth ϕ u ϕ [ u K ( R ∗ ) ] Fig. 9.
The upper panel depicts the mass accretion rate, the middle panelthe radial velocity, and the lower panel the azimuthal velocity. All threequantities are plotted as functions of azimuth at the fixed radius r = . R ∗ after 145 orbits. The upper two curves are shifted in phase withrespect to the lower one such that maxima in ˙ M (nearly) coincide withminima in u ϕ . ment is only true for parallel flows. A normal-mode analysis ofthe linearized Euler-equations of a rotating incompressible flowwith zero mean radial velocity eventually yields( ω + im Ω ) (cid:32) ∂ ∂ r + r ∂∂ r − m r (cid:33) φ − im r ∂ Z ∂ r φ = φ ( r ) is the amplitude of a streamfunction ψ (cid:48) ( r , ϕ, t ) = φ ( r ) exp( ω t + im ϕ ) such that u (cid:48) r = im φ/ r and u (cid:48) ϕ = − ∂φ/∂ r , where the prime denotes perturbed quantities; Z = / r ∂ ( ru ϕ ) /∂ r is the vorticity. Equation (22) is the equiva-lent of Rayleigh’s equation for a rotating fluid. Rayleigh showedthat if φ = R and R enclosing the flow, anecessary condition for instability is that the gradient of Z mustchange sign at least once in the interval R < r < R (Drazin &Reid 2004). We take this analogue of Rayleigh’s inflection pointtheorem as a motivation to study the vorticity and its derivative,though Eq. (22) is clearly an idealization of the situation at hand.However, a detailed stability analysis of the RHD equations of arotating stratified fluid is beyond the scope of this work.Figure 11 displays the derivative of the vorticity, ∂ Z /∂ r , forthree points in time. At orbit 50, the system is in a state of quies-cence (see also Fig. 6). After 150 orbits the second outburst hasjust started to evolve, and 50 orbits later the system is recover-ing from the outburst. While the vorticity derivative of the quietand post-outburst state does not change sign for r > . R ∗ , itdoes change several times for orbit 150. The region r > . R ∗ isof interest here because this is where the rotation profile is flat.There are, of course, multiple changes of sign for r < . R ∗ inall curves. However, those points are associated with the strong r [ R ∗ ] d e r i v a t i v e o f t h e v o r t i c i t y orbit 50orbit 150orbit 200 Fig. 11.
Derivative of the vorticity, ∂ Z /∂ r , as a function of radius fororbit 50 (quiet state), orbit 150 (pre-outburst state), and orbit 200 (post-outburst state). The inset is a zoom-in of the region 1 . ≤ r ≤ . shearing regions that we have already discussed in Sect. 4. Afteranalyzing several time steps before, during, and after outbursts,we draw the following picture: The consecutive shock passagesnear the maximum of Ω cause the flow field to adjust such thatmultiple changes of sign are created in the derivative of the vor-ticity. Since this is a necessary condition for the instability of arotating, hydrodynamical flow, it is probable that a classical KHinstability develops in the flat regions of the flow. The role, how-ever, that the KH instability plays in the context of the enhancedaccretion and AM transport is not yet entirely clear. It might, onthe one hand, act as a trigger for a mode of the sonic instabilitywhich is excited during the outburst and e ffi ciently drives the ac-cretion. Once a su ffi cient amount of the mass in the BL has beenaccreted, the mode can no longer be sustained and the systemreturns to quiescence. One candidate for this mode might be theupper branch, which we normally do not observe in our simu-lations. A scenario similar to this option was also discussed inBelyaev et al. (2013a). If, on the other hand, the KH instabilityitself is responsible for the enhanced accretion, AM should betransported mainly by turbulent stresses, as opposed to the wavetransport through the mode mentioned above. We did not findclear signs of turbulence during the outburst, nor did we find ev-idence of the upper mode. It is, however, plausible that the AMtransport is still accomplished by waves and the KH instabilitymight trigger or enhance a mode of the lower branch.We did find evidence for the existence of a KH instabilityduring the outburst. We consider for this purpose Fig. 12, a 2Dimage of the z -component of the vorticity ∇× u at orbit 172 .
6, i.e.at the peak of the second outburst. It shows patterns that are com-monly associated with the classical KH instability, for instance,the curly arms and the cat’s eyes. By analyzing simulation C wefind that these patterns become more distinct with increasing res-olution (see also Fig. 11 in Belyaev et al. 2012). The KH featuresare, however, masked to some extent by the waves generated bythe main instability. This becomes clear when we investigate theoutbursts at later times when the main instability is decreasing.In addition, the flat region where we expect KH instability towork occurs in the fast rotating part of the disk. Static grid codessu ff er from disadvantages when the fluid is boosted and fluidmixing instabilities can be suppressed (Springel 2010). Article number, page 13 of 18 & A proofs: manuscript no. aa26005-15
Fig. 12.
Image of the z -component of the vorticity, ∇ × u , at the peak ofthe second outburst at t = . r [ R ∗ ] T e ff [ K ] Orbit 160Orbit 170Orbit 173Orbit 180Orbit 190Orbit 210
Fig. 13.
The azimuthally averaged e ff ective temperature of the BL forsix time steps centered around the second outburst (see Fig. 6). T e ff in-creases considerably during the outburst and then declines to approxi-mately the pre-outburst level ∼
30 orbits after the peak of the outburst.The black dashed line corresponds to the surface temperature of thestandard solution of Shakura & Sunyaev (1973).
We now investigate the outburst described in Sect. 5.2 from anobserver’s point of view and look at the radiation emitted by theBL and the disk. In Fig. 10 we have exclusively considered phys-ical quantities, which are measured in the equatorial plane of thedisk, among them the midplane temperature T c . If the medium isoptically thick, however, the emergent spectrum is dominated bythe temperature on the surface of the disk where the vertical op-tical depth is approximately one, and locally resembles the spec-trum of a blackbody radiating with temperature T e ff . Since wehave no knowledge of the vertical structure of the disk, we usethe approximation of Hubeny (1990), which is a generalizationof the gray atmosphere (e.g. Rybicki & Lightman 1986). The ef-fective temperature can then be calculated with an approximated optical depth in vertical direction (Suleymanov 1992): T ff = T /τ e ff τ e ff = τ R + √ + τ P . (23)Here, τ R and τ P ( τ = κρ H = . κ Σ ) are the Rosseland and thePlanck mean optical depth, respectively.Figure 13 shows the e ff ective temperature of simulation Bcalculated according to Eq. (23). At the beginning of the outburst(orbit 160) the e ff ective temperature features only a slight peakin the BL and the temperature in the BL and in the disk di ff er atmost by 500 Kelvin. The reason why the e ff ective temperature ofthe BL and the disk are very similar at the beginning of the out-burst episode is the following. The system, which is just on theverge of an outburst at this point in time, comes out of a preced-ing state of quiescence. Therefore, the energy dissipation of thewaves is small and the density is rather high owing to the ine ffi -cient mass transport. Consequently, both the midplane tempera-ture and the optical depth in the BL are comparable to the diskand T e ff is only slightly larger in the BL. Over the course of thenext 13 orbits, however, the system goes through an episode ofe ff ective mass and AM transport, and the midplane temperatureincreases and the surface density decreases (see Fig. 10). Thus,according to Eq. (23), T e ff increases rapidly in the BL yieldinga maximum of 6000 K at r ≈ .
3, which is twice the temper-ature in the disk. Care should be taken when comparing T e ff inour models with conventional BL models where an α -viscosityhas been applied. The latter typically feature a BL that is con-siderably smaller in radial extent than in the models presentedhere. The luminosity of the BL, however, does only depend onthe physical parameters of the system and must be the same inboth models ( GM ˙ M / R ∗ for a non-rotating star). Since the lumi-nosity is the integral of the flux F = σ SB T ff over the area of theBL, it follows that the thinner the BL, the higher the e ff ectivetemperature must be, and vice versa. The relation between themaximum e ff ective temperature in the BL θ BL and in the disk, θ d , can roughly be estimated by (Lynden-Bell & Pringle 1974) θ BL θ d = . (cid:32) δ R ∗ (cid:33) − / , (24)where δ is the width of the thermal BL (see e.g. Regev & Bertout1995; Popham & Narayan 1995). Equation (24) yields a factorof ∼ δ = . R ∗ , which is in good agreement with whatwe observe in our simulation (see Fig. 13). Also visualized inFig. 13 (black dashed line) is the surface temperature accordingto Shakura & Sunyaev (1973), which is given by T e ff ( r ) = (cid:32) GM ∗ ˙ M π r σ SB (cid:34) − j (cid:18) R ∗ r (cid:19) / (cid:35)(cid:33) / , (25)where j is the normalized AM flux. We took j = . ff ective temperature derived from oursimulation agrees very well with Eq. (25) in the disk. In the BL,of course, the analytically derived formula fails to capture therise of the e ff ective temperature.Another important point for the observational appearance ofthe BL is the azimuthal dependence of the e ff ective tempera-ture, which is visualized in Fig. 14 for orbit number 173. We cansee from this depiction that T e ff reaches temperatures of nearly10 000 K in the simulation domain which is high, consideringthe case of a young star. There is strong variation of T e ff in ϕ -direction, especially for the annulus r ≈ .
3, which is due to
Article number, page 14 of 18. Hertfelder and W. Kley: Wave mediated AM transport in astrophysical BLs
Fig. 14.
2D plot of the e ff ective temperature of the BL at orbit 173,which marks the peak of the second outburst (see Fig. 6). T e ff stronglyvaries in azimuthal direction and reaches temperatures up to 10 000 K.The curly arms visible for r < . the shocks of the radial velocity that create density perturbationsin the surface density and then a ff ect the vertical optical depthand the e ff ective temperature. The ϕ -dependent pattern of T e ff rotates with the pattern speed of the acoustic modes, which isclose to the Keplerian velocity, and will express itself in time-dependent variations of the light curve. This is an importantpoint because many systems show light curve oscillations withperiods that are close to the orbital period at the surface of thestar. In cataclysmic variables, for instance, one observes rapidoscillations called dwarf novae oscillations (DNOs, Warner &Robinson 1972) and quasi-periodic oscillations (QPOs, Patter-son et al. 1977). Both phenomena are not yet fully understoodand it is tempting to see a connection to the BL modes on thebasis of the protostar simulations presented here. However, fur-ther research is necessary to clarify this point. We also note thatFig. 14 shows clear indications of KH instabilities. These fea-tures are even more pronounced in snapshots of T e ff taken beforethe peak of the second outburst and reinforce the point that aKH-like instability is responsible for the frequent outbursts. In Fig. 7 we visualized the total AM flux in the disk, which con-sists of the advective part, the viscous part in the disk, and theReynolds part (see Eq. (17)). It is, however, of great interest tosee how much the instabilities contribute to the total AM fluxand to deduce the e ffi ciency of this transport. As we have de-scribed in Sect. 4, the instability in the BL leads to the devel-opment of acoustic waves that are launched from the interfaceand propagate into the disk and the star. As a consequence, AMtransport in our models is dominated by acoustic radiation, i.e.sound waves carrying AM. We do not observe any turbulence inour simulations and hence rule out turbulent stresses as a ma-jor contribution to the AM flux. Unlike the turbulent stresses,the AM transport by waves is a non-local process. Waves cancarry AM of either sign from one point in the domain, e.g. theinterface in the BL where they are launched to another point inthe disk where they deposit the AM into the fluid through wavedissipation. Thus, the long-range exchange of angular momen- B L t h i c k n e ss i n R ∗ simulation Asimulation Bsimulation Csimulation D Fig. 16.
The width of the BL as a function of time for the four simula-tions mentioned in Sect. 3.3. The width is calculated from the radii ofthe region where 0 . < u ϕ / u K < .
8. All four simulation settle downto a BL width of approximately 20 percent of the stellar radius. Simu-lation A approaches this value more slowly and in a slightly di ff erentmanner, which is due to the very low resolution of this particular simu-lation (512 × tum between regions that are otherwise uncorrelated is possiblethrough the wave transport. This is a major di ff erence from the α -models which yield a local e ff ective viscosity due to a localAM exchange phenomenon, such as turbulence. Therefore, it isquestionable whether the processes in the BL can be describedby a local prescription like the α -model.We now analyze the e ffi ciency of the wave mediated AMtransport and consider Fig. 15 for this purpose. It shows the di-mensionless Reynolds stresses α Re , which has been calculatedaccording to Eq. (14), as a function of radius and time. We canapply this description to the AM transport since waves trans-port angular momentum exclusively by stresses. The Reynoldsstress is evidently very strong with peak values up to the orderof unity and it is largest at the top of the BL, at around r = . Ω reaches its max-imum. Thus, our simulations reveal that the stresses are largestaround the very point where they would vanish according to aconventional α -viscosity model. Strictly speaking, this renders alocal prescription that depends on the shear of the angular veloc-ity inapplicable for the BL. Figure 15 also shows that the AMtransport due to the waves includes a large part of the domain,1 . (cid:46) r (cid:46) .
4, and spreads out farther into the disk during out-bursts. This reinforces the point that the AM transport ultimatelygenerated by the shear in the BL is not confined to the BL, butin fact ranges over long distances. As we saw earlier, the activityrises during outbursts and therefore α Re is enhanced during theseperiods as well. The value of α Re is almost exclusively negative,meaning that angular momentum is transported toward the star.This behavior is characteristic of the lower branch of the acousticmodes (Belyaev et al. 2013a). We refrain from presenting plotsof the parameter α ν that can be calculated according to Eq. (15)because it is artificial to boil the non-local AM transport downto a local parameter. However, for the sake of completeness, itshould be mentioned that Eq. (15) yields values on the order of10 − − ff ects the width of the region where the az- Article number, page 15 of 18 & A proofs: manuscript no. aa26005-15
Fig. 15.
Time-radius image of α Re spanning several hundred orbits. We note that the values of α Re have been multiplied by ( −
1) in order to allowfor the logarithmic scale. For each point in time the values have been time average over ten orbits. The vertical dashed lines denote the locationsof the outbursts mentioned earlier in the text. imuthal velocity rises from stellar rotation to the Keplerian ro-tation in the disk. Figure 16 depicts the BL thickness as a func-tion of time for all four simulations. The thickness of the BL isdefined to be the region in which0 . < u ϕ ( r ) / u K ( r ) < . . (26)The width of the BL increases very rapidly and grows in sizeby a factor of 2 at the beginning of the simulation, i.e. when thesonic instability has grown to establish an e ffi cient AM trans-port. After the initial fast growth the width of the BL increaseson average linearly over the course of approximately 500 orbitsuntil it settles down at a value of ∼ . R ∗ . We have run additionalsimulations in order to test the robustness of the BL thickness of20 % of the stellar radius. In these simulations, the initial dis-tance of the velocity drop from the inner edge of the simulationdomain has been extended by altering both the parameter F (seeSect. 3.2) and the point where the computational domain begins.We thereby test whether the thickness is a ff ected by the BL run-ning up against the inner edge of the domain. However, theseruns are in perfect agreement with the simulations presented hereand saturate at a BL width of ∼ . R ∗ as well.There are several spikes in the temporal evolution of the BLwidth for all three simulations that are clearly associated with thefrequent outbursts. During the high activity state the BL broad-ens considerably in a very short period of time and typically re-sides at a larger width in the aftermath of the outburst. Becauseof the high resolution and the resulting long simulation timeswe can follow the highly resolved simulation only for about 200orbits. During that time it behaves exactly like the other simula-tions with the exception that the width of the BL increases muchmore quickly. This is in good agreement with the fact that thegrowth rate of the sonic instability depends on the resolution insuch a way that it is faster when the the resolution is higher (seeSect. 4). Since the BL widens to a large extent, it extends very farinto the disk and into the star. The latter means that a consistentsimultaneous treatment of the star is required for numerical sim-ulations involving large simulation times. We will discuss thispoint in detail in the next section. We ran and monitored our reference simulation for over 2000orbits (see Fig. 6), and the instability and the periodic outbursts do not die out during that time. There is, however, a change inthe pattern of the acoustic waves that occurs around orbit 1200,which coincides precisely with the moment when the base ofthe BL reaches the inner boundary of the simulation domain.We consider this a problematic moment, since it is not possibleto treat the star correctly in 2D r - ϕ simulations. The geometryof the coordinate system inherently gives the star a cylindricalshape in our simulations. In Sect. 3.2 we describe how we treatthe inner boundary in order to approximate the beginning of thestar and that the radial velocity is set to outflow with a certain ve-locity at this point. Therefore, we might expect unphysical andartificial processes once the vortices of the acoustic modes reachthis point because the radial velocity is fixed to a certain valueand an inflow in the domain is not supported. Furthermore, weforce the azimuthal velocity to zero at the inner boundary. How-ever, the acoustic modes seem to be quite robust since they do notdie out despite the di ffi culties with the artificial inner boundary.The change in the pattern at orbit 1200 might be due to a switch-ing of the modes to the upper branch, since the spiral arms of thewaves are unwinding considerably and changing from leading totrailing. Another 1000 orbits later the process is reversed againand the system switches back to the lower branch.However, in order to capture the physics correctly at theselate stages in time and to still make reliable statements once theBL extends all the way to the surface of the star or even intothe star, it is necessary to attach a realistic model of the star tothe disk. We consider, for instance, a wave launched in the BLthat propagates toward the star. It can penetrate deep into the starbefore it is stopped by the steep density gradient and dissipatesthe angular momentum somewhere in the star. On the other hand,many problems arise when the star is treated consistently. Theproblem must then be simulated in spherical coordinates, i.e. 3Dsimulations are necessary, and the density varies considerablyfrom the interior of the star to the disk. Simulations of this kindare very demanding from a computational point of view, but alsovery tempting, and we plan to investigate this problem in thefuture.
6. Summary and conclusion
We have performed highly resolved, 2D hydrodynamical simu-lations of the BL surrounding a young star in order to investigatethe AM transport driven by instabilities. Extending previous sim-
Article number, page 16 of 18. Hertfelder and W. Kley: Wave mediated AM transport in astrophysical BLs ulations, we have a net mass flow through the disk, we utilizeda realistic equation of state, and (within the framework of one-temperature approximation) we included full radiative transferin the disk plane as well as vertical cooling through an realis-tic estimate of the vertical optical depth, thereby employing aquasi 3D radiation transport. Moreover, those simulations werestarted from state-of-the-art 1D models of the BL and thus com-prise realistic profiles for the density and temperature. We ranthe simulations for over two-thousand orbits in order to studythe long-term behavior.Our simulations confirm that the supersonic velocity drop inthe BL is indeed prone to the sonic instability, a kind of super-sonic shear layer instability whose subsonic counterpart is theKelvin-Helmholty instability (Belyaev & Rafikov 2012). Shortlyafter starting the simulations, the sonic instability sets in andstarts to grow rapidly for about 15 orbits. It then saturates andthe BL is dominated by acoustic waves that propagate both intothe disk and into the star. These acoustic modes manifest them-selves in three distinct patterns that can be related to the threebranches of the instability of a plane parallel vortex sheet and donot die out for the whole simulation time.Additionally, the system repeatedly undergoes outburstswhere the wave activity as well as the AM and mass transportincrease considerably. We argue that these outburst are likelytriggered by a secondary KH instability that develops in the flatregion of the azimuthal velocity profile due to several changes ofsign of the vorticity. Since the e ff ective temperature also showsstrong variations during outbursts, they may play an importantrole in explaining variations in the light curve, such as FU Or-outbursts or DNOs and QPOs in the case of BLs around whitedwarfs. It is tempting to associate these phenomena with the out-bursts we observe in our simulations; however, this topic must belooked into further in order to draw reliable conclusions. Such aninvestigation would certainly involve the study of how the out-burst reacts to a change in parameters, among other things. It isalso vital to perform 3D simulations on this question since theremight be a dependence on dimensionality (Belyaev et al. 2013a).Our main objective in this work was the clarification of theAM transport in unmagnetized astrophysical BLs. We found thatthe acoustic modes indeed transport AM through the BL and thedisk, and that this mechanism is even highly e ffi cient, reachingvalues for the dimensionless Reynolds stress α Re of ∼ .
01 inquiet states and up to unity in high activity states (see Fig. 15).Along with AM, mass is also transported e ffi ciently through thedomain (see Fig. 7). The BL reacts to the enhanced transport bywidening considerably until it reaches a thickness of about 0 . R ∗ (see Fig. 16).One of the most important elements of wave mediated AMtransport is the intrinsic non-locality of this process, i.e. wavescan extract AM from one point in the disk and release it in an-other farther o ff . There are several implications that emerge fromthis: Thus far, α -models for the parametrization of the viscosityhave been used both in the disk and in the BL, possibly withsome corrections for the BL. This prescription is justified in thedisk where, under certain conditions, the system is susceptibleto the MRI. The MRI leads to turbulence in the gas and angularmomentum is transported via turbulent stresses. This is a localprocess that also depends on the shearing at this point and a lo-cal e ff ective viscosity can be derived by the model of Shakura& Sunyaev (1973). Furthermore, it allows for e ffi cient mixing ofthe gas through the turbulence. In contrast, the AM transport wehave investigated in this work does not depend on the local shear-ing, for instance, but is a non-local mechanism, instead. There-fore, it is doubtful whether a conventional α -model is applicable in the BL, and it will be hard to give a simple parametrization forthe AM transport in the BL at all. However, α Re features a rathersmooth behavior and might be utilized for a simplified picture ofthe stresses in the BL. In addition, the mixing in the BL mightnot be as pronounced as in the disk, which could have importantimplications concerning the spectrum of the BL.Connected with the notion of non-local AM transport is thefact that the kinetic energy, which resides in the gas shortly be-fore it is decelerated to stellar rotation, does not need to be re-leased as locally as was previously assumed. Waves also trans-port energy and release it to the fluid where they are damped ordissipate. This might have important implications for the struc-ture and observational appearance of the BL. We consider, forinstance, the alternative BL theory of the spreading layer (Inog-amov & Sunyaev 1999, 2010) where the gas is not deceleratedin the disk midplane but rather in two belts in the northern andsouthern hemisphere of the star. It is di ffi cult to bring a local vis-cosity model into accordance with this theory, since the frictionon the surface of the neutron star is far too small. The non-localAM transport described in this work could, however, shed newlight upon the theory of the spreading layer.Simulations that investigate both the vertical structure andthe non-axisymmetric properties can only be performed in 3D,however. This is a vital extension that should be taken as the nextstep in order to bring the models closer to reality. In a 3D spheri-cal geometry, it would also be possible to simultaneously modelthe star, a point that has been illustrated in Sect. 5.5. Another ma-jor extension is the inclusion of magnetic fields, i.e. performingmagnetohydrodynamical simulations of the BL. It would then bepossible to relinquish any viscosity prescription whatsoever anddirectly simulate the MRI in the disk and the acoustic modes inthe BL. There might be important interactions between the tur-bulent disc and the BL modes and new wave branches might addto the three discussed here. However, the extensions mentionedabove will greatly increase the computation time and possiblyrequire more e ffi cient codes and faster hardware. Acknowledgements.
Marius Hertfelder received financial support from the Ger-man National Academic Foundation (Studienstiftung des deutschen Volkes).This research made use of matplotlib (Hunter 2007), a python module fordata visualization. Most of the simulations were performed on the bwGRiD clus-ter in Tübingen, which is funded by the Ministry for Education and Research ofGermany and the Ministry for Science, Research and Arts of the state Baden-Württemberg, and the cluster of the Collaborative Research Group FOR 759,funded by the DFG. We also thank the referee for his constructive commentswhich helped to improve this paper.
References
Abramowicz, M., Brandenburg, A., & Lasota, J.-P. 1996, MNRAS, 281, L21Balbus, S. A. 2003, ARA&A, 41, 555Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214Balbus, S. A. & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1Balbus, S. A. & Lesa ff re, P. 2008, New A Rev., 51, 814Balbus, S. A. & Papaloizou, J. C. B. 1999, ApJ, 521, 650Baruteau, C. 2008, PhD thesis, CEA Saclay, Service d’Astrophysique, 91191Gif / Yvette Cedex, FranceBell, K. R. & Lin, D. N. C. 1994, ApJ, 427, 987Belyaev, M. A. & Rafikov, R. R. 2012, ApJ, 752, 115Belyaev, M. A., Rafikov, R. R., & Stone, J. M. 2012, ApJ, 760, 22Belyaev, M. A., Rafikov, R. R., & Stone, J. M. 2013a, ApJ, 770, 67Belyaev, M. A., Rafikov, R. R., & Stone, J. M. 2013b, ApJ, 770, 68Bertout, C. & Regev, O. 1992, ApJ, 399, L163Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (PrincetonUniversity Press)Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46,253Cohen, M. & Kuhi, L. V. 1979, ApJS, 41, 743Drazin, P. G. & Reid, W. H. 2004, Hydrodynamic Stability
Article number, page 17 of 18 & A proofs: manuscript no. aa26005-15