Dusty disk winds at the sublimation rim of the highly inclined, low mass YSO SU Aurigae

TL;DR

This study uses high-resolution infrared interferometry and advanced modeling to reveal a dusty disk wind around the young star SU Aurigae, providing new insights into disk structure and dust dynamics near the sublimation rim.

Contribution

First implementation of dusty disk wind in radiative transfer modeling of protoplanetary disks, explaining K-band excess and disk features around SU Aurigae.

Findings

Inner disk rim at ~0.17 au with inclination ~51 degrees

Detection of a dusty disk wind accounting for K-band excess

Inclined, flared disk with specific grain size and scale height

Abstract

T Tauri stars are low-mass young stars whose disks provide the setting for planet formation. Despite this, their structure is poorly understood. We present new infrared interferometric observations of the SU Aurigae circumstellar environment that offer 3 x higher resolution and better baseline position angle coverage over previous observations. We investigate the characteristics of circumstellar material around SU Aur, constrain the disk geometry, composition and inner dust rim structure. The CHARA array offers opportunities for long baseline observations, with baselines up to 331 m. Using the CLIMB 3-telescope combiner in the K-band allows us to measure visibilities as well as closure phase. We undertook image reconstruction for model-independent analysis, and geometric modeling. Additionally, the fitting of radiative transfer models constrains the physical parameters of the disk. For the first time, a dusty disk wind is introduced to the radiative transfer code TORUS to model protoplanetary disks. Our implementation is motivated by theoretical dusty disk winds, where magnetic field lines drive dust above the disk plane close to the sublimation zone. Image reconstruction reveals an inclined disk with slight asymmetry along its minor-axis, likely due to inclination effects obscuring the inner disk rim through absorption of incident star light on the near-side and thermal re-emission/scattering of the far-side. Geometric modelling of a skewed ring finds the inner rim at 0.17+/-0.02 au with an inclination of 50.9+/-1.0 degrees and minor axis position angle 60.8+/-1.2 degrees. Radiative transfer modelling shows a flared disk with an inner radius at 0.18 au which implies a grain size of 0.4 um and a scale height of 15.0 au at 100 au. Among the tested radiative transfer models, only the dusty disk wind successfully accounts for the K-band excess by introducing dust above the mid-plane.

Tables

Table 1. Table 1: Stellar parameters of SU Aurigae

Parameter	Value	Reference
RA (J2000)	$\mathrm{04\; 55\; 59.39}$	(1)
DEC (J2000)	$+\mathrm{30\; 34\; 01.50}$	(1)
Mass	$2.0M_{\odot }$	(2)
Sp. Type	G2 IIIe	(2)
Distance	${157.68}_{-1.48}^{+1.49}\mathrm{pc}$	(1)
Radius	$3.5R_{\odot }$	(2)
$K_{mag}$	$5.99$	(3)
$T_{\mathrm{eff}}$	$5860K$	(2)

Table 2. Table 2: Observing log from 1999 to 2014 from the CHARA, KI and PTI interferometers.

Date	Beam Combiner	Stations	Pointings	Calibrator (UD [mas])
2010-10-02	CHARA/CLIMB	S1-E1-W1	2	HD 29867 ($0.280\pm 0.007$), HD 34499 ($0.257\pm 0.006$)
2010-12-02	CHARA/CLIMB	S2-E1-W2	1	HD 32480 ($0.236\pm 0.006$)
2010-12-03	CHARA/CLIMB	(S2)-E1-W2	1	HD 32480 ($0.236\pm 0.006$), HD 36724 ($0.233\pm 0.006$)
2012-10-18	CHARA/CLIMB	S1-E1-W1	2	HD 27777 ($0.204\pm 0.006$), 34053∗*∗*$\ast$Calibrator found to be a binary, solution used shown in Section A.
2012-10-19	CHARA/CLIMB	E2-S1-W2	4	HD 27777 ($0.204\pm 0.006$), HD 31592∗*∗*$\ast$Calibrator found to be a binary, solution used shown in Section A., HD 34053∗*∗*$\ast$Calibrator found to be a binary, solution used shown in Section A.
2012-10-20	CHARA/CLIMB	S1-W1-W2	3	HD 31592∗*∗*$\ast$Calibrator found to be a binary, solution used shown in Section A., HD 34053∗*∗*$\ast$Calibrator found to be a binary, solution used shown in Section A.
2012-11-27	CHARA/CLIMB	S1-E1-W1	5	HD 32480 ($0.236\pm 0.006$), HD 31706 ($0.219\pm 0.005$)
2012-11-28	CHARA/CLIMB	S1-E1-E1	4	HD 32480, ($0.236\pm 0.006$) HD 31706 ($0.219\pm 0.005$),
				HD 33252 ($0.294\pm 0.007$)
2014-11-25	CHARA/CLIMB	E2-S2-W2	3	HD 33252 ($0.294\pm 0.007$)
2014-11-26	CHARA/CLIMB	E2-S2-W2	3	HD 33252 ($0.294\pm 0.007$)
2009-10-31	CHARA/CLASSIC	S1-E1	4	HD 32480 ($0.236\pm 0.006$)
2009-11-01	CHARA/CLASSIC	S1-E1	3	HD 32480 ($0.236\pm 0.006$), HD 24365 ($0.319\pm 0.008$)
2011-11-07	KI	FT-SEC	5	HD 27777 ($0.204\pm 0.006$)
1999-10-09	PTI	NS	2	HD 30111 ($0.555\pm 0.056$)
1999-11-03	PTI	NS	2	HD 28024 ($0.222\pm 0.026$), HD 27946 ($0.398\pm 0.034$),
				HD 25867 ($0.280\pm 0.007$)
1999-11-04	PTI	NS	1	HD 28024 ($0.222\pm 0.026$), HD 32301 ($0.506\pm 0.054$),
				HD 25867 ($0.280\pm 0.007$)
1999-12-07	PTI	NS	6	HD 27946 ($0.398\pm 0.034$), HD 30111 ($0.555\pm 0.056$)
2000-10-16	PTI	NS	4	HD 30111 ($0.555\pm 0.056$)
2000-11-13	PTI	NW	17	HD 30111 ($0.555\pm 0.056$)
2000-11-14	PTI	NW	9	HD 30111 ($0.555\pm 0.056$)
2003-10-16	PTI	SW	6	HD 28024 ($0.222\pm 0.026$), HD 30111 ($0.555\pm 0.056$)
2003-10-22	PTI	SW	4	HD 28024 ($0.222\pm 0.026$), HD 30111 ($0.555\pm 0.056$),
				HD 27946 ($0.398\pm 0.034$), HD 25867 ($0.508\pm 0.046$),
				HD 29645 ($0.507\pm 0.035$)
2004-10-01	PTI	NW	2	HD 29645 ($0.507\pm 0.035$)
2004-10-05	PTI	SW	1	HD 29645 ($0.507\pm 0.035$)

Table 3. Table 3: Best fit parameters for the simple geometric models investigated.

Parameter	Explored Parameter Space	Gaussian	Ring	Skewed Ring
$R_{\min }$ [mas]	$0.0-10.0$	…	$1.0\pm 0.13$	$1.1\pm 0.12$
FWHM	$0.0-15.0$	$2.01\pm 0.02$	…	…
$w$	$0.01-2.0$	…	…	$0.96\pm 0.02$
$\theta$ [∘] (PA)	$0.0-180.0$	$60.8\pm 1.15$	$61.36\pm 1.2$	$61.0\pm 1.0$
$i$ [∘]	$0.0-90.0$	$51.4\pm 1.04$	$50.91\pm 0.88$	$51.2\pm 1.1$
${\chi }_{\mathrm{red},\mathrm{Vis}}^2$	…	$10.97$	$11.86$	$8.57$
${\chi }_{\mathrm{red},\mathrm{CP}}^2$	…	$1.90$	$1.90$	$1.66$∗*∗*$\ast$The closure phase quoted is the achieved when allowing the skewed ring to become asymmetric. While the software did detect an asymmetry, it failed to constrain its location in the disk. $\theta$ is the minor-axis position angle of the disk.

Table 4. Table 4: Best-fit parameters resulting from SED and visibility fitting.

Parameter	Literature value	Reference	Range explored	Best fit value
$i$	${60}^{\circ }$	(1),(2)	$30-{80}^{\circ }$	${50}^{\circ }$
$R_{\mathrm{inner}}$	$0.18$	(1)(3)	$0.1-0.6\mathrm{au}$	$0.15\mathrm{au}$
$R_{\mathrm{outer}}$	$100\mathrm{au}$	(1)	$20.0-120.0\mathrm{au}$	$100.0\mathrm{au}$
$h_0$	…	…	$7.0-20.0\mathrm{au}$	$15.0\mathrm{au}$
${\alpha }_{\mathrm{disk}}$	…	…	$1.0-3.0$	$2.4$
${\beta }_{\mathrm{disk}}=\left(\alpha -1\right)$	…	…	$0.0-2.0$	$1.4$
$\mathrm{Dust}:\mathrm{Gas}$	$0.01$	(1)	$0.01-0.008$	$0.01$
$a_{\min }$	$0.1\mu m$	(1)	$0.1-1.4\mu m$	$0.39\mu m$
$q_{\mathrm{dist}}$	$3.0$	(1)	$2.00-4.00$	$3.06$
$T_{\mathrm{sub}}$	$1600$	(1)	$1400-2000$	$1600K$
Dusty disk wind parameter	Literature value	Reference	Range explored	Best fit value
$R_{0min}$	…	…	$2.0-10.0R_{☉}$	$4.5R_{☉}$
$T_{\mathrm{wind}\left(\mathrm{near}\mathrm{surface}\right)}$	…	…	$1400-2400K$	$1600K$
$\mathrm{Opening}\mathrm{Angle}$	…	…	$25-{55}^{\circ }$	${45}^{\circ }$
$\dot{M}$	…	…	${10}^{-6}-{10}^{-12}{M}_{\odot }{\mathrm{yr}}^{-1}$	${10}^{-7}{M}_{\odot }y{r}^{-1}$

Table 5. Table 5: χ 2 superscript 𝜒 2 \chi^{2} results of radiative transfer fits to the squared visibilities, the closure phases (CP) and the SEDs.

Model	${\chi }_{V^2}^2$	${\chi }_{CP}^2$	${\chi }_{SED}^2$
IN05	91.9	0.1	151.1
THM07	84.2	0.1	137.7
BK12	35.6	0.1	121.3

Table 6. Table 6: Binary fit parameters for the two calibrator stars HD 31952 and HD 34053.

Star	Obs. Date	Flux Ratio	Separation [mas]	Position Angle [∘]	${\mathrm{UDD}}_{\mathrm{pri}}$ [mas]	${\mathrm{UDD}}_{\sec }$ [mas]
HD 31952	2012-10-(19,20)	$3.76\pm 0.07$	$7.03\pm 0.12$	$95.55\pm 2.35$	$0.15\pm 0.02$	$0.15\pm 0.03$
HD 34053	2012-10-(18,19,20)	$1.17\pm 0.05$	$1.66\pm 0.23$	$85.82\pm 3.27$	$0.14\pm 0.02$	$0.16\pm 0.04$

Table 7. Table 7: Photometric values used to construct the SED of SU Aur

Wavelength [$\mu m$]	Flux [Jy]	Reference
0.15	1.31E-04	Bianchi et al. (2011)
0.23	0.00219	Bianchi et al. (2011)
0.42	0.233	Ammons et al. (2006)
0.44	0.402	Anderson & Francis (2012)
0.53	0.6	Ammons et al. (2006)
0.69	1.32	Morel & Magnenat (1978)
0.79	1.47	Davies et al. (2014)
0.88	1.75	Morel & Magnenat (1978)
1.24	2.08	Röser et al. (2008)
1.25	2.12	Ofek (2008)
1.63	2.47	Ofek (2008)
2.17	2.71	Röser et al. (2008)
2.19	2.62	Ofek (2008)
3.35	2.6	Cutri & et al. (2014)
3.40	2.44	Bourgés et al. (2014)
4.50	1.75	Esplin et al. (2014)
4.60	2.78	Cutri & et al. (2014)
5.03	2.58	Bourgés et al. (2014)
7.88	1.99	Esplin et al. (2014)
8.62	2.36	Abrahamyan et al. (2015)
11.57	2.83	Cutri & et al. (2014)
11.60	3.52	Abrahamyan et al. (2015)
18.40	6.47	Abrahamyan et al. (2015)
22.11	9.24	Cutri & et al. (2014)
23.90	12.8	Abrahamyan et al. (2015)
61.89	12.2	Abrahamyan et al. (2015)
65.04	9.89	Tóth et al. (2014)
90.06	8.8	Tóth et al. (2014)
140.10	10.2	Tóth et al. (2014)
160.11	8.88	Tóth et al. (2014)
849.86	0.074	Mohanty et al. (2013)
887.57	0.071	Andrews et al. (2013)
1300.90	0.03	Mohanty et al. (2013)
1333.33	0.0274	Andrews et al. (2013)

Equations

$$\centering\rho(r,z)=\frac{\Sigma(r)}{h(r)\sqrt{2\pi}}\exp\bigg{\{}-\frac{1}{2}\bigg{[}\frac{z}{h(r)}\bigg{]}^{2}\bigg{\}}.\@add@centering$$

$$T_{sub}=G{\rho }^{\gamma }(r,z),$$

$$\rho=\rho_{1}\bigg{(}\frac{r_{0}}{r_{1}}\bigg{)}^{3/2}\eta(\chi),$$

$$\chi =z/r.$$

$$\rho=1.064\times 10^{-15}\Big{(}\frac{\dot{M}_{out}}{10^{-7}M_{\odot}yr^{-1}}\Big{)}\Big{(}\frac{M_{*}}{0.5M_{\odot}}\Big{)}^{-1/2}\\ \times\frac{1}{\mathrm{ln}(r_{0max}/r_{0min})\psi_{0}(1-h_{0}\xi_{0}^{{}^{\prime}})}\mathrm{\leavevmode\nobreak\ g\leavevmode\nobreak\ cm^{-3}},$$

$${\xi }_0^{{}^{\prime }}=tan({\theta }_0).$$

Full text

*Aims. *We aim to investigate the characteristics of the circumstellar material around SU Aur, constrain the disk geometry, composition and inner dust rim structure.

*Methods. *The CHARA array offers unique opportunities for long baseline observations, with baselines up to $331$ m. Using the CLIMB 3-telescope combiner in the K-band allows us to measure visibilities as well as closure phase. We undertook image reconstruction for model-independent analysis, and fitted geometric models such as Gaussian and ring distributions. Additionally, the fitting of radiative transfer models constrain the physical parameters of the disk. For the first time, a dusty disk wind is introduced to the radiative transfer code TORUS to model protoplanetary disks. Our implementation is motivated by theoretical models of dusty disk winds, where magnetic field lines drive dust above the disk plane close to the sublimation zone.

*Results. * Image reconstruction reveals an inclined disk with slight asymmetry along its minor-axis, likely due to inclination effects obscuring the inner disk rim through absorption of incident star light on the near-side and thermal re-emission/scattering of the far-side. Geometric modelling of a skewed ring finds the inner rim at $0.17\pm 0.02\leavevmode\nobreak\ \mathrm{au}$ with an inclination of $50.9\pm 1.0^{\circ}$ and minor axis position angle $60.8\pm 1.2^{\circ}$. Radiative transfer modelling shows a flared disk with an inner radius at $0.18$ au which implies a grain size of $0.4\,\mu$m assuming astronomical silicates and a scale height of $15.0$ au at $100$ au. Among the tested radiative transfer models, only the dusty disk wind successfully accounts for the K-band excess by introducing dust above the mid-plane.

1 Introduction

Protoplanetary disks are observed across all masses of young stellar objects (Lazareff et al. 2017; Kraus et al. 2017a). Created as a consequence of the conservation of angular momentum during the star formation process, infalling gas and dust form a circumstellar disk around protostellar objects. Circumstellar disks are understood to be the birthplace of planets. The scattering and thermal re-emission of starlight by optically thick dust in the inner disk induces stong near-infrared (NIR) excess emission.

Developments in the field of optical interferometry allowed for the first spatially resolved observations inner astronomical unit-scale of young stellar objects (YSOs) (Millan-Gabet et al. 1999; Akeson et al. 2000). It is within these inner regions that important processes, such as star-disk interactions (accretion and outflows), dust sublimation and planet formation occur. Accretion processes are one of the dominant effects in planet formation in the inner disk. Internal friction, or viscosity within the disk drives accretion onto the central star. In order to preserve angular momentum, some material is lost through outflow processes such as jets (Williams & Cieza 2011). Similarly, outflowing material such as magnetically driven winds can remove angular momentum from the disk, efficiently driving accretion processes (Turner et al. 2014). The current generation of interferometers provide unparalleled opportunities to study the smallest scales of disk structures (Davies et al. 2018; Setterholm et al. 2018), in particular the Centre for High Angular Resolution Astronomy (CHARA) can provide baselines up to $331$ m allowing us to study the very inner disk regions.

The first inner rim models of protoplanetary disks based upon spectral energy distribution (SED) analysis adopted a vertical inner wall of emission (Dullemond et al. 2001; Natta et al. 2001). However, interferometric observations failed to find evidence of the strong asymmetries predicted at high disk inclinations by this model (Monnier et al. 2006a; Kraus et al. 2009). An alternative curved rim model was later proposed, with a curvature arising from a gas-density-dependent sublimation temperature (Isella & Natta 2005). The brightness distribution predicted by these models is much more symmetric in nature. However, these models predict a very pronounced second visibility lobe, which has yet to be observed in other protoplanetary disks such as HD 142666 (Davies et al. 2018), AB Aur, V 1295 Aql and MWC 275 (Tannirkulam et al. 2008; Setterholm et al. 2018). Instead a much flatter second lobe is observed, indicating that a significant fraction of NIR flux originates outside of standard disk models. These rim-only models are as such insufficient to explain many interferometric observations.

Several physical mechanisms have been proposed to explain the observed emission outside of the dust sublimation rim. One prominent explanation is the existence of dusty disk winds which send material along magnetic fields near the dust sublimation zone. This allows for optically thick material to exist close enough to the central star to contribute to the NIR emission above the disk structure. This model has been shown to successfully account for the NIR excess of the SED and the basic visibility features of AB Aur, MWC 275 and RY Tau (Königl & Salmeron 2011; Petrov et al. 2019).

SU Aurigae (SU Aur) is a $5.18\pm 0.13\,\mathrm{Myr}$ old (Bochanski et al. 2018) low mass pre-main-sequence star in the Upper Sco star forming region at a distance of $158^{+1.49}_{-1.48}\,\mathrm{pc}$, obtained from GAIA DR2 parallax measurements (Luri et al. 2018). As a G2-type star it has a similar effective temperature to the sun (DeWarf et al. 2003), but a much higher bolometric luminosity at $12.06\,\mathrm{L_{\odot}}$ (calculated from GAIA DR2, Evans et al. 2018) putting it in the sub-giant star class. The stellar parameters adopted are listed in Table 1. SU Aur is known to be variable in the V band, varying up to $0.5$ mag over a several day cycle (Unruh et al. 2004). On the other hand variability in the K band is minimal (Akeson et al. 2005) allowing us to assume the flux contribution from the star is constant across different epochs of observations. As such, any variation in the visibility is likely geometric.

Spectroscopic and photometric monitoring of SU Aur by Petrov et al. (2019) has revealed that a dusty disk wind is the potential source of the photometric variability in both SU Aur and RY Tau. The characteristic time of change in the disk wind outflow velocity and the stellar brightness indicate that the obscuring dust is located close to the sublimation rim of the disk, in agreement with previous theoretical disk wind models (Bans & Königl 2012; Königl & Salmeron 2011).

Interferometric observations carried out by Akeson et al. (2005) using the Palomar Testbed Interferometer (PTI), a 3-telescope interferometer with baselines up to $110$ m, found a disk inclined at $62^{\circ+4}_{-8}$ with a minor axis position angle of $24\pm 23^{\circ}$ and a sublimation rim at $0.21$ au. The disk was modelled as a flared disk with a vertical inner wall of emission using radiative transfer to fit both interferometric and photometric data. However, optically thick gas close to the central star was needed to fit the SED. The disk geometry is in agreement with values from Eisner et al. (2014) using the Keck interferometer with values derived from an upper limit to the fit of the Br$\gamma$ emission. They find an upper limit on the disk inclination of $50^{\circ}$ and minor axis position angle of $50^{\circ}$. The difference in the derived position angles in previous studies is likely due to the lack of Br$\gamma$ emission and poor signal to noise of Eisner et al. (2014) and the limited baseline range of Akeson et al. (2005). Hence, new observations with significantly better uv coverage on much longer baselines were needed. Additionally, polarimetric imaging campaigns, undertaken by de Leon et al. (2015) and Jeffers et al. (2014) revealed the presence of tails protruding from the disk at both H-band and visible wavelengths, likely associated with an extended reflection nebula and a possible undetected brown dwarf companion encounter. A companion search was undertaken by the SEEDS (Strategic Exploration of Exoplanets and Disks with Subaru) imaging survey and ruled out the presence of a companion down to $10$ Jupiter masses at separations down to $15\,\mathrm{au}$, contradicting the potential brown dwarf encounter theory of de Leon et al. (2015).

This paper presents the lowest mass YSO to be studied with very long baseline ($>110$ m) NIR interferometry to date and is the first study probing the detailed rim structure of SU Aur with interferometric observation on baselines up to $331$ m. Three different modelling methodologies were applied: (i) Image reconstruction was used to obtain a model-independent representation of the data and to derive the basic object morphology. (ii) Following this geometric model fitting allowed us to gain an appreciation for the viewing geometry of the disk by fitting Gaussian and ring models to the data. (iii) Finally, we combine interferometry and photometry to derive physical parameters with radiative transfer analysis, where our particular focus is on the physical characteristics of the inner rim.

Our observations are described in section 2, with image reconstruction detailed in section 3. Our geometrical model fitting approach is described in section 4 and radiative transfer modelling and SED fitting are discussed in section 5. A discussion and analysis of the results can be found in section 6, followed by concluding remarks made in section 7.

2 Observations

The CHARA array is a Y-shaped interferometric facility that comprises six $1\,$m telescopes. It is located at the Mount Wilson Observatory, California, and offers operational baselines between $34$ and $331\,$m (ten Brummelaar et al. 2005). The CLIMB instrument, a three telescope beam combiner (ten Brummelaar et al. 2013), was used to obtain observations in the near-infrared K-band ($\lambda=2.13\,\mu m,\Delta\lambda=0.35\,\mu m$) between October 2010 and November 2014. We obtained 28 independent measurements of SU Aur, using seven different 3-telescope configurations with maximum physical baseline of $331\,$m corresponding to a resolution of $\lambda/(2B)=0.70\,\mathrm{mas}$ [milliarcseconds], where $\lambda$ is the observing wavelength and $B$ is the projected baseline. In addition, a small number of observations were taken in 2009 using the two-telescope CLASSIC beam combiner (ten Brummelaar et al. 2013), also at CHARA in the K-band along the longest ($331\,$m) projected baseline. Details of our observations, and the calibrator(s) observed for the target during each observing session, are summarised in Table 2. The uv plane coverage that we achieved for the target is displayed in Figure 1. Our data covers a relatively wide range of baseline lengths and position angles, making the data set suitable for image reconstruction.

The CLIMB and CLASSIC data was reduced using pipelines developed at the University of Michigan (Davies et al. 2018). This is much better suited to recovering faint fringes from low visibility data than the standard CHARA reduction pipeline of (ten Brummelaar et al. 2012). The measured visibilities and closure phases were calibrated using interferometric calibrator stars observed alongside the target. Their adopted uniform diameters (UDs) were obtained from JMMC SearchCal (Bonneau et al. 2006, 2011), when available, or gcWeb222http://nexsciweb.ipac.caltech.edu/gcWeb/gcWeb.jsp and are listed in Table 2.

During the data reduction it was found that our calibrators HD 31592 and HD 34052 exhibited strong closure phase signals, indicating the presence of close companions around these stars. In order to ensure these binary calibrators could be used to calibrate the primary science target, we cross-calibrated the data on these stars with other calibrators observed during the same nights to determine the binary parameters, as outlined in Section A in the appendix. Based on the fitted binary parameters, we could then correct the transfer function and use the data for the calibration of SU Aur.

In addition to the new observations from CHARA, archival interferometric data from other facilities was included for our analysis. A small amount of data was available from the Keck Interferometer (KI, Colavita et al. 2013; Eisner et al. 2014) from 2011 along a single $84$ m baseline, while a larger amount of data was also available from the Palomar Testbed Interferometer (PTI, Colavita et al. 1999) from 1999 to 2004 using a 2 telescope beam combiner on 3 different physical baselines between $86$ and $110$ m. This data was published in Akeson et al. (2005). These additional measurements complement the CHARA observations in the intermediate baseline range; the full uv coverage is shown in Figure 1.

Both the PTI and KI data were calibrated using the standard method outlined by Boden et al. (1998) using the wbCalib software available from NExcScI333http://http://nexsci.caltech.edu/software/V2calib/wbCalib/. The calibration pipeline works in conjunction with getCal and the Hipparcos catalogue (Perryman et al. 1997) for calibrator star astrometry and diameters. This was the same process used by Akeson et al. (2005) to extract visibilities from the PTI data. Our re-reduction of this data agrees with the results shown in the literature. The fully reduced data from all instruments is shown in Figure 2.

As we combine several years worth of data, care was taken to check for time dependencies in the visibilities of baselines of similar length and position angle. Variability in the K band is known to be minimal, so any time dependencies in the visibility amplitudes is likely geometric. However, no significant time dependencies were discovered.

3 Image Reconstruction

Image reconstruction techniques require broad and circular uv coverage along as many baseline lengths as possible. Fortunately, the data from the observations lends itself to this process as the uv plane has been well sampled, though some small gaps remain in the position angle coverage. By comparing visibilities from different instrument at similar baseline lengths and position angles, we can see there is likely very little extended emission in this system and so no correction for instrument field of view is required. This technique is useful for interpretation of non-zero closure phases, indicative of asymmetric distributions, in a model-independent way. Our closure phase values are shown in Figure 2. There are many different algorithms with which to reconstruct images from interferometric data, but the process described here involved the use of the Polychromatic Image Reconstruction Pipeline (PIRP) which encompasses the $MiRA$ reconstruction algorithm by Thiébaut (2008). The reconstruction procedure and results are described below.

In the $MiRA$ routine, the object is modelled as an unresolved central star with an extended, model-independent, environment (Kluska et al. 2014). Both components have different spectral behaviours and so differing spectral indices. Additionally, the type and weight of the regularisation was explored, $MiRA$ allows for either quadratic smoothing or total variation regularisations to be implemented. The regularisation plays the role of the missing information by promoting a certain type of morphology in the image. The quadratic smoothing algorithm aims for the smallest possible changes between pixels to produce a smoother image, it is particularly useful as it’s quadratic nature means it is less likely to find local minima. On the other hand, total variation aims to minimise the total flux gradient of the image and is useful to describe uniform areas with steep but localised changes. These regularisations are considered to be the best ones for optical interferometric image reconstruction; (Renard et al. 2011). The size and number of pixels also plays an important role in image reconstruction. One cannot simply use the maximum number of pixels of the smallest size to obtain better resolution, they have to be chosen to match uv plane sampling. It was found that a quadratic smoothing regularisation with a weight of $1\times 10^{9}$ and $526\times 526$ pixels of $0.1$ mas in size provides the best-fit image reconstruction when utilising exact Fourier transform methods. The optimal regularisation parameters were determined using the L-curve method. PIRP also allows for bootstrap iterations (Efron & Tibshirani 1994) starting from a previous best image. This involves a random draw of data points within the dataset to determine the reliable features of the image. The process was carried out 500 times allowing for a pixel by pixel error estimation, see Kluska et al. (2016).

The final image is shown in Figure 3 (upper left panel), which also shows the asymmetry in the intensity map (upper right panel). The contours represent the $1\sigma$ (dashed line), $3\sigma$ and $5\sigma$ (solid lines) uncertainties. The asymmetry is calculated by rotating the image through $180^{\circ}$ and subtracting it from the un-rotated image. It is this residual flux that produces the non-zero closure phases, highlighting any areas of greater emission within the disk. Using this technique we can see that the disk has greater intensity in the Eastern regions, which contains $8\%$ more flux that the Western regions with a flux ratio of $1.07$, where the remaining flux is in the central star. By inclining a flared disk with the Western region (bottom right of image) towards the observer the nearside of the inner rim becomes self-shadowed by the near-side disk rim, so the Eastern region inclined away from the observer appears brighter. The best-fit image shows a disk radius of $1.0$ mas, a minor-axis position angle of $41^{\circ}\pm 3$ and an inclination of $51^{\circ}\pm 5$. The fit of the image to the visibility and closure phases in the data is shown in Figure 3, with a combined visibility and closure phase reduced $\mathrm{\chi^{2}_{red}}$ of $4.57$.

4 Geometric Model Fitting

The next step in interpreting our interferometric observations is the fitting of simple geometric models to the observed quantities. The visibility profile (Figure 2) reveals a clear drop in visibility through short and intermediate baselines with a possible plateau/second lobe at the longest baselines. In the visibility profile we do not see any evidence for structures on distinct different spatial scales that might indicate the presence of a binary companion or extended halo emission. The closure phases are all $<20^{\circ}$, indicating weak asymmetric features as evidenced in the image reconstruction above. To explore the viewing geometry of the disk, a series of simple intensity distributions were tested against our data. All the models tested contained an unresolved point source that was used to represent the central star, a reasonable assumption given the expected angular diameter of the star (see Table. 1). The disk component was modelled as one of four intensity distributions: (i) A Gaussian to simulate a disk with unresolved inner rim with a FWHM free parameter. (ii) A ring to simulate emission from a bright inner rim only with a defined fractional width equal to $20\%$ of the radius with an inner radius ($R_{min}$) free parameter. (iii) A skewed ring model with a diffuse radial profile defined by an inner radius ($R_{min}$), a FWHM ($w$) and an azimuthal brightness modulation, where $c_{j}$ and $s_{j}$ are the cosine and sine amplitudes for mode $j$, in an attempt to model disk asymmetries. Lazareff et al. (2017) find that the diffuse skewed ring here can be used to successfully model a wide range of YSOs.

With the exception of the skewed ring, these models are intrinsically axisymmetric, but we project the brightness distribution in order to mimic viewing geometry effects that are parameterised with an inclination angle $i$ (defined with $0^{\circ}$ as face-on) and a disk position angle $\theta$. We measure disk position angles along the minor axis and follow the convention that position angles are measured from North ($\theta=0^{\circ}$) towards East.

The stellar-to-total flux ratio can be calculated by comparing K-band photometry (Cutri et al. 2003b) with the stellar atmosphere models of Castelli & Kurucz (2004), which gives a ratio of $1.17$ (using the stellar parameters listed in Table 1). However, it cannot be fixed as this has been shown to introduce unreasonable large scale components when fitting long baseline data (Ajay & Venkataraman 2013). As such the parameter space of the flux ratio was explored step-wise for all models with a range of model parameters. In this way the stellar-to-total flux ratio was constrained for all models. It was found that a ratio of $1.13\pm 0.01$ provided the best fit to the data for the ring, skewed ring and TGM models. After determining the star-to-disk flux ratio, the parameter space of the geometric models could be fitted to the observed visibilities and closure phases for each of the observed baselines given initial parameter constraints based on the literature values of Akeson et al. (2005) described in Section 1. We used a bootstrap method to explore the parameter space around these initial values and to compute uncertainties on the individual parameters by fitting Gaussian distributions to parameter histograms.

The best-fit parameters and associated errors for each of the geometric models are shown in Table 3. All test models agree with respect to the position angle and inclination of the disk very well (Table 3). However, none of the models provide good fits to the data as evidence by the $\chi^{2}_{red}$ values of between 11.86 (ring model) and $8.57$ (skewed ring). The skewed ring model found a minor axis position angle of $61.0^{\circ}\pm 1.0$ and inclination of $51.2^{\circ}\pm 1.1$. Of the individual models, the skewed ring provides the best fit to both the visibilities and closure phases. Both the Gaussian model and the Skewed Ring are found to be quite diffuse with a FWHM of $2.01\leavevmode\nobreak\ \mathrm{mas}\pm 0.02$ and width of $0.96\leavevmode\nobreak\ \mathrm{mas}\pm 0.02$ respectively.

Overall, we achieved the best-fit with a skewed ring with an inner radius of $1.10\leavevmode\nobreak\ \mathrm{mas}\pm 0.12$ and a ring FWHM of $0.96\leavevmode\nobreak\ \mathrm{mas}\pm 0.02$, thus making the ring very diffuse with only a marginally defined inner radius. There was no evidence discovered for any over-resolved extended emission or ’halo’ found in many other objects (Monnier et al. 2006b; Kraus et al. 2009).

In order to model the observed CP signal, with a maximum of $20\pm 12^{\circ}$, which may indicate a slight asymmetry in the disk, we introduced asymmetries to the skewed ring model (Lazareff et al. 2017). However, this improved the fit only marginally over the standard ring model from Table 3. It was found that zero closure phases produced a reduced chi-squared ($\mathrm{\chi^{2}_{red,CP}}$) fit of $1.90$, while a sinusoidally modulated asymmetric ring only resulted in a $\mathrm{\chi^{2}_{red,CP}}$ of $1.66$. The model visibility curves corresponding to the best-fit skewed ring model are shown in Figure 4, 5 and 6 (red curve). The different panels show visibilities towards different position angle bins.

Of the simple geometric models tested, the skewed ring model provides the best fit. However, non can be said to provide a good fit to the observed data, as evidence by the $\chi^{2}_{red}$ values shown in Table 3. As such, more complex disk structures are required, such as flared disks, different rim morphologies or disk winds. In the next section, these possibilities are explored in detail using radiative transfer techniques. This allows us to not only explore complex geometries, but to derive physical parameters such as radial density profiles and the disk scale heights.

5 Radiative Transfer Modelling

We used the TORUS Monte-Carlo radiative transfer code (Harries 2000), allowing for the simultaneous fitting of visibility, closure phase and photometric data to further constrain the geometry and physical dust properties of the SU Aurigae circumstellar disk.

Starting from the disk properties derived by Akeson et al. (2005, Table 4), we explored radiative transfer models with different scale heights (where scale height is that of the gas parameterised at $100\leavevmode\nobreak\ \mathrm{au}$ with flaring index $\beta$) and inner rim shapes. In our TORUS simulations, the dust was allowed to vertically settle and the dust sublimation radius was left as a free parameter, allowing the inner rim radius to define itself based on well-defined rules of the Lucy (1999) iterative method to determine the location and the temperature structure of the whole disk. This is implemented whereby the temperature is initially calculated for grid cells in an optically thin disk structure, with dust added iteratively to each cell with a temperature lower than that of sublimation, until the appropriate dust to gas ratio is reached ($0.01$). Once TORUS has converged to radiative equilibrium a separate Monte Carlo algorithm is used to compute images and SEDs based on the optical properties of the dust species implemented. We confirmed that stellar photosphere models of Castelli & Kurucz (2004) using these stellar parameters can reproduce the photometry measurements of SU Aur reasonably well across the visible continuum. The grain size distribution used by Akeson et al. (2005) is a distribution of astronomical silicate grains up to $1\,\mathrm{mm}$ in size. We adopt a silicate grain species with dust properties and opacities adopted from Draine & Lee (1984). The initial density structure of the gas is based upon the $\alpha$-disk prescription of Shakura & Sunyaev (1973) where the disk density is given as:

[TABLE]

Here, $z$ is the vertical distance from the midplane while the parameters $h(r)$ and $\Sigma(r)$ describe the scale height and the surface density respectively.

In our radiative transfer models, we represent the stellar photosphere with a Kurucz (1979) model atmosphere using the stellar parameters outlined in Table 1. The photometric data was obtained from a wide range of instruments from the VizieR database and are compiled in Table 7. Where multiple observations in the same waveband were present care was taken to minimise the total number of instruments and keep the number of observation epochs as close as possible to minimise any potential variability effects.

Visibilities were calculated from synthetic images of the disk system (as shown in Figure 7) extracted through application of the van Cittert-Zernicke theorem, applied using a 1D Fourier transforms projected onto the observed baseline position angle. Phases are also extracted from the images and are used to calculate the closure phase, as described in Davies et al. (2018).

The parameter space of the radiative transfer models was explored objectively using the values of Akeson et al. (2005) as a starting point. A range of physically realistic values for each parameter was explored in a broad grid of models (as described in Table 4). A $\chi^{2}$ value was then computed for the visibilties, closure phases and SED fits of each model allowing the grid to be refined around the minimum. The interferometric and photometric data points were fitted simultaneously, with the resulting $\chi^{2}$ values shown in Table 5. The silicate feature at $10\,\mathrm{\mu m}$ allows us to place some constraints on the dust sizes, as larger grains produce smaller features. The growth of dust grains and their effects on observed silicate features is described in a review by Natta et al. (2007). The IR flux is controlled by the morphology of the sublimation rim. As the inner radius increases, the amount of circumstellar material emitting in mid-IR wavelengths is reduced, and the IR emission decreases. The shape of the mid-IR excess also describes the degree of flaring present in the disk, where large excess indicates greater flaring. A larger flaring power in a disk leads to an increasing surface intercepting the starlight, and therefore an increase in reprocessed radiation. We adjusted the total dust mass in the model in order to match the millimeter flux. A detailed description of the effect of disk parameters on the SEDs of protoplanetary disks can be found in Robitaille et al. (2007).

Due to the optical depth of the system the inner-rim of the disk appears as the brightest part of the disk at NIR wavelengths. This is because the rest of the disk is shadowed by the inner rim and only rises out of shadow in cooler regions of longer wavelength emission.

5.1 Sublimation Rim Model

The curved rim of Isella & Natta (2005), henceforth IN05, is based upon a single grain size prescription with a gas density-dependent sublimation temperature. Due to a vertical gas pressure gradient, the sublimation temperature decreases away from the mid-plane, creating a curved rim. The sublimation temperature of the grains follows

[TABLE]

where the constant $G=2000$ K, $\gamma=1.95\times 10^{-2}$, $r$ is the radial distance into the disk and $z$ is the height above the midplane (Pollack et al. 1994).

A curved rim is shown to be a viable disk model by Flock & Turner (2016) and Flock et al. (2016), based upon extensive hydrodynamical simulations. A wide range of disk structure parameters are explored to find the best fit solution to both the interferometric measurements and the SED. Figure 4 shows the results of radiative transfer modelling of the IN05 rim prescription. The model SED shows a clear flux deficit in the NIR, with a K-band flux of just 64% of the 2MASS photometric point, far outside the limited range of variability of SU Aur. This is also clear in the visibility curves, shown in green, where the overall shape is a good fit, but the minimum visibility is too high due to a larger than expected stellar contribution.

An alternative sublimation front geometry is proposed by Tannirkulam et al. (2007), henceforth THM07. This model employs a two-grain scenario, were a mixture of small $0.1\,\mathrm{\mu m}$ grains and large $1.2\,\mathrm{\mu m}$ grains has been adopted, with the mass of larger grains fixed at $9$ times the mass of smaller grains. The smaller grains are not allowed to settle, so the scale height is fixed to that of the gas. The larger grains are allowed to settle to 60% of the scale height of the gas. This combined with the larger grains existing closer to the star, due to more efficient cooling, leads to an elongated and curved sublimation front. Figure 5 shows the results of this modelling. The SED again shows a clear deficit in NIR flux, with a K-band flux of just 68%, comparable with that of the IN05 prescription. The presence of dust closer to star changes the shape of the visibility curve dramatically, with the first lobe now extending to much longer baselines. The two grain THM07 model proves to be a worse fit than the single grain IN05 model.

Our results show that the best fit of the curved-rim disk model of IN05 can be achieved with a single silicate grain species with a single grain size of $0.1\,\mu m$. The addition of larger grain species further reduced near-infrared flux in 1-3 $\mu m$ region, resulting in a poorer fit to the SED in both the shape and magnitude of the NIR excess. The disk is also found to be highly flared and extending from $0.15$ to $100$ au, loosely constrained by the long wavelength photometry, while still solving for vertical hydrostatic equilibrium. Importantly, comparison of the stellar atmosphere, represented by the Kurucz model atmosphere (Castelli & Kurucz 2004), with the photometric data can provide the stellar-to-total flux ratio for each waveband. In the K-band this ratio is found to be $1.27$, meaning the circumstellar environment contributes $44\%$ of the total flux. The dust to gas ratio is fixed to $0.01$, as taken from literature values (Akeson et al. 2005). The NIR flux deficit results also in a poor fit to the K-band visibilities.

5.2 Dusty Disk Wind Model

In an attempt to increase the NIR flux contributions in the model, we explored a dusty disk wind scenario as set out by Bans & Königl (2012), henceforth BK12. This mechanism is based on the presence of a large-scale, ordered magnetic field which threads the disk. The field could originate from the interstellar field that permeates the molecular cloud and is dragged by in-falling gas into the disk. The magnetic field strength required to drive these winds is of the order of kGauss. Fields of this strength have been shown to be present in other T Tauri stars ($2.35\pm 0.15\,\mathrm{kG}$ for T Tau) by Guenther et al. (1999) and Johns-Krull et al. (1999). In this model material is flung out along magnetic field lines highly inclined to the disk surface. The high magnetic pressure gradient above the disk surface accelerates the material which is then collimated through the azimuthal and poloidal field components (Bans & Königl 2012). These centrifugally driven winds are highly efficient at distributing density above and below the plane of the disk, carrying angular momentum away from the disk surface.

The BK12 wind is launched from the disk surface in a region between an inner and outer radius, located at a maximum of twice the sublimation radius. It is FUV radiation that plays a key role (through photoevaporation) in the mass loading of the wind at the inner launching region. The gas density structure is given as

[TABLE]

where $\rho$ is the density along a flow line at a distance $r_{0}$ from the symmetry axis, i.e. along this disk surface. $\rho_{1}$ denotes the density at the fiducial radius of $r_{1}$, which is given as $1$ au while $\chi$ is related to the cylindrical coordinate $z$ through

[TABLE]

$\eta(\chi)$ is obtained from the solution of the MHD wind equations. The into-wind mass transfer is linked to the density through

[TABLE]

where $\dot{M}_{out}$ is the mass outflow rate, $M_{*}$ is the stellar mass, $h_{0}$ is the disk scale height and $\psi_{0}$ is the ratio of vertical speed to Keplarian speed at the disks surface. $\xi_{0}^{{}^{\prime}}$ is related to the angle $(\theta_{0})$ at which the poloidal component of the magnetic field threads the disk, defining the opening angle of the disk wind (Safier 1993)

[TABLE]

This rate generally controls the magnitude of the NIR excess added by the dusty wind.

This prescription is taken from Safier (1993), and assumes a steady, axisymmetric, effectively cold disk outflow. For the dust distribution, a constant dust-to-gas ratio is assumed to match that of the disk. The wind is populated with dust in the same way as the disk and also converges towards radiative equilibrium with each Lucy (1999) iteration. Full details of the disk wind model implemented can be found in Bans & Königl (2012) and Königl & Salmeron (2011). A wide parameter search was undertaken, in order to determine the optimum solution (see Table 4). The stellar parameters were fixed to those shown in Table 1 and the same silicate prescription was adopted for all models, as in the IN05 and THM07 prescriptions The position angle of the disk was fixed to that of the best fit geometric model listed in Table 3.

The results of our radiative transfer modelling of the BK12 wind are shown in Figure 6. The cross section shows the very inner-rim of disk, uplifted dust above and below the mid-plane can clearly be seen. The inner-rim is also curved using the IN05 prescription, although a grain size of $0.4\,\mathrm{\mu m}$ is required to produced the observed excess across the infrared. A range of grain sizes were tested from $0.1$ to $1.5\,\mathrm{\mu m}$ but $0.4\,\mathrm{\mu m}$ grains provided the best fit to the K-band photometric flux values.

An into-wind mass outflow rate of $1\times 10^{-7}\mathrm{M_{\odot}yr^{-1}}$ was required to uplift enough material to reproduce the observed excess. Lower into-wind outflow rates do not allow enough material to exist exterior to the inner-rim to reproduce the observed K-band excess in the photometry and interferometric data. If one assumes an outflow to accretion ratio of $0.1$ one can predict an accretion rate of $1\times 10^{-6}\mathrm{M_{\odot}yr^{-1}}$. This is unphysically high for an object such as SU Aur given the expected age of the system ($5.18\pm 0.13\leavevmode\nobreak\ \mathrm{Myr}$), assuming that the rate stayed constant over full period. In addition, this is in disagreement with non-detection of $\mathrm{Br\gamma}$ emission by Eisner et al. (2014) which is suggestive of a lower outflow/accretion rate. Other free parameters in this parameterisation of a dusty disk wind cannot reproduce the effect of a high into-wind mass accretion rate.

The dust wind scenario has the effect of increasing the amount of the dust close to the star, where temperatures are sufficient for NIR flux contribution. As shown in Figure 6 this model is in agreement with the observed NIR photometry, with a K-band flux of 102% of the photometric value, well within the range of variability of SU Aur. The rest of the SED is still well fitted, as the curvature and shape are based upon the IN05 prescription described above. The visibility curves shown provide a good fit to the data points, matching the lowest visibilites well. However, one cannot negate the potentially unphysically high outflow/accretion rates required to successfully model the data. The introduction of disk wind also had the interesting effect of flattening the ’bump’ in the second visibility lobe present in many disk models, including the IN05 prescription. The $\chi^{2}$ fits for the visibilities and SED for all three models are shown in Table 5.

6 Discussion

In investigating the circumstellar environment of SU Aur we have explored the structure and composition of the disk and greatly improved the constraints on the parameters initially taken from literature. The wide variety of techniques used to analyse the interferometric data allow us to precisely define the disk characteristics.

Image reconstruction shows a disk inclined at $52.8^{\circ}\pm 2.2$, this is in agreement with values of $63^{\circ+4}_{-8}$, ~$60^{\circ}$ and ~$50^{\circ}$ found by Akeson et al. (2005); Unruh et al. (2004); Jeffers et al. (2014) respectively. The minor axis $\theta$ on the other hand was found to be $50.1^{\circ}\pm 0.2$, greater than the literature values of $24^{\circ}\pm 23$ and $15^{\circ}\pm 5$ found by Akeson et al. (2005); Jeffers et al. (2014). This difference is likely due to either: The poor uv coverage and lack of longer baselines in previous interferometric studies, both of which make estimating the position angle and inclination particularly unreliable. Other non-interferometric studies focus on the outer disk, rather than the inner au-scale regions. The image reconstruction also reveals evidence of slight asymmetries within the disk at $3$ sigma significance level, that are consistent with an inner disk rim seen at an intermediate inclination. As the K-band emission primarily traces the very inner region of the rim; if a disk is inclined the near side of the rim will be partially obscured from view, while the far side of the rim will be exposed to observation. This explanation can successfully account for the over-brightness observed in the Eastern disk region in the image reconstruction shown in Fig. 3 and can also be seen in the radiative transfer image shown in Fig. 3. Models of the effect of inclined disk on the observed brightness distribution are described by Jang-Condell & Turner (2013). There are several other scenarios that have been used to explain asymmetries in protoplanetary disks in the past. Two possible scenarios are: Firstly, a shell ejection episode that can carry dust and gas away from the central star can be capable of reproducing the photometric variability in different epochs of observations (Borges Fernandes et al. 2009; Kluska et al. 2018). Also, the presence of a companion embedded within the disk can create dust trapping vortices that capture dust grains. These vortices, however, are known to trap primarily large grains (mm-size), not the small micron-sized grain we observe in the infrared, as shown by Kraus et al. (2017b); van der Marel et al. (2013). We rule out the presence of a companion by undertaking a companion search using our geometric models; a second point source was iterated through the parameter space with a grid size of $100$ mas in steps of $0.1$ mas (see Sec. 3). However, the model fit did not improve significantly by adding an off-center point source. We therefore favour the explanation of asymmetry arising from an inclined disk.

Our geometric model fits were key in understanding the circumstellar environment of SU Aur. It was found that a skewed ring structure is able to fit our data best, which suggests that we trace a diffuse inner disk edge. This finding is consistent with the studies on other YSOs that found bright inner rims, such as the Lazareff et al. (2017) survey of 51 Herbig AeBe stars using the PIONIER instrument at the VLTI. They found that over half of the disks could be successfully modelled using a diffuse ring structure. A high optical depth of the circumstellar material causes a bright inner rim, where most of the radiation is absorbed, scattered or re-emitted. Our best-fit Skewed ring model suggests that the stellar-to-total flux ratio is $1.13$ and achieves a $\mathrm{\chi^{2}_{red}}$ of $8.57$. Akeson et al. (2005) find that $44\pm 9\%$ of the total flux is in the SED K-band excess, with a $4\%$ of flux in an extended envelope, this is in good agreement with the values found from geometric modelling. The skewed ring model fits a ring of radius $0.17\pm 0.02$ au at an inclination of $51.2^{\pm}1.1^{\circ}$. These values are in remarkable agreement with both the values derived from image reconstruction in this study and the literature values of Akeson et al. (2005) of $0.18\pm 0.04$ au and $62^{\circ+4}_{-8}$, where the slight differences are likely due to difference between the diffuse profile, skewed ring and the standard ring structures employed. The minor axis $\theta$ of $61.0\pm 1.0^{\circ}$ is similar to the image reconstruction value, which is significantly larger than literature values. As above, this is most likely due to the poor uv coverage and short baselines available in previous interferometric studies, making the estimates of position angle and inclination particularly unreliable. The skewed ring also introduces modulated asymmetries into the ring profile. The skewed ring geometric model with azimuthal brightness modulation results in an improved fit with $\mathrm{\chi^{2}_{red}}=1.66$, where the contrast of the asymmetry is consistent with the one found using image reconstruction techniques. Similarly this can be attributed to inclination effects of the viewing geometry (Jang-Condell & Turner 2013). However, the position angle of the asymmetry is not well constrained in these models.

Radiative transfer modelling of the disk allowed us to fit a physical disk model to the visibility and photometry data simultaneously, meaning the 3-D density distribution of the disk can be explored. In this paper, three different geometries are considered: The single grain curved rim of IN05, the two grain curved rim of THM07 and the addition of a dusty disk wind of BK12 to the single grain curved rim. In the case of the IN05 prescription, we follow the idea that the sublimation temperature is gas-pressure dependent allowing the rim shape to be defined as described in Sec. 5.1. This model provides good constraints on the characteristic size of the near-infrared emitting region and the flaring in the colder regions, with a sublimation temperature of $1600$ K corresponding to an inner radius of $0.12$ au, slightly smaller than than the literature values of $0.18\pm 0.04$ au (Akeson et al. 2005) and $0.17\pm 0.08$ au (Jeffers et al. 2014). The disk was also found to be at an inclination of ~$50^{\circ}$ and position angle of ~$45^{\circ}$, in agreement with both the literature and above mentioned methods. This is shown to fit the photometry well at both shorter and longer wavelengths. However, there is a clear deficit in the IR excess, which also leads to poorly fitted visibilities due to an over-estimation of the stellar-to-total flux ratio. The same issue is obvious in the THM07 prescription, where larger grains are introduced allowing dust to exist closer to the star, though this model also fails to reproduce the visibility curve of the geometric modelling, as the inner radius of the disk is much smaller. A two grain model does not provide a good fit to these observations.

The dusty disk wind prescription of BK12 was incorporated into the single grain disk model of IN05. Dust flung out from the inner regions of the disk is carried far above and below the mid-plane. This dust is directly exposed to stellar radiation so is hot enough to contribute to NIR flux, whilst also obscuring the direct stellar flux. This was shown to be a physically viable scenario for YSOs, including SU Aur, by Königl & Salmeron (2011) and Petrov et al. (2019). The resulting model reveals a disk with an inner radius of $0.15$ au and a dust-to-gas ratio of $0.01$. The inner radius is in agreement with the literature values discussed above. The flaring parameters $\alpha_{disk}$ and $\beta_{disk}$ were fixed such that $\alpha_{disk}=\beta_{disk}+1$ and found to be $2.4$ and $1.4$ respectively. The dusty disk wind mechanism can directly reproduce the flux ratio in the K-band, allowing for an good visibility fit and an improved SED fit. As the disk wind rises above the mid-plane it also shields the cooler parts of the disk, reducing the longer wavelength flux compared to rim-only models. This was compensated for in our models by increasing the scale height of the disk to $15$ au at a radius of $100$ au. However, the implementation of the dusty disk wind in this scenario is not completely physical, owing to the high into-wind outflow rate of $1\times 10^{-7}\mathrm{M_{\odot}yr^{-1}}$ required. The BK12 model also had the effect of flattening the second visibility lobe, a feature found in other YSOs (Tannirkulam et al. 2008; Setterholm et al. 2018) and potentially opens powerful future modelling pathways for these objects. The grain size of the silicate dust species that produced the best fit was found to be $0.4\,\mathrm{\mu m}$ with no evidence of larger grains, as this addition resulted in a worse fit to the shape and magnitude of the NIR excess, particularly the shape of the Silicate feature around $10$ um. This is in contrast to other inner disk studies where larger $1.2$ um grains are required (Kraus et al. 2009; Davies et al. 2018).

All the inner rim models investigated differ from the model proposed by Akeson et al. (2005) whereby a vertical inner wall was combined with a small optically thick inner gas disk very close to the star aligned with the outer disk. This optically thick gas was implemented through very simple black-body emission models and allowed the author to successfully reproduced the observed NIR bump in excess flux. While TORUS could implement this type of black-body emission, it is unable to simulate the gas emission in a self-consistent physical manner.

7 Conclusions

This interferometric study of SU Aurigae has revealed the complex geometry and composition of the disk around SU Aurigae. We summarise our conclusions as follows:

• •

We reconstruct an interferometric image that confirms the inclined disk described in literature. We see evidence for an asymmetry in the brightness distribution that can be explained by the exposure of the inner-rim on the far side of the disk and its obscuration on the near side due to inclination effects. Our data set does not permit the imaging fidelity that would be needed to detect evidence of ongoing planetary formation within the inner disk, such as small-scale asymmetries, gaps or rings.

• •

We see no evidence for a companion, in either the reconstructed images nor in the geometric model fitting procedures.

• •

Our simple geometric model fits reveal a disk of inclination $51.2\pm 1.2^{\circ}$ along a minor axis position angle of $61.0\pm 1.0^{\circ}$ and an inner radius of $1.12\pm 0.12$ mas ($=0.17\pm 0.02$ au). The disk is best modelled with a skewed ring which has a Gaussian ring width profile and sinusoidally modulated asymmetry. However, the poor $\chi^{2}_{red}$ of this model fit means the uncertainties quoted here are likely not representative of the the true range of values.

• •

Radiative transfer modelling shows that simple curved rim disk geometries of IN05 and THM07 cannot effectively model both the SED and visibility data. A deficit of NIR flux is obvious in the failure to reproduce our K-band observations.

• •

A dusty disk wind scenario can successfully account for both the observed excess in the SED and the observed visibilities. The dusty disk wind scenario described here lifts material above the disk photosphere, thus exposing more dust grains to the higher temperatures close to the star responsible for the NIR excess. However, the high accretion rate required to reproduce the stellar-to-total flux ratio may make this scenario physically invalid.

• •

Our best-fit model (dusty disk wind model) suggests that the dust composition in the disk is dominated by medium sized grains ($0.4\,\mathrm{\mu m}$) with a sublimation temperature of $1600$ K. Introducing larger grains results in a worse fit to the SED shape and NIR excess. The disk is also shown to be highly flared ($15$ au at $100$ au).

• •

The dusty disk wind model predicts a rather flat visibility profile at long baselines. This class of models avoids the pronounced visibility ’bounce’ that are associated with sharp edges in brightness distributions, as predicted by rim-only models. Therefore, these models may also open a pathway to physically model other YSOs that have been observed with $\gtrsim\,300$ m infrared long-baseline interferometry, such as AB Aur, MWC 275, and V1295 Aql (Tannirkulam et al. 2008; Setterholm et al. 2018) which all observe very flat long baseline visibility profiles.

Acknowledgements.

We acknowledge support from an STFC studentship (No. 630008203) and an European Research Council Starting Grant (Grant Agreement No. 639889). This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in Ochsenbein et al. (2000).

Appendix A Appendix A.

A.1 Binary Fit Procedure

The two calibrators HD 31952 and HD 34053 were found to be binary systems, based on strong non-zero closure phase signals. As both of these stars were needed for the calibration for 3 nights of our data, a binary fit was undertaken in order to recalculate the transfer function with which our data was calibrated. The software package LITpro was used to construct the fits within a search radius of $10$ mas. The parameters fitted are: angular separation, position angle and the uniform disk diameters (UDD) of the primary and secondary components. Two uniform disks were used to represent the stars and $\chi^{2}$ maps were constructed to find the best-fit location of the secondary star. The parameters are displayed in Table 6.

Appendix B Appendix B.

Table of photometry used in the SED fitting procedure.