© The Authors 2024

1 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.5c_s (see also the initial studies by Guilloteau & Dutrey 1998 and Dartois et al. 2003). A later study by Flaherty et al. (2020) obtained lower, but still significant, turbulence values of 0.25–0.33c_s 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. 2015, 2017, 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 observational 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 Fig. 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 Fig. 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₂H transitions, measuring the radial profile of line broadening due to nonther-mal motions, and relating this to vertically and radially resolved turbulence. CN and C₂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₂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, 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 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 ~ 100m s⁻¹ 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₂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₂H and comparing them to those of CO isotopologs. Section 4 shows the parametric model we used to extract the linewidth from CN, and Sect. 5 presents our results on the radial turbulence profile as traced by CN. Finally, Sect. 6 discusses the findings of this work in the context of previous constrains, and Sect. 7 summarizes our main conclusions.

2 Observations

This work uses ALMA Band 7 C₂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 line-free spectral windows and channels with the tools of Common Astronomy Software Applications (CASA; version 5.4.0, McMullin 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.5 km). Three rounds of phase and one round of amplitude calibration were performed, the initial time solution interval was 4800 s 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₂H and CN hyperfine group was imaged using tclean with Briggs weighting and Keplerian masking¹. Multiple Keplerian masks, one for each of the expected hyperfine components, were used in the cleaning. A stellar mass of 1.1 M_⊙ (Teague et al. 2021) and a constant z/r value of 0.2 and 0.3 was used for C₂H and CN emission, respectively. To compromise between optimal angular resolution and high S/N, C₂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₂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⁻¹ for C₂H and 100 m s⁻¹ for CN, to optimize spectral resolution and S/N. The final rms of the images for each molecule are 0.6 and 1 mJy beam⁻¹ for C₂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⁻¹. 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.

Fig. 1 ALMA observations of CN and C2H 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).

3 Emission morphology

3.1 Vertical profile

From the channel maps it is possible to visually identify the CN and C2H 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₂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.

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 $$z(r)=z_0\times {\left(\frac{r}{100\text{ au }}\right)}^{\varphi }\times \exp \left[{\left(\frac{-r}{r_{\text{taper }}}\right)}^{\psi }\right].$$(1)

The best-fit values for z₀, ϕ, ψ, 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₂H, together with the best-fit exponentially tapered power-law model.

Comparison with previously studied CO isotopolog emission from IM Lup (Paneque-Carreño et al. 2023) is presented in Fig. 3 demonstrating that the location of CN (N = 3−2) is comparable to that of ¹³CO (J = 2−1). The C₂H (N = 4−3) emission is located slightly below, tracing the same region as C¹⁸O (J = 2−1). We note that C₂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).

Fig. 2 Channel map emission for CN (left) and C2H (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 C2H, respectively).

3.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, $$T_{\text{b}}=\frac{hv}{k}{\ln }^{-1}\left(1+\frac{2h{v}^3}{I{c}^2}\right),$$(2)

where h is the Planck constant, k the Boltzmann constant, c the speed of light, and v the frequency of the emission. To compare the brightness temperatures of CN and C₂H to those of the CO isotopologs previously studied in Paneque-Carrefio et al. (2023), all brightness temperature profiles were calculated using cubes with an equal spectral resolution of 0.2km s⁻¹ 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-Carrefio et al. 2023), are presented in Fig. 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₂H emission is spatially colocated with ¹³CO; however, both transitions show lower brightness temperatures than ¹³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).

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 ΓM Lup, in particular ¹²CO and ¹³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 C2H 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.

Fig. 3 Vertical profiles of each tracer. Panels A and B show the data and best fit model for CN and C2H, 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 C2H (in colors) compared to the averaged CO isotopolog emission surfaces as presented in Paneque-Carrefio et al. (2023).

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.

4 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.

4.1 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, 2022, 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 $$I_{\text{m}}\left(R,z;v_{\text{chan }}\right)=I_p\exp \left[-\frac{{\left(v-v_{\text{chan }}\right)}^2}{\text{Δ}{V}^2}\right],$$(3)

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. 2021, 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 Fig. 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 Fig. 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₂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 Fig. 3 for a comparison on how close the CN hyperfine components are compared to C₂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.

4.2 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 ~200 m s^–1 residual between 200 and 250 au, which coincides with a ~50 m s^–1 negative residual in the linewidth. This velocity variation is also observed in ¹²CO (see Fig. 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.

5 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\sqrt{\log 2}\text{Δ}V$. The linewidth can then be written as $$\text{Δ}{V}^2=V_{\text{th }}^2+V_{\text{non-th }}^2.$$(4)

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 $$\text{Δ}{V}^2=\frac{2k{T}_{\text{kin }}}{m_X}+v_{\text{turb}}^2.$$(5)

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, 𝓜, such that $$ℳ=\frac{v_{\text{turb }}}{\sqrt{2}{c}_s},$$(6)

where c_s is the local sound speed in the studied layer $c_s^2=\frac{k{T}_{\text{kin }}}{\mu m_{\text{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, 2017, 2020; Teague et al. 2016, 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 Fig. 7. Our results indicate turbulent motions of 𝓜 = 0.4–0.6 at radii beyond ~l50au. 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 <5m 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 (v = ac_sH, Shakura & Sunyaev 1973) it is possible to relate the turbulence Mach number 𝓜 to α following, α ~ 𝓜² (Rosotti 2023). Our measurements indicate α ~ 10^–1.

The total measured linewidth (ΔV) and turbulence values increase steeply in the inner ~l50au. 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 ~22l m s^–1. Thermal broadening only accounts for 105-132 m s^–1, using the ¹³CO and ¹²CO brightness temperatures, respectively. Therefore, the nonthermal broadening ranges between 177 and 194m 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).

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–1).

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

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

6 Discussion

6.1 Vertical location of CN and C₂H

Our work traces the CN and C₂H emission originating from a vertical region similar to ¹³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₂H coming from vertically elevated layer in the outer disk, as previous works studying other protoplanetary disks had only traced C₂H inward of 150 au and closer to the midplane (Law et al. 2021b; Paneque-Carrefio et al. 2023).

Both CN and C₂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₂H. IM Lup is the perfect laboratory for these studies.

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 13CO temperature is overlaid in blue.

6.2 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 (<20au, 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 ΓM 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 (~3l3 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 ther-mochemical 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. 2024). 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⁻¹) traced by CN at R >l50au 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). ΓM 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₂CO, and HCN, which have their vertical heights characterized and emit from in-between the midplane and ¹³CO (Paneque-Carreño et al. 2023). Unfortunately, the sensitivity and spectral resolution of the available data for these tracers (≤200 m s^–1, Öberg et al. 2021) is not sufficient to detect nonthermal broadening of ≤90 m 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¹⁸O J = 2–1 data from the MAPS program (Öberg et al. 2021). Figure 3 shows that C¹⁸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¹⁸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¹⁸O are slightly higher, but comparable within the uncertainty to the results obtained from CN assuming the temperature profile of ¹³CO. We note that while the brightness temperature of C¹⁸O is below the assumed excitation temperature (see Fig. 4), the column density of CO in IM Lup (Zhang et al. 2021) and estimates from Elias 2–27 of C¹⁸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¹⁸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¹⁸O considering opacity broadening gives a turbulent dispersion profile below that of CN (using the ¹³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.

Fig. 10 Same as Fig. 7, but including an analysis of the nonthermal component of C18O J = 2–1 data. The solid black line traces the results when applying the same methodology as previously used in CN to C18O. The dashed black line is the turbulence profile obtained after assuming a 25% opacity broadening in the measured linewidth, which may be the case if C18O has optical depths of 1–2.

6.3 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 Fig. 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 ΓM 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 <25m s^–1. The width would have to be broader than ~60au 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 ~220au (Fig. 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.

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.

6.4 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, 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. 2015, 2017, 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 freeze-out 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 𝓜 = 0.18–0.3 (Flaherty et al., in prep.). Within our errors (see Fig. 7) these values are consistent with our determination of 𝓜 = 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 𝓜 ~ 0.1 (Teague et al. 2016, 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.

7 Conclusions

We have presented new Band 7 observations at high spectral and spatial resolution of CN N = 3−2 and C₂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¹³ CO J = 2−1. This is the first time C₂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 beyond radial distances of 150 au, turbulent velocities are between 𝓜 = 0.4 and 0.6 (α ~ 10⁻¹).

• 2.

  Combined with constraints of low turbulence in the mid-plane (Franceschi et al. 2023; Jiang et al. 2024), 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.

Acknowledgements

This paper makes use of ALMA data from project #2019.1.01357.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), and by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101019751 MOLDISK). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 823823 (DUSTBUSTERS).

Appendix A Radial features

Fig. A.1 Azimuthally averaged radial intensity profiles for each of the observed hyperfine groups. The solid line traces the average value and the shaded region the dispersion within each radial bin. The vertical black line indicates the peak of emission, which coincides for all transitions at ~106 au.

The integrated intensity emission maps (moment 0) show a distinct radial feature for all the studied hyperfine groups (Fig. 1). We compute the azimuthally averaged intensity profiles from these maps by accurately deprojecting using the best-fit vertical surfaces and the function radial_profile from GoFish (Teague 2019). Only emission from within a 60° wedge along the semi-major axis is considered, to avoid deprojection effects from emission close to the minor axis (Law et al. 2021a). Our results are shown in Figure A.1, where all tracers display a ringed structure peaking at ~ 106 au.

To characterize the width and any azimuthal variations of the ring feature, a set of two Gaussians are fitted to the radial profiles as done in Law et al. (2021a). A broad component Gaussian is used to trace the emission at further out radii and a second Gaussian traces the ring structure. Figure A.2 shows the fitted Gaussians and FWHM value of the ring-component Gaussian. We used a Gaussian function such that $$I=I_0\exp \left[-\frac{1}{2}\frac{{\left(R-R_0\right)}^2}{{\sigma }^2}\right],$$(A.1)

where I is the measured intensity, R the radial distance, I₀, R₀ and σ are the fitted parameters and the $\text{FWHM}=2\sqrt{2\log (2)}\sigma .$ While the peak of emission was coherent across different hyperfine groups, the ring width has slight variations. Overall the feature shows an average width of ~89 au in CN and ~109 au in C₂H.

Considering the broadest value of ~115 au we calculate the averaged intensity throughout the angular extent of the ring, using azimuthal bins of 5°. Figure A.3 shows the results. CN displays a symmetric ring across azimuth, with slight variations along the semiminor axis likely due to projection effects. C₂H has a lower S/N and may show tentative asymmetries; however, the quality of the data is not sufficient to characterize them in detail.

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. A.3 Azimuthal emission curves along the ring for each transition. The solid line indicates the mean values within each azimuthal bin and the shaded region the dispersion. The panel to the left shows the radial limits used to extract the ring emission.

Appendix B Channel map analysis

To extract the vertical location of the emission layer, we traced the emission maxima from the channel map emission as explained in section 3.1. The channels, masks, and extracted emission maxima are presented in Figures B.1 and B.2 for CN and C₂H, respectively.

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. B.2 Same as Fig. B.1 but for C2H N = 4–3 emission.

Appendix C DISCMINER best-fit parameters

Table C.1 indicates the best-fit parameters obtained by DISCMINER for the CN emission model.

Appendix D Effect of vertical layer width

To estimate the effect of the Keplerian shear due to a nonzero vertical width of the CN emitting layer, we considered the best-fit CN emitting surface as a lower limit and used the geometry of the disk to calculate the line-of-sight path for different vertical widths (see panel A in Figure D.1). The rotation velocity of a parcel of gas located at a position (r, z), in cylindrical coordinates, from the star is given by $$V_{\text{rot }}^2=\frac{G{M}_{*}{r}^2}{{\left(r^2+z^2\right)}^{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₁, z₁) and final (r₂,z₂) 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 the temperature is not expected to change significantly in these upper disk layers for IM Lup (Cleeves et al. 2016), we assumed a constant intrinsic thermal linewidth in the Gaussian profiles for our estimation. The final linewidth broadening caused by a nonzero emission surface is then calculated by subtracting the intrinsic and measured linewidth values.

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 (r1, z1) 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 (r2, z2) 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 (ΔVmeasured) is taken from the resulting average Gaussian, shown as a dashed black curve and compared to the intrinsic value (ΔVgas).

References

All Tables

All Figures

Fig. 1 ALMA observations of CN and C2H 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).

In the text

Fig. 2 Channel map emission for CN (left) and C2H (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 C2H, respectively).

In the text

Fig. 3 Vertical profiles of each tracer. Panels A and B show the data and best fit model for CN and C2H, 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 C2H (in colors) compared to the averaged CO isotopolog emission surfaces as presented in Paneque-Carrefio et al. (2023).

In the text

	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.
In the text

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–1).

In the text

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.

In the text

	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 13CO. The orange curve traces results using the brightness temperature of 12CO. Shaded regions represent the dispersion when considering an error of 50 m s–1 on the measured linewidth.
In the text

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

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 13CO temperature is overlaid in blue.

In the text

	Fig. 10 Same as Fig. 7, but including an analysis of the nonthermal component of C18O J = 2–1 data. The solid black line traces the results when applying the same methodology as previously used in CN to C18O. The dashed black line is the turbulence profile obtained after assuming a 25% opacity broadening in the measured linewidth, which may be the case if C18O has optical depths of 1–2.
In the text

	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.
In the text

	Fig. A.1 Azimuthally averaged radial intensity profiles for each of the observed hyperfine groups. The solid line traces the average value and the shaded region the dispersion within each radial bin. The vertical black line indicates the peak of emission, which coincides for all transitions at ~106 au.
In the text

	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.
In the text

	Fig. A.3 Azimuthal emission curves along the ring for each transition. The solid line indicates the mean values within each azimuthal bin and the shaded region the dispersion. The panel to the left shows the radial limits used to extract the ring emission.
In the text

	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.
In the text

	Fig. B.2 Same as Fig. B.1 but for C2H N = 4–3 emission.
In the text

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 (r1, z1) 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 (r2, z2) 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 (ΔVmeasured) is taken from the resulting average Gaussian, shown as a dashed black curve and compared to the intrinsic value (ΔVgas).

In the text