Accurate measurements of the physical structure of protoplanetary discs are critical inputs for planet formation models. These constraints are traditionally established via complex modelling of continuum and line observations. Instead, we present an empirical framework to locate the CO isotopologue emitting surfaces from high spectral and spatial resolution ALMA observations. We apply this framework to the disc surrounding IM Lupi, where we report the first direct, i.e. model independent, measurements of the radial and vertical gradients of temperature and velocity in a protoplanetary disc. The measured disc structure is consistent with an irradiated self-similar disc structure, where the temperature increases and the velocity decreases towards the disc surface. We also directly map the vertical CO snow line, which is located at about one gas scale height at radii between 150 and 300 au, with a CO freeze-out temperature of 21 ± 2 K. In the outer disc (>300 au), where the gas surface density transitions from a power law to an exponential taper, the velocity rotation field becomes significantly sub-Keplerian, in agreement with the expected steeper pressure gradient. The sub-Keplerian velocities should result in a very efficient inward migration of large dust grains, explaining the lack of millimetre continuum emission outside of 300 au. The sub-Keplerian motions may also be the signature of the base of an externally irradiated photo-evaporative wind. In the same outer region, the measured CO temperature above the snow line decreases to ≈15 K because of the reduced gas density, which can result in a lower CO freeze-out temperature, photo-desorption, or deviations from local thermodynamic equilibrium.

These dynamical processes remain hard to constrain by existing observations, however. Rotational molecular lines are a powerful tool that can be used to probe the disc structure, but constraints rely in most cases on complex model fitting or on indirect methods. For instance, detailed radiative transfer modelling of multiple CO lines can provide an estimate of the vertical temperature structure of the disc (Dartois et al. 2003), indicate freeze-out of the gas phase molecules onto dust grains in the cold midplane (e.g. Qi et al. 2011; Zhang et al. 2017), and suggest a very low level of turbulence (e.g. Flaherty et al. 2015). Indirect mapping via chemical imaging (Qi et al. 2013, with N2H+; and Mathews et al. 2013, with DCO+) has also been used to estimate the location of the CO snow line. However, Aikawa et al. (2015), van’t Hoff et al. (2017), and Huang et al. (2017) have shown that these species are not robust CO snow line tracers, and that detailed chemical modelling is required to infer the midplane CO snow line location from N2H+ or DCO+ observations.

Channel maps of 12CO (left) and 13CO (right) displayed with a channel width of 80 m/s. The data is not continuum subtracted. The dv is relative to the systemic velocity of 4.5 km s-1. The beam is indicated in the bottom right corner.

The combination of high spatial and spectral resolution and sensitivity offered by ALMA opens new avenues to directly map the disc thermal and kinematic structure by resolving the gas disc both radially and vertically. Dutrey et al. (2017) recently introduced a method to map the thermal and density gas structure of discs at close to edge-on inclinations. Discs at intermediate inclinations are also ideal targets as the Keplerian velocities spatially separate the emitting regions, eliminating line of sight confusion (e.g. de Gregorio-Monsalvo et al. 2013; Rosenfeld et al. 2013).

IM Lupi is a M0V T Tauri star (distance of 161 ± 10 pc, Gaia Collaboration 2016) surrounded by a large and massive disc with an inclination of 48° ± 3° (Pinte et al. 2008; Cleeves et al. 2016). The disc is detected in rotational CO emission up to a radius of 970 au, in the millimetre continuum up to 310 au (Panić et al. 2009; Cleeves et al. 2016), and in scattered light up to 320 au with a well-defined disc morphology (Pinte et al. 2008, rescaled to the updated distance), but with an extended component up to 720 au in radius. We present here a framework to directly measure the altitude, velocity, and temperature of the CO layers from new high spectral (0.05 km s-1) and intermediate spatial (0.4′′) resolution ALMA observations of the disc surrounding IM Lupi.

2. Observations and data reduction

IM Lupi was observed with ALMA in band 6 on the night from 9 to 10 June 2015 with a total on-source time of 37.4 min (Program ID 2013.1.00798.S). The array was configured with 37 antennas and baselines ranging from 21.4 to 784 m. Titan was used as a flux calibrator and the quasars J1517-2422 and J1610-3958 were used as bandpass and phase calibrators. Two of the four spectral windows of the correlator were configured in time division mode (TDM) centred at 218.9 GHz and 231.7 GHz, each with 1.875 GHz usable bandwidth. The other two spectral windows were configured in frequency division mode (FDM) to target the 12CO J = 2−1, 13CO J = 2−1, and C18O J = 2–1 transitions, with a channel spacing of 30.5 kHz and a channel averaging of 2. Because of the Hanning smoothing used to reduce spectral ringing, the actual spectral resolution is 35 kHz, corresponding to 45–48 m/s depending on the line. The median partial water vapour was 0.93 mm. The data sets were calibrated using the Common Astronomy Software Applications pipeline (CASA, McMullin et al. 2007, version 4.6.0). We performed three successive rounds of phase self-calibration, the last with solutions calculated for each individual integration (6s), followed by a phase and amplitude self-calibration. The continuum self-calibration solutions were applied to the CO lines. A CLEAN mask was manually defined for the continuum image and each channel map. Channel maps were produced both with and without continuum subtraction (using the CASA task uvcontsub). We estimate the flux calibration uncertainty to 10%.

The disc is spatially resolved in all transitions and channel maps show the typical butterfly pattern of discs in Keplerian rotation. The measured fluxes are presented in Table 1 and the representative channel maps for 12CO and 13CO in Fig. 1. The C18O channel maps and continuum image are shown in Figs. A.1 and C.1, respectively. Our data and results are consistent with those of Cleeves et al. (2016). We focus in the following on the spatial distribution of the emission in individual channel maps.

3. Reconstructing the altitude, velocity, and temperature of the CO emitting layers

In a given channel map, the line emission is concentrated along the isovelocity curve, i.e. the set of points where the projected velocity of the emitting surface is constant. As a result, the emission originating from the upper and lower surfaces (i.e. above and below the midplane as seen from the observer), and from the near and far sides of these surfaces, is well separated (Figs. 2 and 1). By locating the emission in each channel map, we can directly reconstruct, via simple geometrical arguments, the position and velocity of each of the CO layers.

3.1. Altitude

We consider the coordinate system where the x-axis is aligned with the disc major axis, and we denote (x⋆, y⋆) the position of the star which was determined by the peak of the continuum map (Fig. 2). Because the gas is vertically pressure supported, we assume that, at any point in the disc, the gas is rotating on a circular orbit parallel to the disc midplane, i.e. we assume that the bulk of the motion is described by Keplerian rotation and neglect for this analysis any radial or vertical gas circulation. The radius and altitude of this orbit are denoted by r and h, respectively.

For a given offset δx = x−x⋆ in the image plane, a vertical line will intersect the isovelocity curve of the upper disc surface1 at two points (for small enough x), which belong to the same orbit, i.e. at the same distance from the star and at the same altitude. Due to line broadening, and to finite spatial and spectral resolution, the emission is not located exactly on the isovelocity curve, but forms a narrow band around it. We estimate the position of the emission by its maximum, as illustrated in Fig. 2. If we denote yn and yf the ordinates of the two points on the near and far side of the disc, the coordinates of the centre of the projected circular orbit passing through those two points are (xc, yc), where xc = x⋆ and yc = (yf + yn)/2. Any point on this ellipse fulfils (x−x⋆)2 + ((y−yc)/cosi)2 = r2, where the orbital radius r is the length of the semi-major axis of the ellipse and i the disc inclination.

By using the point on the far side of the disc surface, we can for instance derive (1)The altitude h of the orbit is derived by noting that the offset between the centre of the ellipse and the star is simply hsini: (2)

3.2. Gas velocity

In a given velocity channel, the projected radial velocity νobs is known and can be expressed as νobs = νsyst + ν cosθ sini, where ν is the actual gas velocity around the central star, νsyst the systemic velocity, and θ is the true longitude. By noting that cosθ = (x−x⋆) /r, we obtain the actual gas velocity at each point: (3)

3.3. Gas kinetic temperature

For optically thick emission, and as long as the emission fills the beam, the observed brightness temperature Tb is very close to the excitation temperature Tex, Tb = Tex(1−exp(−τ)). Low J CO lines are close to local thermodynamic equilibrium, with an excitation almost equal to the gas kinetic temperature Tkin (e.g. Pavlyuchenkov et al. 2007). For optically thick lines, the peak intensity then provides an accurate estimate of the gas temperature. This is the case for the 12CO and 13CO lines (see discussion in Sect. 4), but the C18O emission is partly optically thin and the brightness temperature is lower than the actual excitation temperature. To reconstruct the disc temperature profile, we measured the brightness temperature at all points extracted above. We also extracted the brightness temperature of points on the disc lower surface in the same way (Fig. 2).

3.4. Reconstructing the full spatial, kinematic, and thermal structure of the CO layers

As described above, by locating the position of the maxima in each isovelocity curve and for each x in a given channel map (where we sampled x every 0.08′′, i.e. ≈1/4 of the beam), we can directly measure r and the corresponding elevation h and velocity ν, and we can estimate Tgas.

Each channel probes a range of distances from the star depending on the velocity. By repeating this procedure for each channel, we can reconstruct the complete 2D velocity and temperature distribution of the CO emitting layers, as shown in Fig. 4.

3.5. Caveats

The central channels of 12CO (within ±0.7 km s-1 from the systemic velocity) suffer from strong foreground extinction (Panić et al. 2009; Cleeves et al. 2016). They were excluded from the analysis for the temperature estimate. The morphology of the channel maps is not affected, however, and we included them to measure the altitude and velocity.

The method presented in this section can only be applied in regions where the emission originates from a vertically thin layer, so that an altitude h can be defined, i.e. a clear maximum can be detected in the channel map. This is only possible when the projected width of the emitting layer is significantly smaller than the beam. This is the case for instance when the molecular layer itself is a narrow layer, or when the optically depth gradient is steep around the τ = 1 surface and we can only see a narrow layer.

These conditions are met for the CO layers of IM Lupi up to a radius 450 au, where we are able to distinctly identify the location of the isovelocity curves, but are no longer valid far from the star even though the emission is detected up to ≈1000 au. As already noted by Cleeves et al. (2016), the CO emission at large scales (> 450 au) becomes very diffuse, making it impossible to locate the isovelocity curve and to estimate an altitude. This is very likely the result of UV photo-desorption of the CO by external irradiation from the interstellar radiation field (see e.g. Fig. 5d), allowing gas CO to exist even in the midplane and to emit from a vertically extended region.

The presented framework also requires sufficient signal in a given channel to accurately locate the emission. For 12CO we were able to use the intrinsic channel width, but for 13CO and C18O we needed to bin the data to channel widths of 80 and 320 m/s, respectively. This introduces spatial smearing due to the disc Keplerian velocity and reduces correspondingly the range over which we can accurately measure the CO layer altitude2. For 13CO and C18O, this prevented a measurement of the altitude at radii smaller than ≈ 150 au and ≈ 220 au, respectively.

Continuum subtraction can significantly affect the measured brightness temperature and apparent morphology of line emission when the line is optically thick, and the line and continuum intensities are comparable (Boehler et al. 2017). At the centre of an optically thick line, the dust continuum is absorbed by the molecule. As the continuum is usually estimated from line-free channels, this results in oversubtraction of the continuum close to the line centre. The authors conclude that using non-continuum-subtracted maps is preferable in order to estimate the line brightness temperature (see e.g. their Fig. 8). This conclusion also holds for the framework we have presented, and using non-continuum-subtracted channel maps will provide a more accurate estimate of the brightness temperature and will prevent artificial spatial shift of the line emission due to oversubtraction of the local continuum emission.

At the spatial resolution of our observations, the continuum subtraction has no impact on the results as the continuum is optically thin in the regions where we can isolate the various isovelocity curves, and the corresponding continuum intensity is negligible compared to the line intensity. As a sanity check, we performed our analysis on both the continuum-subtracted and non-subtracted maps and did not find any significant difference. At higher spatial resolution, however, the effect is very likely to become important as the continuum emission in our current central beam is significant, producing a central “hole” in the continuum-subtracted 13CO and C18O maps (Öberg et al. 2015; Cleeves et al. 2016), which is absent in the non-continuum-subtracted maps (Figs. 1 and A.1).

4. Results and discussion

4.1. CO layers

As expected, we find that the 12CO emitting layer lies higher in the disc than the 13CO, due to its much higher optical depth (Fig. 4, left panel). This is readily visible in Fig. 1, where the two surfaces are geometrically much closer to one another for 13CO than 12CO. While already suggested by previous ALMA observations (e.g. de Gregorio-Monsalvo et al. 2013; Rosenfeld et al. 2013, for HD 163296), this is the first direct, model-independent measurement of the altitude of the CO layers. In particular, we clearly see that CO surfaces are flaring and seem to flatten out outside of ≈ 300 au. At a distance of 200 au, the 12CO emission originates from an altitude of ≈65 au, while the 13CO originates from ≈25 au. This corresponds to about 2.5 and 1 hydrostatic scale height, respectively. Interestingly, these values are very similar to the numbers that Dartois et al. (2003) and Dutrey et al. (2017) have found for DM Tau and the Flying Saucer, respectively.

A linear regression analysis indicates a power-law exponent of 1.8 ± 0.2 and 2.1 ± 0.4 for the altitude of 12CO and 13CO layers at radii smaller than 300 au, respectively. It should be noted that these values are much larger than the actual faring index of the disc pressure scale height (≈1.15–1.2 in the various published disc models). This is expected, however, as the shape of the emitting CO layer is set by the irradiation from the star (responsible for both the photo-dissociation and disc temperature). The emitting molecular layer should roughly follow a layer with a constant optical depth as seen from the star, and not the actual disc scale height. Similarly, large exponents were measured around the Herbig star HD 97048 for the PAH emission (Lagage et al. 2006) and scattered light (Ginski et al. 2016) layers, which are also set by the stellar irradiation.

The C18O emission is more difficult to locate due to the lower signal-to-noise ratio and required velocity binning, but it appears co-spatial with the 13CO emission. Since 13CO and C18O have a typical optical depth ratio of about 8:1, one expects the C18O to be closer to the midplane than 13CO by the same reasoning as for 13CO and 12CO. The fact that we do not see any significant difference in altitude between 13CO and C18O suggests that we are actually probing a physical layer in the disc below which there is no significant CO emission, hence this layer marks the vertical freeze-out line.

Single 12CO channel with a schematic of the quantities we measure. The disc has been rotated so that the semi-major axis is horizontal. The position of the central object is marked by a star. For a given offset δx along the disc major axis, two maxima N and F on the same vertical are located on the near and far sides of the disc’s upper surface. These two points are on the same inclined circular orbit (white ellipse), whose centre is marked by a diamond. Indicated are the geometrical parameters from which the cylindrical radius, height, and velocity of the points F and B can be measured. By repeating the procedure for each δx (green points, sampled every 0.08′′) and each velocity channel, a full mapping of the CO upper surface is performed. The cyan points mark the location of the lower surface (near side), while the magenta points represent the symmetric of the far upper surface relative to the disc major axis.

4.2. CO depleted midplane

The disc lower surface (i.e. on the other side of the disc midplane) is clearly visible, in particular for 12CO, telling us that the CO gas abundance is significantly reduced in the disc midplane. To illustrate this, we generated a representative disc model (Fig. 5) using the radiative transfer code MCFOST (Pinte et al. 2006, 2009). This model qualitatively explains the observed morphological features, but is not in any sense optimized and should only be considered for general guidance. We also produced a few additional models where we modified physical ingredients (UV chemistry, CO condensation) to discuss their impact on the observed channel maps (Fig. 5). The model description and parameters are presented in Appendix B. In particular, we see that the CO needs to be physically depleted in the disc midplane to be able to observe the lower disc surface. The temperature gradient alone (Fig. 5, panel a) results in a strong attenuation by the cold CO in the midplane, masking the lower surface.

Because the 12CO emission on the lower surface is seen through the disc, it should originate from the lowest boundary of the actual CO layer, i.e. exactly at the vertical CO snow line. It can potentially provide a more accurate estimate of the altitude of the CO snow line, and of the CO layer thickness by comparison with measurements from the upper surface. We attempted to locate the 12CO lower surface using the same formalism as for the upper surface. But, surprisingly, the lower surface appears shifted further away from the midplane than the upper surface is (Fig. 2), indicating an apparent higher altitude, while one would expect the opposite (see Fig. 3). The 13CO does not show any shift on the other hand. The absence of shift is consistent with the 13CO emission probing the actual physical layer of CO. We interpret the apparent shift of the lower 12CO surface as a radiative transfer effect.

Measured altitude, velocity, and brightness temperature of the CO layers as a function of radius: 12CO (blue), 13CO (red), and C18O (green). Data points were extracted for all channel maps and binned according to distance. The error bars represent the dispersion within the bins added quadratically to the uncertainty resulting from inclination (48 ± 3°). Left: solid lines display a broken power-law fit (below and above r = 300 au) for the 12CO and 13CO. Centre: black solid line represents the Keplerian velocity for a 1 M⊙ central star, derived from the C18O position-velocity diagram (Appendix A). The shaded area represents an uncertainty of 0.1 M⊙. The red and blue dashed lines show the expected Keplerian velocities taking into account the altitude of the 12CO and 13CO layers. The red and blue solid lines display the velocity of the 12CO and 13CO layers taking into account the altitude and the pressure gradient term, as derived from our MCFOST model. Right: circles and triangles represent the temperature measured on the upper and lower surfaces, respectively. The brightness temperatures are similar on the upper and lower surfaces for 13CO and C18O; for clarity, only the upper surface is plotted.

Representative models of IM Lupi. The disc has been rotated so that the semi-major axis is horizontal. From left to right a): without any freeze-out of CO; b): with complete freeze-out of CO, where T< 21 K; c): with depletion of gaseous CO by 10-4, where T< 21 K; d): same as c) including the interstellar radiation field. The models are designed for qualitative comparison only. They are convolved to the ALMA resolution, but not degraded to the same noise level as the observations.

If the CO is entirely depleted at temperatures below 20 K, the lower disc surface appears very wide and uniformly bright in the model, with no shift relative to the upper surface (Fig. 5b). To reproduce the observed morphology, we need a partial depletion of ≈ 10-4 (i.e. an abundance of ≈ 10-8) of CO where T< 20 K (Fig. 5c). The exact value is not well constrained and depends on other parameters, in particular the local turbulent velocity dispersion. This small fraction of cold CO gas below the snow line attenuates the centre of the emission from the lower surface. Emission can only escape where the velocity shift is large enough, i.e. slightly away from the isovelocity curve, resulting in an apparent spatial shift of the lower surface. In our simple model, this is reproduced by a constant depletion factor. In reality, this means that the vertical snow line is not infinitely sharp, as expected, but with a CO abundance that decreases progressively with altitude, i.e. over a few K. It is potentially a powerful diagnostic to constrain the physics of the CO condensation on grains, and the possible vertical “smearing” of the ice line due to turbulent vertical mixing (Aikawa 2007; Furuya & Aikawa 2014).

4.3. Vertical velocity gradient and sub-Keplerian rotation

The different altitudes probed by the CO isotopologues also allow for the first measurement of the vertical gradient on the rotational velocity of the disc. At high altitude, the 12CO emitting layer is rotating slower than the 13CO/C18O emitting layer, which is located closer to the midplane (Fig. 4, central panel).

Within 300 au, the difference in the measured azimuthal velocities between 12CO and 13CO/C18O is in good agreement with the vertical dependence of the stellar gravitational field (Fig. 4), where the radial projection of the gravitational force per unit mass decreases with altitude: (4)At larger radii, the gas becomes clearly sub-Keplerian. This can be understood by including the pressure gradient term (e.g. Rosenfeld et al. 2013): (5)With both the density and temperature profile generally decreasing with distance, this term is negative, reducing the gas velocity compared to the fiducial Keplerian velocity. This phenomenon is thought to play an important role in the radial migration and growth of solids in these discs (Weidenschilling 1977; Barrière-Fouchet et al. 2005; Brauer et al. 2008).

Interestingly, the radius at which we see significant sub-Keplerian velocities coincides with the radius where the gas surface density transitions from a power law to a tapered profile (Panić et al. 2009; Cleeves et al. 2016), i.e. where the disc pressure gradient is becoming larger and where the velocities should deviate the most from Keplerian. At a radius of 400 au, the measured deviations from the altitude dependent Keplerian rotation are of the order of 0.15 km s-1 for both 13CO and 12CO. This is consistent with the expected pressure gradient as derived from the temperature and density structure of our MCFOST model (Fig. 4), which predicts a slightly smaller correction of ≈ 0.1 km s-1 at 400 au. This value depends on the local density and temperature gradients, which are not well constrained at these distances from the star, and the pressure gradient term might be able to fully explain the observed deviations. We note that the deviations from Keplerian velocities start at slightly different radii for the two isotopologues, ≈300 au for 12CO and ≈350 au for 13CO. The difference is about one beam size and would need to be confirmed with higher spatial resolution observations. This behaviour is not captured by our model, suggesting that an additional process might be at play.

A possible explanation is that we could also be observing the base of a photo-evaporative wind. Haworth et al. (2017) showed that the large-scale halo of CO emission (> 1000 au) described by Cleeves et al. (2016) could be interpreted as the result of a photo-evaporative wind created by irradiation from the local (weak) external radiation field. The CO abundance in the flow remains high enough to be comparable to the observed large-scale emission around IM Lupi. Interestingly, as shown by Facchini et al. (2016) and Haworth et al. (2016), such a photo-evaporative wind can also result in a sub-Keplerian rotation field (see e.g. Fig. 10 of Haworth et al. 2016), between the outer edge of the disc and the H-H2 transition. For the 1D models of IM Lupi presented by Haworth et al. (2017), the expected deviation from Keplerian rotation is at the 0.1–0.2 km s-1 level (Haworth, priv. comm.) and would also be consistent with our observations. Such a photo-evaporative wind could also explain why the deviations from Keplerian rotation seem to start closer in at high altitude where the density is lower, allowing the external radiation to penetrate further in. However, both Facchini et al. (2016) and Haworth et al. (2016) predict that the wind should also be associated with a radially increasing gas temperature profile. This is clearly not the case for IM Lupi, and more detailed modelling (e.g. 2D) might be necessary to test the hypothesis of external photo-evaporation further.

4.4. Two-dimensional temperature structure and vertical CO snow line

In the upper disc surface, we find that the 12CO temperature is decreasing with radius, as expected from direct heating by the star (Fig. 4, right panel and Fig. 6). In the lower surface, on the other hand, the 12CO temperature is virtually constant at ≈21 K, between 130 and 300 au. Because we are looking at this surface through the disc, we are looking at the emitting CO layer closest to the midplane. At these large distances from the star, the dust continuum is optically thin (see Appendix C) and the measured CO brightness temperature is attenuated by at most 3% by continuum absorption. Because the whole CO layer is optically thick, we are observing its top in the upper surface, and its bottom in the lower surface. The difference in brightness temperature between the two sides of the disc is then a direct observation of the vertical temperature gradient.

Comparison of the measured temperature (points) with the temperature structure from our MCFOST model (background). The same colour map is used for the data and the model. The white contour is the T = 21 K level, and the dashed contours correspond to T = 15, 30, 40, and 50 K.

The constant 21 K of the lower surface is consistent with the expected CO condensation temperature for binding onto pure CO ice (e.g. Sandford & Allamandola 1988; Bisschop et al. 2006; Martín-Doménech et al. 2014; Fayolle et al. 2016). We interpret it as a strong indication that we are seeing gaseous CO just above the freeze-out region. This is reinforced by the measured brightness temperature of the 13CO on the upper surface, which is extremely close to the 12CO brightness temperature on the lower surface (Fig. 4, right panel). The combination of 12CO and 13CO allows us to probe both sides of the gaseous CO just above the freeze-out region. In other words, our method allows us to directly map the vertical snow line as a function of radius. The very similar temperature profiles, despite a typical factor ≈ 70 in abundance, indicate that the lines remain optically thick and that we are indeed measuring the gas temperature in this region, hence directly measuring the CO condensation temperature.

The C18O emission has a brightness temperature of ≈8 K on both surfaces between 150 and 300 au. With an excitation temperature of 21 K as for 13CO, this means that the C18O optical depth is of the order of 0.5. With a typical abundance ratio of 8:1, the corresponding 13CO optical depth is ≈ 4, confirming our assumption of Tex = Tkin. This moderate optical depth also explains why we observe the gas 13CO in the dense region just above the CO snow line, but not at higher altitudes where it becomes optically thin (while 12CO remains optically thick) due to the reduced densities.

At distances greater than 300 au, where the gas density drops, the temperature of the CO region above the snow line appears to decrease (as measured by both 12CO and 13CO) with distance down to ≈15 K at 400 au. This is consistent with the expected dependence of the CO condensation temperature on the density (e.g. Collings et al. 2003; Bisschop et al. 2006; Hollenbach & Gorti 2009). Detailed modelling by Cleeves et al. (2016) has shown that UV photo-desorption naturally explains the diffuse CO emission outside 400 au (see also Fig. 5d). We might also observe this effect between 300 and 400 au, where the increasing external radiation field can allow the CO to remain in the gas phase at lower and lower temperatures as the density decreases with radius. Another possibility is that we are starting to see deviations from local thermodynamic equilibrium as the density becomes progressively lower than the critical density, and the measured excitation temperature can become lower than the actual kinetic temperature.

Closer in, the emission along the isovelocity curve becomes very compact and no longer fills the beam, resulting in a drop in the brightness temperature which can no longer be used to estimate the gas temperature. Interestingly, this drop does not occur at the same radius for the high altitude 12CO emitting region (r ≈ 160 au) and for the CO layer just above the snow line (r ≈ 130 au). This is a direct indication that, in this range of radius, the vertical temperature gradient is becoming steeper at high altitude in the disc, as expected from models of irradiated discs (Fig. 6). When the beam becomes larger than the actual emitting region, a fraction of the beam is filled with emission originating from the wings of the local spectral line, i.e. with lower optical depth, meaning that part of the emission originates from further away along the line of sight. If the temperature gradient is shallow, the emission will remain at the same level, as seen for 13CO, C18O, and 12CO on the lower disc surface. This is consistent with the similar temperature we obtain while observing this layer from both sides with 13CO and 12CO.

5. Concluding remarks

We have presented a general framework to measure directly the altitude, velocity, and temperature of the CO emitting layers in protoplanetary discs at intermediate inclination. These simple geometrical considerations were applied to IM Lupi. Because these measurements are performed directly on the ALMA channel maps and do not rely on any modeling of the disc, they put unbiased constraints on the disc structure, which in turn are critical to feed models of disc structure, dust evolution, and early stages of planet formation.

The results are in agreement with the broad picture of how we understand protoplanetary discs, with a tapered-edge disc structure passively heated by the central star. The accuracy of our measurements is currently limited by the spatial resolution of ≈ 0.4′′of our data. The presented framework is generic. We applied it to IM Lupi which holds an unusually large disc. However, the required high spectral resolution (≲ 0.1 km s-1) can be reached by ALMA down to spatial resolutions of ≈ 0.1–0.2′′in a few hours to detect emission with Tb ≈ 20 K, making it possible to reproduce this analysis for a significant number of discs and constrain their thermal and kinematic structures in a consistent way.

At these high spatial resolutions, discs very often display structures, such as dust rings and gaps (ALMA Partnership et al. 2015; Isella et al. 2016; Fedele et al. 2017; van der Plas et al. 2017). The presented framework may offer a way to constrain local variations in altitude, temperature, and/or velocity of the CO layers in the vicinity of the observed structures, potentially pinpointing the mechanisms at the origin of these structures.

1

The same formalism can be used for the lower surface, but for reasons we detail in Sect. 4, it cannot be used here and we only describe the formalism for the upper surface.

2

This mainly affects our ability to measure h, which depends on the small distance between the star’s position and the centre of the gas orbit. The radius and velocity are less affected as they depend on greater lengths in the channel map.

Acknowledgments

We thank T. J. Haworth, J. Lesur and C. McNally for the useful discussions. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00798.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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We acknowledge funding from ANR of France under contract number ANR-16-CE31-0013, from the Australian Research Council (ARC) under the Future Fellowship number FT170100040, from the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation, and from the European Union Seventh Framework Programme FP7-2011 under grant agreement No. 284405.

Appendix A: C18O observational results

The C18O channel maps are presented in Fig. A.1. In Fig. B.1, we compare positions-velocity diagrams along the disc major axis with curves representing Keplerian rotation in a geometrically thin disc with an inclination of 48°. The position–velocity diagrams are nicely reproduced by a central mass of 1 ± 0.1M⊙, assuming a distance of 161 pc. This is in agreement with the value used by Panić et al. (2009; rescaling for the difference in distance) and Cleeves et al. (2016).

Appendix B: Model description and parameters

We use the Monte Carlo radiative transfer code MCFOST (Pinte et al. 2006, 2009) to build a simple model of the CO emission. Our model is an extension of the model presented in Pinte et al. (2008), where the radial power-law surface density was replaced by a tapered-edge structure, following Panić et al. (2009). This model is only meant to provide a qualitative description of the ALMA band 6 data. A complete fit of all the observations of IM Lupi is beyond the scope of this paper, and we refer the readers to the work by Cleeves et al. (2016) for the most recent global model of IM Lupi.

The central object is described by a photosphere at 3900 K, a radius of 2.15 R⊙ and a mass of 1 M⊙. An ultra-violet excess is added to the phosphere, defined by a parameter fUV, which describes the fractional luminosity excess between 91.2 and 250 nm.

An external, isotropic, interstellar radiation field is added in model d presented in Fig. 5. It is composed of a highly diluted UV field, an infrared background radiation field, and the cosmic microwave background. We follow the same prescription as Woitke et al. (2016, see their Appendix A.4 for details). In this work, the UV component is scaled to 4 G°, following the results of Cleeves et al. (2016), who estimated the local radiation field from geometrical arguments based on Hipparcos data.

We consider an axisymmetric flared density structure with a Gaussian vertical profile. The disc scale height is described by a power-law (B.1)and its surface density by (B.2)where r is the radial coordinate in the equatorial plane and h0 is the scale height at the radius r0 = 100 au. The disc extends from rin, that we fixed to 0.3 au, to an outer rout, set to 8 times the critical radius rtap, i.e. where the density is so low that the actual value does not affect any of the synthetic observables.

Dust grains are defined as homogeneous and spherical particles with sizes distributed according to the power law dn(a) ∝ a-3.5da between amin = 0.03 μm and amax. Optical properties are computed using the Mie theory, adopting the same olivine stoichiometry from Dorschner et al. (1995) as in Pinte et al. (2008).

Position–velocity diagram of the continuum subtracted C18O line along the disc major axis. The Δv is relative to the systemic velocity of 4.5 km s-1. The solid line shows the midplane Keplerian velocity for a stellar mass of 1 M⊙ and an inclination of 48°. The dashed lines correspond to stellar masses of 0.9 and 1.1 M⊙. The data have been spectrally smoothed to a resolution of 80 m s-1. The spectro-spatial resolution is indicated in the lower left corner.

For the dust radial migration, we follow the simple prescription introduced in Pinte & Laibe (2014). The maximum expected grain size is set to the optimal size of migration (B.3)where ρg and cs are the gas density and sound speed, ρdust the intrinsic dust density, and ΩK the Keplerian frequency respectively.

CO chemistry is dominated by a few key mechanisms: photo-dissociation by UV radiation in the upper layers, condensation on the dust grains at low temperatures (T ≲ 20 K), and photo-desorption of the CO ice by UV radiation in the outer disc. We have implemented a very simple model of the CO chemistry, which was calibrated from ProDiMo models (Woitke et al. 2016, see e.g. their Fig. C.2). We initially assume that the CO abundance X(CO) is constant everywhere. CO is photo-dissociated, i.e. we set its gas abundance to 0 if (B.4)where n is the number density of hydrogen atoms and χ describes the intensity of the UV radiation field, following Woitke et al. (2009)(B.5)with FDraine = 1.921 × 1012 m-2 s-1, and where uλ is the energy density of the radiation field which was computed by MCFOST at each point in the model.

CO is frozen out on the dust grains if T< 21 K, i.e. its abundance is multiplied by a factor ϵ< 1, except if (B.6)where CO is photo-desorbed.

The dust temperature structure is computed using MCFOST. We then assume that the gas temperature is equal to the dust temperature and that the CO emission is at local thermodynamic equilibrium, and compute the synthetic channel map via a ray-tracing method.

As our aim was simply to build an illustrative model, we only adjusted the Band 6 CO data in the image place via a genetic algorithm (each generation was composed of 100 models and the genetic algorithm was run for 50 generations). As a consequence, our model may not be a unique solution, and the model parameters presented in Table B.1 must be interpreted with care.

Parameters of our disc model, as derived by our genetic algorithm from fitting of the Band 6 ALMA data only.

Appendix C: Continuum optical depth

To assess the continuum optical depth of the disc, a more accurate radiative transfer model of the continuum emission (Fig. C.1) is required. We use the same iterative procedure as in Pinte et al. (2016) for HL Tau. In short, we first extract a brightness profile along the disc major axis in the CLEANed map (average of both sides). We deconvolve this profile by a Gaussian corresponding to the size of the beam in that direction. We convert the deconvolved brightness profile into a surface density profile assuming an initial power-law radial temperature profile (exponent = −0.5), from which we compute a MCFOST model. We compare the resulting synthetic surface brightness profile with the observed one and iteratively correct the surface density until the synthetic brightness profile does not change by more than 1% at any radius. We reach convergence after about ten iterations. The resulting continuum optical depth as a function of radius is presented in Fig. C.2.

Continuum image of IM Lupi at λ = 1.3 mm. Contours begin at 4 σ and step in factors of 2 in intensity up to 1024 σ (green contour), where σ = 0.05 mJy/beam is the RMS measured on the map away from the source.

In the regions we are studying in this paper, i.e. outside 100 au, the maximum optical depth is ≈ 0.03, meaning that the flux on the lower surface of the disc is attenuated by at most 3%, due to continuum absorption by the dust.

Interestingly, the optical depth and corresponding dust surface density show two dips at ≈ 85 and 210 au, suggesting the presence of two gaps, and the corresponding two rings peaking at ≈ 115 and 240 au. They very likely correspond to the rings observed by Cleeves et al. (2016) at ≈150 and 250 au, after applying an unsharp mask.

All Tables

All Figures

Channel maps of 12CO (left) and 13CO (right) displayed with a channel width of 80 m/s. The data is not continuum subtracted. The dv is relative to the systemic velocity of 4.5 km s-1. The beam is indicated in the bottom right corner.

Single 12CO channel with a schematic of the quantities we measure. The disc has been rotated so that the semi-major axis is horizontal. The position of the central object is marked by a star. For a given offset δx along the disc major axis, two maxima N and F on the same vertical are located on the near and far sides of the disc’s upper surface. These two points are on the same inclined circular orbit (white ellipse), whose centre is marked by a diamond. Indicated are the geometrical parameters from which the cylindrical radius, height, and velocity of the points F and B can be measured. By repeating the procedure for each δx (green points, sampled every 0.08′′) and each velocity channel, a full mapping of the CO upper surface is performed. The cyan points mark the location of the lower surface (near side), while the magenta points represent the symmetric of the far upper surface relative to the disc major axis.

Measured altitude, velocity, and brightness temperature of the CO layers as a function of radius: 12CO (blue), 13CO (red), and C18O (green). Data points were extracted for all channel maps and binned according to distance. The error bars represent the dispersion within the bins added quadratically to the uncertainty resulting from inclination (48 ± 3°). Left: solid lines display a broken power-law fit (below and above r = 300 au) for the 12CO and 13CO. Centre: black solid line represents the Keplerian velocity for a 1 M⊙ central star, derived from the C18O position-velocity diagram (Appendix A). The shaded area represents an uncertainty of 0.1 M⊙. The red and blue dashed lines show the expected Keplerian velocities taking into account the altitude of the 12CO and 13CO layers. The red and blue solid lines display the velocity of the 12CO and 13CO layers taking into account the altitude and the pressure gradient term, as derived from our MCFOST model. Right: circles and triangles represent the temperature measured on the upper and lower surfaces, respectively. The brightness temperatures are similar on the upper and lower surfaces for 13CO and C18O; for clarity, only the upper surface is plotted.

Representative models of IM Lupi. The disc has been rotated so that the semi-major axis is horizontal. From left to right a): without any freeze-out of CO; b): with complete freeze-out of CO, where T< 21 K; c): with depletion of gaseous CO by 10-4, where T< 21 K; d): same as c) including the interstellar radiation field. The models are designed for qualitative comparison only. They are convolved to the ALMA resolution, but not degraded to the same noise level as the observations.

Comparison of the measured temperature (points) with the temperature structure from our MCFOST model (background). The same colour map is used for the data and the model. The white contour is the T = 21 K level, and the dashed contours correspond to T = 15, 30, 40, and 50 K.

Position–velocity diagram of the continuum subtracted C18O line along the disc major axis. The Δv is relative to the systemic velocity of 4.5 km s-1. The solid line shows the midplane Keplerian velocity for a stellar mass of 1 M⊙ and an inclination of 48°. The dashed lines correspond to stellar masses of 0.9 and 1.1 M⊙. The data have been spectrally smoothed to a resolution of 80 m s-1. The spectro-spatial resolution is indicated in the lower left corner.

Continuum image of IM Lupi at λ = 1.3 mm. Contours begin at 4 σ and step in factors of 2 in intensity up to 1024 σ (green contour), where σ = 0.05 mJy/beam is the RMS measured on the map away from the source.

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.