Nicole Pawellek
email: pawellek@mpia.de

Attila Moór

Julien Milli

Ágnes Kóspál

Johan Olofsson

Péter Ábrahám

Miriam Keppler

Quentin Kral

Adriana Pohl

Jean-Charles Augereau

Anthony Boccaletti

Gaël Chauvin

Élodie Choquet

Natalia Engler

Thomas Henning

Maud Langlois

Eve J Lee

François Ménard

Philippe Thébault

Alice Zurlo

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

Keywords: infrared, stars -circumstellar matter -stars, individual (49 Cet)

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

INTRODUCTION

Circumstellar discs around young stars are natural byproducts 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 ∼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 ∼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 16L and an effective temperature of 9000 K (Roberge et al. 2013). Its distance is given by the new Gaia data release as 57.0±0.3 pc (Gaia Collaboration et al. 2016Collaboration et al. , 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 discs 1 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 β Pic shows an extent of ∼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;[START_REF] Dent | [END_REF]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 3. 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.

OBSERVATIONS AND DATA REDUCTION

Broad band SPHERE data

We observed the disc around 49 Cet in the programme 198.C-0209(N) on 19 th November 2016 for 1.8 hours on source with the SPHERE/IRDIS instrument of the VLT [START_REF] Beuzit | Circumstellar Dust Disks and Planet Formation[END_REF]Dohlen et al. 2008). The observations were carried out using the broad band Y filter with a central wavelength of 1.04 µ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, [START_REF] Marois | [END_REF]. 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/arcsec 2 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/arcsec 2 , 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 arcsec 2 . The pixel scale of IRDIS is 0.01225 arcsec/pixel [START_REF]Ground-based and Airborne Instrumentation for Astronomy VI[END_REF]. 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.

Narrow band SPHERE data

We also observed 49 Cet on 23 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 [START_REF] Delorme | Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics[END_REF]. The obtained reduced master cubes were then utilised as input for high contrast imaging post-processing performed by the Speckle Calibration (SpeCal) package [START_REF] Delorme | Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics[END_REF]Galicher et al. 2018). The reduced image in ADU is depicted in Fig. 2.

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 (80 au) and ∼ 4.9 (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 • 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 • 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 ∼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 ∼ 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 • ±1.0 • . 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.

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 ∼ 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 ∼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 ∼130 au and ∼170 au in the eastern and western direction in both SPHERE images.

Radial profile of thermal emission

In order to compare the brightness profile of the nearinfrared 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 µm continuum images for 49 Cet with a beam size of 0.47 × 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 ∼280 au and the western side up to ∼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 ∼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.

MODELLING

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

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.

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 [START_REF] Brott | The Three-Dimensional Universe with Gaia[END_REF]. The dust composition is assumed to be astronomical silicate (Draine 2003) with a bulk density of 3.3 g/cm 3 .

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 µ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 powerlaw for both the radial and the size distribution of the dust using the surface number density N(r, s):

N(r, s) ∼ r r 0 -p × s s 0 -q
.

(1)

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, Σ(r, s), by Σ(r, s)ds = πs 2 × N(r, s)ds.

(2)

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 min = 70 au and r max = 298 au to fit the SED.

We assume the grain sizes to lie between a minimum and a maximum value, s min and s 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 min , the size distribution index, q and the amount of dust, M dust , for particles between s min and s max . We infer the amount of dust by using the bulk density and the dust volume V:

M dust = × V = × 4π 3 s max ∫ s min 2π r max ∫ r min N(r, s) r dr s 3 ds. (3)
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 ∼ 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

R disc R blackbody = A × L L sun B .
(4) For pure astronomical silicate Pawellek (2017) give the parameters as A = 6.49 ± 0.86 and B = -0.37 ± 0.05. Using these values and the solar luminosity of 16L we find an estimated "true" location of the warm component of ∼ 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 min = 5.14 ± 0.76 µm and the size distribution index q = 3.77 ± 0.05. These parameters are applied to all semidynamical 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 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 d = 8.8 × 10 -4 .
The error bars give the uncertainties of the fit parameters and are inferred in the following way. We start from the 10000 thermal emission, fixed 15 scattered light, fixed q 3.77 ± 0.05 . . .

f d 8.8 × 10 -4 . . . T dust [K] 65 ± 3 . . . T warm [K]
125 ± 2 . . . χ 2 red 1.03 . . . Notes: The blow out limit is given for pure astronomical silicate as 2.9 µm. The different maximum grain sizes are explained in Section 4.3. For the dust temperature of both components we take the wavelength of the SED peak and calculate the temperature by using Wien's displacement law.

Figure 6. SED of 49 Cet using pure astronomical silicate with a bulk density of 3.3g/cm 3 for the outer and inner component.

The thin black solid line represents the stellar photosphere, the red solid line shows the warm inner ring, the blue lines the outer belt.

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

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 (≥ 20 µm) since the grains are assumed to be compact spheres (e.g., Schuerman et al. 1981;[START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF]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:

α = 90 • -i, (5)
which is ≈ 10 • 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 πs 2 and therefore, the Rayleigh criterion provides the minimum angle observable for the grain at a certain wavelength λ:

α = 1.22 × λ 2 × s = 70 • × λ 2 × s . ( 6
)
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:

s ≈ λ × 70 • 2 × (90 • -i) = 1.04 µm × 70 • 2 × 10 • = 3.5 µm. (7)
This shows, that the scattered light contribution of particles larger than 3.5 µ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π×s/λ > 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 ∼ 5 µ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 µ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 χ 2 -minimisation assuming that an ideal residual image should expose no emission in any pixel. The χ 2 -parameter is then computed for each pixel by

χ 2 = N pixel i=1 F i F noise 2 (8)
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 noise = 0.056 mJy/beam (Hughes et al. 2017). The number of pixels for the thermal emission maps is 1024 × 1024 and for the scattered light images 255 × 255.

ANALYSIS

Thermal emission

The thermal emission image observed by Hughes et al. (2017) with ALMA at 850 µm is compared to the semidynamical 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. 789.

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 χ 2 of χ 2 red = 1.4. The grain parameters derived by SED-fitting are applied to our model, where the minimum grain size is 5.14 µ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 χ 2 red = 1.18. After subtracting the model from the observed data, there are some 2σ residual structures left throughout the disc. However, they do not seem to be of systematic origin.

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 χ 2 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 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 ring = 140 au. The parameter χ 2 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.

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 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 χ 2 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 trans = 110 au, p 1 = -1.0 and p 2 = 1.0 and leads to χ 2 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σ remnants are left.

Comparing the models

By comparing the different radial profiles we find that all models achieve a similar χ 2 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

BIC = χ 2 + J × log e (N), (9)
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 × log e (∆BIC), where ∆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 χ 2 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.

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. 101112. 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 • 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 χ 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 oversubtraction 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 χ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.

DISCUSSION

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 [START_REF] Marois | [END_REF]) 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 χ 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 bestfitting 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.

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 µ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 ± 0.76 µm was found as dominant grain radius. The ratio of the minimum particle size to the inferred blow-out limit of 2.9 µ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 ± 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 µm are expelled from the 49 Cet system [START_REF] Burns | [END_REF][START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF]. Fig. A2 shows that grains smaller than 1 µ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 µ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.

Model images

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 µm) not sensitively traced at submm wavelengths.

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 semimajor 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 µm was assumed while the blow out grain size for 49 Cet lies at 2.9 µm using pure astronomical silicate. The grains between 2.9 µm and 5.14 µ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 µm to the upper cut-off size of 15 µ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 µ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 χ 2 red of 2.27 in comparison to the models using the SED size distribution where we found χ 2 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 µm) and the best-fit size of the SED (5.14 µm). The model using a size of 5 µm delivers the best result with χ 2 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.

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 µ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, β (see Eq. A1), and therefore the dust particle orbits. It could be shown that β 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).

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 µ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 ∼ 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 µ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 (∼ 5 µ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. 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)(2010)(2011)(2012) and grant number 312430 for FP7 (2013)(2014)(2015)(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), PY-THEAS/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.

β is increasing with decreasing grain size. Therefore, the socalled blow-out grain size for 49 Cet is 2.9 µ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:

N(r) = N 0 × r r 0 -p , ( A5
)
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 p , e p , f p ) to calculate the particle orbits (a d , e d ).

A2 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. A2 the absorption efficiency is depicted for a wavelength of 850 µm and the scattering efficiency for 1.04 µm. These are the wavelengths used by SPHERE and ALMA for the 49 Cet observations used in this study. Grains with a size of ∼ 0.5 µm posses the highest scattering efficiency while ∼ 200 µ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. A3 fiducial discs are given for thermal emission and scattered light. We use the 1.04 µm and 850 µ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 µ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 µ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 µm grains possesses the largest emission at 850 µm compared to 5, 10 and 1000 µm. The absorption efficiency decreases for particles smaller than 100 µ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 A1, the radiation pressure parameter, β, is directly proportional to the radiation pressure efficiency, Q pr (s), averaged over the stellar spectrum which is given by

Q pr = ∫ Q pr (s, λ) × F λ dλ ∫ F λ dλ . (B1)
Here, Q pr describes the radiation pressure efficiency depending on grain radius, s, and wavelength, λ, while F λ gives the stellar flux density. For each grain and wavelength, Q pr is calculated by

Q pr (s, λ) = Q ext (s, λ) -Q sca (s, λ) × cos(ϑ) (s, λ), (B2)
where Q ext and Q sca are the extinction and scattering efficiencies and cos(ϑ) the asymmetry parameter, depending on the scattering angle, ϑ [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF]. This parameter is also called g. All the parameters used to calculate Q 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. B1. 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. 2017Engler et al. , 2018)). The values from Mie show strong deviations from this value. For 0.1 µm-sized particles, g decreases to values close to zero for wavelength larger than 1 µm while for larger grains this decrease is moved to longer wavelengths.

In Fig. B2 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 ∼ 3 × 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. B3). 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, Q pr is overestimated for both, g = 0.1 and 0.5 with differences of up to ∼ 100% assuming g = 0.1 and ∼50% assuming g = 0.5. In case of g = 0.9, the deviation of Q pr is close to zero for grains larger than ∼ 2 µm while we underestimate Q pr for smaller particles.

The parameter β changes by the same order of magnitude as Q pr (see Fig. A1). Hence, similar-sized particles possess different eccentricities for both methods . For example, if we assume a grain with a certain size s for which β = 0.1 using Mie theory, we get β for HG twice as large assuming g = 0.1 (see Figs. A1,B3). This leads to an increase of the particle's eccentricity by a factor of 2.25 compared to Mie, using Eq. A3 and a circular parent belt. For larger β this effect is even stronger due to the non-linear relation between β and grain eccentricity (see Fig. A1). Only for the smallest particles (s ∼ 0.01 µ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.

Figure 1 .

 1 Figure 1. Classical ADI-reduced scattered light Y-band image of 49 Cet. The white dashed line shows the inferred disc extent of 280 au, the blue solid line the location of the surface brightness peak at 140 au. The numbers 1 and 2 indicate the residuals . Further explanations are given in the text.

Figure 2 .

 2 Figure 2. Classical ADI-reduced scattered light H-band image of 49 Cet. The white dashed line shows the inferred disc extension of 200 au, the blue solid line the location of the surface brightness peak at ∼140 au.

Figure 3 .

 3 Figure 3. Surface brightness as a function of radial distance to the star for the SPHERE observations. The upper panel shows the Y-band data, the lower the H-band data. The red dashed line depicts the radial profile of the disc, the black solid line the 1σ noise level.Bright blue filled areas in both panels represent the location of symmetric peaks in western and eastern direction. The grey filled area shows the inner region of 80 au where the noise level is comparable to the signal.

Figure 4 .

 4 Figure 4. Surface brightness as a function of radial distance. The thick black solid line shows the observational data for ALMA, the black dashed line for SPHERE H-band normalised to the ALMA data. The blue dotted line shows the disc with single power law radial profile, the green dash-dotted line a two -power law radial profile and the red dash-double dotted line a single power law including a narrow dust ring. Blue shaded areas represent symmetric peaks in western and eastern direction at the same position as in Fig. 3. The grey filled area shows the inner region of 80 au.

Figure 5 .

 5 Figure 5. Surface density profiles for the parent belts inferred from semi-dynamical models. Red shows the single power law; green the power law including a narrow ring and blue the double power law model. The vertical green dotted line and the blue dashed line give the position of the planetesimal ring and the transition radius, respectively. The different radial distribution indices used for the models are depicted at the corresponding parts of the profiles. The noise is caused by the finite number of planetesimals in the models.

Figure 7 .

 7 Figure 7. Modelling results of 49 Cet. Left: ALMA image at 850 µm.Middle: Best fitting model using a single power law radial distribution with p = 0.6 generated with MODERATO. A grain size distribution with s min = 5.14 µm, s max = 1.00cm and q = 3.77 was used to generate the disc. Right: Subtraction of the model disc from the observational data. The contour levels show 2, 4 and 6 times the σ level given as 0.056 mJy/beam inHughes et al. (2017). The dotted contours give the negative scale, the white dashed area the beam size. The data are normalised to their maximum. The value of χ 2 red is 1.18.

Figure 8 .

 8 Figure8. Same as Fig.7, but now using the single power law radial distribution with p = 0.5 including an additional narrow ring at r ring = 140 au. For this model χ 2 red = 1.16.

Figure 9 .

 9 Figure 9. Same as Fig. 7, but now using the double power law radial profile with a transition radius of r trans = 130 au and radial distribution indices of p 1 = -1.0 and p 2 = 1.0. The best fit χ 2 red is 1.15.

Figure 10 .

 10 Figure 10. Modelling results of 49 Cet. Left: ADI-reduced Y-band image as shown in Fig. 1 excluding the region within 80 au. Middle: model generated with MODERATO assuming a single power law as radial distribution with p = 0.6. Right: Subtraction of the model disc from the observational data. The model leads to χ 2 red = 2.11.

Figure 11 .

 11 Figure 11. Same as Fig. 10, but now using the single power law radial distribution with p = 0.5 including an additional narrow ring at r ring = 140 au. For this model χ 2 red = 2.16.

Figure 12 .

 12 Figure 12. Same as Fig.10, but now using the double power law radial profile with a transition radius at r trans = 130 au and radial distribution indices of p 1 = -1.0 and p 2 = 1.0. The best fit χ 2 red is 1.99.

Figure 13 .

 13 Figure13. Same as Fig.10, using the double power law radial profile and a size distribution between 2.9 µm and 15 µm. The radial distribution index is p = 0.6.

Figure 14 .

 14 Figure 14. Same as Fig. 10, using the single power law radial profile with a fixed asymmetry parameter of g = 0.1.

Figure A1 .

 A1 Figure A1. Left: Radiation pressure parameter as a function of grain size calculated for 49 Cet (solid black line). Solid red line: using Mie theory; dash-dotted blue line: using Henyey-Greenstein approach with g = 0.1; black dashed line shows β = 0.5. Right: particle eccentricity as a function of the radiation pressure parameter (solid black line) assuming a circular orbit for the parent body. The dashed lines show β = 0.5 and e = 1.0.

Figure A2 .

 A2 Figure A2. Scattering and absorption efficiency as a function of grain size. The blue dashed line shows the scattering efficiency at a wavelength of 1.04 µm and the red dash dotted line the absorption efficiency at 850 µm. The dust composition is pure astronomical silicate.

Figure A3 .

 A3 Figure A3. Model discs for different grain sizes. From left to right: The dust is made of a single grain size of 5 µm, 10 µm, 100 µm, 1000 µm. The first row gives the scattering emission of the discs of a constant mass at 1.04 µm. The second row shows the same as the first row, but now the discs have the same number of particles. The third row shows the discs with a constant mass in thermal emission at 850 µm. The bottom row shows the same as the third row, but now the discs have a constant number of particles.

Figure B1 .

 B1 Figure B1. Asymmetry parameter as function of wavelength.Comparison of Mie theory and Henyey-Greenstein approach using astronomical silicate and a fixed g of 0.5. Black solid line shows the fixed asymmetry parameter; colour dashed lines show g for different grain radii (red dashed: 0.1 µm, green dotted: 1.0 µm, blue dash-dotted: 10 µm).

Figure B2 .

 B2 Figure B2. Radiation pressure efficiency as function of wavelength for the same grain radii as in Fig. B1. The solid lines shows the result for Mie theory, the dashed lines for Henyey-Greenstein assuming g = 0.5.

Figure B3 .

 B3 Figure B3. Ratio of the averaged radiation pressure efficiencies of Henyey-Greenstein and Mie theory as a function of grain size for different values of g. Red dashed: g = 0.1; green dotted: g = 0.5, blue dash-dotted: g = 0.9; black solid: ratio of 1

Table 1 .

 1 Disc parameters inferred in different wavelength ranges.

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

Table 2 .

 2 Parameters for the radial distribution models. The abbreviation PL stands for power law. The parameter values of the first table part were taken from Hughes et al. (2017), corrected for the new distance of 57.1 pc, and used for the modelling of this work. In the second part of the table the values in column H17 stem from Hughes et al. (2017) and are corrected for the new distance of 57.1 pc. Parameters in the column Fit are the best fit results of this work.

	Parameter		Single PL	Double PL
			no ring	with ring
			H17	Fit H17	Fit H17	Fit
	i [ • ]		79.3		79.3	79.2
	P A [ • ]	108.7	108.8	108.3
	r min [au]	70		58	26
	r 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 ring [au]	range	. . .		. . . 90-150	. . .
		step size	. . .		. . .	10	. . .
		value	. . .		109	140	. . .
	r 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
	χ 2 red	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
	Notes:

Table 3 .

 3 Continuum flux density.

	Wavelength Flux density	Instrument	Reference
	[ µm]	[mJy]
	0.44 22260 ± 350	TYCHO	1
	0.55 21600 ± 200	TYCHO	1
	1.24 10180 ± 190	2MASS	2
	1.65 6300 ± 130	2MASS	2
	2.16 4373 ± 81	2MASS	2
	3.35 1990 ± 120	WISE	3
	4.60 1294 ± 49	WISE	3
	5.86	690 ± 70	IRS	4
	7.07	480 ± 50	IRS	4
	8.97	320 ± 30	IRS	4
	9.0	366 ± 13.5	AKARI/IRC	5
	10.8	250 ± 50	Keck/MIRLIN	6
	11.40	210 ± 20	IRS	4
	11.56	211 ± 21	WISE	3
	12.0	330 ± 66	IRAS	6
	12.50	200 ± 26	Keck/MIRLIN	8
	13.90	180 ± 20	IRS	4
	17.90	186 ± 25	Keck/MIRLIN	8
	18.0	199 ± 16	AKARI/IRC	5
	22.09	238 ± 24	WISE	3
	24.00	259 ± 10	Spitzer/MIPS	9
	25.0	380 ± 76	IRAS	7
	60.0 2000 ± 400	IRAS	7
	63.19 2090 ± 350	Herschel/PACS Spec	9
	70.00 2163 ± 151	Herschel/PACS	10
	71.42 1749 ± 123	Spitzer/MIPS	10
	72.84 1950 ± 320	Herschel/PACS Spec	9
	78.74 1900 ± 310	Herschel/PACS Spec	9
	90.0 1776 ± 295	AKARI/FIS	5
	90.16 1880 ± 320	Herschel/PACS Spec	9
	100.0 1910 ± 380	IRAS	7
	100.0 1919 ± 134	Herschel/PACS	10
	145.54 1160 ± 180	Herschel/PACS Spec	9
	150	750 ± 500	ISO	4,11
	157.68	980 ± 130	Herschel/PACS Spec	9
	160.00 1066 ± 75	Herschel/PACS	10
	170 1100 ± 500	ISO	4,11
	250.00	372 ± 27	Herschel/SPIRE	9
	350.00	180 ± 14	Herschel/SPIRE	9
	450.0	125 ± 18	JCMT/SCUBA-2	12
	500.00	86 ± 9	Herschel/SPIRE	9
	850.0	17 ± 3	ALMA	4
	850.00	13.5 ± 1.5	JCMT/SCUBA-2	12
	1200.00	12.7 ± 2.8	IRAM	13
	1300.0	2.1 ± 0.7	CARMA	4
	1330	5.5 ± 0.7	ALMA/ACA	14
	9000.0 0.0251 ± 0.0055	VLA	15
	References: [1] -Høg et al. (2000); [2] -2MASS All-Sky Catalog
	of Point Sources; [3] -Wright et al. (2010); [4] -Hughes et al.
	(2017); [5] -AKARI All-Sky Survey Bright Source Catalog [6]
	-Jayawardhana et al. (2001); [7] -IRAS Faint Source Catalog;
	[8] -Wahhaj et al. (2007); [9] -Roberge et al. (2013); [10] -Moór
	et al. (2015) [11] -ISO; [12] -Holland et al. (2017); [13] -Bockelée-
	Morvan et al. (1994); [14] -Moór et al. (in prep) [15] -MacGregor
	et al. (2016);

position of the minimum χ 2 in the parameter space. New parameter values are generated and the resulting χ 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

Table 4 .

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

Parameter Best fit values Comment s min

[ µm]

5.14 ± 0.76 . . . s max

[ µm]

Table 5 .

 5 Comparison of thermal emission models

	No. Model	J Free parameters	χ 2 [×10 6 ] BIC [×10 6 ]
	1 SPL	2 scaling, p	1.24	1.24
	2 SPL + ring 3 scaling, p, r ring	1.22	1.22
	3 DPL	4 scaling, p 1 , p 2 , r trans	1.21	1.21
	Notes: SPL is the single power law model, SPL + ring the single
	power law including a narrow ring and DPL the double power
	law model. The number of data points assumed is the number of
		pixels: N = 1024 × 1024.

Table 6 .

 6 Comparison of scattered light models Comparing models χ 2 [×10 9 ] B The model numbers are given in Tab. 5. The number of data points assumed is the number of pixels of the scattered light image: N = 255 × 255.

	1, 2	8.58 4.8
	1, 3	8.58 6.2
	2, 3	8.58 4.8
	Notes:

MNRAS 000, 1-19 (2018)

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000, 1-19(2018)

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,

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.

A1 Radiation pressure

As shown in earlier studies (e.g., [START_REF] Burns | [END_REF]Wyatt et al. 1999) the orbital parameters depend on the radiation pressure parameter, β, given as

Here, L star and M star are the stellar luminosity and mass, G the gravitational constant, c the speed of light, Q pr the radiation pressure efficiency averaged over the stellar spectrum, 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 p , eccentricity, e p and true anomaly f p , then the orbit parameters of the dust grain can be calculated with

where is the longitude of pericentre. To compute the β parameter, we calculate Q pr using Mie theory, where the particles are assumed to be compact spheres [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF].

In Fig. A1 on the left side the β parameter is shown as a function of grain size. To calculate β, we use the stellar properties of 49 Cet and astronomical silicate (Draine 2003) with a bulk density of 3.3 g/cm 3 . On the right side the particle eccentricity is given as a function of β. We use Eq. A3 and assume e p and f p to be equal to zero. In this case, the eccentricity turns 1 when β is 0.5 and thus, grains below a value of β = 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,


Talha Ijaz

Ben Yang

Ruipu Wang

Jiazhe Zhu

Aftab Farrukh

Gong Chen

Gregory Franc

Yang Zhang
email: zhyangnano@ustc.edu.cn

Andre Gourdon
email: andre.gourdon@cemes.fr

Zhenchao Dong
email: zcdong@ustc.edu.cn

Self-decoupled tetrapodal perylene molecules for luminescence studies of isolated emitters on Au(111)

Self-decoupled tetrapodal perylene molecules were designed, synthesized and deposited on the Au(111) surface through electrospray ionization technique. Photoluminescence and lifetime measurements show that the chromophore groups of the designed molecules are well decoupled from the gold substrate. Molecule-specific emissions from both neutral and anionic molecules were observed in the preliminary electroluminescence measurements. The

emergence of significant emission when the tip is positioned at the molecular center suggests that there is a considerable vertical component of the transition dipole of the designed molecule along the tip axial direction. Our results may open up a new route for the realization of nano light sources and plasmonic devices based on organic molecules.

Scanning tunneling microscope (STM) induced molecular emission offers unprecedented opportunities for gaining insights into single-molecule characteristics by exploring optoelectronic phenomena occurring at the single-molecule level. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] A key issue to obtain single molecular electroluminescence in STM is to electronically decouple the molecular emitters from the metallic substrate since the fluorescence of the emitters directly adsorbed on the metal surface would be quenched due to the direct electron transfer between molecules and metal substrate. The most widely used approach to achieve electronic decoupling is to physically introduce an insulating spacer between the molecule and metal substrate, such as oxides, 1 halides, 3 thiol layers, 16 and molecular multilayers. 2 A less explored approach is the chemical route that involves chemical modifications by adding spacer groups to the emitting unit, which would allow studying a wide range of functionally tunable molecules directly adsorbed on metal substrates. [17][18][19][20] Indeed, the electronic self-decoupling effect was achieved previously by chemically adding tripodal anchors to the porphyrin chromophore. 18 However, a majority of porphyrin molecules with tripodal anchors still tend to adopt flat-lying configurations after deposition on a metal substrate and thus still suffer from severe fluorescence quenching. 18 Consequently, how to further optimize the chemical design of anchor supported molecules to achieve efficient decoupling, regardless of the adsorption configurations, remains to be a challenging issue to be solved. One other challenge is how to align the transition dipole of the emitters along the vertical direction that would favor strong plasmonic enhancement? In this work, we shall first present the chemical design, synthesis, and deposition of self-decoupled tetrapodal perylene (TP) molecules, which is better optimized for self-decoupling properties. Then, we shall investigate the optical properties of the TP molecules on Au(111) using photoluminescence (PL) and scanning tunneling microscope induced luminescence(STML) techniques. Figure 1(a) shows the chemical structure of our designed molecule. The tetrahedral geometry is adopted to explore the idea that no matter in what configuration molecules are adsorbed on the surface there would always be one emitter pointing in the upward direction, as shown in the inset of Fig. 1(a). To this end, the tetraphenylmethane is selected to act as a tetrahedral backbone for the targeted tetrapodal molecule. Perylene diimide (PDI) molecules are chosen due to their high photochemical stability and quantum yield. PDI functions as both emitter and spacer, though separately, while alkyl chains (i.e., two hexylheptyl chains for each PDI unit) serve as solubilizing moieties and anchoring groups.

The target TP molecule, namely, tetra-4-[N-(1-Hexylheptyl)perylene-3,4:9,10tetracarboxylicbisimid-N′-yl]phenylmethane, was obtained by a route similar to previous reference, 21 but by direct condensation of four PDI precursors with tetrakis(4aminophenyl)methane in presence of zinc acetate/imidazole rather than by conversion to the tetraformamide. The deposition of isolated TP molecules on metal surfaces is challenging because of its large molecular weight and low thermal stability, making it impossible to transfer by thermal sublimation under ultrahigh vacuum (UHV) conditions. In order to overcome such difficulty, we employ the electrospray ionization (ESI) technique to deposit these molecules onto the Au(111) metal surface. 22,23 As shown in Fig. 1(b), a commercial ESI system (Molecular Spray UHV4) with three-stage pumping is mounted to the loadlock chamber. A solution with a concentration of 34 µM was prepared by dissolving TP molecules in solvents consisting of a mixture of dichloromethane and methanol with a volume ratio of 4: 1. This solution was passed through a stainless steel emitter capillary held at 1.8 kV. This leads to the formation of spray outside the vacuum system. The spray enters the pre-vacuum chamber through an aperture, then passes through differentially pumped vacuum stages, and finally deposits onto atomically flat Au(111) surfaces. 22 The vacuum of the loadlock chamber varies slightly during the ESI deposition for 300 seconds, typically from a base pressure of 1.510 -6 Torr to 2.010 -6 Torr.

All of the STM, PL, and STML experiments were performed in the observation chamber of a low-temperature UHV-STM (Unisoku) with a base pressure of 810 -11 Torr at 78 K or 7 K, unless otherwise noted. STM imaging was performed at constant-current mode with sample biased. A gold film with a thickness of 200 nm was prepared by thermal deposition of gold on a freshly cleaved mica surface. Before ESI deposition, the atomically flat Au(111) surface was achieved by cycles of sputtering and annealing. Electrochemically etched silver tips used in the experiments were cleaned by electron bombardment and argon ion sputtering, followed by tip indentations to achieve desired nanocavity plasmon (NCP) modes for resonance enhancement. The optical setups were detailed in previous references. 7,9,24 For lifetime measurements, a 532 nm pulsed laser with a 35 ps pulse width was used and synchronized with the time-correlated single-photon counting (TCSPC) setup. In order to examine the self-decoupling efficiency of the designed TP molecules on Au(111), we first carried out photoluminescence measurements to reveal their emission feature. Figure 2(a) shows a typical PL spectrum acquired from the TP molecules adsorbed on Au(111), featuring three emission peaks at 540 nm, 585 nm and 637 nm, respectively. These peak positions are similar to the spectrum acquired from TP molecules in solution, as shown in Fig. 2(b), and are also consistent with the emission peak values of PDI monomer molecules reported in the literature. [26][27][28] In addition, according to the mirror-image symmetry between the absorption and fluorescence spectra, the emission peak at 540 nm can be assigned to the 00 line, and the peaks at 585 nm and 637 nm can be assigned to the associated vibronic progression with an energy difference of 0.17 eV. The observation of clear TP molecular emission in PL suggests that the direct electron transfer between the emitter and the substrate is blocked and as a result, the emitter group (top PDI molecular unit) is indeed well decoupled from the underneath gold substrate. In other words, the chemical design strategy via a tetrahedral tetrapod works well.

In order to probe the influence of both the local environment around the emitter and the metal substrate on the exciton decay dynamics of TP molecules, we also performed fluorescence lifetime measurements on the 00 emission of the TP molecules adsorbed on Au(111) and in solution, respectively, as shown in Fig. 2(c) and 2(d). The fluorescence decay of the TP molecules in solution can be well fitted by a single exponential function with a lifetime of 4.31(1) ns, which is very close to the PDI monomer lifetime of about 4 ns reported in previous literature. 27,29 Such an agreement suggests that the tetropodal chemical design does not generate additional nonradiative channels that would otherwise quench the molecular fluorescence. However, the fluorescence of the TP molecule on Au(111) exhibits a biexponential decay feature with the major component (83%) for a longer lifetime of 0.58(2) ns and the minor component (17%) for a shorter lifetime of 0.12(1) ns. Both lifetimes are much shorter than that observed for the TP molecules in solution, which suggests that the presence of the metal substrate opens up fast energy transfer channels. 30,31 The shorter lifetime (0.12 ns) might be assigned to the direct decay channel between the excited chromophore to the metal substrate through dipole-dipole interactions, while the relatively longer lifetime (0.58 ns) might be attributed to the indirect decay channel that is associated with the interaction between the top perylene emitter and the underneath spacer perylene units, the latter are close to the metal substrate and could then dissipate the transferred energy rapidly to the substrate. In this context, the dominant role of the 0. We also performed preliminary STML measurements on the TP molecules directly adsorbed on Au(111) to probe the electroluminescence behavior of individual TP molecules under tunneling electron excitations, as schematically illustrated in Fig. 3(a). Here, we would like to note that, unlike the constant observation of three characteristic emission peaks in the PL measurements for the TP molecules on Au(111) [Fig. 2(a)], STM induced luminescence acquired above individual TP molecules on Au(111) was found not only different from the PL spectrum, but also exhibiting a variety of different spectral profiles. Such distinctive luminescence behavior suggests that tunneling electrons (and tip approaching) can easily modify the molecular configuration or electronic states since TP is a complicated bulky molecule with long alkyl chains. The easily changing nature of the molecular configurations or electronic states under tunneling electron excitations makes it highly challenging to run STML measurements and to understand their emission features as well. For example, as shown in Fig. 3(b), we observed two sharp emission peaks at 593 nm and 641 nm, respectively, which are evidently different from the broad NCP emission measured on the nearby bare Au(111) surface [red curve in Fig. 3(b)] and thus suggest the molecular origin of the emission. In addition, the observed double peaks are found to resemble the 01 and 02 vibronic features of the PL spectrum of TP on Au(111) with a similar energy spacing of 0.16 eV despite small redshifts in the peak positions. Therefore, the double-peak STML spectrum is attributed to the fluorescence from a neutral TP molecule. The redshift of the peak positions in the STML measurement with respect to those in the PL measurements may originate from either configuration variations or the stark effect. [32][33][34] The reason for the absence of 00 emission is unclear, but might be related to the easily changing nature of molecular configurations upon electron excitations, which may favor vibronic emissions.

In addition to the double-peak spectral profile such as Fig. 3(b), STML spectra were often found to also exhibit sharp emission peaks in the long-wavelength region above 700 nm. As exemplified in Fig. 3(c), two sharp emission peaks at 856 nm and 923 nm were observed, which could be related to either anionic emissions from charged PDI units 35,36 or phosphorescence. 13,37 The relatively broad peak at ~700 nm might be related to the 03 emission of a neutral TP molecule in terms of the vibrational energy spacing of ~0.16 eV observed above, though the possibility stemming from molecule-modulated NCP emission cannot be completely excluded. 15,38 The observation of emissions from both neutral and charged TP molecules indicates that tunneling electrons not only can modify molecular configurations, but also can change molecular charging states. 11,13 In an effort to improve the configurational stability of the TP molecules, we further cooled the system to the liquid-helium temperature (~7 K at the sample). However, the configurational stability of the TP molecules does not seem to improve since various types of STML spectral profiles were still observed at such a low temperature even for very small currents down to 2 pA, which suggests that the TP molecular configurations are very susceptible to perturbation by tunneling electrons. As shown in Fig 3(d), two very sharp emission peaks at ~856 nm and 937 nm were observed, which can be attributed again to the anionic emission or phosphorescence of TP molecules. The origin of the broad emission around 700 nm with three minor peaks added on top (ranging from 632 nm to 725 nm) is still unclear but appears to contain mixed contributions from the TP molecule and NCP emission.

Notably, all the blue STML spectra shown in Fig. 3 were acquired when the tip was positioned above the center of individual TP molecules. The observation of significant molecule-specific emission intensity at the center position (stronger than or comparable to other molecular sites) is substantially different from the STML behavior on a flat-lying zinc phthalocyanine (ZnPc) molecule, where the intensity acquired at the ZnPc center is extremely weak because the ZnPc transition dipole is oriented horizontally and the associated total dipole symmetry of the whole tipmoleculesubstrate system tends to cancel out. 4,7 Such a great contrast in the STML intensity at the molecular center between a TP and a ZnPc suggests that the transition dipole of the TP molecule on Au(111) is likely to stand up straight or at least contain considerable vertical component along the tip axial direction. Such a vertical-dipole configuration may offer new opportunities for probing possible strong coupling between the NCP and molecular emitter.

In summary, we have designed and synthesized a tetrahedron-like tetrapodal perylene molecule that can suppress fluorescence quenching from metal substrates through the built-in self-decoupling units, regardless of adsorption configurations. We have also demonstrated that molecular electrospray technique can be used to deposit such large functional but thermally unstable molecules on Au(111) for the luminescence study of isolated molecules.

FIG. 1 .

 1 FIG. 1. (a) Chemical structure of the designed TP molecule, with the inset on the top-right showing the decoupling strategy. (b) Schematic of electrospray ionization for molecular deposition in vacuum. P1, P2, P3 shows differentially pumping stages and S represents the Au substrate. (c) A 3D view of STM topograph of single TP molecules on Au(111) (2.7 V, 2 pA, 27  22 nm 2 ). (d) Typical height profile of a single selfdecoupled TP molecule on Au(111) for the line trace shown in the inset STM image (2.7 V, 2 pA, 20  15 nm 2 ).

Figure 1 (

 1 Figure 1(c) shows the STM image of deposited molecular species on the Au(111) surface.Isolated bright spots, with a surface coverage of 4%, can be readily identified. Typical dimension parameters measured for these bright spots are 3.2±0.2 nm in width and 1.4± 0.1 nm in apparent height, respectively. Based on theoretical calculations with semiempirical quantum mechanical PM6 method,25 typical values for the configuration-optimized TP molecule are 3.0 nm in both width and height. The consistency between the measured width of the bright spots and the calculated width of a TP molecule suggests that each bright spot is likely to correspond to a single TP molecule, as subsequently further confirmed by PL and STML measurements. In addition, at low coverage, the TP molecules tend to remain as isolated single entities instead of forming molecular aggregates, suggesting that the TP molecules are anchored on the Au(111) surface after ESI deposition.

FIG. 2 .

 2 FIG. 2. (a) PL spectrum of self-decoupled TP molecules on Au(111) at 7 K, excited by a 532 nm pulse laser. (b) PL spectrum of TP molecules in a solution mixed with dichloromethane and methanol, measured at ambient temperature. (c) Time-resolved fluorescence decay curve of TP molecules adsorbed on Au(111) at 7 K. The blue curve is a fit convoluted with the instrument response function (IRF, black curve). (d) Time-resolved fluorescence decay curve of self-decoupled TP molecular solution measured at ambient temperature. A bandpass filter was used in the TCSPC lifetime measurements for the 00 emission.

  58 ns component, together with the observation of clear molecule-specific emission, indicates that such self-decoupled tetrapodal chemical design indeed suppresses effectively the fluorescence quenching from the metal substrate. All these PL results also give further justification that the bright individual spots observed in the STM image of Fig. 1(c) are indeed the target self-decoupled TP molecules.

FIG. 3 .

 3 FIG. 3. (a) Schematic of STML setup for a single self-decoupled TP molecule on Au(111). (b) One example of STML spectra (2.7 V, 30 pA, 30 s) at 78 K measured at the labeled positions in the inset STM image (2.7 V, 2 pA, 33  20 nm 2 ). (c) Another example of STML spectra (2.7 V, 30 pA, 30 s) at 78 K. (d) One type of typical STML spectra (2.7 V, 2 pA, 300 s) acquired at 7 K. Solid lines show the smoothed curves with raw data denoted by the gray dots. The red curves are STML spectra measured on the bare Au(111) surface near TP molecules.

  Both steady-state PL and transient lifetime measurements for isolated TP molecules on Au(111) justify the design concept through the observation of clear molecule-specific emissions and shortened lifetime components. Preliminary STML measurements on individual TP molecules on Au(111) indicate that tunneling electrons can readily modify the molecular configuration and even charging states, resulting in various luminescence spectral profiles and emission bands. The observation of significant emissions from the molecular center in STML measurements also suggests a standing-up orientation of the top-emitter transition dipole and further confirms the design concept. Future works should focus on the improvement of the stability of molecular configurations, so that molecular configurations, charging states, and luminescent properties can be better controlled and correlated. Our results provide new routes for controlling transition dipole orientations and for developing single-molecule light sources and organic optoelectronic devices.

This work is supported by the National Key R&D Program of China (grant numbers 2016YFA0200600 and 2017YFA0303500), the National Natural Science Foundation of China, the Chinese Academy of Sciences, and Anhui Initiative in Quantum Information Technologies. Talha Ijaz acknowledges support by the World Academy of Sciences. Yang Zhang acknowledges support by Excellent Young Scientist Foundation of the National Natural Science Foundation of China.