<figure><img src="image_1.png" /></figure>

# Imaging the Inner Astronomical Unit of the Herbig Be Star HD 190073

Nour Ibrahim<fnref n="1" />, John D. Monnier<fnref n="1" />, Stefan Kraus<fnref n="2" />, Jean-Baptiste Le Bouquin<fnref n="3" />, Narsireddy Anugu<fnref n="4" />, Fabien Baron<fnref n="3" />,<fnref n="5" />, Theo Ten Brummelaar<fnref n="4" />, Claire L. Davies<fnref n="2" />, Jacob Ennis<fnref n="1" />, Tyler Gardner<fnref n="1" />, Aaron Labdon<fnref n="6" />, Cyprien Lanthermann<fnref n="4" />, Antoine Mérand<fnref n="7" />, Evan Rich<fnref n="1" />, Gail H. Schaefer<fnref n="4" />, and Benjamin R. Setterholm<fnref n="1" />

## Abstract

The inner regions of protoplanetary disks host many complex physical processes such as star–disk interactions, magnetic fields, planet formation, and the migration of new planets. To study directly this region requires milliarcsecond angular resolution, beyond the diffraction limit of the world's largest optical telescopes and even too small for the millimeter-wave interferometer Atacama Large Millimeter/submillimeter Array (ALMA). However, we can use infrared interferometers to image the inner astronomical unit. Here, we present new results from the CHARA and VLTI arrays for the young and luminous Herbig Be star HD 190073. We detect a sub-astronomical unit (sub-AU) cavity surrounded by a ring-like structure that we interpret as the dust destruction front. We model the shape with six radial profiles, three symmetric and three asymmetric, and present a model-free image reconstruction. All the models are consistent with a near face-on disk with an inclination $\le 20^\circ$, and we measure an average ring radius of $1.4 \pm 0.2$ mas (1.14 au). Around 48% of the total flux comes from the disk with 15% of that emission appearing to emerge from inside the inner rim. The cause of emission is still unclear, perhaps due to different dust grain compositions or gas emission. The skewed models and the imaging point to an off-center star, possibly due to binarity. Our image shows sub-AU structure, which seems to move between the two epochs inconsistently with Keplerian motion and we discuss possible explanations for this apparent change.

*Unified Astronomy Thesaurus concepts: Young stellar objects (1834); Herbig Ae/Be stars (723); Protoplanetary disks (1300); Long baseline interferometry (932)*

## 1. Introduction

As the catalogs of exoplanets grow, we increasingly see evidence of otherworldly environments that we do not observe in our local solar system (see Rein 2012; Cacciapuoti et al. 2022). Observing the early formation stages of planetary systems is key to understanding and eventually predicting the worlds we continually find. The planet formation process is still not well understood, especially around massive stars. Herbig Ae/Be stars are a class of intermediate- to high-mass (2–10 $M_\odot$) pre-main-sequence stars of spectral type A or earlier that exhibit excess near-infrared (NIR) emission, which is associated with circumstellar disks. Disks around the more massive B stars are harder to study because they disappear on shorter timescales, around a few $10^5$ yr, likely due to photoevaporation caused by UV radiation from the star. Herbig Be stars, in particular, are often observed with a factor of 5–10 lower disk masses than Herbig Ae and T Tauri stars, or in some cases, are not detected at all (Alonso-Albi et al. 2009). The inner regions of these disks are particularly interesting because they host many complex processes that go into planet formation, yet they are not well studied. Imaging sub-astronomical unit (sub-AU) scales requires milliarcsecond angular resolution, which is not

achievable using the world's largest optical telescopes and is even too small for the millimeter-wave interferometer Atacama Large Millimeter/submillimeter Array (ALMA). However, IR long-baseline interferometers can probe disks at the sub-AU scale and give us a better look at the inner regions of the disks.

Modeling has been a key tool for studying the structure of circumstellar accretion disks. The early models by Hillenbrand et al. (1992), which assumed a flat and optically thick disk extending a few stellar radii close to the star, were able to reproduce photometric NIR excess emission measurements. However, this theoretical picture was not confirmed by the interferometric observations carried out by Millan-Gabet et al. (1999), which found that the disk sizes had to be much larger than what was predicted by the thin and optically thick disk models. In 2001, Natta et al. (2001) and Dullemond et al. (2001) proposed a new model for the inner disk region, suggesting an optically thin cavity around the star and a “puffed-up” inner rim wall that is truncated at a sub-AU radius where the rim temperature is equal to the dust sublimation temperature. Further observations by Millan-Gabet et al. (2001) and Tuthill et al. (2001) showed that this model, where the majority of the NIR excess originates from the dust sublimation radius, explained the large photometric NIR bump. Monnier & Millan-Gabet (2002) validated the “puffed-up” inner rim with the introduction of a size–luminosity diagram where they also showed that the dust sublimation temperature was between 1500 and 2000 K for their sample. However, more recent submilliarcsecond interferometry observations of AB Aurigae

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Astronomy Department, University of Michigan, Ann Arbor, MI, 48109, USA; inoura@umich.edu</footnote>

<footnote n="2">² Astrophysics Group, Department of Physics & Astronomy, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK</footnote>

<footnote n="3">³ Institut de Planetologie et d'Astrophysique de Grenoble, Grenoble F-38058, France</footnote>

<footnote n="4">⁴ The CHARA Array of Georgia State University, Mount Wilson Observatory, Mount Wilson, CA, 91203, USA</footnote>

<footnote n="5">⁵ Department of Physics and Astronomy, Georgia State University, 25 Park Place NE 605, Atlanta, GA, 30303-2911, USA</footnote>

<footnote n="6">⁶ European Southern Observatory, Casilla 19001, Santiago 19, Chile</footnote>

<footnote n="7">⁷ European Organisation for Astronomical Research in the Southern Hemisphere (ESO), Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany
*Received 2022 August 13; revised 2023 January 17; accepted 2023 January 18; published 2023 April 24*</footnote>
<table><caption>**Table 1**
HD 190073 Stellar Properties from Guzmán-Díaz et al. (2021)</caption><thead><tr><th>Property</th><th>HD 190073</th></tr></thead><tbody><tr><td>α (J2000)</td><td>20h03m02s51</td></tr><tr><td>δ (J2000)</td><td>+5°44′16″66</td></tr><tr><td>Spectral Type</td><td>B9</td></tr><tr><td>T<sub>eff</sub></td><td>9750 ± 125 K</td></tr><tr><td>R<sub>*</sub></td><td>9.68 ± 0.44 R<sub>☉</sub></td></tr><tr><td>Age</td><td>0.30 ± 0.02 Myr</td></tr><tr><td>Distance</td><td>824.16 ± 21.83 pc</td></tr><tr><td>log(L<sub>*</sub>) L<sub>☉</sub></td><td>2.88 ± 0.03</td></tr><tr><td>Mass</td><td>6.0 ± 0.2 M<sub>☉</sub></td></tr><tr><td>V mag</td><td>7.79 ± 0.06</td></tr><tr><td>H mag</td><td>6.61 ± 0.07</td></tr></tbody></table>

and MWC 275 (HD 163296) carried out by Tannirkulam et al. (2008) showed that models in which the dust evaporation rim solely produces the NIR excess fail to explain the data. They found that a significant amount of the inner emission emerges from within the dust sublimation radius which likely does not have a sharp edge. Lazareff et al. (2017)'s Precision Integrated-Optics Near-infrared Imaging ExpeRiment (PIONIER) survey confirmed this general picture in addition to adding 27 new objects to the size-luminosity diagram by Monnier & Millan-Gabet (2002) and found some variations between sources.

Most of what we know about Herbigs comes from observations of the more common and close-by Herbig Ae stars. Due to the rareness of massive stars and the very short disk lifetimes, only a handful of Herbig Be stars with detected disks are close enough for us to study, one of which is the B9 star HD 190073, also known as V1295 Aquila (see the stellar parameters in Table 1) (Rich & Monnier 2022). HD 190073 is a pre-main-sequence star with a spectral type of B9 and a temperature of approximately 9750 K. It has a luminosity of ~760 $L_{☉}$ and a model-derived mass of $6.0 \pm 0.2 M_{☉}$. Based on its mass, this star is expected to contract rapidly onto the main sequence with a spectral type of B5 and a temperature of 15,500 K in ~10⁵ yr (Cox 2000). While HD 190073 shares some similarities with other well-studied pre-main-sequence stars such as Herbigs AB Aur and MWC 275, it has more than twice the mass and will likely have a significantly different experience in terms of its disk evolution. It also has a detected magnetic field which is uncommon among Herbigs (Alelian et al. 2013). The magnetic field is much weaker than those of T Tauri stars, but it could be measured due to the star's narrow spectral lines. The low $v \sin i$ 8.6 km s⁻¹ is partly explained by a face-on geometry (Catala et al. 2007). Unfortunately, ALMA has not observed it yet so we cannot confirm through imaging that the disk is face-on. 2D models from the PIONIER/Very Large Telescope Interferometer (VLTI) survey found that HD 190073 is nearly face-on and symmetric, but they did not have enough resolution to see any inner emission (Lazareff et al. 2017). The results from this PIONIER study revealed that the inclination of the HD 190073 disk is <20°, and therefore it is plausible that the slow $v \sin i$ is due to the low inclination.

Setterholm et al. (2018) recently presented the first Center for High Angular Resolution Astronomy (CHARA) results on HD 190073 with three times higher resolution. They used broadband data taken over many years and combined them as a single epoch. The resulting visibility curves lacked the expected bounce at larger baselines associated with a ring-like structure. The HD 190073 models showed a face-on disk that

favored a Gaussian shape rather than a thin ring. This conclusion was similar to the results from Tannirkulam et al. (2008) for the older, less massive, and less luminous Herbig Ae stars AB Aur and MWC 275.

Here in this work, we will revisit this analysis using new “snapshot” observations taken over a much shorter period of three months and with a much denser (u, v) coverage. We obtained the new H-band interferometry using the updated Michigan InfraRed Combiner—eXeter (MIRC-X) (Anugu et al. 2020) instrument at the CHARA Array (ten Brummelaar et al. 2005) and PIONIER (Le Bouquin et al. 2011) at the VLTI. With improved angular resolution and increased (u, v) coverage, we present new models with six different radial profiles along with a model-free image reconstruction of the sub-AU ring. The enhanced data quality allows us to model finer details in the disk structure and interpret the inner region dynamics more accurately.

We start by describing our observation routines to collect the CHARA and VLTI data, as well as discussing our data reduction techniques in Section 2. Next, we explain our methods for combining the data from the two instruments into two epochs, and we show a basic presentation of the (u, v) coverage, visibility, and closure phase (CP) measurements for each epoch in Section 3. We present our symmetric and asymmetric models in Section 4 and describe them extensively. Then in Section 5 we present model-free images and compare them to the models. Finally, we discuss in Section 6 what the results of the modeling and imaging reveal about the inner disk of HD 190073.

## 2. Observations and Data Reduction

This paper relies on new IR interferometry data from two different facilities, CHARA and VLTI. We coordinated observations between Mt. Wilson in California and Paranal Mountain in Chile using similar wavelengths for one month to get better-quality data. It was also important to get data around the same time due to time variability because of the short rotational period of the sub-AU region of the disk, ~0.8 yr, which we calculated from previous estimations of the inner radius. Interferometric time variability has been reported for other young stellar objects (YSOs), and it was shown to be present in >10% of the accretion disk sample studied by Kobus et al. (2020). The goal of the experiment is to make the best-quality images and models by combining the baselines and (u, v) coverages of the two facilities. We were granted observing time over a three-month period, so we expect to be able to resolve some sub-AU disk motion if present. We have six observations using the CHARA array's MIRC-X instrument, and seven observations using the VLTI array's PIONIER instrument. In this section, we describe the observations at the two facilities.

### 2.1. CHARA

HD 190073 was recorded using the MIRC-X instrument in the summer of 2019 as listed in Table 2. MIRC-X is an IR six-telescope beam combiner at the CHARA Array. The CHARA interferometer is located at Mt. Wilson in California, USA, with an array of six 1 m telescopes disposed in a Y shape, optimizing the imaging capability. The telescopes combine to have baselines of up to 331 meters. The maximum angular resolution is $\lambda/2B$, which corresponds to 0.51 mas in the H band
<table><caption>**Table 2**
CHARA/MIRC-X Observations of HD 190073</caption><thead><tr><th>UT Date</th><th>Configuration</th><th>Band</th><th>no. V<sup>2</sup> (×8)</th><th>No. CP (×8)</th><th>Calibrator(s)</th><th>Diameter (mas)</th><th>Epoch</th></tr></thead><tbody><tr><td>2019-05-09</td><td>E1-W2-W1-S1-E2</td><td>H</td><td>9</td><td>7</td><td>HD 190753</td><td>0.323 ± 0.008</td><td>A<sup>a</sup></td></tr><tr><td>2019-06-05</td><td>E1-W2-W1-S2-S1-E2</td><td>H</td><td>67</td><td>77</td><td>HD 190753<br>HD 191656</td><td>0.323 ± 0.008<br>0.430 ± 0.01</td><td>A</td></tr><tr><td>2019 June 9</td><td>E1-W2-W1-E2</td><td>H</td><td>26</td><td>24</td><td>HD 197551</td><td>0.644 ± 0.02</td><td>A</td></tr><tr><td>2019-07-11</td><td>E1-W2-W1-S2-S1-E2</td><td>H</td><td>66</td><td>74</td><td>HD 190753</td><td>0.323 ± 0.008</td><td>B</td></tr><tr><td>2019-07-12</td><td>E1-W2-W1-S2-E2</td><td>H</td><td>25</td><td>23</td><td>HD 190753</td><td>0.323 ± 0.008</td><td>B</td></tr><tr><td>2019-08-27</td><td>E1-W2-W1-S2-S1-E2</td><td>H</td><td>36</td><td>46</td><td>HD 191656<br>HD 190753</td><td>0.430 ± 0.01<br>0.323 ± 0.008</td><td>B<sup>a</sup></td></tr></tbody></table>

**Note.**

<table><caption>**Table 3**
VLTI/PIONIER Observations of HD 190073</caption><thead><tr><th>UT Date</th><th>Configuration</th><th>Band</th><th>No. V<sup>2</sup> (×6)</th><th>No. CP (×6)</th><th>Calibrator(s)</th><th>Diameter (mas)</th><th>Epoch</th></tr></thead><tbody><tr><td rowspan="2">2019 June 9</td><td rowspan="2">D0-G2-J3-K0 (Medium)</td><td rowspan="2">H</td><td rowspan="2">12</td><td rowspan="2">8</td><td>HD 188385</td><td>0.221 ± 0.01</td><td rowspan="2">A</td></tr><tr><td>HD 189509</td><td>0.247 ± 0.03</td></tr><tr><td rowspan="2">2019 July 10</td><td rowspan="2">D0-G2-J3-K0 (Medium)</td><td rowspan="2">H</td><td rowspan="2">24</td><td rowspan="2">16</td><td>HD 190753</td><td>0.323 ± 0.02</td><td rowspan="2">A</td></tr><tr><td>HD 191840</td><td>0.288 ± 0.01</td></tr><tr><td rowspan="2">2019 July 20</td><td rowspan="2">A0-G1-J2-J3 (Large)</td><td rowspan="2">H</td><td rowspan="2">30</td><td rowspan="2">20</td><td>HD 190753</td><td>0.323 ± 0.02</td><td rowspan="2">A and B</td></tr><tr><td>HD 191840</td><td>0.288 ± 0.01</td></tr><tr><td rowspan="2">2019 July 30</td><td rowspan="2">A0-B2-C1-D0 (Small)</td><td rowspan="2">H</td><td rowspan="2">24</td><td rowspan="2">16</td><td>HD 190753</td><td>0.323 ± 0.02</td><td rowspan="2">A</td></tr><tr><td>HD 191840</td><td>0.288 ± 0.01</td></tr><tr><td rowspan="2">2019 August 4</td><td rowspan="2">A0-B2-C1-D0 (Small)</td><td rowspan="2">H</td><td rowspan="2">12</td><td rowspan="2">8</td><td>HD 190753</td><td>0.221 ± 0.01</td><td rowspan="2">B</td></tr><tr><td>HD 191840</td><td>0.247 ± 0.03</td></tr><tr><td rowspan="2">2019 August 5</td><td rowspan="2">D0-G2-J3-K0 (Medium)</td><td rowspan="2">H</td><td rowspan="2">6</td><td rowspan="2">4</td><td>HD 188385</td><td>0.221 ± 0.01</td><td rowspan="2">B</td></tr><tr><td>HD 189509</td><td>0.247 ± 0.03</td></tr><tr><td rowspan="2">2019 August 6</td><td rowspan="2">D0-G2-J3-K0 (Medium)</td><td rowspan="2">H</td><td rowspan="2">6</td><td rowspan="2">4</td><td>HD 188385</td><td>0.221 ± 0.01</td><td rowspan="2">B</td></tr><tr><td>HD 191840</td><td>0.288 ± 0.01</td></tr></tbody></table>

($\lambda$ = 1.64 µm). MIRC-X uses single-mode fibers to combine coherently the light from the six telescopes simultaneously with an image-plane combination scheme, and typically delivers a visibility precision better than 5% and CP precision better than 1°. The data were collected using the spectral resolution $R$ = 50 mode using a prism dispersive element giving eight spectral channels spread over $\Delta\lambda$ = 0.27 µm. We reduced the data using the MIRC-X reduction pipeline<fnref n="8" /> (version 1.3.3) and calibrated the data using a custom IDL routine.<fnref n="9" /> The MIRC-X data reduction pipeline produces science-ready visibilities and CPs written in OIFITS format (Pauls et al. 2005; Duvert et al. 2017). The raw visibility is calibrated using stars with known diameters. Typical observations execute a standard calibrator-science-calibrator (CAL-SCI-CAL) cycle. Calibrators are usually chosen, using the SearchCal software, to be close to the science target, both in terms of sky position and magnitude, and have a smaller angular diameter so that their visibility on a given baseline is less dependent on the diameter (Chelli et al. 2016).

## 2.2. VLTI

HD 190073 was recorded with the PIONIER instrument on multiple occasions throughout 2019, as listed in Table 3. PIONIER is a four-telescope beam combiner operating in the H band ($\lambda$ = 1.64 µm) at the VLTI. Data were obtained using the

four auxiliary telescopes (ATs) in multiple baseline configurations giving baselines ranging from 11.3 to 140.0 m, allowing for a maximum spatial resolution ($\lambda/B$) of 1.12 mas. Additionally, a prism dispersive element was used giving six spectral channels spread over $\Delta\lambda$ = 0.30 µm. Observations were taken in concatenations of calibrator (CAL) and science (SCI) targets in blocks of either CAL-SCI-CAL or CAL-SCI-CAL-SCI-CAL, to allow for effective monitoring of the transfer function and precise calibration. The data were reduced and calibrated using the standard PNDRS package v3.52 (Le Bouquin et al. 2011).

# 3. Basic Data Presentation

## 3.1. Combining the VLTI and CHARA Data

Using the angular resolution of the instruments and the Keplerian velocity of the inner disk, we calculated the timescale on which each instrument is able to resolve a change in the disk (e.g., rotation). CHARA has a maximum angular resolution of $\lambda/2B$ = 0.51 mas for $\lambda$ = 1.65 µm and $B$ = 331 m. We estimated the Keplerian velocity of the dust rim using $R \sim 1.64$ au from Setterholm et al. (2018) and the star's mass of $M = 6 M_\odot$, which gave us a rotational period of $T$ = 304.5 days or 0.83 yr. That means that it takes 13 days for us to be able to resolve the motion of one resolution element using MIRC-X. Doing the same calculation for the large VLTI baseline, we found that it takes 28 days to resolve a change in one resolution element using PIONIER. With that information, we aimed to

<footnote>a For imaging, this epoch was not included; see Section 5.</footnote>

<footnote n="8">⁸ https://gitlab chara.gsu.edu/lebouquij/mircx_pipeine</footnote>

<footnote n="9">⁹ Contact monnier@umich.edu for the IDL routine.</footnote>
split the data into two epochs, such that the amount of time in
which we expect each epoch to change is minimized. We
grouped the earlier three CHARA nights to make epoch A and
the later three nights to make epoch B. When combining the
VLTI nights, we had to consider the baseline configuration as
well as the dates. The different configurations (small, medium,
and large) correspond to different (u, v) coverages. The small
configuration contributes to the short baselines which show us
the large-scale components of the disk, while the large
configuration, which contributes the long baselines, allows us
to constrain the small-scale components of the inner region. To
get the most coverage, we included the first four nights (2019
June 9, 2019 July 10, 2019 July 20, and 2019 July 30) into
epoch A. Even though the fourth night, 2019 July 30, is more
than 28 days away from the first night, 2019 June 9, we
included it in epoch A because it was the closest night with a
small baseline configuration. The shorter the baseline, the
longer the timescale for resolving a significant change. That
timescale for the small configuration is >1 month, so adding
2019 July 30 to epoch A is justified. The last three nights (2019
August 4, 2019 August 5, and 2019 August 6) went into epoch
B, as well as 2019 July 20, which is shared between the two
epochs as it is the only large baseline configuration we have.
The last column of Table 2 and Table 3 indicates the epoch to
which each night corresponds. We show the combined full (u,
v) coverage, squared visibilities ($V^2$), and CPs from both
instruments for each epoch in Figure 1.

Based on early studies of MIRC (Monnier et al. 2012), we
adopted the same minimum systematic errors. There are two
types of errors associated with $V^2$ measurements. There are
multiplicative errors due to the effects of seeing, and there are
additive errors correcting for biases due to the limitations of the
pipeline. The multiplicative error is 6% while the additive error
is $2 \times 10^{-4}$. These are used to calculate a new minimum error.
If the systematic error is larger than the statistical one, then we
adopt the larger one. The same thing is done for triple
amplitudes but the associated errors are 10% for multiplicative
and $1 \times 10^{-5}$ for additive. CPs have an error floor of $0.1^\circ$.

In Section 3, we discuss our method of reducing the data and
combining the CHARA and VLTI nights into two epochs.

## 3.2. Wavelength Interpolation

Since the wavelength channels vary between MIRC-X and PIONIER, and in fact, even between nights on the same instrument, we fit each the square visibility and CP measurement as a function of wavelength to a quadratic with a linear regression and resampled our data at standard wavelengths in the range 1.5–1.7 µm at 50 nm increments. This wavelength “smoothing” produces evenly spaced visibility measurements in five wavelength channels which we then use for modeling and imaging in the following sections. The amplitude of the error estimates was also fit to a quadratic and resampled at the chosen wavelengths. We then applied the same minimum error calibration as described in Section 3.1.

## 3.3. Squared Visibilities

The combination of CHARA and VLTI allows us to sample
the (u, v) coverage well. Just looking at the non-wavelength-
smoothed visibility curve in Figure 1, we can visually detect a
bounce that was not as clear in previous studies ( Setterholm
et al. 2018). Due to the short time between observation nights,

we are able to avoid some smearing effects that the rotation of
the disk can cause on longer timescales. As mentioned in
Section 3.2 above, due to the mismatch in wavelengths
between nights and instruments, we will be using the
wavelength-smoothed data for modeling in the next section.
An initial qualitative look at the visibility curves can already
point to features that we can predict. The short baseline
measurements from PIONIER appear to be slightly less than
unity at zero baseline, which could indicate the presence of an
over-resolved large-scale halo structure. As the baseline
increases, the visibility drops off rapidly at the short baselines.
The drop-off in visibility tells us about the disk diameter.
Toward longer baselines from MIRC-X, we see a bounce in the
visibility. A bounce in the intermediate baselines indicates the
existence of a cavity in the disk, and the severity of the bounce
characterizes the sharpness of the rim. A shallow bounce points
to a fuzzy rim, while a steep bounce points to a sharper
transition (Dullemond & Monnier 2010). The fits by Lazareff
et al. (2017) of HD 190073 seemed to favor a ring structure like
we expect to see from the visibility curves, while the fits by
 Setterholm et al. (2018) were not able to recover a cavity and
concluded a Gaussian shape instead.

## 3.4. Closure Phases

Closure Phase (CP) measurements from both instruments
and for both epochs can be seen in the bottom two panels of
Figure 1. The majority of the CPs are relatively small, $\lesssim 10^\circ$,
which tells us that the disk is not far from symmetric. On the
shorter baselines, the CP measurements are closest to zero,
which tells us that the larger-scale brightness distribution is
nearly symmetric. However, since the CPs at longer baselines
are nonzero and non-180° we expect the disk to have skewness
on the finer spatial scale, albeit small.

Similar to the visibility plots, we are not showing the error
bars on the measurements. That was done mainly to show a
simplistic visualization of what the data look like. Typical
errors are $\lesssim 5^\circ$.

# 4. Modeling

We split our simple geometric models into two categories in
this section. First, we assume point symmetry and fit the
squared visibilities and not the CPs since the disk is nearly
face-on and we have small CP measurements. The symmetrical
models will allow us to compare our modeling to previous
work, as well as to constrain some physical parameters.
Second, we allow the models to be asymmetric by fitting both
the squared visibilities and CPs and adding more fitting
parameters to test which model best fits the CPs. All modeling
was done using PMOIRED<fnref n="10" /> (Mérand 2022). In order to
combine data from all wavelengths, we modeled the star as a
point source with a fixed power spectrum proportional to $\lambda^{-4}$.
The second component that went into the model is a disk with a
power spectrum $\beta_{\text{disk}} \lambda^{\alpha_{\text{disk}}}$ and we let $\alpha_{\text{disk}}$ and $\beta_{\text{disk}}$ be free
parameters. We added a resolved halo component to the models
to account for any large-scale scattered light which helps us to
constrain the flux contributions better. This halo component
had a similar power spectrum to the disk $\beta_{\text{halo}} \lambda^{\alpha_{\text{halo}}}$ and we let
$\alpha_{\text{halo}}$ and $\beta_{\text{halo}}$ be free parameters. $\alpha_{\text{disk}}$ and $\alpha_{\text{halo}}$ are what we
are going to refer to as spectral slopes in the next section.

<footnote n="10">¹⁰ https://github.com/amerand/PMOIRED</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top row: the (u, v) coverage of the combined data for epochs A (left) and B (right). Data from MIRC-X are shown in black and PIONIER are shown in red. Middle row: unsmoothed combined squared visibility measurements from both instruments. The colors correspond to the H-band wavelengths with the darkest being the longest wavelength. Bottom row: CP measurements from both instruments.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** We introduce a new function to describe the radial profile of the ring. The Double Sigmoid model allows us to constrain the radii and the sharpness of the inner and outer rims separately. The blue curve corresponds to Equation (2), which, by definition, has no flux inside the inner rim. The orange dashed curve corresponds to the α Double Sigmoid from Equation (3), which allows the center to have flux. The black dashed lines indicate the inner $r_1$ and outer $r_2$ radii of the profiles. The red horizontal lines indicate the $\sigma_1$ and $\sigma_2$ scale factors.</figcaption></figure>

### 4.1. Symmetric Modeling

We introduce two geometrical models with different profiles that we refer to as the Doughnut model and the Double Sigmoid model. The Doughnut model is a simple parabolic model that allows us to constrain the flux, inclination, projection angle, which indicates the direction of the major axis, spectral slopes, outer radius of the ring, and thickness. The brightness profile has the following functional form:

$$f(r) = 1 - \left( \frac{r - \bar{r}}{2(r_{\text{max}} - r_{\text{min}})} \right)^2, \quad (1)$$

where $\bar{r}$, $r_{\text{max}}$, and $r_{\text{min}}$ are the mean, maximum, and minimum radii, respectively. The Double Sigmoid model adds more complexity to allow us to fit and scale the inner and outer radii independently, and the profile takes the following functional form:

$$f(r) = \frac{1}{1 + e^{-\frac{(r-R_{\text{in}})}{\sigma_{\text{in}}}}} * \frac{1}{1 + e^{-\frac{(r-R_{\text{out}})}{\sigma_{\text{out}}}}}, \quad (2)$$

where $r_{\text{in}}$ and $r_{\text{out}}$ are the inner and outer radii, respectively, and $\sigma_{\text{in}}$ and $\sigma_{\text{out}}$ are their respective scale factors that set the sharpness of the edges. By definition, the Double Sigmoid model forces the inner cavity to be dark. To investigate the possibility of flux coming from the inner cavity, we modified the Double Sigmoid profile to take the following form:

$$f(r) = \left( \alpha + \frac{1 - \alpha}{1 + e^{-\frac{(r-R_{\text{in}})}{\sigma_{\text{in}}}}} \right) * \frac{1}{1 + e^{-\frac{(r-R_{\text{out}})}{\sigma_{\text{out}}}}}, \quad (3)$$

where the $\alpha$ term is allowed to vary to add flux to the center or go to zero if a completely dark center is the best fit. Figure 2 is an example of what that the intensity profiles might look like using $R_{\text{in}} = 1$, $R_{\text{out}} = 3.4$, $\sigma_{\text{in}} = 0.2$, $\sigma_{\text{out}} = 0.3$, and $\alpha = 0.2$.

The full list of best-fit parameters from the three models are compiled in Table 4. Starting with the Doughnut model, as seen in the top two panels of Figure 3, we notice that nearly all fitting parameters are consistent between epochs A and B. The

disk is nearly face-on with a small inclination angle. The inner radii are very small, indicating a transition that is well beyond the expected dust wall. We found that around half of the flux originates from the disk, 42.5% comes from the star, and the rest is from an over-resolved halo component. We fixed the star's spectrum to be $\propto \lambda^{-4}$, and fit the spectra of the disk and halo as power laws. The spectral slope of the disk was fit to be $-0.10 \pm 0.14$, while the halo's slope was fit to be much steeper $2.08 \pm 1.16$. The goodness of these fits can be seen in Figure 4, where the model is fitted to the squared visibility. The reduced (abbreviated “red”) $\chi^2$ for epoch A was calculated to be $\chi^2_{\text{red,A}} = 1.49$ and for epoch B $\chi^2_{\text{red,B}} = 1.72$, which are also indicated on the plots. Note that the $\chi^2_{\text{red}}$ values are based on $V^2$ fits only. We see in the top two panels of Figure 4 that the model fits the longer baselines well, following the bounce and recovering an inner cavity as we expected, as shown in the top two panels of Figure 3. The shorter baselines, on the other hand, are not fit as well using this model. The Doughnut model relates the outer and inner radii and does not allow us to change the sharpness of the inner and outer rims. Seeing that the shorter baselines were not fit very well using this model points to the larger-scale features not being properly represented by the Doughnut model. To get around this, we introduced the more elaborate Double Sigmoid model as discussed above.

Similar to the Doughnut model, the Double Sigmoid model allows us to fit for the flux, inclination, and projection angle, spectral slopes, as well as inner and outer radii and their sharpness. The middle two panels of Figure 3 show the results of the symmetric Double Sigmoid model for epochs A and B. Similar to the Doughnut model, the disk seems almost face-on, with a clear inner cavity. While the fitted inclination angle for A, $i_A = 10^\circ.1 \pm 2^\circ.9$, matches closely to the angle from the Doughnut model fits, the inclination angle for epoch B, $i_B = 18^\circ.1 \pm 1^\circ.8$, is higher than expected. The inclination of the disk cannot change between epochs so we estimate the inclination to generally be $\lesssim 20^\circ$. The projection angles are significantly lower than the ones from the Doughnut model. That could be attributed to the almost face-on nature of the disk making it harder to constrain the projection angle. The inner
<table><caption>**Table 4**
Best-fit Results from the $V^2$ Modeling</caption><thead><tr><th>Model</th><th>epoch</th><th>incl. (deg)</th><th>PA (deg)</th><th>R<sub>in</sub> (mas)</th><th>R<sub>out</sub> (mas)</th><th>σ<sub>in</sub> (mas)</th><th>σ<sub>out</sub> (mas)</th><th>%f Star</th><th>%f Disk</th><th>%f Halo</th><th>Disk Slope</th><th>Halo Slope</th><th>χ<sup>2</sup><sub>red</sub></th><th>α</th></tr></thead><tbody><tr><td rowspan="2">Doughnut</td><td>A</td><td>10.94±3.0</td><td>42.30±10.0</td><td>0.29±0.01</td><td>2.60±0.02</td><td>...</td><td>...</td><td>43±0.5</td><td>51±0.5</td><td>6±0.5</td><td>-0.10±0.1</td><td>2.08±1.0</td><td>1.49</td><td>...</td></tr><tr><td>B</td><td>10.05±3.0</td><td>55.4±20.0</td><td>0.35±0.01</td><td>2.65±0.02</td><td>...</td><td>...</td><td>42±0.5</td><td>49±0.5</td><td>9±0.5</td><td>-0.32±0.2</td><td>2.10±1.0</td><td>1.72</td><td>...</td></tr><tr><td rowspan="2">Double Sigmoid</td><td>A</td><td>10.14±3.0</td><td>3.4±10.0</td><td>1.26±0.04</td><td>1.25±0.05</td><td>0.26±0.02</td><td>0.43±0.02</td><td>42±0.5</td><td>54±0.5</td><td>4±0.5</td><td>-0.07±0.1</td><td>3.65±2.0</td><td>1.13</td><td>...</td></tr><tr><td>B</td><td>18.1±2.0</td><td>9.58±6.0</td><td>1.28±0.04</td><td>1.27±0.05</td><td>0.27±0.02</td><td>0.45±0.02</td><td>42±0.5</td><td>54±0.5</td><td>4±0.5</td><td>-0.25±0.2</td><td>3.87±2.0</td><td>1.35</td><td>...</td></tr><tr><td rowspan="2">α Double Sigmoid</td><td>A</td><td>9.91±3.0</td><td>0.2±20.0</td><td>1.20±0.1</td><td>1.20±0.2</td><td>0.01±0.07</td><td>0.48±0.02</td><td>41±0.5</td><td>56±0.5</td><td>3±0.5</td><td>0.09±0.1</td><td>4.56±3.0</td><td>1.09</td><td>0.19±0.01</td></tr><tr><td>B</td><td>19.1±2.0</td><td>11.35±5.0</td><td>1.25±0.02</td><td>1.24±0.04</td><td>0.10±0.04</td><td>0.50±0.02</td><td>41±0.5</td><td>55±0.5</td><td>4±0.5</td><td>-0.13±0.2</td><td>4.62±3.0</td><td>1.33</td><td>0.17±0.04</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top row: epoch A and B symmetric Doughnut models. Epoch A: $\chi^2_{red} = 1.49$ and Epoch B: $\chi^2_{red} = 1.72$. Middle row: epoch A and B symmetric Double Sigmoid models. Epoch A: $\chi^2_{red} = 1.13$ and Epoch B: $\chi^2_{red} = 1.35$. Bottom row: epoch A and B symmetric α Double Sigmoid models. Epoch A: $\chi^2_{red} = 1.09$ and Epoch B: $\chi^2_{red} = 1.33$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Global fits of $V^2$ for the Doughnut model (top row), the Double Sigmoid (middle row), and the $\alpha$ Double Sigmoid (bottom row). We see that the Doughnut model does the poorest job of fitting the short baselines while both Double Sigmoid models do a much better job. All three models fit the bounce at the longer baselines well, which was expected since we see an inner cavity in all of them.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Shown in green are the epoch A and B symmetric Doughnut brightness profiles normalized by the total flux. Both the inner and outer edges have some abruptness in the transition. We see a very narrow inner radius and a wide disk overall. Shown in orange are epochs A and B symmetric Double Sigmoid brightness profiles normalized by the total flux. Overplotted in blue is the normalized profile of the α Double Sigmoid model. We see that the outer part of the profiles roughly matches while the inner parts differ. The α Double Sigmoid model shows that 15% of the disk flux comes from the inner region and then shows a sharp inner edge transition.</figcaption></figure>

rim being less diffuse and the outer rim being more diffuse in this model leads to a larger inner radius and smaller outer radius compared to the previous model. It is important to note that the $R_{in}$ and $R_{out}$ parameters do not necessarily reflect the size of the physical inner and outer radii. They are simply the best-fit parameters that reproduce the best brightness profile, and from those profiles, we estimate the physical radii based on the intensity transition. We show the normalized brightness profiles in Figure 5 for both models and both epochs to get a better idea of how they compare. The Double Sigmoid profile is shown in orange while the Doughnut profile is shown in green. Indeed, we see the much smoother inner and outer rims of the Double Sigmoid model compared to the sharp drops of the Doughnut model. The Double Sigmoid model shows a slightly sharper inner edge compared to the outer edge. While the inner radius of the Double Sigmoid model seems larger, and the outer seems smaller, due to the smoother edges, we can see from the profiles that the physical size of the disk from both models does not vary drastically, and the radii are not as different as the fitting parameters might lead us to believe. Other notable differences include the disk flux contribution increasing in the Double Sigmoid model, while the spectral slope decreases indicating an even smaller dependence on color. On the other hand, the halo flux contributions decreased, while the spectral slope grew redder, showing a stronger dependence on color. The reduced $\chi^2$ for epoch A was calculated to be $\chi^2_{red,A} = 1.13$ and for epoch B $\chi^2_{red,B} = 1.35$, which are also indicated on the plots. We see in the middle two panels of Figure 4 that the model fits the longer baselines well, following the bounce and recovering an inner cavity like we expected, as shown in the top middle panels of Figure 3. Unlike the Doughnut model, the Double Sigmoid model fits the shorter baselines much better.

We clearly see that the $\chi^2_{red}$ significantly improves with the Double Sigmoid model, where we see a bigger inner cavity and sharper inner rim.

The α Double Sigmoid models can be seen in the bottom row of Figure 3. The initial global fit showed that adding the α term does improve the fit as evident by $\chi^2_{red}$ decreasing for both epochs, $\chi^2_{red,A} = 1.09$ and $\chi^2_{red,B} = 1.33$, which can be seen on the bottom row of Figure 4. This model tells us that ~15% of the light of the disk emerges from inside the rim. The normalized brightness profile can be seen in blue in Figure 5.

Next, to investigate the inner flux dependence on color, we fit each of the five wavelength channels independently using the same modified Double Sigmoid profile. We see in Figure 6 that all five wavelengths take a similar overall shape with varying flux contributions to the center. We normalized the profiles by the total flux and calculated the fractions of flux contributed by α to the total flux of the disk, which we show as the filled areas under the profile curves. These fractions allow us to calculate the spectrum of the inner disk independently of the spectrum of the outer disk. To calculate the inner disk spectrum, we used the inner flux fractions to calculate how much flux they contribute to the total light from the model that consists of the star, disk, and halo. We also calculated the spectra of the disk excluding the inner emission and the halo. We assumed the star's spectrum to be $\lambda^{-4}$. To test the plausibility of the spectra, we added the flux contributions from the star, disk, inner emission, and halo and compared the total to the previously obtained photometry which we show in Figure 7. The spectral energy distribution (SED), which was built using SEDBYS<fnref n="11" /> (Davies 2021), shows the photometric data in the *J*, *H* (shaded), and *K* bands with the spectra that we

<footnote n="11">¹¹ https://gitlab.com/clairedavies/sedbys</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Normalized flux profile of disk emission. The filled areas correspond to the fraction of flux contributed by the inner emission we define as $\alpha$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** SED in the J, H, (shaded in gray), and K bands shown as the black dots. Overplotted are the fitted spectra from epoch A. The star was assumed to have a $\lambda^{-4}$ spectrum, which is shown in orange. The lighter orange line shows the extension of the slope across the J and K bands. The blue line indicates the flux contribution from the ring, not including the inner emission. The inner emission flux is shown in purple. The over-resolved halo flux contribution is shown in red. The green line indicates the sum of all four contributions and it matches the observed photometry flux (the SED data and references are in Table 5 in Appendix A).</figcaption></figure>

calculated for the four components of HD 190073 overplotted. We see the slope of the total spectrum matches the data closely. The spectrum of the disk emission is rather flat which is more consistent with dust emission rather than the steeper free-free spectrum we would expect from gas emission.

## 4.2. Asymmetric Modeling

So far, we have been assuming that the disk is symmetric and therefore only fitting the squared visibilities. However, as mentioned before, the CPs are nonzero and non-180°, which means that there is asymmetry in the disk. In this section, we use three different simple geometric models to attempt to characterize the asymmetry.

These simple models are no longer valid in describing the geometry of the disk so we have to introduce asymmetry parameters. Azimuthal asymmetry is expected from flared disks

seen at an inclination. Since the Double Sigmoid model performed better than the Doughnut model in the symmetrical case, we used a Double Sigmoid skewed ring model with azimuthal asymmetry (Monnier et al. 2006) as the base of our modeling and added different asymmetries to it. We tested multiple azimuthal asymmetries such as adding a point-source free parameter, adding a higher-order multipole for azimuthal asymmetry, considering an off-center star, and more.

These skewed models allow us to vary the brightness distribution across the disk. In this paper, we present three of the models we tested and show how their CP $\chi^2_{red,CP}$ values compare to give us a better picture of the disk. A full list of the best-fit parameters from these three models can be found compiled in Table 6 in Appendix B. Starting with a simple skewed Double Sigmoid ring, we see in the top two panels of Figure 8 that one side of the disk tends to be brighter. This
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Top row: epoch A and B Skewed Double Sigmoid ring model. Second row: epoch A and B higher-order Skewed Double Sigmoid ring model. Third row: epoch A and B off-center star Double Sigmoid model. Bottom row: epoch A and B model-free images. The white circle in the left bottom corner of the left panel corresponds to the effective beam size $\lambda/(2B_{max}) = 0.55$ mas, where $\lambda = 1.65$ µm and $B_{max} = 300$ m.</figcaption></figure>
model introduces two additional fitting parameters to represent the harmonic azimuthal variation. The amplitude is referred to as “Az Amp₁” in Table 6 in Appendix B, and the projection angle, which is defined with respect to the global projection angle, is denoted as “Az PA₁.” Both parameters have a subscript that represents their order. The amplitude of variation is stronger in epoch A compared to B. The inclination of the flared disk causes the side farther to us to appear brighter but this skewness does not produce a good fit, as indicated by $\chi^2_{red,CP,A} = 7.33$ and $\chi^2_{red,CP,B} = 4.68$. This model is not complex enough to model the asymmetries, so in the next model, we introduced a higher order of asymmetry, $\sin(2\theta)$, to the ring. The higher order allows us to model simple structures that might be in the disk. This adds an extra pair of fitting parameters, Az Amp₂ and Az PA₂, which we list in Table 6 as well. Models of the disk, which can be seen in the second row of Figure 8, produce two distinct bright regions, and the CP $\chi^2_{red,CP}$ has decreased, but not by much. We find $\chi^2_{red,CP,A} = 7.24$ and $\chi^2_{red,CP,B} = 4.49$ for epochs A and B, respectively. The first-order azimuthal amplitudes stay consistent with the previous model. Interestingly enough, the second-order amplitudes have the opposite strengths for the two epochs. Epoch A Az Amp$_{A,2}$ is lower than Az Amp$_{A,1}$, and epoch B Az Amp$_{B,2}$ is higher than Az Amp$_{B,1}$. We are not giving significant physical meanings to these models because they are meant to describe a complex structure with simple parameters. Instead of introducing even higher orders, we went with a less complex approach by testing a skewed Double Sigmoid model with an off-center star, which can be seen in the third row of Figure 8. This model favored a slightly off-center star, which dropped the CP $\chi^2_{red,CP}$ significantly. The simple skewed ring started with a $\chi^2_{red,CP,A} = 7.33$ and $\chi^2_{red,CP,B} = 4.68$, which dropped slightly by adding the second azimuthal parameter to be $\chi^2_{red,CP,A} = 7.24$ and $\chi^2_{red,CP,B} = 4.49$, and then dropped sharply with the addition of the off-center star to be $\chi^2_{red,CP,A} = 1.54$ and $\chi^2_{red,CP,B} = 3.44$. The amplitudes of skewness of the two epochs almost match in this model. Two extra parameters were introduced to characterize how much the star was fit away from the center, an x offset, and a y offset. The x offsets are similar for both epochs and they place the star $0.1 \pm 0.01$ mas west of the center. The y offsets both placed the star $0.1 \pm 0.01$ mas south of the center. While none of the asymmetric models fit to $\chi^2_{red} \sim 1$, we are inclined to favor the off-center star as being at least a plausible component. The reason for the star’s off-center location has not been established yet. One of our leading theories is the presence of a binary in the center, which will need further studies.

## 5. Imaging

Since we were not getting a good fit using simple models, we decided to try imaging in the hopes of unlocking some of the mystery. As mentioned in Section 4, the CP data show nonzero and non-180° measurements that indicate some form of asymmetry in the disk. However, our asymmetric models did not succeed in characterizing the asymmetry. One possibility is an off-center star, but the physical reason behind that is still not established. Since the goal of imaging is to search for finer-scale details, we wanted to fine-tune the data binning in the two epochs. As mentioned in Section 3.1, CHARA can resolve the

motion of one resolution element every 13 days. However, the first night in epoch A is more than 13 days away from the other two nights. Similarly, the last night in epoch B is more than a month away from the first two nights. Therefore, we removed the 2019 May 09 night from epoch A and 2019 August 27 from epoch B.

Images were reconstructed using conventional regularized maximum likelihood and the open-source OITOOLS.jl package<fnref n="12" />. Only the power spectra and CPs were used for the reconstructions. In line with previous mixed modeling and imaging works (e.g., SPARCO in Kluska et al. 2014), the star was modeled as a point source in the center of the field of view, with the spectrum proportional to $\lambda^{-4}$, while the environment was assumed to follow a $\lambda^\alpha$ dependency, with $\alpha$ free to vary. In addition, we modeled potential extended or background emission as a zero visibility component, also following $\lambda^{-4}$. Thus there were three free parameters in addition to the image pixels: $\alpha$ and the flux ratios disk/star and extended emission/star.

Minimization was ensured by the gradient-based VMLM-B (Thiebaut 2002) algorithm from OptimPack, enforcing positivity of the disk intensity. Three other priors were employed in addition to positivity: compactness to regularize large-scale emission, total variation to handle small-scale pixel correlations (see, e.g., Thiebaut & Young 2017, for the mathematical expressions of these), as well as a novel azimuthal regularization to “circularize” YSO disks at medium scales. Given four hyperparameters (the centroid, inclination, and position angles of the disk), this regularizer simply computes the location of concentric elliptical rings and then sums the azimuthal pixel variances along them. Minimizing the regularizer enforces relative smoothness along the rings, unlike total variation. Since the expression can be written as a squared $l_2$ norm of a linear operator on the image, it lends itself well to the gradient-based approach. Hyperpriors for the four hyperparameters were uniform and an initial hyperparameter range was roughly determined from the data by directly minimizing the variance of the de-rotated power spectra. The Nelder-Mead simplex (Nelder & Mead 1965) was used to optimize the four hyperparameters, and new reconstructions were run with different sets of hyperparameters until a global minimum was found. Values for the other regularization hyperparameters were chosen based on quick simulations (copying the (u, v) coverage and signal-to-noise ratio of the actual data and reconstructed disks), within the expected range (Renard et al. 2011).

We see a face-on thin ring with some inner emission, similar to the $\alpha$ Double Sigmoid model, and an off-center star. We clearly recover a cavity in the center where we see some fine-scale structure that seems to change. The last row of Figure 8 shows the images that we produced without priors from modeling. The bright spots in the ring seem to correspond with the location of the bright areas seen in a combination of the skewed models. Note that in the images, the point source representing the star is shown in the center of the image while the disk’s center is off-center. The final reduced $\chi^2$ for the images is $\chi^2_{v2} \simeq 1.09$ and $\chi^2_{t3phi} \simeq 0.72$ for epoch A and $\chi^2_{v2} \simeq 1.01$ and $\chi^2_{t3phi} \simeq 1.19$ for epoch B. A closer look at the images can be seen in Figure 9. We explored if the changes might be due to rotation. We performed a qualitative analysis

<footnote n="12">¹² https://github.com/fabienbaron/OITOOLS.jl</footnote>
Epoch A (left), Epoch B (right)

<figure><img src="image_11.png" /><figcaption>Figure 9. Epochs A and B images. The dashed crosshairs show the center of the image. While the star is in the center, the disk is not which shows that this matches our off-center star model. The green arrow indicates a 27° rotation. The spot starts at the bottom of the green arrow in epoch A and moves to the top in epoch B. The white circle in the left bottom corner of the left panel corresponds to the effective beam size $\lambda/(2B_{max}) = 0.55$ mas, where $\lambda = 1.65$ µm and $B_{max} = 300$ m.</figcaption></figure>

of the apparent rotation between epochs A and B by manually rotating the disk to match the bright feature in the bottom right quadrant in epoch A to the bright spot on the middle right edge of the ring in epoch B, which is indicated in Figure 9. The arrows represent a 27° counterclockwise rotation of the disk. Through modeling, we found that the $\alpha$ Double Sigmoid intensity peaked at a radius of $R_A = 1.23 \pm 0.02$ mas and $R_B = 1.5 \pm 0.02$ mas. The Keplerian rotational period at this smaller radius that we average to be $R = 1.35$ mas is 182.7 days. The two epochs are 32 days apart which corresponds to a 63° rotation based on the calculated period. However, this motion is not consistent with the 27° rotation we estimate from the images, which is more than twice as slow as Keplerian rotation. The origin of such features is not well-known but we speculate that the motion could be caused by interactions in the outer disk versus an object embedded in the inner disk, for example. This structure could be related to planet formation, instabilities in the accretion flow, or even magnetic fields. This change could also be an artifact of the (u, v) coverage and not a physical process at all, so additional epochs are needed to confirm the origin of these changes.

## 6. Discussion

The symmetric $V^2$-only models revealed a clear inner cavity in the disk, along with rough radii estimates. We are taking the radius at which the intensity peaks in the Double Sigmoid and $\alpha$ Double Sigmoid models to be representative of the evaporation front. Since the inner radius is set by the dust evaporation temperature, we can use the average radius from the models to estimate the temperature of the dust wall using the following equation from Dullemond & Monnier (2010):

$$T_{\text{dust}} = T_{*} \frac{1}{\epsilon^{1/4}} \sqrt{\frac{R_{*}}{2R}} \quad (4)$$

This equation treats the inner rim as an optically thick wall that radiates like a blackbody, assumes that the gas inside the rim is transparent, and includes the “backwarming” of dust grains. The $\epsilon$ term is the ratio of the effectiveness of emission at the wavelength at which the dust radiates away its heat and absorption at stellar wavelengths. Assuming large dust grains, we can set $\epsilon = 1$ and solve for the temperature of the inner edge using values from Table 1 and $R = 1.35 \pm 0.2$ mas for the radius of the inner rim. We calculate an inner rim temperature of $T_{\text{dust}} = 1367.13 \pm 107$ K, which is a reasonable estimate. If the dust grains have a larger $\epsilon$, meaning that they cool more efficiently, then the grains can exist closer to the star. Our modified $\alpha$ Double Sigmoid model showed that ~15% of the disk flux comes from the inside of the ring. The sub-AU inner emission could be due to such grains with $\epsilon > 1$, gas emission, or grains that can survive at higher temperatures.

We have nonzero CPs which indicates that we have some asymmetry in the ring. The CPs are relatively small so the asymmetry is not huge. From our fits, it seems that the asymmetry is best fit by including an off-center star. While not certain about the physical meaning of that fit, one speculation is that we are looking at a binary star system.

One reassuring confirmation that the model is probable, is that the imaging produced an off-center star as well. The imaging also confirmed the cavity in the center of the ring and showed new evidence of small-scale structures with rotation slightly slower than Keplerian at that radius. There is already evidence of this temporal variation on the inner few AU scale in protoplanetary disks. Kobus et al. (2020) showed evidence of interferometric temporal variation occurring in >10% of their sample of 68 accretion disks. The asymmetric variations could be attributed to companions, forming exoplanets, asymmetries in the dust density distribution, asymmetric illuminations of the circumstellar material as a result of stellar
spots or obscuring, or artifacts. There are a few inner disk dynamics that have been observed that could help unlock the mystery of the variations we see in our disk. Some possible inner disk dynamics include self-shadowing by the inner disk on the outer disk (Garufi et al. 2022), dippers caused by turbulent accretion (Alencar et al. 2010), or even a small planet actively accreting close to the inner rim with accompanying spirals, vortices, or misalignment (see Benisty et al. 2017; Marr & Dong 2022). Characterizing these structures is going to be essential to furthering our understanding of planet formation.

## 7. Conclusions

We tested multiple simple geometric models, both symmetric and asymmetric, in hopes of getting a better glance at the young Herbig Be star HD 190073's accretion disk. With more sensitive data, better (u, v) coverage, and a shorter time span to avoid smearing effects, we were able to model and image the disk with unprecedented accuracy. CHARA's long baselines probe the finer structures of the disk, while the short baselines of VLTI look at the larger-scale structures. For a rapidly changing object like HD 190073, it is important to get near-simultaneous observations with CHARA and VLTI to get as close to an instant snapshot of the disk as possible. While there is still some mystery surrounding HD 190073, we were able to produce convincing evidence of a ring-like structure with an inner cavity, and some evidence of changes. The modified α Double Sigmoid model showed that the inner and outer rims are relatively smooth, with the inner rim having a sharper profile due to the dust destruction front. We also see that the disk is almost face-on with inclination angles of $i \lesssim 20^{\circ}$. The modeling revealed the fractions of flux contributed by each component. The star contributes 41%, the halo 4.6%, and the disk ~56%, with ~15% of that disk flux coming from inside the inner rim. We found that the best simple model to fit the CP measurements is a skewed Double Sigmoid ring with an off-center star, which was backed up by the imaging. The imaging revealed possible evidence of rotating sub-AU features. They seem to rotate 27° counterclockwise in the span of the 32 days between epochs A and B, which would correspond to a rotational period that is more than two times slower than Keplerian. This change could be caused by dynamics in the outer disk, or it could be an artifact.

We plan on continuing the study of HD 190073 in the future. We have new MYSTIC (Monnier et al. 2018) K-band data that we will be adding to get a wider look at the disk. We also plan on quantitatively analyzing the temporal rotation to get a better estimate of the rotation angle and to measure the correlation between the two epochs. The correlation will be able to confirm further that the rotation is not an artifact of the (u, v) coverage. Future data will also reveal if the star is off-center due to an inner binary.

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. MIRC-X received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant No. 639889). J.D.M. acknowledges funding for the development of MIRC-X (NASA-XRP NNX16AD43G and NSF-AST 1909165) and MYSTIC (NSF-ATI 1506540 and NSF-AST 1909165). This research has made use of the Jean-Marie Mariotti Center Aspro and SearchCal services. S.K., N. A., and C.L.D. acknowledge support from an ERC Starting Grant ("ImagePlanetFormDiscs," grant agreement No. 639889), ERC Consolidator Grant ("GAIA-BIFROST", grant agreement No. 101003096), and STFC Consolidated Grant (ST/V000721/1). A.L. received funding from STFC student-ship No. 630008203. Observing travel support was provided by STFC PATT grant ST/S005293/1. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programs 0101.C-0896(A) and 0103.C-0915(A,B,C). E.A.R. and J.D.M. acknowledge support from NSF-AST 1830728.

**Facilities:** CHARA, VLTI.

**Software:** Astropy (Astropy Collaboration et al. 2013, 2018, 2022), PMOIRED (Mérand 2022), SEDBYS (Davies 2021), PNDRS (v3.52; Le Bouquin et al. 2011), OITOOLS.jl (https://github.com/fabienbaron/OITOOLS.jl).

## Appendix A
SED Data Table

The photometry data used in the spectral energy distribution in Figure 7 are displayed in Table 5.

<table><caption>**Table 5**
HD 190073 Photometry Data</caption><thead><tr><th>Wavelength<br>μm</th><th>λF<sub>λ</sub><br>10<sup>-13</sup> W m<sup>-2</sup></th><th>Date</th><th>Reference</th></tr></thead><tbody><tr><td>1.25</td><td>53.036 ± 3.908</td><td>2006 Jun 2</td><td>Tannirkulam et al. (2008)</td></tr><tr><td>1.25</td><td>48.370</td><td>2000 Jun 2</td><td>Lazareff et al. (2017)</td></tr><tr><td>1.25</td><td>50.676 ± 0.887</td><td>2000 Jun 4</td><td>Cutri et al. (2003)</td></tr><tr><td>1.60</td><td>42.956 ± 2.769</td><td>2006 Jun 2</td><td>Tannirkulam et al. (2008)</td></tr><tr><td>1.60</td><td>43.353</td><td>2000 Jun 2</td><td>Lazareff et al. (2017)</td></tr><tr><td>1.65</td><td>40.817 ± 0.639</td><td>2000 Jun 4</td><td>Cutri et al. (2003)</td></tr><tr><td>2.15</td><td>42.304 ± 1.052</td><td>2000 Jun 4</td><td>Cutri et al. (2003)</td></tr><tr><td>2.18</td><td>42.629</td><td>2000 Jun 2</td><td>Lazareff et al. (2017)</td></tr><tr><td>2.18</td><td>43.421 ± 3.199</td><td>2006 Jun 2</td><td>Tannirkulam et al. (2008)</td></tr></tbody></table>
## Appendix B
Asymmetric Modeling Fitting Parameters

Best-fit parameters from the asymmetrical modeling can be seen compiled in Table 6.

<table><caption>**Table 6**
Fitting Parameters Results from Asymmetric Models</caption><thead><tr><th rowspan="2">Model<br>Epoch</th><th colspan="2">Skewed Ring</th><th colspan="2">2nd Az Parameter</th><th colspan="2">Off-center</th></tr><tr><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>incl. (deg)</td><td>14.18 ± 3.0</td><td>16.4 ± 3.0</td><td>25.2 ± 2.0</td><td>30.0 ± 2.0</td><td>10.7 ± 3.0</td><td>19.7 ± 2.0</td></tr><tr><td>PA (deg)</td><td>101.7 ± 10.0</td><td>162.1 ± 10.0</td><td>118.2 ± 5.0</td><td>157.9 ± 5.0</td><td>3 ± 10.0</td><td>157.7 ± 6.0</td></tr><tr><td>R<sub>in</sub> (mas)</td><td>1.38 ± 0.07</td><td>1.27 ± 0.06</td><td>1.39 ± 0.07</td><td>1.26 ± 0.05</td><td>1.16 ± 0.02</td><td>1.11 ± 0.02</td></tr><tr><td>R<sub>out</sub> (mas)</td><td>1.37 ± 0.08</td><td>1.26 ± 0.06</td><td>1.39 ± 0.07</td><td>1.25 ± 0.06</td><td>1.14 ± 0.04</td><td>1.10 ± 0.04</td></tr><tr><td>σ<sub>in</sub> (mas)</td><td>0.30 ± 0.03</td><td>0.28 ± 0.03</td><td>0.31 ± 0.03</td><td>0.28 ± 0.03</td><td>0.18 ± 0.01</td><td>0.17 ± 0.02</td></tr><tr><td>σ<sub>out</sub> (mas)</td><td>0.40 ± 0.02</td><td>0.46 ± 0.03</td><td>0.43 ± 0.03</td><td>0.52 ± 0.03</td><td>0.44 ± 0.01</td><td>0.53 ± 0.02</td></tr><tr><td>% f<sub>star</sub></td><td>42 ± 0.5</td><td>42 ± 0.5</td><td>42 ± 0.5</td><td>42 ± 0.5</td><td>42 ± 0.5</td><td>42 ± 0.5</td></tr><tr><td>% f<sub>disk</sub></td><td>54 ± 0.5</td><td>54 ± 0.5</td><td>54 ± 0.5</td><td>55 ± 0.5</td><td>54 ± 0.5</td><td>54 ± 0.5</td></tr><tr><td>% f<sub>halo</sub></td><td>4 ± 0.5</td><td>4 ± 0.5</td><td>4 ± 0.5</td><td>3 ± 0.5</td><td>4 ± 0.5</td><td>4 ± 0.5</td></tr><tr><td>disk slope</td><td>-0.14 ± 0.2</td><td>-0.27 ± 0.2</td><td>-0.12 ± 0.2</td><td>-0.22 ± 0.2</td><td>-0.13 ± 0.1</td><td>-0.20 ± 0.2</td></tr><tr><td>halo slope</td><td>3.60 ± 3.0</td><td>4.1 ± 3.0</td><td>3.58 ± 3.0</td><td>4.48 ± 4.0</td><td>4.10 ± 2.0</td><td>3.4 ± 3.0</td></tr><tr><td>Az Amp<sub>1</sub></td><td>0.14 ± 0.01</td><td>0.04 ± 0.01</td><td>0.14 ± 0.01</td><td>0.05 ± 0.01</td><td>0.08 ± 0.02</td><td>0.07 ± 0.01</td></tr><tr><td>Az PA<sub>1</sub> (deg)</td><td>-153.2 ± 10.0</td><td>57.1 ± 20.0</td><td>-170.9 ± 6.0</td><td>38.6 ± 10.0</td><td>-170.4 ± 20.0</td><td>39.57 ± 9.0</td></tr><tr><td>Az Amp<sub>2</sub></td><td>...</td><td>...</td><td>0.09 ± 0.02</td><td>0.15 ± 0.02</td><td>...</td><td>...</td></tr><tr><td>Az PA<sub>2</sub> (deg)</td><td>...</td><td>...</td><td>80.05 ± 4.0</td><td>-92.4 ± 3.0</td><td>...</td><td>...</td></tr><tr><td>x offset (mas)</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.12 ± 0.01</td><td>0.09 ± 0.01</td></tr><tr><td>y offset (mas)</td><td>...</td><td>...</td><td>...</td><td>...</td><td>-0.15 ± 0.01</td><td>-0.09 ± 0.01</td></tr><tr><td>Overall χ<sup>2</sup><sub>red</sub></td><td>2.10</td><td>2.91</td><td>2.08</td><td>2.83</td><td>1.35</td><td>2.28</td></tr><tr><td>V<sup>2</sup> χ<sup>2</sup><sub>red</sub></td><td>1.93</td><td>1.74</td><td>1.86</td><td>1.65</td><td>1.96</td><td>1.88</td></tr><tr><td>CP χ<sup>2</sup><sub>red</sub></td><td>7.33</td><td>4.68</td><td>7.24</td><td>4.49</td><td>1.54</td><td>3.44</td></tr></tbody></table>
ORCID iDs

Nour Ibrahim https://orcid.org/0000-0002-1788-9366

John D. Monnier https://orcid.org/0000-0002-3380-3307

Stefan Kraus https://orcid.org/0000-0001-6017-8773

Jean-Baptiste Le Bouquin https://orcid.org/0000-0002-0493-4674

Narsireddy Anugu https://orcid.org/0000-0002-2208-6541

Fabien Baron https://orcid.org/0000-0002-8376-8941

Theo Ten Brummelaar https://orcid.org/0000-0002-0114-7915

Claire L. Davies https://orcid.org/0000-0001-9764-2357

Jacob Ennis https://orcid.org/0000-0002-1575-4310

Tyler Gardner https://orcid.org/0000-0002-3003-3183

Aaron Labdon https://orcid.org/0000-0001-8837-7045

Cyprien Lanthermann https://orcid.org/0000-0001-9745-5834

Antoine Mérand https://orcid.org/0000-0003-2125-0183

Evan Rich https://orcid.org/0000-0002-1779-8181

Gail H. Schaefer https://orcid.org/0000-0001-5415-9189

Benjamin R. Setterholm https://orcid.org/0000-0001-5980-0246

References

Alecian, E., Neiner, C., Mathis, S., et al. 2013, A&A, 549, L8

Alencar, S. H. P., Teixeira, P. S., Guimarães, M. M., et al. 2010, A&A, 519, A88

Alonso-Albi, T., Fuente, A., Bachiller, R., et al. 2009, A&A, 497, 117

Anugu, N., Le Bouquin, J.-B., Monnier, J. D., et al. 2020, AJ, 160, 158

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42

Cacciapuoti, L., Kostov, V. B., Kuchner, M., et al. 2022, MNRAS, 513, 102

Catala, C., Alecian, E., Donati, J.-F., et al. 2007, A&A, 462, 293

Chelli, A., Duvert, G., Bourgès, L., et al. 2016, A&A, 589, A112

Cox, A. N. 2000, Allen's Astrophysical Quantities (4th ed.; New York: AIP Press, Springer)

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, 2246, 0

Davies, C. L. 2021, SoftX, 14, 100687

Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957

Dullemond, C. P., & Monnier, J. D. 2010, ARA&A, 48, 205

Duvert, G., Young, J., & Hummel, C. A. 2017, A&A, 597, A8

Garufi, A., Dominik, C., Ginski, C., et al. 2022, A&A, 658, A137

Guzmán-Díaz, J., Mendigutía, I., Montesinos, B., et al. 2021, A&A, 650, A182

Hillenbrand, L. A., Strom, S. E., Vrba, F. J., & Keene, J. 1992, ApJ, 397, 613

Kluska, J., Malbet, F., Berger, J.-P., et al. 2014, A&A, 564, A80

Kobus, J., Wolf, S., Ratzka, T., & Brunngräber, R. 2020, A&A, 642, A104

Lazareff, B., Berger, J. P., Kluska, J., et al. 2017, A&A, 599, A85

Le Bouquin, J.-B., Berger, J.-P., Lazareff, B., et al. 2011, A&A, 535, A67

Marr, M., & Dong, R. 2022, ApJ, 930, 80

Mérand, A. 2022, Proc. SPIE, 12183, 121831N

Millan-Gabet, R., Schloerb, F. P., & Traub, W. A. 2001, ApJ, 546, 358

Millan-Gabet, R., Schloerb, F. P., Traub, W. A., et al. 1999, ApJL, 513, L131

Monnier, J. D., & Millan-Gabet, R. 2002, ApJ, 579, 694

Monnier, J. D., Berger, J.-P., Millan-Gabet, R., et al. 2006, ApJ, 647, 444

Monnier, J. D., Che, X., Zhao, M., et al. 2012, ApJL, 761, L3

Monnier, J. D., Le Bouquin, J.-B., Anugu, N., et al. 2018, Proc. SPIE, 10701, 1070122

Natta, A., Prusti, T., Neri, R., et al. 2001, A&A, 371, 186

Nelder, J. A., & Mead, R. 1965, CompJ, 7, 308

Pauls, T. A., Young, J. S., Cotton, W. D., & Monnier, J. D. 2005, PASP, 117, 1255

Rein, H. 2012, arXiv:1211.7121

Rich, E. A., Monnier, J. D., Aarnio, A., et al. 2022, AJ, 164, 109

Renard, S., Thiébaut, E., & Malbet, F. 2011, A&A, 533, A64

Setterholm, B. R., Monnier, J. D., Davies, C. L., et al. 2018, ApJ, 869, 164

Tannirkulam, A., Monnier, J. D., Harries, T. J., et al. 2008, ApJ, 689, 513

ten Brummelaar, T. A., McAlister, H. A., Ridgway, S. T., et al. 2005, ApJ, 628, 453

Thiébaut, É. 2002, Proc. SPIE, 4847, 174

Thiébaut, É., & Young, J. 2017, JOSAA, 34, 904

Tuthill, P. G., Monnier, J. D., & Danchi, W. C. 2001, Natur, 409, 1012