An ALMA/NOEMA study of gas dissipation and dust evolution in the 5 Myr-old HD 141569A hybrid disc

Context. The study of gas-rich debris discs is fundamental to characterising the transition between protoplanetary discs and debris discs. Aims. We determine the physical parameters of the brightest gas-rich debris disc orbiting HD 141569A. Methods. We analyse images from the NOrthern Extended Millimeter Array (NOEMA) 1 and the Atacama Large Millimeter/ submillimeter Array (ALMA) in 12 CO, 13 CO J =2 → 1, and 13 CO J =1 → 0 transitions. We incorporate ALMA archival data of the 12 CO J =3 → 2 transition and present continuum maps at 0.87, 1.3, and 2.8mm. We use simple parametric laws with the Diskﬁt code and MCMC exploration to characterise the gas disc parameters and report a ﬁrst attempt to characterise its chemical content with IRAM-30 m. Results. The continuum emission is equally shared between a compact ( (cid:46) 50au) and a smooth, extended dust component ( ∼ 350au). Large millimetre grains seem to dominate the inner regions, while the dust spectral index is marginally larger in the outer region. The 12 CO is optically thick, while 13 CO is optically thin with τ 13 CO ∼ 0 . 15 (C 18 O is not detected). The 13 CO surface density is constrained to be one order of magnitude smaller than around younger Herbig Ae stars, and we derive a gas mass M 12 CO = 10 − 1 M ⊕ . We conﬁrm the presence of a small CO cavity ( R CO = 17 ± 3au), and ﬁnd a possibly larger radius for the optically thin 13 CO J =2 → 1 transition (35 ± 5au). We show that the observed CO brightness asymmetry is coincident with the complex ring structures discovered with VLT/SPHERE in the inner 90au. The 12 CO temperature T 0 (100au) ∼ 30K is lower than expected for a Herbig A0 star, and could be indicative of subthermal excitation. Conclusions. With the largest amount of dust and gas among hybrid discs, HD 141569A shows coincident characteristics of both pro- toplanetary discs (central regions), and debris discs at large distance. Together with its morphological characteristics and young age, it appears to be a good candidate to witness the transient phase of gas dissipation, with an apparently large gas-to-dust ratio ( G / D > 100) favouring a faster evolution of dust grains.


Introduction
The evolution of circumstellar discs from the initial reservoir of the interstellar medium (ISM)-dust and gas mixture towards planetary systems is governed by the growth and coagulation of dust grains, viscous spreading and accretion onto the central star (or on proto-planets), and dissipation of the primordial gas (Williams & Cieza 2011). Although a large amount of high-resolution images and spectroscopic data have been collected during recent decades, the various stages of disc evolution remain poorly constrained by observations . The main accretion phase corresponds to the Class II phases for T Tauri (TTS) and Herbig Ae (HAe) stars. The final phase is represented by debris discs, when the primordial reservoir has been fully dissipated and/or incorporated into planetesimals Reduced datacubes and images are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/635/A94 and planets, and the detected emission from second-generation dust comes from collisions between the leftovers of planet formation or comet evaporation (Wyatt 2008;Hughes et al. 2018). Evidence for structural changes has been collected for a large number of pre-main sequence stars at various stages of the disc evolution (e.g. Andrews et al. 2018;Avenhaus et al. 2018;Garufi et al. 2018). This includes the presence of spirals, rings and gaps, central cavities, or azimuthal asymmetries, where different processes are at play such as the protoplanet(s)-disc gravitational or hydrodynamical interaction, viscous evolution, grain growth and opacity change, gas-grain hydrodynamical interaction, the development of vortices and magneto-rotational instabilities, and the photo-evaporation of gas. Transitional discs with large dust cavities (e.g. Espaillat et al. 2014) were found to represent almost 30% of the millimetre(mm)-bright discs in Taurus and Ophiuchus regions (Andrews et al. 2011). These systems show a lack of continuum mm emission with dust depletion factors in the range of 100-1000 in the disc central regions, while a variety of gas depletion factors have been reported within these cavities (e.g. van der Marel et al. 2016van der Marel et al. , 2018. They may witness gravitational interaction between gas and sub-stellar/planetary companion(s) in the disc cavity, and drag forces between dust grains and gas may cause the mm grains to be trapped by pressure bumps at large radial distances (e.g. Zhu et al. 2011;Pinilla et al. 2012).
While the timescale for disc dissipation has been constrained to about 5 Myr from the fading of infrared (IR)-excess (for dust) and of accretion tracers or molecular lines (for gas) in young clusters (Haisch et al. 2001;Fedele et al. 2010), very few systems are known in the transition phase between protoplanetary and debris discs. A new category of discs has emerged in recent years, which simultaneously display the typical characteristics of debris discs (with in particular a weak fractional excess of IR emission f d = L IR /L < 10 −2 originating from the optically thin dust component of second-generation grains), and the presence of gas tracers (mostly CO lines) that reveal a nonnegligible residual gaseous component (Zuckerman et al. 1995;Moór et al. 2011). Although only a handful of such gas-rich debris discs were known when we started the present project (namely HD 21997, 49 Ceti, β Pictoris), this population has now grown to almost twenty members with typical ages in the range of 10-40 Myr (e.g. Lieman-Sifry et al. 2016;Moór et al. 2017). The origin of the remnant gas is still a subject of controversy: the age of the system must be compared to the lifetime of the detected molecular species, which itself depends on the UV/ X-ray stellar/interstellar emission and on the (complex) shielding effects (Kral et al. 2017). Those systems harbouring supposedly primordial gas are called hybrid discs (Kóspál et al. 2013). The gas-rich debris discs revealed by CO lines simultaneously display a weak dust mm continuum and an IR excess characteristic of more evolved, optically thin debris discs. These systems are thought to witness a short-lived transient phase at the interface between proto-planetary and debris discs, when the material dissipates through photo-evaporation and accretion. In Péricaud et al. (2017), we reported a correlation between the CO flux density and the mm continuum emission for a broad variety of systems including T Tauri, Herbig Ae, transition discs, and even the proto-typical young debris disc β Pic. The position of hybrid discs in the diagram, systematically above the correlation line, suggests a faster evolution of the dust in these systems, leading to a transient and unusual increase of the gas-to-dust ratio.
In this population, the Herbig B9.5V/A0Ve star HD 141569A holds a very peculiar and interesting position. The star is located at 110.6 ± 0.5 pc (Gaia Collaboration 2018), and has two M-stars companions whose gravitational link remains unclear (Reche et al. 2009). With an age of about 5 ± 3 Myr (Weinberger et al. 2000;Merín et al. 2004), it is the youngest debris disc in this population with an IR brightness twice as large as that of β Pic. With the new Gaia distance, the stellar luminosity determined in Merín et al. (2004) turns to 27.0 ± 3.6 L . In IR colour-magnitude diagrams, HD 141569A lies at the exact interface between proto-planetary and debris discs (see, e.g. Wyatt et al. 2015, and Fig. 1 therein). Its cold molecular content was discovered through its CO emission line in the spectroscopic survey of Zuckerman et al. (1995), whereas warmer CO gas was later detected in the inner system by Goto et al. (2006). The debris disc of HD 141569A (L IR /L = 8.4 × 10 −3 , Sylvester et al. 1996) revealed by optical/IR dust emission presents a complex multiple-ring architecture with spirals, arcs, and gaps imaged in scattered light by the HST (Augereau et al. 1999;Weinberger et al. 1999;Mouillet et al. 2001;Boccaletti et al. 2003;Clampin et al. 2003;Konishi et al. 2016). More recently, a series of concentric rings was discovered with VLT/SPHERE with a range of radii of 40-90 au (Perrot et al. 2016), partially corroborating the emission detected around 40 au in the L band by Currie et al. (2016) with Keck/NIRC2.
While the small-grain component has been well studied in the near-infrared (NIR), a thorough analysis of the larger grains observed at mm wavelengths and a comparative physical and morphological analysis of the CO gas distribution is crucial to understanding its evolution status as well as the nature of the observed dust and gas components. Our group obtained observations of HD 141569A in 1999-2000 with the Plateau de Bure Interferometer (PdBI) of Institut de Radio-Astronomie Millimétrique (IRAM), where the gas disc was already partly resolved (briefly reported in Dutrey et al. 2004). These early observations have motivated the present study, but their quality and sensitivity is outmatched by the current dataset and will therefore not be presented here. Observations of this system by the Submilliter Array (SMA), Combined Array for Research in Millimeter-wave Astronomy (CARMA), and the Atacama Large Millimeter/submillimeter Array (ALMA) have been reported, revealing an extended gas disc (from the 12 CO J = 3→2 transition) out to the 250 au main dust ring highlighted in early HST images, and a dust disc extending out to 56 AU from the 870 µm continuum observations (White et al. 2016;Flaherty et al. 2016). The total mass of gas was revised to values much larger than initially derived thanks to the optically thin line 13 CO J = 2→1 (Miley et al. 2018). From the mm continuum emission, a double component with spectral indices consistent with the range of values observed in both debris and primordial discs was reported by , and an unresolved dust ring peaking at 220 au was detected from 1.3 mm ALMA observations by Miley et al. (2018). In this article, we present complementary, multifrequency observations with ALMA and the NOrthern Extended Millimeter Array (NOEMA) 1 (Sect. 2). We characterise the dust mm continuum and the gas emission, taking advantage of the spectral information to derive constraints on the grain growth and emission lines opacities (Sect. 3). We finally model the 12 CO and 13 CO lines with the radiative transfer code Diskfit (Piétu et al. 2007) to characterise the temperature and density laws in this disc (Sect. 4). In Sect. 5, we finally compare the derived morphology at mm wavelengths with the dust disc seen in optical/IR domain and use the new constraints on the gas and dust mass to discuss the origin of the gas and the evolution status of this peculiar system.

PdBI/NOEMA observations
HD 141569 was observed with NOEMA during three runs from December 2014 to April 2016. The 12 CO J = 2→1 transition at 230.538 GHz was mapped with the 7D and 7C configurations (7 antennas) of the IRAM array, with an angular resolution of 2.4 × 1.2 (baseline lengths from 15 to 192 m; all subsequent angular resolution values correspond to natural weighting). The spectral setup used a narrow band (20 MHz bandwidth, resolution 0.20 km s −1 ). Excellent conditions are noteworthy for one track, with precipitable water vapour inferior to 1.5 mm, and a 15 • mean phase rms.
In winter 2014/2015, we also obtained 13 CO J = 2→1 (220.399 GHz) observations. The 6C and 7D configuration antennas were used with a total on-source time of 12.7 h. We reached a 2.3 × 1.2 angular resolution (baseline lengths from 15 to 176 m) and used the same spectral resolution as for 12 CO J = 2→1. For both transitions, 1546+027 and 1508-055 were used for amplitude and phase calibration and we bootstrapped their fluxes using MWC 349 observations (with a flux model of 1.86 Jy at 219 GHz and 1.91 Jy at 230.5 GHz).
The WIDEX correlator provided 4 GHz bandwidth during the observations of the 12 CO and 13 CO. Combining the two observations at 1.3 mm, we obtained a 2.1 × 1.0 beam for the continuum map.

ALMA observations
We obtained ALMA observations of HD 141569 in Band 3 (84−116 GHz). The data were recorded in August 2015 for the 13 CO J = 1→0 (110.201 GHz) and C 18 O J = 1→0 (115.271 GHz) transitions during two runs with 39 and 44 antennas, respectively (project 2013.1.00883, PI. Péricaud). The spatial resolution was 0.76 × 0.56 (baselines up to 1.55 km), and the spectral resolution 0.2 km s −1 . The data were collected with relatively poor weather conditions (PWV∼ 4 mm). We used Ceres and 1550+054 for the flux and amplitude calibration, 1550+0527 and 1517-2422 for bandpass calibration and 1550+0527 for the phase calibration.
In addition, we collected archival raw data from the project 2012.1.00698.S (PI. Boley) to complement our analysis in Band 7 (275−373 GHz), which includes the 12 CO J = 3→2 line at 345.796 GHz. A preliminary analysis of this data set was reported in White et al. (2016). Only a subset of the observations could be properly calibrated. The observations carried out with 32 antennas provide a resolution of 0.43 × 0.35 (baselines up to 650 m), and a spectral resolution of 0.8 km s −1 . The flux and amplitude calibrator is Titan, the phase and bandpass calibrator is 1550+0527. We used the CASA package to calibrate the ALMA data, and we carried out the imaging and analysis of all NOEMA and ALMA observations with the GILDAS software. The imaging was performed with natural weighting for both ALMA and NOEMA.

A deep search for extra molecular lines with IRAM-30 m/EMIR
We also used the EMIR detector on the IRAM 30 m antenna to investigate the molecular content of this hybrid disc (Project 049-15). The observations were performed in Summer 2015 with a total time of 20 h on source. Two frequency bands were used in this project: 87.720-93.500 GHz and 140.520-148.300 GHz. We used Wobbler switching, resulting in very flat baselines over the observed bandwidths. The frequency coverage included lines of HCN J = 1−0, HCO + J = 1−0, and CS J = 3−2 that are usually bright in discs around classical T Tauri (of order 0.5-1.0 Jy km s −1 , see Dutrey et al. 1997) but slightly fainter in Herbig Ae stars (e.g. MWC 480 and AB Aur, Piétu et al. 2007;Schreyer et al. 2008). No emission was detected. Assuming a similar line width to that of the CO lines (7 km s −1 FWHM), the 3σ detection threshold on the integrated line flux is 45 mJy km.s −1 around 90 GHz, and 60 mJy km.s −1 around 145 GHz.

Continuum emission
The continuum emission is resolved at the three frequencies, although with some differences due to the different instrumental Zero-baseline flux F ν (u, v 0, mJy) 6.00 ± 0.25 4.00 ± 0.25 0.70± 0.05 Notes. F ν (u, v 0) corresponds to the maximum flux density measured at the shortest baseline for each frequency. ( * ) To compute the spectral index, we add quadratically an arbitrary uncertainty of 10% of the flux values to account for the instrumental calibration.
resolutions. The images presented in Fig.1 were produced using the CLEAN deconvolution algorithm with natural weighting. The emission appears compact, and centrally peaked, with a possible secondary component which looks more extended in the NOEMA map only (a low level plateau in the 1.3 mm image). The coarser resolution of the NOEMA map (231 × 110 au) is not the only reason: the visibility profiles presented in Fig. 1 clearly show a steep increase at the shortest baselines. This visibility rise at u, v distance shorter than 50-100 m (present at the three frequencies) is characteristic of extended emission (at few arcsec scales), while the more compact, inner disc component is resolved at intermediate baselines (e.g. B ∼ 150−400 m at 0.87 mm). The spatial extent and flux values of the continuum emission were determined by fitting Gaussian functions in the u, v plane. We tested a single and a two-component model. In Table 1 we present the best-fit parameters for the two models, and the measured maximum flux density corresponding to the shortest baseline points (so called zero-baseline flux). The single-component model is dominated by the compact disc component (0.3-0.5 ) at 0.87 and 2.8 mm, while at 1.3 mm only the broad component (2-3 ) is constrained. Therefore, the flux values for this model cannot be used together to estimate the dust spectral index. Our maximum flux values from the shortest baselines are consistent with those reported by single-dish observations or compact arrays (see Table 2). At 1.3 mm, we measure F 230 GHz (u, v 0) = 4.00 ± 0.25 mJy, in agreement with the JCMT/SCUBA measurement (5.4 ± 1.0 mJy, A94, page 3 of 16 A&A 635, A94 (2020) a c e Fig. 1. Left column: continuum emission maps, all contours are set with 3σ spacing, panel a: ALMA data at 0.87 mm, σ = 0.051 mJy beam −1 , contours; (c) NOEMA data at 1.3 mm, σ = 0.078 mJy beam −1 ; (e) ALMA data at 2.8 mm, σ = 0.008 mJy beam −1 . We note that the scale of panel c is changed to include the more extended emission. Right column: real part of the visibility (with spatial binning) as a function of the u, v radius in the Fourier plane, panels b, d, and f: correspond to data in panels a, c, and e, respectively.   Sylvester et al. 2001), but our value is about three times larger than the ALMA measurement of 1.7 mJy reported in Miley et al. (2018). Typical calibration uncertainties would amount to ∼10% for both instruments, but cannot justify this discrepancy. At 2.8 mm (and 0.87 mm resp.) our measurements are consistent with the CARMA (resp. SMA) values F 110 GHz = 0.78 ± 0.30 mJy (resp. F 230 GHz = 8.2 ± 2.4 mJy, Flaherty et al. 2016).
Because the inner component is very compact, the disc orientation and aspect ratio are only roughly reproduced by the elliptical Gaussian fit in the single component model (see Table 1). In the more realistic double component model below, we therefore first deprojected the visibilities assuming the geometric parameters derived at optical wavelength (PA= 356 • , i = 53 • , Augereau et al. 1999), and then fit in the u, v plane with a circular Gaussian function.
The second model suggests that the emission is about equally shared between a compact component, which is concentrated within about 50 au (FWHM int ∼ 30−35 au, a value comparable to the interferometric beams at 0.87 and 2.8 mm), and a much broader component that extends out to about 350 au. The u, v plane coverage does not allow precise determination of the outer radius, and most of the flux of this extended component is filtered by the interferometer at 0.87 and 2.8 mm. It is marginally detected as a low-level plateau in the 1.3 mm map thanks to the larger beam (although its shape and radial extent cannot be fully determined in our data). From higher-resolution ALMA data (0.65 beam) at 1.3 mm, Miley et al. (2018) revealed a lowlevel ring-like emission peaking at 220 ± 10 au that cannot be seen in our NOEMA map. This unresolved ring is different from the extended component that we report from our analysis of the u, v plane. In our higher-resolution maps at 0.87 and 2.8 mm, no such ring-like emission is detected. Assuming a spectral index α mm = 2 (see below), and extrapolating the 220 au ring peak brightness and our brightness sensitivity at 0.87 mm (0.4 beam) and 2.8 mm (0.66 beam), it would be marginally detectable at the ∼3 and 5 σ level, respectively. Such a two-component model was also proposed by  from u, v plane fitting, but our results disagree with their total flux values (e.g. F 345 GHz = 16.8 ± 1.7 mJy), which are systematically larger than both the zero-baseline flux values, and the earlier SMA, APEX, CARMA or JCMT observations (Flaherty et al. 2016;Nilsson et al. 2010;Sylvester et al. 2001). The observed differences between the flux values summarized in Table 2 can be explained firstly by the assumed shape of the model components (Gaussian, uniform discs, rings, power-laws), and the number of emission components taken into account. Secondly, the veryshort-baseline information is not always taken into account in the data analysis (e.g. White et al. 2016;Miley et al. 2018). Thirdly, the u, v coverage and the field of view can vary significantly between instruments, and the signal-to-noise ratio as a function of (u, v) radius can be very different from one dataset to another (e.g. at 1.3 mm).
We estimate the dust spectral index β from the slope of the spectral flux density, and we first add quadratically a typical calibration uncertainty of 10% to the flux density values. The spectral index is defined as β = α − 2, with α = log(F ν 1 /F ν 2 )/ log(ν 1 /ν 2 ), F ν . A fit to the three frequencies yields the slope α. We derive the spectral index β in two cases: (i) for the global disc using the maximum flux values measured at the shortest baselines (u, v) ∼ 0, (ii) for the double component model, which provides a more realistic description of the continuum emission than the single-component model. The values reported in Table 1 suggest that the spectral index is significantly smaller than the canonical ISM value with β ≤ 0.7 at the 3σ level. Whether this is indicative of grain growth or optically thick emission can only be determined with more resolved observations. In the disc modelling by Thi et al. (2014), the dust component is optically thin from optical up to mm wavelength both in the radial and vertical directions, and we therefore expect the constraints on the spectral index to primarily reflect a change in grain properties. Radial variations of this spectral index may also be expected, as already demonstrated in a sample of TTS and HAe stars (Guilloteau et al. 2011;Pérez et al. 2012). The double-component model suggests a larger value of β for the more extended disc component, although the difference with the inner 35 au region is not statistically significant. A similar trend has been claimed by , although their spectral index value for the outer component β = 0.28 ± 0.18 may be biased by an overly large flux value at 345 GHz, as explained above. A dense u, v coverage at both short and very long baselines would be required to firmly establish a change of dust properties between the inner compact and the extended dust disc components.

CO gas emission
The 12 CO and 13 CO J = 2→1 transitions are detected with NOEMA at a high signal-to-noise ratio (S /N 150, and S /N 20 respectively), whereas the 13 CO J = 1→0 detection with ALMA is more marginally significant (S /N ∼ 5). Integrated intensities for all CO transitions are presented in Table 3 with their associated spatial resolution. The intensity we derive for the ALMA 12 CO J = 3→2 observations 18.4 ± 0.3 Jy km s −1 is slightly larger than (but consistent within 2σ with) the value reported by White et al. (2016) and by Flaherty et al. (2016) with the SMA. The intensity maps (zeroth moment) are displayed in Fig. 2 for the three transitions (natural weighting). The 12 CO emission appears slightly more extended in the NOEMA map because of the convolution by a larger beam. The 13 CO J = 2→1 emission appears to be more compact than A94, page 5 of 16 A&A 635, A94 (2020)  the 12 CO J = 2→1 transition which extends beyond 2", likely as a result of a difference in opacity since the sensitivity and resolution are equivalent for these two lines. We also display the ALMA 12 CO J = 3→2 emission for comparison, using uniform weighting in order to show evidence for the tiny inner gas cavity which is readily resolved in this high-resolution image. Given its very low level, the continuum emission was not subtracted in Fig. 2, meaning that this CO cavity cannot artificially result from over-subtraction of the continuum. The separation of the central emission peaks is about 0.8 , which translates into a maximum cavity radius of about 45 au. This value is refined in the subsequent modelling presented in the following sections. The channel maps in Figs. A.1 and A.2 display the typical velocity pattern of rotating discs. In Fig. 3, the CO integrated area maps have been superimposed to the scattered light emission from the HST observations. These figures confirm that most of the molecular emission is apparently less extended than the optical dust disc. The 13 CO J = 1→0 transition has a rather poor  Table 3.

Opacity of CO lines
Using the line ratio of the 12 CO and 13 CO isotopologues for the same transition level, the line opacity can be estimated with the following formula: R 12/13 = S 12 CO S 13 CO = ν 12 ν 13 2 1 − e τ 12 where τ 12 and τ 13 are the opacity of the 12 CO and 13 CO lines, respectively. If we further assume that the opacity ratio of 12 CO over 13 CO is equal to the elemental abundance ratio (Wilson & Rood 1994), then τ 12 /τ 13 = 77 and we can solve for Eq. (1). For the J = 2→1 transition, we measure from NOEMA observations a line ratio R 12/13 = 7.1 ± 1.7, and derive τ 12 = 14 ± 4 and τ 13 = 0.18 ± 0.05. It should be noted that here we derived the CO opacity from the 12 CO / 13 CO flux ratio in velocity channels of 1 km s −1 to minimise the smearing effect due to rotation (Guilloteau et al. 2006) and averaged the results. The same analysis for the J = 1→0 transition, based on our new ALMA-Band 3 measurement and the integrated intensity S ( 13 CO) = 1.6 ± 0.2 Jy km s −1 reported by Flaherty et al. (2016) with CARMA yields τ 12 = 8 ± 2 and τ 13 = 0.11 ± 0.02, in good agreement with the J = 2→1 transition estimate. Even if these values could be biased by a factor of a few due to the uncertainties on the 12 C/ 13 C abundance ratio, as discussed in Kóspál et al. (2013), the 12 CO lines appears to be optically thick and the 13 CO optically thin. The 13 CO transition is therefore a more robust gas mass tracer than 12 CO, and this consideration leads us to reconsider the mass of gas derived by Flaherty et al. (2016) and White et al. (2016) who assumed optically thin emission for the 12 CO lines. Within the uncertainties, our upper limit on C 18 O at 115.271 GHz and 219.560 GHz is consistent with ISM CO isotopologue abundance ratios.

Modelling
We independently model the 12 CO J = 2→1, 13 CO J = 2→1, and 12 CO J = 3→2 lines with the Diskfit code (Piétu et al. 2007). We assume power laws for the CO surface density (Σ CO (r) = Σ 0 (r/100 au) −p ), the temperature T ex (r) = T 0 (r/100 au) −q , and the gas velocity V(r) = V 0 (r/100 au) −a v . In this study, we fixed the velocity exponent so that a v = 0.5 (Keplerian rotation). The disc is assumed to have a sharp outer edge at R out . The vertical density is assumed to have a Gaussian profile (see Eq. (1), Piétu et al. 2007), with the scale height being a (free) power law of radius H(r) = H 0 (r/100 au) −h . The line emission is computed assuming a (total) local line width dV independent of the radius and local thermodynamic equilibrium (LTE; i.e. T 0 represents the rotation temperature of the level population distribution). In addition to the above intrinsic parameters, the model also includes geometric parameters: the source position (x 0 ; y 0 ), the inclination i, and the position angle (PA) of the rotation axis ; and the source systemic velocity V sys relative to the local standard of rest (LSR) frame. The modelling is done in the u, v plane to avoid non-linear effects due to deconvolution. The minimization is completed by a Bayesian analysis using a Markov chain Monte Carlo (MCMC) code, emcee (Foreman-Mackey et al. 2013). Our approach differs from the analysis of White et al. (2016), which did not include the physical disc parameters in the modelling and focussed on the geometrical and kinematic description of the 12 CO J = 3→2 emission. Our results are presented in Table 4. Comparison with observations are shown for the integrated spectra in Fig. 4, and for the channel maps in Fig. 5.

Spatial distribution and kinematics of the CO gas
The spatial extent of the gas disc is constrained to consistent values around R out = 275 ± 3 au by the two transitions of 12 CO. This value is slightly larger (within 3σ) than the 224 ± 20 au reported by Flaherty et al. (2016) from a simultaneous fit to 12 CO J = 3→2 and 12 CO J = 1→0, because these latter authors fixed the density exponent to a smaller value: p = 1.5. The 13 CO line emission is slightly more compact with R out = 240 ± 20 au. Although the discrepancy does not appear to be statistically significant, it may result from the combination of opacity and sensitivity effects. More interestingly, an inner cavity is confirmed with a radius R in = 17 ± 3 au from 12 CO J = 3→2, as first suggested in lowerresolution maps by Flaherty et al. (2016). The spatial resolution and sensitivity of the NOEMA observations are not sufficient to provide a robust constraint, but point towards consistent values of the inner disc edge. The 13 CO J = 2→1 data suggest a larger value R in = 35 ± 5 au. This discrepancy should be taken with care because the beam size of these NOEMA data is more than four to five times larger than the ALMA beam for 12 CO J = 3→2. If confirmed, this result could be explained by the difference in opacity between the two isotopologues: since the 12 CO transition is 70 times more optically thick than the 13 CO transition, a larger inner radius for 13 CO could be indicative of a shallow increase of the density at the inner disc edge, while a sharp edge is expected to produce similar inner radii for both lines. A similar effect has been reported in the 15-70 au gas cavity of transitional disc J160421.7-213028 (Dong et al. 2017), although the nature of the two discs is different, because HD 141569A does not show any sign of having a dust cavity (a compact continuum emission was reported at all frequencies in Sect.3.1). At the reported resolution, the presence of unresolved gaps in the gas distribution in the central 50 au may also result in a similar discrepancy A94, page 7 of 16 A&A 635, A94 (2020) Fig. 4. Observed line profiles and superposed best models obtained with Diskfit (colour dashed lines) for, from left to right: 12 CO J = 2→1, 13 CO J = 2→1 NOEMA data, and 12 CO J = 3→2 ALMA data. between the CO isotopologues. The inner rim value (17 au) for 12 CO is coherent with the spectroscopic detection of warm CO through IR ro-vibrational lines by van der Plas et al. (2015), with a modelled distribution of the gas between 14 and 45 au (after correcting for the revised stellar distance). The inclination and position angle are found to be in the range i = 56−58 • and PA = 86 ± 1 • for 12 CO lines, and i = 53 ± 2 • and PA = 90 ± 2 • for 13 CO, which is in good agreement with the optical values of i = 52.5 ± 4.5 • and PA = 85.4 ± 1 • respectively (Augereau et al. 1999). The disc is flared, although the flaring exponent h is not well constrained, a value consistent with h = 1.25 is found for both 12 CO transitions, while it is naturally not constrained for the optically thin transition of 13 CO. For 12 CO, a scale height H 0 = 17 au is found at 100 au. This is twice the value of the hydrostatic scale height (H hydro = 8 au) for a disc with T 0 = 30 K, as measured from the optically thick transitions of 12 CO. This is consistent with direct measurements of the 12 CO layer location in HD 163296 (de Gregorio-Monsalvo et al. 2013) and the Flying Saucer (Guilloteau et al. 2016), where CO is located at 2−3 scale heights .
In our preliminary fits, the velocity exponent proved to be very close to the Keplerian value a v = 0.5, and we adopted this fixed value for the MCMC exploration reported in Table 4. At 100 au, the rotation velocity is better constrained for 12 CO lines, with a value V 0 = 4.44 ± 0.01 km s −1 that can be converted into a central stellar mass of 2.22 ± 0.01 M . This value is slightly smaller than the mass of 2.39 ± 0.05 M reported by White et al. (2016), but is very consistent with the spectral analysis of Merín et al. (2004) that yielded M = 2.00 +0.18 −0.15 M . The 13 CO analysis leads to a slightly larger mass M = 2.35 ± 0.05 M .

Temperature and density laws
The modelling of the two optically thick transitions 12 CO J = 2→1 and 12 CO J = 3→2 provides very consistent estimates for the excitation temperature, T 0 ∼ 30 K with a relatively flat radial profile q < 0.4 (q = 0.28 ± 0.01 for 12 CO J = 3→2 and q = 0.10 ± 0.09 for 12 CO J = 2→1). The temperature is ill constrained for the optically thin 13 CO transition. A value of 30 K is consistent with the analysis of the 12 CO J = 3→2 detection with SMA by Flaherty et al. (2016, who used a fixed temperature exponent q = 0.5 and found T 0 = 33 +11 −4 K). It is somewhat smaller than the gas temperature of about 50 K derived from 12 CO line for a similar, albeit younger and more optically thick disc around the Herbig Ae star MWC 480 (Piétu et al. 2007). Given the star luminosity of 30 L , blackbody grains in an optically thin disc are predicted to have a temperature of 60 K at a distance of 100 au, and would be much hotter than the detected gas. The derived gas temperature is also significantly smaller than the range of kinetic temperatures computed by Thi et al. (2014) with their Prodimo modelling based on unresolved detections of [OI], C[II] and 12 CO J = 3→2 lines around HD 141569. In their modelling, the gas is heated by Polycyclic Aromatic Hydrocarbons (PAH) and cooled by [OI] and CO lines, and the computed temperature reaches about 100 K at 100 au in the inner disc and drops down to 40 K at 500 au (i.e. well beyond the outer edge of the emission we detect with ALMA/NOEMA). The excitation temperature we  Table 4), and the residual emission. The latter is obtained by subtraction in the u, v plane before imaging and cleaning. Contour spacing is 5σ for the 12 CO transitions (σ = 16 mJy beam −1 for J = 2→1 and σ = 6.3 mJy beam −1 for J = 3→2), and 3σ for the 13 CO J = 2→1 transition (σ = 10 mJy beam −1 ). Residual contour spacing is 3σ for all transitions.
derive from 12 CO may still be smaller than the actual gas kinetic temperature. Low values of the gas temperature (although with a much more dramatic deviation) have also been reported in the hybrid disc HD 21997 by Kóspál et al. (2013), where the excitation temperature derived from CO line ratios drops to values as low as 6−9 K. Flaherty et al. (2016) propose that this low gas temperature could be explained by enhanced cooling from CO (following Hollenbach & Tielens 1999), and/or NLTE effects (following Matrà et al. 2015).
In the absence of C 18 O emission, we use the more optically thin emission of the 13 CO J = 2→1 line to constrain the CO surface density. Our modelling yields a best-fit value Σ 0 = 5 +2 −1 × 10 15 cm −2 at 100 au with an exponent p = 1.8 ± 0.5. We note that the two other transitions of 12 CO lead to much steeper A94, page 9 of 16 A&A 635, A94 (2020) Fig. 6. Global scheme of dust and gas distribution in HD 141569 summarizing the most recent resolved observations. profiles with p ∼ 3 and require a surface density of the order of 6 × 10 16 cm −2 . These values of the surface density are at least one order of magnitude smaller than those found for the younger, gas-rich disc around the Herbig Ae star MWC 480: Piétu et al. (2007) report, at a distance of 100 au, Σ 0 = 8 × 10 16 cm −2 for 13 CO and 6 × 10 18 cm −2 for 12 CO. This system is taken here as a reference as it is one of the very few where the surface density of 13 CO has been precisely constrained through similar disc emission modelling.
We derived a total mass of CO gas based on Σ 0 and used the minimum outer disc edge value R out = 232 au inferred from the 13 CO line (hence a minimum mass value). Because most of the mass lies in the inner regions with a power law distribution, the computed value is sensitive to the inner disc edge, and we take into account the inner cavity detected with the 12 CO emission in the ALMA high-resolution maps (R in = 17 au). Assuming a standard isotopic ratio 12 CO/ 13 CO = 77, these parameters yield M12 CO = 0.10 +0.05 −0.02 M ⊕ (=4 × 10 23 kg), where we take into account the uncertainties on R in , R out , Σ 0 and p. This value is 50 times larger than the CO mass reported by White et al. (2016) based on the 12 CO emission.
The conversion to a total mass of gas is very uncertain. A low CO content compared to the dust emission was already noted in the early detection of CO isotopologues in DM Tau (Dutrey et al. 1997). Several recent studies (Bergin et al. 2013;Reboussin et al. 2015) favour much lower values for the [CO/H 2 ] ratio than the canonical ISM value of 10 −4 (see also the discussion on carbon depletion in e.g. Miotello et al. 2017). We discuss this point again in Sect. 5.2. Here, in order to compare with other works, we nevertheless convert this minimum mass of CO into M H 2 = 71 +36 −14 M ⊕ or 0.21 +0.11 −0.04 M Jup 2.10 −4 M using the standard ISM abundance. Using the grid of models developed by Miotello et al. (2016), which include the effect of isotope selective photodissociation, the luminosity of the 13 CO J = 2→1 line and the non-detection of C 18 O J = 2→1 favours a total mass in the range 1−5 × 10 −4 M , in good agreement with the previous estimate. This is one order of magnitude larger than the value reported by Flaherty et al. (2016) from the SMA detection of 12 CO J = 3→2, but is in good agreement with the earlier modelling of Thi et al. (2014) and Jonkheid et al. (2006), based on spatially unresolved observations of the same line, who found a total gas mass in the range 67−164 M ⊕ . Our value is a factor of three to five smaller than the total gas mass of 6−9 × 10 −4 M reported by Miley et al. (2018) based on the integrated flux of the optically thin transition 13 CO J = 2→1, but these remain in reasonable agreement due to the different methodologies used. This total gas mass value is clearly on the lower bound of the range of disc masses found for Herbig Ae stars, and usually derived from mm dust flux converted to gas mass with an assumed gas-todust ratio of 100 (see, e.g. Williams & Cieza 2011). It is instead typical of discs around very low-mass TTS (for instance, TTS disc masses obtained from 13 CO and C 18 O J = 3→2 in Lupus range between 10 −5 and 10 −3 M , based on models assuming isotope selective photo-dissociation and CO freeze-out on grains Miotello et al. 2017). HD 141569 also harbours one of the most massive molecular gas components among the gas-rich debris discs reported in Moór et al. (2017) based on supposedly optically thin CO isotopologues. Even if we underestimate the CO mass due to 13 CO optical depth effects or the uncertain [CO/H 2 ] abundance, such a low mass supports the process of ongoing gas dissipation in this 5 Myr-old disc.

Comparison with optical/NIR emission
The distribution of matter is strikingly different at optical/NIR and mm wavelengths. A global sketch summarizing and comparing the respective radial extent of the main disc components detected at optical/NIR and mm wavelengths is presented in Fig. 6. The continuum emission at mm wavelength is clearly much more compact than the 250-400 au broad rings of the debris disc revealed at optical/NIR wavelengths, and is also more compact and centrally peaked than the gas emission characterised in the following section. Furthermore, it is coincident with the innermost component recently reported by various authors in the NIR domain. Indeed, in the inner 100 au recent observations with VLT/SPHERE have revealed the presence of a series of concentric rings and arcs with semi-major axes in the range 45−90 au (Perrot et al. 2016). The link between the µm grains responsible for the IR emission and the mm grains detected with ALMA/NOEMA cannot be clarified because of the limited resolution of our data. VLT/VISIR observations in the mid-IR also suggest the presence of an even more internal component (Thi et al. 2014), which was recently detected at L band with the Keck vortex coronagraph within 70 au (Mawet et al. 2017) and may extend up to 100 au. According to the extended scattered-light emission reported with HST/STIS by Konishi et al. (2016), this inner warm component may also be spatially associated with the inner rim of the gas disc that our modelling locates around 15−20 au from 12 CO J = 3→2 maps (similarly to the warm CO detected in spectroscopy in the 14−60 au region by Goto et al. 2006;van der Plas et al. 2015). The disc is not classified as a transitional disc because of the lack of a mm dust cavity. A small cavity of no more than 15 au in radius is observed in the CO maps, consistent with the earlier spectroscopic detections of a warm CO ring, but the continuum emission is strikingly centrally peaked. Inside this gas cavity, there is probably not much material: VLTI-PIONIER observations in the H-band reported a completely unresolved structure (Lazareff et al. 2017), which might be dominated by the stellar emission itself (accounting for 86% of the total flux), without suggestion of an inner ring near the sublimation distance of small grains.  The scattered light images show many asymmetries at various spatial scales, with spirals, arcs, and multiple ring structures. The outer limit of the gas disc is precisely coincident with the so-called "inner ring" near 250 au reported previously. Among the possible ring formation mechanisms, Takeuchi & Artymowicz (2001) proposed a complex mechanism based on the outward migration of small grains at the outer edge of a gas disc. HD 141569 is the only resolved disc so far that displays this combination of distant dust rings and a gas-rich inner disc.
Large-scale asymmetries are also suggested in both wavelength domains. In Fig. 7, we superimpose the SPHERE/ Ks image on the map of our best-fit residual of 12 CO J = 3→2 ALMA emission. This residual emission, which was already reported by White et al. (2016), is significant (up to 6σ) at all velocity channels and peaks between 6 and 9 km s −1 (see Fig. 5). We present here the integrated emission, which extends on the western side with an elliptical shape. It displays two brightness peaks centred at 60-70 au of deprojected distance at azimuth angles of about 220 • and 290 • , in the region of the three concentric rings detected with SPHERE. Since the 12 CO emission is optically thick and the feature is not detected in the isotopologues due to insufficient sensitivity and angular resolution, we cannot disentangle between a local temperature offset and a density feature. However, the significance of this 12 CO brightness increase, the coincidence with the location of the inner rings, and the similar brightness asymmetry raise questions about the interplay between gas and dust in this region. Lyra & Kuchner (2013) proposed a mechanism based on the development of a clumping instability to create eccentric rings from the coupling of gas and grains, but this appears to be efficient in a regime of low gas-to-dust ratio, which does not seem to be favoured in the inner region (see discussion in Sect. 5.2). The 250 au dust ring (and possibly also the outermost 400 au ring) at the edge of the detectable gaseous component (see Fig. 3) may instead be an ideal location to reach the necessary very low value of gas-to-dust ratio (G/D < 1) for such an instability to develop. So far, no planetary companion candidate has been reported in the inner 100 au region with a sensitivity limit of the order of 1−3 M Jup (Perrot et al. 2016). The destruction of planetesimals (by collision or evaporation) could produce both large quantities of small grains and release substantial amounts of gas, as has been suggested by the detection of an asymmetrical emission of CO in β Pictoris (Dent et al. 2014). In HD 141569, such a collisional event or enhanced planetesimal evaporation might be more difficult to demonstrate because of the large residual amount of primordial gas in the region of the NIR rings.

Gas-to-dust ratio
The global gas-to-dust mass ratio in the disc is expected to trace the evolution of the system, but getting an accurate value is challenging because of the numerous assumptions that must be made especially in the computation of the total gas mass, as discussed in the previous section. The mass of dust is no less uncertain as it relies on the unknown value of the dust mass opacity κ ν , and we will follow the prescription: κ ν = 0.02(λ/1mm) −β cm 2 g −1 , which is similar to that adopted by Beckwith et al. (1990) (with β = 1). M dust can be estimated on the basis of the mm dust flux (F ν ), which is constrained in our observations at three different frequencies. From the modelling of the SED, earlier models by Jonkheid et al. (2006) and Thi et al. (2014) showed that the dust emission is optically thin from optical to mm wavelengths, meaning that the mm flux can be used to provide a robust estimate of the dust mass. Using the mm dust fluxes reported in Table 1, the extreme values of β in the range [0,0.7], and typical values of the temperature for dust grains located at 50-100 au from this 30 L star (T d = 65−90 K), M dust is computed with the standard equation: We obtain values for M dust in the range 0.03−0.52 M ⊕ . This value is significantly smaller than for other Herbig Ae systems. In order to compare with dust masses derived for other gas-rich debris discs reported in Moór et al. (2017), we alternatively used the same dust opacity (κ ν=230 GHz = 2.3 cm 2 g −1 ) and their assumed temperature (T d = 50 K), which leads to M d = 0.46 M ⊕ , which makes it the most massive gas-rich debris disc with HD 131488 in this study. With our gas mass estimate, 70 M ⊕ , and our first choice of dust opacity, yielding M dust ∼ 0.03−0.52 M ⊕ , we can obtain a gas-to-dust mass ratio G/D in the range 135−2370. This is most likely a severe lower limit, since we assumed a high CO/H 2 abundance ratio to derive the gas mass. In T Tauri discs, it is in general necessary to invoke an overall depletion of CO by a factor of ten or so to bring the derived G/D ratio to the ISM value of 100 (see, Piétu et al. 2007;Williams & Best 2014), although the warmer discs around HAeBe stars have higher CO/dust ratios, and can sometimes be directly interpreted by G/D = 100 and a normal CO/H 2 ratio (e.g. AB Aur Piétu et al. 2007, or also HD 163296 in Williams & Best 2014. The estimated gas-to-dust mass ratio appears at least one order of magnitude larger than the values reported in Miotello et al. (2017) for a sample of sources in Lupus, where gas mass estimates were also derived from optically thin CO isotopologues its extremely low continuum which suggests a faster evolution of the dust component as already highlighted in Péricaud et al. (2017).