Turbulent convection in protoplanetary discs and its role in angular momentum transfer

We present a model for the transport of anisotropic turbulence in an accretion disc. The model uses the Reynolds stress tensor approach in the mean field approximation. To study the role of convection in a protoplanetary disc, we combine the turbulence model with a radiative transfer calculation, and also include convection using the mixing length approximation. We find that the turbulence generated by convection causes the angular momentum of the accretion disc to be directed outwards. We also confirm the conclusions of other authors that turbulent convection is unable to provide the observed disc accretion rates as well as a heat source sufficient for the convection to be self-sustaining. The reasons for the latter are the strong anisotropy of the turbulence together with the low efficiency of the energy transfer from the background velocity shear to the turbulent stress tensor.


INTRODUCTION
The theory of disc accretion is used in astrophysics to explain a wide range of observed sources and phenomena: active galactic nuclei, the evolution and variability of close binary systems, the formation of jets and bipolar outflows, the structure of protoplanetary discs, the formation of planetary systems, and many others, see for example Shakura (2018); Hartmann (2009); Armitage (2015).In all these objects, accretion takes place under different physical conditions, varying in temperature, density, degree of ionisation, magnetic induction, radiation field, presence of dust, and so on.However, the processes that influence accretion have common consequences: this is the redistribution and removal of angular momentum, which allows matter to accrete from the disc to the central object.
Turbulence is thought to play an important role in many accretion processes.The classical approach to describe turbulent accretion relies on the formalism of turbulent viscosity, which is mathematically equivalent to molecular viscosity.In the Shakura & Sunyaev (1973) model, this formalism is reduced to the setting of an alpha parameter that relates the coefficient of turbulent viscosity to the speed of sound and the scale height of the disc.This phenomenological approach has proved extremely useful for describing the structure and evolution of astrophysical discs, but the question ⋆ E-mail: kurbatov@inasan.ru of the causes and properties of the turbulence itself remains beyond its scope.
Currently, one of the most active areas of astrophysical research is the study of protoplanetary discs (PPDs) around young stars.This interest is stimulated by the progress in observational techniques providing the means to obtain direct images of the discs at different wavelengths (see the review by Andrews 2020).The high angular and spectral resolution allows, among other things, the reconstruction of the detailed distribution of the turbulent gas velocity across the disc, with the estimates varying greatly for different sources (see Flaherty et al. 2017;Guilloteau et al. 2012).In this context, protoplanetary discs can be considered as a convenient natural laboratory for studying the physics of accretion and turbulence in general.
It is widely accepted that turbulence in the disc arises due to some instability.In protoplanetary discs, the possible triggers of turbulence could be gravitational, thermal, magneto-rotational, baroclinic, streaming, vertical shear and other instabilities (see the reviews by Armitage 2015;Bae et al. 2022;Lesur et al. 2022).These instabilities appear at different dynamical and thermodynamical conditions in the disc, (see, e.g., Pfeil & Klahr 2019).Each instability is the subject of extensive research.For example, in Klahr & Hubbard (2014) the authors considered the effects of radial buoyancy in discs.It was shown that in a rotating flow, radial buoyancy together with centrifugal force can cause epicyclic oscillations with increasing amplitude (Latter 2016;Volponi 2016), the phenomenon named convective oversta-bility.In the nonlinear regime, this instability can lead to a subcritical baroclinic instability (Lyra 2014), as well as to the growth of large-scale vortices, which may play a role in planet formation (Raettig et al. 2021).
In the present work we are interested in the convective instability.The link between convection, turbulence and angular momentum redistribution has been investigated in many studies.The idea that convection in protoplanetary discs can not only transfer heat but also provide viscosity and thus influence the evolution of the disc was formulated by Cameron (1978) and Lin & Papaloizou (1980).This idea has generated a lot of interest, but after several decades of research the role of convection in the transfer of angular momentum is still controversial, see a detailed historical review in Klahr (2007), and also the recent papers by Held & Latter (2018, 2021).A representative example is that in early numerical models, convection was found to cause the transfer of angular momentum towards the accretor (Stone & Balbus 1996), which would correspond to a negative alpha parameter.In later work, using high-resolution numerical schemes, it was shown that the angular momentum of the accreting matter is transferred outwards (see Held & Latter 2021, and discussion therein).Held & Latter (2018) presented the results of 3D modelling of convection in a disc, illustrating the emergence of convective cells, eddies and other coherent structures upon the initiation of convection.At the same time, they noted that they could not obtain a self-sustaining convection regime in the disc.Held & Latter (2021) showed that the interaction of convective and magneto-rotational instabilities ensures the periodic nature of accretion in the disc.Pavlyuchenkov et al. (2020) and Maksimova et al. (2020) also showed that convective instability in a protoplanetary disc can lead to irregular accretion onto a star.This result is relevant in the context of the search for physical mechanisms to confirm the scenario of episodic accretion in protoplanetary discs (Hartmann 2009), which is important for solving the problem of observed accretion luminosities and for explaining the nature of young stellar objects with luminosity outbursts, such as FU Ori and EX Lup type stars.However, the key approximation of the model presented in Pavlyuchenkov et al. (2020) is the assumption that the emerging convection is accompanied by high turbulent viscosity.
The relationship between convection, turbulence and accretion in the disc layer is pictured in Fig. 1, which is based on the energy circulation scheme.Let us suppose that there is an initial heating source in the medium.Thermal energy is transferred by radiative diffusion and is eventually emitted as infrared radiation.Certain conditions can lead to the development of convection, which not only transfers some of the heat, but also excites the turbulence.The dissipation of the turbulence eventually converts the kinetic energy back into thermal energy.In addition to this cycle, the turbulence may be intensified by the background shear flow, which dissipates its kinetic energy and replenishes the heat budget, providing another energy source for convection.This way, the turbulence converts the gravitational energy of the gas into the thermal energy.In this scheme, the fundamental question is how significant is the contribution of the energy of the differential rotation to the turbulence strength.vection, but its underlying physical mechanisms may vary.
As an example, one can suggest the heating due to density waves excited by the disc self-gravity (Cossins, Lodato, & Clarke 2009).Another possible sources are the magnetorotational instability (Held & Latter 2018) or dissipation of large scale magnetic field (Béthune & Latter 2020).Cosmic rays, which can penetrate quite deep into the disc, can also be considered as a source of heating (D'Alessio et al. 1998).In addition, there are factors that directly prevent convection.For example, the heating of the disc by stellar and interstellar radiation helps to establish a positive temperature gradient in the upper layers of the disc, making it stable against convection.
Although direct numerical simulations can produce very realistic results, a complete self-consistent three-dimensional calculation of the accretion disc evolution with a sufficiently high spatial and temporal resolution remains a challenging task.Even if the full-scale numerical model is assumed to be sufficiently resolved and accurate, there remains the problem of interpretation, i.e. the assessment of the importance of one or another physical factor that influences the gas dynamics, transport processes, etc.In this paper, we implement a non-isotropic turbulent transport model (in the mean-field approximation), together with a calculation of the convective flow (in the mixing-length approximation) and radiative transport to study the problem of turbulent convection and its role in the redistribution of angular momentum.The mean-field approach makes the calculations simpler than the full three-dimensional hydrodynamic calculations, while allowing the most important physical processes to be emphasized.Thus we can explicitly determine the contribution of different factors involved in the energy cycle to the dynamics of turbulence.
We base our modelling on the mean field turbulence model proposed by Canuto (1992Canuto ( , 1993Canuto ( , 1997) ) for stellar atmospheres.This model relies on the momentum representation to describe the fields of velocity, density and pressure fluctuations up to the fourth order moments.It is quite complex and has never been fully implemented for astrophysical applications.We formulate a reduced version of the Canuto model, where only the dynamics of the turbulent stress tensor is calculated explicitly, while the remaining closures are implemented in the gradient or algebraic approximations.As in the original model of Canuto, the only seed of the turbulence is the convective flux.However, the turbulence can grow or decay due to interaction with the background shear flow.The second component of our model is infrared (IR) radiative transfer for PPDs by Pavlyuchenkov et al. (2020).It uses temperature-dependent opacities for a mixture of graphite and silicate dust grains, which makes it possible to realistically simulate the conditions for the development of the convective instability.This model also takes into account the absorption of the radiation from the central star and the interstellar medium.
In Section 2, we present the mean-field turbulence transfer equations, as well as the closures and the convective flux.In Section 3, the full turbulent convection and the radiative transfer model is formulated in the one-dimensional cylindrical frame.We also perform the test calculations and implement the model for a vertical column in a protoplanetary disc.Discussion and conclusions are presented in Sections 4 and 5, respectively.

Mean field and turbulence transfer equations
As mentioned in the introduction, in this paper we use the mean-field approach to model the turbulence.In this approach, we can distinguish two methods for the description of the turbulence: the filter method and the statistical method.The former uses a spatial and temporal filter, then the details of the sub-scale flow are only approximated by averaging on the filter scale.This approach is implemented in the class of subgrid models and in the Large Eddy Simulation (LES, Leonard 1975;Meneveau et al. 1996).Another approach is to introduce a statistical ensemble for the turbulent fluctuations, assuming that they are stochastic.The properties of the fluctuations are then formulated in terms of the statistical moments of this ensemble (see e.g. the references in Canuto 1997).Since there is only one realisation of the flow in any given problem, the LES method may seem more physically justified.In addition, this method can explicitly describe non-local effects, such as the interaction of subgrid and supergrid scale structures (Leonard 1975;Stewart 1976).The statistical approach, in turn, greatly simplifies the computations, as it allows the use of empirical information on the amplitudes of the fluctuations and their mutual correlations.We are interested in the effect of the turbulence on the mean flow rather than in the detailed spatial and temporal structure of the turbulence, so we describe its properties in terms of the statistical moments of the velocity, density and pressure fluctuations.
We will only provide the final expressions for the turbulent transfer model of Canuto (1992Canuto ( , 1997) (see also references therein), without the detailed derivation.We write the equations in arbitrary curvilinear coordinates with the metric tensor κ ij .Later on, the model will be implemented in cylindrical frame.Let us define the variables that characterise turbulent flow: volume density ρ, velocity v i , pressure p, internal energy e, and the Reynolds stress tensor wij.These quantities satisfy the dynamic equations: (4) In addition to the quantities listed above, these equations also include gravitational acceleration g i , heat source q, convective flow F i conv , and several closures to the equation for the Reynolds stress tensor (w ijk , B ij , Π ij , and ϵ), which will be defined later.
Thermodynamical variables are related to each other via the ideal gas equation of state: where T is the temperature; µ is the weight of a gas particle in hydrogen atom mass units, mH; R = kB/mH = 8.25 × 10 7 erg g −1 K −1 is the gas constant; γ is the adiabatic index; cp and cv are the specific heat at constant pressure and volume, respectively.

Closures
The closures in the r.h.s. of the Eqs.( 3) and ( 4) are 2nd and 3rd-order statistical moments of turbulent fluctuations.It is also possible to formulate dynamic equations for these moments.However, since the model already contains a large number of parameters, we will use algebraic closures.Making the model more complex will only make the results more difficult to interpret.A valuable quantity in the algebraic closures is the turbulence correlation time.Under accretion disc conditions, Keplerian time appears to be a natural time scale for turbulence correlation.An estimate obtained by Stewart (1976) by analysing an equation similar to Eq. (2) leads to an expression for the correlation time of the form where |Ω| is the angular velocity of the gas rotation; MT is a turbulent Mach number, it depends on the local speed of sound c 2 s : Mach numbers estimated from the non-thermal broadening of spectral lines in protoplanetary discs range from 0.06 (Flaherty et al. 2017) to 0.5 (Guilloteau et al. 2012), which gives estimates for turbulence correlation time within 0.3 ≲ |Ω| tT/(2π) ≲ 2.7.The B ij tensor is responsible for buoyancy effects.It can be expressed in terms of the convective flux: The buoyancy tensor is the seed of turbulence in the present model.Note that the r.h.s of Eq. ( 3) contain the convective source B k k .In an accretion disc, one can expect the pressure gradient to be directed towards the disc mid-plane, while the convective flow is directed away from it.Hence we can conclude that B k k ⩾ 0, which means that convection takes away the thermal energy of the mean flow.
It is widely accepted that pressure fluctuations play a key role in the development of turbulence.In the Π ij tensor these effects manifest themselves in isotropising the turbulence ("return-to-isotropy"), buoyancy and interacting with the background flow.The following form of this tensor was derived from symmetry and dimension considerations (see references in Launder 1974;Speziale 1991;Canuto 1997): where b ij characterises the deviation from isotropy, U ij and V ij are the symmetric and asymmetric parts of the strain rate tensor of the background flow, respectively, Note that the expression (10) depends linearly on the components of the Reynolds tensor, except for the last term.
The variable ϵ is related to the dissipation of turbulent energy.This process takes place on small scales where the anisotropy is on average weak.For this reason, we use the classical isotropic closure and K = (1/2) w k k is the turbulent kinetic energy volume density.
To write out a closed expression for the third-order velocity momentum w ijk , we use the gradient approximation: The factor νT is interpreted as the turbulent kinematic viscosity coefficient.It can be written as follows The dimensionless factor Cν is usually set equal to 0.09 (Launder 1974;Speziale 1991).As can be seen, the tensor w ijk describes the diffusion of turbulence.
The constants in the expressions ( 10) and ( 16) were obtained experimentally (see papers by Launder 1974;Speziale 1991;Canuto 1992Canuto , 1993, and references therein): Finally, let us consider the gas energy balance.Taking the trace of Eq. ( 4), one can get Here one can see that the B k k source enters the r.h.s. with a positive sign.By the comparison with Eq. ( 3), we can conclude that convection converts the thermal energy of the mean flow into the turbulent energy.
Projecting Eq. ( 2) onto the velocity vector gives the equation for the kinetic energy of the mean flow: Combining Eqs. ( 3), ( 18) and ( 19), we get the equation for the total energy: This equation is completely conservative except for external sources.This means that the thermal, kinetic and turbulent energy of the medium can only transform to each other.

Convective flux and conditions for instability
In Mixing Length Theory (MLT), heat is transferred by convective elements formed by convective instability.It is usually assumed that the elements move with a characteristic velocity vconv under the effect of the buoyancy force and that they transfer the excess heat ρcp∆T to the surroundings.The convective flux can be written as ρcp∆T vconv.Different versions of the theory differ in the way the quantities ∆T and vconv are estimated.In this way, the radiative heat losses of the convective element along its path, its viscous deceleration, and the spreading of convection beyond the convective zone (overshooting) (Canuto 1992) can be considered.
The conditions for convective instability are fulfilled in the regions where the temperature gradient exceeds the adiabatic (more precisely, the isentropic) gradient in the direction opposite to the gravitational acceleration.Let us denote the excess temperature gradient as where n i is the unit vector in the direction of the gravitational acceleration g i (note that the centrifugal acceleration also contributes to the value of g i ).The adiabatic gradient is also expressed in terms of the acceleration vector, The instability condition is −g k β k > 0.
Let us write down the results given in Hansen & Turbulent convection in protoplanetary discs 5 Kawaler (1994) for the problem of stellar convection, without going into a detailed derivation.The characteristic parameter of the theory is the length of the mixing path ℓ.This is the distance traveled by the convective element before it mixes with the surrounding matter.Another characteristic parameter is the growth rate of the convective instability (Hansen & Kawaler 1994, Chap. 5), where ν mol and ν rad are the molecular (collisional) and radiative thermometric conductivities, respectively, cm 2 /s cm 2 s −1 ; N is the Brunt-Väisälä frequency, it is defined as In a convectively stable medium (g k β k > 0), N is the frequency at which the gas element oscillates due to buoyancy and gravity.In the limit of weak convection, |N | ≪ (ν mol +ν rad )/(2ℓ 2 ), the increment is ω ≈ ℓ 2 |N | 2 /(ν mol +ν rad ).
In the opposite limit the increment is saturated as ω ≈ |N |.
The collisional thermometric conductivity coefficient for neutral atoms is where v th = (3RT /µ) 1/2 is the mean thermal velocity of the molecules; σnn = 3 × 10 −16 cm 2 is the collision cross section for neutral hydrogen.The radiative thermometric conductivity coefficient has the form here c is the speed of light; a rad = 7.56×10 −15 ergcm −3 K −4 is the radiation density constant; κR is the Rosseland mean opacity.It should be noted that in many astrophysical applications, the collisional mechanism of the heat conduction can be neglected.The speed of the convective element and the excess temperature are estimated as follows: where |β| is the magnitude of the vector β i , Eq. ( 21).As a result, the convective energy flux takes the form The MLT has one free parameter, the mixing path length ℓ.The convective flux ( 29) is very sensitive to this parameter: from ℓ 8 in weak convection to ℓ 2 in strong convection.The pressure scale height of a star is usually adopted as the mixing length in stellar convective shells models.Under the accretion disc conditions, the thermal scale height of the disc can be taken as the mixing length.

Final system of equations
Consider a cylindrical coordinate system 1 (r, ϕ, z) and a rotational axisymmetric flow in a narrow radial annulus of radius r.Our requirement is that all quantities are independent of the azimuthal angle.It will be assumed that the disc is close to a mechanical equilibrium.In this case, the gradients can be estimated as where H is the vertical thermal scale in the disc; Ω = v ϕ /r is the angular velocity of the gas rotation.If no luminosity outbursts are considered in the disc, then the radial and vertical velocities of the gas in typical accreting discs become essentially subsonic: where M ≪ 1 is the Mach number for the radial and vertical background gas velocities.In this approximation the advection terms in the expressions ( 2)-( 4) can be neglected, except for the centrifugal force in the Euler equation, as well as the corresponding components in the turbulence transport equation.For all quantities except the background gas angular velocity, the radial dependence is neglected.We also neglect the self-gravity of the disc and consider only the gravity of the star.Finally, the system of equations ( 2)-( 4) takes the following form: In Eq. ( 33), ΩK is the Keplerian angular velocity at the radius r; νT is the turbulent viscosity coefficient ( 16).Due to the chosen approximation, only the z-component of the convective flux vector (29) remains non-zero providing the only non-zero component of the buoyancy tensor: Fconv . (41) The components of the isotropisation tensor are now where ( As can be seen from the expressions ( 35)-( 40), the convective heat flux is the seed of turbulence.The background flow is only involved in the amplification or weakening of different wij components.For the MLT flux (29) we take the mixing length to be equal to the vertical thermal scale of the disc, ℓ ≡ H = cs/|Ω|.
There are several sources on the r.h.s. of the heat balance equation ( 34): the first one is responsible for convective heat transfer, the second one describes the consumption of thermal energy for convective motions, the third one provides energy input due to turbulence dissipation.The last source, q, accounts for heating by stellar and interstellar radiation, exchanging energy with its own IR radiation, and may also include an additional heat source.The model for q is given in Section 3.3.
The Eqs. ( 35)-( 40) require boundary conditions.At the upper boundary of the disc, it is natural to set each of wij to zero.In the disc mid-plane (z = 0), the boundary conditions are as follows: It can be seen that in the system ( 35)-( 40), the last two equations do not contain any energy source or sink and only govern the redistribution of the turbulent energy between the components wrz and w ϕz .Given the boundary conditions, this means that if these components are initially zero, they will remain zero in the future.Thus, the equations for the wrz and w ϕz are not considered further.

Testing the model of turbulent transfer
Here we will briefly analyse the turbulent transfer model and try to reveal the role of the free parameters CΠ1-CΠ3 defined in Eq. ( 17).From Eqs. ( 42)-( 45) one can see that CΠ1/tT is the inverse characteristic time of Reynolds stress tensor isotropisation.The parameters CΠ2 and CΠ3, in turn, are the coupling constants between the turbulence and the background shear flow (its symmetric and asymmetric parts, respectively, see Eqs. ( 12) and ( 13)).Finally, CBBzz/2 is the heat loss per unit time due to the excitation of turbulence by the convection channel.
Let us convert the field variables to the dimensionless form: We simplify the model by neglecting turbulent diffusion, which is applicable to the conditions deep within the convective zone.In this case the equations ( 35)-( 38) are reduced to a system of ODEs.The local background shear velocity profile is assumed to be in general form, rather than Keplerian, Ω ∝ r −q .It is important to note that Ω is the projection of the angular velocity vector onto the OZ axis.It can be shown that the w rϕ component always enters the equations in a combination sign(Ω) w rϕ .Thus, the inversion of the angular velocity leads to a change of the sign of w rϕ (and also of wrz and w ϕz , see Eqs. ( 38)-( 40) and the corresponding equations for Πij).Without loss of generality we assume Ω > 0.
It is useful to regroup the dimensionless equations to the following representation: It is seen here that the component wrϕ is coupled to wrr + wϕϕ by CΠ2 and to wrr − wϕϕ by CΠ3.
In a steady state limit, the energy equation gives The wrϕ component is therefore determined by the difference between the amount of dissipated turbulent energy and the energy supplied by the convective source.Despite that convection is the seed of energy for turbulence in our model, the background velocity shear is also important.Once the convection is able to produce the turbulence, the velocity shear starts to amplify the components of the Reynolds stress tensor (through the coupling constants CΠ2 and CΠ3), making the r.h.s. of the Eq. ( 65) non-zero.Since the wrϕ component is responsible for the angular momentum transfer in the disc, the efficiency of the transfer is clearly related to the energy balance.When there is no background shear, q = 0, the parameter wrϕ is decoupled from the turbulent energy K, so there is no transfer of angular momentum.
Let us assume wrr = wϕϕ .Then it can be shown from Eqs. ( 61)-( 63) that in the steady-state limit all the turbulence components are zero.Since this derivation is independent of the source of the turbulence (in r.h.s. of the Eq. ( 64)), this is true not only for convection-generated turbulence, but also in the general case of rotational shear flows.The latter means that turbulence is always anisotropic in accretion discs (at least in the no-diffusion approximation).
We have performed calculations of the model ( 61)-( 64) with a fixed convective source Bzz = 0.01 and fixed values of the constants CB, CΠ1, CΠ2 and CΠ3 from the Eq. ( 17).Fig. 2 shows that after a monotonic growth phase lasting 1-3 disc periods, the turbulence reaches a steady state.Varying the velocity profile index q by 10%-15% significantly affects the turbulence intensity, in particular the wrϕ component changes by a factor of two.Simulations with high velocity profile indices revealed that q wrϕ ∼ K/τT ∼ K for q ≳ 3 (not shown in Fig. 2).Negative values of q leads to the negative values of wrϕ , though the dependence on q is weaker.
It is interesting to see how the choice of constant parameter values affects turbulence.For random sets of the constant parameters and various velocity profile indices, we ran a ensemble of test calculations, starting from zero initial conditions.The convective source was fixed to Bzz = 0.01 as it only determines the value of the turbulence energy in the steady state limit, see Eq. ( 65).The constants were chosen randomly from the intervals 0 ≤ CΠ1 ≤ 5, 0 ≤ CΠ2 ≤ 2, and 0 ≤ CΠ3 ≤ 2 (blue dots in Fig. 3).The velocity profile index q was set to 1.5.Two qualitative indicators of the solutions are shown in Fig. 3.The first indicator, min{ wrr, wϕϕ , wzz} declares physical constraints: the quadratic velocity correlators should not be negative.Solutions that satisfy this constraint leave this indicator at zero value.It is seen that the physically allowed values of CΠ1 cannot be lower than ∼ 1.The allowed values of CΠ2 are bounded in a quite narrow range around the experimental value (17).Varying CΠ3 within the considered limits does not violate the physical constraints.The second indicator in Fig. 3, min wrf , shows under which conditions the direction of the angular momentum flux changes.Note that only the physically allowed solutions are shown here, i.e. the solutions with min{ wrr, wϕϕ , wzz} = 0.The sign of the off-diagonal component wrϕ is not sensitive to CΠ1 and CΠ3, however, it is sensitive to CΠ2. Surprisingly, the experimental value of ).Different values of the velocity profile indices were explored: q = 1.5 (solid curves), q = 1.65 (dashed), q = 1.25 (dot-dashed).
CΠ2 from Eq. ( 17) only slightly exceeds the lower limit of the range where min wrϕ ≥ 0. The above boundaries have been estimated from the considered ensemble runs and are shown in more detail in Fig. 4. Fig. 5 shows the effect of the outlying values of the constant parameters on the qualitative behavior of the model.Assigning CΠ1 a value below the acceptable range leads to the negative steady-state limit in wrr and wϕϕ .On the other hand, the assignment of CΠ2 to the outlying value results in oscillating (including negative) solutions but positive limits.Interestingly, in these calculations the turbulence energy is always positive and quite stable.

Combining the turbulent convection model with the radiative transfer model
We calculate radiative transfer and heat balance using the thermal model of Vorobyov & Pavlyuchenkov (2017) and Pavlyuchenkov et al. (2020).The model considers the energy exchange between the gas and the IR radiation, the UV heating by stellar and interstellar radiation (SUV), and an additional heat source (Sext): where E rad is the energy volume density of the IR radiation; κP is the Planck mean absorption coefficient.The Sext source can be associated with some dissipation processes that are not explicitly included in our model (see below in this section).The heating due to the stellar and interstellar irradiation, SUV, is calculated by a direct integration of the radiative transfer (Vorobyov & Pavlyuchenkov 2017).The IR radiative transfer is simulated by means of the Eddington approximation, which is reduced to the following system of equations: where F rad is the IR radiation flux; κR is the Rosseland mean absorption coefficient.The radiation flux boundary conditions in the mid-plane and at the upper boundary zmax are of the form where TCMB = 2.73 K.An important feature of the disc thermal model is the use of temperature dependent Rosseland and Planck mean opacities, as it has been found that an increase in opacity with temperature is necessary for the onset of convection (Lin & Papaloizou 1980).The opacities have  17).Coloured lines: One of the constants is varied (see the legend in the first plot).
been taken from Pavlyuchenkov et al. (2020) where they were obtained for a mixture of graphite and silicate dust grains.The solution of the subsystem of the radiative transfer equations ( 66)-( 68) is found by an implicit method.The details and tests are described in the appendix of Vorobyov & Pavlyuchenkov (2017).
We noted above that in our model, convection is the seed of turbulence.To initiate convection, an external heat source is required.It is included in the heat energy equation as Sext in Eq. ( 66).We define the external energy injection rate by using the accretion rate Ṁ as a measure (Shakura & Sunyaev 1973): where Σ is the local surface density of the disc.
Figure 6 shows the thermal structure of a protoplanetary disc obtained in the 1+1-dimensional framework using the radiation transfer model without convection and turbulence.By 1+1 formalism, we mean that self-consistent density and temperature distributions in vertical direction (the first '1'-dimension) are obtained separately for each radial position (the second '1'-dimension ) of the disk.In our calculations, the surface density distribution and stellar heating depend on the radial position, but the columns do not affect each other.The purpose of this calculation was to locate the strongest convective instability.This model includes all the heating sources from the one-dimensional model above, along with the external heating Sext.The parameters of the model are: the stellar mass 1 M⊙, the star effective temperature 5780 K, the accretion rate onto the star 10 −7 M⊙ yr −1 , the density 10 3 g cm −2 at a radius of 1 au, assuming a power law of the surface density Σ ∝ r −1 .The convective flow was included in this calculation, but was not involved in the formation of the disc thermal structure.The last two plots in the Fig. 6 show the region of convective instability (∇ − ∇ ad > 0 2 ).The convection is concentrated in a rather shallow torus near the mid-plane (not reaching it due to symmetry constraints).The maximum value of the convective flux is reached in the inner part of the computational domain, r ≈ 1 au, and drops by three orders of magnitude at r ≈ 3 au.Note that the flux distribution has two peaks, (cf.with Pfeil & Klahr 2019, Fig. 10).
We have performed several runs considering different radial positions of the disc column.We have chosen the column at a radius of 3.34 au for two reasons.First, as the average temperature of the disc decreases with distance from the star, the role of radiative heat transfer within the convective zone is reduced in favour of convective heat transfer.Second, the thickness of the convective zone decreases with radial distance.This choice allows us to see the effect of turbulent diffusion in detail.

Numerical solution scheme
The solution of the full system of equations ( 33)-(40) for each time step was divided in two stages.In the first stage, the heat balance, radiative transfer and hydrostatic equations were solved jointly with the given sources SUV, Sext and wij.Within a time step, the heat balance equation was linearised in temperature and the hydrostatic equation was linearised in density.In the radiative transfer equation, the spatial derivative operator was discretised using a standard scheme and expressed by a tridiagonal matrix.This system of equations was solved using a completely implicit iterative scheme.The detailed scheme of the solution is described in Vorobyov & Pavlyuchenkov (2017) and Pavlyuchenkov et al. (2020).
In the second stage, the turbulence transfer equations were solved.This was done by discretisation the spatial derivative operator and solving the entire system of equations on the spatial grid as a system of ODEs in time.The solution was performed using an explicit-implicit scheme LSODA.3

Disc column with turbulent convection
The vertical structure of the disc column was modelled in the one-dimensional turbulent disc convection approach.We set the values of the external heat source Sext by parameterising it with the accretion rate, Eq. ( 71): Ṁ = 10 −7 M⊙ yr −1 (Model A) and Ṁ = 10 −4 M⊙ yr −1 (Model B).These options correspond to a quiescent and outburst state of the disc, respectively (Audard et al. 2014;Fischer et al. 2022).
Previously, we ran a 1+1-dimensional disc model in order to find the location of the convective zone and to select the most interesting conditions for a full calculations including turbulent convection.However, that model did not take the convective heat transfer into account.As a result, the density and temperature distributions were unstable and therefore not physical.As will be seen later, the account of the convective heat transfer significantly changes the steady state of the disc.We start the full model with the same initial conditions as in the previous model, namely assuming an isothermal disc with a temperature of 100 K.The main parameters of the model are listed in Table 1.Calculations were continued until steady states were reached.In total, Model A was run for 127 Keplerian orbits, and the Model B was run for eight orbits.We also run each model without convection (hence, without turbulence) to have a baseline to compare our convective models to.The model without convection consists of the hydrostatic (33), heat transfer (34), and radiative transfer (66)-( 68) equations, only: all the convection and turbulence terms (wzz, Fconv, Bzz, and ϵ) have been omitted in this model.Using the proposed model, it is interesting to evaluate how effectively the convective turbulence utilizes external heating and background shear flow.It is also important to investigate whether the energy cycle (Fig. 1) can be selfsustaining without external heating, leading to steady convection and turbulence.The purpose of this modelling is to quantify the effective accretion rate and the dissipation rate of turbulence.

Model
The results for Model A are shown in Fig. 7.As can be seen, the convective zone extends upward to approximately one and a half thermal scale heights.In a significant part of the convective zone, the convective flow turns out to be comparable in magnitude to the radiative flow.Thus, convection provides about half of the total heat flow.The total energy flux is only a few percent higher than the energy flux in the calculation without convection (dashed lines).This, however, is enough to reduce the temperature of the inner layers of the disc, within the thermal scale, by about 25 K.In our model, the convective flow is directed along the OZ axis, so the convection only excites the wzz component of the Reynolds tensor directly.The remaining components of this tensor arise from wzz due to the pressure tensor Πij, then amplify through interaction with the background flow.Note that wzz dominates in magnitude, it exceeds wrr and w ϕϕ by 3-5 times and exceeds w rϕ by two or more orders of magnitude.The intensity of turbulence and the corresponding contribution to the width of the observed spectral lines is determined by the diagonal sum of the Reynolds tensor, k w kk .Due to the turbulent diffusion, the region of developed turbulence is almost twice as thick as the convective zone.In the bulk of the disc material, the turbulent Mach number MT varies in the range 0.03-0.08,which is in good agreement with the estimates obtained from observations (Flaherty et al. 2017).
The rate of angular momentum removal due to the turbulence and the corresponding accretion rate (which we will call the effective accretion rate) depend on the off-diagonal component w rϕ (Shakura & Sunyaev 1973): The effective accretion rate can also be estimated in terms of the alpha parameter: due to the fact that w rϕ ≪ k w kk (see Fig. 7 and the Table 2), one cannot draw a correct conclusion about the accretion rate from the turbulent kinetic energy, or, equivalently, from the width of the spectral lines.In the Model A, the effective accretion rate is Ṁeff = 1.7 × 10 −10 M⊙ yr −1 , which is 1000 times less than the accretion rate Ṁ that determines the heating of the layer.
In the inner layers of the disc, z < 0.5 au, the external source ρSext dominates among all sources, except for a narrow region near the top boundary of the convective zone, where the convective source ρSconv is more important.Heating from the convective source is comparable to external heating in absolute value, but its role is to redistribute and slightly reduce thermal energy.The source ρS turb associated with turbulence dissipation is smaller than ρSext by one or two orders of magnitude everywhere in the disc.The integrated values of each of the sources (Q = 2 zmax 0 ρS dz) are given in Table 2. Also note that in the Model A even more energy is spent to start convection than is returned to the heat budget from dissipation: |Qconv| > Q turb .

Model B (
The Model B results are shown in Fig. 8. Steady state was reached in eight Keplerian orbits, i.e. much faster than in the Model A. This is due to higher temperature of the gas in the Model B (∼ 2000 K) comparing to the Model A (∼ 200 K), hence the shorter thermal time scale.The radiative time scale may be estimated using a well known approximation derived for an optically thick layer, see e.g.Wu & Lithwick (2021): Since κR ∝ T (Pavlyuchenkov et al. 2020), an order of magnitude increase in temperature results in two order of magnitude reduction of the characteristic radiation time.
In case of higher heating power Qext, which corresponds to the accretion rate Ṁ = 10 −4 M⊙ yr −1 , the convection is also developed.The magnitude of the convective flow in the Model B is much higher than in the Model A, but is much lower than the radiative flow and does not affect the thermal structure of the layer.The maximum amplitude of turbulent fluctuations, which is reached near the upper boundary of the vertical thermal scale, exceeds 0.2 cs.However, the w rϕ tensor component responsible for the transfer of angular momentum is two orders of magnitude lower than the wzz component, as in the quiescent model.In the Model A Ṁeff / Ṁ = 1.7 × 10 −3 , while in the Model B Ṁeff / Ṁ = 1.7 × 10 −4 .

DISCUSSION
In this paper we have used the idea that convection in protoplanetary discs is turbulent.In general, this is not necessarily the case.In order to quantify the transition of the convective flow between the laminar and turbulent regimes, let us use an empirical criterion based on the Reynolds number Re = V L/ν mol , where V is the flow velocity, L is its characteristic scale, and ν mol is the molecular kinematic viscosity 3.5 × 10 −3 7.9 × 10 −3 α eff 1.6 × 10 −5 1.5 × 10 −4 coefficient.When Re exceeds a critical value Recr, the flow becomes turbulent.Different values of the critical Reynolds number correspond to different types of flow: from tens for rotational flows to 10 3 and even 10 5 in the special cases of the flow in a tube (Landau & Lifshitz 1959).We can apply these empirical considerations to the parameters of the convective flow.Let us determine V as the convective element velocity ωℓ, Eq. ( 27), and L to be equal to the mixing length ℓ.After substitution of the weak convection limit for ω, see Eq. ( 23), the Reynolds number becomes This expression is actually the Rayleigh number Ra (Canuto 1992;Held & Latter 2021).We may therefore speculate that for the laminar-turbulent transition of the convective flow, the critical Rayleigh number should be of the order of the critical Reynolds number.Typically we had Ra ≳ 10 9 in our calculations, see Figures 7 and 8, which is high enough for turbulence to develop.The efficiency of convection in driving accretion has long been debated.At one point, a number of papers appeared which argued that the convection that develops in an accretion disc could not lead to an outward transfer of angular momentum.Stone & Balbus (1996) proposed some analytical arguments to support the idea that in the accretion disc with convective turbulence the angular momentum flux is directed inward.Their argument was based on the assumption that angular variations of pressure in axisymmetric turbulent flows are subdominant compared to the radial variations (see their comments below the Eq. ( 11)).This is not the case in our model, since the radial and angular components of the isotropisation tensor Πij are comparable in magnitude.Previously, Held & Latter (2021) noted that in the high-resolution numerical model the hydrodynamic convection can transport angular momentum outwards (Lesur & Ogilvie 2010;Held & Latter 2018).
The models based on the mean-field approach inevitably depend on free parameters.Given the conventional, experiment-based values of the parameters, Eq. ( 17), the angular momentum in our model is transferred outwards.The numerical experiments with a local model (Sec.3.2) have shown that the direction of the angular momentum trans- port depends crucially on the value of the parameter CΠ2.The assumed value of CΠ2 is quite close to a critical below which the angular momentum flux changes its sign.
Our motivation for studying convection in protoplanetary discs was particularly inspired by the idea that turbulent convection could arrange irregular accretion from the disc to a star, as proposed by Pavlyuchenkov et al. (2020); Maksimova et al. (2020).The importance of convection has also been mentioned in other studies, e.g. by Hirose (2015); Held & Latter (2021), when considering high accretion states of MRI active protoplanetary discs.However, our results indicate that turbulent convection is a weak mechanism for angular momentum transfer in the protoplanetary disc.In fact, we support the results of Lesur & Ogilvie (2010) and Held & Latter (2018) that the turbulence generated by convection does not provide the observed disc accretion rates and sufficient heat influx for convection to be self-sustaining.There are two reasons for this: the anisotropy of the turbulence, and the fact that convection is too weak a source of turbulence.
The first reason is that the wzz element of the isotropic part of the Reynolds tensor is the only element excited by convection, while w rϕ is the only element responsible for the removal of angular momentum from the disc.The energy exchange between Reynolds stress tensor elements is not efficient enough, so w rϕ is more than two orders of magnitude smaller than wzz.
To estimate the importance of the second factor, we can look at the turbulence in the steady state near the maximum of wij in Fig. 7 or 8 (bottom row, left).Under these conditions, the diffusion term disappears from the equations ( 35)-( 38), while the system of linear algebraic equations remains.The only inhomogeneous term in this system is Bzz.Thus, the solution of the resulting system is proportional to the value of this convective source, see Section 3.2.To make the effective accretion rate Ṁeff formally equal to the given Ṁ , it is necessary to increase the convective flux by three orders of magnitude in the Model A and by four orders of magnitude in the Model B. This is hardly possible, even considering the uncertainty of the mixing length ℓ.
We note that convection might still play an important role in facilitating angular momentum transport if turbulence is excited not by convection alone, but by the collective effects of different instabilities.For example, the papers by Hirose et al. (2014); Coleman et al. (2018); Scepi et al. (2018) and Held & Latter (2021) present calculations of 3D MHD models showing that the joint action of convection and magneto-rotational instability can increase α to the observed values.
It would be interesting to apply our mean field model to other instabilities, such as vertical shear instability, streaming instability, etc. (see the reviews by Bae et al. 2022;Lesur et al. 2022).For example, Stoll et al. (2017) used numerical simulations to show that vertical shear instability leads to the appearance of anisotropic turbulence.

CONCLUSIONS
In this study, we have presented a model for the transport of anisotropic turbulence in a protoplanetary disc.The model includes time-dependent heat transfer by radiative diffusion and convection, developing on the background of hydrostatic equilibrium.The time-dependent turbulent transport model is based on the mean-field approach formulated in terms of Reynolds stresses.The seed of turbulence in our model is the convective instability, hence the convection flux.In addition, the turbulence interacts with the background shear flow, which increases the amplitude and anisotropy of the turbulence.The advantage of this model is that it allows to explicitly measure the contribution of different factors involved to the cycle of thermal and turbulent energy (Fig. 1).At the same time, it should be noted that this approach does not allow the study of the detailed spatial and temporal structure of turbulence, but only its effect on the mean flow.
The aforementioned model was used to study turbulence driven by convection in accretion discs.Two models of protoplanetary discs have been examined, one for the quiescent state and one for the outburst state.We can agree with the results presented in Held & Latter (2018) that convection-induced turbulence results in the outward transfer of angular momentum.The amplitude of the turbulence (Mach number ∼ 0.1) is in agreement with the estimates from molecular line observations (Flaherty et al. 2017(Flaherty et al. , 2018)).However, the turbulence is found to be too weak to either reproduce the heat flux through the dissipation channel (Q turb /Qext = 3.7 × 10 −2 in Model A) to support self-sustaining convection, or to generate sufficient torque to drive the accretion rate and the associated heating source ( Ṁeff / Ṁ = 1.7 × 10 −3 in Model A).A possible explanation for this is that convection only excites the vertical component of the Reynolds stress tensor, wzz, directly, while angular momentum transfer is controlled by the mixed component w rϕ .The latter is weak because of weak coupling of the turbulence and background velocity shear.It would be worth investigating instabilities occurring in both the radial and azimuthal directions as a potential source of mixed excitation.For example, subcritical baroclinic instability, magneto-rotational instability and others.

Fig. 1 Figure 1 .
Figure 1.Energy circulation in an accretion disc with a convective turbulence source.

Figure 3 .Figure 4 .
Figure 3. Qualitative indicators of the solutions to system (61)-(64), for q = 1.5.Top row: Minimum values of the components of the Reynolds tensor diagonal over the whole simulation time.Bottom row: Minimum values of the off-diagonal component wrϕ over the simulation time.Each dot is a simulation for some set of the constant parameters.The pictures in the bottom row show only the points corresponding to solutions with a non-negative diagonal, i.e. physically allowed solutions.The set of blue dots is obtained by uniformly sampling all three parameters: 0≤ C Π1 ≤ 5, 0 ≤ C Π2 ≤ 2, 0 ≤ C Π3 ≤ 2.The orange dots are obtained by sampling only the constant parameter that labels the horizontal axis, while the other two parameters are fixed at their experimental values (17).The latter are also marked with the thick dots.

Figure 5 .
Figure 5. Local turbulent convection model with q = 1.5, where different constant parameters are tested.Black lines: The experimental values of the constants, see Eq. (17).Coloured lines: One of the constants is varied (see the legend in the first plot).

Figure 6 .
Figure6.The structure of a protoplanetary disc with no convective energy transport and no turbulence.From left to right: gas number density, temperature, temperature gradient excess, and the magnitude of the convective flux.

Figure 7 .
Figure 7.The Model A results, Ṁ = 10 −7 M ⊙ yr −1 .Top row, from left to right: density, temperature, bulk radiation density.Middle row: Rayleigh number, energy flux, volumetric heat source (yellow dashed line denotes negative ρSconv values).Bottom row: turbulent stress tensor components, turbulent Mach number and volumetric force density.The blue dashed lines indicate the results of the model without taking into account convection.Dotted orange vertical lines mark the top boundary of the convective zone.Dotted black vertical lines mark the thermal scale of the disc.

Table 1 .
Parameters of disc models with turbulent convection

Table 2 .
Main results of Models A and B. Ṁ is the accretion rate parameterising the external heat source Qext; Ṁeff is the effective accretion rate; Q UV is stellar UV radiation source; Qconv is the source associated with the convection; Q turb is the turbulent energy dissipation rate; W ij are the z-integrated components of the Reynolds stress tensor; α eff is the effective Shakura-Sunyaev parameter.