A multi-wavelength study of the debris disc around 49 Cet

TL;DR

This study combines multi-wavelength observations and modeling to analyze the dust properties and structure of the debris disc around 49 Cet, revealing insights into grain size distribution and disc morphology.

Contribution

It presents new scattered light images and semi-dynamical dust models that unify thermal emission and scattered light data for the 49 Cet debris disc.

Findings

Disc with maximum surface density at 110 au fits observations

Grains near blow-out limit and large grains originate from same planetesimal population

Radiation pressure significantly influences dust grain distribution

Abstract

In a multi-wavelength study of thermal emission and scattered light images we analyse the dust properties and structure of the debris disc around the A1-type main sequence star 49~Cet. As a basis for this study, we present new scattered light images of the debris disc known to possess both a high amount of dust and gas. The outer region of the disc is revealed in former coronagraphic H-band and our new Y-band images from the Very Large Telescope SPHERE instrument. We use the knowledge of the disc's radial extent inferred from ALMA observations and the grain size distribution found by SED fitting to generate semi-dynamical dust models of the disc. We compare the models to scattered light and thermal emission data and find that a disc with a maximum of the surface density at 110~au and shallow edges can describe both thermal emission and scattered light observations. This suggests that grains close to the blow-out limit and large grains stem from the same planetesimal population and are mainly influenced by radiation pressure. The influence of inwards transport processes could not be analysed in this study.

Tables

Table 1. Table 1: Disc parameters inferred in different wavelength ranges.

Wavelength [$\mu \text{m}$]	Radius [au]	Inclination [∘]	PA [∘]	Comments / Instrument
NIR	65 - 250	73	110	VLT/SPHERE, Choquet et al. (2017)
NIR	80 - 280	…	$106\pm 1$	VLT/SPHERE, this work
12.5, 17.9	30 - 60	$60\pm 15$	$125\pm 10$	Keck/MIRLIN, inner disc component, Wahhaj et al. (2007)
70	200	$>44$	$105\pm 1$	Herschel/PACS, Roberge et al. (2013)
70	192	$67.4\pm 2.7$	$109.0\pm 3.9$	Herschel/PACS, Moór et al. (2015)
100	196	$67.2\pm 2.5$	$109.4\pm 4.8$	Herschel/PACS, Moór et al. (2015)
160	209	$56.7\pm 15.5$	$93.4\pm 13.9$	Herschel/PACS, Moór et al. (2015)
450	$421\pm 16$	$74\pm 13$	$130\pm 10$	SCUBA2, Holland et al. (2017)
850	117	$80.6\pm 0.4$	$109.1\pm 0.4$	ALMA, Hughes et al. (2017)

Table 2. Table 2: Parameters for the radial distribution models.

Parameter		Single PL				Double PL
		no ring		with ring
		H17	Fit	H17	Fit	H17	Fit
$i$ [∘]		79.3		79.3		79.2
$PA$ [∘]		108.7		108.8		108.3
$r_{\text{min}}$ [au]		70		58		26
$r_{\text{max}}$ [au]		298		291		302
$p$	range	…	0.0-2.0	…	0.0-1.5	…
	step size	…	0.2	…	0.5	…
	value	1.27	0.6	0.75	0.5	…
$r_{\text{ring}}$ [au]	range	…		…	90-150	…
	step size	…		…	10	…
	value	…		109	140	…
$r_{\text{trans}}$ [au]	range	…		…		…	50-170
	step size	…		…		…	20
	value	…		…		92	110
$p_1$	range	…		…		…	-3.0-0.0
	step size	…		…		…	0.5
	value	…		…		-2.7	-1.0
$p_2$	range	…		…		…	0.0-2.0
	step size	…		…		…	0.5
	value	…		…		1.50	1.0
${\chi }_{\text{red}}^2$	thermal	1.40	1.18	1.30	1.16	1.32	1.15
	scattered	2.38	2.11	2.37	2.16	2.37	1.99

Table 3. Table 3: Continuum flux density.

Wavelength	Flux density			Instrument	Reference
$\left[\mu \text{m}\right]$	[mJy]
0.44	22260	$\pm$	350	TYCHO	1
0.55	21600	$\pm$	200	TYCHO	1
1.24	10180	$\pm$	190	2MASS	2
1.65	6300	$\pm$	130	2MASS	2
2.16	4373	$\pm$	81	2MASS	2
3.35	1990	$\pm$	120	WISE	3
4.60	1294	$\pm$	49	WISE	3
5.86	690	$\pm$	70	IRS	4
7.07	480	$\pm$	50	IRS	4
8.97	320	$\pm$	30	IRS	4
9.0	366	$\pm$	13.5	AKARI/IRC	5
10.8	250	$\pm$	50	Keck/MIRLIN	6
11.40	210	$\pm$	20	IRS	4
11.56	211	$\pm$	21	WISE	3
12.0	330	$\pm$	66	IRAS	6
12.50	200	$\pm$	26	Keck/MIRLIN	8
13.90	180	$\pm$	20	IRS	4
17.90	186	$\pm$	25	Keck/MIRLIN	8
18.0	199	$\pm$	16	AKARI/IRC	5
22.09	238	$\pm$	24	WISE	3
24.00	259	$\pm$	10	Spitzer/MIPS	9
25.0	380	$\pm$	76	IRAS	7
60.0	2000	$\pm$	400	IRAS	7
63.19	2090	$\pm$	350	Herschel/PACS Spec	9
70.00	2163	$\pm$	151	Herschel/PACS	10
71.42	1749	$\pm$	123	Spitzer/MIPS	10
72.84	1950	$\pm$	320	Herschel/PACS Spec	9
78.74	1900	$\pm$	310	Herschel/PACS Spec	9
90.0	1776	$\pm$	295	AKARI/FIS	5
90.16	1880	$\pm$	320	Herschel/PACS Spec	9
100.0	1910	$\pm$	380	IRAS	7
100.0	1919	$\pm$	134	Herschel/PACS	10
145.54	1160	$\pm$	180	Herschel/PACS Spec	9
150	750	$\pm$	500	ISO	4,11
157.68	980	$\pm$	130	Herschel/PACS Spec	9
160.00	1066	$\pm$	75	Herschel/PACS	10
170	1100	$\pm$	500	ISO	4,11
250.00	372	$\pm$	27	Herschel/SPIRE	9
350.00	180	$\pm$	14	Herschel/SPIRE	9
450.0	125	$\pm$	18	JCMT/SCUBA-2	12
500.00	86	$\pm$	9	Herschel/SPIRE	9
850.0	17	$\pm$	3	ALMA	4
850.00	13.5	$\pm$	1.5	JCMT/SCUBA-2	12
1200.00	12.7	$\pm$	2.8	IRAM	13
1300.0	2.1	$\pm$	0.7	CARMA	4
1330	5.5	$\pm$	0.7	ALMA/ACA	14
9000.0	0.0251	$\pm$	0.0055	VLA	15

Table 4. Table 4: Parameters for the grain size distribution assuming a single power law radial profile.

Parameter	Best fit values	Comment
$s_{\text{min}}$ [$\mu \text{m}$]	$5.14\pm 0.76$	…
$s_{\text{max}}$ [$\mu \text{m}$]	10000	thermal emission, fixed
	15	scattered light, fixed
$q$	$3.77\pm 0.05$	…
$f_{\text{d}}$	$8.8\times {10}^{-4}$	…
$T_{\text{dust}}$ [K]	$65\pm 3$	…
$T_{\text{warm}}$ [K]	$125\pm 2$	…
${\chi }_{\text{red}}^2$	1.03	…

Table 5. Table 5: Comparison of thermal emission models

No.	Model	$J$	Free parameters	${\chi }^2$ [$\times {10}^6$]	$BIC$ [$\times {10}^6$]
1	SPL	2	scaling, $p$	$1.24$	1.24
2	SPL + ring	3	scaling, $p$, $r_{\text{ring}}$	$1.22$	1.22
3	DPL	4	scaling, $p_1$, $p_2$, $r_{\text{trans}}$	$1.21$	1.21

Table 6. Table 6: Comparison of scattered light models

Comparing models	${\chi }^2$ [$\times {10}^9$]	$B$
1, 2	8.58	4.8
1, 3	8.58	6.2
2, 3	8.58	4.8

Equations

$$N(r,s)\sim {\left(\frac{r}{r_0}\right)}^{-p}\times {\left(\frac{s}{s_0}\right)}^{-q}.$$

$$\Sigma (r,s)\text{d}s=\pi s^2\times N(r,s)\text{d}s.$$

$$M_{\text{dust}}=\varrho \times V=\varrho \times \frac{4\pi }{3}\underset{s_{\text{min}}}{\overset{s_{\text{max}}}{\int }}2\pi \underset{r_{\text{min}}}{\overset{r_{\text{max}}}{\int }}N(r,s)rdr{s}^{3}ds.$$

$$\frac{R_{\text{disc}}}{R_{\text{blackbody}}}=A\times {\left(\frac{L}{L_{\text{sun}}}\right)}^B.$$

$$\alpha =90^{\circ }-i,$$

$$\alpha =1.22\times \frac{\lambda }{2\times s}=70^{\circ }\times \frac{\lambda }{2\times s}.$$

$$s\approx \frac{\lambda \times 70^{\circ }}{2\times (90^{\circ }-i)}=\frac{1.04\mu m\times 70^{\circ }}{2\times 10^{\circ }}=3.5\mu m.$$

$${\chi }^2=\sum _{i=1}^{N_{\text{pixel}}}{\left(\frac{F_{\text{i}}}{F_{\text{noise}}}\right)}^2$$

$$BIC={\chi }^2+J\times {\log }_e(N),$$

$$\beta \equiv \frac{\left|\overrightarrow{F_{\text{rad}}}\right|}{\left|\overrightarrow{F_{\text{G}}}\right|}=\frac{3L_{\text{star}}}{16\pi Gc{M}_{\text{star}}}\frac{{\overline{Q}}_{\text{pr}}}{\varrho s}.$$

$$a_{\text{d}}=\frac{a_{\text{p}}(1-\beta )(1-e_{\text{p}}^2)}{1-e_{\text{p}}^2-2\beta (1+e_{\text{p}}\cos (f_{\text{p}}))}$$

$$e_{\text{d}}^2=\frac{{\beta }^2+e_{\text{p}}^2+2\beta e_{\text{p}}\cos (f_{\text{p}})}{(1-\beta {)}^2}$$

$$\tan (\varpi -{\varpi }_{\text{p}})=\frac{\beta \sin (f_{\text{p}})}{\beta \cos (f_{\text{p}})+e_{\text{p}}},$$

$$N(r)=N_0\times {\left(\frac{r}{r_0}\right)}^{-p},$$

$${\overline{Q}}_{\text{pr}}=\frac{\int Q_{\text{pr}}(s,\lambda )\times F_{\lambda }\text{d}\lambda }{\int F_{\lambda }\text{d}\lambda }.$$

$$Q_{\text{pr}}(s,\lambda )=Q_{\text{ext}}(s,\lambda )-Q_{\text{sca}}(s,\lambda )\times ⟨\cos (\vartheta )⟩(s,\lambda ),$$

Full text

††pubyear: 2018††pagerange: A multi-wavelength study of the debris disc around 49 Cet††thanks: Based on observations collected at the European Southern Observatory under ESO programmes 198.C-0209(N) and 097.C-0747(A)–B

1 Introduction

Circumstellar discs around young stars are natural by-products of star formation. They serve as reservoir for mass accretion when protostars form and afterwards transform to places where planets can form. At the beginning of their evolution, primordial discs are mostly composed of gas and only a minor mass fraction is present in small solid dust particles. The gas plays a major role in controlling the disc dust dynamics (Beckwith et al., 2000). Due to viscous accretion (Lynden-Bell & Pringle, 1974) and photo-evaporation (Alexander et al., 2006) the gas is removed during the discs evolution process and is mostly depleted during the first $\sim$10 Myr (e.g., Zuckerman et al., 1995; Fedele et al., 2010).

After this evolutionary phase, the dynamics of dust particles, being no longer governed or stabilised by the gas, are strongly influenced by stellar radiation pressure and Poynting-Robertson drag. In such environment the lifetime of grains is much shorter than the lifetime of the host star (e.g., Dominik & Decin, 2003; Wyatt, 2005). In the last decades observations at infrared wavelengths revealed hundreds of dust dominated discs around stars with a wide range of ages (Hughes et al., 2018). Considering the limited lifetime, the dust material of these debris discs could not be leftover from the primordial stage but rather be comprised of second generation particles continuously replenishing from collisions and evaporation of previously formed larger planetesimals (e.g., Wyatt, 2008; Krivov, 2010). Besides dust, recent observations revealed gas, mostly CO molecules, in emission in $\sim$20 of these debris systems (Hughes et al., 2018). Similarly to dust particles the detected gas is also thought to be derived from larger planetesimals and thus having secondary origin (Kral et al., 2017), though in some very gas-rich, young systems a primordial, residual origin also cannot be ruled out (Kóspál et al., 2013).

49 Cet (HD 9672) is one of the most prominent examples for young, dust-rich gas-bearing debris discs. The A1V-type host star has a stellar luminosity of 16$L_{\odot}$ and an effective temperature of 9000 K (Roberge et al., 2013). Its distance is given by the new Gaia data release as $57.0\pm 0.3$ pc (Gaia Collaboration et al., 2016, 2018; Bailer-Jones et al., 2018). 49 Cet is found to be a member of the 40 Myr old Argus association (Torres et al., 2008; Zuckerman et al., 2012; Zuckerman, 2018). The dust disc around 49 Cet was discovered by the Infrared Astronomical Satellite (IRAS) (Neugebauer et al., 1984; Sadakane & Nishida, 1986) and found to be one of the brightest discs with a fractional luminosity above $10^{-3}$ (Jura et al., 1993). It was first spatially resolved in the mid-infrared with the MIRLIN instrument of the Keck II 10 m telescope (e.g., Wahhaj et al., 2007). However, in the mid-infrared regime the disc exhibits only a moderate amount of emission, suggesting that the dust needs to have a low temperature. 49 Cet is one of those rare gaseous debris discs, where besides CO gas (Zuckerman et al., 1995) several other gas compounds have been detected either in emission (Donaldson et al., 2013; Higuchi et al., 2017) or in absorption (Roberge et al., 2014). The origin of this gas is still debated but most current evidences lean towards a second-generation origin (e.g., Hughes et al., 2017).

Recently, the results of the numerous studies of the 49 Cet debris disc were summarised in Choquet et al. (2017), who presented the first analysis of scattered light images of this disc observed by Hubble Space Telescope/NICMOS and Very Large Telescope/SPHERE. The authors generated a schematic view of the system and furthermore, they investigated the possibility of existing but yet unseen planets around the 49 Cet host star. In addition to that, Hughes et al. (2017) presented ALMA images with a spatially resolved surface density distribution of the molecular gas and the dust continuum emission.

Thanks to the wide wavelength coverage of observations of 49 Cet, it is possible to compare its system properties inferred by studies at different wavelengths. An interesting question is whether we see similar disc structures for different dust grain populations. By combining near-infrared scattered light emission, tracing smaller grains, and thermal emission in the sub-mm regime, sensitive to large particles, we are able to address this issue. Furthermore, by combining different kind of observations at multiple wavelengths, it is possible to constrain system parameters which remained unconstrained in former studies focusing on one wavelength range or type of emission. Due to the developments in observational techniques and instruments this approach became feasible for debris discs only recently, but at the time of this study a handful of analyses addressed this topic already (e.g., Augereau et al., 2001; Ertel et al., 2011; Donaldson et al., 2013; Schüppler et al., 2015; Lebreton et al., 2016; Olofsson et al., 2016, Olofsson et al., in prep.). Nevertheless, most resolved discs (146, List of resolved debris discs111https://www.astro.uni-jena.de/index.php/theory/catalog-of-resolved-debris-disks.html Hughes et al., 2018) do not have data in both wavelength regimes so that a detailed analysis of the disc extent in scattered light and thermal emission remains a future project.

Large dust particles are most effectively traced by observations of thermal emission at long wavelengths (e.g., far-infrared or sub-millimetre). Since these grains are not sensitive to radiation pressure forces we assume their position to be close to the planetesimal belt invisible for us. On the other hand, scattered light in the near-infrared is dominated by small particles which are highly affected by radiation pressure and can have highly eccentric orbits as a consequence. As a result we would expect that by observing small grains we should see a more extended dust disc than by tracing large particles at mm-wavelengths. Indeed, for example the debris disc around $\beta$ Pic shows an extent of $\sim$150 au at sub-mm wavelengths while in near-infrared a halo of small particles becomes visible up to 1800 au (Ballering et al., 2016). However, for some discs such as Fomalhaut comparable disc radii are found in thermal emission and scattered light (Holland et al., 1998; Dent et al., 2000; Kalas et al., 2005). In the case of 49 Cet, we are fortunately able to compare former ALMA studies (Hughes et al., 2017) with observations in scattered light (this work, see also Choquet et al., 2017) to address the question of the disc extent as well. An overview of the theoretical background for the dust dynamics in debris discs is presented in Appendix A.

In this paper, we will present new Y-band data taken with VLT/SPHERE and compare them to previous SPHERE H-band data from Choquet et al. (2017). The observations and data reduction of the SPHERE data are described in Section 2. We will analyse the radial extent of the disc at different wavelengths in Section 1. The modelling of the disc is described in Section 4. Furthermore, we will investigate the question of whether it is possible to find a model which can fit the disc parameters inferred from both thermal emission and scattered light images. Section 5 shows the analysis of the different models investigated, while the results are discussed in section 6. A summary is given in section 7.

2 Observations and data reduction

2.1 Broad band SPHERE data

We observed the disc around 49 Cet in the programme 198.C-0209(N) on 19${}^{\text{th}}$ November 2016 for 1.8 hours on source with the SPHERE/IRDIS instrument of the VLT (Beuzit et al., 2019; Dohlen et al., 2008). The observations were carried out using the broad band Y filter with a central wavelength of 1.04$\,\mu\hbox{m}$, a width of 139 nm and the coronagraph N_ALC_YJ_S (Martinez et al., 2009; Carbillet et al., 2011) with a diameter of 185 mas. The observations were performed in pupil tracking mode to allow for angular differential imaging (ADI, Marois et al., 2006). The total on-sky rotation during our observations was 70*∘*.

After basic reduction steps (flat fielding, bad-pixel correction, background subtraction, frame registration, frame sorting) we processed the data with a classical ADI reduction technique, which consisted of building a model PSF from the mean of all pupil-stabilised images, which was then subtracted from each frame before de-rotating and stacking the images.

In Fig. 1 the result of the ADI-reduction of the 49 Cet debris disc is shown. The image was normalised to mJy/arcsec2 in the following way. On the non-coronagraphic image, we measured the flux density encircled within a circle of radius 0.1 arcsec, encompassing the PSF core, wings and diffraction spikes from the spiders. Then this flux density is corrected by the transmission of the neutral density filter used to obtain the non-coronagraphic image, and by the ratio between the Detector Integration Time (DIT) of the coronagraphic and non-coronagraphic images, to obtain a reference conversion value. To convert the coronagraphic image from ADU to mJy/arcsec2, the coronagraphic image is divided by the reference conversion value, multiplied by the stellar flux density found to be 11.9 Jy at the central wavelength of the Y band and divided by the pixel surface area in arcsec2. The pixel scale of IRDIS is 0.01225 arcsec/pixel (Maire et al., 2016). The ADI reduction algorithm induces self-subtraction of any extended astrophysical signal (Milli et al., 2012). The flux displayed in Fig. 1 did not take this self-subtraction effect into account, which requires a forward-modelling approach and is described in Section 4.3.

2.2 Narrow band SPHERE data

We also observed 49 Cet on 23${}^{\text{rd}}$ July 2016, with VLT/SPHERE in the framework of an open time programme 097.C-0747(A). We used pupil-stabilised imaging with the N_ALC_YJH_S coronagraph, having a coronagraphic mask with a diameter of 185 mas. We used the IRDIFS observing mode, which provided IFS observations in the Y-J range and IRDIS data in the H23 dual band. The data reduction was carried out by the SPHERE Data Centre using their pipeline (Delorme et al., 2017). The obtained reduced master cubes were then utilised as input for high contrast imaging post-processing performed by the Speckle Calibration (SpeCal) package (Delorme et al., 2017; Galicher et al., 2018).

The reduced image in ADU is depicted in Fig. 2.

3 Disc geometry and radial profiles

In our ADI-reduced Y-band image (Fig. 1) we clearly detect the debris disc of 49 Cet between a radial distance of $1.4\arcsec$ (80 au) and $\sim 4.9\arcsec$ (280 au). We note that these values are detection limits and thus, might not resemble the true, possibly larger, disc extent. At an angle of $45^{\circ}$ to the semi-major axis we see dark patterns as residuals from the reduction process (shown as No. 1 in Fig. 1.). In addition to that, there seems to be emission at a $90^{\circ}$ angle which is identified as diffraction pattern residuals stemming from the star (shown as No. 2 in Fig. 1.). In the reduced H23-band image the disc is detected between 80 and $\sim$200 au where the smaller extent is caused by a weaker detection of the disc. Nevertheless, in contrast to the Y-band data, neither the dark pattern, nor the diffraction spikes are visible in the H23-band.

In the following study we will focus on the Y-band data due to the stronger detection of the debris disc. Furthermore, we will concentrate on the resolved outer disc component and will not take a possible inner ring into account which was proposed at a location of $\sim 10$ au by former studies (e.g., Wahhaj et al., 2007; Chen et al., 2014; Pawellek et al., 2014).

We derived the PA of the disc in Y-band, and found a best value of $106.2^{\circ}\pm 1.0^{\circ}$. To do so, we de-rotated the disc to align the major axis with the horizontal of the de-rotator. Then we fitted the vertical profiles of the disc between a separation of 45 and 105 pixels with a Gaussian profile. We iterated with the de-rotation angle until the centroids of the Gaussian had a slope of zero (e.g., Lagrange et al., 2012; Milli et al., 2014). Comparing our PA to values from former studies which lie between 93 and 130 degrees (see Tab. 1) we find a good agreement with the Herschel/PACS measurements. Comparing directly to other SPHERE observations is difficult since the value of 110∘ inferred by Choquet et al. (2017) was not fitted but fixed.

3.1 Radial profiles of scattered light images

We extracted the radial profile of the surface brightness using the same method as described in Choquet et al. (2017) where the authors analyse H-band data of 49 Cet. In this method we produce slices along the semi-major axis with a length of 93 pixels above and under the semi-major axis and a width of 2 pixels. Then we calculate the mean value of the flux density for each slice. With the sizes chosen for the slices we are able to directly compare our Y-band profile to the H-band data of Choquet et al. (2017). We estimate the noise level of the images by generating similar slices as for the radial profile itself. Then we rotate these slices by 90∘ to get the perpendicular direction to the disc. Finally, we calculate the standard deviation of each slice.

The result is shown in Fig. 3, where we normalise the surface brightness to the flux density of the host star given as 6.4 Jy for H-band and 11.9 Jy for Y-band.

Our Y-band observations reach a disc signal-to-noise level of 7, which is stronger than for the H-band data due to a longer exposure time. The noise that dominates at short separation is the residuals from the subtraction process of the Point Spread Function (PSF). Beyond $\sim 3$ arcsec, we are in the background limited regime where the background noise dominates (sky background + readout noise of the detector). The noise level is in the same order of magnitude as the disc signal within a region of 80 au and thus, we will exclude the inner region from further analyses. Furthermore, the disc is detected within $\sim$280 au based on the noise level as well. In case of the H-band data the eastern side of the disc is well detected while the western side remains close to the noise level. In addition, the disc is found to be of a similar extent as suggested by our Y-band image. Analysing the radial profile we find local peaks of the surface brightness in an area between $\sim$130 au and $\sim$170 au in the eastern and western direction in both SPHERE images.

3.2 Radial profile of thermal emission

In order to compare the brightness profile of the near-infrared scattered light with the thermal dust emission in the sub-millimetre, we use ALMA data directly taken from Hughes et al. (2017) who presented 850$\,\mu\hbox{m}$ continuum images for 49 Cet with a beam size of $0.47\times 0.39$ arcsec. To obtain the radial profile we apply the same method as described for the scattered light observations, but adapt the size of the slices according to the ALMA pixel scale of 9.76562 mas/pixel. The resulting radial profile is depicted in Fig. 4. In contrast to the calculated noise level as a function of radial distance for SPHERE data we use the noise level of 0.056 mJy/beam as given in Hughes et al. (2017) to infer the disc radius.

We find that the eastern wing extends to $\sim$280 au and the western side up to $\sim$250 au which is in agreement with the scattered light data. However, in case of the thermal emission measurements the radial profile shows a clear drop of the surface brightness around $\sim$200 au down to values comparable to the noise level of 0.056 mJy/beam so that we assume to see the real disc extent rather than a radial detection limit.

Comparing our inferred disc extent with values derived by former studies (see Tab. 1) we find comparable disc radii for different Herschel/PACS observations as well as for the SPHERE H-band observation.

4 Modelling

4.1 Modelling strategy

The 49 Cet debris disc system was subject to several modelling projects, either concentrating on the dust (e.g., Wahhaj et al., 2007; Pawellek et al., 2014; Choquet et al., 2017), the gas (e.g., Roberge et al., 2014; Miles et al., 2016) or both components (e.g., Zuckerman & Song, 2012; Roberge et al., 2013; Nhung et al., 2017; Hughes et al., 2017). In this study we focus on the dust component. Our goal is to investigate whether disc models inferred from thermal emission images can explain the scattered light observations obtained in near-infrared.

We start with a fit of the spectral energy distribution (SED) of 49 Cet to get information on the dust grain size parameters (see Section 4.2 for a detailed description). Then we study the radial distribution of the particles. In a former study Hughes et al. (2017) investigated three different surface density profiles for the grains visible at sub-mm wavelengths: a single power law, a power law including a narrow dust ring, and a double-power law model. We assume that large dust particles trace the planetesimals best and thus, use the parameters of the three surface density profiles mentioned before to generate the plantesimal belts for our so-called semi-dynamical disc models. These models assume that dust grains are released from the planetesimals and that their orbits are altered by stellar radiation pressure. We use the term “semi-dynamical” since no time-dependency of moving particles is taken into account. The actual profiles of the planetesimal belts from which the dust is produced are depicted in Fig. 5 while their parameters are specified in Tab. 2. The noise visible in Fig. 5 stems from the number of planetesimals used for the model setting. The sizes of the dust grains follow the distribution inferred by the SED modelling. From the semi-dynamical models we generate thermal emission and scattered light images. Details for this procedure are described in Section 4.3. In a final step the resulting images are compared to the observational images by forward modelling. At the time of the ALMA study (Hughes et al., 2017), the new Gaia distances were not available yet, so we correct the inner and outer disc radii for them.

4.2 SED modelling

The photometric data are collected from the literature and summarised in Table 3. We use the SONATA code (Pawellek et al., 2014; Pawellek & Krivov, 2015) to fit the SED and apply the same stellar photospheric model used in Pawellek et al. (2014) to determine the influence of the host star. Here, the stellar temperature, metallicity, and surface gravity are taken into account to generate a synthetic spectrum by interpolating the PHOENIX-GAIA model grid (Brott & Hauschildt, 2005). The dust composition is assumed to be astronomical silicate (Draine, 2003) with a bulk density of 3.3 g/cm3.

The SONATA code calculates the temperature and the thermal emission of each dust particle at a given distance to the star, where the number of particles is determined by the dust mass. Then it sums up the emission of all particles to generate the SED. The flux densities given for wavelengths shorter than 10$\,\mu\hbox{m}$ are not used to fit the dust disc since in this wavelength regime the stellar photosphere rather than the dust dominates the emission. The code assumes a power-law for both the radial and the size distribution of the dust using the surface number density $N(r,s)$:

[TABLE]

Here, $r$ represents the disc radius, $s$ the grain radius, $s_{0}$ and $r_{0}$ normalisation factors and $p$ and $q$ the power-law indices for the radial and size distribution. The surface number density is directly connected to the surface density, $\Sigma(r,s)$, by

[TABLE]

Due to the fact that SEDs use photometric data integrated over the whole disc area, they cannot differentiate between different radial profiles. Therefore, we assume the same single power law model given by Hughes et al. (2017) with $p=1.27$, $r_{\text{min}}=70$ au and $r_{\text{max}}=298$ au to fit the SED.

We assume the grain sizes to lie between a minimum and a maximum value, $s_{\text{min}}$ and $s_{\text{max}}$. We fix the maximum grain size to 1 cm, because larger grains do not contribute any more to the SED in the wavelength range observed. Furthermore, we fix the radial parameters and hence, we are left with three free parameters to fit, namely the minimum grain size, $s_{\text{min}}$, the size distribution index, $q$ and the amount of dust, $M_{\text{dust}}$, for particles between $s_{\text{min}}$ and $s_{\text{max}}$. We infer the amount of dust by using the bulk density $\varrho$ and the dust volume $V$:

[TABLE]

Roughly 2/3 of all discs investigated so far show evidence for an asteroid belt analogue (also called warm component) in addition to an Edgeworth-Kuiper-belt analogue (cold component, e.g., Ballering et al., 2013; Chen et al., 2014; Pawellek et al., 2014). We apply the criteria given in Pawellek et al. (2014) to check for the presence of a warm component for 49 Cet and indeed, our modelling results suggest the existence of an inner dust belt as well. Since there is a large degeneracy between the outer and inner disc parameters, we have to make several assumptions for our SED model. At first we have to identify the location of the warm component which is not spatially resolved in the images. To do so, we fit the warm dust component with a pure blackbody model to infer its blackbody radius. Here we assume the same radial distribution index ($p=1.27$) as for the outer component. The blackbody radius is found to be $\sim 8$ au. Now, we estimate the “true” disc radius by applying the method presented in Pawellek & Krivov (2015) and Pawellek (2017). Here the authors found a relation between the true disc radius and the blackbody radius in the shape of

[TABLE]

For pure astronomical silicate Pawellek (2017) give the parameters as $A=6.49\pm 0.86$ and $B=-0.37\pm 0.05$. Using these values and the solar luminosity of 16$L_{\odot}$ we find an estimated “true” location of the warm component of $\sim 16$ au. We assume the same minimum grain size and size distribution index for both the warm and cold component and therefore increase the number of free parameters from three to four by adding the amount of dust of the warm component. The SONATA code uses the simulated annealing approach (Pawellek et al., 2014) to fit the SED. Considering the scattered light images, we are not sensitive to the region within 80 au. Thus, we focus on the cold dust component of the debris disc.

In the best fitting SED model the minimum grain size is $s_{\text{min}}=5.14\pm 0.76\,\mu\hbox{m}$ and the size distribution index $q=3.77\pm 0.05$. These parameters are applied to all semi-dynamical models with the different radial profiles analysed in this work. We checked the validity of our assumption by using the double power law as radial profile for the SED (Tab. 2) and found the grain size parameters to be close to the values of the single power law model lying within their given confidence intervals. The resulting SED is shown in Figure 6 while the grain size parameters are given in Table 4. The amount of dust, $M_{\text{dust}}$ , in debris discs depends on the SED model settings, i.e., the grain size and radial distribution and is therefore difficult to compare for different studies. Hence, we focus on the fractional luminosity which is only determined by photometric data points. It is given as $f_{\text{d}}=8.8\times 10^{-4}$.

The error bars give the uncertainties of the fit parameters and are inferred in the following way. We start from the position of the minimum $\chi^{2}$ in the parameter space. New parameter values are generated and the resulting $\chi^{2}$ of the model is compared to its minimum value. There is a probability, that a new minimum can be found in the direction of the new parameter values. If this probability reaches a critical value then the fit parameters are saved. In the end, it is counted how often the code reaches a certain parameter value. The resulting distribution in the parameter space represents an estimate for the probability distribution of the parameters and thus, allows us to calculate the confidence levels for the parameters assuming that the values follow a normal distribution (simulated annealing, e.g., Pawellek, 2017).

4.3 Generating images

After fitting the SED we generate the semi-dynamical disc models using the MODERATO-code (Olofsson et al., in prep., for similar approaches see Wyatt et al., 1999; Lee & Chiang, 2016). The code assumes surface density profiles for the parent belt where each planetesimal releases grains of different sizes following a power law (see Eq. 1). Then the orbits of the individual dust particles altered by stellar radiation pressure are calculated and from their position the scattered light and thermal emission are inferred. We correct the possible overabundance of highly eccentric small particles in the way described by Strubbe & Chiang (2006) and Thébault & Wu (2008).

The radial and grain size distribution assumed for the models are given in Tabs. 2 and 4. For each radial profile a small grid of parameter values is assumed to find the best fitting model. Due to the high calculation time for one run, Monte-Carlo fitting is not applicable. In all cases we fix the inclination, PA, and the minimum and maximum disc radius to the values inferred by Hughes et al. (2017) to generate the thermal emission maps. Considering the scattered light maps we use the PA of 106.2 inferred in this work.

In the first scattered light study of 49 Cet, Choquet et al. (2017) use the anisotropic scattering approach by Henyey & Greenstein (1941) where a single scattering asymmetry factor,$|g|$, is assumed for the whole material without a differentiation between grain sizes. In scattered light studies of debris discs, this approach is often used to get a general idea of the scattering properties of the disc material (e.g., Schneider et al., 2006; Thalmann et al., 2013; Millar-Blanchaer et al., 2015; Olofsson et al., 2016; Engler et al., 2018). Since we generate a physical model with a distribution of particles inferred by SED fitting, we use Mie theory instead to stay consistent and assume pure astronomical silicate (Draine, 2003) as dust composition similar to the SED modelling. However, we additionally use the Henyey-Greenstein approach to compare our results to previous studies.

We are aware that using Mie-theory to calculate the optical properties for scattered light models leads to an overestimation of forward scattering, especially for larger particles ($\geq 20\,\mu\hbox{m}$) since the grains are assumed to be compact spheres (e.g., Schuerman et al., 1981; Bohren & Huffman, 1983; Weiss-Wrana, 1983; Mugnai & Wiscombe, 1986; McGuire & Hapke, 1995).

To avoid this effect we assess the maximum grain size still contributing to the scattered light. This is done by estimating the width of the forward scattering peak. The minimum scattering angle achievable for a disc is given by the disc inclination:

[TABLE]

which is $\approx 10^{\circ}$ in the case of 49 Cet. On the other hand, analogue to a telescope, a dust particle of radius $s$ can “detect” the incoming light (Babinet’s principle) with the geometric area $\pi s^{2}$ and therefore, the Rayleigh criterion provides the minimum angle observable for the grain at a certain wavelength $\lambda$:

[TABLE]

Equalising the angles of Eqs. 5 and 6, it is possible to infer a grain radius for which the scattered light observations are still sensitive. We get an estimate for the maximum particle size of the 49 Cet debris disc as:

[TABLE]

This shows, that the scattered light contribution of particles larger than $3.5\,\mu\hbox{m}$ is minor, which is in agreement with Zubko (2013). In this study, the authors analysed the contribution of large particles in cometary dust by using the Discrete Dipole Approximation method (DDA). They found that for a dust composition with moderate or high absorption, such as carbonaceous material, grains with a size parameter $x=2\pi\times s/\lambda>15$ are not needed to model the scattered light observations. However, adding large grains to the model did not change the results significantly, although the computation time was much higher.

With SED fitting we found a minimum grain size of $\sim 5\,\mu\hbox{m}$ which is not contradictory to the above mentioned estimate. It just states that the larger particles are not as effectively traced as smaller grains. Due to possible estimation uncertainties we set the maximum particle size for scattered light observations to a value of $15\,\mu\hbox{m}$.

The semi-dynamical models are compared to scattered light observations by subtracting the convolved disc model from our pupil-stabilised cube of frames. Then we perform ADI on this cube to get the final image. In case of thermal emission, the models are convolved with the assumed ALMA beam and then subtracted from the observational image. We use a scaling factor to adapt the flux density of the model to our observations for scattered light and thermal emission. This is inevitable since calculating the flux density of all dust particles necessary to equal the measured emission exceeds the computational resources available. Due to the different subtraction processes we note that the scaling factors are different for scattered light and thermal emission. To estimate the quality of the model matching the observations we use $\chi^{2}$-minimisation assuming that an ideal residual image should expose no emission in any pixel. The $\chi^{2}$-parameter is then computed for each pixel by

[TABLE]

In case of the scattered light data, the noise (or error) is estimated in the following way. We use the noise level inferred from the observational image (see Fig. 3). Then we rotate this profile around the centre of the image to generate a noise map. In case of the thermal emission image we use the noise level given as $F_{\text{noise}}=0.056$ mJy/beam (Hughes et al., 2017). The number of pixels for the thermal emission maps is $1024\times 1024$ and for the scattered light images $255\times 255$.

5 Analysis

5.1 Thermal emission

The thermal emission image observed by Hughes et al. (2017) with ALMA at 850$\,\mu\hbox{m}$ is compared to the semi-dynamical models based on the three different radial profiles given in Tab. 2 and illustrated in Fig. 5. The best fit results are shown in Figs. 7-9.

5.1.1 Single power law distribution

At first a single power law model is assumed as radial distribution. Using the parameter values from Hughes et al. (2017), where the disc is confined between 70 and 298 au and possesses a surface density power law exponent of 1.27, we get a reduced $\chi^{2}$ of $\chi^{2}_{\text{red}}=1.4$.

The grain parameters derived by SED-fitting are applied to our model, where the minimum grain size is 5.14$\,\mu\hbox{m}$, the maximum grain size 1.0 cm and the distribution index $q=3.77$.

Next, we analyse a small parameter grid. Here, we fix the minimum and maximum disc radius as well as the inclination and PA, while only one parameter, $p$, is varied between 0.0 and 2.0 in steps of 0.2 . The result of this fit is shown in Fig. 7, where $p=0.6$ is the best fit value. We could achieve $\chi^{2}_{\text{red}}=1.18$. After subtracting the model from the observed data, there are some $2\sigma$ residual structures left throughout the disc. However, they do not seem to be of systematic origin.

5.1.2 Single power law distribution with a narrow ring

In a second approach a narrow dust ring is added to the single power law of the radial distribution. The parameter values can be found in Tab. 2. Here, $N(r,s)$ of the extended disc is calculated at the location of the ring. Then the surface number density of the ring is set to a value twice as high which is in agreement with the parameters given by Hughes et al. (2017). In the final step, the grains of both model parts (disc + ring) are summed up to derive the thermal emission map.

The parameters given by Hughes et al. (2017) lead to $\chi^{2}_{\text{red}}=1.30$ which is an improvement to their pure single power law model. The ring is confined between 108 and 110 au and possesses the same grain size parameter values as the extended disc between 58 and 291 au.

Fitting a grid of models, we vary two parameters: the distribution index $p$ between 0.0 and 1.5 in steps of 0.5 and the position of the included ring ($r_{\text{ring}}$) between 90 and 150 au in steps of 10 au. The results are shown in Fig. 8 where we found a best fitting model with $p=0.5$ and a ring position of $r_{\text{ring}}=140$ au. The parameter $\chi^{2}_{\text{red}}$ is 1.16 and thus, slightly smaller than for the pure single power law model. Analysing the residual image (right panel) there are less residuals found in the inner part of the disc in contrast to the single power law model while in the outer region some larger residuals are left.

5.1.3 Double power law distribution

The third approach contains a double power law distribution as radial profile. To implement this setting, a similar approach is used as for the former model. We calculate $N(r,s)$ at the position of the transition radius, $r_{\text{trans}}$. At this position, one power law changes into the other. We demand that $N(r,s)$ has to be the same value for both parts of the radial profile to ensure a continuous surface density profile.

Hughes et al. (2017) inferred a transition radius between both profile parts at 92 au while the inner disc region is given with $p_{1}=-2.7$ and the outer fixed to $p_{2}=1.5$. Using these values we get a $\chi^{2}_{\text{red}}$ of 1.32.

Now we vary the transition radius between 50 and 170 au in steps of 20 au and the radial distribution indices between -3.0 and 0.0 for the inner disc region and 0.0 and 2.0 for the outer one in steps of 0.5 each, leading to three free parameters. The best fitting model yields $r_{\text{trans}}=110$ au, $p_{1}=-1.0$ and $p_{2}=1.0$ and leads to $\chi^{2}_{\text{red}}=1.15$. The result is depicted in Fig. 9. Comparable to the single power law model including a narrow ring the inner disc region shows no systematic residuals while in the outer part some $2\sigma$ remnants are left.

5.1.4 Comparing the models

By comparing the different radial profiles we find that all models achieve a similar $\chi^{2}_{\text{red}}$ value between 1.15 and 1.18. Thus, a conclusion of which model gives the best result is hardly possible. Therefore, we apply the Bayesian information criterion ($BIC$) being defined as

[TABLE]

where, $J$ gives the number of free parameters and $N$ the number of data points. The resulting values for each model are given in Tab. 5.

Following the classification given in Kass & Raftery (1995) we compare each model by calculating the parameter $B=2\times\log_{e}(\Delta BIC)$, where $\Delta BIC$ represents the difference of $BIC$ between the models. If $B$ lies between 0 and 2 the more complex model is not significantly better than the simple model, by values between 2 and 6 it is possible that the more complex model is better, and between 6 and 10 the more complex model is more likely. For values larger than 10 the probability of the more complex model is much higher than for the simple one.

We now compare the models 2 and 3 (see Tab. 5) with the single power law model. Model 2 reaches a value of $B=19.9$, while model 3 gives $B=20.1$. This indicates that the single power law model including a ring and the double power law model have a much higher probability than the single power law model. In the final step we compare model 2 with model 3 and find $B=18.5$. Thus, although the difference in $\chi^{2}_{\text{red}}$ is only minor, the information criterion clearly states that the double power law model is preferable to the power law model including a narrow ring.

Our fitting results deliver slightly different values than provided by the model of Hughes et al. (2017). The main difference here is that in the aforementioned study the applied surface density models for the dust grains do not consider a grain size distribution or any forces altering the particle orbits. We assume a dust orbital model including radiation pressure and a particle size distribution to generate thermal emission maps. Nevertheless, both approaches show that a single power law distribution might not be the appropriate way to simulate the debris disc of 49 Cet.

5.2 Scattered light

In this section the best fitting models inferred from the ALMA data are compared to the scattered light data presented in this work. For reasons of comparability we also analyse the models with parameters stemming from Hughes et al. (2017). We use the grain size values given in Tab. 4 and the radial parameters listed in Tab. 2. The results are shown in Figs. 10-12.

The surface density distribution of each planetesimal belt is calculated in the same way as for the thermal emission maps. Each planetesimal then releases dust particles which are influenced by radiation pressure. From the position of the grain the appropriate scattered light is derived and summed up into a map. In contrast to the thermal emission images the scattered light maps possess a complex artefact structure inclined by 45∘and $90^{\circ}$ to the major axis which makes the interpretation of the modelling results difficult. So far there is no criterion to distinguish the artefact areas from the actual scattered light so that we cannot exclude these parts from $\chi^{2}$ calculations.

Nevertheless, we find that all three surface density profiles can describe the scattered light observations to a certain extent, but we see residuals along the major-axis and over-subtraction in the inner region close to the central star for all three radial profiles. However, the residuals in the outer disc region are at the same order as the estimated noise level of the image.

Due to the artefact structures in the image in combination with the method to find the best fit (see Eq. 8), the $\chi 2$ values inferred by forward modelling are very large and thus, we use again the $BIC$ introduced in the former section to compare the models with each other. The results are listed in Tab. 6.

By following the classification of (Kass & Raftery, 1995) again, we find that the probability of a double power law profile is higher than the single power law and the power law including a narrow ring. Hence, the scattered light data lead to similar results as the thermal emission data showing that a double power law profile as depicted in Fig. 5 is more likely to occur in the debris disc around 49 Cet than the other profiles investigated in this study.

6 Discussion

6.1 Radial profiles

We found that the radial extent of the 49 Ceti debris disc is the same for both ALMA and SPHERE observations as shown in Fig. 4. As theory predicts, the orbits of small grains should be highly altered by radiation pressure and thus, we would expect a larger disc extent for our scattered light image than for the thermal emission map. A reason could be the rather noisy scattered light detection of the disc only showing the brightest parts of it so that we underestimate its outer radius in scattered light.

In comparison to other disc detections with SPHERE (e.g., Olofsson et al., 2016; Wahhaj et al., 2016; Engler et al., 2018) the disc of 49 Cet while possessing a large fractional luminosity is rather faint in its surface brightness. On the one hand this might indicate a dust composition strongly absorbing, such as carbon. On the other hand it might indicate the absence of small dust particles which is supported by the results of the SED fit and other studies, focusing on thermal emission and reporting that the disc is faint even in the mid-infrared (Wahhaj et al., 2007). Another explanation is that the disc possesses a broad extent from at least 80 au up to 280 au (see Fig. 3). Thus, a high number of small dust grains might be present in the disc, but due to its large on-sky area the particles are widely distributed so that the surface brightness stays low. In addition to that, the ADI reduction process (Marois et al., 2006) leads to a stronger subtraction for broad discs and so the low detection might be an observational bias.

While the single power law profile was already excluded as the best-fit model by the aforementioned ALMA study (Hughes et al., 2017), the power law including a ring and the double power law profile remained indistinguishable for the analysis of thermal emission. To a certain extent we could reproduce this result with our semi-dynamical models, but we found that all three profiles are capable of describing the observational data (Fig. 4). Nevertheless, we found different parameter values for the best fitting fiducial discs. Applying the Bayesian information criterion we found that despite the comparable $\chi^{2}$ values for all three profiles, the double power law has the highest probability.

Summarising the results we infer the following architecture of the outer region of the debris disc around 49 Cet. From observed surface brightness profiles we find that the Kuiper-belt analogue is located between 80 au and 280 au and shows an increase of small particles between 130 au and 170 au visible in scattered light (Fig 3). Assuming the best-fitting double power law model, a maximum for large grains, seen at sub-mm wavelengths, is found around 110 au. The edges inwards and outwards of the ring are shallow.

6.2 Size distribution

In Pawellek et al. (2014) the authors found a weak relation between the minimum grain size and the stellar luminosity of the host star. However, they stated that this trend is also consistent with being independent of the stellar luminosity and found an average value of roughly 5 $\mu$m to be valid for the majority of debris discs investigated. This is in good agreement with our SED-fit results where a value of $5.14\pm 0.76\,\mu\hbox{m}$ was found as dominant grain radius. The ratio of the minimum particle size to the inferred blow-out limit of 2.9$\,\mu\hbox{m}$ for astrosilicate is 1.7. While it is directly connected to the dynamical excitation of dust-producing planetesimals in the disc, the low value of the ratio argues for a collisionally active disc for which a value around 2 is suggested (e.g., Krivov et al., 2006; Thébault & Augereau, 2007).

As shown in former studies (e.g., Pawellek, 2017), the size distribution index, $q$, is mostly influencing the steepness of the SED longwards of the far-infrared wavelengths. Hence, the now available VLA data at 9 mm (MacGregor et al., 2016) help to strengthen the constraint. Compared to an idealised collisional cascade, where the index lies at $q=3.5$ (Dohnanyi, 1969), we inferred a value of $3.77\pm 0.05$. This is somewhat larger, but still in agreement with collisional models (e.g., Löhne et al., 2008; Gáspár et al., 2012; Kral et al., 2013; Löhne et al., 2017). The result supports the assumption of a dynamically excited disc producing more grains closer to the blow-out limit.

Besides a higher excitation level which is expected for discs around earlier-type stars (Pawellek & Krivov, 2015), the disc around 49 Cet is known to host a high amount of gas (e.g., Zuckerman et al., 1995; Hughes et al., 2017, Moor et al., submitted). The presence of gas might lead to a longer residence time of sub blow-out grains produced by collisions so that we would expect an increased number of small grains (Kral et al., 2018). However, we found the surface brightness to be rather weak in scattered light speaking against this assumption.

A reason for a possible lack of small particles might be given by the blow-out limit for pure astronomical silicate, showing that grains smaller than 3$\,\mu\hbox{m}$ are expelled from the 49 Cet system (Burns et al., 1979; Bohren & Huffman, 1983). Fig. 16 shows that grains smaller than 1$\,\mu\hbox{m}$ possess the highest scattering efficiency. For larger grains this parameter stays nearly constant and thus, the scattered light is mainly determined by the particle cross section which is decreasing with increasing particle size, assuming a power law size distribution (e.g., Dohnanyi, 1969). Furthermore, considering the inclined disc of 49 Cet, grains larger than 3.5$\,\mu\hbox{m}$ are not as effectively contributing to the scattered light as smaller particles at the observational wavelength used (Zubko, 2013).

Another possibility for the lack of small grains is a dynamically “cold” disc where an imbalance between the particles’ production and destruction rate leads to a natural depletion of grains up to a few times the blow-out size (see Fig. 7 of Thébault & Wu, 2008). Furthermore small-grain depletion might be caused by the surface energy conservation criterion suggested by Krijt & Kama (2014). Here, the minimum collisional fragment size is determined by the energy necessary to form the small grains. Nevertheless, the aforementioned mechanisms seem to be negligible for discs around A-type host stars (Thebault, 2016).

Thus, the explanation of a low surface brightness due to a broad on-sky area of the disc might be more appropriate than a lack of small grains although further studies are needed to confirm this suspicion.

6.3 Model images

6.3.1 Thermal emission

Considering the thermal emission maps we found elongated residuals in the outer disc region for the single power law profile including a narrow ring and the double power law profile which can be explained in different ways. Firstly, they might be modelling artefacts since for the single power law model without a ring these extended residuals are not visible. Secondly, the observations might suggest that there is a higher amount of grains in several parts of the outer disc than provided by the models. We do not see a symmetric distribution such as ring shaped residuals and thus, it is possible that the higher amount of particles is caused by recent collisions between massive bodies (Olofsson et al., 2016) or stirring effects, such as self-stirring (e.g., Wyatt, 2008; Kennedy & Wyatt, 2010) or planetary stirring (e.g., Mustill & Wyatt, 2009), although there were no planets found in the 49 Cet system so far (Choquet et al., 2017).

Besides the dust stirring mechanisms the particles might be dragged outward by present gas (e.g., Takeuchi & Artymowicz, 2001; Thébault & Augereau, 2005; Krivov et al., 2009; Kral et al., 2018, Moór et al., submitted). Nevertheless, the influence of the gas on dust particles may be limited to the smallest grains ($<5\,\mu\hbox{m}$) not sensitively traced at sub-mm wavelengths.

6.3.2 Scattered light

Comparing the scattered light maps to our ADI-reduced image we found that the double power law model reproduces the observations best. However, residuals along the semi-major axis can be found in all residual images, but due to the complex artefact structure the remaining scattered light is of the order of the noise level that was estimated for the disc.

In all semi-dynamical models applied to the scattered light observations, the grain size distribution inferred from SED fitting was used. Here, a minimum grain size of 5.14$\,\mu\hbox{m}$ was assumed while the blow out grain size for 49 Cet lies at 2.9$\,\mu\hbox{m}$ using pure astronomical silicate. The grains between 2.9$\,\mu\hbox{m}$ and 5.14$\,\mu\hbox{m}$ move on bound orbits around their host star and would therefore significantly contribute to the flux density in scattered light if they were present in the disc.

Thus, we generated an additional model assuming the best fitting parameters of the single power law radial profile ($p=0.6$) and a size distribution ranging from the blow out grain size of 2.9$\,\mu\hbox{m}$ to the upper cut-off size of 15$\,\mu\hbox{m}$. The result can be seen in Fig. 13.

As suspected, the highly eccentric grains lead to a broader dust distribution and thus, the flux density in the outer disc regions increases comparable to the model using $5.14\,\mu\hbox{m}$ as minimum grain size. Furthermore, the forward scattering becomes stronger so that the disc signal is overestimated close to the star leading to an increased $\chi^{2}_{\text{red}}$ of $2.27$ in comparison to the models using the SED size distribution where we found $\chi^{2}_{\text{red}}=2.11$. In an extended study we analysed scattered light models using a grid of different minimum grain sizes between the blow-out limit (2.9$\,\mu\hbox{m}$) and the best-fit size of the SED (5.14$\,\mu\hbox{m}$). The model using a size of 5$\,\mu\hbox{m}$ delivers the best result with $\chi^{2}_{\text{red}}=2.11$. Thus, the semi-dynamical models give us an additional opportunity to estimate the minimum dust grain size besides a pure SED fit.

6.4 Scattering properties and impact on our results

An explanation for the residuals, although being close to the estimated noise level, might be given by the simplified model assumption considering the dust grains to be compact spheres in order to use Mie-theory. It is well established that the optical properties, such as absorption and scattering efficiencies or asymmetry parameter, strongly depend on the particle shape which can lead to significant deviations from the Mie values assumed (e.g., Schuerman et al., 1981; Weiss-Wrana, 1983; Mugnai & Wiscombe, 1986). Hence, in comparison to non-spherical grains, the Mie theory approach overestimates the fraction of forward scattering for larger particles (e.g., Arnold et al., 2018). In order to stay in an appropriate grain size regime, we estimated an upper cut-off grain size of $15\,\mu\hbox{m}$ for particles still significantly contributing to the measured scattered light and omitted all larger particles (see Sec. 4.3, Zubko, 2013).

To weaken the effect of forward scattering, many debris disc studies use the Henyey-Greenstein approach (Henyey & Greenstein, 1941) where the asymmetry parameter, $g$, is fixed to a specific value between -1 and 1 to give an integrated scattering phase function independent of particle composition and sizes.

In the first scattered light study of 49 Cet, Choquet et al. (2016) gave a best fit value of $g=0.1$. To check whether the remaining emission is caused by the overestimation of forward scattering we calculated a set of models based on the single power law radial distribution and fixed the asymmetry parameter for each fiducial disc. In this approach, the $g$ parameter is varied between -0.9 and 0.9 in steps of 0.1. We found a best fit model with $g=0.1$ which is in agreement with Choquet et al. (2017). The result is shown in Fig. 14.

Since the amount of forward scattering is reduced in the Henyey-Greenstein approach, we could reduce the overestimation of scattered light close to the star. Nevertheless, due to the higher particle eccentricities (see Appendix B for more details) in comparison to the discs assuming Mie theory, we get an overestimation of flux density in the outer disc region by assuming the radial distribution index to be $p=0.6$.

At this point, we have to give a word of caution about the usage of the Henyey-Greenstein approach together with our semi-dynamical dust modelling since it is an approach often used in former studies (e.g., Lee & Chiang, 2016; Esposito et al., 2016; Olofsson et al., 2016). The MODERATO code uses Mie theory to calculate the optical parameters of the generated dust grains. However, it also serves as basis to compute the radiation pressure parameter, $\beta$ (see Eq. 10), and therefore the dust particle orbits. It could be shown that $\beta$ depends as strongly on the shape of the dust grains as the optical parameters and thus, can easily differ by 75% comparing equivalent surface areas of non-spherical and spherical grains (e.g., Schuerman et al., 1981). By fixing the asymmetry parameter, the optical properties of the dust material assumed are tampered and hence, using the Henyey-Greenstein approach is inconsistent with the dust population inferred from Mie theory. Furthermore, to calculate the flux density of the scattered light not only the asymmetry parameter, but the scattering and extinction efficiencies are needed. These parameters are not provided by the Henyey-Greenstein approach and thus, have to be taken from Mie calculations or comparable methods which enlarges the already existing inconsistency. A more detailed comparison between the Henyey-Greenstein approach and Mie theory can be found in Appendix B.

In contrast to Mie and Henyey-Greenstein theory, a more sophisticated approach to model the scattered light might be given in the Discrete Dipole Approximation (DDA, e.g., Draine & Flatau, 2013).

7 Summary

We performed a multi-wavelength study of the debris disc around 49 Cet focusing on the dust component. We presented new scattered light data obtained with the SPHERE/IRDIS instrument on the VLT in H23 dual band and Y band filters and used thermal emission data from a former ALMA study at $850\,\mu\hbox{m}$. The H23-band detection was weak in comparison to the Y-band data so we focused our analysis on the Y-band data. For future SPHERE observations of comparable debris discs we therefore suggest using a broad band filter for stronger detections.

In both wavelength regimes the disc radius was found to be $\sim 280$ au. While for the ALMA detection a clear decrease in surface brightness could be observed, the extent of the disc in scattered light is limited by the noise level and might therefore be underestimated.

In a first modelling step, we fitted the SED of 49 Cet to obtain information on the grain size distribution. We inferred a dominant particle radius of 5.14$\,\mu\hbox{m}$ and a size distribution index of 3.77, suggesting that the disc is dynamically excited.

We computed semi-dynamical models using three different radial surface density profiles. By comparing our models to the observations in thermal emission and scattered light we found that all three can reproduce both SPHERE and ALMA observations, but that the double power law profile achieves the best fit results. Hence, we assume that the small grains ($\sim 5\,\mu\hbox{m}$) dominant in scattered light and large particles dominant in thermal emission stem from the same planetesimal belt and are mainly influenced by radiation pressure and collisions. Since transport processes, such as Poynting-Robertson drag, have greater effect on small particles we would expect differences in the radial dust distribution of small and large grains. However, these differences might only be observable in the inner disc region between the planetesimal belt and the star (Kennedy & Piette, 2015). Since our observations are not sensitive for this inner region of the 49 Cet debris disc we cannot draw any conclusions on such inward directed transport processes.

Acknowledgements

NP is thankful for fruitful discussions with Torsten Löhne and Alexander Krivov. AM acknowledges support from the Hungarian National Research, Development and Innovation Office NKFIH Grant KH-130526. JO acknowledges financial support from the ICM (Iniciativa Científica Milenio) via the Núcleo Milenio de Formación Planetaria grant, from the Universidad de Valparaíso, and from Fondecyt (grant 1180395) This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 716155 (SACCRED). FMe acknowledges funding from ANR of France under contract number ANR-16-CE31-0013. AZ acknowledges support from the CONICYT + PAI/ Convocatoria nacional subvención a la instalación en la academia, convocatoria 2017 + Folio PAI77170087. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF–Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). We also acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU in France. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). It has also been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). MRM, HMS, and SD are pleased to acknowledge this financial support of the SNSF. Finally, this work has made use of the the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice) and Observatoire de Paris/LESIA (Paris) and is supported by a grant from Labex OSUG@2020 (Investissements d’avenir - ANR10 LABX56). We thank P. Delorme and E. Lagadec (SPHERE Data Centre) for their efficient help during the data reduction process.

Appendix A Theoretical background

The orbits of dust particles are altered by different mechanisms, such as collisions, Poynting-Robertson drag, and stellar radiation pressure is one of the strongest processes.

A.1 Radiation pressure

As shown in earlier studies (e.g., Burns et al., 1979; Wyatt et al., 1999) the orbital parameters depend on the radiation pressure parameter, $\beta$, given as

[TABLE]

Here, $L_{\text{star}}$ and $M_{\text{star}}$ are the stellar luminosity and mass, $G$ the gravitational constant, $c$ the speed of light, $Q_{\text{pr}}$ the radiation pressure efficiency averaged over the stellar spectrum, $\varrho$ the bulk density of the dust material and $s$ the grain radius. If we assume a dust particle is released from a planetesimal possessing the orbital parameters semi-major axis, $a_{\text{p}}$, eccentricity, $e_{\text{p}}$ and true anomaly $f_{\text{p}}$, then the orbit parameters of the dust grain can be calculated with

[TABLE]

[TABLE]

[TABLE]

where $\varpi$ is the longitude of pericentre. To compute the $\beta$ parameter, we calculate $Q_{\text{pr}}$ using Mie theory, where the particles are assumed to be compact spheres (Bohren & Huffman, 1983).

In Fig. 15 on the left side the $\beta$ parameter is shown as a function of grain size. To calculate $\beta$, we use the stellar properties of 49 Cet and astronomical silicate (Draine, 2003) with a bulk density of 3.3 g/cm3. On the right side the particle eccentricity is given as a function of $\beta$. We use Eq. 12 and assume $e_{\text{p}}$ and $f_{\text{p}}$ to be equal to zero. In this case, the eccentricity turns 1 when $\beta$ is 0.5 and thus, grains below a value of $\beta=0.5$ stay in bound orbits around the star, while for values larger than 0.5 the particles are expelled from the stellar system on either parabola or hyperbola orbits. In general, $\beta$ is increasing with decreasing grain size. Therefore, the so-called blow-out grain size for 49 Cet is 2.9$\,\mu\hbox{m}$ meaning, that grains smaller than this value are expelled from the stellar system.

Furthermore, the smaller the grains the larger the eccentricity gets. As a consequence, small grains spend more time in the apocentre region compared to large particles and thus, it is expected that by tracing small grains in observations we would see larger disc radii than by tracing large grains.

To calculate the dust disc we generate at first the parent belt using the surface number density. In many cases a size independent power law is assumed:

[TABLE]

where the parameters are the same as described in Sec. 4.2. Now we have a population of planetesimals at different radii and angles releasing the dust grains. The planetesimals give us the parameters necessary ($a_{\text{p}}$, $e_{\text{p}}$, $f_{\text{p}}$) to calculate the particle orbits ($a_{\text{d}}$, $e_{\text{d}}$).

A.2 Tracing different grain sizes

It is necessary to investigate the absorption and scattering efficiencies for each grain size and wavelength to get information on the emission of the particles. For thermal emission the absorption efficiency of the grains is of importance while for scattered light emission it is the scattering efficiency. In Fig. 16 the absorption efficiency is depicted for a wavelength of 850$\,\mu\hbox{m}$ and the scattering efficiency for 1.04$\,\mu\hbox{m}$. These are the wavelengths used by SPHERE and ALMA for the 49 Cet observations used in this study.

Grains with a size of $\sim 0.5\,\mu\hbox{m}$ posses the highest scattering efficiency while $\sim 200\,\mu\hbox{m}$ sized particles show the largest absorption efficiency at the observed wavelengths. This means, we would expect to trace the small grains best with scattered light emission in the near infrared and the large grains with thermal emission in the sub-mm region. However, taking the radiation pressure into account, we expect to see grains larger than the blow-out size.

In Fig. 17 fiducial discs are given for thermal emission and scattered light. We use the 1.04$\,\mu\hbox{m}$ and 850$\,\mu\hbox{m}$ as observational wavelengths. Furthermore, we assume two cases. The first being that all discs possess the same mass and the second that all discs contain the same number of particles.

The scattered light images at 1.04$\,\mu\hbox{m}$ are shown in the first two rows, where the first row represents discs of the same mass and the second row discs with the same particle number. We see in the first case that the discs get fainter for an increasing grain size and decrease in radius. The latter fact is mainly caused by stellar radiation pressure. In the second case, where the particle number stays constant the discs get brighter.

One reason is the decreasing scattering efficiency. However, it stays nearly the same for grains larger than 10$\,\mu\hbox{m}$ and thus, the second reason is the constant disc mass leading to a smaller number of large grains. Another point is the wavy structure of the emission is caused by the phase function which depends on the grain size, the observation wavelength and the scattering angle. Considering the thermal emission images it is obvious that the disc with dust made of 100$\,\mu\hbox{m}$ grains possesses the largest emission at 850$\,\mu\hbox{m}$ compared to 5, 10 and 1000$\,\mu\hbox{m}$. The absorption efficiency decreases for particles smaller than 100$\,\mu\hbox{m}$. Hence, even the larger number of small particles cannot compensate the lower absorption.

Appendix B Comparison of Mie theory and Henyey-Greenstein approach

As shown by Equation 10, the radiation pressure parameter, $\beta$, is directly proportional to the radiation pressure efficiency, $\overline{Q}_{\text{pr}}(s)$, averaged over the stellar spectrum which is given by

[TABLE]

Here, $Q_{\text{pr}}$ describes the radiation pressure efficiency depending on grain radius, $s$, and wavelength, $\lambda$, while $F_{\lambda}$ gives the stellar flux density. For each grain and wavelength, $Q_{\text{pr}}$ is calculated by

[TABLE]

where $Q_{\text{ext}}$ and $Q_{\text{sca}}$ are the extinction and scattering efficiencies and $\langle\cos(\vartheta)\rangle$ the asymmetry parameter, depending on the scattering angle, $\vartheta$ (Bohren & Huffman, 1983). This parameter is also called $g$. All the parameters used to calculate $Q_{\text{pr}}$ depend on the grain size and the wavelength.

Using the Henyey-Greenstein approach (HG), the asymmetry parameter is fixed and therefore independent of grain size and wavelength, while all other parameters are still calculated by Mie theory or comparable methods. Since the extinction and scattering efficiencies as well as the $g$ parameter depend on the optical constants, i.e. refractive indices, of the dust composition used, the fixation of $g$ resembles a variation of the dust material.

A comparison of $g$ values given by Mie and HG is presented in Fig. 18. Here, $g$ is set to 0.5 resembling a best-fit value for different studies using the HG approach (e.g., Millar-Blanchaer et al., 2015; Olofsson et al., 2016; Engler et al., 2017, 2018). The values from Mie show strong deviations from this value. For 0.1$\,\mu\hbox{m}$-sized particles, $g$ decreases to values close to zero for wavelength larger than $1\,\mu\hbox{m}$ while for larger grains this decrease is moved to longer wavelengths.

In Fig. 19 we present the radiation pressure efficiency as a function of wavelength for different grain radii. Depending on the particle size, HG and Mie lead to similar results for wavelengths roughly longer than $\sim 3\times s$. For shorter wavelengths, the differences between the methods can easily reach 50%, although these deviations decrease for smaller particles (assuming $g=0.5$).

In the next step we compare the radiation pressure efficiencies averaged over the stellar spectrum for both, HG and Mie (Fig. 20). Here, the value calculated by HG is divided by the value inferred with Mie theory and given as a function of grain radius. Using three different $g$ values between 0.1 and 0.9 for HG, $\overline{Q}_{\text{pr}}$ is overestimated for both, $g=0.1$ and 0.5 with differences of up to $\sim 100\%$ assuming $g=0.1$ and $\sim$50% assuming $g=0.5$. In case of $g=0.9$, the deviation of $\overline{Q}_{\text{pr}}$ is close to zero for grains larger than $\sim 2\,\mu\hbox{m}$ while we underestimate $\overline{Q}_{\text{pr}}$ for smaller particles.

The parameter $\beta$ changes by the same order of magnitude as $\overline{Q}_{\text{pr}}$ (see Fig. 15). Hence, similar-sized particles possess different eccentricities for both methods . For example, if we assume a grain with a certain size $s$ for which $\beta=0.1$ using Mie theory, we get $\beta$ for HG twice as large assuming $g=0.1$ (see Figs. 15, 20). This leads to an increase of the particle’s eccentricity by a factor of 2.25 compared to Mie, using Eq. 12 and a circular parent belt. For larger $\beta$ this effect is even stronger due to the non-linear relation between $\beta$ and grain eccentricity (see Fig. 15). Only for the smallest particles ($s\sim 0.01\,\mu\hbox{m}$) HG and Mie reach similar values. Hence, by applying the Henyey-Greenstein approach the dust population inferred with Mie theory is altered by HG.