Tiny Grains Shining Bright in the Gaps of Herbig Ae Transitional Discs

TL;DR

This study uses advanced imaging and modeling to analyze the inner regions of two Herbig Ae transitional discs, revealing asymmetric features and challenging previous interpretations of companion candidates.

Contribution

Introduces a forward modeling approach combining radiative transfer and MCMC to analyze disc structures in high-resolution imaging data.

Findings

Detected emission interior to known disc walls at ~15AU and ~7AU.

Identified asymmetries in both discs.

Suggests companion candidates may be disc features rather than objects.

Abstract

This work presents a study of two Herbig Ae transitional discs, Oph IRS 48 and HD 169142; which both have reported rings in their dust density distributions. We use Keck-II/NIRC2 adaptive optics imaging observations in the L' filter (3.8 micron) to probe the regions of these discs inwards of ~20AU from the star. We introduce our method for investigating these transitional discs, which takes a forward modelling approach: making a model of the disc (using the Monte Carlo radiative transfer code RADMC), convolving it with point-spread functions of calibrator stars, and comparing the convolved models with the observational data. The disc surface density parameters are explored with a Monte Carlo Markov Chain technique. Our analysis recovers emission from both of the discs interior to the well known optically thick walls, modelled as a ring of emission at ~15AU in Oph IRS 48, and ~7AU for HD 169142, and identifies asymmetries in both discs. Given the brightness of the near-symmetric rings compared to the reported companion candidates, we suggest that the reported companion candidates can be interpreted as slightly asymmetric disc emission or illumination.

Tables

Table 1. Table 1 : Parameters of Oph IRS 48 and HD 169142

Oph IRS 48
Right Ascension	16h 27m 37.190s	1
Declination	$-{24}^{\circ }30\prime 35.03″$	1
Alternate Names	WLY 2-48,	1
	2MASS J16273718-2430350,	1
	YLW 46	1
Stellar Type	A0	2, 3
Distance to object	$134.4\pm 2.2\mathrm{pc}$	4
W1 magnitude	5.786	5
HD 169142
Right Ascension	18h 24m 29.779s	1
Declination	$-{29}^{\circ }46\prime 49.37″$	1
Alternate Names	MWC 925	1
Stellar Type	A5	6
Distance to object	$114.0\pm 0.8\mathrm{pc}$	4
W1 magnitude	6.203	5

Table 2. Table 3 : Tested radii and fluxes for rings in Oph IRS 48 and HD 169142, where F 𝐹 F refers to the fraction of flux that comes from each ring of radius r 𝑟 r . The visibility curves generated by these model systems are shown in Figure 5

Oph IRS 48
$F_1$	$r_1\left(\mathrm{AU}\right)$	$F_2$	$r_2\left(\mathrm{AU}\right)$	$F_3$	$r_3\left(\mathrm{AU}\right)$
0.71	-	-	-	-	-
0.66	14.3	-	-	-	-
0.38	12.7	0.29	25.2	-	-
0.22	9.8	0.26	17.0	0.20	26.9
HD 169142
$F_1$	$r_1\left(\mathrm{AU}\right)$	$F_2$	$r_2\left(\mathrm{AU}\right)$	$F_3$	$r_3\left(\mathrm{AU}\right)$
0.40	-	-	-	-	-
0.37	11.0	-	-	-	-
0.26	10.7	0.14	20.7	-	-
0.24	10.2	0.12	18.6	0.03	28.4

Table 3. Table 4 : Fixed parameters for Oph IRS 48 and HD 169142.

Oph IRS 48
Set Parameter	Value	Reference
Stellar Type	A0	1, 2
Star Temperature	$9000K$	2
Star Mass	$2.0M_{\odot }$	2
Star Radius	$2.24R_{\odot }$	Assumed
		(young A0 star)
Disc Mass	${10}^{-4}{M}_{\odot }$	3
$r_{\mathrm{wall}}$	$67\mathrm{AU}$	4, adjusted to new
		distance
${\delta }_{\mathrm{wall}}$	0.1	3
Dust Type	5.6 nm C & PAH	Assumed, 5
Distance to object	$134.4\mathrm{pc}$	6
Flux Ratio	7.8	This work
HD 169142
Set Parameter	Value	Reference
Stellar Type	A5	7
Star Temperature	$8250K$	7
Star Mass	$1.65M_{\odot }$	7
Star Radius	$1.56R_{\odot }$	Assumed, 7
Disc Mass	${10}^{-3}{M}_{\odot }$	Assumed, 7
$r_{\mathrm{wall}}$	$55\mathrm{AU}$	8, 9
${\delta }_{\mathrm{wall}}$	0.1	Assumed
Dust Type	5.6 nm C & PAH	Assumed, 5
Distance to object	$114\mathrm{pc}$	6
Flux Ratio	5.2	This work

Table 4. Table 5 : General parameters for the synthetic data set.

Parameters that will remain set.
Parameter	Value
Star Temperature	$10000K$
Star Mass	$2.0M_{\odot }$
Star Radius	$2.0R_{\odot }$
Disc Mass	${10}^{-3}{M}_{\odot }$
$r_{\mathrm{wall}}$	$50\mathrm{AU}$
${\delta }_{\mathrm{wall}}$	0.1
Dust Type	5.6 nm CP dust
Distance to object	$100\mathrm{pc}$

Table 5. Table 6 : Results for the synthetic companion data.

Synthetic Companion Data
Parameter	Target	Symmetric Model	Asymmetric Model	Companion Model
Dust to gas ratio	$1\times {10}^{-2}$	1.51$\pm 0.47\times {10}^{-2}$	6.56$\pm 1.80\times {10}^{-3}$	1.12$\pm 0.11\times {10}^{-2}$
${\delta }_d$	$1\times {10}^{-4}$	7.87$\pm 1.61\times {10}^{-5}$	1.27$\pm 0.20\times {10}^{-4}$	9.83$\pm 0.60\times {10}^{-5}$
${\delta }_1$	$1\times {10}^{-2}$	8.89$\pm 1.42\times {10}^{-3}$	1.22$\pm 0.14\times {10}^{-2}$	1.11$\pm 0.08\times {10}^{-2}$
$r_d$ $\left(\mathrm{AU}\right)$	0.5	0.39$\pm 0.07$	0.64$\pm 0.08$	0.42$\pm 0.03$
$r_d$ $\left(\mathrm{mas}\right)$	5	3.9$\pm 0.7$	6.4$\pm 0.8$	4.2$\pm 0.3$
$r_1$ $\left(\mathrm{AU}\right)$	10.0	10.0$\pm 0.1$	10.2$\pm 0.1$	9.95$\pm 0.10$
$r_1$ $\left(\mathrm{mas}\right)$	100	100$\pm 1$	102$\pm 1$	99.5$\pm 1.0$
$r_2$ $\left(\mathrm{AU}\right)$	25.0	23.9$\pm 0.6$	23.3$\pm 0.4$	23.6$\pm 0.5$
$r_2$ $\left(\mathrm{mas}\right)$	250	239$\pm 6$	233$\pm 4$	236$\pm 5$
Inclination ${(}^{\circ })$	30.0	29.7$\pm 1.2$	31.9$\pm 1.1$	30.7$\pm 1.2$
Position angle ${(}^{\circ })$	50.0	54.1$\pm 2.9$	57.4$\pm 2.2$	56.4$\pm 2.0$
Star offset $\left(\mathrm{mas}\right)$	1.3	-	1.6$\pm 0.3$	1.1$\pm 0.1$
Star position angle ${(}^{\circ })$	90.9	-	79$\pm 15$	75$\pm 13$
Companion offset $\left(\mathrm{mas}\right)$	180	-	-	186$\pm 2$
Companion position angle ${(}^{\circ })$	232.8	-	-	235.6$\pm 2.2$
Companion contrast (mag)	5.4	-	-	5.7$\pm 0.2$
Flux Ratio	30.31	30.33	30.32	30.32
reduced ${\chi }_{\mathrm{shot}}^2$	-	67.5	27.3	24.4

Table 6. Table 7 : Best fit model parameters for 2015 Oph IRS 48 data and 2014 HD 169142 data, using mean and standard deviation. Best-fit models were not calculated for the 2016 data of both objects, however the best-fit models from the other years were tested on the 2016 data.

Oph IRS 48
Parameter	Symmetric	Asymmetric	Companion
Dust to gas ratio	3.10$\pm 0.66\times {10}^{-3}$	3.20$\pm 0.26\times {10}^{-3}$	2.06$\pm 0.09\times {10}^{-3}$
${\delta }_d$	2.48$\pm 0.43\times {10}^{-4}$	2.11$\pm 0.23\times {10}^{-4}$	4.20$\pm 0.35\times {10}^{-4}$
${\delta }_1$	1.01$\pm 0.10\times {10}^{-2}$	9.40$\pm 0.61\times {10}^{-3}$	1.24$\pm 0.06\times {10}^{-2}$
$r_d$ $\left(\mathrm{AU}\right)$	1.6$\pm 0.2$	1.4$\pm 0.1$	2.0$\pm 0.1$
$r_d$ $\left(\mathrm{mas}\right)$	$12\pm 1$	$10\pm 1$	$15\pm 1$
$r_1$ $\left(\mathrm{AU}\right)$	15.0$\pm 0.3$	14.6$\pm 0.1$	15.2$\pm 0.2$
$r_1$ $\left(\mathrm{mas}\right)$	$112\pm 2$	$108\pm 1$	$113\pm 1$
$r_2$ $\left(\mathrm{AU}\right)$	32.1$\pm 0.9$	31.0$\pm 0.3$	31.3$\pm 0.3$
$r_2$ $\left(\mathrm{mas}\right)$	$238\pm 7$	$231\pm 2$	$233\pm 2$
Inclination ${(}^{\circ })$	54.8$\pm 0.5$	54.2$\pm 0.3$	52.4$\pm 0.3$
Position angle ${(}^{\circ })$	99.0$\pm 0.6$	98.8$\pm 0.3$	96.5$\pm 0.4$
Star offset ($\mathrm{mas}$)	-	5.0$\pm 0.4$	5.3$\pm 0.3$
Star offset position angle ${(}^{\circ })$	-	118.7$\pm 7.8$	122.9$\pm 8.2$
Companion offset ($\mathrm{mas}$)	-	-	104$\pm 2$
Companion position angle ${(}^{\circ })$	-	-	288.6$\pm 4.6$
Companion contrast (mag)	-	-	3.97$\pm 0.05$
Flux Ratio (7.8)	7.8	7.8	7.8
Reduced ${\chi }_{\mathrm{shot}}^2$	65.1	60.0	48.0
Reduced ${\chi }_{\mathrm{shot}}^2$ 2016 data	68.3	61.0	52.6
HD 169142
Parameter	Symmetric	Asymmetric	Companion∗
Dust to gas ratio	8.88$\pm 1.34\times {10}^{-4}$	7.16$\pm 0.63\times {10}^{-4}$	6.81$\pm 0.64\times {10}^{-4}$
${\delta }_d$	6.14$\pm 0.57\times {10}^{-4}$	7.14$\pm 0.47\times {10}^{-4}$	7.17$\pm 0.42\times {10}^{-4}$
${\delta }_1$	2.77$\pm 0.38\times {10}^{-2}$	2.02$\pm 0.41\times {10}^{-2}$	2.02$\pm 0.38\times {10}^{-2}$
$r_d$ $\left(\mathrm{AU}\right)$	0.11$\pm 0.01$	0.11$\pm 0.01$	0.11$\pm 0.01$
$r_d$ $\left(\mathrm{mas}\right)$	$0.96\pm 0.08$	$0.96\pm 0.08$	$0.96\pm 0.08$
$r_1$ $\left(\mathrm{AU}\right)$	8.4$\pm 0.4$	7.3$\pm 0.5$	7.3$\pm 0.4$
$r_1$ $\left(\mathrm{mas}\right)$	$74\pm 4$	$64\pm 4$	$64\pm 4$
$r_2$ $\left(\mathrm{AU}\right)$	15.4$\pm 1.3$	13.5$\pm 0.6$	13.4$\pm 0.6$
$r_2$ $\left(\mathrm{mas}\right)$	$135\pm 11$	$118\pm 5$	$118\pm 5$
Inclination ${(}^{\circ })$	30.3$\pm 1.0$	29.8$\pm 0.4$	29.5$\pm 0.6$
Position angle ${(}^{\circ })$	12.6$\pm 2.1$	15.0$\pm 0.9$	15.7$\pm 1.5$
Star offset ($\mathrm{mas}$)	-	0.09$\pm 0.03$	0.11$\pm 0.03$
Star offset position angle ${(}^{\circ })$	-	15.0$\pm 9.1$	7.4$\pm 8.9$
Companion offset ($\mathrm{mas}$)	-	-	131 (fixed)
Companion position angle ${(}^{\circ })$	-	-	2.6 (fixed)
Companion contrast (mag)	-	-	5.7$\pm 0.7$
Flux Ratio (5.2)	5.2	5.2	5.2
Reduced ${\chi }_{\mathrm{shot}}^2$	5.91	5.58	5.51
Reduced ${\chi }_{\mathrm{shot}}^2$ 2016 data	44.1	42.8	45.7

Table 7. Table 8 : Values used for the spatial grid in RADMC3D .

Variable	Grid Limits	Number of points between Grid Limits
$r$ ($\mathrm{AU}$)	$r_{\mathrm{dust}}$, $r_1$, $r_1+0.1$, $r_2$, $r_2\times 1.1$, 100	5, 20, 30, 20, 40
$\theta$	0, $\pi /3.$, $\pi /2.$, $2\pi /3.$, $\pi$	10, 30, 30, 10
$\varphi$	0, 2$\pi$	60

Table 8. Table 9 : Set parameters for all models, most are defaults that RADMC3D  chooses, from the ‘problem_params.inp’ file.

Parameter	Value
Continuous Stellar Source	False
Discrete Stellar Source	True
Coordinate System	Spherical
Number of Points for Wavelength grid	19, 50, 30
Number of Points for Wavelength grid (SED)	100, 100, 30
Bounds of wavelength grid (µm)	0.1, 7.0, 25.0, 10000
Bounds of wavelength grid (µm) (SED)	0.1, 1.5, 25.0, 10000
Number of refinement levels	3
Number of the original grid cells to refine	3
Number of grid cells to create in a refinement level	3
Bulk Density of materials (${\mathrm{gcm}}^{-3}$)	3.6, 1.8
Grain size distribution power exponent	-3.5
Maximum grain size	10.0
Minimum grain size	0.1
Mass fractions of the dust components to be mixed	0.75, 0.25
Number of grain sizes	1
Use finite size of star	no - take as point source
Modified Random Walk	Off
Number of photons for image generation	$2\times {10}^4$
Number of photons for SED generation	$1\times {10}^5$
Output format for RADMC3D files	ASCII
Scattering mode	isotropic
Dust Temperature equal Gas Temperature	yes
Background Density (${\mathrm{gcm}}^{-3}$)	$1\times {10}^{-30}$
Pressure scale height at innermost radius ($\mathrm{AU}$)	0.0
Reference radius at which $H_p$/R is taken ($\mathrm{AU}$)	100.
Ratio of the pressure scale height over radius at reference radius for $H_p$/R	0.1
Flaring index	1./7.
Power exponent of the surface density distribution as a function of radius	-1.0
Outer boundary of the puffed-up inner rim in terms of innermost radius	0.0
Outer radius of the disc $\mathrm{AU}$	100.
Surface density at outer radius of the disc	0.0
Surface density type	polynomial
Outer boundary of the smoothing in the inner rim in terms of innermost radius	1.0
Power exponent of the density reduction inside of the the inner rim smoothing	0.0

Table 9. Table 10 : References for SED data for Oph IRS 48 and HD 169142.

Oph IRS 48
Wavelength ($\mathrm{\mu m}$)	Reference
0.43, 0.64	Zacharias et al. (2005)
0.64, 0.79	Erickson et al. (2011)
1.24, 1.66, 2.16	Cutri et al. (2003)
3.4, 4.6, 12, 22	Wright et al. (2010)
3.6, 4.5, 5.8, 8, 70	van Kempen et al. (2009)
12, 25, 60 ,100	Helou & Walker (1988)
18.7	Yamamura et al. (2010)
18.7	Geers et al. (2007a)
70	Fedele et al. (2013)
HD 169142
Wavelength (µm)	Reference
0.15, 0.18, 0.22, 0.25, 0.33	IUE archival data
0.36, 0.44, 0.55, 0.64,	Sylvester et al. (1996)
0.79, 3.77, 4.78
1.24, 1.65, 2.16	2MASS All-sky Point
	Source Catalog
3.35, 4.6, 11.6, 22.1	WISE All-sky Data Release
	Catalog
10.8, 18,2	Jayawardhana et al. (2001)
11.7, 18.3	Mariñas et al. (2011)
12, 25, 60, 100	IRAS Point Source Catalog
18	AKARI/IRC All-sky Survey Point
	Source Catalog
	(Ishihara et al., 2010)
18.8, 24.5	Honda et al. (2012)
65, 90, 140, 160	AKARI/FIS All-sky Survey
	Point Source Catalog
	(Version 1.0)
70, 160	Meeus et al. (2010)

Equations

$$r_{uv}^{\prime }$$

$${\chi }_{\mathrm{shot},\mathrm{i}}^2={\min }_j\sum _k\frac{(D_{i,k}-(M\otimes P_j{)}_k{)}^2}{{\sigma }_k^2}.$$

$${\chi }_{\mathrm{tot}}^2=\sum _i{\chi }_{\mathrm{shot},i}^2.$$

$$ln(\mathrm{likelihood})=\frac{-{\chi }_{\mathrm{tot}}^2}{2T_{\mathrm{MC}}},$$

$${\chi }^2=\frac{-ln(\mathrm{likelihood})\times 2\times T_{\mathrm{MC}}}{n_{\mathrm{pixels}}\times n_{\mathrm{images}}},$$

Full text

††pubyear: 2019††pagerange: Tiny Grains Shining Bright in the Gaps of Herbig Ae Transitional Discs–Tiny Grains Shining Bright in the Gaps of Herbig Ae Transitional Discs

1 Introduction

Transitional discs are a subset of protoplanetary discs that have a region which is depleted in dust (but not necessarily depleted in gas), known as a gap or hole. These objects are a subset of Lada (1987) Class II objects, with flat or declining mid-infrared (mid-IR) excesses in their spectral energy distributions (SEDs). Giant planets that have formed in protoplanetary discs are predicted to carve out gaps or holes, producing structures that might be observed in transitional discs (Strom et al., 1989; Skrutskie et al., 1990; Marsh & Mahoney, 1992, 1993). It has also been suggested that the gaps in protoplanetary discs could be caused by multi-planet systems (Dodson-Robinson & Salyk, 2011), however the prevalence of systems with 3 or more Jupiter mass planets is very low ($<1$%, Han et al., 2014), so this is unlikely to be the only cause of transitional disc structure. By direct imaging of transitional discs, we can look for signs of planet formation or other disc evolution such as dust asymmetries, dust depleted regions, gaps, rings, and holes.

The transitional discs studied in this work are Oph IRS 48 and HD 169142. Both are Herbig Ae stars known to have polycyclic aromatic hydrocarbons (PAHs) in their discs. Information on these objects is summarised in Table 1. Both objects are known to have structure in their outer discs with rings of emission at $50-60\,\mathrm{AU}$ identified from observations at mm wavelengths (e.g. Geers et al., 2007a; Brown et al., 2012b; van der Marel et al., 2013; Quanz et al., 2013; Maaskant et al., 2014).

There are several objects for which there are detections of dust within the regions thought to be depleted of CO. This is the case for V1247 Orionis, for which the presence of carbon-rich dust inside the gap region has been found through near-IR observations (Kraus et al., 2013). Previous modelling of infrared spectra of both Oph IRS 48 and HD 169142 have indicated the likely presence of PAHs within the gas gap regions (regions depleted in CO; Geers et al., 2007b; Maaskant et al., 2014; Seok & Li, 2017).

Oph IRS 48 is located in the $\rho$ Ophiuchus star forming region, as catalogued by Elias (1978) and confirmed by Wilking et al. (1989). The average distance to the cloud’s core was $120.0^{+4.5}_{-4.2}\,\mathrm{pc}$ found by Loinard et al. (2008) using Very Long Baseline Array (VLBA) data. The parallax for Oph IRS 48 as given by Gaia DR2 (Gaia Collaboration et al., 2018) is $7.44\pm 0.12$ millisecond of arc (mas), corresponding to a distance of $\sim 134.4\pm 2.2\,\mathrm{pc}$, and this is the distance adopted in this work.

The outer disc of Oph IRS 48 is known to have a strong asymmetric feature in the millimetre-sized grains. The asymmetry in the dust of the outer disc (at a deprojected distance of $67\,\mathrm{AU}$ from the star) was discovered with ALMA sub-millimetre observations and indicates a separation of the micron- and millimetre-sized dust grains (van der Marel et al., 2013). There is also a separation of the very small grains that do not follow the larger micron sized dust, as reported as unresolved emission in Geers et al. (2007a). Later observations found that centimetre-sized dust grains are further concentrated in the region of the millimetre-grains (van der Marel et al., 2015), consistent with being caused by a vortex induced by a pressure or density gradient generating a dust trap, which may have been generated by a companion (van der Marel et al., 2013, 2015).

Inside the $67\,\mathrm{AU}$ asymmetric ring, the structure of the Oph IRS 48 disc becomes more complicated. It is thought that there is a depletion of dust within $67\,\mathrm{AU}$, given that there is a no excess at wavelengths short of 10 µm in the SED (Maaskant et al., 2013). However, where these depletions are seen can depend on the wavelength of the observations, as different materials (gas, dust, PAHs, large or small grains) are traced by different wavelengths (e.g. Brown et al., 2012b). The disc is considered to be mostly depleted of dust within $\sim 23\,\mathrm{AU}$ (Bruderer et al., 2014), where there is a wall, but the flux deficit caused by depletion depends strongly on the wavelength of observation. Oph IRS 48 is still accreting (accretion rate $\sim 10^{-9}\,\mathrm{M}_{\odot}\,\mathrm{yr^{-1}}$, Salyk et al., 2013), and SED modelling also confirms the presence of inner-disc material (e.g. Bruderer et al., 2014).

HD 169142 was previously thought from optical photometry to be located at a distance of $145\,\mathrm{pc}$ (Sylvester et al., 1996), but the parallax was measured with Gaia (Gaia Collaboration et al., 2018) to be $8.77\pm 0.06\,$mas, corresponding to a distance of $114.0\pm 0.8\,\mathrm{pc}$. We adopt the new distance in this work, and have adjusted the linear radii of the previously noted features to correspond to the new distance.

HD 169142 was reported as having circumstellar material by Walker & Wolstencroft (1988), and subsequent studies have found that this disc has structure, including an inner cavity and rings (e.g. Panić et al., 2008; Honda et al., 2012; Quanz et al., 2013; Osorio et al., 2014; Seok & Li, 2016; Monnier et al., 2017; Fedele et al., 2017; Macías et al., 2017). The most commonly modelled structure for HD 169142 thus far includes a small inner rim (Panić et al., 2008; Honda et al., 2012), with an optically thick ring at $\sim 20\,\mathrm{AU}$ (Osorio et al., 2014; Honda et al., 2012; Panić et al., 2008), and a gap from $31\,\mathrm{AU}$ to an optically thick wall at $55\,\mathrm{AU}$ (Quanz et al., 2013).

It is possible that the gaps in the disc are caused by planet-disc interactions. Biller et al. (2014) and Reggiani et al. (2014) suggested from L’ (3.8 µm) observations using a coronagraph that a potential companion is located at a distance of either $\sim 13\,\mathrm{AU}$ ($110\pm{30}\,$mas, position angle of $0\pm{14}\,^{\circ}$, Biller et al., 2014) or $\sim 18\,\mathrm{AU}$ (156$\pm{32}\,$mas, position angle of $7.4\pm{11.3}\,^{\circ}$, Reggiani et al., 2014) from the star (interior to the $\sim 20\,\mathrm{AU}$ ring). The potential companion has an L’ apparent Vega magnitude of $12.2\pm 0.5$ mag (Reggiani et al., 2014), which corresponds to a contrast of 6.4 mag (Biller et al., 2014), or $6.5\pm 0.5$ mag (Reggiani et al., 2014). These studies both used observations taken with the Very Large Telescope using NACO, and the coronograph available for this instrument. Both studies used similar methods to detect the candidate companion, which included angular differential imaging (ADI) that self-subtracted any azimuthally symmetric structures. Other studies investigate whether this possible companion, and potentially a second candidate companion further out in the disc are causing the structural features seen in the disc (Fedele et al., 2017; Kanagawa et al., 2015).

There have been few previous imaging studies of these discs at wavelengths shorter than $\sim$8 $\mu$m – wavelengths most sensitive to structures inside the known rings and where planet formation signatures, including those of circumplanetary accretion discs (Zhu, 2015), could be seen. The key reason for this gap is the higher contrasts and angular resolutions needed for shorter wavelengths, meaning that specialised techniques are required, such as ADI (Reggiani et al., 2014), or aperture mask interferometry (Schworer et al., 2017).

It is possible that there is a companion in each of the discs, clearing the gaps and driving the asymmetries (e.g. van der Marel et al., 2013; Biller et al., 2014; Reggiani et al., 2014). To investigate whether there is a planet inside the gas gaps of these discs at $\sim 20\,\mathrm{AU}$ we present our observations and analysis of both discs as follows. Section 2 discusses the observations used for this work, which differs from previous studies as we are more sensitive to extended structures than aperture mask interferometry, but are still sensitive to circularly symmetric features, unlike ADI. Section 3 discusses evidence for significant emission inside the previously inferred walls. Sections 4 and 5 outlines the computational methods used in this paper and tests them on a synthetic data set. Section 6 describes the resulting best fit physical disc models. Section 7 summarises our results in the context of the field, and discusses possible future work in this area.

2 Observations and Data Reduction

2.1 Observations

Observations were acquired with the NIRC2 instrument on Keck II over three observing runs (June 2014, 2015, 2016). All observations used in this work were taken with the L’ filter of the NIRC2 camera, using a ‘large hex pupil’, a $512\times 512$ subarray, and a two-point dither mode separated by $\sim$3.5”, utilising the top left and bottom right quadrants of the detector. These observations did not use aperture masking or coronography. Table 2 summarises these observations, including total exposure times and point spread function (PSF) calibrators. We use both natural and laser guide star adaptive optics (AO) in this work (because HD 169142 is bright enough to not need the laser, and Oph IRS 48 is fainter). Our 2014 and 2016 observations of HD 169142 used natural guide star AO and the 2015 and 2016 observations of Oph IRS 48 used laser guide star AO.

The calibrators for HD 169142 were chosen specifically for that object, while Oph IRS 48 was calibrated against other Class II objects that formed a survey of the Ophiuchus star forming region for accreting exoplanets. These calibrators were expected to have bright, unresolved inner discs, with cross-calibration between different calibrators used to check for any measurably resolved objects. No potential calibrators were eliminated.

2.2 Data Reduction

Each image was corrected for detector nonlinearity using the algorithm from the IDL program linearize_nirc2.pro 111http://www.astro.sunysb.edu/metchev/ao.html, and then divided by a mean dome flat field. Images taken with the star at the opposite dither position were used as an estimate of sky, which was subtracted from each image. Bad pixels were primarily identified by searching for outliers in mean pixel value or variance in a sequence of dark frames or flat frames. Additional bad pixels or cosmic rays were found by Fourier transforming each data frame, masking out spatial frequencies below Nyquist, inverse transforming and searching for significant peaks. The bad pixels were corrected using the algorithm in Ireland (2013), which involves finding the pixel value that minimises Fourier power above the Nyquist frequency. Images were also centred to the peak value (which is typically the centre of the star) and cropped to $128\times 128$ pixels, to remove edge effects.

When making the set of cleaned images to be used for subsequent analysis, images extincted by cloud by more than a factor of 2 or with Strehl more than $\sim$20% below the maximum were rejected. For the 2014 HD 169142 data, only one of the 12 target images was rejected. For the 2015 Oph IRS 48 data, three of the 18 target images were rejected. For the 2016 HD 169142 data, two of the 10 target images were rejected. For the 2016 Oph IRS 48 data, six of the 18 target images were rejected. The 2016 data were taken through cloud, and for all data sets we began data acquisition before the low-bandwidth wavefront sensor showed low errors. All raw data are publicly available in the Keck observatory archive, and we have made the cleaned data cubes222http://www.mso.anu.edu.au/$\sim$mireland/Birchall18/ available.

3 Resolved Disc Emission

3.1 Richardson-Lucy Deconvolution

The Richardson-Lucy method of deconvolution is an iterative deconvolution method based on finding the mostly likely image that fits a data set given a known PSF. It was developed independently by Richardson (1972) and Lucy (1974). We performed a Richardson-Lucy deconvolution on our reduced data, in order to isolate the resolved structure in the discs, and as a simple first-pass imaging algorithm to inspire physical model-fitting.

We used 50 iterations of the Richardson-Lucy algorithm, and our initial model was a point-source. We used this algorithm on all target-image – calibrator-image pairs. However, we only considered the calibrator image that, when convolved, had the smallest root-mean-square (RMS) with respect to the target star image (i.e. the one that returned the image most similar to each target image). We chose to use 50 iterations because at this point the RMS between the target and calibrator images is stable, and not decreasing dramatically with more iterations. Also, with more iterations, the deconvolution noticeably over fits the data.

One should note that the point-source initial model for our implementation of the Richardson-Lucy algorithm has an effect of regularising the final model to be a point source with resolved structure, rather than a marginally-resolved central source with structure. A point source with a disc is a suitable physical approximation for our deconvolved images as the star and its inner ring of material ($<1\,\mathrm{AU}$) is unresolved.

The deconvolved images for Oph IRS 48 and HD 169142 are shown in the right panels of Figures 1 and 2. In both cases, a resolved ring of emission is clearly visible at radii within the previously published transitional disc holes ($\sim 23\,\mathrm{AU}$ for Oph IRS 48, and $\sim 20\,\mathrm{AU}$ for HD 169142).

3.2 Fourier Analysis

The presence of significant rings interior to the $\sim\,20\,\mathrm{AU}$ wall inferred from previous data sets (e.g. Bruderer et al., 2014; Osorio et al., 2014, for Oph IRS 48 and HD 169142 respectively) in the Richardson-Lucy deconvolved images is something that could in principle also arise from instability in the adaptive optics system. This possibility was ameliorated in our data collection by having multiple epochs on the targets and multiple calibrators at each epoch. In addition, Strehl ratios were reasonably high (0.5–0.7) and stable (within $\sim$10% in the used images) given the long wavelength of our observation.

In order to quantitatively examine the effects of different calibrators, we performed a power spectrum analysis, as pioneered in speckle interferometry (Labeyrie, 1970), where spatial power spectra were calculated for both individual calibrator stars and target stars, and then divided by the mean of the PSF calibrator power spectra in two dimensions. These calibrated power spectra are then equal to the power spectrum of the target itself, with the dispersion amongst the calibrator power spectra representing the uncertainty in this technique.

We plot the azimuthal averages of the visibilities (i.e. the square root of the power spectra) from each non-rejected frame of the target and calibrator in Figures 3 and 4 for the 2014 and 2015 epochs, divided by the mean azimuthally averaged visibilities. It can be seen that for both HD 169142 and Oph IRS 48, the target visibility curves fall well outside the range of the calibrators, indicative of dominant structures that are fully resolved on 0.2–0.3" scales. The asymptotic visibilities of $\sim$0.26 for Oph IRS 48 and $\sim$0.57 for HD 169142 measures directly that only $\sim$26% and $\sim$57% of the L’ flux from Oph IRS 48 and HD 169142 come from the sum of the central star and inner disc, in disagreement with previous models constrained by the SED (e.g. Seok & Li, 2016; Bruderer et al., 2014, for HD 169142 and Oph IRS 48 respectively).

In order to further analyse these visibilities, we used a model where each target was represented by an azimuthally symmetric inclined thin disc. The de-projected spatial frequency co-ordinate $r^{\prime}_{uv}$ is derived from the $(u,v)$ coordinate using:

[TABLE]

where $\theta_{uv}$ is the position angle of the $(u,v)$ coordinate, $r_{uv}$ its magnitude, $\theta_{\rm maj}$ the position angle of the disc major axis and $i$ the disc inclination. The Fourier power was then binned and averaged over this de-projected Fourier coordinate $r^{\prime}_{uv}$ to enable robust data visualisation and azimuthally symmetric fitting. For Oph IRS 48, an inclination of $55^{\circ}$ and a major axis position angle of $99^{\circ}$ were used, while for HD 169142, an inclination of $30^{\circ}$ and position angle of $13^{\circ}$ were used (based on literature values and our findings in Section 6 - where they differ from literature values). After azimuthally averaging the deprojected power spectra and taking the square root, we plot the visibility curves in Figure 5. Rising visibilities to mid or high spatial frequencies demonstrate that fits of ring-like or sharp-edged flux distributions are appropriate, rather than e.g. Gaussian or power-law dust distributions. We fitted several parameterized models to these distributions in order to determine the required model complexity (Table 3 and Figure 5). One ring at the previous location of the wall inferred from the SED or long wavelength observations is clearly inconsistent with the data, which requires a more continuous distribution of emission, including an emission component inside the previous ($\sim 20\,\mathrm{AU}$) gap.

4 Theoretical disc modelling

In order to interpret the observational data, radiative transfer models of the discs were generated. We chose to follow previous literature on these discs, (e.g. Bruderer et al., 2014, Oph IRS 48) and to fit power-law surface density distributions with multiple gaps. We recognise that alternative models with smoother distributions of dust and gas may fit our data equally well, and could be attempted as an extension to this work, constrained by hydrodynamical models or more detailed SED modelling, for example. Limitations of our work are discussed in Section 4.2.5. The radiative transfer disc models used in this work are made using RADMC3D (Dullemond, 2012)333Code is available at: http://www.ita.uni-heidelberg.de/$\sim$dullemond/software/radmc-3d/, which is a Monte Carlo radiative transfer code. A forward modelling approach was taken; the models are generated by this code, rotated to match the observed position angle at each epoch, and convolved with a PSF so that a comparison can be made between the model and the observational data.

The radiative transfer models are discussed in Section 4.1, and the details of the parameters of the models are in Section 4.2.

4.1 Radiative transfer

The protoplanetary disc function of RADMC3DPy, the companion Python library to RADMC3D (Dullemond, 2012), was used to create three dimensional disc density distributions with structures such as gaps and rings, which were then fed into RADMC3D.

Three model types are tested; a symmetric disc, an asymmetric disc and an asymmetric disc with a companion, with each model type being an extension on the previous one. These three model types are used to identify the presence of asymmetries in the data and attempt to characterise the asymmetry to see if it is point-like. To simulate an asymmetric disc, the central star can be moved slightly which gives an asymmetric dust illumination. In some cases companions can cause the disc to become eccentric (Thalmann et al., 2010), meaning that the star is no longer at the centre of the disc. Sub-stellar companions are simulated by adding photospheric emission from a second point source to the RADMC3D disc model. The companion source is less massive, smaller in radius and cooler than the central star. The parameters involved in generating these models are discussed in Section 4.2.

4.2 Parameters

4.2.1 Globally constant parameters

The stellar emission type, dust model, and stellar source type were held constant for both objects. We use Kurucz444http://kurucz.harvard.edu/grids.html (Castelli & Kurucz, 2004) stellar models, with solar abundances and micro-turbulence of $2\,\mathrm{km\,s^{-1}}$ (an arbitrary choice, as we are not attempting to model the spectrum of the star in detail). The RADMC3D emission source used for either object is a point source.

The chosen dust model used for the final modelling of both objects is a fine-grained carbon and PAH dust mixture (henceforth CP dust). The dust model chosen was the Draine & Li (2001, 2007) dust, which includes carbon in graphitic form and also grain sizes down to nm scale which are mostly PAHs. The discs of both Oph IRS 48 and HD 169142 are known to contain PAHs (e.g. Geers et al., 2007b; Maaskant et al., 2014; Seok & Li, 2016, Oph IRS 48, both, HD 169142, respectively). The reasoning behind the choice of dust is explained further in Section 6.1. The dust best suited from the Draine & Li (2001, 2007) set of neutral dust was a grain-size of $5.6\times 10^{-3}$ µm (5.6 nm). At 5.6 nm the dust mixture can give strong PAH features, but also has a contribution from a more general graphitic source (Draine & Li, 2001; Draine, 2003; Draine & Li, 2007). This 5.6 nm CP dust is used for the determination of the other disc parameters.

Other parameters that are the same for both objects include the number of photons (set to 106), and the spatial resolution. These were both set based on results of convergence tests, ensuring a reproducible output. The interested reader is directed to Appendix C for information on some of the other RADMC3D parameters that were used and convergence studies.

4.2.2 Object constant parameters

The stellar temperature and mass are fixed for each object based on literature spectral type and age for each object, with the radius coarsely adjusted to match the observed visible flux (Table 4). As we model the small-grained dust only, we can not fit directly to the gas density. We choose to fix the mass of the disc (most of the mass is in the outer disc), to a value from the literature. Instead of a varying disc mass, we have a dust-to-gas mass ratio that is independent of radius, with density (of dust and gas) that is able to deviate from dust-to-gas ratio at various radii, as discussed in the next section. Both of the objects studied in this work are known to have either outer rings or asymmetries observed in much longer wavelengths (e.g. mm wavelengths), so we include these known parameters as the outer wall of the disc, $r_{\mathrm{wall}}$ and outer disc dust depletion $\delta_{\mathrm{wall}}$. The flux ratio refers to the ratio of flux between the star and the disc, and is discussed in more detail in Sections 4.2.4 and 6.1.

4.2.3 Parameters varied for each object

The symmetric model uses the following eight parameters: the overall dust-to-gas ratio; the inner radius of the disc $r_{d}$; the depletion factor of the inner disc, $\delta_{d}$; the radius of the first wall, $r_{1}$; the depletion of the first gap, $\delta_{1}$; the radius of the second wall, $r_{2}$; the inclination of the disc; and the position angle of the object on the sky. The first six parameters are illustrated in a schematic plot of the midplane density of a disc in Figure 6, and a schematic indicating the radii and other parameters is shown in Figure 7. The inclination is whether the object is seen face-on, edge-on or somewhere in-between. The position angle quoted is the position angle of the disc major axis from North towards East, with the near side of the modelled disc at this position angle plus 90 degrees. Note, however, that it is uncertain which minor axis corresponds to the near side.

The use of two walls in the disc at $r_{1}$ and $r_{2}$ allows us to explore the position of the known wall at $r_{2}$, while introducing a new wall at a closer radius, $r_{1}$. Both the objects investigated here are known to have rings of emission at a radius of $\sim 20\,\mathrm{AU}$, and so to recover the position of that ring, while also fitting our data we use two radii here.

The asymmetric model case moves the central star. Moving the central star simulates the appearance of an asymmetric disc, without making an asymmetric dust distribution, which could add an arbitrary degree of complexity. For YSOs, it is possible for the star to be slightly offset from the centre of the disc, for example this was one explanation for the peri-centre offset of LkCa15 (in e.g., Thalmann et al., 2010). The eccentricity of the disc in systems where the star has a peri-centre offset is expected to be driven by the presence of a companion. The off-centre star placement means that one side of the disc is more strongly illuminated than the other, and is indicated by the position of the star in Figure 7. The asymmetry parameters come from changing the $x$ and $y$ positions of the central star (units of this movement are in $\mathrm{AU}$). The positive $x$ direction is to the right along the semi-major axis, and the positive $y$ direction is to the right along the semi-minor axis. These parameters are converted to a separation in mas and position angle offset when presented in Section 6.

The companion model is constructed by adding a second point source to the asymmetric disc model. We do not explicitly suggest that a point-like asymmetry is likely to be an accreting exoplanet at our poor angular resolution, but include a point source for simplicity of parametric modelling. The parameters that are varied in the model for the companion to ensure a better fit are the radius and $x$ and $y$ positions of the companion (where the $x$ and $y$ directions are the same as for the asymmetry). As with the star parameters, the position of the companion is converted to a separation given in mas and a position angle, so that it is more easily compared to literature values for the position of potential companions. The radius of the planet was used as a proxy for the brightness within our simulation setup and it is converted to a contrast of the companion to the total model, given in magnitudes. Please see Appendix B for details on the companion parameters. An example of a companion is denoted by the red circle in Figure 7.

4.2.4 Priors

For some of the parameters that will be explored, Bayesian priors were set. The priors for the radii make sure that they stay in the order of $r_{d}<r_{1}<r_{2}<r_{\rm wall}$. As minimum inner disc radius, we set 0.1 AU, however, this prior is not important in the determination of the likelihood peak.

When the star is placed off-centre, there is a prior on the displacement so that the displacement does not exceed a radius of $1\,\mathrm{AU}$ from its central position, which we find does not affect the likelihood peak.

There is also a prior on the flux ratio of the star to the total flux (star+disc). The flux ratio prior is used so that without actually modelling the full SED (which would require multiple dust models), we can still constrain the amount of dust in the unresolved inner disc, and to have a self-consistent disc model that fully takes into account shadowing of the outer disc by the inner disc.

For both Oph IRS 48 and HD 169142, we fit the flux ratio of the entire disc and star system to the star. We calculate the flux ratio to fit by taking the stellar photosphere models of the objects from our simulations and literature observations at 3.4 (WISE) or 3.6µm (Spitzer/IRAC) photometry. To find our model flux ratio, we calculate the total intensity of a model with just a star, and for a model of a star and a disc, and then calculate the ratio of these intensities.

For Oph IRS 48 the flux ratio used was $\sim 7.8$ and for HD 169142 it was $\sim 5.2$. Uncertainties for each of the flux ratios are estimated for use in the emcee part of the method. The flux ratio places a constraint on the brightness of the disc in the model, and in particular, the unresolved inner disc.

4.2.5 Limitations

There are some limitations to our method and the way our models are constructed. Our models have a similar power-law density structure to those used in other works (e.g. Bruderer et al., 2014; Maaskant et al., 2014, Oph IRS 48, both, respectively), as seen in Figure 6. Alternative density structures (e.g. without discontinuities) would also fit the data, however we choose to analyse one class of models with sufficient parameters. We consider a phenomenological model, and do not test whether the vertical structure is in equilibrium. We also assume that the distances to the objects are fixed, and do not use the uncertainties associated with them in our modelling. The interested reader is directed to Appendix C for more details on the parameters not discussed in Section 4.2.

4.3 Convolving the models

The final step before the model image can be compared to the observational data is to convolve it with a PSF, chosen from our PSF library generated by the observation of calibrator objects. The model image is rotated to the same position angle as the observational data, then convolved with each of the PSFs in the relevant PSF library, shifted on a sub-pixel scale to match the position of the target data, and then the best PSF for convolution is chosen for each model image. The best PSF is found for each of the rotated images by comparing with observational data and calculating the $\chi^{2}_{\rm shot}$ for each of the models, and the best convolved model is the one with the lowest $\chi^{2}_{\rm shot}$.

The process of choosing the best convolved image is described by Equation 2,

[TABLE]

The $\chi^{2}_{\rm shot}$ for each data image $D$, of which there are $i$, is calculated by finding the minimum of $\chi^{2}_{\rm shot}$ value over the set of PSFs $P$, of which there are $j$, over all of the pixels with the index $k$; where $\sigma_{k}^{2}$ is the variance over those pixels, accounting for readout noise, target shot noise and background shot noise. Note that we deliberately choose not to attempt to model speckles and adaptive optics phase noise as additional uncertainties, as it is in principle not a fundamental noise source in a well characterised and calibrated adaptive optics imager. We do, however, scale our final reported uncertainties appropriately, as described in Section 6.2.

The rotated model image is convolved with all of the PSFs, and the best of these is chosen to calculate the $\chi^{2}_{\rm shot}$ for the data image $D_{i}$.

These best $\chi^{2}_{\rm shot}$ values for each image are then summed to find the total $\chi^{2}_{\rm shot}$ value, $\chi^{2}_{\mathrm{tot}}$,

[TABLE]

This process of choosing the single best calibrator can in principle under-fit the data in a similar way that the locally optimized combination of images (LOCI) algorithm under-estimates the brightness of detected companions (Soummer et al., 2012). We checked that this type of systematic was not severe by verifying that many different calibrator images were chosen for differing target images. Given how flat the squared visibility versus baseline curve of the calibrators are when calibrated against each other (see Figures 3 and  4), this systematic is limited to missing flux in our model at up to the few percent level in the outer portion of the disc. A more robust algorithm would need to develop a significantly larger set of model PSFs (e.g. from a linear combination of calibrator images) and marginalise over all possible PSFs, as described in Ireland (2016).

4.4 Model parameter estimation via MC in MC

All of the aforementioned steps and the parameter estimation are completed using our MCinMC technique. MCinMC is the use of the Monte Carlo radiative transfer code RADMC3D in combination with the Markov Chain Monte Carlo (MCMC) code emcee (Foreman-Mackey et al., 2013). Given the nature of the large parameter space, and the risk that some of the parameters may be correlated, an MCMC approach was taken to ensure that the best-fit parameters were the global minimum, rather than a local minimum. emcee is an affine-invariant ensemble sampler, meaning that it is not affected by covariances between parameters. emcee uses a number of walkers to explore a parameter space. In this work the walkers are initialised from a group of points in the N-dimensional parameter space and the parameter space is explored. emcee can be easily made parallel onto different compute cores to make the computations more efficient.

Once the best $\chi^{2}_{\rm shot}$ for each individual data image and its corresponding rotated, convolved model image is calculated, the $\chi^{2}_{\rm shot}$ values are summed and used to calculate the log likelihood of the model;

[TABLE]

where the total $\chi^{2}_{\rm shot}$ is that calculated in Equation 3, and $T_{\mathrm{MC}}$ is the emcee temperature, set so that scaling the error bars by $\sqrt{T_{\mathrm{MC}}}$ would give a reduced $\chi^{2}$ of 1. We choose this approach rather than adding additional uncertainty estimates to our data variance, in order to show how much better the fits could be if PSF models were improved. emcee uses the log likelihood to find the best model for the data. The reduced $\chi^{2}_{\rm shot}$  is calculated from this log likelihood, given by:

[TABLE]

where the number of pixels in the image $n_{\mathrm{pixels}}$, and the number of images $n_{\mathrm{images}}$ normalise over the degrees of freedom, making this the reduced $\chi^{2}_{\rm shot}$. The number of pixels in the image ($n_{\mathrm{pixels}}$), is actually the number used in the calculation, rather than the total number of pixels in the image. The number of images ($n_{\mathrm{images}}$) is the number of target star images.

5 Synthetic Data Set

Our method was tested by generating a synthetic data set within MCinMC and then running a full parameter exploration on the synthetic data. The synthetic data set was generated by convolving a RADMC3D model, with parameters as listed in Table 5, and the values in the column labelled ‘Target’ in Table 6 with a sequence of calibrator (PSF) images from HD 167666. When this data set was tested a different sequence of calibrator images (HD 170768) were used for generating the convolved model for comparison with the data. It was found that the method could recover the parameters of the initial model to within 5 standard deviations for all parameters, with our approximate method of accounting for point-spread function fitting uncertainties as described in Section 6.2.

5.1 Testing the Synthetic Data Set with MCinMC

A model with a disc asymmetry and companion was used to generate the synthetic data set and then this data set was tested in the MCinMC code. The data set was tested as a symmetric, asymmetric and asymmetric with companion model. The results of this exploration are shown in Table 6, with figures of the results in Figure 8. The target parameters are those used to construct the data set, listed as Target in Table 6.

The symmetric model is the worst fit for the synthetic data, at a reduced $\chi^{2}_{\rm shot}$ of 67.5. For the symmetric model only the inner wall radius and the inclination are the same (within uncertainties) as those used to generate the data, however most were still within a few standard deviations of the true values. Since the data is asymmetric with a companion, it is not surprising that the symmetric model does not fit the data well. The $\chi^{2}_{\rm shot}$ for the symmetric model is more than twice that of the asymmetric model.

Most parameters for the asymmetric model are closer to those used to generate the data than the symmetric results are, and those closest to the parameters used are the two depletion factors, the radius of the inner disc, the radius of the first wall, the inclination, and the asymmetry parameters. The asymmetric model $\chi^{2}_{\rm shot}$ of 27.3 is a large improvement over the symmetric fit.

The companion model has an improved, but similar $\chi^{2}_{\rm shot}$ to the asymmetric model. The resulting fits to the 8–13 parameters for the companion fit are mostly within three standard deviations of the data. When the model is convolved with the PSF to generate the synthetic data, the partly-resolved companion gives a signal similar to a low-order aberration, and our simplified error model does not take this into account adequately. Note that in this case, the companion separation is dominated by the $x$ position (causing the slight deviations in the companion offset and position angle), which agrees reasonably well with the input model.

The residual images for both the asymmetric and companion models have obviously lower residuals than the symmetric model. The symmetric model residual shows that there is an asymmetry in the data that is not being fit with this model. The asymmetric and companion model residuals mostly show the noise in the PSFs, although there is a slight difference between the two at the companion position, showing that the companion model fits a feature where the asymmetric model cannot. The small uncertainty on the brightness of the companion and its position indicate that there is a strong asymmetry in that part of the disc.

5.2 Summary of Synthetic Data Set Results

In this test our method is able to distinguish between a symmetric or asymmetric disc, but not as well between an asymmetry caused by an eccentric disc or a point-like asymmetry. The symmetric model will not fit well to a disc with obvious asymmetries. It is difficult to determine the cause of an asymmetry with our method, however we are sensitive to whether an asymmetry is present, and the strength of the asymmetry.

Our method is limited by the uncertainties in the PSFs, making it difficult to determine the difference between a disc asymmetry and a point-like asymmetry in the disc. Given the faint nature of planetary companions it can be difficult to determine the difference between noise, disc asymmetries and a point-like asymmetry.

6 Results

6.1 Dust

It was initially found that carbon dust gave a better result (lower reduced $\chi^{2}_{\rm shot}$) than silicate dust; and that dust containing carbon and PAHs gave a better result than plain carbon dust. Carbon produced better models for Oph IRS 48 than silicate dust could, because silicates were too cool to produce the emission at $\sim$10 AU radii (Figure 9). Relative to the absorption of stellar photons at $\la 1$ $\mu$m wavelengths, silicate dust has a large emission feature at 11 µm which effectively cools the dust (Draine, 2003).

Using carbon dust, it was not possible to simultaneously have the $\sim$13 AU ring bright enough to fit our adaptive optics data and still have an inner disc bright enough to fit the SED. The amount of carbon dust needed in the inner disc to fit the SED would cause too much shadowing on the 13 AU ring. Therefore, a dust which could more effectively absorb blue and ultraviolet wavelengths was needed.

Both of the discs studied in this work have clear PAH features in their SEDs (e.g. Bruderer et al., 2014; Seok & Li, 2016; Maaskant et al., 2014, Oph IRS 48, HD 169142, both, respectively). In our attempts to fit the near-IR excesses seen in the SED, CP dust was used, and the same dust was used for both objects due to the similarities in their SEDs. Given that we are not trying to fit the whole SED, we chose a dust model from the Draine & Li (2001, 2007) dust, which includes fine grained graphitic carbon, and PAHs. From inspection of the SEDs the dust chosen from the Draine & Li (2001, 2007) set of neutral CP dust has a grain-size of 5.6 nm. Even though it is known that some of the PAHs in these discs are ionised (Maaskant et al., 2014), we chose the neutral dust to use in both cases because it was similar to the PAH emission of the objects and able to replicate the near-IR excess. For an investigation of how differently sized dust grains from the Draine & Li (2001, 2007) set of neutral CP dust change the fit to the SED, see the Appendix (D).

The reason that the smaller grains (CP dust) better recover the near-IR excess is that they are generally warmer than larger grains (see Figure 9), due to their relatively higher UV absorption. There is also no need to invoke quantum heating in this case, as $\sim 5$ nm radius dust has a sufficient number of atoms to not be significantly warmed by single photon absorption events. Thus quantum heating in the disc is ignored here, though it may become required to model images at all wavelengths (especially even shorter wavelengths) simultaneously.

One consideration to be able to have both a bright inner disc and bright walls at 13 and $30\,\mathrm{AU}$ (as is needed for Oph IRS 48, and similar to what is needed in HD 169142) would be to have in inclined inner disc. With an inclined inner disc, carbon dust could still be used and the disc walls would still be bright. However, for both Oph IRS 48 and HD 169142, the deconvolved images show symmetric features, with no evidence for the shadows that would be expected at opposite sides of the disc if the inner disc was inclined. Therefore an inclined inner disc with carbon dust is not considered here and we instead use CP dust.

A mix of different dust types would benefit the modelling of the disc and the SED simultaneously. More freedom in the dust distribution and composition would allow for better modelling of disc features, and possibly determine the difference between a disc feature or a companion. This was not justified for our current method, given the high speckle residuals, which were at the same level for our data sets as for synthetic data sets using multiple calibrator stars.

6.2 Fits to Oph IRS 48 and HD 169142

The best-fit model parameters for both objects are shown in Table 7. The primary fitting for HD 169142 was done for the 2014 data and for Oph IRS 48 it was done for the 2015 data, the fit for the best fit model for each of these data sets is then tested on the 2016 data. We do not report independent fits to the 2016 data because as mentioned in section 2.2, some observations were taken through cloud.

In the following sections the results of each model type for each object will be discussed. For each object there was a symmetric, an asymmetric and a companion model tested. The model setups are described in Section 4.2. Figures 10 and 11 show each of these models; their convolved counterpart; the residual of the data and the model; and the ratio of the data and the model, for Oph IRS 48 and HD 169142 respectively, and are discussed in Sections 6.3, and 6.4.

The values and uncertainties quoted in Table 7 are the mean and standard deviation of the results found with emcee. The data used in this work has uncertainties with several components, including background and speckle noise.

We choose to only consider the target and background shot noise to give the reader a clear indication on the scope of improvement in model fitting that would be possible with an improved point-spread function model (e.g., interpolation between PSFs). However, we also need to robustly estimate statistical uncertainties from our current fitting method, and choose to do this by setting the temperature in the Monte-Carlo Markov Chain.

Naively, emcee temperature would be increased by the value of reduced $\chi^{2}_{\rm shot}$, which is equivalent to scaling error bars on the data pixels by $\sqrt{\chi^{2}_{\rm shot}}$ to give $\chi^{2}=1$ due to uncertainties in addition to the modelled shot noise. However, this does not take into account the highly correlated speckle noise, and would underestimate our final uncertainties. The speckle noise has a characteristic solid angle of correlated noise which is $\sim(\lambda/D)^{2}$, or $\sim N_{\rm ind}=70$ pixels. Therefore, to account for the speckle noise we set the temperature in emcee ($T_{\mathrm{MC}}$) as $\sim N_{\rm ind}\times\chi^{2}_{\rm shot}$ ($\sim 10000$), where $\chi^{2}_{\rm shot}$ is the reduced chi-squared of the best fitting model. Thus the $\chi^{2}_{\rm shot}$ of $\gg 1$ we have in our results means that photon statistics of the target and background do not dominate the uncertainty. We have not studied sufficient calibrator-calibrator pairs to re-define our uncertainties reliably (based on the noise from the speckles).

The main purpose of this work was to fit the inner regions of the disc, rather than the whole disc. The interested reader is directed to Appendix D for an investigation of the SEDs of the discs.

6.3 Oph IRS 48

For Oph IRS 48 the aim was to study the disc structure within $\sim 20\,\mathrm{AU}$ and also to see if a companion could be detected. In our analysis of Oph IRS 48, the 2015 data were used, with the best model for each of the model types for the 2015 data tested on the 2016 data. Figure 10 shows the resulting models, convolved images, residuals and ratio images.

One commonality shared by the Oph IRS 48 models in their residual and ratio images in Figure 10 is the two features (one on the left and one on the right of the disc) that are consistently present. These ‘ear’ features are green in the residuals and light green in the ratio images. They are present in the same location as where the model appears brighter at around $30\,\mathrm{AU}$ ($\sim\pm 0.25\,$″), because of the overlap of the bright $\sim 30\,\mathrm{AU}$ ring.

This region of the disc is expected to become optically thick with opacity dominated by large grains, but our model only includes very small grains as we focus on the inner regions. This issue with the ‘ears’ in the residual and ratio images could possibly be rectified by using a mixture of dust grains with varying sizes and compositions. There are some other features in the residual images that are consistent across all of the models, which may be due to an inadequate PSF library. As mentioned above when uncertainties were discussed, this noise is also the reason for the high reduced $\chi^{2}_{\rm shot}$ values.

All three of the Oph IRS 48 models have similar values for the parameters. The symmetric and asymmetric models are most similar with nearly all parameters matching within uncertainties. The companion model is slightly different; as it must allow for the presence of the companion changing the brightness of the disc, thus adjusting the dust parameters and stellar position on accordingly. The similarities and differences between the models are discussed in the following sections.

6.3.1 Symmetric Disc

The first column of Figure 10 shows the results for the symmetric model. The residual and ratio images show how well the model fits the data. The L’ band flux ratio of 7.8 is a match for the value we were fitting to. The symmetric model of Oph IRS 48 has a bright inner disc, $r_{1}$ at $\sim 112\,$mas ($\sim 15\,\mathrm{AU}$) and $r_{2}$ at $\sim 238\,$mas ($\sim 32.1\,\mathrm{AU}$). The location of $r_{1}$ is consistent with the location of the PAH emission modelled by Geers et al. (2007b, $\sim 16\,\mathrm{AU}$), and the rings suggested by Brown et al. (2012b, $\sim 15\,\mathrm{AU}$ and $\sim 34\,\mathrm{AU}$).

6.3.2 Asymmetric Disc

As mentioned before, the asymmetry in the disc was introduced by moving the star. This introduces two new parameters that the symmetric model did not have, the star offset and position angle (generated from the $x$ and $y$ position of the star). Because there are additional parameters a better fit is expected, and in this case a slightly better fit is achieved for the 2015 data.

Aside from the two new parameters, most of the parameter values are the same within uncertainties as for the symmetric model case.

6.3.3 Asymmetric Disc with Companion

As mentioned earlier, a companion is added to the asymmetric model by including a second RADMC3D point source. The results for the planet model differ slightly from the asymmetric and symmetric model for the parameters relating to the dust structure. This is because having the second source causes some changes to the brightness of the disc. The values of $r_{1}$ and $r_{2}$ are similar to those of the symmetric and asymmetric models, but $r_{d}$ has moved slightly further out.

The companion has a brightness contrast of 3.97$\pm{0.05}\,$mag with respect to the rest of the image. Note that this is a factor of $\sim$30 fainter than the disc emission or a factor of $\sim$5 fainter than the local disc emission within a diffraction limit ($\sim$80 mas). The position of our point-like asymmetry at $\sim 104\,$mas and $\sim 288\,^{\circ}$ is consistent with one of the point sources found by Schworer et al. (2017) at $\sim 105\,$mas to the west of the star. However, the Schworer et al. (2017) object has a contrast of $\sim 3.3\,$mag.

The reduced $\chi^{2}_{\rm shot}$ for this more complex model is the lowest of the three models for both 2015 and 2016, as reported in Table 7.

6.3.4 Summary of Oph IRS 48 Results

The disc parameters are similar across all models for Oph IRS 48, however there is improvement with the addition of an asymmetry and a companion. Given this improvement, the small uncertainty in the companion radius and the narrow brightness contrast range, it is likely that there is an asymmetry in this region of the Oph IRS 48 disc. However, we are unable to identify whether this is due to a disc asymmetry or a point-like companion at this stage.

In all cases the inner disc radius ($r_{d}$) is greater than what is considered to be the canonical inner disc for Oph IRS 48. Schworer et al. (2017) find using imaging and SED fitting that there is likely no emission from the disc between 0.4 and $1\,\mathrm{AU}$, and that the very small particles in the disc are present from $11\,\mathrm{AU}$. We find that for the disc to be bright enough to match our observations, we do still require the inner disc, but present beyond $1\,\mathrm{AU}$, with a second ring of emission present at $\sim 15\,\mathrm{AU}$ ($\sim 112\,$mas).

Generally for Oph IRS 48 the $\chi^{2}_{\rm shot}$ values are quite high. Some part of this high $\chi^{2}_{\rm shot}$ is likely not only due to an inadequate PSF library, but also due to insufficient complexity in the model, as evidenced by the negative residuals to the east and west at $\sim 0.25\,\arcsec$ separation.

6.4 HD 169142

One goal of the study on HD 169142 was to see if we could detect the candidate companion proposed to be in the inner disc. For HD 169142, the 2014 data were used for the analysis, with the best model for each of the model types for the 2014 data tested on the 2016 data. The results of each model type are discussed below, and the figures corresponding to the model are shown in Figure 11. The same model types and figures are included here as were in the previous section (6.3) for Oph IRS 48.

The inclination and position angle in all cases for the HD 169142 disc are different to those in the literature. The inclination we find here is $\sim 30\,^{\circ}$, rather than the $\sim 13\,^{\circ}$ reported in Panić et al. (2008), and our position angle is $\sim 15\,^{\circ}$ compared to their $\sim 30\,^{\circ}$. This difference may arise from the large difference in the spatial scales probed: the Panić et al. (2008) observations using rotational lines of CO probed $>1\,$″separations compared to our $\sim 0.1\,$″separations. A disc warp, for example, could produce these slightly differing inner and outer disc geometries.

6.4.1 Symmetric Disc

For the symmetric model, the best fit has $r_{1}$ at $74\,$mas ($\sim 8.4\,\mathrm{AU}$), with $r_{2}$ at $135\,$mas ($15.4\,\mathrm{AU}$). Previous studies of the disc found a ring at $\sim 175\,$mas ($\sim 20\,\mathrm{AU}$), using modelling of the SED based on mid-IR observations (Honda et al., 2012) and imaging of large grains from 7mm observations (Osorio et al., 2014). Our deconvolved image of HD 169142 suggested a single ring at $\sim 88\,$mas ($\sim 10\,\mathrm{AU}$), similar to that marginally detected by Ligi et al. (2018, $\sim 100\,$mas). We attempted modelling HD 169142 with one wall, but found that the $\chi^{2}_{\rm shot}$ values were significantly lower when we used two walls.

The reduced $\chi^{2}_{\rm shot}$ for the 2014 data was 5.91 and for the 2016 data it was 44.1. There is more of a discrepancy between the results than there was for Oph IRS 48, likely due to the very significant difference in the observing conditions in this case (flux varying by a factor of $\sim$2 during observations even after removing the most cloud-affected data) and this $\chi^{2}$ metric only taking shot noise into account (Section 4.3).

6.4.2 Asymmetric Disc

The results for the asymmetric model are similar to the symmetric model. The results fit the symmetric ones within a few standard deviations. For this model $r_{1}$ is at $64\,$mas ($\sim 7.3\,\mathrm{AU}$), with $r_{2}$ at $118\,$mas ($\sim 13.5\,\mathrm{AU}$). The location of the star to generate the asymmetry is $3-\sigma$ from the origin, suggesting that the star may not be off centre in HD 169142. The reduced $\chi^{2}_{\rm shot}$ value for the asymmetric model is 5.58, which is a slight improvement over the symmetric model. The reduced $\chi^{2}_{\rm shot}$ of the 2016 data is 42.8 using the asymmetric model.

6.4.3 Asymmetric Disc with Companion

The companion model has a similar density structure to both the symmetric and asymmetric models. There is a wall at $64\,$mas ($7.3\,\mathrm{AU}$), similar to the wall found in the symmetric and asymmetric models, as well as a wall at 118 mas ($\sim 13.4\,\mathrm{AU}$). When the companion position was freely explored, it was not detected in any location with certainty. However, when the companion is placed at a location consistent with those found in the Biller et al. (2014) and Reggiani et al. (2014) studies (these are consistent with each other within the uncertainties), the $\chi^{2}_{\rm shot}$ was improved and the brightness contrast was found to be $5.7\pm{0.7}\,$mag. Our companion is at a lower significance than those papers, and not as significant as a general disc asymmetry. Our brightness contrast is roughly consistent with the $6.5\pm 0.5\,$mag found by Reggiani et al. (2014), and the $\sim 6.4\,$mag reported by Biller et al. (2014). This point-like asymmetry is a factor of $\sim$100 fainter than the disc emission or a factor of $\sim$20 fainter than the local disc emission within a diffraction limit ($\sim$80 mas).

For the 2016 data the reduced $\chi^{2}_{\rm shot}$ is 45.7, which is slightly higher than for the asymmetric model.

6.4.4 Summary of HD 169142 Results

In summary the results of the HD 169142 models are consistent with each other, and indicate the presence of an asymmetry in the disc. We are able to recover the previously detected companion-like asymmetry, but are unable to determine the cause of the asymmetry at this stage. We also detect symmetric rings of dust at $\sim 61\,$mas ($\sim 7\,\mathrm{AU}$) and $\sim 114\,$mas ($\sim 13\,\mathrm{AU}$). Our rings are at a similar location to the $\sim 100\,$mas ($\sim 11\,\mathrm{AU}$) ring and asymmetry found by Ligi et al. (2018). The location of our inner disc ($r_{d}$) is at $0.11\,\mathrm{AU}$ (0.96 mas), which is slightly beyond the radius of $0.07\,\mathrm{AU}$ found by Chen et al. (2018) from their modelling of the inner disc.

The substantial difference between the $\chi^{2}_{\rm shot}$ values for the 2014 and 2016 epochs for HD 169142 is due to the large variation in observing conditions between epochs. As discussed in Section 2.2, the observing conditions for the 2016 epoch were affected by weather, particularly HD 169142.

6.5 Additional Companions

We attempted to fit additional wide companions to our model residuals (i.e. companions present in addition to the companions fit in Sections 6.3.3 and 6.4.3), assessing the significance of companion fits by comparing fitted companion brightness to the root-mean square azimuthally averaged contrast as a function of separation. No additional companions were detected, and the resulting contrast limits are shown in Figure 12, with contrasts between 5 and 7 magnitudes with respect to the system absolute L’ magnitudes of $\sim-0.6\,$mag for Oph IRS 48 and $\sim 0.6\,$mag for HD 169142. The relatively poor contrast for Oph IRS 48 is likely due to the inadequacy of our model in fitting the disc structure, resulting in relatively high residuals. Additional companions would be detectable if they exist in the space below the lines in Figure 12.

The possible 3.5 $\rm{M_{J}}$ object at $40\,\mathrm{AU}$ ($0.3\,\arcsec$) discussed in Schworer et al. (2017), deduced from disc structure and not from direct planetary emission, would then not be detectable in our data, unless it had a very high accretion rate of $\sim 3\times 10^{-5}\,\rm{M_{J}}\,\rm{yr}^{-1}$ according to the models of Zhu (2015). Note that our companion model is consistent with the asymmetry seen by Schworer et al. (2017) which is interpreted as an asymmetric ring. Potential companions aside from those near the location of the Biller et al. (2014), Reggiani et al. (2014) and Ligi et al. (2018) candidates in HD 169142 include those suggested by Pérez et al. (2019), which is at beyond $0.5\,$″; Gratton et al. (2019), at a distance of $0.335\,$″, but a contrast of $10.1\,$mag; and Pohl et al. (2017), one of which is at our inner limit and the other beyond our limit; thus none of these would be detectable with this method.

7 Conclusions

Here we introduced our method for the study of transitional discs. The method utilises images of target and calibrator stars taken with a thermal IR filter behind adaptive optics to reveal detail about the structures of their discs. We use our new MCinMC (see Section 4) code to make models of the observed discs, then convolve the models with calibrator images, and finally compare the convolved models with the data to determine how well the model fits our data. The key conclusions we draw from this work are summarised here:

In contrast to the historical view of transitional discs, we find that the “gap” region in two Herbig Ae stars is radiatively dominated by emission from very small ($\la 5\,$nm) grains and PAHs. This has been suggested before for both objects, and this paper confirms their presence with spatially resolved detections. 2. 2.

We confirm previously reported brightness asymmetries (Biller et al., 2014; Reggiani et al., 2014; Schworer et al., 2017) in both HD 169142 and Oph IRS 48 using our complementary technique. However, the detection of overall disc asymmetry (modelled as an offset central star) was more significant than a point source companion model. For both objects, the co-located disc emission (within a diffraction limit) was brighter than previously reported “companions” by a factor of $\sim$5 or more, leading us to conclude that there is no need to invoke companions as anything more than a modelling convenience to explain the asymmetric emission.

The structure we detect in these discs suggests that the most common explanation of transitional discs having cleared inner holes may be too simple. To better understand these transitional discs and the mechanisms within them, we need more complete models of the disc geometry and composition and an improved understanding of the PSF uncertainties in the observational data.

An improvement that could be made to our method would be to interpolate between the PSFs and to use a larger library, possibly incorporating multiple nights of observations or models of optical aberrations, in a similar way to LOCI (Lafrenière et al., 2007) or KL eigenimages (Soummer et al., 2012). These methods will ultimately be limited by the angular resolution of a single telescope, and improved spatial resolution with MATISSE (Lopez et al., 2014) or future concepts such as the Planet Formation Imager (PFI, Ireland et al., 2016) may be required to definitively distinguish between disc features and signs of exoplanets.

Acknowledgements

We would like to thank the anonymous reviewer for their comments, which greatly improved the quality of this paper.

The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

E.K.B. would like to thank the Australian Government for their support through the Australian Government Research Training Program Stipend Scholarship and the Research School of Astronomy and Astrophysics at the Australian National University for the Masters of Astronomy and Astrophysics (Advanced) Scholarship.

M.J.I. gratefully acknowledges funding provided by the Australian Research Council’s Future Fellowship (FT130100235).

C.F. gratefully acknowledges funding provided by the Australian Research Council’s Discovery Projects (grants DP150104329 and DP170100603) and Future Fellowship Scheme (grant FT180100495), as well as the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). The Monte Carlo simulations and data analyses presented in this work used high performance computing resources provided by the Australian National Computational Infrastructure (grant ek9), and the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, in the framework of the National Computational Merit Allocation Scheme and the ANU Allocation Scheme. We further thank for supercomputing resources at the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391) and the Partnership for Advanced Computing in Europe (PRACE grant pr89mu).

S.K. acknowledges support from an ERC Starting Grant (Grant Agreement No. 639889) and STFC Rutherford Fellowship (ST/J004030/1).

These data were collected thanks to support from NASA KPDA grants (JPL-1452321, 1474717, 1485953, 1496788) and J.D.M. acknowledges NSF AST. 1311698.

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.

Appendix A Parameter Convergence Study

A convergence test was completed to determine the optimal number of photon packets to use in the model, so that it was stable, but not too computationally expensive. It was found that $10^{6}$ photon packets was good choice for the RADMC3D thermal Monte Carlo to converge. As seen in the plot in Figure 13 there is a minimum for reduced $\chi^{2}_{\rm shot}$ for the best model found with $10^{4}$ or $10^{5}$ photon packets, because increasing the number of photon packets changed the best fit parameters.

The convergence test used half and double the number of points in each parameter for the number of grid points. The result was that these different grids give a very similar $\chi^{2}_{\rm shot}$ values to the one shown in Table 8.

Appendix B Parameters for the companion brightness

In all cases the companion had a temperature of 3500K and a mass of $10^{-3}\,\mathrm{M}_{\odot}$, these were chosen due to modelling constraints. The brightness of the companion was adjusted by changing its radius, with a larger radius making the companion brighter. There was a lower limit on the radius of the companion, $0.02\,\mathrm{R}_{\odot}$, which had a negligible amount of flux, and we found that this was not important to the determination of the log likelihood.

Appendix C Additional RADMC3D Parameters

Geometric parameters for the spatial grid are shown in Table 8, and additional parameters that were used (many of which did not deviate from the default RADMC3D values) are shown in Table 9.

Appendix D Spectral Energy Distributions and Dust Types

The SED of the symmetric disc model for Oph IRS 48 (left) and HD 169142 (right) are shown in Figure 14. The SED models are reddened to match the photometry of each of the objects, and are included to show that without explicitly trying to fit the SED we are able to replicate the shape. The references for the data used in the SEDs can be found in Table 10.

For Oph IRS 48 it is difficult to de-redden the photometry or redden the model, because there is so much extinction towards it. We use the Cardelli et al. (1989) reddening laws with $R_{V}=6.5$ and $A_{V}=12.9$, as was found in Schworer et al. (2017). To calculate the level of reddening that was to be applied we used the extinction package for python by Barbary (2016).

We note that the SED is not fit very well by our model, and this is due to the constraints we have placed on ourselves for modelling the disc, as well as the large amount of extinction towards Oph IRS 48. A better fit could be achieved by differing dust types with radius, or with a different model configuration. As fitting the SED was not the key goal of this work, but fitting the Keck L’ filter data was, we leave detailed fitting of the SED for future work.

The SED of the symmetric model for HD 169142 (right of Figure 14) is able to recover the flux at our target wavelength, but does not do well beyond this, particularly not in the un-modelled cool outer disc region.

The model was reddened using the same methods as for Oph IRS 48, but with $R_{V}=3.1$ and $A_{V}=0.31$. A better fit could be achieved with a different model configuration, specifically one that allows for varying the dust chemistry and distribution of the disc.

We also tested HD 169142 with a few different dust sizes to investigate how grain-size changed the fit to the SED. The results of these tests are found in Figure 15. Each dust type was started from the best fit with the $5.6\,$nm dust and then the MCinMC code was run to find a best fit model for the new dust type. The grain-sizes investigated were $5.0\,$nm, $6.3\,$nm, and $10.0\,$nm, with the $5.6\,$nm dust to compare to. The composition of the dust is all the same carbon and PAH dust mix as discussed above from the Draine & Li (2001, 2007) set of neutral carbon and PAH dust.

Dust smaller that the one used for our analysis has very strong emission from PAHs, whereas dust of a larger grain size has weaker PAH emission. Larger dust is too cool to recover the features we would like to replicate. The chosen dust size of $5.6\,$nm allows us to have PAH features with some smoothing due to graphite. Using a significantly smaller dust would move into the realm of needing quantum heating of the dust particles, which would be difficult to parameterise for our model, so we chose to use dust models readily available from the literature.