© The Authors 2023

1 Introduction

Outflows from protoplanetary systems are one of the key mass-loss mechanisms during the planet formation process. They remove both excess mass and angular momentum from the system, which is a crucial process as the final masses and rotation rates of stars are known to be significantly less than the initial mass of protostellar cores.

Within young stellar objects (YSOs), there are several different mechanisms of outflow. Firstly, accretion and magnetospherically driven jets can emerge from the poles of the star. Secondly, photoevaporative winds caused by the ultra-violet (UV) disassociation of molecules in the upper layers of the outer disk can cause significant mass-loss. Such photoevaporative winds are usually associated with higher-mass, hotter objects such as Herbig Ae/Be stars. Finally, magnetospherically driven dusty disk winds can originate from the inner disk whereby material is lifted from the disk plane along inclined magnetic field lines. Magnetospheric winds require the presence of a strong magnetic field, usually associated with T Tauri stars with convective envelopes rather than fully radiative interiors. This allows for optically thick material to exist close enough to the central star to contribute to the NIR emission exterior to the main disk structure. This model has been shown to successfully account for the NIR excess of the spectral energy distribution (SED) and the basic visibility features of AB Aur, MWC 275, and RY Tau (Königl & Salmeron 2011; Petrov et al. 2019). While all of these mechanisms have been observed, it is not fully understood why some YSOs only appear to exhibit a sub-selection of outflow mechanisms. One of the first stars observed to have an inner dusty disk wind was the T Tauri star SU Aurigae (Petrov et al. 2019).

For a full description of the literature regarding SU Aurigae and the basic stellar properties, readers can refer to the previous paper by us (Labdon et al. 2019, henceforth LA19). In this previous work, we studied the circumstellar environment of SU Aur using interferometric observations from the centre for high angular resolution astronomy (CHARA) and Palomar testbed interferometer (PTI). The disk was found to be inclined at 51.2 ± 1.2° with a position angle of 61.0 ± 1.0° and was best modelled with a ring-like geometry with a radius of 0.17 ± 0.02 au. Additionally, radiative transfer modelling of visibilities and the spectral energy distribution (SED) allowed us to find that the NIR excess could only be reproduced in the presence of a dusty disk wind, where material is lifted from the disk along magnetic field lines allowing the reprocessing of additional stellar radiation.

Since the publication of LA19, additional relevant pieces of literature have been published. 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 at visible wavelengths. The characteristic time of change in the disk wind outflow velocity and the stellar brightness indicates 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). Recent Atacama Large Millimeter/submillimeter Array (ALMA) and spectro-polarimetic high contrast imager for exoplanets research (SPHERE) observations by Ginski et al. (2021) have revealed a significant disk warp between the inner and out disks of ~70°. This misalignment is shown to cause large shadows on the outer disk as it blocks light from the central star. Their observations also suggest that SU Aur is currently undergoing a late infall event with significant amounts of material falling inwards from the outermost regions of the disk. Such events have the opportunity to significantly impact the evolution of the disk.

This paper presents one of the first six-telescope optical inter-ferometric studies of a YSO to date utilising state-of-the-art observations covering a wider range of baseline position angles and lengths (up to 331 m) (other firsts include Kraus et al. 2020; Davies et al. 2022). Three different modelling methodologies were used to interpret our data and to provide direct comparisons to LA19. Firstly, Image reconstruction was used to obtain a model-independent representation of the data and to derive the basic object morphology. Secondly, following this geometric model fitting allowed us to gain an appreciation of the viewing geometry of the disk by fitting Gaussian and ring models to the data. In addition, more complex geometric modelling was used to explore the chromaticity of the data. Finally, we combine inter-ferometry and photometry to derive physical parameters with radiative transfer analysis, where our focus was on confirming the presence of a dusty disk wind.

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 Michigan InfraRed Combiner-eXeter (MIRC-X) instrument (Anugu et al. 2020; Kraus et al. 2018), a six-telescope beam combiner, was used to obtain observations in the near-infrared H-band (λ = 1.63 μm, Δλ = 0.35 μm) between September and October 2018. We obtained 11 independent pointings of SU Aur, using a mixture of 5 and 6-telescope configurations due to the short delay line limitations of CHARA. We obtained a maximum physical baseline of 331 m corresponding to a resolution of λ/(2Β) = 0.70 mas [milliarcseconds], where λ is the observing wavelength and B is the projected baseline. Details of our observations and the calibrator(s) observed for the target during each observing session are summarised in Table 1. The uv plane coverage that we achieved for the target is displayed in Fig. 1. Our data covers an exceptionally wide range of baseline lengths and position angles, making the data ideally suited for image reconstruction.

The MIRC-X data were reduced using the standard Python pipeline developed at the University of Michigan by (J.B. le Bouquin, N. Anugu, T. Gardner). The measured visibilities and closure phases were calibrated using interferometric calibrator stars observed alongside the target. Their adopted uniform diameters (UDs) were obtained from SearchCal (Bonneau et al. 2006, 2011), and are listed in Table 1.

Considering the short timescale over which the observations were taken, the effect of time dependencies/variability of the object is considered minimal. However, care was taken to check for time dependencies in the visibilities of baselines of similar length and position angle. Variability in the H band is known to be minimal, so any time dependencies in the visibility amplitudes are likely geometric. However, no significant time dependencies were discovered.

Fig. 1 Coverage of the uv plane of the interferometric MIRC-X observations obtained with the CHARA array.

Fig. 2 Visibilities and closure phases of the image reconstruction. Black triangles with error bars are the original calibrated observables (squared visibilities on the left and closure phases on the right), over-plotted as blue markers are the model observables of the reconstructed image. Below each plot are the fit residuals normalised by the standard deviation as black circles.

3 Image reconstruction

Image reconstruction techniques require broad and roughly circular uv coverage along as many baseline lengths as possible. Fortunately, the data from these observations lends itself to this process as the uv plane has been well sampled, though some small gaps remain in the position angle coverage. This technique is useful for the 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 S QUEEZE algorithm (Baron et al. 2010). S QUEEZE employs an MCMC approach to image reconstruction and was chosen due to the wide range of available regularisation options and its ability to implement SPARCO, a semi-parametric approach for image reconstruction of chromatic objects (Kluska et al. 2014).

In the S QUEEZE/SPARCO 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 were explored, S QUEEZE allows for a very wide range of regularisation algorithms to be implemented. Regularisation aids in completing the missing information caused by a sparse uv coverage by promoting a certain type of morphology in the image.

Total variation (TV) was found to most reliably reproduce the best image, TV 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 play an important role in image reconstruction. One cannot simply use the maximum number of pixels of the smallest size to obtain a better resolution, and 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 252 × 252 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.

The final image is shown in Fig. 3 (top left panel) with the 1,3 and 5σ significance levels shown in white, green and blue contours respectively. The inclination of the disk appears to be greater than that found by LA19 with a similar minor-axis position angle. There also appears to be a central bulge along the minor disk axis, likely caused by the over-brightness of the star along this axis. The brightness distribution shows a brighter structure along the north-west of the outer disk, parallel to the major axis of the disk. This is consistent with the asymmetry found by LA19 and indicates a highly inclined disk where the far side of the inner rim is directly exposed to the observer, while the nearside is obscured by flaring in the outer disk. There are smaller significant structures to the south-east of the disk also, we interpret these to be the shadowed near side of the rim due to their smaller extent than the northern features. We are not confident in the exact shape of these objects, given their irregularity.

In order to highlight the radial brightness distribution across the disk, Fig. 4 shows a flattened profile, the elongation of the bulge can be seen in the NW and SE directions with the extended rim material at radii out to 1 mas.

To quantitatively measure the size and orientation of the emitting region, a simple ellipse was fitted to 3, 4, and 5σ flux significance contours. The averaged results find an inclination of 46° ± 6 and a position angle of 53° ± 4. Fitting an ellipse to lower significance levels is not possible due to their irregular shape, the position angle of the ellipses is in good agreement with the values derived in Sect. 4. The inclination shows slight deviations from these values.

The chromaticity of the object is measured using two variables in the S PARCO implementation. the stellar-to-total-flux ratio at the central wavelength; and d_env, the spectral index of the extended environment. Only the spectral index of the extended environment is needed as interferometric data is only sensitive to the relative difference in spectral index between the star en circumstellar material. d_env is found to be 1.7 + 0.8, which corresponds to a temperature of ${1257}_{-231}^{+234}$ assuming the objects NIR emission is in a Rayleigh-Jean regime. This is within the range of sublimation temperatures of typical disk astronomical silicates, as expected at such small disk radii, is found to be 0.55 ± 0.26, consistent with values measured from the SED and those found from geometric modelling (see Sect. 4). The large errors associated with and d_env result from some degeneracy between the parameters; see Kluska et al. (2014) for a full description of the procedure.

The visibility and closure phase fits of the image reconstruction are shown in Fig. 2 (top panels) along with the residuals of the fit (bottom panels). The combined visibility and closure phase reduced chi-squared ${\chi }_{\text{red}}^2$ of the image reconstruction was found to be 4.38.

Fig. 3 Reconstructed and Radiative transfer images of SU Aurigae. Top Left: image reconstruction resultant bootstrapped image, including beam size and orientation. The coloured contours represent significance flux levels of 1σ (white), 3σ (green) and 5σ (blue). Top right: radiative transfer image produced using TORUS, including a dusty disk wind. Bottom middle: reconstruction of simulated data created using the best fit TORUS image above, reconstruction parameters are equivalent to those of the top-left image. Colours are normalised intensity.

Fig. 4 Azimuthal profile corrected for the inclination of the reconstructed image, created from angular slices through the image from the central star. In the full image of Fig. 3, north is up, and east is left. R is the radial distance from the central star, Θ is the polar angular direction.

4 Geometric modelling

To understand the geometry of the system one must consider the application of simple geometric models. In this section, we explore several approaches to modelling our data with both non-chromatic 'grey' models and techniques that explore chromaticity.

4.1 Basic geometric models

The fitting of Gaussian and ring-like distributions to the interferometric variables allows highly accurate estimations of the characteristic size, inclination and position angle of the object. In all models, the central star is modelled as a point source, which is an acceptable assumption given the expected angular diameter of the star. The disk parameters are then fitted in the RAPIDO (Radiative transfer and Analytic modelling Pipeline for Interferometric Disk Observations) framework (Kreplin et al. 2018), available in-house at the University of Exeter. RAPIDO utilises the Markov chain Monte Carlo (MCMC) sampler emcee to produce a fit and error estimate (Foreman-Mackey 2016). Three disk models were employed, a standard Gaussian brightness distribution which is characterised by its full-width-half-maximum (FWHM), along with two ring models, a sharp ring with a width fixed to 20% of the disk radius (R) and a 'skewed' ring with a more diffuse radial profile produced by convolving with a Gaussian with a FWHM. The skewed ring is also capable of modelling azimuthal modulation or disk asymmetries, a detailed description of this model can be found in Lazareff et al. (2017). In addition to the model-specific parameters, we also fitted the inclination (INC), minor-axis position angle (PA) and disk-to-total flux ratio (f_disk)- As we see no evidence of time variability in the data we are able to fit all data simultaneously. The results from the simple geometric model fitting are shown in Table 2.

Out of the geometric models tested, the Gaussian model is considered to be the best fit. Even though the skewed ring produced a slightly small χ² value for the closure phase and visibility measurements, we do not consider this significant, given the additional degrees of freedom in the model. In addition, as shown in Fig. 5, the best fit skewed ring tends towards a Gaussian distribution, with a very small inner radius 0.17 mas with a wide FWHM ring width 0.75 mas. The skewed ring fails to reproduce the small but complex asymmetries seen in the image reconstruction. We attribute this to a failure of the underlying ring model rather than a lack of asymmetries within the disk; if a simple ring model does not represent the data well, introducing asymmetries cannot be expected to improve the result.

The best fit Gaussian model finds a disk of FWHM 1.52 ± 0.01 mas which is inclined at 56.9 ± 0.4° and a minor-axis position angle of 55.9 ± 0.5°. In addition, we find that 43 ± 1% of the total flux originates from the disk in the Η band. This is consistent with measurements based on the infrared excess of the spectral energy distribution (SED) in LA 19 and also with the flux ratio found in the image reconstruction algorithm.

The primary limitation of the simple geometric models described above is that they are intrinsically 'grey' in nature. Meaning they contain no spectral information, hence the large χ² values obtained in the fitting process. In order to better model the spectral dependency of the visibility, more complex temperature gradient models that are able to account for observing wavelength must be employed.

4.2 Temperature gradient models

A physically correct model can be applied by considering the temperature gradient of the disk. A temperature gradient model (TGM) allows for the simultaneous fitting of interferometric and photometric observables. The origin of the photometric data used is described in Table A.1. It is built up by several rings extending from an inner radius R_in to an outer radius R_out. Each ring is associated with temperature and hence flux. Therefore, a model SED can be computed by integrating over the resulting blackbody distributions for each of the concentric rings. Such a model allows us to not only build up a picture of the temperature profile but also approximate the position of the inner radius. The TGM is based upon a T_R = T₀(R/R₀)^−Q profile where T₀ is the temperature at the inner radius of the disk R₀, and Q is the exponent of the temperature gradient (Kreplin et al. 2020; Eisner & Hillenbrand 2011). A TGM represents an intrinsically geometrically thin disk. A point source is used at the centre of each model to represent an unresolved star, which is a reasonable approximation given the expected angular diameter of 0.05 mas (Pérez et al. 2020). Also included in the fitting of photometric parameters is a treatment of interstellar extinction based on Fitzpatrick (1999) with an E_(B−V) = 0.5 (Bertout et al. 2007).

The inclination and position angle of the disk is maintained at fixed values of 56.9 ± 0.4° and 55.9 ± 0.5° respectively, from the fitting of the Gaussian distribution. This was done to reduce the number of free parameters in the model. The fitting was undertaken using all of the visibility data shown in Fig. 3 and all the SED points simultaneously. The fitting and error computation was once again done using emcee (Foreman-Mackey 2016).

The results of the temperature gradient modelling are shown in Fig. 7. We find an inner disk radius of 0.15 ± 0.04 au where the temperature is equivalent to 2100 ± 200 Κ and decreases with an exponent of Q = 0.62 ± 0.02. The inner disk radius is the point at which the dusty disk is truncated due to the sublimation of material, in contrast to more extreme objects such as FU Ori where the inner disk radius is equivalent to that of the stellar radius indicating boundary layer accretion (Labdon et al. 2021). Interior to the sublimation radius is expected a hot dust-free inner disk from which material can be magnetospherically accreted. However, our low spectral resolution continuum observations are not sensitive to these regions. The outer edge of the temperature gradient model is a fitted parameter found to be 0.2 ± 0.03 au. This does not correspond to a physical outer edge of the disk but to a potential temperature discontinuity, such as where a viscously heated region ends. This is also the approximate region at which the NIT visibilities are no longer sensitive.

An inner disk radii temperature of 2100 ± 200 Κ may be considered high, but is broadly consistent with laboratory sublimation temperatures for silicate grains. An exponent of Q = 0.62 ± 0.02 is slightly larger than that expected from a disk heated by stellar radiation alone and may indicate the presence of additional disk heating mechanisms, such as viscous heating (Pringle 1981; Kenyon & Hartmann 1987; Dullemond & Dominik 2004). A higher-than-expected sublimation temperature can also be explained by the presence of refractory grains which sublimate at a high temperature than standard astronomical silicates (Benisty et al. 2010)

The resultant SED of the inner disk is shown in Fig. 6 as the orange curve. Beyond 5–6 μm the TGM fails to fit the shape of the disk accurately. This is likely due to the flared nature of the disk in contrast to the 'flat' TGM model which will have the strongest effect and longer wavelengths. In addition, the strong disk warp reported by Ginski et al. (2021) would also introduce some temperature discontinuity, the modelling of which is beyond the scope of this work.

Fig. 5 Geometric model images, corresponding to best-fit parameters described in Table 2. Top left: Gaussian model brightness distribution. Top right: ring model brightness distribution. Bottom left: skewed Ring model brightness distribution. Bottom right: full reconstructed image, repeated from Fig. 3.

Fig. 6 Spectral energy distribution of SU Aurigae. Green points are photometric data from a variety of instruments. Red line is Spitzer IR data. The Black dashed line is direct radiation from the stellar photosphere. The blue line is the best TORUS computed radiative transfer model inclined at 56°. The orange line is the SED computed from the simple temperature gradient models described in Sect. 4.2.

Fig. 7 Temperature gradient profile across the inner au of SU Aur. The dusty disk extends down an inner radius, consistent with the expected sublimation radius and temperature. Interior to this expected a hot dust-free inner disk which is not detected in the continuum observations.

5 Radiative transfer

In order to provide a more physical model and directly compare to LA 19, we used the Transport Of Radiation Using Stokes intensities (TORUS) Monte-Carlo radiative transfer code (Harries et al. 2019) to simultaneously fit the visibility, closure phase and photometric data of the SU Aurigae system. The models adopted here are based on the disk models used by LA 19, adapted to account for the higher inclination and different observing wavelengths. In these TORUS simulations, the dust was allowed to vertically settle to the scale height of the gas component 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). 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. We adopt a silicate grain species with dust properties and opacities adopted from Draine (2003). For a more detailed description of TORUS and the algorithms used, see Davies et al. (2018), Labdon et al. (2019).

The dusty disk wind model is adapted from Bans & Königl (2012). This mechanism is based on the presence of a large-scale, ordered magnetic field which threads the disk along which disk material is flung out. 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. A full description of the implementation within the TORUS radiative transfer code can be found in LA 19.

The disk model adopted follows the curved inner rim prescription of Isella & Natta (2005) with a density-dependent sublimation radius whereby grains located in the disk midplane are better shielded due to higher densities and so can exist closer to the central star than grains in the less dense upper layers. A full summary of the disk parameters can be found in Table3. The only parameters we explored with respect to LA19 were the disk scale height (h₀), the grain size (a) and the flaring index β_disk. The key difference in the models described here compared to LA19 is the grain size adopted. Here we adopt a smaller grain size of 0.14 μm, which in turn leads to a slightly smaller inner radius of 0.16 au. Additionally, in order to improve the SED fit at longer wavelengths we also adopt a more modest scale height of 9.0 au. The resultant SED from the radiative transfer model is shown in Fig. 6, in addition to the model stellar photosphere, also calculated within TORUS. The SED fits well across the optical and IR and longer mm wavelengths, however, is a relatively poor fit across the 8–40 μm range. We attribute this to the disk warp reported in Ginski et al. (2021) which would result in a physical disk break and temperature discontinuity. The geometry governing disk warps is little understood, in particular, how it would affect the inner and outer edges of the warp. The shapes of these rims, and how is the vertical structure of the disk affected at this point, are not known. Such a model is well beyond the scope of this paper for 2 reasons: firstly we lack high spatial resolution data at M/F-IR wavelengths which might cover the location of the disk warp. The complex geometry of the warp would require extensive modelling, which is difficult given the limitations of the radiative transfer code used.

The presence of a disk wind is once again required to fit both the visibilities and the SED. The absence of a disk wind fails to reproduce IR excess across both the H and K bands, with insufficient NIR disk flux. This significantly impacts both the SED and interferometric fit. A disk wind is required to eject more hot dust above the midplane of the disk where it is directly exposed to stellar radiation which is reprocessed as an IR excess.

This is shown in Fig. 8, where the squared visibilities are shown for both a disk model with and without the dusty wind environment. The disk wind model provides a far superior fit to the observations, being able to successfully reproduce the NIR excess. Beginning from the wind parameters found in LA19, we run a new grid of models exploring the disk wind parameter space. The parameters explored and their ranges are described in Table 3, along with the resulting best-fit values. The results are broadly similar to those found in LA19, with the exception of a slightly higher temperature of material in the disk wind, closer to dust temperatures found in the inner disk. Of particular note is the into-wind accretion rate of 10⁻⁷ M_⊙. Considering the historically accepted on-to-star accretion rate to into-wind accretion ratio of 0.1, this level of transport is perhaps unphysically high given the age of the star, but a discussion of possible mechanisms are included in Sect. 6.

The final computed image is shown in Fig. 3 (top left panel) and shows the clear asymmetry originating from the inclination. In order to approximate what this computed image would look like if observed at the same resolution as the original observations, we computed synthetic visibility and closure phases based on the radiative transfer images (as described in Davies et al. 2018). Artificial noise and error bars were computed to be representative of the original data and to ensure an accurate representation. These synthetic observables were then reconstructed in the same manner as the original data, as described in Sect. 3. Care was taken to ensure consistency in the reconstruction parameters for both the real and synthetic observables. The reconstructed TORUS image is also shown in Fig. 3 (bottom right panel), and shows clear similarities with both the original TORUS image and the image reconstructed from the original data.

6 Discussion

Our extensive observations and analysis of the circumstellar environment of SU Aurigae have revealed the details of the inner disk in unprecedented detail. The wide variety of techniques used to analyse our interferometric data allows us to precisely define the disk characteristics.

Image reconstruction is a crucial, model-independent, method of analysis which is ideally suited to our dataset with extensive uv and baseline coverage. Our analysis reveals an elliptical shape, indicative of an object with a high inclination as shown in Fig. 3. There appears to be a central bulge to the disk, however, this feature is not thought to be physical but rather a manifestation of the brightness of the central star combined with the width of the disk at this point. The more extended material in the image is thought to be a depiction of the far side of the disk rim, which is unobscured by the outer disk, this is reinforced by the fact that the extended material is bright in the northeast of the image. There is a significant asymmetric feature in the form of a thin brightness on the north-eastern edge of the disk. The unique shape of this feature indicates that this is again caused by the high inclination obscuring the nearside disk rim. The effect of an inclined disk on the observed brightness distribution is described extensively by Jang-Condell & Turner (2013) and accurately describes the observations here. By fitting an ellipse to the significant features in the image, we gain a quantitative measure of the disk geometry. The results find an inclination of 46° ± 6 and a position angle of 53° ± 4, very similar to those derived from more reliable geometric model fits described below. The scale and shape of the image are similar to that of LA19, with a slightly higher inclined viewing angle. We consider the images of this work to be the more accurate depiction of SU Aur given the higher quality observations, taken over a much shorter timescale, with the added detail and resolution this entails.

Although geometric modelling is much more constrained in the geometries it can explore, it provides a more quantitative view of the disk. It was found that the model which best fit our data was a simple Gaussian distribution with a point source representing the star. A Gaussian model is consistent with other work, both on this object by LA19, but also in other YSO studies such as the survey by Lazareff et al. (2017), who find that little under half of their 51 objects can be modelled by a Gaussian structure. The Gaussian fitted in this work has an FWHM of 1.52 ± 0.01 mas (0.239 ± 0.002 au) at an inclination of 56.9° ± 0.4 at a minor axis position angle of 55.9° ± 0.5 with a stellar-to-total flux ratio of 0.57 ± 0.01. The reduced; χ² value for the visibilities is 11.63 and 6.05 for the closure phases, which are equivalent to 0^° for this centrosymmetric model. These values are in agreement with the literature values of Akeson et al. (2005) who find a K band radius of 0.18 ± 0.04 au and an inclination of $62{°}_{-8}^{+4}$. Similar values for the inclination in literature are ~60° and ~50° found by Unruh et al. (2004), Jeffers et al. (2014) respectively. The minor axis position angle derived here is significantly greater than the literature values of 24° ± 23 and 15° ± 5 found by Akeson et al. (2005), Jeffers et al. (2014). This difference is likely due to either: The poor uv coverage and the 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, which is shown by Ginski et al. (2021) to be misaligned compared to the inner disk.

The geometric modelling results are broadly similar to those presented in our previous work LA19, where an inclination of 50.9 ± 1.0° and minor axis position angle of 60.8 ± 1.2° were found, and the data were marginally better described by a ringlike brightness distribution. The values and models presented here are considered to be more accurate due to high precision observations and significantly smaller potential for temporal variations, as the data of LA19 was coalesced over 14 yr. These values are also consistent with observations of the outer disk by Ginski et al. (2021) where dark shadows are observed in scattered light originating from a significant disk warp between the inner and outer regions. On larger scales, the near side of the disk is seen to the northeast; in our observations, it is seen to the south-west.

In modelling the temperature gradient of SU Aur, we can gain an appreciation for the spectral dependence of our inter-ferometric variables across the 6 spectral channels of MIRC-X. Our modelling finds a disk which extends down to 0.15 ± 0.04 au where the temperature is equivalent to 2100 ± 200 K and decreases with an exponent of Q = 0.62 ± 0.02. The outer edge of this temperature regime was found to be 0.20 ± 0.03 au, showing this prescription only covers the innermost regions of the disk. Modelling outer regions of the disk are beyond the scope of this paper, as our NIR interferometric data does not cover emission from these regions and larger radii are likely affected by the significant disk warp causing a temperature gradient discontinuity. A temperature gradient exponent of Q = 0.62 ± 0.02 found here lies between two established models from the literature. That of Pringle (1981), who find that a steady state, optically thick accretion disk heated by viscous processes will exhibit an exponent of 0.75 and of Kenyon & Hartmann (1987), Dullemond & Dominik (2004), who show that a flared disk heated by reprocessed stellar radiation alone will exhibit an exponent of ≤0.5. As such, it is difficult to comment on the heating of the circumstellar environment of SU Aurigae. In particular, given the strong dusty disk wind environment, which likely impacts the observed temperature gradient. But it is possible that the disk is not heated by stellar radiation alone, but additional heating processes, such as viscous heating, may also be present.

The radiative transfer modelling presented in this paper is heavily based on LA19. For a detailed discussion of the motivation behind certain choices, particularly in relation to the shape of the inner rim, we recommend the reader see Sect. 5 of that paper. In this work, we can achieve a similar SED fit to LA19, including the adoption of a smaller 0.14 μm grain size. It is noted that the smaller grain size is in line with older radiative transfer work of SU Aur by Akeson et al. (2005). The inner rim extends down to 0.13 au at a sublimation temperature of 2000 K, compared to 0.15 au from LA19. This is within the uncertainties of values predicted by the temperature gradient modelling and is roughly consistent with older literature values of 0.18 ± 0.04 au and 0.17 ± 0.08 au by Akeson et al. (2005) and Jeffers et al. (2014) respectively. The relationship between grain size and inner rim radius is within the TORUS radiative transfer framework is discussed in Davies et al. (2020). The flaring parameters α_disk and β_disk were fixed such that α_disk = β_disk + 1 and found to be 2.3 and 1.3 respectively, a more physical representation, than the values depicted in LA19. Similarly to LA19, a dusty disk wind is required to fit both the SED across the NIR and visibilities, as shown in Fig. 8.

The TORUS implementation of s dusty wind does not depend on the underlying launching mechanism but just prescribes a geometry above and below the disk, which is populated by dust grains, where it can reprocess stellar radiation to contribute to the NIR excess. This was first proposed by Bans & Königl (2012) in the context of a magnetospherically driven disk wind. Figure 8 highlights how a standard (no-wind) model cannot sufficiently fit the interferometric data and how adding a dusty disk wind can provide a significantly better fit to the data. The maximum temperature of the dust in the wind is similar to the temperature of the dust at the sublimation temperature of the disk 1900 and 2000 K respectively, as expected given the dust launches from close to to the sublimation rim. In Appendix B, each baseline is plotted in a separate panel to explore the chromatic (temperature) gradient of the data. The simplified heating description can explain the discrepancy in the gradient and level of some baselines within the radiative transfer model. TORUS produces a disk heated by reprocessed stellar radiation alone, with no internal disk heating such as viscous heating. We see from the specific temperature gradient modelling that we might expect some viscous heating within the disk.

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. If one assumes the historically accepted outflow-to-accretion ratio of 0.1, the resulting onto-star-accretion rate is greater than those typically found in T Tauri stars. In particular, this contrasts with the measured accretion rate of SU Aur (Pérez et al. 2020). However, the inflow-to-outflow ratio is the subject of some discussion, with recent works suggesting the ratio may be closer to unity in some situations (Pascucci et al. 2023), particularly by invoking magnetically driven winds originating from a dead zone close to the sublimation radius. Based on our data, we cannot differentiate such specific wind launching mechanisms. In addition, the suggestion that SU Aur is undergoing a late-stage infall event Ginski et al. (2021) could be a potential explanation for such a high level of mass transport through the system.

Late infall events typically occur after the depletion of the protostellar envelope, when a large amount of material exterior to the disk falls inwards. The triggering mechanisms behind such events are poorly understood as there are few observational examples. From simulations, it is thought that a large infall event can not only directly increase the accretion rate (Dullemond et al. 2019), but also dramatically re-sculpt the disk, creating gaps, spiral arms and strong misalignments (Kuffmeier et al. 2021). In the case of SU Aur a large infall of material is thought to have caused a strong disk warp between the inner and outer disk (Ginski et al. 2021). It is also not unreasonable to assume this lead to an enhanced accretion rate, which would allow for the strong disk winds needed in the model presented here.

The simplest way to directly compare the analytical methods employed in this work is through the images produced. Figure 3 shows a collage image reconstruction and radiative transfer methods, while Fig. 5 compares the various geometric images. They clearly show the similarities between the images with very similar inclinations and position angles with asymmetries oriented in the same direction. The asymmetries all appear to result from the high inclination of the disk, causing the near side inner rim to be shadowed by the flared outer disk. The lowest panel of Fig. 3 shows an image reconstructed from synthetic visibilities and closure phases obtained from the radiative transfer image. Care was taken to ensure the observables were matched in baseline length and position angle and detector mode, while the image reconstruction process was exactly that described in Sect. 3. The remarkable similarity between the different images, in particular the real image reconstruction and the synthetic reconstruction, provides strong reinforcement that the models adopted are truly representative of the disk structure.

Fig. 8 Coloured circles are the MIRC-X squared visibilities, where the colour represents the wavelength across the H-band. Black crosses are the corresponding squared visibilities extracted from TORUS radiative transfer model images. Left: best fit TORUS model, including a strong dusty disk wind component. Right: the same base disk model, excluding any dusty disk wind component. Lower panels show the residuals in the fits between data and models.

Fig. 9 Black crosses are the closure phase data obtained with MIRC-X. Overlaid as coloured points are the TORUS model closure phases; the colours represent the wavelength of the spectral channels and follows the same convention as other plots in this work. Below in black points are the normalised residual errors of the fit.

7 Future observations

Recent improvements at CHARA include the commissioning of the MYSTIC (Michigan Young STar Imager at CHARA) beam combiner (Monnier et al. 2018), a 6-telescope K-band combiner capable of working in parallel with MIRC-X for simultaneous observations. The combination of MIRC-X and MYSTIC would allow us to investigate the temperature structure of the inner disk in greater detail while also obtaining greater precision due to the higher sensitivity of MYSTIC. In addition, longer wavelength mid-infrared (MIR) interferometry with the MATISSE (Multi AperTure mid-Infrared SpectroScopic Experiment) instrument at the Very Large Telescope interferometer (VLTI) would potentially allow us to resolve the warped region of the disk, to confirm the nature of the geometry in this region for direct comparison to late infall event simulations.

8 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 the literature. We see evidence of 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 simple geometric model fits show that the circumstel-lar environment is best modelled as a Gaussian distribution with a disk of inclination 56.9 ± 0.4° along a minor axis position angle of 55.9.0 ± 0.5° and an FWHM of 1.52 ± 0.01 mas (0.239 ± 0.002 au). Such geometry is consistent with strong disk shadows observed in the outer disk originating from a disk warp.

• We model the radial temperature profile of the inner disk and find a disk which extends down to 0.15 ± 0.04 au where the temperature is equivalent to 2100 ± 200 K and decreases with an exponent of 0.62 ± 0.02.

• A dusty disk wind scenario is still required to 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. The high accretion rate required to reproduce the stellar-to-total flux ratio could be explained by a late infall event.

• Our best-fit model (dusty disk wind model) suggests that the dust composition in the disk is dominated by small grains (0.14 μm) with a sublimation temperature of 2000 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 (9 at 100 au).

We have shown that dusty disk winds and late infall events have profound impacts on protoplanetary disks. These processes will have implications for the efficiency of planet formation within the disk, the study of which will require more observations of similar objects.

Acknowledgements

We acknowledge support from an STFC studentship (No. 630008203) and a European Research Council Starting Grant (Grant Agreement No. 639889). MIRC-X received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant No. 639889). J.D.M. acknowledges funding for the development of MIRC-X (NASA-XRP NNX16AD43G, NSF-AST 1909165) and MYSTIC (NSF-ATI 1506540, NSF-AST 1909165). This work is based upon observations obtained with the Georgia State University Center for High Angular Resolution Astronomy Array at Mount Wilson Observatory. The CHARA Array is supported by the National Science Foundation under Grant No. AST-1636624 and AST-2034336. Institutional support has been provided by the GSU College of Arts and Sciences and the GSU Office of the Vice President for Research and Economic Development. 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 Additional table of photometry

Table of photometry used in the SED fitting procedure.

Appendix B Additional figure

Fig. B.1 Squared visibilities split by baseline/pointing against wavelength. Black crosses are the measured visibilities. Red lines represent synthetic visibilities from the final TORUS radiate transfer model. The legend of each plot is the length of the individual baseline.

References

All Tables

All Figures

	Fig. 1 Coverage of the uv plane of the interferometric MIRC-X observations obtained with the CHARA array.
In the text

	Fig. 2 Visibilities and closure phases of the image reconstruction. Black triangles with error bars are the original calibrated observables (squared visibilities on the left and closure phases on the right), over-plotted as blue markers are the model observables of the reconstructed image. Below each plot are the fit residuals normalised by the standard deviation as black circles.
In the text

Fig. 3 Reconstructed and Radiative transfer images of SU Aurigae. Top Left: image reconstruction resultant bootstrapped image, including beam size and orientation. The coloured contours represent significance flux levels of 1σ (white), 3σ (green) and 5σ (blue). Top right: radiative transfer image produced using TORUS, including a dusty disk wind. Bottom middle: reconstruction of simulated data created using the best fit TORUS image above, reconstruction parameters are equivalent to those of the top-left image. Colours are normalised intensity.

In the text

	Fig. 4 Azimuthal profile corrected for the inclination of the reconstructed image, created from angular slices through the image from the central star. In the full image of Fig. 3, north is up, and east is left. R is the radial distance from the central star, Θ is the polar angular direction.
In the text

	Fig. 5 Geometric model images, corresponding to best-fit parameters described in Table 2. Top left: Gaussian model brightness distribution. Top right: ring model brightness distribution. Bottom left: skewed Ring model brightness distribution. Bottom right: full reconstructed image, repeated from Fig. 3.
In the text

Fig. 6 Spectral energy distribution of SU Aurigae. Green points are photometric data from a variety of instruments. Red line is Spitzer IR data. The Black dashed line is direct radiation from the stellar photosphere. The blue line is the best TORUS computed radiative transfer model inclined at 56°. The orange line is the SED computed from the simple temperature gradient models described in Sect. 4.2.

In the text

	Fig. 7 Temperature gradient profile across the inner au of SU Aur. The dusty disk extends down an inner radius, consistent with the expected sublimation radius and temperature. Interior to this expected a hot dust-free inner disk which is not detected in the continuum observations.
In the text

Fig. 8 Coloured circles are the MIRC-X squared visibilities, where the colour represents the wavelength across the H-band. Black crosses are the corresponding squared visibilities extracted from TORUS radiative transfer model images. Left: best fit TORUS model, including a strong dusty disk wind component. Right: the same base disk model, excluding any dusty disk wind component. Lower panels show the residuals in the fits between data and models.

In the text

	Fig. 9 Black crosses are the closure phase data obtained with MIRC-X. Overlaid as coloured points are the TORUS model closure phases; the colours represent the wavelength of the spectral channels and follows the same convention as other plots in this work. Below in black points are the normalised residual errors of the fit.
In the text

	Fig. B.1 Squared visibilities split by baseline/pointing against wavelength. Black crosses are the measured visibilities. Red lines represent synthetic visibilities from the final TORUS radiate transfer model. The legend of each plot is the length of the individual baseline.
In the text