Asymmetric mid-plane gas in ALMA images of HD~100546

TL;DR

This study presents new ALMA observations of HD 100546's transitional disc, revealing azimuthal asymmetries, mass estimates, and grain growth evidence, advancing understanding of planet-forming environments.

Contribution

First detection of C$^{18}$O in HD 100546's disc and detailed analysis of asymmetries, mass, and grain growth in a proto-planet hosting disc.

Findings

Azimuthal asymmetry in gas and dust regions.

Mass estimates of dust and gas in the disc.

Evidence of grain growth to cm scales inward.

Abstract

In this paper we present new ALMA observations towards the proto-planet hosting transitional disc of Herbig Ae/Be star HD 100546. This includes resolved 1.3 mm continuum, $^{13}$CO and the first detection of C$^{18}$O in this disc, which displays azimuthal asymmetry in regions spatially coincident with structures previously identified in HST images related to spiral arms. The lower limit on the mass of the dust disc is calculated to be 9.6x10$^{-4}$M$_\odot$. A firm lower-limit on the total gas mass calculated from optically thin, mid-plane tracing C$^{18}$O (2-1) emission is 0.018M$_\odot$ assuming ISM abundances. These mass estimates provide an estimate of gas-to-dust ratio in the disc of 19, the ratio will increase if C$^{18}$O is relatively under-abundant in the disc compared to CO and H2. Through deprojection and azimuthal averaging of the image plane we detect 1.3 mm continuum emission out to 290+/-10 au,$^{13}$CO to 390+/-10 au and C$^{18}$O to 300+/-10au. We measure a radially increasing millimetre spectral index between wavelengths of 867$\mu$m and 1.3 mm, which shows that grain sizes increase towards the star, with solid particles growing to cm scales in the inner disc.

Figures

Tables

Table 1. Table 1: Sexadecimal positions of central emission in the 1.3 mm continuum and CO isotopologues as determined by fitting a Gaussian to integrated emission maps in the image plane (‘Im.’ above), and a uniform disc in the uv plane (‘uv’). Positions are rounded to the nearest milliarcsecond. The angular separation between these positions and the GAIA position of HD 100546, corrected for proper motion between GAIA and ALMA observations, is given in the final column.

	Coordinates of best	Offset from Stellar
	fitting model	Position ("), P.A.
Cont. (Im.)	11:33:25.305 -70.11.41.258	0.082, 79∘
13CO (Im.)	11:33:25.310 -70.11.41.291	0.073, 49∘
C18O (Im.)	11:33:25.305 -70.11.41.300	0.098, 55∘
Cont. (uv)	11:33:25.305 -70.11.41.257	0.082, 80∘
13CO (uv)	11:33:25.308 -70.11.41.276	0.072, 63∘
C18O (uv)	11:33:25.300 -70.11.41.280	0.110, 71∘

Equations

$$F_{\nu }=\frac{{\kappa }_{\nu }{B}_{\nu }(T)M_{dust}}{D^2},$$

$$M_{\text{gas}}=\frac{4\pi }{h{\nu }_{21}}\frac{F_{21}m{d}^2}{A_{21}{x}_2},$$

$$I=\frac{h\nu A_{J,J-1}N(J)}{4\pi },$$

$$N_{J}=N(total)\frac{2J+1}{Z}\frac{1}{exp\big{(}\frac{hB_{e}J(J+1)}{kT}\big{)}},$$

$$Z=\sum_{J=0}^{\infty}(2J+1)~{}exp\Big{(}\frac{-hB_{e}J(J+1)}{kT}\Big{)},$$

$$I=\frac{h\nu A_{J,J-1}}{4\pi}\frac{N_{\rm total}~{}(2J+1)}{exp\Big{(}\frac{hB_{e}J(J+1)}{kT}\Big{)}}\frac{1}{\sum_{J=0}^{\infty}(2J+1)~{}exp\Big{(}\frac{-hB_{e}J(J+1)}{kT}\Big{)}}.$$

$$\frac{S{(}^{13}\mathrm{CO})}{S({\mathrm{C}}^{18}\mathrm{O})}=\frac{{\nu }_{{}^{13}\mathrm{CO}}}{{\nu }_{{\mathrm{C}}^{18}\mathrm{O}}}\frac{T_{ex{,}^{13}\mathrm{CO}}}{T_{ex,{\mathrm{C}}^{18}\mathrm{O}}}\frac{1-e^{{\tau }_{{}^{13}\mathrm{CO}}}}{1-e^{{\tau }_{{\mathrm{C}}^{18}\mathrm{O}}}}.$$

$${\beta }_{\mathrm{mm}}(R)=\frac{\mathrm{dlog}{\text{F}}_{\nu }\text{(R)}-\mathrm{dlog}{\text{B}}_{\nu }}{\mathrm{dlog}\nu }.$$

Full text

††pubyear: 2018

1 Introduction

Interferometric observations at millimetre and sub-millimetre wavelengths often reveal millimetre-sized particles in the outer regions of Herbig and T-Tauri discs. Radial drift timescales for particles in the outer regions of the disc ($>$100 au) are of order $10^{5}$years, however disc lifetimes are of order $10^{6}-10^{7}$ years, meaning the particles in the outer disc should have already drifted onto the star. Retention of large grains in the outer disc, or alternatively some form of continuous replenishment, therefore seems to be essential in order to explain the observational data (Klahr & Bodenheimer, 2006; Birnstiel & Andrews, 2013). The consequences of mechanisms that concentrate disc material locally are now becoming frequently observed. Results from ALMA have revealed highly asymmetric, so called ‘lop-sided discs’ such as SAO 206462 and SR 21 (Perez et al., 2014), IRS 48 (van der Marel et al., 2015) and HD 142527 (Casassus et al., 2015). Radial gas pressure bumps can also result in dust traps that form symmetric, concentric rings in the dust disc (Pinilla et al., 2012b; Dullemond et al., 2018).

Continuum emission at mm-wavelengths is also a diagnostic tool for tracing the dust in a protoplanetary disc, specifically the coolest dust contained in the mid-plane and in the outer disc, where the majority of solid mass is expected to be contained. Measured flux can be related to the mass of emitting particles by assuming optically thin emission and using

[TABLE]

where Fν is the continuum flux measured at frequency $\nu$, $\kappa_{\nu}$ is the specific opacity and D is the distance to the source from the observer. Continuum emission also depends on particle size; millimetre continuum emission comes most efficiently from a narrow range of grain sizes between 1mm and 1cm (Takeuchi & Lin, 2005). If grains have grown beyond $a_{max}\approx 3\lambda$, where $\lambda$ is the wavelength of observation, they emit less efficiently (Draine, 2006). Radial drift and grain growth predict the largest solid particles to be concentrated toward the star, leaving the smaller ones behind. As emission brightness is related to grain size we can therefore expect a variation in observed disc size at different observing wavelengths, but also a variation in spectral index, $\alpha_{\rm mm}$, across the disc. Spatially resolved continuum images of discs at multiple wavelengths allow for a calculation of a position-dependent $\alpha_{\rm mm}$, where, for optically thin emission, S${}_{\nu}\propto\nu^{\alpha_{\rm mm}}$ (e.g. HD 163296, Guidi et al. 2016). Surveys of star forming regions find that disc-averaged spectral indices can lie in the range 1.5 $\leq\alpha_{\rm mm}\leq$ 3.5 (Ricci et al., 2012; Pinilla et al., 2014). If the size of grains in the disc varies with radius, then so will the spectral index.

Draine (2006) computed the frequency dependent opacity of different grain size distributions. A calculation of the opacity index $\beta$, where $\kappa_{\nu}\propto\nu^{\beta}$, between two spatially resolved observations at different frequency can therefore provide a quantitative estimate as to how large grains have grown in different parts of the disc by a comparison to the calculations of Draine (2006). $\beta$ values evolve with grain growth; $\beta\approx$ 1.5-2 is associated with the ISM (Mathis et al., 1977), and for the vast majority of disc masses, these ‘lower’ values of $\beta\leq 1$ are likely to be as a result of grain growth rather than other potential explanations such as optically thick emission (Ricci et al., 2012).

In order to calculate a reliable total mass of the disc we must have a robust gas mass measurement (Williams et al., 2014) as well as one for the dust. The mass of the dust disc is canonically related to the gas mass by an empirical ratio taken from the ISM; g/d (gas-to-dust ratio) = 100. Modelling of gas emission has increasingly suggested this may not be the case however. A complication is that modelling of the gaseous spectral lines is generally degenerately sensitive to g/d - that there is no unique g/d that satisfies the observations (Bruderer et al., 2012; Boneberg et al., 2016). Using ALMA data, a recent survey of the Lupus star forming region gives evidence that supports low observed g/d in protoplanetary discs; Miotello et al. (2017) find that 23 out of their 34 discs with gas and dust detections have g/d < 10. Chemical reprocessing can leave CO depleted in discs, with CO locked into larger icy bodies (Du et al., 2015) or converted into complex molecules (Kama et al., 2016; Yu et al., 2016). The carbon chemistry that controls the abundances is highly sensitive to the temperature profiles of the dust within the disc (Reboussin et al., 2015).

A complication for inferring gas properties such as the total mass is that gas emission is typically either optically thick or requires high sensitivity to detect. In this study we observe C18O, a less optically thick alternative to other CO isotopologues (Van Zadelhoff et al., 2001; Dartois et al., 2003). The trade-off here is that high sensitivity is required to detect a molecule much less abundant than other gas tracers. For example in the ALMA survey of Lupus, Ansdell et al. (2016) detect C18O (3-2) emission in 11 out of their 89 targets with integrations of between 30-60 seconds per source. Detection rates are higher for intermediate mass Herbig stars, where the brighter star supports a warmer disk meaning more C18O is found in the gas phase. To date there is only a handful of C18O detections in the discs of Herbig stars, including MWC 480 (Guilloteau et al., 2013), HD 169142 (Panić et al., 2008; Fedele et al., 2017), HD 163296 (Rosenfeld et al., 2013) and AB Aur (Pacheco-Vazquez et al., 2016). When the optically thin C18O and other isotopologues are observed it permits us to compute the CO disc mass and temperature. Emission from C18O traces the dense disc mid-plane (Van Zadelhoff et al., 2001; Dartois et al., 2003), where we expect grain growth and the formation of giant planets. Resolved detections of emission from this molecule will therefore allow us to study a part of the disc that is crucial for planet formation.

Much recent work has focused on transition discs. Transition discs are a class of protoplanetary discs that have a lesser flux in the mid-IR due to the opening of a gap or cavity in the inner disc which is generally explained by some combination of mechanisms such as grain-growth (Tanaka et al., 2005), photo evaporation (Owen, 2016; Ercolano & Pascucci, 2017), clearing by a forming exoplanet or binary companion (Baruteau et al., 2014; Price et al., 2018) and dead-zones (Dzyurkevich et al., 2010). This clearing eventually leads to gas-poor, optically thin, debris discs (Espaillat et al., 2014; Wyatt et al., 2015). There is now an increasing number of observational examples of resolved inner cavities of transition discs (Williams & Cieza, 2011; van der Marel, 2017; Pinilla et al., 2018). Characterising the distribution of gas and dust in such systems may be able to shed light on how this transition occurs.

In this paper we present new ALMA observations towards Herbig Ae/Be star HD 100546 in order to characterise the dust and gas in its disc. In the next section we briefly review recent and related published work on HD 100546. In Section 3 we describe the data taken and then present the initial results and measurements in Section 4. In Section 5.1 we calculate the gas and dust mass in the disc and assess apparent asymmetry in the interferometric images produced. We use the deprojected, azimuthally averaged flux profiles to determine the radial extent of the emission for each of the observed tracers and calculate $\beta_{\rm mm}$ to infer levels of grain growth. Our main results are summarised in Section 6.

2 HD 100546

This study focuses on a fascinating example of planet formation; the well-studied disc around Herbig Ae/Be star HD 100546. The star is located at a distance of 110pc (Lindegren et al., 2018) with a mass reported by van den Ancker et al. (1997) of 2.4 M*⊙* from pre-main sequence (PMS) track fitting and more recently 1.9 M*⊙* by Fairlamb et al. (2015) derived through X-shooter spectroscopy and PARSEC PMS tracks (Bressan et al., 2012). Although previously thought to be as old as $\geq$10Myr (van den Ancker et al., 1997), Fairlamb et al. (2015) calculate an age for the star of 7.02$\pm$1.49 Myr and using the updated parallax from GAIA data release 2, Vioque et al. (2018) find a younger age of 5.5${}^{+1.4}_{-0.7}$ Myr.

2.1 Evidence for proto-planet(s)

The system has been long associated with active planet formation and so is well studied. A protoplanet has been confirmed at $\sim$50 au via direct imaging (Quanz et al., 2015) and there has been some evidence of an inner companion as well. This includes SED modelling and mid-IR interferometry (Mulders et al., 2013; Panić et al., 2014) , ALMA observations (Walsh et al., 2014) , near IR spectroscopy (Brittain et al., 2013, 2014) and a point source in GEMINI Planet Imager (GPI) data at 10 au (Currie et al., 2015) support the potential HD 100546‘c’. However this predicted companion has alluded detection in more recent GPI and Magellan Adaptive Optics System observations by Follette et al. (2017). The position of the potential inner companion coincides with the position of the rim of the inner gap at $\sim 11\pm 1$ au, interior to which is a small inner dust disc of $\leq$0.7 au (Panić et al., 2014) and a recently detected bar-like structure reaching across the inner gap in both continuum and H$\alpha$ emission (Mendigutía et al., 2017).

2.2 Structure in the disc

Observing with ATCA, 7mm continuum emission can be fit by a Gaussian with FWHM between 50 and 60 au (Wright et al., 2015), and the authors assert that grains in the disc are well-processed and have grown as large as 5cm and possibly further. Walsh et al. (2014) are able to detect a significant amount of weak, extended continuum emission beyond this and find a best-fit model for the dust in the uv plane comprised of two rings with a Gaussian profile; one centred at 26 au with a FWHM of 21 au and one at 190 au with FWHM 75 au. Pineda et al. (2014) find that mm-sized grains are trapped within a similar inner ring and note brightness asymmetries in the dust between the South-East and North-West of the disc. All three of these studies at (sub-)millimetre wavelengths favour arguments of at least one giant planet orbiting within the disc. The disc does not appear to be axi-symmetric; scattered light studies routinely observe multiple asymmetric features and spiral arms within the dust disc (Garufi et al., 2016; Follette et al., 2017), including a ‘dark lane’ observed by Avenhaus et al. (2014) in the same direction as the horse-shoe shaped asymmetry identified in 7 mm images of (Wright et al., 2015). The exact nature of this asymmetry at millimetre wavelengths has until now remained largely undefined.

The gas in the disc is much more spatially extended than the dust, CO has been measured out to $\sim$400au (Walsh et al., 2014). Very little is known about the gas mass of this uniquely interesting disc due to a lack of data observing optically thin gas lines. Gas mass calculations based on optically thick 12CO emission can lead to underestimations of as large as two orders of magnitude in comparison to those using optically thin CO isotopologues (Moór et al., 2011; Miley et al., 2018).

3 Observations & Data Reduction

We present new ALMA observations of HD 100546 with band 6 receivers across 4 unique spectral windows (SPWs); project ID 2015.1.01600.S, PI: O. Panić. The two continuum windows, SPW0 and SPW1, both contain 128 channels of width 15625 kHz to give a total bandwidth of 2 GHz, SPW0 has a central frequency of 232 GHz, with SPW1 at 234 GHz. SPW2 contains 1920 channels of width 61.035 kHz, which at central frequency 220398 MHz gives resolution of 0.083kms*-1* for the window containing the 13CO (2-1) emission line. The C18O (2-1) line is in SPW3, this window has 3840 channels of 0.166 kms*-1* at central frequency 219560 MHz. Observations were made on the 28th March 2016 for a total time of 625s on source, employing 43 antennas with a baseline range of 15-460 metres giving a maximum angular scale of 10$\aas@@fstack{\prime\prime}$6. The flux calibrator was J1107-4449 and the phase calibrator was J1136-6827, the latter is found at the position 11:36:02.09787 -68:27:05821, which is (0$\aas@@fstack{\prime\prime}$22, 1$\aas@@fstack{\prime\prime}$73) degrees offset from the target.

Calibration and reduction was carried out using casa version 4.7.0. and imaging was performed using casa tasks and the clean algorithm (Rau & Cornwell, 2011). Line free continuum channels were used to construct image models to be used for self-calibration, the solutions were then applied to all channels. After applying self-calibration, images were created with briggs weighting and robust=0.5, giving a synthesized beam of 1$\aas@@fstack{\prime\prime}$05 x 0$\aas@@fstack{\prime\prime}$77, P.A. = -20*∘* for the continuum image. Source emission is contained within the inner third of the primary beam and so primary beam correction was not applied to images. Data cubes were created for the line emission by subtracting the continuum from SPWs containing the line emission using casa task uvcontsub. Images of the lines have beam sizes of 1$\aas@@fstack{\prime\prime}$08 x 0$\aas@@fstack{\prime\prime}$82, P.A. = -19*∘* for 13CO and 1$\aas@@fstack{\prime\prime}$10 x 0$\aas@@fstack{\prime\prime}$83, P.A. = -20*∘* for C18O. Integrated intensity, from which fluxes are measured, and intensity weighted velocity maps were created with casa task immoments by summing all emission in channels containing significant (>3$\sigma$) signal; in the case of 13CO this corresponds to a range between -0.62 and 11.4 kms*-1* and for C18O between -2.0 and 10.8 kms*-1*.

We also utilise archival cycle 0 ALMA observations of HD 100546 in band 7 (Project ID: 2011.0.00863.S, PI: C. Walsh), as described by Walsh et al. (2014). We re-reduce and image this data, applying self calibration as described above, to obtain continuum emission observed at a wavelength of 867$\mu$m. The cycle 0 data used J1147-6753 as a phase calibrator, offset from the primary target by (1$\aas@@fstack{\prime\prime}$2, 2$\aas@@fstack{\prime\prime}$3). The continuum image has a synthesized beam of dimensions 0$\aas@@fstack{\prime\prime}$95 x 0$\aas@@fstack{\prime\prime}$42, 38*∘*. The minimum baseline of 21 m results in a maximum angular scale of $5\aas@@fstack{\prime\prime}$3.

4 Results

The integrated continuum flux in the 1.3 mm image of 492 mJy is measured with an image rms of 0.5 mJy/beam. The 13CO (2-1) integrated flux is measured at 12.872 Jy kms*-1* and the first detection of C18O in the disc gives an integrated flux of 2.948 Jy kms*-1* for the (2-1) transition, with rms noise 0.033 Jy/beam kms*-1* and 0.027 Jy/beam kms*-1* respectively. The 2$\sigma$ contour reaches a maximum separation of 2$\aas@@fstack{\prime\prime}$5 in the 1.3 mm image, 3$\aas@@fstack{\prime\prime}$6 in 13CO and 3$\aas@@fstack{\prime\prime}$0 in C18O. We discuss asymmetry and the radial extent of the emission maps thoroughly in Section 5.2.

To find a position angle for the source we use the data cubes from the CO isotopologue emission to identify the minor axis of the disc. We find the minor axis by locating the channel at the systemic velocity of the source, identified by the symmetry of the map as here it will show no Doppler shifting. Channel maps show the minor axis of 13CO in the channel at 5.51$\pm$0.08 kms*-1* and for C18O at 5.33$\pm$0.17 kms*-1*. We adopt an average of these two velocities as $v_{\rm sys}$ from this data set, giving 5.42 kms*-1*, shown in Figure 2 along with the minor axis channel of each isotopologue. This average is consistent with the two $v_{\rm sys}$ determined from the channel maps within the uncertainties (which we assume to be represented by the channel widths). Figure 3 plots the observed spectrum of each of the two lines. Both lines display a double peaked spectrum, where each of the C18O peaks is broader than the peaks of 13CO. The mid-point of the 13CO spectrum is consistent with the average $v_{\rm sys}$. The mid-point of the C18O spectrum is harder to determine due to its shape. The local minimum at the centre of the double peak profile is offset from the adopted $v_{\rm sys}$ by 0.33 kms*-1*. Walsh et al. (2014) constrain a system velocity of 5.7 kms*-1* from their ALMA observations of 12CO.

The orientation of the disc minor axis is used to determine the position angle of the disc. We use our constrained position angle as an a priori estimate for a Gaussian fit of the data, in order to determine the inclination of the source. Any asymmetry in the emission map will influence the fitting of a Gaussian, by using this observationally determined a priori estimate for position angle, we can ensure our fitted inclination is more accurate.

For 13CO this method gives a position angle of the major axis, measured East of North, of $142\pm 0.8^{\circ}$ and for C18O, $146\pm 2.3^{\circ}$. The average value of the two is $144\pm 2.4^{\circ}$ which we adopt henceforth. Uncertainties give the variation in position angle when measured from the adjacent spectral channels. Keeping this position angle fixed, we then fit an elliptical Gaussian to each of the CO isotopologue integrated emission maps to determine inclinations of $42.6^{\circ}$ and $43.7^{\circ}$ for 13CO and C18O respectively. We adopt an average inclination of $43^{\circ}$ for all further analysis.

These values are consistent with most reported position angles from previous sub-mm observations of the disc with ALMA; Walsh et al. (2014) find a position angle of 146${}^{\circ}\pm 4^{\circ}$ and inclination 44${}^{\circ}\pm 3^{\circ}$, whilst Pineda et al. (2014) find 145.14${}^{\circ}\pm 0.04^{\circ}$ and 41.94${}^{\circ}\pm 0.03^{\circ}$. Panić et al. (2014) observed the disc with the MID-infrared Interferometric instrument (MIDI) on the VLT and find a position angle of 145$\pm 5^{\circ}$, their inclination is considerably higher however because their observations probe the inner wall of the disc at 11 au. There is growing evidence for a warp in the disc (Quillen, 2006; Panić et al., 2010; Walsh et al., 2017) which would cause inclination to deviate radially.

Asymmetry in the 12CO line profile is detected in both APEX and ALMA observations (Panić & Hogerheijde, 2009; Walsh et al., 2017). Images of the ALMA observations of 12CO achieve a synthesizing beam of 0$\aas@@fstack{\prime\prime}92$x0$\aas@@fstack{\prime\prime}$38 (Pineda et al., 2014). We do not detect this asymmetry in the less abundant 13CO or C18O lines, finding instead double peaked line profiles symmetric to within uncertainties (Figure 3). 12CO emission traces higher scale heights in the disc, meaning that the asymmetry could be restricted to the upper vertical layers of the disc. There already exists evidence for asymmetry of the disc surface in the scattered light images that identify spiral arms and other structures (Ardila et al., 2007; Avenhaus et al., 2014; Garufi et al., 2016; Follette et al., 2017). The APEX observations have a larger beam and therefore the results may be affected by the extended envelope of the disc as detected by scattered light as far as $\sim 1000$ au (Grady et al., 2001; Ardila et al., 2007).

We adjusted the GAIA data release 2 (2018) position of HD 100546 for proper motion between the dates of the ALMA and GAIA observations, giving a position of 11h33m25.3209s -70d11m41.2432s, within an accuracy of (0.04, 0.05) mas (Lindegren et al., 2018). We find that the centre of emission for each of the tracers we observe is offset from this position.

Table 1 gives the coordinates of the centre of emission for the continuum and CO isotopologues, accompanied by its offset from the GAIA coordinates. A 2D Gaussian is fitted to the emission maps in Figure 1, the positions of these fits are taken to be the centre of emission. We also undertake a similar analysis, this time fitting a uniform, elliptical disc geometry in the uv plane using CASA task uvcontsub. We repeat this for each tracer, giving the position of the best-fit disc model in Table 1. The benefit of uv-analysis is that it avoids being misled by spurious artefacts that can be introduced into images during the image reconstruction process. In each case the uv-fitted position is very similar to that determined in the image plane. The uniform-disc fits give a position angle of 146*∘* and inclination of 42*∘*, values which are within a few degrees of those we constrain above, and those from previous ALMA observations (Pineda et al., 2014; Walsh et al., 2014).

The uncertainty in an astrometric position in an ALMA image depends on a number of factors, which we address here. The signal-to-noise based target position-fitting error is 1.3 mas. This is consistent with the ALMA technical handbook Section 10.6.6 . The dominating factor in positional uncertainty is introduced when the phase corrections are transferred from the phase calibrator to the target. The phase calibrator is 1.7 degrees in separation from HD 100546, and from the phase drift between scans we estimate an total uncertainty on the astrometric position of $\sim$10 mas. This is relatively low as a result of a stable phase in the observations. Other factors such as antenna position errors can degrade astrometry; the ALMA technical handbook cautions that with reasonable phase stability, positional accuracy is $\sim$1/20 of resolution, which in this case would correspond to a maximum of 50 mas. We adopt this greater value of 50 mas as a conservative uncertainty on coordinate positions in the image. Each of the best-fit centres of emission in Table 1 have a separation from the central star that exceeds the 50 mas uncertainty in ALMA image positions. An offset of emission from the stellar position may indicate an over-density of disc material towards one side of the disc. We discuss asymmetry in the disc in more detail in Section 5.2.1.

5 Analysis and Discussion

5.1 Disc Mass

5.1.1 Dust Mass

We measure the flux within a 3$\sigma$ contour of the 1.3 mm continuum image. The opacity depends on the maximum grain size, for which we only have loose constraints. A minimum dust mass estimation can be made by assuming a dust grain size distribution $dn/da\propto a^{-3.5}$ (Mathis et al., 1977). Assuming a maximum grain size $a_{\rm max}$ = 1 mm yields $\kappa_{1.3\rm mm}$ = 1.15 $\rm cm^{2}g^{-1}$ (Draine, 2006). Mid-plane temperature of the dust is difficult to constrain from observations, however we can get a firm handle on the mid-plane disc mass by assuming a sensible temperature range for this dense, cool environment, we therefore adopt a range of 20K – 40K. In this transition disc, gas and dust are likely to be well coupled, and our strong gas line detections suggest that the disc temperature must be over the CO freeze-out temperature of $\sim$20 K. Herbigs host more luminous, warmer discs than their T Tauri counterparts, and so the temperature is likely to be above this 20 K lower bound. We adopt an upper bound of 40 K. In the radiative transfer models of A type stars by Panić & Min (2017), gas mass in the disc must be increased by over two orders of magnitude to increase mid-plane disc temperatures from 20 to 40 K at short separations from the star, and so we consider this to be a generous upper limit.

Using Equation 1 with the assumptions above, this gives a minimum total dust mass of 4.1$\times 10^{-4}~{}$M*⊙* for $T$= 40 K, and 9.6$\times 10^{-4}~{}$M*⊙* for $T$= 20 K. In Figure 4 we show how this dust mass can vary with mid-plane temperature or if the assumed level of grain growth varies, thereby changing $a_{\rm max}$ and the emissivity of the dust. Between $a_{\rm max}$= 1 mm and 10 mm, the most likely values in a protoplanetary disc at this stage of evolution, the masses remain almost constant and are less sensitive to temperature. This indicates that for reasonable assumptions of $a_{\rm max}$ and mid-plane temperature, the dust mass is $\approx$10${}^{-3}~{}$M*⊙*.

The mass range calculated is comparable with the estimate of Wright et al. (2015) who employ the accretion disc models of d ’Alessio et al. (2005) to find Mdust of order 1.0$\times 10^{-3}~{}$M*⊙. Repeating the calculation above using instead the ALMA measured flux at 867$\mu$m from Walsh et al. (2014) with a temperature of 30K and $a_{\rm max}$ = 1 mm, a disc mass of 0.7$\times 10^{-3}~{}$M⊙* is found, similar to the results between 1mm and 10mm in Figure 4.

This is more massive than most discs around low-mass stars, which fall in the range $10^{-6}-10^{-4}\rm M_{\odot}$ (Manara et al., 2018). This dust mass is also at the higher end of the range found for Herbig Ae stars $10^{-4}<M_{\rm dust}(M_{\odot})<10^{-3}$ (Panić at al., in prep). Given the age of the star this is quite exceptional, suggesting a prime object to study longevity of dust in discs.

5.1.2 Gas Mass

From the integrated line flux of C18O above 3$\sigma$, assuming optically thin emission and local thermodynamic equilibrium, we calculate the total gas mass similarly to Hughes et al. (2008); Matrà et al. (2015); White et al. (2016),

[TABLE]

where F21 is the integrated line flux, $\nu_{21}$ is the rest frequency of the transition, d the distance to the source, m is the mass of the CO molecule, $A_{21}$ is the appropriate Einstein coefficient, and $x$ is the fractional population of the upper level. We adopt ISM fractional abundances for CO and its isotopologues, specifically CO/H${}_{2}\sim 10^{-4}$, and then adopt the relative fraction of C18O in the ISM, CO/C${}^{18}\rm O\sim$ 550 (Wilson & Rood, 1994). A gas mass is calculated in this way over a range of temperatures. The minimum gas mass was calculated with a temperature of 20 K, giving a minimum disc gas mass of 0.018 M*⊙*. This estimated gas mass is around 1% of the stellar mass, corresponding to around 18 Jupiter masses. For comparison, in an ALMA survey of Lupus, most low-mass discs have a gas disc mass of $\leq 10^{-3}\rm M_{\odot}$. The gas mass found here is large even for Herbig stars, where the majority have discs within the range (1-10)$\times 10^{-3}M_{\odot}$ when calculated from mid-plane tracing CO isotopologues (Panić at al., in prep). HD100546 therefore appears to still possess enough disc mass to form Jupiter mass planets, a very significant amount for a disc that has reached the transition disc stage of its evolution and shows detections of active planet formation (Quanz et al., 2015).

The calculated gas mass assumes ISM abundances of C18O relative to CO. Isotope-selective photodissociation can alter these ratios significantly in some discs, which can lead to underestimates of disc mass by up to an order of magnitude (Miotello et al., 2014). The CO line luminosities we measure in this data are unable to conclusively determine the extent of any isotope-selective photodissociation when compared to the the models of (Miotello et al., 2016), and so our calculation represents a lower limit. Carbon freeze-out beyond the snowline can also result in low measured gas masses; T Tauri discs appear to have significant levels of carbon freeze-out, for example in TW Hya (Schwarz et al., 2016). Herbig discs such as HD 100546 on the other hand remain warmer for a greater radial extent due to the more massive host star (Panić & Min, 2017). This pushes the CO snowline to a larger radial separation allowing a larger amount of carbon in the disc to remain in the gas phase. In the disc of HD 100546 we do not expect freeze-out to play a significant role. Fitting physical-chemical models to observations, Kama et al. (2016) find that within the disc of TW Hya, C and O are strongly under-abundant, whereas in HD 100546, the depletion of gas-phase carbon appears quite low.

In light of the discussion above, our calculated gas estimate represents a firm lower limit, but a very large one in comparison to other Herbig discs. The gas mass estimated here is similar to that achieved through line modelling of another bright Herbig disc, HD 163296 (Boneberg et al., 2016), and shows there is still a large budget remaining in the disc for giant planet building in comparison to discs that are in a more advanced stage of their evolution towards a debris disc, e.g. HD 141569 (White et al., 2016; Miley et al., 2018).

5.1.3 Gas To Dust Ratio

Taking the dust mass calculated with $T$=20K and the minimum gas mass calculated at the same temperature, we calculate a g/d for HD 100546 of $\approx$ 19. Low g/d (< 100) is observed in the majority of sources in recent surveys of discs in star forming regions (Ansdell et al., 2016; Long et al., 2017). Studies of individual Herbig discs suggest g/d is generally low in these larger, brighter discs too (Meeus et al., 2010; Boneberg et al., 2016). The g/d affects disc evolution because gas dynamics dominate the behaviour of the dust in gas-rich protoplanetary discs. The uncoupling of these two populations when gas densities drop sufficiently low is a signpost for the transition to a debris disc.

Low observed g/d can mean a great amount of gas has been lost, or it can be the result of carbon depletion in the disc through sequestration into large icy bodies or the conversion of CO into more complex molecules. This ratio is indicative only of regions within the CO snowline where freeze out is avoided. In the similar Herbig disc HD 163296 the snowline is calculated to be at around 90 au (Boneberg et al., 2016), in cooler discs the snowline will be at a lesser separation, e.g. TW Hya, 40 au (Schwarz et al., 2016).

5.2 Disc Asymmetry

5.2.1 Evidence in ALMA Images

Asymmetry in the disc has previously been reported with evidence in asymmetric 12CO line profiles (Panić & Hogerheijde, 2009), horse-shoe structures in cm-wavelength observations (Wright et al., 2015) and the identification of multiple scattered light features (e.g. Ardila et al., 2007; Garufi et al., 2016; Follette et al., 2017). Asymmetric features can form in the disc as a result of mechanisms such as planet-disc interactions (Baruteau et al., 2014), interactions between binary companions (Price et al., 2018) or particle traps within the disc (van der Marel et al., 2013; Van Der Marel et al., 2016). Such features are therefore particularly interesting in a disc such as that of HD 100546 which is known to host an embedded protoplanet, and exhibits a range of observed structures across different spatial scales.

In the data presented here we see a striking match with the extended emission in the C18O map and the structures in HST scattered light observations of the disc. The low-level bulges in C18O emission (Figure 5) persist in images created using different weighting parameters in the image reconstruction process, and are present both before and after applying self-calibration. Using the Advanced Camera for Surveys (ACS) on board the Hubble Space Telescope (HST) Ardila et al. (2007) identify 3 structures in the outer disc which are attributed to the known spiral arms in the disc. The authors name these structures 1, 2 and 3 and we adopt this convention here as well. Figure 5 displays the HST image in the F435W band presented by (Ardila et al., 2007). Each of the low-level bulges in the extended C18O (2-1) emission is found in the same region of the disc as one of the 3 scattered light structures. This provides tentative evidence of a mid-plane gas equivalent to the scattered light features. The C18O emission is optically thin, and can therefore be interpreted as an increase in local mid-plane density. An increased local gas density will increase the scale height of the disc, the subsequent elevation of small particles towards the scattering surface will result in a stronger scattered light flux, as we see in the HST data. The C18O emission bulges towards the outer disc may therefore be signalling the mid-plane density features that correspond with the disc’s spiral arms. In this section we make a general analysis of any observational evidence for asymmetry in the disc, within the limits of the angular resolution of our ALMA observations.

In Section 4 we accurately measure the inclination and position angle, both of which are very close to values derived from previous sub-millimetre ALMA observations of the disc (Pineda et al., 2014; Walsh et al., 2014). We take this geometry to construct an ellipse with a semi-major axis equivalent to 100 au in the disc. We can make an initial visual assessment of asymmetry in the disc by plotting this ellipse over the emission maps in Figure 1 centred at the stellar position. The mid-plane tracing 1.3 mm continuum and C18O (2-1) emission in particular show evidnece on an offset disc in Figure 6, on scales larger than that the offsets listed in Table 1. In both panels we omit high level contours internal to the inner ellipse for clarity. Towards the northwest in the panels of Figure 6, the 75$\sigma$ contour (in the 1.3 mm image) and the 15$\sigma$ contour (C18O) overlap with the 100 au ellipse, whereas to the southeast there is a significant offset. C18O in particular shows azimuthal asymmetry in the outer disc, with low level emission extending much outward as previously highlighted in Figure 5.

The largest asymmetry in the disc is along the minor axis (Figure 6) with the brighter region found towards the position of Structure 2 in Figure 5. We can also see this in Figure 7, which shows slices from the emission maps taken along the minor axis in both directions, centered at the stellar position.

C18O (bottom panel) is brighter in the southwest for the majority of the disc, and emission extends further by $\sim$0$\aas@@fstack{\prime\prime}$3 in comparison to that in the northeast. In the continuum (top panel) the southwest side is brighter than the northeast throughout the disc. The greatest difference between the two sides is at an angular separation of 0$\aas@@fstack{\prime\prime}$7, where the southwest is brighter by 50$\sigma$. A low-level bump in the continuum emission just beyond 2$\aas@@fstack{\prime\prime}$0 occurs at the same radial location as the extended C18O emission and scattered light emission in Figure 5. 13CO emission (middle panel) is also brighter in the southwest for as far as $\sim 2\aas@@fstack{\prime\prime}$8, beyond which the northeast of the disc is brighter.

5.2.2 Comparison to previous asymmetries in HD 100456

Figure 5 displays how the low-level, extended emission towards the edge of the gas disc is spatially coincident with structures identified by Ardila et al. (2007), potentially linking observed scattered light features with the mid-plane gas above which the small, scattering grains are suspended. In other scattered light observations, Avenhaus et al. (2014) see a dark lane that extends from $\sim 0\aas@@fstack{\prime\prime}2-0\aas@@fstack{\prime\prime}6$ along the minor axis. The authors attribute this to effects due to the polarised scattering function, rather than any physical feature in the disc. However Wright et al. (2015) also observe an arc of millimetre emission in the same direction, at the same separation as the scattered light dark lane, but not quite as extended.

Mid-IR interferometry is able to resolve the inner disc, where asymmetry is found close to the inner gap of the disc at 9 au (Panić et al., 2014). Features in the the inner disc can have a strong effect on the outer disc. For example there has been recent evidence of warps from ALMA observations, e.g. this disc (Walsh et al., 2017) and T Tauri star AA Tau (Loomis et al., 2017). A warp of the inner disc can lead not only to shadowing of the outer disc but potentially even a break of the disc into two distinct planes. Substructure can then develop as a result of this geometry; a misaligned disc can result in shadows that create asymmetric illumination of the disc, whilst a ‘broken’ disc can contain dust streamers that cross gaps in the disc (Dutrey et al., 2014; Loomis et al., 2017). Recent very high angular resolution (0$\aas@@fstack{\prime\prime}05$) ALMA observations of the disc resolve the very inner regions of the disc, revealing an asymmetric circumstellar ring in the continuum, brightest towards the North East (Pineda et al., 2018). The authors attribute this asymmetry to another scattered light spiral feature, this time identified by SPHERE (Garufi et al., 2016), on angular scales much less than that of the synthesized beam achieved in the observations presented here.

Many of these features can arise through interactions with a central companion as shown by 3D hydrodyanamic simulations of the transition disc HD 142527 (Price et al., 2018), this is particularly interesting in regard to any asymmetries found in the HD 100546 disc when bearing in mind the evidence found in recent years for a further companion at smaller separation.

5.3 Asymmetry Analysis

Asymmetry in the surface brightness of the gas emission can occur from increased temperature, increased density or a combination of the two.

In 3D hydrodynamic simulations, spiral arms launched by planets observed at 1.3 mm have been shown to induce azimuthal variation in gas surface density of approximately 50-60% (Juhász & Rosotti, 2018), and should also be observable in optically thin gas observations. In this scenario, a large spiral arm could increase gas density locally, thereby elevating the scale height of micron-sized grains to produce the scattered light features. Alternatively a local density increase of 30% over a length scale$\approx$ disc scale height is enough to provide the conditions for dust trapping of larger particles in the outer disc (Pinilla et al., 2012a). Another explanation for local increase in density is a gravitationally unstable disc, where clumping occurs due to gravitational forces if the free fall time scale is short compared to sound speed crossing time scales and shearing time scales (Toomre, 1964). This is generally found in discs at an earlier stage in their evolution whilst discs are still very massive $(\rm M_{disc}/M_{\star}>10\%$), which is not the case here.

Brightness asymmetry can also result from varying temperature in the disc. Azimuthal variation in disc temperature may result from an inner warp that obscures sections of the disc or azimuthal variation of density.

5.3.1 Density Scenario

In this scenario we assume that the difference in observed brightness is due solely to changes in gas density. We split each half of the disc, bisected along the major axis, into three 60∘ degree segments and sum the flux in the segment. By calculating the brightness ratio between two segments on opposite sides of the disc, we can calculate the increase in surface density required to produce the difference in flux. The intensity of integrated line emission is given by

[TABLE]

where AJ,J-1 is the Einstein coefficient for spontaneous emission from level J to level J-1, and N(J) is the column density in a specific level. The intensity is therefore proportional to the density of molecules in the N(J) level.

The largest brightness difference is between the two segments along the minor axis; the southwest is brighter than the northeast by a factor of 1.10$\pm$0.15 . The asymmetry identified in Section 5.2.1 is most obvious in the outer disc. We also compare the brightness ratio of the outer disc, we do this in the same fashion, except only summing the emission beyond a deprojected radius of 100 au. In this case, the brightness ratio is 1.27$\pm$0.18 . We are limited to comparisons that average over regions of the disc due to the angular resolution of our observations, but local density increases on smaller scales can alert us to certain physical mechanisms in the disc, some of which we discuss below.

Given the 1$\aas@@fstack{\prime\prime}$0 resolution of the data, it is feasible that there are much larger azimuthal density variations for a given radial separation that are diluted by the beam. This might arise due to a local increase in density due to, for example, a spiral arm. Ardila et al. (2007) associate the three structures in Figure 5 with spiral arms, and it is towards each of these structures that the outer disc of C18O emission seems least symmetric.

From this data we cannot rule out the possibility of an unresolved vortex capable of trapping particles at large separation in the HD 100546 disc, as our density variation is similar to that predicted by Pinilla et al. (2012a). High resolution data across multiple wavelengths would help to identify any such feature (e.g. Cazzoletti et al., 2018).

Calculation of the Toomre parameter, Q, based on the radial profile of C18O rules out gravitational instability in the disc. Assuming C18O traces total disc mass, Q(R) for this disc reaches a minimum $\sim 24$ at $\sim$80 au confirming that the disc is stable against gravitational collapse. However a clump in the outer disc will experience relatively long orbital timescales, meaning it will take longer to break up. Fragments of the disc that were previously unstable might linger at the outer edges of the disc as they dissipate.

5.3.2 Temperature Scenario

To test the possibility of temperature variation as the cause of asymmetry in the outer C18O disc, we measure the average flux in different segments of the disc as in the previous section, and then calculate the temperature variation that would be required to reproduce such a flux difference. Through a comparison of the average fluxes in opposite segments, we can calculate the increase in temperature required if this was to be due to temperature alone (i.e. no density increase).

We calculate the column density of C18O gas in the J=2 level using

[TABLE]

where T is the midplane temperature and Be is the rotation constant. The partition function, Z, is calculated using

[TABLE]

Substituting for Z and N(J), Equation 3 becomes

[TABLE]

Assuming a range of mid-plane temperatures, we then calculate $\Delta$T required to produce the observed brightness ratio of 1.27 from Section 5.3.1. We plot the results against the range of assumed mid-plane temperatures in Figure 8.

Figure 8 shows that even for a mid-plane temperature as cool as 20K (below which CO is no longer in the gas phase), a 75% increase in average temperature is required. From 24K and above, the average temperature of the outer disc mid-plane must increase by greater than 200%. There are no known physical processes expected within the disc mid-plane which could produce such a large change in temperature, especially when averaged over such a large proportion of the disc. If we were to assume that inner parts of the disc dominate the flux in each segment, temperatures may be larger than in the cool outer disc. However, the radiative transfer models of Panić & Min (2017) show that in reality an increase in disc temperature from 20 to 32K (comparable to the 20-35K increase in Figure 8), requires an increase in gas density of order 10, far greater than what we detect here. We can therefore rule out temperature as a sole cause of any flux asymmetry. In order to explain the difference in brightness between the two regions of the disc we must invoke some sort of density increase as well.

5.4 Radial Flux Profiles

Ignoring azimuthal variations in disc emission, we can consider how the flux changes in the radial direction by de-projecting and then azimuthally averaging, around the mm continuum centre, the emission from images in Figure 1. The resultant average profiles are shown in Figure 9.

Taking advantage of the high signal to noise in the data we adopt the following method to obtain an accurate uncertainty estimate for the radial flux profiles.

Each image is averaged off-source in ‘empty’ sky locations at 45*∘* position angle intervals at a separation that does not include any significant emission from the central source, but is also well within the primary beam. The resultant radial profile gauges the level of fluctuations due to background noise in that part of the image (shown as grey lines in Figure 9, more visible in the insets). The amplitudes of these fluctuations as a function of radius are fitted with polynomials, the polynomials are then averaged over the different off-source locations. The resulting final polynomial gives the average level of off-source fluctuations in the image as a result of the azimuthal averaging process. We use this for error bars on the average flux curve. Our method also allows for a determination of average outermost detection radius that does not rely on a pre-defined function or model, we take this to be the point at which the average flux becomes less than 2 times the average noise fluctuations at that radius; in Figure 9 this corresponds to the error bars reaching the grey shaded area.

Figure 9 plots the azimuthally averaged continuum flux and integrated line emission against deprojected radius (using inclination=43*∘* and PA=144*∘* as determined in Section 4). The 1.3 mm continuum for as far as an average, deprojected radial separation of 286$\pm$11 au. The 13CO profile shows us that the gas in the disc is significantly more extended; detected as far out as 385$\pm$11 au, a comparable distance to the 390$\pm$20 au extent of 12CO in the disc reported by Walsh et al. (2014).

The less abundant molecule C18O is detected out to a shorter separation of 297$\pm$11 au. Found at lower scale-heights, this molecule traces the cool, dense mid-plane Van Zadelhoff et al. (2001); Dartois et al. (2003) in which giant planets are expected to form.

A large amount of mass spread over large radial separations is important for planet formation in the disc. Giant planets grow more massive in discs of higher initial disc (gas) mass (Mordasini et al., 2012). Furthermore, population synthesis models show that significant amounts of giant planets migrate into the central star after formation (See review Mordasini, 2018). In order to survive, some cores accrete late in the disc lifetime and grow as they migrate inward (Baruteau et al., 2014). Our results suggest the disc of HD 100546 is well placed to support such late survivors, because a large amount of planet building material ($\sim$18 MJ) still exists over a large region of the disc (297 au outermost detection of C18O). This is particularly relevant in the case of this disc, where a confirmed proto-planet lies towards the outer disc at a radial separation of $\sim$50 au (Quanz et al., 2015).

The difference in the outer radius at which the discs are detected may be a true physical difference, or could be as a result of the different sensitivity between the continuum and the line emission. We explore this by calculating surface density profiles for each tracer, which are plotted in Figure 10.

In Figure 10 we plot the surface density profile of each tracer in order to assess how it changes with radial separation. The 13CO curve represents a lower bound because the emission is optically thick for most of the disc. The continuum surface density profile is scaled by an assumption that g/d in the disc is the canonical value of 100, giving an upper bound with which we can compare to the gas. The profiles in Figure 10 tend to show a steady decline towards the outer disc, rather than a sharp fall that might signify a definitive outer edge of the disc.

The surface density of mm particles, as calculated from the 1.3 mm emission, is fairly large to begin with but begins to drop sharply after $\sim$100 au. The C18O surface density declines over the extent of the disc on the log scale in Figure 10, converging with the g/d scaled continuum curve. The mass calculations in Section 4 result in a disc g/d lower than 100. If values of <100 are used to scale the 1.3 mm surface density curve, in the outer disc C18O would soon have a greater surface density than the scaled dust, suggesting we have enough sensitivity in our gas line observations to confirm that the gas is more extended than the dust. Determining the outer radius of a disc in either gas or dust is complex, but Facchini et al. (2017) use thermo-chemical models to find that the observed difference between dust and gas disc radii is largely down to optical depth effects. The effect of radial drift is shown by a sharp outer edge to the dust disc. Such a feature is absent from the dust surface density profile in Figure 10 due to the $\sim$100 au resolution in these observations.

We see evidence of optical depth effects in this data set in the CO isotopologues; throughout the disc 13CO has a lower surface density than the much less abundant C18O. In order to estimate the optical depth of these species we use Equation 7,

[TABLE]

This is a method adopted from Schwarz et al. (2016) where we have used the ratio of fluxes instead of brightness temperatures. We assume 13CO to be optically thick and the isotopic ratio 13CO/C18O=8 based on ISM values (Wilson, 1999). The ratio of excitation temperatures of the two isotopologues is assumed to be $\approx$ 1. Figure 11 shows the resulting $\tau$(R) as well as the flux ratio from the two isotopologues.

At the outer edge of C18O detection, the ratio of isotopologue emission rises up to $\sim$7, with a steep increase towards the ISM ratio of 8 from $\sim$250 au as the edge of the C18O disc is approached. C18O is found to be optically thin for the full radial distance over which we detect it, with $\tau\leq$ 0.68 as shown in Figure 11. This optical depth result can help to explain the isotopic ratio curve, which rises as the optical depth decreases. Assuming the intrinsic isotopic abundance ratio within the disc remains constant, the 13CO optical depth will decrease as the C18O does with separation from the star. As 13CO becomes increasingly optically thin, the observed abundance ratio will tend towards its true value.

5.4.1 Spectral Index

We now utilise the 867 $\mu$m continuum data from the ALMA band 7 data on this source (Walsh et al., 2014). The continuum from this data is imaged using the same beam as the band 6 1.3mm continuum. Using azimuthal averaged profiles of each continuum image we can directly calculate an average, radially dependent spectral index for the disc, plotted as the central green line in Figure 12. In the upper panel of Figure 12, $\alpha_{\rm mm}$(R) begins in the inner disc with a value of 2.1 and steadily rises with radius by $\Delta\alpha_{\rm mm}\approx$ 0.3 .

$\alpha_{\rm mm}$(R) remains within the range typically seen in protoplanetary discs (Ricci et al., 2012; Pinilla et al., 2014), covering a range of values that are low in comparison to that of the ISM. Low spectral index at millimetre wavelengths can be down to one of two reasons; grain growth has lead to larger dust particles ($a_{\rm max}\geq 3\lambda$) which emit less efficiently and reduce the calculated $\alpha$, or optically thick emission. As we do not have robust information on the sizes of grains present in the disc, we calculate the optical depth of 1.3 mm continuum emission as a function of radius using a range of emissivities at 1.3 mm corresponding to different particle sizes from Draine (2006) . The results are shown in Figure 13 and suggest that over all emissivities tested, the 1.3 mm emission is always optically thin for the entirety of the disc. Through comparing spectral indices calculated from previous mm observations, Walsh et al. (2014) find that at sub-mm wavelengths the continuum emission may be approaching the optically thick regime however.

In order to constrain what the size of the emitting grains are we calculate the opacity index $\beta_{\rm mm}$. In an optically thin disc under the Rayleigh-Jeans approximation, this can be related to the spectral index by $\alpha_{\rm mm}=\beta_{\rm mm}+2$. At $\lambda\approx$ 1 mm we trace the cooler dust, where the Rayleigh-Jeans approximation does not hold, and so we take an approach similar to that of Guidi et al. (2016) and calculate $\beta_{\rm mm}$ with

[TABLE]

We adopt a temperature profile that obeys a power law, $T(R)=T_{100}~{}(R/R_{100})^{0.5}$, and take $T_{100}$=25K. Figure 14 demonstrates how temperature affects the calculated $\beta_{\rm mm}$(R). Although temperatures in the mid-plane and the outer disc are cool and relatively stable, Figure 14 shows that an inaccurate temperature prescription can alter the $\beta$ significantly; a difference of 20K in $T_{100}$ changes $\beta$ at $\sim$200 au by 0.32. For comparison, this difference is larger than the difference between the characteristic $\beta$ for $a_{\rm max}$=1mm and $a_{\rm max}$=1cm. Under-estimating disc temperature increases the inferred size of grains in the cooler outer disc, meaning the radial size sorting of grain sizes is of a lesser degree and that grain sizes thoughout the disc are of large ($\sim$mm) sizes.

In the lower panel of Figure 12 we show $\beta_{\rm mm}$ as a function of radius compared with opacity indexes calculated by Draine (2006). Coloured areas on the plot show radii at which $\beta_{\rm mm}$ indicates a particular $a_{\rm max}$ has been achieved through grain growth, determined by when the lower bound of the calculation crosses the indicative horizontal lines. The curve takes a similar form to that in the top panel, except for a larger variation of $\Delta\beta_{\rm mm}(R)\approx$ 0.8 over the extent of the disc.

The main result from Figure 12 is that $\beta_{\rm mm}(R)$ calculated between 867$\mu$m and 1.3 mm indicates grain growth to cm-sizes throughout the disc. In fact, the characteristic $\beta$ for $a_{\rm max}=1,10,10^{2}$ cm, at millimetre wavelengths are all $\approx 0.8$ (see Draine, 2006, Figure 2), so dust could be growing even larger than cm sizes without us being able to detect it. Within the uncertainties due to flux calibration, this result is consistent with grain growth to cm-size for the full extent of the disc detected at both continuum wavelengths, where we calculate $\alpha$ over the shorter $\sim 230$ au band 7 continuum radius (Walsh et al., 2014). $\beta_{\rm mm}(R)$ for the disc inward of $\sim$200 au is consistent with pebble (>cm) sized dust. The upper curve remains below the top dashed line for the full extent of the disc, confirming that grains have evolved from an ISM-like $\beta$ throughout the disc. Overall, mid-plane dust has grown substantially from ISM levels within the disc.

Radially increasing $\beta_{\rm mm}$ in the outward direction has been predicted by dust evolution models that include realistic coagulation and fragmentation (Brauer et al., 2008; Birnstiel et al., 2010). Radially increasing $\beta$ also agrees with past studies that fit multi-wavelength observations with parametric models of discs in the uv-plane, where only those with a radially increasing $\beta$ could reproduce observed emission (Banzatti et al., 2011; Pérez et al., 2012; Trotta et al., 2013; Tazzari et al., 2016).

Guidi et al. (2016) use high resolution data from ALMA and the VLA to observationally confirm a radially increasing dust opacity in the disc of intermediate mass star HD 163296, a bright disc around a Herbig star that has been found to have a series of dust-depleted rings (Isella et al., 2016). We have been able to confirm this now for HD 100546, adding further evidence for the prediction of Wright et al. (2015) who suggest grains may have grown to at least 5 cm. Walsh et al. (2014) compare their measured flux with previous mm observations by Henning et al. (1994); Wilner et al. (2003) in order to calculate disc-averaged $\beta_{\rm 3mm-1mm}$=0.7-0.8, and $\beta_{\rm 1mm-870\mu m}$=-0.4, values that lie within the range we see over the extent of the disc in Figure 12.

We have assumed in our analysis that dust grains in the disc have a size distribution of $dn/da\propto a^{-p}$, where p=3.5, an index typical of the ISM (Mathis et al., 1977) that is canonically adopted when modelling discs. As grains grow, the overall distribution of mass between grain sizes is likely to evolve. Birnstiel et al. (2011) showed that when both coagulation and fragmentation are included within a disc grain growth model, the size distribution will deviate from p=3.5; growth of grains tends to lower p, whereas fragmentation will re-introduce small particles and raise p. Draine (2006) demonstrates how $\beta$ decreases with index p for a constant $a_{\rm max}$ in Figure 6 of that paper. We would expect ISM-like grains at the outer edge of the disc, yet this is not shown in Figure 12. Either we lack the sensitivity to detect the full extent of the dust, or the dust in the disc would be better suited to an alternative value of p. The latter would corroborate with the results of Wright et al. (2015) when studying HD 100546, where the authors find that only d ’Alessio et al. (2005) disc models with p=2.5 fit their data. Tighter constraints could be achieved by comparing against an extensive grid of opacity calculations at millimetre wavelengths over a range of p values.

6 Conclusions

We have presented new ALMA observations of the HD 100546 protoplanetary disc, including the first detection of C18O in the disc.

We measure an integrated flux of 492 mJy from the 1.3 mm continuum image. We derive an outer-most detection of the dust disc at a radial separation from the star of 290$\pm$10 au at 1.3 mm and the lower limit on dust disc mass was calculated to be 4.1-9.6 $\times 10^{-4}~{}$M*⊙*.

We detect line emission from CO isotopologue transitions 13CO (2-1) and C18O (2-1). 13CO has an integrated flux of 12.872 Jy kms*-1* and extends out to 390$\pm$10 au, a size comparable with previous observations of 12CO in the disc. The optically thinner, mid-plane tracing isotopologue C18O is detected for the first time in the disc, giving an integrated flux of 2.948 Jy kms*-1* with emission detected as far as 300$\pm$10 au.

From these detections we calculate a considerable gas mass of 0.018 M*⊙, around 1% of the stellar mass. We do not expect freeze-out in this warm Herbig disc. If carbon sequestration is efficient the total gas mass may be larger, but not lower, by our calculations. The ALMA observations of continuum and CO isotopologue emission presented here suggest the lower limit of the gas to dust ratio for the HD 100546 disc is 19. We have assumed ISM abundances of C18O relative to H2 in the disc and so the gas to dust ratio in the disc will be larger in the event that C18*O is significantly depleted relative to other species.

We observe low-level extended C18O emission that is spatially coincident with spiral arm structures identified in scattered light observations by the HST. The extended C18O emission may be tracing the mid-plane density variation that creates the scattered light features. In each of our observed tracers, the disc is seen to have a puzzling offset from the known stellar position of between $\approx$70-110 milliarcseconds.

Spectral index, $\alpha_{\rm mm}$, is derived as a function of radius in the disc and rises with separation from the star suggesting a decrease of maximum dust grain size in the outward direction. We directly relate this to the size of grains through calculating the opacity index, $\beta_{\rm mm}$. This gives results consistent with the growth of grains beyond mm-sizes throughout the disc, and to larger cm-sized ‘pebbles’ within $\sim$200 au. A shallower dust grain size distribution is the alternative explanation, suggesting a larger proportion of the solid mass is contained within the largest grains in comparison to the ISM, or indeed optically thick dust at 867$\mu$m.

Acknowledgements

JM is funded through the University Research Scholarship from the University of Leeds. The work of OP is supported by the Royal Society Dorothy Hodgkin Fellowship. TJH is funded by an Imperial College Junior Research Fellowship. We thank Catherine Walsh for useful advice and discussion. We thank David Ardila for providing the HST images of HD 100546. This paper makes use of the following ALMA data: ADS/JAO.ALMA$\#$2015.1.01600.S, ADS/JAO.ALMA$\#$2011.0.00863.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. This research made use of APLpy, an open-source plotting package for Python hosted at http://aplpy.github.com, of Astropy, a community-developed core Python package for Astronomy (The Astropy Collaboration et al., 2013), of NASA’s Astrophysics Data System and the SIMBAD database, operated at CDS, Strasbourg, France.