High turbulence in the IM Lup protoplanetary disk Direct observational constraints from CN and C 2 H emission

,


Introduction
The study of planet formation requires high sensitivity observations that can spatially and spectrally resolve the motions of the material in protoplanetary disks around young stars.In this context, the Atacama Large Millimeter/submillimeter Array (ALMA) is a unique facility, able to characterize these objects and their evolution mechanisms, which depend on the transportation of angular momentum (Lynden-Bell & Pringle 1974).Turbulence is the mechanism proposed to shape the surface density of disks, by transporting angular momentum through velocity fluctuations and providing an effective viscosity at microscopic scales (Shakura & Sunyaev 1973).There are many hydrodynamic processes that could be the origin of turbulent motions (see Lesur et al. 2023, and references therein); however, direct measurements of turbulence in protoplanetary disks have been hard to obtain, with only a few upper limits (e.g., Hughes et al. 2009;Flaherty et al. 2015;Teague et al. 2016) and one candidate with measured turbulence, DM Tau.In this system, high spatial and spectral resolution ALMA observations of CS allowed Guilloteau et al. (2012) to initially measure turbulence values of 0.3-0.5cs (see also the initial studies by Guilloteau &Dutrey 1998 andDartois et al. 2003).A later study by Flaherty et al. (2020) obtained lower, but still significant, turbulence values of 0.25-0.33cs through modeling of CO emission in DM Tau.
Current methods used to constrain turbulence mostly depend on comparing observations to complex parametric and radiative transfer models, which depend on fundamental physical disk properties such as surface density and temperature structure (e.g., Hughes et al. 2009;Guilloteau et al. 2012;Flaherty et al. 2015Flaherty et al. , 2017Flaherty et al. , 2020)).While these methods constrain their best-fit parameters with state-of-the-art analysis, there are multiple effects that are hard to account for since precisely resolving the gas surface density profile requires knowledge of the dust grain distribution and excitation conditions of the observed molecule (Zhang et al. 2021;Miotello et al. 2023).Furthermore, the presence of substructure in the disk may or may not cause gaps and rings in the surface density, depending on their origin, which is also not trivial to determine (Rosotti et al. 2021;Zhang et al. 2021;Teague et al. 2018a).
Due to their complexity, previous uses of these models to extract disk turbulence have assumed a single turbulence value for the whole protoplanetary disk, neglecting the possibility of a radial or vertical dependence.However, turbulence could have vertical variations, depending on the launching mechanism (Forgan et al. 2012;Simon et al. 2015;Shi & Chiang 2014), and while various molecular tracers are used to probe the vertical structure in searches for vertical turbulence variations (Flaherty et al. 2017), locating the exact position of the molecular emission in the disk is not straightforward in the modeling process and in some cases presents strong discrepancies with observa-tional constraints.A successful case is the best-fit model to extract turbulence in DM Tau, where the CO (J = 2 − 1) emission is expected to arise from a layer of z/r ∼0.3-0.4 (see Figure 2 in Flaherty et al. 2020), in good agreement with observational constraints from Law et al. (2023).However, analysis of HD 163296, where no significant turbulence was detected, predicts that CO (J = 2 − 1) emission is coming from z/r <0.1 (see Figure 16 in Flaherty et al. 2017), which is many times lower than the measured emission surface of CO (J = 2 − 1), located at z/r ∼0.3 (Law et al. 2021b;Paneque-Carreño et al. 2023).These differences may be due to the limitations of the method, a lack of spatial resolution, or errors in the assumed density and temperature structure.
To alleviate these issues, this work focuses on directly tracing the vertical location of optically thin CN and C 2 H transitions, measuring the radial profile of line broadening due to nonthermal motions, and relating this to vertically and radially resolved turbulence.CN and C 2 H are expected to originate from chemical processes requiring energetic UV-photon radiation (Visser et al. 2018;Cazzoletti et al. 2018;Bergin et al. 2016;Teague & Loomis 2020), and therefore their emission should arise from a sweet spot located at high elevation from the midplane (Schreyer et al. 2008;Cazzoletti et al. 2018), in a zone known as the photon-dominated region (Aikawa et al. 2002).As the emission should come from a vertically thin slab, it is possible to accurately determine the location of the emission surfaces of CN and C 2 H despite the emission being optically thin (see Paneque-Carreño et al. 2022).In parallel, optically thick CO isotopologs can be used as tracers of the disk structure, with the resolved vertical emission location (Law et al. 2021b;Paneque-Carreño et al. 2022) allowing us to provide a global three-dimensional characterization of the temperature structure, molecular emission, and turbulence in the upper disk layers.This is ultimately achieved by measuring the emission linewidth and separating the thermal from the nonthermal component.As we assume optically thin molecular emission, the nonthermal component is assumed to be mainly turbulent (Hacar et al. 2016).A key tool in our analysis is the use of DISCMINER (Izquierdo et al. 2021(Izquierdo et al. , 2023)), a disk analysis package that fits the spectra per pixel of our data to precisely extract and parametrize the radial linewidth profile of the data, accounting for the effects of beam convolution and upper and lower emission surfaces.
Our work focuses on the protoplanetary disk around the young K5 star IM Lup, located at a distance of 158 pc (Gaia Collaboration et al. 2023).This system has been extensively studied using multiple tracers.Infrared scattered light shows a large, vertically extended disk, with radial rings and gaps (Avenhaus et al. 2018).Millimeter continuum traces tightly wound spiral arms and an outer ring (Cleeves et al. 2016;Huang et al. 2018a,b), and gas line emission shows a rich molecular reservoir (Pinte et al. 2018a;Öberg et al. 2021).Recent observational and theoretical works have proposed that IM Lup has an inner turbulent disk (Bosman et al. 2023) and possible turbulent velocities of ∼100 m s −1 in the outer disk (Cleeves et al. 2016).Due to its inclination and vertical extent, it is possible to trace the vertical disk structure and molecular layering of various lines; however, previous studies have been unable to extract CN or C 2 H vertical profiles due the low signal-to-noise ratio (S/N) of the lower rotational transitions (Paneque-Carreño et al. 2023).New, higher-energy ALMA Band 7 transitions have sufficient S/N for an adequate characterization of their radial, azimuthal, and vertical features.
This paper is structured as follows.Section 2 presents the data and explains the calibration procedures.Section 3 studies the emission morphology, characterizing the vertical extent of the disk and the brightness temperature of CN and C 2 H and comparing them to those of CO isotopologs.Section 4 shows the parametric model we used to extract the linewidth from CN, and Section 5 presents our results on the radial turbulence profile as traced by CN.Finally, Section 6 discusses the findings of this work in the context of previous constrains, and Section 7 summarizes our main conclusions.

Observations
This work uses ALMA Band 7 C 2 H (N = 4 − 3) and CN (N = 3 − 2) observations of IM Lup (Project code: #2019.1.01357.S, PI: R.Teague) taken during the second semester of 2019 and first semester of 2021.After initial calibration through the ALMA pipeline (Hunter et al. 2023), the data were flux-calibrated and self-calibrated in phase and amplitude using only the linefree spectral windows and channels with the tools of Common Astronomy Software Applications (CASA;version 5.4.0, Mc-Mullin et al. 2007) as described in the following.
The observations from each execution block (EB; seven in total) were aligned to a common phase center, considering the peak of the continuum emission of each EB as fitted by CASA task imfit.After alignment, an amplitude scaling was applied to correct for the relative flux scale variations between EBs, which were <10%, following the procedures of ALMA large program DSHARP (Andrews et al. 2018).Self-calibration was then applied to the short-baseline observations (two EBs, with a longest baseline of 0.3 km).Five rounds of phase and one round of amplitude calibration were performed, until there was no further improvement of the S/N, starting with a solution interval of 3360s and taking half the time step in each subsequent round.The S/N improvement in the short-baseline data was 726%, from 190 to 1380.Finally, self-calibration was performed on the joint dataset of the previously self-calibrated short baselines, concatenated to the flux-calibrated long baselines (five EBs, with a maximum baseline of 1.5km).Three rounds of phase and one round of amplitude calibration were performed, the initial time solution interval was 4800s and the final S/N increase was 428%, from 301 to 1289.
Each of the calibration solutions from the continuum analysis was then applied to the whole dataset, to correct spectral windows and channels with line emission.Using the CASA task uvcontsub, we produced continuum-subtracted measurement sets for each of our target emission lines.Each C 2 H and CN hyperfine group was imaged using tclean with Briggs weighting and Keplerian masking 1 .Multiple Keplerian masks, one for each of the expected hyperfine components, were used in the cleaning.A stellar mass of 1.1M ⊙ (Teague et al. 2021) and a constant z/r value of 0.2 and 0.3 was used for C 2 H and CN emission, respectively.To compromise between optimal angular resolution and high S/N, C 2 H emission was imaged using a robust parameter of 2.0, resulting in a 0.27 ′′ × 0.25 ′′ beam (PA: 64.83 • ).CN is imaged using a 0.5 robust parameter with a final beam of 0.22 ′′ × 0.19 ′′ (PA: 57.56 • ).JvM correction (Czekala et al. 2021) was applied to our final images to account for the noise caused by having non-Gaussian beams.The epsilon correction factors are 0.24 and 0.52 or C 2 H and CN, respectively (for details on this procedure and the epsilon factor, refer to Czekala et al. 2021).The images used in the analysis where obtained using a channel spacing of 200 m s −1 for C 2 H and 100 m s −1 for CN, to optimize spectral resolution and S/N.The final rms of the images for each  molecule are 0.6 mJy/beam and 1 mJy/beam for C 2 H and CN, respectively.
Figure 1 shows the integrated intensity maps (moment 0) and stacked spectra for the data imaged at the native spectral resolution, 50 m s −1 .The spectra were computed using the GoFish (Teague 2019) function, integrated_spectrum, over the whole disk.Each of the detected hyperfine components are marked with an orange vertical line and their properties, together with the integrated fluxes indicated in Table 1.A clear ring is observed in both tracers for all hyperfine groups, a detailed analysis of this feature can be found in Appendix A.

Vertical profile
From the channel maps it is possible to visually identify the CN and C 2 H emission originating from the upper and lower sides of the disk.Tracing the vertical location of each molecule allows us to characterize the disk structure and physical conditions of the emitting gas in defined regions.To extract the vertical profiles, we used the code ALFAHOR (Paneque-Carreño et al. 2023), which traces the maxima of emission within a given masked region.The masks are visually determined and encircle the far and near sides of the disk upper surface (for details on the method and terminology, see Pinte et al. 2018a;Paneque-Carreño et al. 2022, 2023).Only channels where both near and far sides of the upper surface can be confidently identified are used for this analysis.Figures B.1 and B.2 present the masks and location of the emission maxima for each of the occupied transitions.
A key assumption of this method is that the intensity maxima trace a single isovelocity along each channel.Due to the strong blending of hyperfine components in the CN emission and the low S/N of some detected transitions, we only used the unblended and strongest component CN N = 3−2J = 5/2−3/2F = 7/2 − 5/2.C 2 H also has blending in both transitions; however, it has a sufficient separation that we can identify both components in each group (see Fig. 2).However, we extracted the surface only using the J = 9/2 − 7/2 hyperfine group, as it has a higher S/N.We expect that these transitions are representative of the spatial location of all the detected transitions of the same molecule, as they all have similar upper state energies (see Table 1) and thus, originate from the same layer.
Table 1.Molecular line data from the LAMBDA database (Schöier et al. 2005).

Species
J Transition  The resulting vertical distribution is parametrized with an exponentially tapered power law to capture the inner flared surface, plateau and turnover regions at larger radii (Teague et al. 2019) following The best-fit values for z 0 , ϕ, ψ, and r taper were found using a Markov chain Monte Carlo sampler as implemented by emcee (Foreman-Mackey et al. 2019) and shown in Table 2. Figure 3 shows the data considering all of the extracted maxima for CN and C 2 H, together with the best-fit exponentially tapered powerlaw model.
sion is located slightly below, tracing the same region as C 18 O (J = 2 − 1).We note that C 2 H presents a large scatter, likely due to the lower S/N of the observations; however, our result on the average vertical profile (panel C of Fig. 3) is consistent with z/r ∼0.2 expected from the visual channel map analysis (Fig. 2).

Brightness temperature
Considering the best-fit vertical surfaces, the radial peak brightness temperature profiles were extracted from the peak intensity maps (moment 8).The conversion from peak intensity (I) to brightness temperature (T b ) was done using Planck's law,   12 CO J=2-1 13 CO J=2-1 C 2 H N=4-3, J=9/2-7/2 3 1.5 0 -1.5 -3 ["] where h is the Planck constant, k the Boltzmann constant, c the speed of light, and ν the frequency of the emission.To compare the brightness temperatures of CN and C 2 H to those of the CO isotopologs previously studied in Paneque-Carreño et al. ( 2023), all brightness temperature profiles were calculated using cubes with an equal spectral resolution of 0.2 km s −1 to match the resolution of the Molecules with ALMA at Planet-forming Scales (MAPS) Program data (Öberg et al. 2021).
The peak intensity maps and brightness temperature profiles, compared to those of CO isotopologs (Paneque-Carreño et al. 2023), are presented in Figure 4.The radially extended bright feature along the major axis seen in the peak intensity maps for all tracers is a radiative transfer effect related to the higher density of isovelocity contours and does not indicate a physical feature.Based on our previous analysis we know that CN and C 2 H emission is spatially colocated with 13 CO; however, both transitions show lower brightness temperatures than 13 CO for all hyperfine groups.The drop in temperature at the innermost radii is dominated by beam dilution as the emitting area becomes comparable to or smaller than the spatial resolution; however, this drop extends out to ∼150 au, which is beyond the beam size.Dust absorption is likely a big contributor to the low temperatures, as several works have shown that the inner disk of IM Lup is optically thick (Cleeves et al. 2016;Sierra et al. 2021).
["] Residuals The brightness temperature may be a tracer of the excitation temperature of the region, but only in the cases where the emission is optically thick (τ >>1), as it is proportional to T ex (1−e −τ ).CO is expected to be abundant in IM Lup, in particular 12 CO and 13 CO should be optically thick (Cleeves et al. 2016) and therefore their brightness temperature may be taken as a tracer of the kinetic temperature in the region.On the other hand, CN and C 2 H are expected to have lower column densities and be optically thin tracers (Cazzoletti et al. 2018;Bergner et al. 2021;Bergin et al. 2016), which is corroborated by the lower brightness temperatures measured in this work.We note that the temperature difference between both CN hyperfine groups is due to the strong blending in the J = 7/2 − 5/2 transition (see Fig. 1), which means that the peak flux for a given pixel is the sum of multiple hyperfine components.This is the same situation that affects the integrated flux values shown in Table 1.

Parametric emission model: DISCMINER fit
Knowing that the molecular emission is optically thin enables us to directly determine turbulence from the nonthermal contributions to the linewidth (Hacar et al. 2016).This is not possible in optically thick tracers (such as CO isotopologs) due to the additional opacity broadening that affects the linewidth in those cases, which is why in these cases radiative transfer models are needed to account for the opacity effect.Extracting the linewidth is not straightforward, it requires a parametric model that accounts for the vertically elevated emission, the upper and lower sides as well as radial intensity and velocity gradients.To this end, we used DISCMINER (Izquierdo et al. 2021), which allowed us to find the best-fit parametrization of our emission.

Model setup
DISCMINER determines the best orientation, velocity, surface and line profile parameters for the observed emission by creating a model intensity, I m as a function of the disk cylindrical coordinates (R, z) and the velocity.The emission kernel has two components, one for the upper and another for the lower surface, which are well separated in the case of our studied emission.
Contrary to previous works with DISCMINER, which focus on CO isotopolog emission (Izquierdo et al. 2021(Izquierdo et al. , 2022(Izquierdo et al. , 2023)), this work traces optically thin molecules that will be well reproduced by a Gaussian line profile (Hacar et al. 2016).Therefore, we used a Gaussian kernel for each surface such that with I p being the peak intensity and ∆V the linewidth at a given radial and vertical distance from the star.Since we assumed fully optically thin emission, the final model intensity is given by the direct sum of the upper and lower disk surface intensity contributions and compared against the emission spectra extracted from each pixel of our observational data cube (for specific details, see Izquierdo et al. 2021Izquierdo et al. , 2023)).Table 3 indicates the equations used to model each of the disk attributes.Each emission surface is modeled independently as an exponentially tapered power law, the intensity is considered to be a power law with a Gaussian component, to account for the ring in the integrated emission maps (see Figure 1).The linewidth is assumed to follow a simple radial and vertical power law (as indicated in Table 3).DISCMINER is not built to fit emission with multiple hyperfine components; therefore, we created a cube that only covers the emission of the strongest unblended hyperfine component of CN 3-2, J = 5/2 -3/2 F= 7/2-5/2 (see Figure 1).It is not possible to analyze the other CN hyperfine components, due to the strong blending of the high S/N lines, and the low signal of the other detections.C 2 H is also excluded from this analysis due to its low S/N and the close blending of the lines, which does not allow us to isolate them as we can for CN (see Figure 3 for a comparison on how close the CN hyperfine components are compared to C 2 H in the channel map emission).Future work incorporating hyperfine components to the DISCMINER kernel could allow us to optimize our model using the whole dataset, but it is beyond the scope of this study.Comparison between various radial profiles between the data (magenta curves) and the DISCMINER model (gray curves).All panels except the height profile indicate the profile dispersion as a shaded region.Velocity and linewidth panels have a noticeable difference located at ∼220 au, a blue residual curve characterizes this feature in the upper right corner of each panel.The linewidth panel shows the measured linewidth extracted from the data and the beam-convolved model by fitting the spectra.Additionally, the dashed black line indicates the intrinsic parametric form of the linewidth, as fitted by DISCMINER.The beam-convolved model will always have a broader linewidth than the parametric form and is used to directly compare our model to the data, which has a finite sptial resolution.However, the intrinsic linewidth is what we used in our calculations, as this is expected to be the actual line broadening due to thermal and nonthermal motions.
Notes.G is the gravitational constant and D 0 = 100 au is a normalization factor.R is the cylindrical radius and z is the height from the midplane, always positive, r is the spherical radius.The other variables are free parameters, even if named alike they are modeled independently for each attribute.

Comparison to CN observations
The best-fit model parameters are found in Appendix C. Figure 5 shows the resulting model channel maps, convolved with a Gaussian beam to be at the same spatial resolution as the observations.The 5σ residuals are presented in the bottom row, no coherent or systematic pattern is identified.Figure 6 shows further comparison between the model and data.There is good agreement between the derived vertical and azimuthally averaged intensity profiles and the DISCMINER parametric forms.
The rotation curve shows a ∼200m s −1 residual between 200 and 250 au, which coincides with a ∼50m s −1 negative residual in the linewidth.This velocity variation is also observed in 12 CO (see Figure B.18 in Izquierdo et al. 2023) and the location of the feature is close to a dust gap at 209 au and a dust ring at 220 au (Huang et al. 2018a), which may be related to these variations.
We note that the linewidth values extracted from the channel map images (red and gray curves, with uncertainty regions in lower right panel) are higher than the best-fit DISCMINER parametric form (traced by dashed black line).This is because the profiles extracted from the image correspond to direct measurements of the linewidth from the spectra extracted at each pixel in the image plane.These measurements will always be broader (have higher values) than the underlying molecular emission linewidth due to the linewidth broadening caused by the finite spatial resolution of the observations (usually called beam broadening and discussed also in Teague & Loomis 2020).When we refer to the DISCMINER parametric form, we are describing what is expected to be the native molecular linewidth that, convolved with a Gaussian beam equal to the spatial resolution of the ALMA images (∼0.22 ′′ ), will match the broader values extracted directly from the data.

Turbulence as traced by CN
The measured linewidth of an optically thin tracer will be dominated by doppler broadening caused by thermal and nonthermal motions (Hacar et al. 2016).As we used a Gaussian kernel (see Eq. 3), our measured linewidth ∆V corresponds to the velocity dispersion such that the full width at half maximum (FWHM) of the kernel is 2 log2 ∆V.The linewidth can then be written as The thermal component is well characterized and proportional to the ratio between the kinetic temperature and the molecular mass of the observed tracer (Myers et al. 1991).The nonthermal component can originate from non-Keplerian motion such as photo-evaporative winds, planet-disk interactions, or turbulence (e.g., Fuller & Myers 1992;Hacar et al. 2016;Izquierdo et al. 2023).For this analysis we considered the case where nonthermal motion is fully due to turbulence; the implications and possibilities of other origins are reviewed in the discussion section.Therefore, the two components of the measured linewidth can be decomposed as Here k is the Boltzmann constant and m X the molecular mass, in this case, for CN, m X = 26.0(Klisch et al. 1995;Schöier et al. 2005).The nonthermal broadening is assumed to be completely due to turbulence, turb .Turbulence can then be expressed in terms of the Mach number, M, such that where c s is the local sound speed in the studied layer c 2 s = kT kin µm H (Balbus & Hawley 1998;Flaherty et al. 2015).We generally refer to turbulence in terms of the Mach number so it may be directly compared to previous works and estimations (Flaherty et al. 2015(Flaherty et al. , 2017(Flaherty et al. , 2020;;Teague et al. 2016Teague et al. , 2018b)).By spatially localizing the CN emission vertically and radially close to the optically thick CO isotopes, we have a tool to directly access the kinetic temperature without need to use parametric models or additional molecular transitions.Using the linewidth parametrization found by DISCMINER we can immediately calculate the turbulence radial profile, presented in Figure 7.Our results indicate turbulent motions of M =0.4-0.6 at radii beyond ∼150 au.The shaded regions of the plot represent the variation of our results considering a 50 m s −1 error in the linewidth measurement.While the statistical error of our fit is <5 m s −1 , we took the native resolution of our data as the largest possible deviation, in case effects such as Hanning smoothing are significantly affecting our linewidth extraction (Cortes et al. 2023).Under the α-prescription for turbulent viscosity (ν = αc s H, Shakura & Sunyaev 1973) it is possible to relate the turbulence Mach number M to α following, α ∼ M 2 (Rosotti 2023).Our measurements indicate α ∼ 10 −1 .
The total measured linewidth (∆V) and turbulence values increase steeply in the inner ∼150 au.While IM Lup has been proposed to be extremely turbulent in the inner disk (Cleeves et al. 2016;Bosman et al. 2023) it is unclear if our analysis is sensitive to this effect, or if the high turbulence should extend beyond the innermost 20 au, given the lower estimated scale height for millimeter dust in this system (Bosman et al. 2023).Linewidth measurements are strongly affected by overlapping velocities in the inner disk (Teague & Loomis 2020;Izquierdo et al. 2022); therefore, we only considered the values beyond 150 au as accurate measurements and do not refer to the high values inward of 150 au.
The average total linewidth measured outward of 150 au is ∼221 m s −1 .Thermal broadening only accounts for 105-132 m s −1 , using the 13 CO and 12 CO brightness temperatures, respectively.Therefore, the nonthermal broadening ranges between 177 and 194 m s −1 and would require temperatures ∼40 K higher than measured to account for it as thermal motion.Figures 8 and 9 illustrate the need of a high nonthermal component in order to explain the measured total linewidth (see Fig. 8) and the data spectra (see Fig. 9).6).It is compared to the thermal component obtained using the 13 CO (blue) or 12 CO (orange) brightness temperatures.

Vertical location of CN and C 2 H
Our work traces the CN and C 2 H emission originating from a vertical region similar to 13 CO, at z/r =0.2-0.3 constrained radially between 50 and 350 au.These results are in agreement with chemical models of IM Lup that locate the tracer at z/r =0.2 (Cleeves et al. 2018).However, while in agreement with theoretical models, this is the first direct observational evidence of C 2 H coming from vertically elevated layer in the outer disk, as previous works studying other protoplanetary disks had only traced C 2 H inward of 150 au and closer to the midplane (Law et al. 2021b;Paneque-Carreño et al. 2023).
Both CN and C 2 H are formed via UV-driven chemistry (Nagy et al. 2015;Bergin et al. 2016;Visser et al. 2018;Cazzoletti et al. 2018;Miotello et al. 2019;Guzmán et al. 2021;Bergner et al. 2021) and tracing their vertical distribution in various transitions will allow us to further understand the range of penetration that these energetic photons have toward the midplane, or the mixing processes that may result in some transitions tracing different vertical zones.The distribution of the small grains may have a particularly strong effect on how the UV field affects the disk (Bergin et al. 2016;Bosman et al. 2021).Combining scattered light observations with high resolution ALMA data of UV-sensitive molecular lines will give insight into these effects.Overall, future work should focus on characterizing differences in the vertical location of various transitions for a single system, particularly of CN and C 2 H. IM Lup is the perfect laboratory for these studies.

Origin of high turbulence in IM Lup
The physical structure of IM Lup has been thoroughly studied with models (e.g., Cleeves et al. 2016) that propose the existence of a turbulent inner disk (<20 au Cleeves et al. 2018;Bosman et al. 2023) with high dust optical depth (Huang et al. 2018a;Sierra et al. 2021) that may be causing the cavity seen in molecular emission.Additionally, a series of instabilities may be ongoing in IM Lup.Bosman et al. (2023) proposed that the turbulent inner disk could be related to vertical shear instability (Urpin & Brandenburg 1998).The high disk-to-star mass ratio measured in multiple works (0.1-0.2, Lodato et al. 2023;Cleeves et al. 2018;Zhang et al. 2021;Sierra et al. 2021) suggests gravitational instabilities (Kratter & Lodato 2016) throughout some portion of the system.Finally, analysis of the ionization structure suggests a radial gradient in the cosmic ray ionization rate and adequate conditions for magneto-rotational instabilities (MRIs; Balbus & Hawley 1991) in disk regions beyond the millimeter continuum disk edge (∼313 au) or above 20 au or z/r ∼ 0.25 from the midplane (Seifert et al. 2021).
Indeed, high turbulence in this system is not surprising.Any of these instabilities, or their combined effect, could be causing turbulent motions.There have already been suggestions of nonthermal broadening from low-resolution CO data, as thermochemical models are able to explain the CO morphology using 100 m s −1 turbulent velocities (Cleeves et al. 2016) a value within the order of our derived range of 177-194 m s −1 .In contrast with this, low turbulence values (α ∼ 10 −3 ) are derived from simultaneous modeling of millimeter (ALMA) and micrometer (SPHERE) dust observations (Franceschi et al. 2023;Jiang et al. 2023).This low turbulence is mainly driven by the millimeter dust distribution, which traces the cold midplane.Through our analysis, we constrain high turbulence values (α ∼ 10 −1 ) traced by CN at R >150 au and in the vertically elevated disk layers (above the midplane) at z/r = 0.25 − 0.3, which would suggest a vertical gradient in the turbulence, with higher values in the upper disk and lower values toward the midplane.Theoretical MRI studies have shown that high α values may correlate with higher accretion rates (Simon et al. 2015).IM Lup has a measured disk accretion rate of 10 −8 M ⊙ yr −1 (Alcalá et al. 2017), which is typical of lower α.This is reminiscent of DM Tau, where the system showed measurable turbulence, but a low accretion rate (Flaherty et al. 2020).It is important to note that the measured accretion rate is taken from the instantaneous accretion rate onto the star; however, the α values represent disk motions throughout the outer, and in the case of IM Lup the elevated, regions of the disk.Future work should look into this with dedicated models.
A vertical gradient in the turbulence, with high values approaching the sound speed above 2H, with H being the disk hydrodynamical scale height, is a key prediction of MRIs (Simon et al. 2015;Flock et al. 2015).The scale height in IM Lup is expected to be H ∼ 0.1 (Zhang et al. 2021;Paneque-Carreño et al. 2023); therefore, CN emission is originating above 2H.To fully confirm a vertical gradient in turbulence, studying molecules in the intermediate layers is necessary, such as HCO + , H 2 CO, and HCN, which have their vertical heights characterized and emit from in-between the midplane and 13 CO (Paneque-Carreño et al. 2023).Unfortunately, the sensitivity and spectral resolution of the available data for these tracers (≤200m s −1 , Öberg et al. 2021) is not sufficient to detect nonthermal broadening of ≤90m s −1 , such as detected for CN.Follow up observations at higher spectral resolution would allow us to replicate this analysis and confirm the existence of a vertical gradient.
As a first attempt to compare our results with other less abundant tracers, we repeated our analysis on C 18 O J = 2 − 1 data from the MAPS program (Öberg et al. 2021).Figure 3 shows that C 18 O J = 2 − 1 comes from an elevated layer, comparable to that of CN, but slightly below.This emission is expected to come close to the CO freeze-out region and previous works have assumed an excitation temperature of 21 K for the C 18 O gas (Pinte et al. 2018a;Paneque-Carreño et al. 2022).Assuming that the nonthermal contribution is fully turbulent, the extracted turbulent dispersion profile of C 18 O are slightly higher, but comparable within the uncertainty to the results obtained from CN assuming the temperature profile of 13 CO.We note that while the brightness temperature of C 18 O is below the assumed excitation temperature (see Figure 4), the column density of CO in IM Lup (Zhang et al. 2021) and estimates from Elias 2-27 of C 18 O under similar conditions (Paneque-Carreño et al. 2022), show that this isotopolog may have optical depths of 1-2.In this case, opacity broadening could account for up to 25% of the measured linewidth value (Hacar et al. 2016), biasing our result to higher turbulence profiles.Figure 10 shows the comparison between the CN results and the turbulent dispersion profiles of C 18 O assuming a fully turbulent nonthermal contribution in the same way as our CN analysis (solid line) and considering an opacity broadening of 25% to the measured linewidth (dotted line).While the analysis of C 18 O considering opacity broadening gives a turbulent dispersion profile below that of CN (using the 13 CO temperature), all profiles are comparable within their uncertainty.This is likely due to the similarity in the emitting region, which overlap within the dispersion (for details, see Paneque-Carreño et al. 2023) and highlights the need for a follow up study using other optically thin tracers that clearly probe different innermost disk layers.

Possible overestimation of the nonthermal component
We considered the possibility that we overestimated the turbulence in our analysis, for example by underestimating the temperatures in the CN emitting region or not taking into account alternative origins of the nonthermal broadening.Regarding the temperature profiles, we are confident that CO traces optically thick gas, even toward the outer disk regions.We confirmed this by comparing our temperature profiles to the gas temperature from dedicated thermochemical models of IM Lup (see Figs. 8 and 9 of Cleeves et al. 2016), which coincide with the values used in this work; therefore, we have not underestimated the temperature.Determining other contributors to nonthermal broadening is harder: it could be that CN emission is not as optically thin as expected, that it is vertically extended and the velocity dispersion along the line of sight is considerable, or that some strong non-Keplerian motion is affecting the system.Below, we analyze each of these scenarios.At optical depths of τ ∼0.1 the opacity contribution to the linewidth is negligible, but at slightly higher values, even while τ <1 as expected for CN in IM Lup (Bergner et al. 2021), opacity broadening could contribute up to ∼15% of the measured linewidth (Hacar et al. 2016).For our measured linewidth of ∼221 m s −1 outside of 150 au we considered a systematic error of 50 m s −1 to account for Hanning smoothing effects, which corresponds to a 23% error in the linewidth measurement and also accounts for the possibility of having CN emission with an optical depth up to τ ∼2.Our results would still be consistent with high turbulence, even if the emission were marginally optically thick and opacity was contributing to the nonthermal broadening.
CN emission is expected to arise from a geometrically thin layer (Cazzoletti et al. 2018;Paneque-Carreño et al. 2022); however, at the peak of emission this layer may extend vertically up to ∼40 au (see Figure 3 from Cazzoletti et al. 2018).It is not possible to estimate the vertical width of the CN emission with our current analysis, but we can test the impact on our measured linewidth caused by Keplerian shear along the line of sight for IM Lup.By considering the inclination of the system and taking the best-fit CN emission layer as the position of the lowest parcel of material, we measured the linewidth broadening caused by having different vertical widths.The detail of this analysis can be found in Appendix D. Figure 11 comprises the results, showing that at radial distances beyond 150 au the variation in the linewidth caused by vertical widths of up to 50 au is always ≤25 m s −1 .The width would have to be broader than ∼60 au to cause significant variations out to further radii, which is not expected for CN (Cazzoletti et al. 2018).This indicates that even in the case of CN being vertically extended, our results beyond 150 au will not vary beyond our assumed uncertainty of 50 m s −1 (Fig. 7).There are no immediate signs of non-Keplerian motion (as seen in other systems, for example HD 163296; Pinte et al. 2018b;Teague et al. 2018a;Izquierdo et al. 2022) beyond the strong residual at ∼220 au (Figure 6), which is in any case radially localized and would not change our main results across the complete radial extent.However, IM Lup has been found to be prone to strong photo-evaporative winds, even in the presence of a weak external radiation field.Theoretical work by Haworth et al. (2017) shows that, if present, these winds would originate in the upper disk layers and possibly cause mass loss rates comparable to the disk accretion rate (10 −8 M ⊙ yr −1 , Alcalá et al. 2017), due to the large vertical extent of the disk, which makes the surface material weakly gravitationally bound.Winds are expected to broaden the line (Fuller & Myers 1992); however, as we do not detect any indication of them in the spectra or image plane, we did not include this uncertainty in our analysis.

Comparison with other methods and results
We have presented a new methodology to directly extract turbulence, a fundamental disk property, from molecular line emission.By tracing the vertical location of an optically thin tracer and comparing it to that of optically thick CO lines we can immediately access the kinetic temperature of the region (Law et al. 2021b(Law et al. , 2023;;Paneque-Carreño et al. 2022) and from the measured linewidth calculate the nonthermal contribution attributed to turbulent motions.This requires a precise parametric model of the spectra to extract the linewidth accounting for the upper and lower disk emission layers and the broadening caused by beam convolution, which is why we use DISCMINER (Izquierdo et al. 2021).
Previous works have mostly put upper limits on turbulence by modeling the emission using physical disk models (Hughes et al. 2009;Guilloteau et al. 2012;Flaherty et al. 2015Flaherty et al. , 2017Flaherty et al. , 2020)).This approach is required when the observed tracer is optically thick, as the linewidth will be affected by an opacity broadening due to the column density of the molecule (Hacar et al. 2016) and the comparison between model and data must be done taking into consideration radiative transfer calculations (Flaherty et al. 2015).However, this has the additional difficulty of determining surface density and temperature structure parameters (e.g., Hughes et al. 2009;Guilloteau et al. 2012), which are hard to constrain (see Miotello et al. 2023, and references within).Models of CO emission must also account for freezeout and desorption processes that may be relevant (Öberg et al. 2015), adding further complexity in this kind of analysis to obtain turbulence (Flaherty et al. 2020).
Analyzing CO isotopolog data in IM Lup using these complex models also results in high nonthermal motions associated with turbulence of M = 0.18-0.3(Flaherty et al. in prep).Within our errors (see Figure 7) these values are consistent with our determination of M = 0.4-0.6;however, they are a factor of 2 lower.With only one system where it is possible compare both methods, we are not able to conclude if this is a systematic difference or how it would affect the determination of nonthermal motions in other disks.We note that previous best-fit disk radiative transfer models used for studying turbulence in HD 163296 show strong differences with the disk vertical structure that is more recently constrained by Paneque-Carreño et al. (2023).Flaherty et al. (2017) predict that the expected location for CO (J = 2 − 1) in HD 163296, based on their models used for setting upper limits on the turbulence, is z/r < 0.1 (see Fig. 16 of Flaherty et al. 2017).Direct observational measurements, however, trace the CO (J = 2 − 1) emission surface in HD 163296 at z/r ∼0.3 (Law et al. 2021b;Paneque-Carreño et al. 2023), indicating that some of the underlying model assumptions may be incorrect and highlighting the need for high resolution data as well as the difficulty and possible errors of such methods.
Another approach of directly determining the thermal structure and turbulence is to obtain the temperature using multiple transitions, which is particularly useful for face-on disks where the vertical dimension is not accessible.This has been applied to the TW Hya disk to constrain an upper limit on turbulence of M ∼0.1 (Teague et al. 2016(Teague et al. , 2018b)).This method is complementary to our approach, as it relies on direct observational constraints rather than physical models; however, it carries the uncertainty on the exact emission height of different transitions.

Conclusions
We have presented new Band 7 observations at high spectral and spatial resolution of CN N = 3 − 2 and C 2 H N = 4 − 3 in IM Lup.The emission from both tracers seems to be optically thin and comes from a vertically elevated layer colocated with 13 CO J = 2 − 1.This is the first time C 2 H is seen to be vertically elevated as predicted by thermochemical models.
Through the use of a parametric model, we traced the emission morphology and constrained the linewidth of CN, which we then decomposed into its thermal and nonthermal components.Assuming that the nonthermal component corresponds to turbulent broadening, we conclude the following: 1. IM Lup presents strong turbulence in the upper layers (z/r ∼0.25), as traced by CN emission.In this layer and be-yond radial distances of 150 au, turbulent velocities are between M = 0.4 and 0.6 (α ∼10 −1 ). 2. Combined with constraints of low turbulence in the midplane (Franceschi et al. 2023;Jiang et al. 2023), IM Lup shows evidence for a vertical gradient of turbulence, with higher values in the upper layers.This is a key prediction of MRI (Simon et al. 2015), which has been proposed to be ongoing in IM Lup due to its ionization structure (Seifert et al. 2021).3. The use of optically thin tracers, with emission constrained to a localized vertical region, together with the temperature profiles extracted from CO isotopologs can be leveraged to directly measure turbulence through simple parametric models of the line intensity, without the uncertainties imposed by the density and temperature structure of protoplanetary disks.
Further use of ALMA's exquisite spectral and spatial resolution will allow us to expand this study to other sources and directly confirm if IM Lup is an exception, or if other disks have turbulent motions comparable to the sound speed in the upper layers.V 2 rot = GM * r 2 (r 2 + z 2 ) 3/2 , (D.1) where G is the gravitational constant and M * the stellar mass.Using this formula, we calculated the velocity for a parcel of material along the line-of-sight path, between the initial (r 1 , z 1 ) and final (r 2 , z 2 ) positions, as shown in Figure D.1.The velocity variation is not linear, as it depends on both radial and vertical differences; to be sensitive to these variations, we sampled the line-of-sight path every 1 au.Using a series of Gaussian profiles with line centers at each of the sampled velocities, we computed the resulting average profile, which is of Gaussian form, and estimated the measured linewidth from it (panel B in Figure D.1).As mJy km s 1 beam 1 ]

Fig. 1 .
Fig. 1.ALMA observations of CN and C 2 H line emission from IM Lup.Left column: Stacked, integrated spectra extracted from the whole radial extent of the disk, for each transition, obtained using GoFish.The corresponding hyperfine groups are indicated in the top left of each panel.Vertical lines mark the location of the central frequency for each hyperfine component within the spectral range.The channel spacing is 50 m s −1 , matching the native spectral resolution of the data.Right column: Integrated intensity maps (moment 0).

Fig. 2 .
Fig. 2. Channel map emission for CN (left) and C 2 H (right) at equal relative velocity to systemic motion (-1.2 km s −1 ).Continuous line traces the upper emission surface and dashed line the lower surface.Different colored contours indicate different hyperfine components.Contours are not fit to the data, they represent the visual identification of the emission surfaces and follow constant z/r (0.3 and 0.2 for CN and C 2 H, respectively).
CN N = 3 − 2 C 2 H N = 4 − 3 studied CO isotopolog emission from IM Lup (Paneque-Carreño et al. 2023) is presented in Figure 3 demonstrating that the location of CN

Fig. 3 .
Fig. 3. Vertical profiles of each tracer.Panels (A) and (B) show the data and best fit model for CN and C 2 H, respectively (the molecule and hyperfine group used is shown in the upper left corner).Each colored dot represents the vertical location extracted at a given radius from the maxima of emission traced in the channel maps.Overlaid is an exponentially tapered model with the best-fit parameters and the statistical uncertainty of the profile.Panel (C) presents a comparison between the averaged values and standard deviation within radial bins for CN and C 2 H (in colors) compared to the averaged CO isotopolog emission surfaces as presented in Paneque-Carreño et al. (2023).

Fig. 4 .
Fig. 4. Brightness temperature analysis for each tracer.Left panels: Peak intensity (moment 8) maps for each hyperfine group under study.The extended feature toward the semimajor axis results from the radiative transferred effects due to higher velocity density in this region.Right panel: Brightness temperature profiles of each hyperfine group, compared to those of CO (J = 2 − 1) isotopologs.

Fig. 5 .
Fig. 5. DISCMINER results and comparison to CN data.First row: Selection of channel maps from the CN datacube.Velocities are relative to the systemic velocity, and the spectral resolution of each channel is 0.1 km s −1 .Middle row: Channel maps from the best-fit DISCMINER model, at the same spatial resolution and intensity scale as the observations.Bottom row: Residual (observations minus model) emission.Contours indicate 5σ positive (red) and negative (purple) residuals, where σ is the data rms (1 mJy/beam).
Fig. 6.Comparison between various radial profiles between the data (magenta curves) and the DISCMINER model (gray curves).All panels except the height profile indicate the profile dispersion as a shaded region.Velocity and linewidth panels have a noticeable difference located at ∼220 au, a blue residual curve characterizes this feature in the upper right corner of each panel.The linewidth panel shows the measured linewidth extracted from the data and the beam-convolved model by fitting the spectra.Additionally, the dashed black line indicates the intrinsic parametric form of the linewidth, as fitted by DISCMINER.The beam-convolved model will always have a broader linewidth than the parametric form and is used to directly compare our model to the data, which has a finite sptial resolution.However, the intrinsic linewidth is what we used in our calculations, as this is expected to be the actual line broadening due to thermal and nonthermal motions.

Fig. 7 .
Fig. 7. Turbulence extracted from the nonthermal component of the CN (N = 3 − 2) linewidth.The blue curve indicates results when considering the brightness temperature of 13 CO.The orange curve traces results using the brightness temperature of 12 CO.Shaded regions represent the dispersion when considering an error of 50 m s −1 on the measured linewidth.

Fig. 8 .
Fig.8.Total intrinsic linewidth determined from the best-fit DIS-CMINER model (solid black line; this is the same as the dashed black line in the bottom-right panel of Figure6).It is compared to the thermal component obtained using the 13 CO (blue) or 12 CO (orange) brightness temperatures.

Fig. 9 .
Fig. 9. Spectra of four different pixels in the CN disk emission (in gray), indicated on the central integrated moment map.Overlaid are the DISCMINER model spectra from the best-fit parameter fit, considering the total measured linewidth (thermal and nonthermal components) in red.For comparison, the best-fit DISCMINER model considering a linewidth from only the thermal component calculated using the 13 CO temperature is overlaid* in blue.

Fig. 11 .
Fig. 11.Linewidth broadening due to the vertical Keplerian shear, for different vertical widths.Each colored line indicates a specific emission width.The vertical dashed line at 150 au indicates the radial distance beyond which we consider our nonthermal linewidth component measurements robust.Details on the method are found in Appendix D.

Fig
Fig. A.2. Gaussian fits to the radial intensity profiles of each transition.The data are shown in black, pink dashed lines trace the two individual Gaussian components, and the dashed red line presents the model curve, obtained by summing both Gaussians.The FWHM of the inner Gaussian, which traces the ring component, is indicated for each panel.

Fig
Fig. B.1.CN N = 3 − 2 emission, showing the channel maps used to extract the vertical emission layer with ALFAHOR.The masked regions are indicated in the black contours and black dots within the masks trace the emission maxima.

Fig. D. 1 .
Fig. D.1.Sketch of the method used to estimate the effect of a nonzero vertical emitting width for CN.Panel (A) shows the geometry of CN emission in IM Lup, considering an inclination of 47.8• with respect to the observer.For the vertical structure, the solid curve indicates the best-fit vertical profile of CN emission and the dashed curve the extent considering a fixed vertical width.In each case, the blue dot is the observed parcel of gas at position (r 1 , z 1 ) that has to cross the line-of-sight path as indicated by the red line, leaving the zone of CN emission at the red dot location (r 2 , z 2 ).For different radial distances, the path varies.Panel (B) shows the velocity variation between the blue dot and the final red dot positions.Each solid-line Gaussian is centered at a determined velocity, depending on their (r, z) position along the line-of-sight path.The measured linewidth (∆V measured ) is taken from the resulting average Gaussian, shown as a dashed black curve and compared to the intrinsic value (∆V gas ).

Table 2 .
Best-fit parameters of vertical profiles.