Characterisation of global flow and local fluctuations in 3D SPH simulations of protoplanetary discs
Abstract
A complete and detailed knowledge of the structure of the gaseous component in protoplanetary discs is essential to the study of dust evolution during the early phases of preplanetesimal formation. The aim of this paper is to determine if threedimensional accretion discs simulated by the Smoothed Particle Hydrodynamics (SPH) method can reproduce the observational data now available and the expected turbulent nature of protoplanetary discs. The investigation is carried out by setting up a suite of diagnostic tools specifically designed to characterise both the global flow and the fluctuations of the gaseous disc. The main result concerns the role of the artificial viscosity implementation in the SPH method: in addition to the already known ability of SPH artificial viscosity to mimic a physicallike viscosity under specific conditions, we show how the same artificial viscosity prescription behaves like an implicit turbulence model. In fact, we identify a threshold for the parameters in the standard artificial viscosity above which SPH disc models present a cascade in the power spectrum of velocity fluctuations, turbulent diffusion and a mass accretion rate of the same order of magnitude as measured in observations. Furthermore, the turbulence properties observed locally in SPH disc models are accompanied by meridional circulation in the global flow of the gas, proving that the two mechanisms can coexist.
keywords:
accretion, accretion discs – hydrodynamics – turbulence – methods: numerical – protoplanetary discs.1 Introduction
Protoplanetary discs (abbreviated PPD) are discs composed mainly of gas and dust in quasiKeplerian rotation around young stars. They are believed to be the birth place of planets. The theory of planet formation is complex because several scales (from m to hundreds of thousands of km), several forces (e.g. electrostatic, magnetic, gravitational, drag) and several processes (e.g. turbulence, chemical and thermodynamical transformations, instabilities) are involved. One of the more poorly understood stages is the evolution of msize dust grains into kmsize objects, called planetesimals.
Among the several mechanisms the dust is subject to (e.g. radial drift and vertical settling), possible large scale motion and turbulence are expected to be of relevant importance. Large scale motion can lead dust to travel to very different locations in the disc. Turbulence can have two competing effects mediated by gas drag: stirring up and diffusing dust particles, impeding their agglomeration, or trapping them inside eddies, favoring their agglomeration (see e.g. Cuzzi et al., 1993; Carballido et al., 2008; Cuzzi et al., 2008).
Here we focus on the global flow and on the expected turbulent behaviour of only the gaseous component of PPD, which represents the medium in which dust evolves.
In the literature, two large classes of models of turbulent discs are widely used: continuous and discrete models. Continuous models of viscous discs are directly derived from the NavierStokes equations in presence of the gravitational potential of the central star, therefore they include the ingredients of physical viscosity and star gravity. Among the most popular models of this type we find the analytic 1D models by Pringle (1981) and LyndenBell & Pringle (1974). Since the physical molecular viscosity of the gas in PPD is very low (see e.g. Armitage, 2007), high Reynolds numbers are expected and therefore the gas is believed to be turbulent. The most commonly used model for such discs is the Shakura & Sunyaev (1973) model (hereafter SS73), which is a Prandtl model for turbulence: the NavierStokes equations are averaged and the Reynolds stress tensor is modeled by a viscosity called turbulent viscosity. The resulting total viscosity is the sum of the molecular and the turbulent viscosities, the former being negligible in PPD. Continuous models are usually applied to turbulent discs whose source of turbulence is not known or not explicitly declared; in such cases the SS73 relation is used to express the turbulent viscosity parameter. The dimensionless parameter , usually taken as a constant, collects all the ignorance of the source of turbulence, is the gas sound speed and the scale height of the disc. Generally, discs are good models for turbulence and capture its main effects. An example is given by the twodimensional models calculated by Takeuchi & Lin (2002, hereafter TL02) that show how the prescription is able to describe the mechanism of meridional circulation in discs. However, turbulent fluctuations are not directly reproduced in such disc models, only their effect on the dynamics of the global discs can be studied. In addition, if necessary, turbulent diffusion has to be added to the basic equations.
A different line of research starts from a proposed explicit hypothesis concerning the source of turbulence and derive the corresponding disc evolution. Given the complex dynamics, numerical simulations are often required, therefore these consist of mainly discrete models. Turbulent fluctuations can be resolved down to the resolution scale, which limits the size of the smallest structures. A first example of such models are discs produced by Magnetohydrodynamics (MHD) simulations where the source of turbulence is the MagnetoRotational Instability (MRI) (e.g. Fromang & Nelson, 2006). The equations at the base of such models are the Euler equations in presence of an adequate magnetic field and of the gravitational potential of the central star. Therefore they include the ingredients of magnetic field and star gravity, however physical viscosity (included in continuous models) is not present. In addition, turbulent diffusion (not added a priori) has been found in the outcome of simulations. A second example of discrete models are systems produced by simulations of selfgravitating discs where the source of turbulence is the Gravitational Instability (GI) (e.g Rice et al., 2005). The equations at the base of such models are the Euler equations in presence of both selfgravity and the gravitational potential of the central star. In addition, the numerical scheme is completed by an artificial viscosity term. Therefore they include the ingredients of selfgravity, star gravity and a ‘physicallike’ viscosity.
The nature of the source of turbulence in accretion discs and particularly in PPD is still an open issue. Even in the absence of magnetic fields and selfgravity, other mechanisms can be sources of turbulence, like interactions with external objects or pure hydrodynamic and thermal instabilities, such as the Rossby instability (Lovelace et al., 1999) or the baroclinic instability (Klahr & Bodenheimer, 2003). The possibility to have a hydrodynamic instability in accretion discs is still debated. In fact, due to the Rayleigh criterion (see e.g. Balbus & Hawley, 1991; Armitage, 2007), Keplerian discs are believed to be stable with respect to hydrodynamic instabilities. However this is proved only for 1D discs, if the third dimension is considered the behaviour could change and there may be the possibility to trigger instabilities (see e.g. Nelson et al., 2012).
As done in most continuous analytic models, we prefer to model the effects of turbulence (and not its cause) by the use of an effective viscosity. Here we consider threedimensional accretion discs modeled by means of standard Smoothed Particle Hydrodynamics (SPH, for a review see Monaghan, 2005), which we call SPH disc models, they are discrete models, such as those of the second class. The SPH scheme includes an artificial viscosity term (a Von Neumannlike viscosity) that has been introduced to treat shocks correctly (see e.g. Cullen & Dehnen, 2010). Far from the shock, the artificial viscosity is usually turned off by a switch. Different flavours of the artificial viscosity term are discussed in Sect. 2.1.
The aim of the present paper is to answer the following question: can SPH disc models correctly reproduce both the observed properties of PPD and the expected effect of turbulence? In order to clarify these two points we present a detailed characterisation of the properties of the global flow and particularly of the density and velocity fluctuations of SPH disc models.
The global flow of SPH models gives us information about the average properties of the gaseous disc, while the divergence from the average values defines the SPH fluctuations associated to the considered field (e.g. velocity and density fields). In particular, SPH fluctuations result from the combination of standard numerical noise and physical fluctuations. The former, present in all numerical schemes, is due in the case of the SPH method both to the discretisation of the continuum equations and to the SPH approximations; it depends on the number of particles and on the SPH kernel used in the simulations (see Sect. 2.1). The latter is present if the gas flow is turbulent; it depends on the Reynolds number of the flow which in turn is related to the physical viscosity of the gas. Both the numerical noise and physical fluctuations are related to the artificial viscosity term in the specific way we will show in this paper.
Turbulence modeling with the SPH method has recently become an active field of research (see Monaghan, 2005, and references therein). The first effort of implementing turbulence models in the SPH equations (e.g. Violeau & Issa, 2007) is now moving in the direction of quantifying the ability of the SPH method to intrinsically reproduce turbulence (e.g. Ellero et al., 2010; Monaghan, 2011). A comparison between SPH and grid methods is presented in Price & Federrath (2010) and a detailed and clear relation between turbulence and resolution in the SPH method is explained in Price (2012). In those studies, isotropic homogeneous turbulence in 2D boxes with periodic boundary conditions and with negligible gravity effects are considered. However, the case of anisotropic and inhomogeneous turbulence with free boundary condition such as in 3D gaseous discs, where also the action of gravity from the central star is relevant, has been addressed only in a few works (Murray, 1996; Lodato & Price, 2010), but without considering the global flow and the statistics of the velocity and density field and the match with the available observations for PPD. The present work aims to contribute in this direction.
In Sect. 2 we describe the main features of the SPH code we have used to model the accretion discs considered in the present work, the adopted reference disc model and the simulations of the different models we performed. In Sect. 3 we present the diagnostics applied to characterise both the global gas flow and the fluctuations of the velocity field of the gas, along with the known reference values for PPD. We then show the results of the applied diagnostics respectively for the global flow (Sect. 4), the magnitude of SPH fluctuations (Sect. 5) and the structure of SPH fluctuations (Sect. 6). The results are discussed in Sect. 7 where we also draw our conclusions.
2 SPH models of accretion discs
In this work we focus on a PPD described by given analytic initial density and locally isothermal sound speed profiles. The system is not homogeneous nor isotropic. We calculate the quasistationary state of such a disc using different values of the three numerical parameters (number of particles) and and (artificial viscosity parameters) controlling the SPH fluctuations, in order to characterise their effects.
In this section we present the code and the adopted units, the reference disc model and the sets of simulations we performed.
2.1 The SPH code
We use the twophase SPH code described in BarrièreFouchet et al. (2005). The two phases represent gas and dust that interact via aerodynamic drag. The gas is described by the Euler equations and artificial viscosity.
Here we consider only the gas phase because we are interested in the characterisation of the global gas flow and its velocity and density fluctuations. Thus, the momentum conservation equation takes the form:
(1) 
where the usual SPH approximation of replacing integrals with sums over a finite number of particles has been performed. The term is the velocity of SPH particle , the time, the SPH particle mass, the pressure, the mass density, the position of particle with respect to the central star of mass , a parameter used to prevent singularities and the SPH smoothing kernel, with the distance between the particle and its neighbour , and the smoothing length. Here we use a cubicspline kernel truncated at :
(2) 
where , is the dimension of the simulation (in this work = ) and = is the normalization constant respectively in 1, 2 and 3 dimensions. The effect of different kernels will be addressed in a future work.
The smoothing length is variable and is derived from the density :
(3) 
with . This choice of guarantees a roughly constant number of neighbours (for 3D simulations: ). The smoothing length defines the resolution of the SPH simulations: for a larger number of particles the target number of neighbours is reached inside a smaller volume (of radius from Eq. 2), implying a smaller smoothing length. The gas is described by a locally isothermal equation of state. The artificial viscosity term is described in Sect. 2.1.1.
2.1.1 The artificial viscosity term
Two different implementations of the artificial viscosity term are considered. The first one, originally introduced by Monaghan & Gingold (1983) (hereafter MG83) and subsequently refined by Lattanzio et al. (1985); Monaghan & Lattanzio (1985); Monaghan (1992), is defined by:
(4) 
where and are respectively the relative distance and relative velocity of particles and . The overlined quantities are averages between particle and its neighbouring particle : , , . Finally,
(5) 
with . Different combinations of the two artificial viscosity parameters have been used so far in different applications, in particular the combination has been claimed to give good results (Monaghan, 1992, 2005). We will refer to and as ‘standard’ values for the artificial viscosity parameters.
The second implementation of artificial viscosity (hereafter identified as Mu96) is that adopted by Murray (1996) and Lodato & Price (2010):
(6) 
In contrast to the MG83 implementation, here the artificial viscosity is applied to all particles (the switch present in the MG83 version is not applied) and the term is set to zero. Therefore, discs with the standard MG83 artificial viscosity implementation are naturally half as viscous as those with the Mu96 implementation.
2.1.2 The link to physical viscosity
The reason of the Mu96 choice is to better match the conditions under which the artificial viscosity can represent a physical viscosity. In fact, Meglicki et al. (1993) showed that in the 3D continuum limit the Mu96 artificial viscosity has the form of a physical (shear and bulk) viscous force. The shear term is equivalent to an viscosity (that we call , where the subscript ‘cont’ refers to continuum) given by:
(7) 
where the average on is taken along both the azimuthal and vertical directions and the semithickness of the disc is defined in Sect. 2.2 (note that this expression is valid for the cubic spline kernel). However, depends on the resolution, because of the presence of the smoothing length (see Lodato & Price, 2010).
The SPH equations with the Mu96 artificial viscosity and with a very large number of particles (high resolution) are therefore equivalent to those of a viscous fluid with an given by Eq. 7. Thus, artificial viscosity can be used to control the effective viscosity in the disc and it is expected to be responsible of physical fluctuations of the flow when the gas is in the turbulent regime.
The 3D continuum limit of the MG83 artificial viscosity is much more complex than the one presented for the Mu96 case and at the moment a simple relation such as that in Eq. 7 is only available for the term (Meru & Bate, 2012):
(8) 
The factor of two between the numerical coefficients of Eqs. 7 and 8 naturally arises from the fact that the MG83 artificial viscosity is applied only to half the particles because of the presence of the switch (see Eq. 4). Since the Mu96 formula (developed much earlier than the MG83 formula) is usually adopted to estimate the effective viscosity in SPH discs independently of the particular AV implementation, SPH discs are generally less viscous than what has been considered so far. In Sect. 4.4 we present a method to determine the effective viscosity of an SPH disc with a generic artificial viscosity implementation. It is tested using Eq. 7 in the Mu96 case. It can therefore be used for sampling numerically the  relation in the full () MG83 case, through highresolution simulations.
We note that in the majority of SPH simulations the use of the artificial viscosity is still preferred to the direct implementation of the viscous stress tensor to model a pure shear viscosity. The main reason is that the latter method does not conserve the total angular momentum exactly (Schäfer, 2005; Lodato & Price, 2010) as the artificial viscosity term does. Secondarily, it is more computationally demanding because of the presence of the second spatial derivative of the velocities.
2.1.3 The units
The internal units of the code are chosen to be 1 M for mass, 100 au (Astronomical Units) for length and to give a gravitational constant . With these values the time unit is therefore yr.
2.2 The reference disc model
We consider a typical T Tauri disc of mass orbiting around a one solar mass star ( M). It extends from to au and is characterised by a density profile given, in cylindrical coordinates, by
(9) 
(expansion at small , see Laibe et al. 2012 for the rigorous expression) with , and .
(10) 
is the semithickness of the disc, related to the sound speed and the angular velocity . The disc is locally isothermal with a temperature radial profile given by
(11) 
with , which leads to the sound speed profile
(12) 
Note that the sound speed coefficient and the sound speed exponent determine respectively the semithickness of the disc and its radial dependence:
(13) 
with . The resulting radial profile of the surface density is
(14) 
with the surface density at and .
When only the gas phase is considered, all models are selfsimilar and different physical scales correspond to the same dimensionless model. Therefore, in the following, results are mainly expressed in code units.
The reference values we adopt in the following are and au. At this location, the disc is slightly flared with , kg m and m s ( K). For reference, the corresponding value at 1 au are: kg m and K.
We follow the evolution of the reference disc model, after pressure equilibrium has been reached, up to time in code units (corresponding to 15.9 orbits at 100 au). All figures are plotted for that time, unless stated otherwise.
2.3 The simulations
The simulations we have performed are listed in Table 1. The first column gives the name of the simulation, the second and third columns display the values of the two artificial viscosity (AV) parameters and , the fourth the number of SPH particles used for sampling the disc, the fifth the kind of artificial viscosity. Simulations can be divided in five sets. Each simulation can belong to more than one set. Sets are displayed in the last column of the table. Simulations in Set and Set allow to study the effect of changing the artificial viscosity parameters for the MG83 and for the Mu96 artificial viscosity respectively. With simulations in Set it is possible to study the effect of the different artificial viscosity implementations. In Set are collected simulations designed for studying the effect of resolution. Finally, simulations in Set are used for testing the  relations presented in Eqs. 7 and 8. The role of each set is summarised in Table 2.
Name  AV  Set  
S1  0.1  0.0  MG83  ,  
S2  0.1  0.2  MG83  
S3  0.1  0.5  MG83  
S4  0.1  2.0  MG83  
S5  0.1  10.0  MG83  
S6  1.0  0.0  MG83  ,  
S7  1.0  2.0  MG83  ,  
S8  1.0  5.0  MG83  ,  
S9  1.0  10.0  MG83  
S10  2.0  0.0  MG83  
S11  2.0  4.0  MG83  
S12  2.0  10.0  MG83  
S13  5.0  0.0  MG83  
S14  5.0  2.0  MG83  
S15  5.0  10.0  MG83  
S16  0.1  0.0  Mu96  ,  
S17  1.0  0.0  Mu96  ,,  
S18  2.0  0.0  Mu96  
S19  1.0  2.0  MG83  
S20  1.0  2.0  MG83  ,  
S21  1.0  5.0  MG83  
S22  1.0  5.0  MG83  
S23  1.0  5.0  MG83  
S24  1.0  5.0  MG83  
S25  1.0  0.0  Mu96  ,,  
S26  5.0  0.0  Mu96  
S27  5.0  0.0  Mu96 
Set  Description 

Changing and in MG83 artificial viscosity  
Changing in Mu96 artificial viscosity  
Changing the artificial viscosity model (MG83, Mu96)  
Changing  
Testing the  relations (Eqs. 7 and 8) 
Midplane and vertical cross sections of the density profiles of four of the performed simulations are shown in Fig. 1.
3 Diagnostics and reference values for PPD
We consider two kinds of diagnostics which refer to the global flow and to fluctuations of several quantities, with particular focus on density and velocity field. We first define the diagnostic and then give the expected values in some known cases. In the following, the symbol represents average quantities. In this paper we perform both the standard azimuthal and vertical averages and the azimuthal averages only , in order to characterise also the vertical extension.
3.1 Diagnostic for the global disc flow
Here, we present the selected quantities used to characterise the global flow.

Mass distribution. We look at:

the radial profile of the surface density ,

the density distribution in the direction perpendicular to the disc midplane, for a selected radius .


Velocity structure. We focus on:

the radial velocity maps , where the radial component of the velocity field is azimuthally averaged.

the local Mach number of the flow in the midplane Ma, where is the modulus of the local velocity field.


Macroscopic turbulence signatures. Four items are analysed.

The mass accretion rate onto the central star : we are particularly interested in this quantity for two reasons: (1) it gives indirect information about the turbulent viscosity coefficient and (2) it is constrained by observations. For PPD around T Tauri stars the measured values of mass accretion rates onto the central star are M (Hartmann et al., 1998; Andrews et al., 2009), they can be reproduced by a value (Hartmann et al., 1998; King et al., 2007).

Since gives an information restricted to the very inner region of the disc (where our inner boundary condition is free), we also analyse the radial profile of the local mass accretion rate: , averaged in the azimuthal and vertical direction.

The effective viscosity of the gaseous disc (we express the viscosity in terms of the dimensionless parameter of the SS73 parametrization) that characterises the simulated disc.
The effective viscosity in the simulated disc model is estimated by means of fits to the twodimensional analytic models of accretion discs (see details in Appendix B). We use 2D analytic models because they are well suited for the azimuthal symmetry present in our simulated discs. For this reason we call this way of estimating the effective viscosity . In particular, we derive by comparing the vertical profile of the radial velocity of the gas flow to the expression derived by analytic disc models for different radial positions. We rewrite the vertical profile of the radial velocity, given by TL02 and Fromang et al. 2011, highlighting the scale height of the disc at radial position :
(15) Lodato & Price (2010) estimated by fit of the surface density profile of 1D viscous accretion discs (Pringle, 1981) and found good agreement with in SPH simulations of warped discs when the necessary equivalence conditions (large number of particles and Mu96 artificial viscosity, see Sect. 2.1.2) are met.

Once we know the effective viscosity of the disc, we can derive the corresponding Reynolds number:
(16)

3.2 Diagnostic for fluctuations
For each simulated disc we have studied both the magnitude and the structure of the fluctuations present in the disc with particular attention to the velocity field. Since current observations have not yet reached the resolution necessary to directly detect turbulence and study its features in protoplanetary discs (some progress has been recently claimed by Hughes et al. 2011 and Guilloteau et al. 2012, who measured some turbulent velocities in mm observations) the resulting properties of density and velocity fluctuations have been compared to the typical behaviour observed in turbulence experiments and to results from grid and particlebased simulations available in the literature.
The components of the velocity field are identified by and those of the fluctuating velocity field are with .
3.2.1 Magnitude of velocity fluctuations
Concerning the magnitude of velocity fluctuations, we focus on the turbulent viscosity coefficient (point M1), the diffusion coefficients (point M2) and the Mach number of velocity fluctuations (point M3).

The turbulent viscosity coefficient is often derived from Reynolds averages of NavierStokes equations with the turbulent viscosity hypothesis (e.g. Pope, 2000) and is thus related to Reynolds stresses:
(17) (see Appendix A.1).
Here we call the corresponding SS73 parameter for a disc in quasiKeplerian rotation, it is given by:
(18) where Eq. 17 has been combined with the SS73 relation = (see Appendix A.2). For each radial and vertical position, averages are performed spatially along the azimuthal direction in order to be able to derive the radial and vertical dependence of .
The corresponding turbulent Reynolds number is:
(19) 
Turbulent flows are often approximated by means of models where turbulence is described as a diffusive process. Here we want to determine if diffusion is present in our simulations. To this end the method used by Fromang & Papaloizou (2006) in order to calculate the turbulent diffusion coefficient is well suited for Lagrangian codes such as SPH:
(20) where is the vertical velocity correlation function:
(21) with representing the vertical position of a given particle at time and the vertical position of the same particle at time ; in a similar way is the velocity of the same particle at time and at time . In order to compute the diffusion coefficient at a selected radial location , we average the diffusion coefficient of single particles that were at at the beginning of the simulation.
A diffusion coefficient is found by Fromang & Papaloizou (2006) in MHD local shearing box stratified simulations.

The strength of velocity fluctuations and therefore their subsonic or supersonic behaviour can be determined by the corresponding Mach number, defined as the ratio of the modulus of velocity fluctuations to the gas sound speed:
(22) We have computed azimuthalonly and azimuthal and vertical averages in order to present – maps and radial profiles.
3.2.2 Structure of velocity fluctuations
Concerning the structure of velocity fluctuations, we focus on the anisotropy vector (point S1), on the power spectrum (point S2) and on the possible presence of intermittency (point S3).

In order to determine the isotropic or anisotropic character of velocity fluctuations, we define the anisotropy vector = in terms of the velocity dispersion in the three spatial directions, in analogy with galactic dynamics (see e.g. Bertin, 2000):
(23) where the square of the velocity dispersion = = = and , and refer to the three cylindrical components. Systems with anisotropy in the direction are characterised by a velocity dispersion of the component that is larger than that of the other two components: and . Therefore, the component of the anisotropy vector will be . Averages are performed along both the vertical and the azimuthal direction.
From the plots of Fromang & Nelson (2006) and Flock et al. (2011) we can deduce that , therefore MHD discs are radially anisotropic and the corresponding anisotropy vector is . In contrast, in experiments of turbulence produced by round jets a corresponding azimuthal anisotropy has been observed (e.g. Pope, 2000).

The power spectrum is computed along a ring of selected radius , centred on the star and located at the selected vertical position (see Appendix A.3):
(24) with the Fourier transform of the component of the velocity fluctuation vector and the wavenumber corresponding to length scale . We analyse the power spectrum of the velocity field in order to determine if an energy cascade, characteristic of turbulent systems, is present and if differences related to the position inside the disc are relevant.

Highly turbulent flows are characterised by intermittency (e.g. Frisch, 1996). This feature implies nonGaussian probability distribution functions (PDFs) with corresponding nonGaussian higher order moments. Here we focus on the density and acceleration PDFs and on their (skewness, noted ) and (kurtosis, noted ) order moments (see Appendix A.3).
Data are not available for the power spectrum of fluctuations and for the phenomenon of intermittency in accretion discs. In the classical Kolmogorov theory of incompressible homogeneous turbulence the power spectrum of velocity fluctuations has a power law form with = . Recent simulations of supersonic compressible gas found a slope of (Price & Federrath, 2010, and references therein). Intermittency, highlighted by nonGaussian PDFs and higher order moments, is experimentally found in very high Reynolds number incompressible flows. In SPH simulations of a simple weakly compressible periodic shear flow, Ellero et al. (2010) find that the PDF of the particle acceleration is in good agreement with nonGaussian statistics observed experimentally for incompressible flows.
In simulations of supersonic compressible turbulence (see e.g. Price & Federrath, 2010), lognormal distributions of the density PDF have been observed:
(25) 
with and a width controlled by the Mach number of the fluctuations
(26) 
where . The expressions of the skewness and kurtosis in terms of for the lognormal distribution are
(27) 
When the lognormal distribution tends to a Gaussian distribution and the skewness and kurtosis tend to their gaussian values (0 and 3 respectively).
4 The disc global flow
We now consider how the main properties of the global disc flow (mass distribution, velocity structure, mass accretion rate and effective disc viscosity) depend on the strength of artificial viscosity. The convergence of each quantity is also studied.
4.1 Mass distribution
The position of the peak of the surface density profile moves towards larger radial distances and its maximum decreases (see left panel of Fig. 2), when the two AV parameters and are increased: such behaviour is in agreement with the viscous spreading mechanism that is expected here due to the existing correlation between artificial and physical viscosity (see Sect. 2.1.2).
An increase in resolution (see right panel of Fig. 2) has a qualitatively similar effect to a decrease in viscosity: the surface density peak moves towards the centre of the disc, becoming higher (note also that the curve becomes less noisy). In the central and external regions the projected density profile converges already for a number of particles as small as . In the inner region (), a trend towards convergence is present, but a higher number of particles, , is necessary.
The vertical density profile at a selected radius (see insets in Fig. 2 for ) follows the expected Gaussian distribution of Eq. 9, showing that the desired profile is well reproduced by the sampling procedure and by numerical thermalisation. It is not significantly affected by changing , or .
4.2 Velocity structure
4.2.1 Radial velocity maps
Maps of the azimuthally averaged radial velocity in the  plane are displayed in Figs. 3 and 4 for changing artificial viscosity and resolution respectively.
Effect of artificial viscosity parameters. A transition from a chaotic to a regular accretion pattern is clearly observed in Fig. 3 when is increased for a given (panels from top to bottom) and as well as when is increased for a given (panels from left to right).
At large radii, a well extended decretion region dominates in all models, as expected from the conservation of angular momentum. At smaller radii, instead, the structure of the gas flow strongly depends on the values of the AV parameters. A decretion flow located around the midplane, associated with an accretion flow in thin surface layers, appears when are increased. Such accretion layers become thicker for larger AV parameters. Therefore, the gas flow is characterised by the phenomenon of meridional circulation. This mechanism has been found in some twodimensional viscous disc models (e.g. TL02 that extends the standard one dimensional SS73 models to two dimensions: radial and vertical) but is not reproduced by MHD disc simulations where turbulence is induced by the MRI (Fromang et al., 2011; Jacquet, 2013). The presence of meridional circulation has important implications for the topic of particle mixing in PPD and is one of the possible mechanisms able to explain the presence of crystalline solid particles in the outer regions of T Tauri stars recently observed by Spitzer (e.g. Sargent et al., 2009). Finally, in the very inner region of the disc a small and well defined accretion region is present in all models.
Effect of artificial viscosity implementation. The two panels in the bottom row of Fig. 4 show the structure of the radial velocity flow of two discs at different resolution simulated with the Mu96 artificial viscosity. The resulting flow structure can be directly compared with the panel above, which represents the flow structure in an equivalent disc at the same resolution but simulated with the MG83 artificial viscosity. In the simulations at lower resolution (), the meridional circulation pattern is better reproduced by the MG83 implementation. However, in the more resolved discs () the two patterns are very similar. The two small accretion regions present in the outer Mu96 discs are due to a longer relaxation timescale for this particular artificial viscosity implementation (in longer evolution simulations, not shown here, we have observed a net outward gas flow in agreement with the MG83 case and the expected viscous spreading).
Effect of resolution. The first row in Fig. 4 shows how the radial velocity structure has already converged after particles.
TL02 found that the outflow zone around the midplane shrinks as the radial volume density gradient is reduced () and when the condition , that guarantees a net accretion of the gas onto the star, is violated. Therefore, we expect that the thickness of the accretion layers depends on the values (, ) of the surface density and temperature radial profiles of the gaseous disc model. However, a detailed study of the dependence of the structure of the accretion layers on the gas density and temperature profiles and the correlation with the central mass accretion rate is out of the scope of this paper and will be addressed in a future work. For the initial profile of the disc model used in this work we have , which implies net decretion but is very close to the limit between net accretion and decretion. The surface density of the simulated discs after numerical relaxation reaches a smooth profile that, in the inner region, is slightly shallower (implying a slightly smaller ) than the initially imposed power law, which is characterised by an inner sharp edge (see Sect. 2.2). This moves the model in the net accretion region and explains the presence of the observed thin accretion layers.
4.2.2 Mach number of the flow
The gas in the disc is supersonic. In fact, since the gas moves on quasiKeplerian orbits, the Mach number is approximated by Ma . With the adopted sound speed profile it becomes: Ma . For the values used here Ma , which is approximately in the range of values 15–25.
4.3 Mass accretion rate
The mass accreted onto the central star (Fig. 5) and the radial profiles of the mass accretion rate (Figs. 6 and 7) are two macroscopic signatures of turbulence.
4.3.1 Central mass accretion rate
Effect of artificial viscosity. The left panel in Fig. 5 shows that increases with larger and/or larger and in the case of the Mu96 artificial viscosity implementation. The effect of is larger, as expected since it controls the first order in the artificial viscosity term.
Effect of resolution. The right panel in Fig. 5 shows that decreases with increasing resolution. Such an effect is due to a decrease of particle noise at higher resolution, which decreases the numerical dissipation. This is linked to the dependence of the effective viscosity on resolution (see Sect. 4.4). In fact, as will be shown in Sect. 6.3.1, the PDFs of several physical quantities are narrower at higher resolution, due to the reduced noise. Concerning the central mass accretion rate we look at the part of the radial velocity distribution corresponding to negative velocities, since only fluid elements with negative radial velocity are responsible for accretion. The radial velocity distribution at higher resolution presents a smaller spread than in the lower resolution case. In addition, SPH particles in low resolution simulations have a larger mass than in higher resolution simulations. Therefore, at lower resolution there is a larger number of more massive SPH particles with higher inward radial velocity, implying higher accretion rates. The trend of the curve shows that the central accretion rate is converging, however convergence is not yet reached at one million particles.
For all the considered combinations of (,) and , we find mass accretion rates consistent with the values observed for PPD. In addition, we observe a correlation between the trend of with respect to , and the increase of thickness of the accreting layers (net accreting mass flux) observed in Fig. 3. Similarly, the trend with respect to correlates with the decrease of thickness of the accreting layers at higher resolution, observed in Fig. 4.
4.3.2 Radial profile of the mass accretion rate
The central mass accretion rate gives information only about the region of the disc that is close to the star. In order to have a global picture of the mass flux in the disc, it is interesting to look at the radial profile of the vertically integrated mass accretion rate . As an example, the profile observed in simulation S8 is displayed in the bottom panel of Fig. 6. The flow is characterised by an inner accreting region separated, at a radius we call , from an outer decreting region with small rates at intermediate radii and larger rates at large radii. It can be understood by a comparison to the surface density profile (top panel of Fig. 6). In the organized flow of model S8, corresponds to the radial location of the surface density peak, which moves outwards with time due to viscous spreading, as the density maximum decreases. The density peak behaves as a reservoir of gas that supplies mass both to the inwards and outwards flows, resulting in an apparent null mass accretion rate at . The varying rate in the decretion region is due to the vertical structure of the flow: Fig. 3 shows accretion throughout the disc height for the inner disc, accretion in the midplane with decretion in the surface layers at intermediate radii resulting in a slightly negative radial mass flux, and decretion throughout the disc height in the very outer regions.
Effect of artificial viscosity. For low values of and fluctuations are so large that the identification of an organized flow of the fluid is not possible. As shown by the top panel of Fig. 7, the mass accretion rate profile is very noisy with a narrow peak in the inner disc, an average value close to zero in the intermediate regions and negative in the outer disc. The peak of the surface density profile decreases with time but less than in the higher case and it does not change position with time. In the map (see the top left corner of Fig. 3) negative and positive regions are clearly mixed.
Increasing the two AV parameters, the profile becomes more regular and smooth, following the structure described above for the specific simulation S8. For = the profile tends to converge when is increased (see the middle panel). In the Mu96 case we observe similar radial profiles for the mass accretion rate (not shown).
Effect of resolution. As shown in the bottom panel of Fig. 7 for the simulations with and increasing resolution, the peak of the radial profile of the mass accretion rate moves inwards and converges for , in agreement with the convergence observed for the radial velocity maps. The reason of the decrease of the peak in higher resolution simulations is the same that explains the similar decreases observed for the central mass accretion rate (see Sect. 4.3.1).
4.4 The effective viscosity of the SPH disc:
In order to estimate the effective viscosity of the gaseous disc, we determine the parameter, which is derived by fitting the vertical profile of the radial velocity expected for a viscous axisymmetric accretion disc, given by Eq. 15, to that derived from the simulations (see Appendix B for the fitting procedure). Two examples, characterised by adequately large AV parameters, are shown in Fig. 8, where the vertical profile of the radial velocity at radial position is displayed for simulation S7 (top panel) and for simulation S20 (bottom panel).
In Fig. 9 we present the values of obtained by averaging the result of the fit at =1 and =2 for the set of simulated discs. The values are averaged in the radial range . In simulations with the vertical profile of radial velocity cannot be described by the profile in Eq. 15 because the high noise level makes the fit impossible. We find that simulation profiles can be fitted by a value of which is roughly the same for the different radial locations we considered. This shows that the structure of SPH models is consistent with the 2D version of SS73 accretion disc models, characterised by a constant value.
Effect of artificial viscosity. In the range of AV parameters were the determinations of the effective viscosity is possible (), the scaling of with the AV parameters (increase with larger and/or in the MG83 implementation and increase with larger for the Mu96 implementation) has the same qualitative behaviour as for the central mass accretion rate (see Fig. 5). Such a match is what we expect, since a more viscous disc is characterised by a larger central accretion rate. The weaker dependence on is in agreement with the fact that is the parameter controlling the second order term of the MG83 artificial viscosity (see Eq. 4).
Effect of resolution. The trend of the curves in the right panel of Fig. 9 is qualitatively similar to that of the central mass accretion rate (see Fig. 5) and shows that the effective viscosity is converging. However one million particles is still not enough in order to reach convergence. The slightly higher resolution required for the convergence of and with respect to that required by the radial velocity maps and by the mass accretion rate radial profile () is due to the fact that the former two parameters are calculated from a smaller subset of simulation particles, the numerical noise is therefore larger.
The result of this subsection is that the effective viscosity in all the SPH models considered here equals a few , in agreement with the order of magnitude deduced from observations.
4.4.1 Reynolds number
Once the effective viscosity of the disc is determined (), it becomes possible to derive the corresponding Reynolds number , shown in Fig. 10. In accordance with the trend of with respect to , and , described in the last paragraph, we found larger Re numbers in less viscous and more resolved discs.
Simulations with one million particles and the MG83 artificial viscosity are closer to, but sill lower than the value , which is considered as a limit for the onset of the turbulent regime in a box with periodic boundaries (see Price, 2012). Therefore, the dynamics of the discs presented in this paper is not expected to be dominated by fullydeveloped turbulence.
4.4.2 The  connection
Figure 11 shows the relation between the effective viscosity (computed in the simulations) and the viscosity (expected when the continuum limit of the SPH equations is taken) for eight of the performed simulations with (S6, S10, S13, S17, S18, S25, S27). Simulations S1 and S16 (with ) are too noisy to allow the determination of . Since the ratio in the expression of depends on the radial location inside the disc (see Eqs. 7 and 8), we have evaluated locally at the radial position both and , whose values are given in Table 3.
We find that both the points concerning the Mu96 case and those concerning the MG83 case follow the appropriate analytic relation. We therefore confirm the necessity of reducing the coefficient in the relation in Eq. 8 for the MG83 case as shown by Meru & Bate (2012) and motivated by the fact that in the MG83 implementation the artificial viscosity is applied only to approaching particles.
Finally, no significant difference is observed when resolution is changed from to particles, suggesting that already with particles the continuum approximation can be considered valid.
Simulation  

S6  0.0301  0.0490 0.0063 
S10  0.0597  0.0668 0.0040 
S13  0.1476  0.1453 0.0023 
S17  0.0601  0.0759 0.0043 
S18  0.1204  0.1249 0.0033 
S27  0.3009  0.2949 0.0075 
S25  0.0372  0.0500 0.0013 
S26  0.1852  0.1776 0.0009 
5 The magnitude of SPH fluctuations
5.1 The Reynolds stress contribution to the SPH disc viscosity:
We now proceed with the determination of , which is the parameter derived from the Reynolds stress, in analogy with turbulence studies and with MHD or selfgravitating disc simulations. Radial and azimuthal fluctuations are computed from the average velocity field determined locally using a – grid.
The radial profile of the vertically averaged is shown in the left panel of Fig. 12, where the effect of changing and is highlighted, and of Fig. 13, where the effect of resolution is shown. In all simulations, the coefficient is characterised by a similar profile: it tends to be approximately constant throughout the radial extension of the disc, with an increase both in the central and in the external radial region.
The right panel of both figures shows the radially averaged vertical profile (as in Fig. 8, the semithickness of the disc is computed from =, see Sect. 2.2): is approximately constant around the midplane and then increases with height. This behaviour qualitatively agrees with that observed in MHD discs by Flock et al. (2011), even if we observe a broader minimum. In addition, we have found that the vertical profile is approximately the same at different radial positions, as shown in Fig. 14 for simulation S24.
Effect of artificial viscosity. The radial and vertical profiles of are very sensitive to the AV parameters and . In fact there is a sharp change in the average value, which decreases from positive to negative values, when the two AV parameters are increased. This transition is clearly visible in both panels of Fig. 12 for : the increase of from 0 to 2 shifts the profile of towards smaller values, moving it in the negative region. A further increase of leads to smaller changes, but the behaviour inverts: higher artificial viscosity now tends to increase the value of the coefficient (which corresponds to a decrease of its absolute value, see left panel) and to extend the constant region around the midplane (see right panel). We find the same trend for other values of (not shown).
The change of sign of is due to the mechanism of meridional circulation that is correctly resolved when adequately high AV parameters are used.
Effect of resolution. For a given combination of AV parameters, the radial and vertical profiles of the coefficient do not change with resolution. The only effect of increasing the number of particles of the simulation is a reduction of the noise, which implies smoother profiles.
We note that the computation of is very sensitive to the correct determination of the average velocity field of the flow, therefore only in simulations with adequately high values of and , where numerical noise and particle disorder are lower, can it be trusted. The absolute value of for the present simulations is . This is comparable to standard values found in MHD discs (see e.g. Fromang & Nelson, 2006). The difference is that meridional circulation is present here and the negative value indicates decretion on average in the midplane, in contrast to the positive value of MHD discs (where meridional circulation is not seen), which is expected for accretion in the midplane.
In conclusion, more viscous discs are characterised by a smaller modulus of and, for all the simulations we performed, we find the ratio of the Reynolds stress to the effective viscosity of the gaseous disc to be . This means that velocity fluctuations are not the main component of the effective viscosity of the disc, in agreement with the low Reynolds number of the global flow in in the present simulations (see Sect. 4.4.1). However, the turbulent Reynolds number, which is associated to the Reynolds stress, amounts to (since Ma and ), close to the turbulent limit. Therefore, physical fluctuations, quantified by the small contribution of to the effective viscosity, can be expected to exhibit turbulent signatures. We will show in Sects. 5.2 and 6.2 that it is indeed the case.
An important point to be highlighted is that the determination of is not a sufficient condition in order to determine the mass accretion rate of the disc, since other sources can be present in the model.
5.2 Diffusion coefficient
For each simulation, we have computed the diffusion coefficient from Eq. 20 at several radial positions in the disc. In Fig. 15 we show the evolution with time of the diffusion coefficient , in units of , measured at the intermediate radial position for some simulations belonging to Set A. Except for simulations from S1 to S3, characterised by larger fluctuations, the diffusion coefficient regularly increases with time and converges to a welldefined value.
For the standard combinations of parameters , is in agreement with the value found by Fromang & Papaloizou (2006) in their MHD local shearing box stratified simulations.
The asymptotic values of the diffusion coefficient for each performed simulation are displayed in Fig. 16.
Effect of artificial viscosity. The left panel of Fig. 16 shows that the coefficient decreases with increasing and (except simulations from S1 to S3 which present a diffusion coefficient without a clear trend, due to the higher noise level) and in the Mu96 case (the artificial viscosity switch on approaching/receding particles is off, resulting in more viscous discs). This trend is due to a reduction of the noise that allows the onset of meridional circulation, which forces the fluid to be more localised in the vertical direction.
Effect of resolution. From the right panel of Fig.16 we observe that the diffusion coefficient decreases with increasing resolution. The reason of this behaviour is the same that explains the qualitatively similar trend we have observed in Sect. 4.3.1 for the central mass accretion rate. In this case, the PDF of the vertical velocity is broader in low resolution simulations with respect to more resolved simulations. Therefore, larger vertical velocity are possible in less resolved discs with resulting higher values (see Eq. 21) of the diffusion coefficient. The curve shows a converging trend, however convergence is not yet reached at one million particles. As for the central mass accretion rate and for the effective viscosity fit, convergence requires more resolution because the diffusion coefficient is computed from a smaller subset of particles than in the case of the radial velocity maps and the mass accretion rate profiles, where an average on the azimuthal component is taken.
We conclude that the diffusive mechanism is correctly represented by intermediate/large AV parameters : more viscous discs are less diffusive, due to the onset of meridional circulation and the decrease of . In addition, for a given combination, less resolved discs are more diffusive: such numerical dependence is of the same type as the dependence of the central mass accretion rate and of the effective viscosity on the resolution of the simulation (see Sect. 4.3.1.)
5.3 Mach number of velocity fluctuations
The azimuthally averaged Mach number of the modulus of velocity fluctuations is displayed in the – plane in the left panel of Fig. 17 for simulation S23, as an example. The right panel shows the corresponding radial profile for each velocity component, vertically and azimuthally averaged. All simulations show qualitatively similar structures and profiles.
Effect of artificial viscosity. The average value of in the radial range is displayed in Fig. 18. From the left panel of the figure, we see that, except for simulation S1, whose fluctuations are close to the sonic regime, all simulations show subsonic fluctuations with Mach number significantly smaller than unity. We observe a decrease of the average Mach number with larger AV parameters. For relatively high and , at each given radius the Mach number increases from a minimum value in the midplane () to a maximum value in the surface layers (), as displayed by the left panel of Fig. 17. This behaviour is qualitatively in agreement with what is found in MHD simulations (e.g. Fromang & Nelson, 2006; Flock et al., 2011), however we observe smaller values. A global maximum is reached in the inner region, where and the two accreting layers meet. In addition, two layers with locally higher fluctuations with are also present at the interface between the surface accreting layers and the midplane outwards flow.
Effect of resolution. Similarly, as a consequence of noise reduction, higher resolution leads to a decrease both of the Mach number and of the extension of the surface layers where it reaches the local maximum for a given radial position (with the same qualitative behaviour of the accreting layers in the maps, see Figs. 3 and 4), as expected. We have found a convergence of the vertically and azimuthally averaged Mach number for a number of particles (see right panel of Fig.18), in agreement with the value found for radial velocity maps (which are connected to radial velocity fluctuations and therefore to the radial part of the Mach number).
Fluctuations tends to be dominant in the regions of the disc where accreting and decreting flows meet: the conelike regions between the accreting surface layers and the decreting midplane layer, and the region around the inner disc edge.
6 The structure of SPH fluctuations
6.1 Anisotropy of velocity fluctuations
We recall that when the system is totally radially anisotropic, while for negative values the system is dominated by anisotropy in different directions. The same holds for the azimuthal and vertical components. The kind of anisotropy of the disc is displayed in Fig. 19 for simulations with changing and . We observe that none of the simulated discs is totally anisotropic in one particular direction, but there is always a component that dominates the other two.
Effect of artificial viscosity. Discs with smaller AV parameters are dominated by radial anisotropy. For example, in the case of simulation S1 with (,)=(0.1,0), we see from Fig. 19 that , and , it follows that is significantly larger both than and than . When and/or are increased, decreases and becomes negative while increases. The vertical anisotropy is always negligible with respect to the radial or azimuthal anisotropy (it is always negative in all performed simulations, as shown by the right panel in Fig. 19).
Effect of resolution More resolved simulations are slightly more azimuthally anisotropic. As for some of the previously considered quantities, we observe that the components of the anisotropy vector, averaged in the vertical, azimuthal and in the usual radial range , converge at (not shown), for the same reason already outlined several times.
6.2 Power spectrum
For each simulation we have computed the power spectrum of velocity fluctuations in a ring located at different radii from the central star and at several altitudes, as explained in Sect. 3.2. As an example, in Figs. 20 and 21 we present compensated spectra (divided by the Kolmogorov slope of , for easy comparison) evaluated at the midplane location = Different locations in the discs are considered later on.
Since we are dealing with an anisotropic system, we look at the onedimensional power spectrum in the three directions: radial, azimuthal and vertical, which are given respectively in the first, second and third column of each row of the two figures.
At least three regions can be identified in the power spectra extracted from simulations. These regions are defined by two wavenumbers: . The smallest scale corresponds to the resolution limit of the simulation that is given by the smoothing length at the considered position (since the discs are axisymmetric there is not dependence on ) and defines the largest wavenumber (for each curve, this is marked by a circle in the plots). The region beyond corresponds to length scales inside the smoothing length, which are below the resolution. Therefore, the resolved region of the spectrum is to the left of the circle. The other scale marks the boundary between regions of the spectrum described by power laws with a different slope: for the slope is close to zero, for the slope is called and is always negative.
In the case of isotropic turbulence one could identify with the ‘forcing scale’ and the region between and with the ‘inertial range’. However, in the present anisotropic case, it is not possible to find a direct link between and the forcing scale in the disc because of the aliasing phenomenon present in onedimensional spectra of turbulent shearstress fluids (see e.g. Pope, 2000), since the power corresponding to wavenumber also contains contribution from larger wavenumbers.
In the following we focus on the qualitative evolution of the power spectrum with changing artificial viscosity and resolution.
Effect of artificial viscosity. What is interesting is the dependence of the shape of the power spectrum on the AV parameters. Simulations with and show a positive slope of the compensated spectra with a decrease just before the resolution scale (), which mimics the decay corresponding to the dissipation scale. In this case a cascade is not clearly identified in any of the three directions. An example is given by simulation S1 with (,)=(0.1,0) in the top row of Fig. 20.
For higher AV parameters the behaviour of the power spectrum depends on the component of the velocity fluctuations, confirming that we are dealing with an anisotropic system also at intermediate scales. In particular, for radial velocity fluctuations in the case of and (e.g. middle and bottom left panels of Fig. 20) and for the azimuthal component when and (bottom central panel of Fig. 20), a progressively more extended flat region (indicating a power spectrum with a slope close to the Kolmogorov value) is visible between the positive slope region on the left and the decay immediately before the resolution scale on the right. The other azimuthal spectra show a positive slope (e.g. top and middle central panels), indicating a subKolmogorov slope. Concerning the vertical component of the velocity, flat regions start to appear for values of close to 1, just above the resolution limit. However, they are less extended than in the case of the other two components. For higher , we observe that larger leads to a steeper slope.
In correspondence to the onset of meridional circulation, we observe the appearance of an approximately flat region in the compensated spectra for all velocity components and at the smallest resolved scales. This result suggests that the system is more isotropic at the small scales above the resolution limit.
Effect of resolution In agreement with previous analyses, we observe a convergence of the power spectrum of velocity fluctuations when a number of particles is used. In addition, higher resolution discs are characterised by a cascade that is globally more extended in the wavenumber space than in lower resolution discs because they resolve smaller length scales, corresponding to larger wavenumbers. The resolution scale therefore shifts towards larger values (see Fig. 21).
In numerical simulations of turbulence a pileup of energy (positive slope of the spectrum) is often observed in the high wavenumber region of the power spectrum (corresponding to the predissipative region, close to grid resolution) and is probably due to numerical dissipation. We observe a similar feature only in the power spectrum of radial velocity fluctuations in intermediate and large viscosity discs (see left panels in Fig. 20 and 21). The pileup is reduced by increasing resolution (see left panels in the top two rows in Fig. 21), supporting the numerical origin of the phenomenon.
Simulations are observed to converge towards a power spectrum with a cascade close to the Kolmogorov one for the radial and azimuthal components. The system reflects its anisotropy in a cascade that is always more extended for radial fluctuations and very short for vertical velocity fluctuations. Such a particular feature of the cascade is due to the different physical scales present in the disc, which is more radially than vertically extended (the disc is thin).
Changing location in the disc. The trends of the power spectrum observed with changing artificial viscosity and resolution are qualitatively similar at different radial locations in the disc (not shown).
In Fig. 22 spectra at different heights above the midplane are compared for each of the three velocity components in the case of simulation S8 where : moving away from the midplane, the slope of the cascade becomes slightly smaller and the cascade disappears due to decreasing resolution. In fact, at higher altitude the gas density is lower, implying a larger smoothing length (i.e. worse resolution). The wavenumber associated to the larger smoothing length is smaller and therefore it is closer to the wavenumber of the ‘forcing scale’, shortening the cascade.
We conclude that the power spectrum preserves its properties throughout the vertical and radial extension of the disc in the resolved regions of the wavenumber space.
6.3 Intermittency
In order to determine if intermittency is present in the simulations we analysed the PDF of the density and of the acceleration in the azimuthal direction (the direction of the shear) and the and order moment of the PDF of several quantities.
6.3.1 Density and azimuthal acceleration PDF
Effect of artificial viscosity. The effects of the AV parameters on the shape of the PDF of the natural logarithm of the density are shown in the left panel of Fig. 23. The distribution of ln is described by a Gaussian distribution. It follows that the density is characterised by a lognormal distribution, as observed in simulations of supersonic compressible turbulence (see e.g. Price & Federrath, 2010, and Sect. 3.2.2).
However, in our simulations the standard deviation of the density logarithm , which controls the width of the distribution, is much smaller, since we are dealing with subsonic physical fluctuations (as shown in Sect. 5.3, Fig. 17 and Eq. 26). In addition, we observe that larger values of and/or lead to a narrower distribution. This behaviour is explained by the same reason: the Mach number of physical fluctuations decreases and leads to a smaller , which in turn produces a PDF close to a Gaussian distribution and gives values of the skewness and kurtosis close to Gaussian values (as we will show in Sect. 6.3.2 and as expected from Eq. 27). The effect of changing for a given is smaller than that of changing for a given , due to the fact that controls the second order term in the MG83 artificial viscosity.
It should be noted here that even for high and/or , although the effective viscosity becomes larger, the Reynolds stress is always non zero and still contributes to it. Therefore physical fluctuations are present and the disc is not laminar. For example, in our discs when the effective viscosity increases to and the Reynolds stress amounts to , with a corresponding turbulent Reynolds number , see Sect. 5.1 and Figs. 12 and 13. If the disc were laminar, one would not expect a distribution of densities but a single value, and therefore a delta function.
In the right column of Fig. 23 we show the PDF of the acceleration in the azimuthal direction since it presents an interesting feature: the distribution shifts towards higher values when both and are increased. This is an effect of the onset of meridional circulation in discs with high enough artificial viscosity. In fact, the presence of meridional circulation imposes a net positive radial velocity (outward flow) to the gas around the midplane that, combined with the disc rotation, implies a spirallike flow. The acceleration is thus characterised by a positive average azimuthal component (in contrast to a pure circular flow where the acceleration is only radial). The peak location moves towards the value =1 for larger AV parameters, however it does not reach unity even for the larger AV combination considered here ()=(5,10) since a fraction of the particles still presents a small but negative azimuthal component of the acceleration. Finally, the PDF of the radial and vertical components of the acceleration (not shown here) present the standard features of a pure circular flow: both of them have a symmetric shape with the peak located around the nonzero Keplerian value or around zero, respectively.
Effect of resolution. The effect of resolution is shown in Figure 24. For all the combinations of AV parameters, more resolved discs are characterised by a narrower density PDF (left column), due to the reduction in the numerical noise. The azimuthal acceleration PDF (right column) is only slightly affected, with the peak shifted towards smaller acceleration for higher resolution, due to a correspondingly lower strength of meridional circulation (see Sect. 4.2). As already found in previous analysed quantities, a number of particles as large as guarantees a good degree of convergence.
In conclusion, the PDFs of density fluctuations reproduce those expected for non intermittent subsonic turbulent flows, confirming that the discs are not laminar.
6.3.2 Higher order moments: S and K
Effect of artificial viscosity. For a given , larger values of result in values of the skewness of the density PDF closer to the Gaussian value (top left plot in Fig. 25). The same effect is observed when is kept constant and increased. All other quantities show Gaussian values independently of the AV parameter used, an example is given by the azimuthal component of the acceleration (top right plot in Fig. 25). Deviations of the kurtosis from the Gaussian value are qualitatively the same as for the skewness.
The simulations with the Mu96 artificial viscosity show systematically lower values of and with respect to simulations with the MG83 artificial viscosity and with the same and . The reason is that for the same , Mu96 discs are more viscous than MG83 discs since artificial viscosity is applied to all particles in contrast to the MG83 case. This trend is in agreement with that observed in MG83 simulations with changing AV parameters: more viscous discs present lower and .
Effect of resolution. The values of S and K only show small fluctuations when the resolution is increased (not shown), suggesting that convergence for these global quantities is present even for a number of particle as low as .
In conclusion, no significant deviation from a Gaussian distribution is observed: intermittency is not present. This result agrees with the intermediate turbulent Reynolds number proper to the present simulations (see Sect. 5.1). In fact intermittency is usually observed in very high Reynolds number flows (e.g. Frisch, 1996).
7 Discussion and conclusion
We have presented the characterisation of the global flow and of the statistical properties of the fluctuations in the velocity and density field of the gas in SPH disc models, which are gaseous discs simulated by means of the SPH method.
An essential tool for most of our results is the new method we have introduced for the determination of the effective viscosity of three dimensional axisymmetric discs. It consists in fitting the analytic expression for the radial velocity derived from two dimensional models (see e.g. TL02, ), depending both on and , to the vertical velocity profile extracted from the simulation at different radial positions in the disc.
We have focused on the effects of changing the number of particles and the values of the two AV parameters and . The relevant results are summarised in the following points:

We have confirmed and quantified the numerical role of the number of particles , which contributes to the numerical component of the SPH fluctuations (see its definition in Sect. 1). In fact, for all studied quantities, an increase in the resolution of the simulation (increase in the number of particles) has the standard numerical role of convergence towards the physical solution, and also leads to a more extended cascade in the power spectrum of velocity fluctuations (smaller scales and therefore higher wavenumbers are resolved).

More significantly, we have found that the artificial viscosity (through the two AV parameters and ) contributes both to the numerical and to the physical component of the SPH fluctuations, in contrast to the behaviour of .
There exists a relationship between artificial viscosity and the following three quantities: (a) physical fluctuations, quantified by , (b) the effective viscosity of the disc, quantified by and (c) the numerical noise. In the hypothesis of an ideal noisefree numerical scheme, we expect that at low the effective viscosity is due to turbulent viscosity, directly produced by physical fluctuations. At intermediate the effective viscosity has a contribution not only from physical fluctuations but also from the AV term, which reproduces the effects of turbulence without directly modeling eddies and vortices in the spirit of the disc model. Finally, in the high range, physical fluctuations are negligible and the effective viscosity is dominated by the AV term. However, numerical schemes are affected by numerical noise, which masks physical fluctuations.
In our case, noise is related to the AV parameters: only for above a threshold is the noise reduced enough to allow physical fluctuations to emerge. We have identified the threshold at (1,2). In fact, in the discs simulated above this threshold, we have observed a sharp transition both of the global and of the local behaviour.

Concerning the global behaviour, we have found that with increasing artificial viscosity the global flow of the gaseous disc evolves from a chaotic radial velocity structure to a more ordered meridional circulation pattern characterised by larger accretion rates and with the Reynolds number , associated to the estimated effective viscosity , below the turbulent limit for all simulations (however we observe that some of the high resolution simulations are approaching it).

Concerning the local behaviour, we observe that, in parallel to the onset of meridional circulation, the average Reynolds stress switches abruptly from positive to negative values and converges to a value that contributes at the 10 per cent level to the effective viscosity (), with an associated turbulent Reynolds number around the turbulent limit. In particular, we have found that with increasing artificial viscosity, physical fluctuations have the following properties:

velocity fluctuations become more and more subsonic and shift from radial to azimuthal anisotropy;

a cascade appears in the power spectrum of velocity fluctuations, which tends towards a Kolmogorovlike spectrum, particularly for the radial component;

turbulent diffusion appears at the same time as the cascade and then the diffusion coefficient tends to decrease. Diffusion tends to be damped by the onset of meridional circulation;

the presence of meridional circulation leads to a shift of the peak of the azimuthal acceleration PDF;

none of the studied SPH disc models presents fluctuations characterised by intermittency.

With these results we can give a first answer to the question raised in the Introduction: can SPH disc models correctly reproduce both the observed properties of PPD and the expected effect of turbulence?
In the case of SPH disc models with the MG83 artificial viscosity we have found that for AV parameters above a threshold that approximately corresponds to the standard values ( and ), the effective viscosity mainly represents viscosity due to turbulence (of unknown origin) in the spirit of the SS73 discs where intermediate and large scale turbulent eddies and vortices are not resolved. However, at the smallest resolved scales, fluctuations with properties similar to those present in turbulent flows have been identified. A meridional circulation pattern and turbulentlike fluctuations coexist in SPH disc models. The former shapes the global scale of the flow, reproducing the observed mass accretion rate onto the star while the latter is present locally at the smallest resolved scales and contributes to only 10 per cent of the effective viscosity of the disc. The resolved physical fluctuations reproduce the main effects expected in a turbulent 3D disc. In fact, the simulated discs are characterised by: a subsonic fluctuating velocity field presenting a turbulentlike diffusion and a small Kolmogorovlike cascade in the power spectrum of velocity fluctuations at the smallest resolved scales (see e.g. simulations S7 and S8). These simulations present values of the Reynolds stress and of the diffusion coefficient that are very close to those observed in MHD accretion discs (where the source of turbulence is the MRI). For values of and below the threshold, the numerical noise dominates and masks physical fluctuations, these models should be avoided.
A direct detection of turbulent eddies would require a reduction of the AV term in order to increase the contribution of the Reynolds stress. However, we have shown that for the two AV implementations considered here the numerical noise dominates and impedes the growth of eddies in the velocity field when the levels of AV are below the threshold: in such cases the positive effects of the AV terms (particularly of the term) in correctly describing SPH particles trajectories and in avoiding particle interpenetration are missing. Therefore, in order to reach higher Reynolds numbers it is necessary both to reduce the effective viscosity of the disc without increasing the numerical noise and to increase the resolution. This is not possible with the two artificial viscosity implementation we are analysing in this work (since reducing the AV coefficients leads to an increase of the numerical noise). Recently developed artificial viscosity switches (e.g Morris & Monaghan, 1997; Cullen & Dehnen, 2010) could help in this direction. However, it is still not possible to simulate fluids at the high Reynolds numbers (Re ) expected in astrophysical flows with the current tools and computer power.
We should stress that simulating very high Reynolds number discs would be necessary if one wanted to investigate the possibility of the existence of pure hydrodynamics turbulence in 3D accretion discs. (We remind the reader that 1D accretion discs are known to be stable with respect to hydrodynamics instabilities due to the Rayleigh criterion, see e.g. Armitage 2007.) However, for our purpose of reproducing the effects of turbulence produced by an unknown source, the effective viscosity combined with a low and possibly intermediate Reynols stress is enough. In fact, these models represent a starting point for the study of dust dynamics, which is affected by physical fluctuations at small scales and meridional circulation at large scales.
In conclusion, in SPH disc models, which a priori include only the ingredients of star gravity and physicallikeviscosity, we have found that the artificial viscosity term, in addition to modeling a physicallike bulk and shear viscosity, can also play the role of an implicit turbulence model. In our SPH disc models, where we do not add any initial turbulent velocity field, the implicit turbulence model sustains and organises the random fluctuations present in the initial conditions of the disc. This result is in agreement with recent indications that a turbulence model is implicitly present in the SPH scheme used to simulate homogeneous and isotropic turbulence with periodic boundary conditions (see e.g. Shi et al., 2012; Monaghan, 2011; Ellero et al., 2010). Here we have analysed the more complex case of systems without boundary conditions and where anisotropic turbulence is expected.
The effects of the density and sound speed profiles on the global and statistical properties of the gas as well as those of the initial set up, of different SPH kernels and of different artificial viscosity implementations will be addressed in more detail in a future work.
Acknowledgments
The authors wish to thank Roland Speith and Guillaume Laibe for fruitful discussions. This research was supported by the Agence Nationale de la Recherche (ANR) of France through contract ANR07BLAN0221 and was conducted within the Lyon Institute of Origins under grant ANR10LABX66. Simulations presented in this work were run at the Service Commun de Calcul Intensif (SCCI) de l’Observatoire de Grenoble, France. Figure 1 was made using the splash software package (Price, 2007).
References
 Andrews et al. (2009) Andrews S. M., Wilner D. J., Hughes A. M., Qi C., Dullemond C. P., 2009, ApJ, 700, 1502
 Armitage (2007) Armitage P. J., 2007, arXiv:0701485
 Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
 BarrièreFouchet et al. (2005) BarrièreFouchet L., Gonzalez J.F., Murray J. R., Humble R. J., Maddison S. T., 2005, A&A, 443, 185
 Bertin (2000) Bertin G., 2000, Dynamics of Galaxies
 Carballido et al. (2008) Carballido A., Stone J. M., Turner N. J., 2008, MNRAS, 386, 145
 Cullen & Dehnen (2010) Cullen L., Dehnen W., 2010, MNRAS, 408, 669
 Cuzzi et al. (1993) Cuzzi J. N., Dobrovolskis A. R., Champney J. M., 1993, Icarus, 106, 102
 Cuzzi et al. (2008) Cuzzi J. N., Hogan R. C., Shariff K., 2008, ApJ, 687, 1432
 Ellero et al. (2010) Ellero M., Español P., Adams N. A., 2010, Phys. Rev. E, 82, 046702
 Flock et al. (2012) Flock M., Dzyurkevich N., Klahr H., Turner N., Henning T., 2012, ApJ, 744, 144
 Flock et al. (2011) Flock M., Dzyurkevich N., Klahr H., Turner N. J., Henning T., 2011, ApJ, 735, 122
 Frisch (1996) Frisch U., 1996, Turbulence
 Fromang et al. (2011) Fromang S., Lyra W., Masset F., 2011, A&A, 534, A107
 Fromang & Nelson (2006) Fromang S., Nelson R. P., 2006, A&A, 457, 343
 Fromang & Nelson (2009) Fromang S., Nelson R. P., 2009, A&A, 496, 597
 Fromang & Papaloizou (2006) Fromang S., Papaloizou J., 2006, A&A, 452, 751
 Guilloteau et al. (2012) Guilloteau S., Dutrey A., Wakelam V., Hersant F., Semenov D., Chapillon E., Henning T., Piétu V., 2012, A&A, 548, A70
 Hartmann et al. (1998) Hartmann L., Calvet N., Gullbring E., D’Alessio P., 1998, ApJ, 495, 385
 Hughes et al. (2011) Hughes A. M., Wilner D. J., Andrews S. M., Qi C., Hogerheijde M. R., 2011, ApJ, 727, 85
 Jacquet (2013) Jacquet E., 2013, A&A, in press (arXiv:1301.5817)
 King et al. (2007) King A. R., Pringle J. E., Livio M., 2007, MNRAS, 376, 1740
 Klahr & Bodenheimer (2003) Klahr H. H., Bodenheimer P., 2003, ApJ, 582, 869
 Laibe et al. (2012) Laibe G., Gonzalez J.F., Maddison S. T., 2012, A&A, 537, A61
 Lattanzio et al. (1985) Lattanzio J. C., Monaghan J. J., Pongracic H., Schwarz M. P., 1985, MNRAS, 215, 125
 Lodato (2008) Lodato G., 2008, New Astronomy Reviews, 52, 21
 Lodato & Price (2010) Lodato G., Price D. J., 2010, MNRAS, 405, 1212
 Lovelace et al. (1999) Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805
 LyndenBell & Pringle (1974) LyndenBell D., Pringle J. E., 1974, MNRAS, 168, 603
 Meglicki et al. (1993) Meglicki Z., Wickramasinghe D., Bicknell G. V., 1993, MNRAS, 264, 691
 Meru & Bate (2012) Meru F., Bate M. R., 2012, MNRAS, 427, 2022
 Monaghan (1992) Monaghan J. J., 1992, ARA&A, 30, 543
 Monaghan (2005) Monaghan J. J., 2005, Reports on Progress in Physics, 68, 1703
 Monaghan (2011) Monaghan J. J., 2011, European Journal of Mechanics  B/Fluids, 30, 360
 Monaghan & Gingold (1983) Monaghan J. J., Gingold R. A., 1983, Journal of Computational Physics, 52, 374
 Monaghan & Lattanzio (1985) Monaghan J. J., Lattanzio J. C., 1985, A&A, 149, 135
 Morris & Monaghan (1997) Morris J. P., Monaghan J. J., 1997, Journal of Computational Physics, 136, 41
 Murray (1996) Murray J. R., 1996, MNRAS, 279, 402
 Nelson et al. (2012) Nelson R. P., Gressel O., Umurhan O. M., 2012, ArXiv eprints
 Pope (2000) Pope S. B., 2000, Turbulent Flows. Cambridge University Press
 Price (2007) Price D. J., 2007, PASA, 24, 159
 Price (2012) Price D. J., 2012, MNRAS, 420, L33
 Price & Federrath (2010) Price D. J., Federrath C., 2010, MNRAS, 406, 1659
 Pringle (1981) Pringle J. E., 1981, Annual review of astronomy and astrophysics, 19, 137
 Rice et al. (2011) Rice W. K. M., Armitage P. J., Mamatsashvili G. R., Lodato G., Clarke C. J., 2011, MNRAS, pp 1535–+
 Rice et al. (2005) Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS, 364, L56
 Sargent et al. (2009) Sargent B. A., Forrest W. J., Tayrien C., McClure M. K., Watson D. M., Sloan G. C., Li A., Manoj P., Bohac C. J., Furlan E., Kim K. H., Green J. D., 2009, ApJS, 182, 477
 Schäfer (2005) Schäfer C., 2005, PhD thesis, EberhardKarls Universität Tübingen
 Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, in H. Bradt & R. Giacconi ed., X and GammaRay Astronomy Vol. 55 of IAU Symposium, Black Holes in Binary Systems: Observational Appearances. pp 155–+
 Shi et al. (2012) Shi Y., Ellero M., Adams N. A., 2012, Phys. Rev. E, 85, 036708
 Takeuchi & Lin (2002) Takeuchi T., Lin D. N. C. ., 2002, ApJ, 581, 1344
 Violeau & Issa (2007) Violeau D., Issa R., 2007, International Journal for Numerical Methods in Fluids, 53, 277
Appendix A Diagnostic details
a.1 The turbulent viscosity coefficient
From the turbulent viscosity hypothesis (Pope, 2000), the Reynolds stresses are written in Cartesian coordinates as follows:
(28) 
where is the rate of strain:
(29) 
and the turbulent kinetic energy:
(30) 
In classical accretion disc theory, only shear viscosity is relevant. Therefore the only nonvanishing component of the stress tensor is the component (see e.g Lodato, 2008), which in cylindrical coordinates becomes:
(31) 
Given the axisymmetry of the disc, the second term on the right hand side of Eq. 31 vanishes and the turbulent viscosity coefficient becomes:
(32) 
with .
a.2 The coefficient
The SS73 coefficient corresponding to the viscosity of Sect. A.1 is defined by the relation that combined with Eq. 32 gives:
(33) 
In the case of quasiKeplerian discs we have: = , with the Keplerian velocity and the corresponding angular velocity. For thin discs, the approximation holds, with the radial spherical component and the radial cylindrical component.
Remembering that the scale height of the disc is related to the sound speed and to the Keplerian angular velocity by , the coefficient takes the form: