Evidence for Non-zero Turbulence in the Protoplanetary disc around IM Lup

The amount of turbulence in protoplanetary discs around young stars is critical for determining the efficiency, timeline, and outcomes of planet formation. It is also difficult to measure. Observations are still limited, but direct measurements of the non-thermal, turbulent gas motion are possible with the Atacama Large Millimeter/submillimeter Array (ALMA). Using CO(2-1)/$^{13}$CO(2-1)/C$^{18}$O(2-1) ALMA observations of the disc around IM Lup at ~0.4"(~60 au) resolution we find evidence of significant turbulence, at the level of $\delta v_{\rm turb}=(0.18-0.30)$c$_s$. This result is robust against systematic uncertainties (e.g., amplitude flux calibration, midplane gas temperature, disc self-gravity). We find that gravito-turbulence as the source of the gas motion is unlikely based on the lack of an imprint on the rotation curve from a massive disc, while magneto-rotational instabilities and hydrodynamic instabilities are still possible, depending on the unknown magnetic field strength and the cooling timescale in the outer disc.


INTRODUCTION
Turbulence is a primary factor in the planet formation process, affecting chemical evolution (Semenov & Wiebe 2011;Furuya & Aikawa 2014;Xu et al. 2017), collisional growth of dust grains (Ormel & Cuzzi 2007;Birnstiel et al. 2010), and the vertical/radial concentration of dust grains (Dullemond et al. 2018;Rosotti et al. 2020;Jennings et al. 2022).Turbulence affects the structures generated in a disc by a newly formed giant planet (Bae et al. 2018), and any uncertainty in the turbulence level influences our ability to interpret the myriad of complex structures that have been revealed among many protoplanetary discs (e.g., Andrews et al. 2018;Long et al. 2018;Öberg et al. 2021).As such, understanding the strength of turbulence, and how it varies from system to system, is important in understanding the conditions under which planets form and evolve.
Often quantified as  in the context of the -disc model of viscosity ( =    where  is the viscosity,   is the sound speed, and  is the pressure scale height Shakura & Sunyaev 1973), early studies of protoplanetary disc evolution inferred a global value of  ∼ 10 −2 based on the measured accretion rate onto the central star (Hartmann ★ E-mail: kmf4@williams.edu(KMF) et al. 1998).More recent observations of dust vertical/radial diffusion (Pinte et al. 2016;Dullemond et al. 2018;Rosotti et al. 2020;Jennings et al. 2022;Ohashi & Kataoka 2019;Ueda et al. 2021;Doi & Kataoka 2021;Franceschi et al. 2022;Villenave et al. 2022), disc sizes (Najita & Bergin 2018;Trapman et al. 2020;Long et al. 2022), and the relation between accretion rate and disc mass/radius as predicted by viscous evolution (Ansdell et al. 2018;Ribas et al. 2020) have found  = 10 −4 − 10 −3 , suggesting more modest levels of turbulence, at least in the outer parts of the disc 1 .
Dust diffusion and angular momentum exchange are valuable but indirect methods for measuring turbulent gas motion.The Doppler shifts of molecular line emission directly traces the kinematics of the gas.Early results (Guilloteau et al. 2012;Hughes et al. 2011) were limited in spatial resolution and sensitivity, but ALMA has dramatically improved on this capability.Observations of the discs around HD 163296 (Flaherty et al. 2017), TW Hya (Teague et al. 2018;Flaherty et al. 2018), and V4046 Sgr and MWC 480 (Flaherty et al. 2020) placed stringent upper limits on the amount of turbulence, constraining the non-thermal velocity component to <0.05c  -<0.15c  between different targets ( less than a few ×10 −3 , assuming  ∼ ( turb /  ) 2 and that  turb /  is constant throughout the disc).The exception to the general trend of weak turbulence was DM Tau, which displayed turbulence of (0.18 -0.28)c  (Flaherty et al. 2020).
Why some systems are turbulent while others are not depends on the physical mechanism(s) that can drive turbulence within a protoplanetary disc.Gravito-turbulence (Shi & Chiang 2014;Forgan et al. 2012) relies on gravitational instabilities within the disc, and hence different turbulence velocities could reflect different disc masses.The magneto-rotational instability (MRI, Balbus & Hawley 1998) relies on a coupling between magnetic fields and partially ionized gas, and variations in the turbulent velocities may reflect either differences in the magnetic field geometry and/or strength or in the ionization of the disc gas (Simon et al. 2018).Hydrodynamic instabilities, e.g., the Vertical Shear Instability (Nelson et al. 2013) which relies on the change in orbital velocity with height within the disc, depend on the cooling timescale (Lyra & Umurhan 2021;Lesur et al. 2022), which in turn may reflect the properties of the dust population in the outer disc.
Molecular line studies are limited in size, and often focus on the brightest discs given the high SNR required to measure the ∼100 m s −1 non-thermal motions characteristic of turbulence in protoplanetary discs.This limits the number of sources in which we can distinguish between turbulence and no/weak turbulence.With this caveat in mind, DM Tau does stand out from the rest of the sources in the molecular line sample as being the youngest system, hinting at a connection between turbulence and age.A change in turbulence with age may reflect changes in disc mass (Manara et al. 2022), magnetic field strength (Simon et al. 2018;Weiss et al. 2021), or cooling timescale, depending on the model for turbulence.A decrease in turbulence with age could affect the turbulence that the planets encounter as they migrate through the gas rich disc (Hühn et al. 2021), with weaker turbulence at later times leading to more compact planetary systems with orbits close to resonances.
Here we examine the young source IM Lup, around which we find evidence of non-zero turbulence.At an age of ∼1 Myr (Alcalá et al. 2017;Mawet et al. 2012) it is comparable in age to DM Tau.It also has a large CO disc; Long et al. (2022) find that DM Tau (876 au) and IM Lup (803 au) have two of the three largest CO discs among the 44 targets they study, with the majority of discs having radii between 50 and 350 au.Such large discs, when combined with the young ages, are consistent with substantial viscous evolution (Trapman et al. 2020;Long et al. 2022), which may be a sign of turbulence.More recently, Paneque-Carreño et al. (2023) examine the non-thermal motion with the CN emission line and find evidence of turbulence at (0.4-0.6)c  .
In section 2 we describe the data and models used to constrain turbulence, while we discuss the finding of non-zero turbulence around IM Lup in section 3.In section 4 we tackle the question of what could be driving the turbulence around DM Tau and IM Lup, and what makes these systems different from the other sources, in order to better understand the expected range of turbulence levels among protoplanetary discs, and to constrain the physical mechanism driving the turbulence.

DATA & MODEL
We use data from project 2013.1.00798.S (PI: C. Pinte), originally presented in Pinte et al. (2018).These data consist of six spectral windows, and in this work we focus on the spectral windows centered on CO(2-1), 13 CO(2-1), and C 18 O(2-1).These spectral windows are centered at 230.54 GHz,220.394 GHz,and 219.556GHz respectively, with a channel spacing of 30.518 kHz, corresponding to a velocity spacing of 39 m s −1 , 41 m s −1 , and 42 m s −1 for the three emission lines.
The raw data were run through the CASA pipeline v4.7.2, while analysis was performed in CASA v5.4.0.Self-calibration was performed using a continuum spectral window with two rounds of phase calibration.For CO(2-1), images are generated using Briggs weighting with robust=0.5,resulting in a beam size of 0. ′′ 3×0.′′ 5 and an rms of 8.6 mJy bm −1 channel −1 .Natural weighting was used for 13 CO and C 18 O resulting in a beam size of 0. ′′ 4×0.′′ 7 for both lines, with an rms of 7.9 mJy bm −1 channel −1 and 5.6 mJy bm −1 channel −1 for 13 CO and C 18 O respectively.
To model the molecular line emission, we use the parametric surface density and disc temperature structure as described in Flaherty et al. (2020) and references therein.To briefly summarize, the surface density is assumed to follow a power law with an exponential tail: where  gas ,   , and  are the gas mass (in M ⊙ ), critical radius (in au) and power law index respectively.The disc extends from  in to 1000 au.The temperature structure is parameterized as a power law with radius, with a smooth function connecting the cold midplane and the warm atmosphere. (2) The parameter   is the height above the midplane at which the gas temperature reaches its maximum value, and  0 was set to twice the local pressure scale height (Dartois et al. 2003), with the pressure scale height defined as: A hydrostatic equilibrium calculation is performed at each radius to derive the gas volume density at each height above the midplane, based on the temperature and surface density structure.CO is spread throughout the disc assuming a constant abundance.In regions with   2 < 1.3 × 10 21 cm −2 , where   2 is the column density of H 2 , CO is assumed to be photo-dissociated (Qi et al. 2011) and the abundance in this region is dropped by eight orders of magnitude.Similarly, the CO abundance is dropped by a factor of five in regions with  gas < 19 K to account for freeze out.We assume isotope abundances of 18 O/ 16 O=557 and 13 C/ 12 C=69 (Wilson 1999) when modeling 13 CO and C 18 O emission.
Extended CO emission has been detected around IM Lup, at a distance where the disc gas temperatures would suggest that CO should be frozen out onto dust grains, indicating that photodesorption returns a substantial mass of CO to the gas phase (Cleeves et al. 2016;Seifert et al. 2021).To accommodate this additional CO, we include CO in the outer disc following the prescription in Flaherty et al. (2020): in the region of the disc where the gas temperature is below the CO freeze-out temperature, we include CO at a constant abundance in the region with 1.3 × 10 21 cm −2 <   2 < 4.8 × 10 21 cm −2 .This prescription maintains CO freeze-out at the midplane, while allowing CO gas in the upper layers of the outer disc where high-energy photons can dislodge CO from ice mantles on dust grains.Cleeves et al. (2016) find evidence for a factor of 20 depletion in CO abundance around IM Lup relative to the standard ISM value, while Zhang et al. (2019) find that the CO depletion is a factor of 100 in the inner disc (<100 au) and is a factor of 5 in the outer disc.We assume a uniform CO depletion of a factor of 20 across the entire disc in our models of the CO emission (i.e., CO/H 2 = 5×10 −6 ), although we also consider models in which the CO has an abundance equal to the ISM value when fitting all three lines together.
The velocity profile is assumed to be Keplerian, with additional corrections for the pressure support of the gas and the height above the midplane.The line profile is assumed to be a Gaussian whose width is set by thermal and non-thermal motions.We assume the nonthermal term is proportional to the local isothermal sound speed: , where  turb is in units of c  .Throughout much of this paper we assume that the non-thermal line broadening is associated with turbulence, although we cannot rule out non-Keplerian non-thermal motion below the spatial resolution limit.In section 4 we discuss such non-Keplerian non-thermal terms in the context of specific physical models (e.g., vertical motions associated with the VSI), and the effect of these types of motion on our turbulence measurement.
We assume M * = 1.1 M ⊙ (Alcalá et al. 2017;Öberg et al. 2021).We use a distance of 158.5 pc from Gaia DR2 measurements of the parallax (Bailer-Jones et al. 2018) and a gas mass of 0.175 M ⊙ (Pinte et al. 2018).During each trial we vary ,   ,  turb ,  atm0 ,  mid0 , inclination (),  in , systemic velocity ( sys ), RA offset from phase centre ( off ), Dec offset from phase centre ( off ), and position angle (PA).For a given set of model parameters, and the resulting density, temperature, velocity structure, model images are generated by solving the radiative transfer equation, as described in Rosenfeld et al. (2013a).Images are generated at the same velocities as the data, with Hanning smoothing applied to the model images.
The posterior distributions for each parameter are estimated using the MCMC routine EMCEE (Foreman-Mackey et al. 2013), based on the Affine-Invariant algorithm originally proposed in Goodman & Weare (2010).This methodology allows us to simultaneously constrain the temperature, density, and turbulence in the disk, accounting for any degeneracies between the parameters.The likelihood of a given model is calculated using the visibilities, with the model visibilities derived from model images using galario (Tazzari et al. 2018).
where the summation is performed over every  and  position and over every channel, and   , is a visibility point in the data at a given  and  and channel, while  , is the visibility point at the same , , and channel derived from the model image.Galario has previously been applied to continuum images, and to extend it to spectral data we applied galario to each channel individually to derive the visibilities.The uncertainties in the data are estimated based on the calculated dispersion among baselines of similar distances in line-free channels.The region that is subject to absorption by the molecular cloud, and is excluded from the MCMC process, is marked with a grey band.Despite the limited spectral range, we are able to place strong constraints on the non-thermal linewidth (  turb =0.237 +0.017 −0.012 c  ), with a significantly better fit to the data than with zero turbulence (green dotted line).
Absorption from the molecular cloud obscures some of the CO emission (van Kempen et al. 2007;Panić et al. 2009) and when fitting to CO(2-1) we exclude channels with  LSR between 4 and 6 km s −1 (Figure 1).Previous studies of turbulence using molecular line emission have focused on high SNR, un-contaminated targets, but expanding the sample requires the consideration of less ideal data, including objects with emission contaminated by cloud absorption similar to what is seen around IM Lup.Part of the goal of our study is to understand how this cloud absorption affects our ability to constrain turbulence in the disk around IM Lup.This absorption is not evident in 13 CO or C 18 O and we use the entire line profile when fitting these lines.Typical MCMC chains consist of 50 walkers and 1000 steps, with convergence on the final solution occurring within 300 steps.The first 500 steps are removed as burn-in, much longer than the typical auto-correlation time, which converges towards ∼60 steps towards the end of the chains.

RESULTS
With our fiducial model, fit to the CO(2-1) data, excluding the region with  LSR between 4 and 6 km s −1 , we find turbulence is nonzero, at  turb =0.237 +0.017 −0.012 c  (Table 1).Figure 1 shows the spectrum of the data and the model defined by the median of the posterior distribution functions (PDFs), while Figure 2 shows the residuals between the model and the data.The model is an excellent match to the data, with some small positive residuals towards the midplane of the disc.We find that we are able to constrain turbulence despite the exclusion of some of the channels.This is because each channel contains information about the turbulence, which we are able to utilize by modeling the full three-dimensional data set.
Figure 3 shows the temperature and height of the CO emitting regions as derived in our fiducial model, as compared to the height and temperature derived by Law et al. (2021).The temperature and height of the CO emitting regions were derived by considering the  distribution of temperatures and heights that contribute to the lines of sight that pass through a given radial bin; the regions in Figure 3 show the 25th and 75th percentile points on this distribution at each radius, representing the boundaries within which the majority of the emission at a given radius arises.The breadth of this region represents the range of optical depths for the lines of sight at a given radius, as well as the front-and back-sides of the disk.As seen in the bottommiddle and bottom-right panels of Figure 4, while the majority of the emission is optically thick and comes from a consistent  above the midplane, some emission arises from deeper in the disk.We are consistent with the Law et al. ( 2021) results, although we derive a higher temperature in the inner disc, where the Law et al. measurements are biased by beam dilution.A comparison of temperatures derived with these different data sets will also be subject to the 10-20% uncertainty in the amplitude calibration (Butler 2012).Given the temperatures of the CO emitting region, the turbulence constraint ( turb =0.237 +0.017 −0.012 c  , c  = √︁ 2/ ℎ ) corresponds to 87 m s −1 at 100 au and 65 m s −1 at 500 au2 .These velocities are consistent with the models of Cleeves et al. (2016), which find that turbulence is weaker than 200 m s −1 (0.5 -0.7 c  ).Our results are slightly weaker than the (0.4-0.2023) assume a gas temperature from the intensity of optically thick CO or 13 CO and compare this to the broadening observed in CN.
While the statistical uncertainty on the turbulence is small within the fiducial model, systematic effects amplify this uncertainty.Accounting for amplitude calibration uncertainty, by modeling the emission with a ±20% offset in total flux (models High sys and Low sys in Table 1), we find a statistically significant difference in  atm0 (29.2 -39.1 K) and  turb (0.18 -0.28 c  ).Increasing   from 2 times the pressure scale height to 4 times the pressure scale height (  = 4 in Table 1) also results in a turbulence that is statistically significantly different from the fiducial model ( turb = 0.30 +0.01 −0.02 c  ) but also does not change our overall interpretation that there is non-zero turbulence around IM Lup.The stellar mass also influences the velocity pattern, and including  * as a free parameter (Stellar Mass in Table 1) does not substantially change the turbulence ( turb = 0.245 +0.010 −0.018 c  ), while finding a stellar mass ( * = 1.088 +0.007 −0.006 M ⊙ ) that is consistent with analysis of the velocity profile ( * ≈ 1.02 M ⊙ Lodato et al. 2023).
Previous work has found that the ratio of the peak flux to the flux at line centre varies with the turbulence level, with smaller peakto-trough ratios associated with stronger turbulence (Simon et al. 2015;Pinte et al. 2022).While the cloud contamination prevents us from accurately measuring the peak-to-trough ratio, we are still able to robustly measure the turbulence, albeit with larger systematic uncertainties.In the next two subsections we examine in detail the influence of midplane temperature and CO abundance on turbulence, but the general result still holds: we find robust evidence for non-zero turbulence in the disc around IM Lup, at the level of (0.18 -0.30)c  .

Midplane Temperature
Our result predicts a CO snowline (T gas =19 K) at the midplane at 154 +15 −12 au, while Qi et al. (2019) find an upper limit on the CO snowline radial location of 59 au based on the innermost edge of the N 2 H + emission.Seifert et al. (2021), building on models from Cleeves et al. (2016), also find a CO snow line that is smaller than our model predictions, while Zhang et al. (2021) put the midplane CO snowline at 15±5 au.The CO emission studied here is optically thick and does not directly trace the midplane temperature.Even the emission from the far side of the disc arises from slightly above the midplane (Figure 3).In this way, the midplane temperature does not directly affect the CO emission.But the midplane temperature does indirectly affect our models through the gas pressure scale height, which roughly scales with the square root of the midplane temperature.Since the height of the emitting region is constrained by the data, the midplane temperature is indirectly constrained and the large height of the CO emitting surface (Pinte et al. 2018;Law et al. 2021) results in a large midplane temperature, which pushes the CO snowline in our model to large radii.The potentially over-estimated midplane temperature may also reflect missing features in the parametric model, such as a jump in the midplane temperature beyond the outer edge of the dust disc (Cleeves 2016), or an increasing CO abundance with radius (Zhang et al. 2021).Flaherty et al. (2018) find that the choice of midplane temperature can bias the turbulence measurement, with an over-estimate of the midplane temperature leading to an underestimate of the turbulence.To understand this effect in the disc around IM Lup, we run a trial with T mid0 set such that the snow line occurs at 59 au (i.e., T mid0 = 19×(59/150) − ), and a second trial where the midplane snowline is at 19 au.As with the fiducial model, we find non-zero turbulence, at a level of  turb =0.25 +0.01 −0.02 c  and 0.256 +0.018 −0.010 c  for the models with the snowline at 59 au and 15 au respectively(R  = 59 and R  = 15 in Table 1).This is consistent with the findings of Flaherty et al. (2018) that a lower midplane temperature is associated with stronger turbulence, although here the effect is at the level of the uncertainty in the turbulent velocity.This indicates that midplane temperature does not strongly affect our estimate of the turbulence.

CO abundance
While midplane temperature does not substantially affect our measurement of non-zero turbulence, we do find that the global CO abundance plays an important role in the measured turbulence level.When running a trial with no CO depletion (i.e., CO/H 2 =10 −4 ) outside of CO freeze-out, we find  turb <0.03c  , much weaker than when CO depletion is included.This anti-correlation between turbulence and CO abundance is similar to that found when using more detailed chemical models of CO depletion (Yu et al. 2017).The connection between turbulence and CO abundance exists because both contribute to emission in the same regions of the channel maps.Figure 4 shows an individual channel from models with (1) low CO abundance and high turbulence, (2) low CO abundance and no turbulence, and (3) high CO abundance with no turbulence.Both the low abundance and high turbulence model (upper left panel) as well as the high abundance and low turbulence model (upper right panel) generate a broadening of the emission in the channel maps.This broadening is not present in the low abundance, zero turbulence model image (upper centre panel).This arises because even though CO is, overall, highly optically thick, there are lines of sight that are optically thin, making the emission along these lines of sight sensitive to the abundance.These lines of sight are the same lines of sight that are most strongly influenced by the turbulence, as they come from regions that have bulk motions at the edge of the velocity channel, which get brought into the channel by turbulent motion.This effect is larger than that seen around TW Hya, where increasing the CO abundance by a factor of 10 led to a change in the upper limit from <0.13c  to <0.08c  (Flaherty et al. 2018).This may be because the cloud contamination for IM Lup excludes the central channels.As shown previously (Simon et al. 2015), the peakto-trough ratio is sensitive to the turbulence level, through its impact on the central velocity channel, and without this information we are more susceptible to the degeneracy with CO abundance.
Without access to the central channels, other data is needed to break the degeneracy between turbulence and CO abundance.The CO(2-1) emission by itself cannot constrain the CO abundance, but when modeled in concert with more optically thin emission lines like 13 CO(2-1) and C 18 O(2-1) we can constrain the turbulence and CO abundance simultaneously.To this end, we simultaneously model CO(2-1), 13 CO(2-1) and C 18 O(2-1), while adding the CO abundance as a free parameter.For simplicity we assume a single turbulence scaling relative to the local sound speed, although models of MRI (e.g.Simon et al. 2015Simon et al. , 2018) ) and VSI (Flock et al. 2017) indicate that turbulence likely varies with height within the disc, beyond a simple scaling with the local sound speed.We use a common set of physical parameters to create a separate set of model visibilities for all three isotopologues and sum the log-likehoods from all three (ln   = ln   + ln  13 + ln  18 ).No additional weighting is applied to the log-likelihood values from each isotopologue, although the result is likely weighted toward CO(2-1) given its higher SNR.Results are shown in Table 1 (Multi-line) and Figure 5; we find non-zero turbulence ( turb =0.20 +0.01 −0.02 c  ), with a strong CO depletion (CO/H 2 = 4.1 +0.2 −0.3 × 10 −6 ).The temperatures and heights of the emitting regions for these molecules are consistent with those values derived by Law et al. (2021) (Figure 6).The CO abundance derived here is consistent with the factor of 20 depletion (i.e., CO/H 2 = 5×10 −6 ) derived by Cleeves et al. (2016) and the factor of 10 -100 depletion derived by Zhang et al. (2021) for this system.
For simplicity we assumed a constant CO depletion across the disc, although observations and models suggest that the CO depletion may vary with radius (Zhang et al. 2019(Zhang et al. , 2021)), with height (Krĳt et al. 2018), and with isotopologue (Miotello et al. 2014).Zhang et al. (2021) find that the CO depletion is a factor of 100 at radii less than 100 au, while the CO depletion is only a factor of 5 at larger radii, bracketing the CO depletion we find, suggesting that our result is best interpreted as an intensity-weighted average CO depletion throughout the disk.Given the degeneracy between CO abundance and turbulence, a spatially-varying CO abundance will most strongly affect efforts to constrain spatially-varying turbulence, which are beyond the scope of this paper.
The multi-line fit converges on a turbulence level that is smaller than the single line fit.This may be due to a decrease in the turbulence level towards the midplane; if the 13 CO(2-1) and C 18 O(2-1) lines are probing regions with weaker turbulence (relative to the local sound speed) then fitting these lines will drag down the average turbulence There is a degeneracy between turbulence and CO abundance, such that low turbulence can be matched with a higher CO abundance.This can be seen in comparing the morphology of the left and right panels, which is similar, despite the differences in abundance and turbulence.Bottom Row: IM Lup data, followed by maps of optical depth (dotted line corresponds to  = 1) and height of the  = 1 surface (dotted line corresponds to z=0) for the fiducial model ([CO/H 2 ]=5×10 −6 ,  turb = 0.24c  ).CO(2-1) is sensitive to the CO abundance, despite being mostly optically thick, because there are lines of sight that are optically thin (i.e., outside of the region marked by the dotted lines in the central panel) that reach deep into the disc.level.As noted above, this behavior is expected from both MRI and VSI.Despite the similar vertical behavior, VSI may generate unique and detectable radial structure in the channel maps, as discussed in more detail below.Modeling the 13 CO(2-1) line by itself ( 13  only in Table 1), or the C 18 O(2-1) line by itself ( 18  only in Table 1), results in turbulence levels consistent with that derived from CO(2-1) by itself (0.237 +0.017 −0.012 c  for CO(2-1), 0.23±0.04c  for 13 CO(2-1), 0.25 +0.07 −0.15 for C 18 O(2-1)), albeit with larger uncertainties even before accounting for systematic effects.Further analysis is needed to fully characterize any evidence for a vertical gradient in the turbulence from the line emission.

Scattered light
We can also characterize the turbulence level using the scattered light images of Avenhaus et al. (2018).The height of the scattered light surface is related to the dust scale height, which in turn is set by turbulence.The  value is related to the dust scale height   , and the gas pressure scale height, , via (e.g., Dullemond et al. 2018) 3 : where  = √︁ /.The Stokes number, , is defined as where  is the particle size,   is the bulk density of the material (=2 g cm −2 ), and Σ  is the gas density.
To estimate   /  from the SPHERE observations we generate synthetic scattered light images using RADMC-3D (Dullemond et al. 2012).We utilize the same density and temperature structure as in the fiducial CO model (Fiducial in Table 1).Initially we consider the dust as a population of 0.1m amorphous silicate grains spread uniformly throughout the radial extent of the disk.The gas to dust mass ratio is assumed to be 100, and the volume density of the dust has a Gaussian profile in the vertical direction with a scale height that is a fraction of the gas scale height (i.e., H  /H is constant throughout the disk), where the gas scale height is approximated by H  =c  /Ω, where c  is defined by the midplane temperature.Scattered light images are generated assuming isotropic scattering at the wavelength of the observations from Avenhaus et al. (2018).The scattered light images are deprojected using diskmap (Stolker et al. 2016) which assumes that the height of the scattering surface follows a power law with radius (z scatter = z 0 (r/ 1 au)  ).Under the assumption of isotropic scattering, a deprojection that uses a combination of z 0 and  that accurately reflects the height of the scattered light surface will produce emission that is axisymmetric about the center of the disk (i.e., radial profiles in any direction away from the center of the disk will be identical).For a given H  /H  , which we vary between 0.3 and 1, we search for a combination of z 0 and  that minimizes the asymmetry of the emission, and exclude values of H  /H  for which the combination of z 0 and  that minimize the asymmetry are inconsistent with the height of the rings observed by Avenhaus et al. (2018).
We find that   /  ≈ 0.7 − 1 (Figure 7), indicating that  ⪆ 1.While our methodology includes some limiting assumptions, they are consistent with a large H  /H  .Dust coagulation models, accounting for growth, fragmentation, and transport, predict a continuous distribution of grain sizes (Birnstiel et al. 2012).Extending our modeling to include a mix of 0.1 and 100m grains, with the large grains have a scale height ten times smaller than the small dust grains, reduces the height of the scattered light surface for a given H  /H  .As a result, including a more realistic grain size distribution would push our result towards higher H  /H  , making our lower limit of 0.7 a conservative estimate.
Given the asymptotic behavior of  as H  /H  approaches 1 ( → ∞ as H  /H  →1), taking   /  = 0.7 provides a lower limit on , which itself scales linearly with the Stokes number.For a broad range of Stokes numbers, the lower limit on  is consistent with the turbulence constraint from the CO emission.
When dust is highly coupled to the gas ( → 0 or  → ∞) the scale height of the dust will depend on the turbulence throughout the vertical extent of the disc.When dust is not perfectly coupled to the gas, dust settling is more strongly dependent on the turbulence at the midplane (Ciesla 2007), even though we are observing light scattering off of small dust grains in the atmosphere of the disk.CO, on the other hand, traces turbulence in the surface layers (Figure 3).Models of MRI and VSI indicate that turbulent velocities can decrease toward the midplane by up to a factor of ten (Simon et al. 2015;Flock et al. 2017).A factor of ten lower turbulence is indicated in Figure 7, and the consistency between this velocity and the scattered light image depends on the size of the dust grains in the disk atmosphere that are responsible for the scattering.Qi et al. (2019), in modeling the SED of IM Lup, find a maximum grain size of 3m for the disc atmosphere, while Tazaki et al. (2023) find that the colour and polarization of scattered light images are consistent with fractal dust aggregates larger than 2m.Franceschi et al. (2022) use the scattered light emission in combination with the continuum emission from ALMA to constrain  and the maximum grain size, finding results consistent with the level of turbulence we derive from CO, under the assumption that turbulence decreases towards the midplane.

DISCUSSION: WHAT IS DRIVING TURBULENCE?
In Flaherty et al. (2020) we found nonthermal gas motion, between 0.25 c  and 0.33 c  , around DM Tau.Here we have examined IM Table 1.2022), who simultaneously constrain  and the maximum grain size.These constraints on the dust settling are also consistent with a decrease in the turbulence between the surface layers probed by CO, and the midplane probed by dust settling.
Lup and found similar non-zero nonthermal gas motion at  turb = 0.18−0.30c  .These results stand in contrast to the non-detections of turbulence around MWC 480, V4046 Sgr, TW Hya, and HD 163296, with upper limits ranging from  turb < 0.15 c  down to  turb <0.05 c  (Flaherty et al. 2018;Teague et al. 2018;Flaherty et al. 2020)  directly derived from CO or CN emission are tracing the behavior in the upper layers of the outer disc rather than the midplane traced by CO diffusion, where turbulence is likely lower, or the inner disk traced by CO obscuration.This dichotomy in turbulent velocities raises the question of what is unique about DM Tau and IM Lup, which in turn is related to the question about what could be driving the turbulence around these sources.E.g., if the turbulence is driven by the magneto-rotational instabilities then the magnetic field and ionization structure may be different between the turbulent and non-turbulent systems, while if turbulence is instead driven by the vertical shear instability, then the density and temperature structure, and the cooling timescale, are important factors.With a sample of six sources, these questions are difficult to fully answer, but we can survey our current knowledge.

Gravito-turbulence
If the disc is sufficiently cold and massive then gravito-turbulence may drive turbulent gas motions.Models of gravito-turbulence predict velocities of 0.2 -0.4c  (Forgan et al. 2012;Shi & Chiang 2014), consistent with the turbulent motion around IM Lup.High resolution dust continuum images of IM Lup reveal spiral arms between 25 and 110 au, which are consistent with gravitational instabilities (Huang et al. 2018;Cadman et al. 2020), although they may alternatively be due to an embedded massive planet (Verrios et al. 2022).
The Toomre Q (=  Ω/Σ) parameter can be used to estimate the susceptibility of the disc to gravitational collapse/instability (Kratter & Lodato 2016).With sufficiently short cooling timescales, the disk will fragment, otherwise gravito-turbulence is more likely.Given that disk fragments have not been observed around IM Lup, we assume the cooling timescale is sufficiently long that the disk would not fragment if Q∼1. Figure 8 shows the Toomre Q value based on the derived disc parameters from the different models.The disc reaches a minimum of Q∼3, which is above the boundary where instabilities begin to set in ( ≲ 1.5, Shi & Chiang 2014).This is consistent with Cleeves et al. (2016) who found that Q reaches a minimum of ∼4 at 70 au.Assuming the temperature profile and surface density shape derived from CO, the disc mass would need to be ≳0.7M⊙ in order for the disc to be gravitationally unstable.
Our ability to rule out gravito-turbulence depends on the assumed radial profile of the disc mass, and the temperature profile.As noted earlier, the midplane temperature may be smaller than derived in our fiducial model, but since the Toomre Q depends on

√
the midplane temperature would need to be an order of magnitude smaller for the disc to be gravitationally unstable.Sierra et al. (2021) find that the radii <20 au are gravitationally unstable, based on a dust mass derived from multi-frequency analysis of radio continuum images, and assuming a gas-to-dust ratio of 100.Lodato et al. ( 2023) use the Table 3. CO velocity profile to derive Q∼1-2 throughout much of the disc, a factor of a few lower than our result despite similar disc masses, due to a much smaller value of   (66 vs 248 au here) and a much steeper temperature profile ( = -0.66 vs. -0.355here).This raises the density in the inner disc, and decreases the temperature in the outer disc, both of which contribute to lowering the Toomre Q throughout the disc.While we have narrowed the region of parameter space (long cooling timescale, to prevent fragmentation, compact disk, steep temperature profile) needed for gravitational instabilities to be present, more direct tracers of the surface density and midplane temperature profile are needed to confidently determine whether or the disc around IM Lup is gravitationally unstable.
In general, disc mass (Table 3) does not stand out as a substantial difference between the turbulent sources (IM Lup: M gas =0.175 M ⊙ , DM Tau: M gas =0.04 M ⊙ ) and the non-turbulent sources (MWC 480: M gas = 0.046 M ⊙ , V4046: M gas = 0.09 M ⊙ , TW Hya: M gas = 0.05 M ⊙ , HD 163296: M gas = 0.09 M ⊙ ).This suggests that gravitoturbulence is not a substantial source of turbulence at these ages even in many of the brightest systems.
This result depends on the assumed disc mass, as well as the surface density and temperature profiles through the disc.While disc masses are difficult to directly estimate from e.g.dust continuum emission (Andrews 2020), a massive disc will leave an imprint on the orbital velocities within the disc, which may be detectable (Verrios et al. 2022).The lower panel of Figure 8 shows the velocity profile, based on our fiducial model, when including the velocity corrections for the pressure gradient and self-gravity for a gravitationally unstable disc (M gas = 0.7 M ⊙ ).The effect of self-gravity is calculated using the method in Veronesi et al. (2021), based on the derivations by Bertin & Lodato (1999) and Lodato (2007).Under self-gravity the orbital frequency becomes: where the last term represents the contribution from disk self-gravity.
Here Φ  is the disk contribution to the gravitational potential, which is calculated via: where  2 = 4 ′ /[( =  ′ ) 2 +  2 and  () and  () are complete elliptic integrals of the first find.If the disc were massive enough to be gravitationally unstable, it would exhibit orbital velocities ∼0.4 km s −1 greater than the Keplerian motion.Pinte et al. (2018) measure sub-Keplerian velocities, deviating from Keplerian motion by ∼0.15 km s −1 , that are inconsistent with the profile expected for a massive disc.The observed velocity profile could be matched with a smaller stellar mass, although this would require M * ∼0.6M ⊙ with M disc ∼0.7 M ⊙ , in contrast to M * =1.1M ⊙ and M disc =0.175 M ⊙ used in our models.We can confidently rule out such a disc to star mass ratio because in that case the signatures of gravitational instability would not be subtle (Kratter & Lodato 2016), and such a low stellar mass has been ruled out by modeling of the disc kinematics (Teague et al. 2021;Lodato et al. 2023;Izquierdo et al. 2023).This indicates that gravito-turbulence is unlikely to be the driver of turbulence around IM Lup.Accounting for disc self-gravity in our modeling of the CO emission, following the velocity prescription of Veronesi et al. (2021), we find that disc self-gravity does not substantially bias our estimate of the turbulence (Table 1).2023) find lower Q values, near instability, with both results having a higher concentration of mass towards the center of the disk.(Bottom:) Orbital velocity in CO emitting region, based on Keplerian motion (grey band, accounting for the varying CO emission heights) and Keplerian motion plus the gas pressure gradient (red band).The blue band includes Keplerian motion, the pressure gradient, and self gravity, with a disc that is massive enough that Q∼1.If the disc were massive enough to be gravitationally unstable, its orbital velocity would show a detectable deviation from Keplerian motion, which is not observed (Pinte et al. 2018, black dots).The drop in the observed velocities below the predictions at  >300 au is possibly a sign of a photoevaporative wind (Haworth et al. 2017).

Magneto-Rotational Instability
In the context of the magneto-rotational instability, the key factors are sufficient ionization and sufficiently strong vertical magnetic fields (Simon et al. 2018).The exact conditions under which MRI is active (e.g.Bai & Stone 2011;Flock et al. 2012) depend on the particulars of the system (e.g., density and temperature) as well as the region of the disk, and the associated non-ideal MHD effects operating (e.g.Ohmic dissipation, Hall effect, Ambipolar diffusion).As an example of the conditions under which MRI could operate, Simon et al. (2018) find that for magnetic field strengths greater than 5-10 G, in the presence of ionization from FUV photons, the MRI would produce detectable (>0.1c  ) turbulence in the ambipolar dominated outer disk.The ionization comes from a combination of FUV, X-ray, and cosmic ray flux with models suggesting that FUV flux plays a larger role in the CO emitting layer (Perez-Becker & Chiang 2011;Simon et al. 2018).Typical FUV luminosities of young stellar objects range from 2×10 −6 − 0.2 L ⊙ (France et al. 2014).There is no clear trend in FUV luminosity among our sample (Table 3), with the nonturbulent HD 163296 having the highest FUV luminosity (Meeus et al. 2012), followed by turbulent IM Lup (Cleeves et al. 2016), non turbulent TW Hya and V4046 Sgr (France et al. 2014), and then turbulent DM Tau having the weakest FUV luminosity (Yang et al. 2012b;France et al. 2014).One caveat in this comparison is that the FUV flux comes from accretion shocks (France et al. 2014), which are highly variable (e.g., the accretion rate onto HD 163296 increased by a factor of 10 in the 15 years between observations Mendigutía et al. (2013)).This, in turn, may lead to time-variable disc ionization (Cleeves et al. 2017).The FUV observations are snapshots in time, which may not reflect the ionizing emission at the time when the turbulence was measured.
Another caveat is that these observations measure the FUV luminosity at the stellar surface, while the FUV flux that reaches the outer disc is more important for influencing turbulence.High-energy photons traveling outwards could be blocked by an inner disc wind (Pascucci et al. 2020), preventing a high stellar FUV flux from ionizing the outer disc.Similarly, a stellar wind can exclude cosmic rays, decreasing the ionization fraction in the disc (Cleeves et al. 2015).Assuming they are unimpeded in their journey to the outer disc, FUV photons penetrate down to Σ = 0.01 − 0.1 g cm −2 , depending on the FUV luminosity, grain abundance, and gas composition, while X-ray photons and cosmic rays can ionize gas closer to the midplane (Perez-Becker & Chiang 2011).The CO emitting heights for our sample are between 0.01 g cm −2 and 0.1 g cm −1 (Figure 9), which, assuming similar levels of irradiating flux, will result in similar ionization structures.As mentioned above, these systems exhibit a range of FUV luminosities and if the depth of ionized region follows the FUV luminosity we would expect ionization to extend deeper into the disks around HD 163296 and IM Lup, and be confined to the largest heights around DM Tau.
Ionized molecules, such as HCO + , N 2 H + and their isotopologues, have been detected around all of our sources (Teague et al. 2015;Cataldi et al. 2021;Aikawa et al. 2021) indicating that at some level the outer discs around these targets contain ionized gas.Whether the disc is sufficiently ionized for MRI requires more detailed analysis of multiple molecular species.Seifert et al. (2021) find high cosmic ray ionization (≳ 10 −17 s −1 ) beyond 80 -100 au in the disc around IM Lup such that the region beyond 100 au is feasibly MRI active, consistent with our detection of non-zero turbulence.Aikawa et al. (2021) find that the ionization rate is higher around IM Lup than around HD 163296 and MWC 480, consistent with the relative turbulence levels between these targets.Cleeves et al. (2015) find a cosmic ray ionization rate ≲ 10 −19 s −1 around TW Hya, again consistent with the weak turbulence around this target (Table 3).
The magnetic field strength remains unknown for planet-forming discs, and could vary enough to explain the difference in turbulence strengths.Simulations of ideal MHD (Hawley et al. 1995;Simon et al. 2009) and ambipolar diffusion in an unstratified disk (Bai & Stone 2011) find turbulent velocity to be correlated with magnetic field strength, but this trend is not as clear in stratified non-ideal MHD simulations (Simon et al. 2018).Broadly speaking, a sufficiently strong magnetic field is needed to activate the MRI, and in this way turbulence is related to magnetic field strength.Recent ALMA observations of line polarization find upper limits of a few mG on the background vertical magnetic field (Vlemmings et al. 2019;Harrison et al. 2021), well above the ∼ 10G regime where MRI operates (Simon et al. 2018).DM Tau and IM Lup, at ∼1 Myr, are younger than the other systems (Table 3) and the magnetic field strength may change over time if the magnetic field is diffused outwards, and/or advected inwards (Bai & Stone 2017;Simon et al. 2018;Cui & Bai 2021).Paleomagnetic measurements within our solar system indicate magnetic fields, of unknown orientation, of ∼0.5 -1 G between 1 and 3 au and > 0.06 G in the outer solar system when our solar system was 2 Myr old, with non-detections with an upper limit of <0.006 G and <0.003 G after about 3.9 Myr and 4.9 Myr for the inner/outer solar system (Weiss et al. 2021).This is consistent with an age dependence for turbulence, assuming all young stars start with similar magnetic field strengths.Recent studies of the vertical extent of mm dust in class 0 and I sources (Villenave et al. 2023;Lin et al. 2023;Guerra-Alvarado et al. 2023) find that these young systems have thicker disks than class II disks, consistent with stronger turbulence at younger ages (but see Lim et al. (2023) for a discussion of how smaller dust scale heights may not necessarily equate to weaker turbulence).
Age is certainly not completely deterministic in setting turbulence.HL Tau is, as a class I object, much younger than DM Tau and IM Lup, but shows weak turbulence (Pinte et al. 2016), as does the class I/II source DG Tau (Ohashi et al. 2023).Similarly, many class II discs with ages similar to DM Tau and IM Lup have gas disc radii consistent with  = 10 −3 − 10 −4 (Najita & Bergin 2018;Trapman et al. 2020;Long et al. 2022).Source-to-source variations in the initial magnetic field strength could result in some systems being MRI active while others are not; observations of magnetic field in protostars are still limited, but they find a variety of strengths (Tsukamoto et al. 2023) and orientations (Hull et al. 2014;Yen et al. 2021).More work is needed to understand the strength of magnetic fields in planetforming discs, and how it evolves with time.

Hydrodynamic Instabilities
If the magnetic field or ionization is sufficiently weak then hydrodynamic instabilities may play a larger role.The Vertical Shear Instability (VSI) is a hydrodynamic instability that is related to the change in orbital motion with height within the disc (Ω ∼ ( 2 +  2 ) −3/2 ), and is expected to play a larger role relative to other hydrodynamic instabilities at the large radii that ALMA observations are sensitive to (Lyra & Umurhan 2021;Lesur et al. 2022).Velocities driven by VSI, which are dominated by vertical motions (e.g., Flores-Rivera et al. 2020), can reach up to ∼100 m s −1 in the upper layers of the outer disc (Flock et al. 2017;Barraza-Alfaro et al. 2021), increasing by a factor of 10 between the midplane and / ∼0.3 (Flock et al. 2017;Pfeil & Klahr 2021).
VSI produces velocities whose amplitude oscillates with a characteristic length scale of ∼0.1R (Flock et al. 2020;Barraza-Alfaro et al. 2021).At small radii, and with low-resolution data, these features are unresolved and would contribute to line broadening in a way that is indistinguishable from isotropic random non-thermal motion, but at larger radii these structures are potentially observable (Barraza-Alfaro et al. 2021) with high-resolution ALMA observations.For the data used here, with a resolution of ∼60 au, radii less than ∼600 au would likely encompass VSI features moving in opposite directions, while higher resolution data from MAPS (Öberg et al. 2021), with a spatial resolution of ∼10au, would only smear together VSI features at radii smaller than ∼100 au.In figure 10 we show the MAPS CO channel maps (Öberg et al. 2021), along with simple models in which a vertical velocity component has been added, that varies as a sinusoidal function of radius with a wavelength of 60 au (=  sin(/(60au))).At a spatial resolution of ∼10 au, the MAPS data can spatially resolve VSI-like features in the outer disc and our simple VSI models produce corrugated structures in the channel maps, similar to those seen in more detailed models (Barraza-Alfaro et al. 2021).Corrugated features are not strongly evident in the MAPS data, but they are difficult to observe at the small velocities (∼50-100 m s −1 ) expected for IM Lup.No corrugated features are evident in the residuals either (Figure 2).
The highly anisotropic nature of VSI (  ≫   ) is very efficient at lofting dust grains in to the disc atmosphere (Flock et al. 2017), consistent with the large dust scale height derived from the IM Lup scattered light image, but may be challenged by the small scale height seen in the millimeter continuum (Franceschi et al. 2022).While larger dust grains are expected to more effectively settle towards the midplane, the VSI may still be effective at lofting up these grains (Dullemond et al. 2022) The dust scale height is directly dependent on the vertical component of the motion, while our observations of the gas motion, under our assumption of isotropic turbulence, would only be sensitive to the line of sight component of this vertical motion.
VSI as the source of the turbulence is challenged by the large disc radius around IM Lup (Table 3), which is consistent with a viscous stress of  ∼ 10 −2 (Trapman et al. 2020;Long et al. 2022), much larger than produced by the VSI ( visc = 1.5 × 10 −4 , Flock et al. 2020).VSI may not exist in isolation, and may operate in concert with MRI and a wind (Cui & Bai 2022).MRI, as well as a wind (Yang & Bai 2021), may contribute to the radial spreading of the disc in the presence of VSI, leading to the observed large disc radii for IM Lup and DM Tau.
In the context of the VSI, turbulence is generated when the cooling timescale is much shorter than the orbital period (Lesur et al. 2022).The cooling timescale in turn depends on the dust population (Umurhan et al. 2017; Lyra & Umurhan 2021), since the large opacity of the dust relative to the gas makes it efficient at cooling the disc.If e.g., dust were present at large radii around DM Tau and IM Lup and not the other sources then this would, at least qualitatively, indicate a shorter cooling timescale around DM Tau and IM Lup.Scattered and polarized light images, sensitive to the small dust grains present in the outer disc, find emission on size scales (Table 3) comparable to the CO emission for all of our sources (Wisniewski et al. 2008;Grady et al. 2010;van Boekel et al. 2017;Avenhaus et al. 2018;Garufi et al. 2024).This indicates that there is no obvious difference in the radial dust distribution between the turbulent and non-turbulent discs, although detailed work is needed to calculate the exact cooling timescales in the outer disc, in order to confirm whether or not the conditions are sufficient for the VSI.
The exact reason for the different turbulence levels among these sources, and the physical origin of the turbulence, is not entirely clear, although there are some hints in this small sample.We can rule out gravito-turbulence in the case of IM Lup, assuming the density and temperature structure derived from the CO observations, as it would impart a detectable signal on the orbital velocities of the CO gas around IM Lup.The disc around IM Lup is sufficiently ionized  2021)) as compared to simple models of VSI-like motion (2nd, 3rd, and 4th rows, Δ  sin(/(60au) )), and isotropic turbulence (bottom row).The VSI-like models produce corrugated features that are visible in the outer disc.These become increasing difficult to detect as the velocity scale approaches that of the observed turbulence around IM Lup (  ∼ 50 − 100m s −1 ).There is no strong evidence for VSI-like behavior in the disc around IM Lup, although it may be below the detection limit of the data.
to be susceptible to the MRI (Seifert et al. 2021), although the unknown magnetic field strength makes it difficult to assess whether or not the MRI is active in this system.The young ages of the turbulent DM Tau and IM Lup discs are consistent with a scenario in which the magnetic field strength evolves with time (Simon et al. 2018), similar to that seen in our solar system (Weiss et al. 2021), although more work is needed to confirm this conclusion.High resolution, deep observations could reveal disc structures driven by the VSI (Barraza-Alfaro et al. 2021), or spatially varying turbulence (e.g., Bosman et al. 2023).The high-resolution MAPS observations have the potential to reveal such structures, as well as to explore e.g.radial variations in turbulence.Such an analysis requires a detailed accounting for variations in the density, temperature, and CO abundance structure beyond the parametric forms employed here.While further work is needed to understand the mechanism driving turbulence, and the demographics of turbulence among a broader sample, the our measurement of turbulence around IM Lup provides valuable clues as to the conditions under which planets can form.ity of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 716155 (SACCRED).This work was also supported by the NKIFIH excellence frant TKP2021-NKTA-64.
This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Figure 1 .
Figure 1.CO(2-1) spectra of the disc around IM Lup (black line) and the median of the PDFs from the fiducial model (red dashed line).The region that is subject to absorption by the molecular cloud, and is excluded from the MCMC process, is marked with a grey band.Despite the limited spectral range, we are able to place strong constraints on the non-thermal linewidth (  turb =0.237 +0.017 −0.012 c  ), with a significantly better fit to the data than with zero turbulence (green dotted line).

Figure 2 .
Figure 2. Channel maps of the CO(2-1) data (top row), the fiducial model (middle row), and residuals (bottom row).In the data and model channel maps the contour indicates the 5 level ( = 9 mJy/beam).In the residual map, positive (black) and negative (red) contours are at levels in multiples of 5.The lower left panel includes a 100 au scale bar and the beam shape.The model reproduces much of the emission, although it does underestimate the emission near the midplane.

Figure 3 .
Figure 3. Temperature (top panel) and height of the emitting region (bottom panel) vs radius for our fiducial model (red band) and as derived directly from the data by Law et al. (2021) (points).In the top panel the dashed line indicates a CO freeze-out temperature of 19 K, while in the bottom panel the dashed line indicates the location of the CO condensation front in our fiducial model.The extended emitting height for CO(2-1) comes from the fact that we see both the near and far side of the disc; the top of this region is seen on the near side of the disc, while the bottom is seen on the far side of the disc.The red Gaussian indicates the beam size in the data we analyze.

Figure 4 .
Figure 4. Top row: Individual channel maps for three different models, with varying levels of turbulence and CO abundance (the white contour indicates the 5 level).There is a degeneracy between turbulence and CO abundance, such that low turbulence can be matched with a higher CO abundance.This can be seen in comparing the morphology of the left and right panels, which is similar, despite the differences in abundance and turbulence.Bottom Row: IM Lup data, followed by maps of optical depth (dotted line corresponds to  = 1) and height of the  = 1 surface (dotted line corresponds to z=0) for the fiducial model ([CO/H 2 ]=5×10 −6 ,  turb = 0.24c  ).CO(2-1) is sensitive to the CO abundance, despite being mostly optically thick, because there are lines of sight that are optically thin (i.e., outside of the region marked by the dotted lines in the central panel) that reach deep into the disc.

Figure 5 .
Figure5.CO(2-1) spectrum (top panels), 13 CO(2-1) spectrum (middle panels) and C 18 O(2-1) spectrum (bottom panels) of the disc around IM Lup (black line) and the median of the PDFs from the multi-line fit (red dashed line).We are able to fit all three spectral lines with non-zero turbulence and with significant CO depletion (left panels), but significantly overestimate the 13 CO and C 18 O flux when assuming no CO depletion (right panels.

Figure 6 .
Figure 6.Temperature (top panel) and height of the emitting region (bottom panel) vs radius for our multi-line model fit.The left panels show CO(2-1) in red, while the right panels show 13 CO(2-1) in blue and C 18 O(2-1) in cyan, with temperatures and emitting heights derived by Law et al. (2021) for CO(2-1) and 13 CO(2-1) shown as circles.Our multi-line fit reproduces the temperature and height of the emitting region for these isotopologues.
CO snowline at the midplane at 15 au.d Within the multi-line fit, the CO abundance is also a free parameter: log([reproduces the CO(2-1) emission, it strongly over-predicts the 13 CO and C 18 O emission.

Figure 7 .
Figure 7. (Top): Height of the scattered light surface, normalized to the radius, for the dust rings identified in the scattered light image from Avenhaus et al. (2018) (black dots), along with models with different dust scale heights.Ratios of dust to gas scale height between 0.7 and 1.0 match the height of the scattered light rings.(Bottom): Constraints on  based on the dust scale height, as a function of the Stokes number.The dark grey band is the constraint on turbulence from the CO observations with CO depletion.The constraint on  from the dust scale height (solid line with arrows) is a lower limit given the asymptotic behavior of  as H  /H  approaches 1, and is consistent with the results fromFranceschi et al. (2022), who simultaneously constrain  and the maximum grain size.These constraints on the dust settling are also consistent with a decrease in the turbulence between the surface layers probed by CO, and the midplane probed by dust settling.

Figure 8 .
Figure 8. (Top:) Toomre Q values, based on disc models derived from the CO observations.The solid curved line represents the fiducial model, while the gray band around this line demonstrates the range of Q values among the different model fits, illustrating the uncertainty in the Q estimate.The radii over which the spiral arms are seen in the dust emission are indicated by the gray shaded region.The horizontal solid line indicates the instability boundary of Q=1.Our derived disc structure lies comfortably above this boundary, indicating that the disc is likely gravitationally stable.Sierra et al. (2021) and Lodato et al. (2023) find lower Q values, near instability, with both results having a higher concentration of mass towards the center of the disk.(Bottom:) Orbital velocity in CO emitting region, based on Keplerian motion (grey band, accounting for the varying CO emission heights) and Keplerian motion plus the gas pressure gradient (red band).The blue band includes Keplerian motion, the pressure gradient, and self gravity, with a disc that is massive enough that Q∼1.If the disc were massive enough to be gravitationally unstable, its orbital velocity would show a detectable deviation from Keplerian motion, which is not observed(Pinte et al. 2018, black dots).The drop in the observed velocities below the predictions at  >300 au is possibly a sign of a photoevaporative wind(Haworth et al. 2017).

Figure 9 .
Figure 9. Grey bands represent the heights between which 50% of the flux for CO(2-1) arises at a given radius (as described in Section 3), for the different systems (the blue band in each panel is the CO(2-1) emitting region around IM Lup).The dashed lines indicate the pressure scale height, , and three times the pressure scale height.The dotted lines indicate surface density boundaries of Σ = 0.01 g cm −2 and 0.1 g cm −1 .While the emitting regions around DM Tau and IM Lup are systematically higher relative to the sources with weak turbulence, they are not systematically higher relative to the pressure scale height, or the column density within the disc.

Figure 10 .
Figure 10.Observed CO(2-1) channel maps (top row, Öberg et al. (2021)) as compared to simple models of VSI-like motion (2nd, 3rd, and 4th rows, Δ  sin(/(60au) )), and isotropic turbulence (bottom row).The VSI-like models produce corrugated features that are visible in the outer disc.These become increasing difficult to detect as the velocity scale approaches that of the observed turbulence around IM Lup (  ∼ 50 − 100m s −1 ).There is no strong evidence for VSI-like behavior in the disc around IM Lup, although it may be below the detection limit of the data.

Figure A1 .
Figure A1.Channel maps of the data, the fiducial model, and the residuals.In the data and model channel maps the contour indicates the 5 level ( = 9 mJy/beam).In the residual maps, positive (black) and negative (red dashed) contours are at levels in multiples of 5.Lower left panels include a 100 au scale bar and the beam shape.Channels between 4 and 6 km s −1 have been removed due to cloud contamination.

Figure B1 .
Figure B1.Channel maps of the data, the fiducial model, and the residuals.In the data and model channel maps the contour indicates the 5 level ( = 9 mJy/beam).In the residual maps, positive (black) and negative (red dashed) contours are at levels in multiples of 5.Lower left panels include a 100 au scale

Table 2 .
Powell et al. (2022)l.(2023)inding of strong non-thermal motion around IM Lup byPaneque-Carreño et al. (2023).In the context of an  disc model, with  ∼ ( turb /  ) 2 , and under the assumption that  turb /  is constant throughout the disc, our result corresponds to  = 0.03−0.08aroundIMLup,  = 0.06−0.10aroundDMTau, and  ≲0.01 around the other sources.Powell et al. (2022)examine the CO abundance in the context of gas diffusion and dust grain surface Bosman et al. (2021)igh gas diffusion, consistent with consistent with  = 0.003 − 0.015, around IM Lup and DM Tau and low diffusion, consistent with  ∼ 10 −4 , around HD 163296 and TW Hya.Bosman et al. (2021)find evidence for high turbulence ( > 4 × 10 −3 ) in the inner 20 au of the disk around IM Lup, based on the dust obscuration of CO emission from the inner disk.These constraints, listed in Table2, are in line with our results, with the caveat that  values