© ESO 2019

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 enabled 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 ofemission (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. 2006; 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 ± 0.13 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.48}^{+1.49}\text{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 L_⊙ (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 three-telescope interferometer with baselines up to 110 m, found a disk inclined at ${62}_{-8}^{°+4}$ with a minor axis position angle of 24 ± 23° 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γ emission. They find an upper limit on the disk inclination of 50° and minor axis position angle of 50°. The difference in the derived position angles in previous studies is likely due to the lack of Brγ emission and poor signal to noise of Eisner et al. (2014) and the limited baseline range of Akeson et al. (2005). Hence, new observationswith 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 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 Sect. 2, with image reconstruction detailed in Sect. 3. Our geometrical model fitting approach is described in Sect. 4 and radiative transfer modelling and SED fitting are discussed in Sect. 5. A discussion and analysis of the results can be found in Sect. 6, followed by concluding remarks made in Sect. 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 (λ = 2.13 μm, Δλ = 0.35 μm) between October 2010 and November 2014. We obtained 28 independent measurements of SU Aur, using seven different two-telescope configurations with maximum physical baseline of 331 m corresponding to a resolution of λ∕(2B) = 0.70 mas (milliarcseconds), where λ 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 Fig. 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 gcWeb¹ 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 Appendix A. 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 two-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 Fig. 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 NExcScI². 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 Fig. 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.

Fig. 1 Coverage of the uv plane of the interferometric observations. Blue points represent observations using the CLIMB instrument, while red stars represent the CLASSIC instrument at the CHARA array. Yellow triangles represent observations using the Keck interferometer in 2011 and green triangles using the PTI with dates between 1999 and 2004.

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 Fig. 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 × 10⁹ and 526 × 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 Fig. 3 (upper left panel), which also shows the asymmetry in the intensity map (upper right panel). The contours represent the 1σ (dashed line), 3σ and 5σ (solid lines) uncertainties. The asymmetry is calculated by rotating the image through 180° 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° ± 3 and an inclination of 51° ± 5. The fit of the image to the visibility and closure phases in the data is shown in Fig. 3, with a combined visibility and closure phase reduced ${\chi }_{\text{red}}^2$ 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 (Fig. 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°, 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° as face-on) and a disk position angle θ. We measure disk position angles along the minor axis and follow the convention that position angles are measured from north (θ = 0°) towards east.

The stellar-to-total flux ratio can be calculated by comparing K-band photometry (Cutri et al. 2003) with the stellar atmosphere models of Castelli & Kurucz (2003), 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 ± 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 inSect. 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 }_{\text{red}}^2$ 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° ± 1.0 and inclination of 51.2° ± 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 mas ± 0.02 and width of 0.96 mas ± 0.02 respectively.

Overall, we achieved the best-fit with a skewed ring with an inner radius of 1.10 mas ± 0.12 and a ring FWHM of 0.96 mas ± 0.02, thus makingthe 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. 2006; Kraus et al. 2009).

In order to model the observed CP signal, with a maximum of 20 ± 12°, 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 (${\chi }_{\text{red},\text{CP}}^2$) fit of 1.90, while a sinusoidally modulated asymmetric ring only resulted in a ${\chi }_{\text{red},\text{CP}}^2$ of 1.66. The model visibility curves corresponding to the best-fit skewed ring model are shown in Figs. 4–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, none can be said to provide a good fit to the observed data, as evidence by the ${\chi }_{\text{red}}^2$ 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.

Fig. 2 Squared visibility and closure phase measurements against the de-projected baseline length of our interferometric data. Blue data points are those from the CLIMB/CHARA instrument. Red data points are from the CLASSIC/CHARA instrument, green data is that obtained from the PTI instrument and orange is that from the KI instrument.

5 Radiative transfer modelling

We used theTORUS 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 au with flaring index β) 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 (2003) 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 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 α-disk prescription of Shakura & Sunyaev (1973) where the disk density is given as: $$\rho (r,z)=\frac{\Sigma (r)}{h(r)\sqrt{2\pi }}\exp \{-\frac{1}{2}{[\frac{z}{h(r)}]}^2\}.$$(1)

Here, z is the vertical distance from the midplane while the parameters h(r) and Σ(r) describe thescale height and the surface density respectively. In our radiative transfer models, we represent the stellar photospherewith a Kurucz (1979) model atmosphere using the stellar parameters outlined in Table 1. The photometricdata was obtained from a wide range of instruments from the VizieR database and are compiled in Table B.1. 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 Fig. 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 χ² 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 χ² values shown in Table 5. The silicate feature at 10 μ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.

Fig. 3 Top left: resultant reconstructed image showing the 3 and 5σ significance levels as solid black lines and the 1σ level as dashed lines. The beam size is shown in the bottom right. The colour bar is the same for both maps, with the maximum intensity normalised to 1 for the sake of readability. Top right: extracted asymmetry in the intensity map shown with 3σ significance level. Bottom left: fitted visibilities, data shown in black with model visibility in blue. Bottom right: fitted closure phases, data shown in black with model CP in blue. The residuals normalised to the standard deviation are plotted in the bottom of each graph.

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 $$T_{\text{sub}}=G{\rho }^{\gamma }(r,z),$$(2)

where the constant G = 2000 K, γ = 1.95 × 10⁻², 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 μm grains and large 1.2 μ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 μm. The addition of larger grain species further reduced near-infrared flux in 1–3 μ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 2003), 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 environmentcontributes 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.

Fig. 4 Results of radiative transfer and geometric modelling with the IN05 disk rim prescription (Isella & Natta 2005) of a gas density-dependent sublimation temperature with a single grain species size 0.1 μm . Top left: disk density cross section of the inner rim. A logarithmic colour scale is used with a minimum density of 1.00 × 10−30 to a maximum of 4.64 × 10−10. Top right:SED computed with our radiative transfer model. Dark blue curve is the simulated blackbody emission of the central star. Green points are photometric observations while the short red line is the Spitzer spectrum. The coloured curves represent the SED at the different inclinations of 50, 57, 64 and 70°. Bottom: visibility data binned by position angle of observation. The data points are split by instrument consistently with Figs. 1 and 2, where blue circles are from CHARA/CLIMB, red stars are from CHARA/CLASSIC, green triangles from PTI and orange diamonds from KI. The red curves are the results of the best fit geometricskewed ring model and the green curves are calculated from the radiative transfer image at an inclination of 50° and a position angle of 60°. The dashed bounding lines indicate the minimum and maximum model visibilities for that position angle bin. The very bottom right panel shows the observed CP measurements (black) and the CP computed from the radiative transfer image (green).

Fig. 5 Results of radiative transfer and geometric modelling with the THM07 disk rim prescription (Tannirkulam et al. 2007) of a gas density-dependent sublimation temperature with two grain species, a majority larger grains at 1.2 μm with fewer smaller grains at 0.1μm. The mass of larger grains is fixed at 9 times the mass of smaller grains. Top left: disk density cross section of the inner rim. A logarithmic colour scale is used with a minimum density of 1.00 × 10−30 to a maximum of 1.23 × 10−09. Top right:SED computed with our radiative transfer model. Dark blue curve is the simulated blackbody emission of the central star. Green points are photometric observations while the short red line is the Spitzer spectrum data. The coloured curves represent the SED at the different inclinations of 50, 57, 64 and 70°. Bottom: visibility data binned by position angle of observation. The data points are split by instrument consistently with Figs. 1 and 2, where blue circles are from CHARA/CLIMB, red stars are from CHARA/CLASSIC, green triangles from PTI and orange diamonds from KI. The red curves are the results of the best fit geometricskewed ring model and the green curves are calculated from the radiative transfer image at an inclination of 50° and a position angle of 60°. The dashed bounding lines indicate the minimum and maximum model visibilities for that position angle bin. The very bottom right panel shows the observed CP measurements (black) and the CP computed from the radiative transfer image (green).

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 ± 0.15 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 $$\rho ={\rho }_1{(\frac{r_0}{r_1})}^{3/2}\eta (\chi ),$$(3)

where ρ is the density along a flow line at a distancer₀ from the symmetry axis, that is, along this disk surface. ρ₁ denotes the density at the fiducial radius of r₁, which is given as 1 au while χ is related to the cylindrical coordinate z through $$\chi =z/r.$$(4)

η(χ) is obtained from the solution of the MHD wind equations. The into-wind mass transfer is linked to the density through $$\text{begin{align}}\rho =\text{\&}1.064\times {10}^{-15}\left(\frac{{\dot{M}}_{\text{out}}}{{10}^{-7}{M}_{\odot }{\text{yr}}^{-1}}\right){\left(\frac{M_{*}}{0.5M_{\odot }}\right)}^{-1/2}\text{nn \&}\times \text{ frac}1\text{{ln}(r_{0\text{max}}/r_{0\text{min}}){\psi }_0(1-h_0\xi \text{\_}0\text{^{'}})\text{}} \text{g} {\text{cm}}^{\text{-3}},\text{nn end{align}}$$(5)

where Ṁ_out is the mass outflow rate, M_* is the stellar mass, h₀ is the disk scale height and ψ₀ is the ratio of vertical speed to Keplarian speed at the disks surface. $\xi \text{\_}0\text{^{'}}$ is related to the angle (θ₀) at which thepoloidal component of the magnetic field threads the disk, defining the opening angle of the disk wind (Safier 1993) $$\xi \text{\_}0\text{^{'}}=\text{tan}({\theta }_0).$$(6)

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. Forthe dust distribution, a constant dust-to-gas ratio is assumed to match that of the disk. The wind is populatedwith 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 Fig. 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 μm is required to produce the observed excess across the infrared. A range of grain sizes were tested from 0.1 to 1.5 μm but 0.4 μm grains provided the best fit to the K-band photometricflux values.

An into-wind mass outflow rate of 1 × 10⁻⁷ M_⊙ yr⁻¹ 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 × 10⁻⁶ M_⊙ yr⁻¹. This is unphysically high for an object such as SU Aur given the expected age of the system (5.18 ± 0.13 Myr), assuming that the rate stayed constant over full period. In addition, this is in disagreement with non-detection of Brγ 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 Fig. 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 χ² fits for the visibilities and SED for all three models are shown in Table 5.

Fig. 6 Results of radiative transfer and geometric modelling with the IN05 disk rim prescription (Isella & Natta 2005) and an additional dusty disk wind (BK12) implemented following Bans & Königl (2012). Top left: disk density cross section of the inner rim. A logarithmic colour scale is used with a minimum density of 1.00 × 10−27 to a maximum of 1.54 × 10−09. Top right: SED output of radiative transfer analysis. Dark blue curve is the simulated blackbody emission of the central star. Green points are photometric observations while the short red line is the Spitzer spectrum data. The coloured curves represent the SED at the different inclinations of 50°, 55° and 60°. Bottom: visibility data binned by position angle of observation. The data points are split by instrument consistently with Figs. 1 and 2, where blue circles are from CHARA/CLIMB, red stars are from CHARA/CLASSIC, green triangles from PTI and orange diamonds from KI. The red curves are the results of the best fit geometricskewed ring model and the green curves are calculated from the radiative transfer image at an inclination of 50° and a position angle of 60°. The dashed bounding lines indicate the minimum and maximum model visibilities for that position angle bin. The very bottom right panel shows the observed CP measurements (black) and the CP computed from the radiative transfer image (green).

Fig. 7 Computed synthetic image from TORUS following the BK12 dusty disk wind prescription. The colour indicates intensity in the K-band.

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° ± 2.2, this is in agreement with values of ${63}_{-8}^{°+4}$, ~60° and ~50° found by Akeson et al. (2005), Unruh et al. (2004), and Jeffers et al. (2014), respectively. The minor axis θ on the other hand was found to be 50.1° ± 0.2, greater than the literature values of 24° ± 23 and 15° ± 5 found by Akeson et al. (2005) and 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) and 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 Sect. 3). However, the model fit did not improve significantly by adding an off-centre 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. Ahigh 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 ${\chi }_{\text{red}}^{\text{2}}$ of 8.57. Akeson et al. (2005) find that 44 ± 9% of the total flux is in theSED 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 ± 0.02 au at an inclination of 51.2 ± 1.1°. 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 ± 0.04 au and ${62}_{-8}^{°+4}$, where the slight differences are likely due to difference between the diffuse profile, skewed ring and the standard ring structures employed. The minor axis θ of 61.0 ± 1.0° 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 ${\chi }_{\text{red}}^{\text{2}}=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 3D 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 Sect. 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 ± 0.04 au (Akeson et al. 2005) and 0.17 ± 0.08 au (Jeffers et al. 2014). The disk was also found to be at an inclination of ~50° and position angle of ~45°, 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 α_disk and β_disk were fixed such that α_disk = β_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 thelonger 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 × 10⁻⁷ M_⊙ yr⁻¹ 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 μ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 μm. This is in contrast to other inner disk studies where larger 1.2 μm 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 andits 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 ± 1.2° along a minor axis position angle of 61.0 ± 1.0° and an inner radius of 1.12 ± 0.12 mas (=0.17 ± 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 }_{\text{red}}^2$ of this model fit means the uncertainties quoted here are likely not representative of 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 μ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 observedwith ≳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: 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 χ² maps were constructed to find the best-fit location of the secondary star. The parameters are displayed in Table 6.

Appendix B: Table of photometry used in the SED fitting procedure

References

All Tables

All Figures

	Fig. 1 Coverage of the uv plane of the interferometric observations. Blue points represent observations using the CLIMB instrument, while red stars represent the CLASSIC instrument at the CHARA array. Yellow triangles represent observations using the Keck interferometer in 2011 and green triangles using the PTI with dates between 1999 and 2004.
In the text

	Fig. 2 Squared visibility and closure phase measurements against the de-projected baseline length of our interferometric data. Blue data points are those from the CLIMB/CHARA instrument. Red data points are from the CLASSIC/CHARA instrument, green data is that obtained from the PTI instrument and orange is that from the KI instrument.
In the text

Fig. 3 Top left: resultant reconstructed image showing the 3 and 5σ significance levels as solid black lines and the 1σ level as dashed lines. The beam size is shown in the bottom right. The colour bar is the same for both maps, with the maximum intensity normalised to 1 for the sake of readability. Top right: extracted asymmetry in the intensity map shown with 3σ significance level. Bottom left: fitted visibilities, data shown in black with model visibility in blue. Bottom right: fitted closure phases, data shown in black with model CP in blue. The residuals normalised to the standard deviation are plotted in the bottom of each graph.

In the text

Fig. 4 Results of radiative transfer and geometric modelling with the IN05 disk rim prescription (Isella & Natta 2005) of a gas density-dependent sublimation temperature with a single grain species size 0.1 μm . Top left: disk density cross section of the inner rim. A logarithmic colour scale is used with a minimum density of 1.00 × 10−30 to a maximum of 4.64 × 10−10. Top right:SED computed with our radiative transfer model. Dark blue curve is the simulated blackbody emission of the central star. Green points are photometric observations while the short red line is the Spitzer spectrum. The coloured curves represent the SED at the different inclinations of 50, 57, 64 and 70°. Bottom: visibility data binned by position angle of observation. The data points are split by instrument consistently with Figs. 1 and 2, where blue circles are from CHARA/CLIMB, red stars are from CHARA/CLASSIC, green triangles from PTI and orange diamonds from KI. The red curves are the results of the best fit geometricskewed ring model and the green curves are calculated from the radiative transfer image at an inclination of 50° and a position angle of 60°. The dashed bounding lines indicate the minimum and maximum model visibilities for that position angle bin. The very bottom right panel shows the observed CP measurements (black) and the CP computed from the radiative transfer image (green).

In the text

Fig. 5 Results of radiative transfer and geometric modelling with the THM07 disk rim prescription (Tannirkulam et al. 2007) of a gas density-dependent sublimation temperature with two grain species, a majority larger grains at 1.2 μm with fewer smaller grains at 0.1μm. The mass of larger grains is fixed at 9 times the mass of smaller grains. Top left: disk density cross section of the inner rim. A logarithmic colour scale is used with a minimum density of 1.00 × 10−30 to a maximum of 1.23 × 10−09. Top right:SED computed with our radiative transfer model. Dark blue curve is the simulated blackbody emission of the central star. Green points are photometric observations while the short red line is the Spitzer spectrum data. The coloured curves represent the SED at the different inclinations of 50, 57, 64 and 70°. Bottom: visibility data binned by position angle of observation. The data points are split by instrument consistently with Figs. 1 and 2, where blue circles are from CHARA/CLIMB, red stars are from CHARA/CLASSIC, green triangles from PTI and orange diamonds from KI. The red curves are the results of the best fit geometricskewed ring model and the green curves are calculated from the radiative transfer image at an inclination of 50° and a position angle of 60°. The dashed bounding lines indicate the minimum and maximum model visibilities for that position angle bin. The very bottom right panel shows the observed CP measurements (black) and the CP computed from the radiative transfer image (green).

In the text

Fig. 6 Results of radiative transfer and geometric modelling with the IN05 disk rim prescription (Isella & Natta 2005) and an additional dusty disk wind (BK12) implemented following Bans & Königl (2012). Top left: disk density cross section of the inner rim. A logarithmic colour scale is used with a minimum density of 1.00 × 10−27 to a maximum of 1.54 × 10−09. Top right: SED output of radiative transfer analysis. Dark blue curve is the simulated blackbody emission of the central star. Green points are photometric observations while the short red line is the Spitzer spectrum data. The coloured curves represent the SED at the different inclinations of 50°, 55° and 60°. Bottom: visibility data binned by position angle of observation. The data points are split by instrument consistently with Figs. 1 and 2, where blue circles are from CHARA/CLIMB, red stars are from CHARA/CLASSIC, green triangles from PTI and orange diamonds from KI. The red curves are the results of the best fit geometricskewed ring model and the green curves are calculated from the radiative transfer image at an inclination of 50° and a position angle of 60°. The dashed bounding lines indicate the minimum and maximum model visibilities for that position angle bin. The very bottom right panel shows the observed CP measurements (black) and the CP computed from the radiative transfer image (green).

In the text

	Fig. 7 Computed synthetic image from TORUS following the BK12 dusty disk wind prescription. The colour indicates intensity in the K-band.
In the text