© The Authors 2024

1 Introduction

In the last decade, our knowledge of protoplanetary disks around young stars has grown considerably thanks to the drastic improvement of observing facilities. It is nowadays established that such disks show different substructures across the optical to submillimeter wavelength range, such as rings, gaps, spiral arms, vortices, warps, and shadows on scales of tens of au (Huang et al. (2018), Garufi et al. (2018) and references therein). Thanks to long-baseline infrared interferometry, the morphology of inner disks at scales of less than 1 au has been revealed (e.g., Monnier et al. 2005; Eisner et al. 2014; Lazareff et al. 2017; GRAVITY Collaboration 2019; Kluska et al. 2020). Temporal photometric variability is a common property of YSOs (e.g., Kóspál et al. 2012; Rice et al. 2015; Wolk et al. 2018; Guarcello et al. 2019; Robinson & Espaillat 2019) and may originate, for instance, from the presence of cool or hot spots on the stellar surface, variable accretion, or changes in the inner disk structure leading to partial occultation of the central star and variable dimming of the system as in the case of “dippers” or UX-Ori type objects. The question of time-variable morphology of the inner disk has been tackled for a handful of sources (e.g., Kluska et al. 2016; Kobus et al. 2020; GRAVITY Collaboration 2021b). The scarcity of such studies is mainly due to the fact that the long-period high-cadence observations of the same object required for such analyses are seldom available. In this context, observations with a large temporal baseline using the PIONIER (Le Bouquin et al. 2011) and GRAVITY instruments (GRAVITY Collaboration 2017) at the VLTI provide a unique opportunity to probe the origin of the variability in the brightness distribution of the innermost regions of YSOs.

In the present paper, we study HD 98922, a B9Ve/A2III Herbig star (Hales et al. 2014; Caratti o Garatti et al. 2015) characterized by a spectral energy distribution (SED) with a high near-infrared (NIR) excess and a low far-infrared (FIR) excess (Garufi et al. 2022). Classifications of the star and estimates of its parameters have shown discrepancies over the years. Indeed, the distance estimates were drastically divergent before the Gaia era, with values going from ~450 pc (Caratti o Garatti et al. 2015) to ~1150 pc (van Leeuwen 2007), which affected the derivation of parameters such as luminosity, mass, and age. For instance, Lee et al. (2016) suggested classification as a post-main sequence giant, but accurate Gaia parallax information supports the classification of HD 98922 as a Herbig Be star (Vioque et al. 2018; Arun et al. 2019). Using VLT/SPHERE, Garufi et al. (2022) set a lower limit of at least 200 au in radius for the physical extent of the dust disk. However, comparison with an ALMA image suggests an even larger radius of ≲500au (Garufi et al. 2022). At scales of smaller than 10 au, Menu et al. (2015) constrained the N-band-emitting dust disk radius to ~7.2 au using VLTI/MIDI interferometric observations.

Furthermore, the system shows a CO-rich circumstellar disk, which Hales et al. (2014) suggest is geometrically flat with an inner radius of ~1 au, extending to ~200 au. Other authors instead suggested a system with a flared CO disk with an inner radius of ~5 au and a flattened dust disk (van der Plas et al. 2015). This transitional disk scenario was suggested because of the relatively strong polycyclic aromatic hydrocarbon (PAH) emission of HD 98922 and its similarities with two other objects, HD 101412 and HD 95881, for which such a disk structure was already proposed (Fedele et al. 2008; Verhoeff et al. 2010). At smaller spatial scales, a rotating gaseous disk inside the dust sublimation radius was proposed to explain the [OI] emission-line profiles (Acke et al. 2005), while a wind or outflow was suggested to explain the P Cygni profiles in the Hα, Si II, and He I lines (Grady et al. 1996; Oudmaijer et al. 2011). Finally, the Brγ emission line was found to arise from a compact region of ~0.65 au in radius, possibly tracing magnetospheric accretion (Kraus et al. 2008) or a disk wind (Caratti o Garatti et al. 2015).

Previous interferometric works (Lazareff et al. 2017; GRAVITY Collaboration 2019; Kluska et al. 2020) revealed the noncentrosymmetric brightness distribution of the inner disk of HD 98922. Here, we exploit a VLTI/GRAVITY and VLTI/PIONIER multi-epoch data set to study the temporal morphological variability of the innermost circumstellar environment of this source, looking at the NIR continuum emission and the hot hydrogen gas. In Sect. 2, we present the observations and describe the interferometric data. In Sect. 3, we present the methodology used to analyze the continuum data and our results. In Sect. 4, we present the methodology used to analyze the Brγ-emitting gas data and the ensuing results. In Sect. 5, we present a new radiative transfer model for the source. In Sect. 6, we discuss some possible interpretations of our results in detail, and finally, in Sect. 7, we summarize our main findings.

2 Observations and data reduction

2.1 The source

HD 98922 is a ~6 M_⊙ pre-main sequence star of less than 1 Myr of age that shows a high accretion rate. We adopt the distance of 650.9 ± 8.8 pc from Gaia DR3 (Gaia Collaboration 2022), which we assume throughout the paper. The most up-to-date stellar parameters are listed in Table 1. HD 98922 is a group II source (Juhász et al. 2010) in the Meeus classification (Meeus et al. 2001), which is typically associated with a flat disk morphology. The K band emission was measured to be inside ~1.5 au (GRAVITY Collaboration 2019) and the H band emission inside ~1.2 au (Lazareff et al. 2017).

2.2 Observations

HD 98922 was observed with VLTI/PIONIER (Le Bouquin et al. 2011) using the four 1.8 m Auxiliary Telescopes (ATs) in 21 different epochs between 2011 and 2016. Data were obtained using small- to large-baseline configurations for different epochs. The data consist of low-spectral-resolution (R ≈ 40) interferometric observables in the H band. The observations span a spatial frequency range between about 5 Mλ and 90 Mλ with a maximal angular resolution of λ/2B ~1.25 mas for the longest baseline of 138.7 m, which corresponds to 0.81 au at 650.9 pc. In total, 45 files were acquired and four files were discarded due to bad weather conditions. The description of the data per epoch can be found in Table A.1 along with the observation logs.

HD 98922 was observed at 13 different epochs between 2017 and 2022 using GRAVITY, with the so-called astrometric, large-, medium-, and small-baseline configurations. The data consist of high-spectral-resolution (R ≈ 4000) observables recorded by the science channel (SC) detector across the K-band with individual integration times of 30 s, as well as low-spectral-resolution (R ≈ 20) observables recorded using the fringe tracker (FT) detector. The spatial frequency ranges between about 5 Mλ and 65 Mλ with a maximal angular resolution of λ/2B ~ 1.72 mas for the longest baseline of 129 m, which corresponds to 1.12 au at the distance of HD 98922. Each observation block corresponds to 5 min of observing time on the object. In total, 97 files were acquired. Detailed information on the GRAVITY observations is given Table A.2. In summary, the PIONIER and GRAVITY dataset spans an 11 yr period from 2011 to 2022.

2.3 Data reduction

The PIONIER data are archival reduced and calibrated data retrieved from the JMMC Optical interferometry DataBase¹. The errors bars were calculated using the reduction and calibration pipeline. These range from 0.25° to 10.1° for the closure phases and from 0.001 to 0.18 for the squared visibilities, depending on the atmospheric conditions.

The GRAVITY data were reduced and calibrated using the GRAVITY data reduction software (Lapeyrere et al. 2014). For the low-resolution FT data, we discarded the first spectral channel, which is typically affected by the metrology laser operating at 1.908 µm. Following GRAVITY Collaboration (2019), we applied a floor value on the error bars of 2% for the squared visibilities and 1° on the closure phases as the error bars computed by the pipeline might be underestimated.

For the high-resolution science (SC) data, the pipeline produces single files containing the spectrum, the calibrated visibilities, and the calibrated differential phases. The observables for one epoch result from the averaging of N single files as listed in the fourth column of Table A.2. Each averaged file is wavelength-calibrated using the position of the telluric absorption lines bracketing the Brγ emission line and then corrected for the radial velocity of the star and the Earth’s motion with respect to the LSR. The wavelength calibration is discussed in more detail in Appendix C. The error bars are calculated by the reduction pipeline for the single files and are propagated through the averaging process. These range from 0.3% to 1% for the spectrum, from 0.001 to 0.01 for the squared visibilities, and from 0.5° to 2° for the differential phases, depending on the observation epoch.

3 Results derived from the continuum interferometric data

The reduced data are illustrated in Fig. 1 for a few selected epochs and are shown in full in Figs. A.1 and A.2. The plots show the calibrated squared visibilities, the closure phases, and the (u,v) coverage. PIONIER and GRAVITY spatially resolves the H band and K band continuum emission of HD 98922 at all baselines and epochs, with visibilities ranging between 0 and 0.8. Clear closure-phase signals up to 20° with PIONIER and 40° with GRAVITY are detected for all epochs. The data therefore clearly suggest asymmetries in the brightness distribution of HD 98922 at spatial scales probed by the VLTI. Importantly, for the configurations with a similar (u,v) coverage, we detect significant variations in the closure phases across the epochs. On the contrary, we observe that the change in the corresponding visibilities is not very strong, typically below V² ~ 0.05.

3.1 Modeling methodology

The observational results show that, for the epochs using VLTI configurations with a similar (u,v) plane coverage, we see clear variations in the interferometric quantities and in particular in the closure phase signal. This indicates a noncentrosymmetric brightness distribution that is temporally variable. We adopt a classical approach in long-baseline interferometry based on the parametric fit of geometrical models to the squared visibilities and closure phase signals. As HD 98922 is known from previous works to host a circumstellar disk (e.g., Kluska et al. 2020), we focus our methodology on the analysis of a disk-like parametric model.

3.1.1 Choice of the model

A geometrical disk model with an azimuthal modulation is one possible solution to account for the asymmetric brightness distribution revealed by the nonzero closure phases (Lazareff et al. 2017); we adopt such a model here. We used chromatic geometric models that consist of a point-like central star, a scattered light component (called halo), and a circumstellar environment to fit the continuum interferometric quantities. The star is assumed to be unresolved and the halo to be fully resolved with a visibility of zero. The circumstellar emission is modeled through an azimuthally modulated wireframe (Lazareff et al. 2017) with a radial brightness distribution given by $$ℱ(r)=\frac{1}{2\pi }\delta \left(r-{\text{a}}_r\right)\cdot \left(1+{\displaystyle \sum _{j=1}^m\left(c_j\cos j\varphi +s_j\sin j\varphi \right)}\right),$$(1)

where ϕ is the polar angle. The order of the azimuthal modulation is taken in our case to be m = 1. The wireframe is convolved by an ellipsoid kernel that regulates the width of the ring-like emission, and whose visibility is given by Eq. (9) of Lazareff et al. (2017). Hence, the model can describe from infinitesi-mally thin rings to very wide rings tending to ellipsoids. It is described by ten parameters: the star flux contribution F_s; the halo flux contribution F_h; the spectral index of the circumstellar disk k_c; the radius of the wireframe a_r and of the Kernel a_k; the inclination i; the position angle (PA); the weighted contribution of a Gaussian or Lorentzian distribution fLor; and the azimuthal modulation parameters c₁ and s₁. More details on the relation between the geometrical, physical, and fitted parameters is provided in Sect. 3.6 of Lazareff et al. (2017).

The total complex visibility of the system at spatial frequencies (u, v) and at wavelength λ is therefore described by a linear combination of the three components as $$V(u,v,\lambda )=\frac{F_{\text{s}}{\left(\lambda /{\lambda }_0\right)}^{k_{\text{s}}}+F_{\text{c}}{\left(\lambda /{\lambda }_0\right)}^{k_{\text{c}}}{V}_{\text{c}}(u,v)}{\left(F_{\text{s}}+F_{\text{h}}\right){\left(\lambda /{\lambda }_0\right)}^{k_{\text{s}}}+F_{\text{c}}{\left(\lambda /{\lambda }_0\right)}^{k_{\text{c}}}},$$(2)

where V_c is the complex visibility of the circumstellar environment, and F_s, F_h, and F_c are the specific fractional flux contributions of the star, of the halo, and of the circum-stellar environment, respectively, at λ₀ = λ_K,0 =2.15 µm (the wavelength of the central spectral channel of the GRAVITY FT), or at λ₀ = λ_H,0 =1.68 µm (for PIONIER). The parameter k_s = d log F_λ,s/d log λ is the spectral index of the star assumed to be a black body at T_eff=10500 K. This translates into a spectral index of k_s =−3.645 at λ_K,0 and −3.523 at λ_H,0. The parameter k_c = d log F_λ,c/d log λ is the spectral index of the circumstellar environment.

At this point, we must reiterate an important convention of the azimuthal modulation relevant for the correct interpretation of our results. The azimuthal modulation, which is parametrized with the variables c1 and s1, has a PA of the peak emission given by the argument of the complex number c1 + j.s1. The origin of the angular position of the azimuthal modulation is the PA of the disk measured from north to east, with east to the left. As an example, an azimuthal modulation described by c1=1 and s1=1 in a disk with a PA = 45° will show a visual rendering where the azimuthal modulation peak emission appears at 90° towards east.

The model fitting consists in an initial minimization procedure with scipy.optimize.minimize using a sequential least squares programming method to obtain an initial guess of the free parameters, followed by a procedure based on a Markov chain Monte Carlo (MCMC; Foreman-Mackey et al. 2013) numerical approach, which is robust against trapping in local minima. We also report, when applicable, the uncertainty on the reduced chi-square ${\chi }_{\text{r}}^2$ by computing the error on the mean of the stochastic variable T_i described in Eq. (8) in GRAVITY Collaboration (2021a).

Fig. 1 HD 98922 PIONIER (first two rows) and GRAVITY (last two rows) data, squared visibilities (left panel), closure phases (central panel), and (u,v) plan coverage (right panel) for two epochs of the complete data set shown in Figs. A.1 and A.2. Colors refer to the different spectral channels.

3.1.2 Global fit modeling

As the epochs are constituted by observations using different array configurations, different epochs probe different spatial frequencies of the source. This may lead to the detection of spurious variability in the model parameters, as shown in our test analysis in Appendix E. Therefore, instead of fitting each single epoch with our ten-parameter model and searching for individual variability trends, we implement a so-called global fit of the data. In this approach, we take advantage of the large temporal baseline of our data set to decipher which of the parameters tested could be the dominant parameter causing the variability visible in the data of Fig. 1. The global fit is obtained by forcing all the model parameters to be nonvariable over time, except for the parameter – or set of parameters – describing the feature to be tested for variability.

3.2 Temporal variability

Using the global fit approach on the continuum data, we explore which among the seven model parameters (star flux contribution F_s, ring inclination i, PA, characteristic size (a_r, a_k), and azimuthal modulation (c₁, s₁) is most prone to cause the temporal variability of the data visible in Fig. 1. The parameters describing the halo contribution F_h, the spectral index k_c, and the weight parameter fLor were not tested for temporal variability. We also report for comparison the analysis result for a fully nonvariable system model. The results are shown in Table 2.

For the GRAVITY data set, we observe from the ${\chi }_{\text{r}}^2$ analysis that the time-variable azimuthal modulation model #6 leads to the smallest ${\chi }_{\text{r}}^2$ value. The nonvariable system (#1) gives the poorest fit. This is less marked for the PIONIER data set for which the ${\chi }_{\text{r}}^2$ values are larger. However, the trend in terms of decreasing ${\chi }_{\text{r}}^2$ suggests an analogous behavior for both the PIO-NIER and GRAVITY data sets. This analysis indicates that, among the tested models, the one with a time-variable disk azimuthal modulation (model #6) best describes our data, and that such asymmetry is probably the dominant effect in the variability of the system, as opposed to other geometrical effects. The parameter values for the best-fit result of model #6 are shown in Table 3 and the corresponding marginal posterior distribution is presented in Appendices B.2 and B.3. The fit results for each individual epoch of our PIONIER and GRAVITY data sets are shown in Appendix A.1.

3.3 Disk azimuthal asymmetry

Our best-fit model exhibits a crescent-like asymmetric feature in the disk resulting from the azimuthal modulation that varies in PA through the epochs and revolves around the central star.

In Fig. 2, we present a subset of the continuum model images corresponding to our fitted models and from which the geometry, extent, and location of the asymmetry can be followed as a function of time. The complete time sequence is presented in Figs. F.1 and F.2. The fit of the variable azimuthal modulation appears quite robust in Fig. B.4, with one single global minimum identified in the c₁,s₁ diagrams. Figure 2 displays the azimuthal uncertainties in the form a white-line cone as derived from the 3σ error bars on c₁ and s₁. The azimuthal uncertainty is generally smaller (up to ~10–20°, depending on the configuration) in the K-band than in the H-band. In addition, the azimuthal position is most poorly constrained with the small configuration (e.g., P2–P5–P6 with PIONIER) because the closure phase signal is marginal for the shortest baselines For epochs only separated by a few days at most and taken with the same array configuration, the azimuthal locations of the asymmetric feature are consistent in most cases within the error bars (e.g., P7–P8–P9–P10–P11 and P19–P20 with PIONIER, or G7-G8 with GRAVITY). We also generally observe that the continuum emission appears azimuthally more compact in the H-band than in the K-band. The brightness contrast between the crescent-like feature and the corresponding centro-symmetric position in the disk is estimated from Figs. F.1 and F.2. The contrast is found to be ~4 in the K-band (ranging from ~ 1.6 to 10 across the epochs, with σ~2.8) and ~2.4 in the H-band (ranging from ~1.3 to 3.3, with σ ~0.7). No specific trend is found in the temporal evolution of the contrast. The dynamical properties of the inner disk feature revealed by our observations are further discussed in Sect. 3.4.

3.4 Orbital period of the crescent-like feature

The revealed azimuthal asymmetry in the inner disk of HD 98922 shows orbital motion around the central star. Assuming this is the same asymmetric feature that is monitored with the VLTI over the 11 yr period in both the H and K bands, we attempt to investigate the orbital properties of the emission feature. We show in Fig. 3 the distribution of the time-variable PA of the emission feature. The small-configuration epochs P2, P5, P6, G6, and G10 are not taken into account as the small baselines provide marginally weaker constraints on the azimuthal modulation because of the low level of closure phase signal detected.

The time origin corresponds to the first PIONIER observation of June 2011. The gray line is a sine function corresponding to a uniform circular motion fitted to the data to derive a period estimate. The best sine fit gives a period of 12.6±0.1 months. The ${\chi }_{\text{r}}^2$ value is large (111 ±51), which is due to the group of epochs at ~60 months that deviate from the best fit. If these epochs were found to correspond to outliers and were removed, the fit would give a ${\chi }_{\text{r}}^2$ value of 38±17, while the estimated period would remain the same (see Table 4). However, at this stage, there is no clear justification for the removal of these epochs.

Using Kepler’s law, we estimate the central mass of our target star for a separation of the azimuthal asymmetry ranging from 0.6 to 1.3 au based on the fitted ring annular radius a_r, as well as for a range from 1.2 to 1.7 au based on the half-light radius a, where $a={\left(a_r^2+a_k^2\right)}^{1/2}$ following Lazareff et al. (2017). From the derived orbital period of the azimuthal feature, we estimate a central mass of ~1–3 M_⊙ depending on the separation reported above. This is significantly lower than the literature value of ~6 M_⊙, which was robustly established using UVES high-spectral-resolution observations (Caratti o Garatti et al. 2015). For this higher mass, the estimated Keplerian orbit for a circular motion would be 9.2±0.4 months, which is shorter than that found with our measurement. It is noteworthy that our sine fit also shows a local minimum for a similar period of around ~8.5 months, although with a poorer ${\chi }_{\text{r}}^2$ value of ~350. We show this sine fit for comparison as a dashed line in Fig. 3.

Some caveats are worth mentioning here. The coverage of the orbital motion is relatively sparse and a finer temporal sampling would allow us to estimate the period with greater confidence. Furthermore, it should be noted that we have implemented here the simplest case of a circular orbit, because our model #6 of Table 2 foresees a nonvariable ring radius a_r. Accounting for possible eccentricity of the orbit of the crescent-like feature may impact the derived orbital period, which is not explored in this work. The orbital motion is further discussed in Sect. 6.1.3.

Fig. 2 Peak-normalized GRAVITY (top row) and PIONIER (bottom row) continuum model images. The dashed white lines represent the ±3σ uncertainty on the PA of the azimuthal modulation. The central object is not displayed but is marked with a star to enhance the circumstellar emission. North is up, east is to the left. See Appendix F for the full data set.

Fig. 3 Sine of the azimuthal modulation PAs as a function of time. The markers represent the sine of the variable dusty feature PAs depicted in Fig. 2, green for the PIONIER data and blue for the GRAVITY one. The full gray line represents the best fit (12.6 ± 0.1 months) of the uniform circular motion PA expressed as a sine function. The dashed gray line represents the fit solution corresponding to a period of 8.50 ± 0.25 months.

3.5 Fractional flux and disk characteristic size

In H-band, we find fractional flux ratios F_s, F_h, and F_c comparable to those found by Lazareff et al. (2017). In the K-band, we measure a stellar flux contribution of ~23% (see Table 3). Different values for the K-band fractional stellar contribution based on photometry and SED analysis are reported in the literature, ranging from ~15% (Kraus et al. 2008; Hales et al. 2014; GRAVITY Collaboration 2019) to ~23% (Caratti o Garatti et al. 2015). However, when comparing these values, we must keep in mind the different distance estimates used in these previous studies. We also find a stronger halo contribution than GRAVITY Collaboration (2019), which can be explained by the fact that these latter authors constrained the circumstellar emission of HD 98922 in the K-band only with the astrometric configuration, whereas the short baselines of the small configuration are usually needed to better constrain the halo.

The radius of the ring disk emission in K-band is found to be a_r=2.02±0.1 mas (or 1.31 ±0.07 au), in good agreement with the estimates of Kraus et al. (2008) and Caratti o Garatti et al. (2015), of namely 2.2 mas and 1.6 mas, respectively. Our result is slightly larger than the estimate of a_r=1.8 mas of GRAVITY Collaboration (2019), who only used two snapshots with the astrometric configuration. In the H-band, we find the circumstellar emission to be more compact than in K, with a characteristic radius of a_r=0.93±0.1mas (or 0.60±0.07 au), which is in line with the estimate of 0.87 mas by Lazareff et al. (2017). Kluska et al. (2020) report a larger half-flux radius of 2.1 mas based on image reconstruction, which evidences the impact of different modeling approaches. The K-to-H size ratio is further discussed in Sect. 6. The ratio a_k/a_r being close to or larger than unity suggests a wide, smooth ring emission as opposed to a sharp edge. The disk is known to have a low inclination, which is therefore more difficult to accurately constrain in the small angle range. This applies to the PA as well (Fig. B.1).

4 Results on the Brγ-line interferometric data

The GRAVITY line spectrum, the calibrated squared visibilities, and the differential phases are shown in Figs. 4 and A.3 for 12 of the 13 epochs after zooming into the spectral region between 2.164 and 2.168 µm. No high-spectral-resolution data could be acquired on June 15, 2018 (epoch G3). From the high-spectral-resolution K band data, HD 98922 shows a slightly blueshifted (≈−23.5 km s⁻¹) single-peaked Brγ emission line in all epochs. Considering the 3 Å spectral resolution, we can consider the peak’s position of the continuum-normalized line (${\mathrm{21\; 659.7}}_{-1.3}^{+0.8}$ Å) to be constant through the epochs at our spectral resolution. The normalized peak flux varies between 1.17 and 1.25 depending on the epoch, while the line width measured at the peak’s 10% flux level ranges from 14 to 15 Å. The total squared visibilities in the Brγ region vary between ≈0 and ≈0.4 depending on the epoch and baseline configurations (medium or large), while they reach ≈0.7 for the small configuration. Finally, the differential phase signals vary between −15° and 25° and have significantly different shapes (flat, single-peaked, double-peaked, or S-shape) for the different epochs and baselines.

Fig. 4 HD 98922 GRAVITY SC data for two different epochs (one per column) with similar (u,v) plan coverage. From top to bottom: Wavelength-calibrated and continuum-normalized spectrum, total squared visibilities (two panels), and total differential phases (two panels). Colors refer to the different baselines.

4.1 Modeling methodology

To estimate the Brγ gas region size, kinematics, and displacement with respect to the continuum emission from the GRAVITY SC visibilities and differential phases, we extrapolated the pure-line contribution (marked with subscript L) from the interferometric observables. Following Weigelt et al. (2011), the pure-line interferometric quantities characterizing the gas-emitting region, the visibility V_L and differential phase ϕ_L, are related as follows: $$F_{\text{L}}^2{V}_{\text{L}}^2=F_{\text{tot}}^2{V}_{\text{tot}}^2+F_{\text{cont}}^2{V}_{\text{cont}}^2-2F_{\text{tot}}{V}_{\text{tot}}{F}_{\text{cont}}{V}_{\text{cont}}\cdot \cos {\varphi }_{\text{tot}},$$(3) $$\sin {\varphi }_{\text{L}}=\sin {\varphi }_{\text{tot}}\frac{F_{\text{tot}}{V}_{\text{tot}}}{F_L{V}_L}.$$(4)

The quantities reported in Eqs. (3) and (4) refers to values inside the Brγ-line spectral region. As the continuum quantities F_cont and V_cont within the line are not directly measurable, they are estimated from the continuum near to the line region. However, hot Herbig stars exhibit a strong Brγ photospheric absorption feature, which needs to be accounted for in order to retrieve the correct pure-line quantities. In the case where photospheric absorption is present, knowledge from a model for the continuum emission is required. One can show that $$\begin{array}{l}{V}_{\text{L}}=(\frac{{\left(V_{\text{tot\hspace{0.17em}}}{F}_{\text{L}/\text{C}}^{\prime }\right)}^2+{\left(\frac{F_s^{\prime }}{F_{\text{cont\hspace{0.17em}}}^{\prime }}(\alpha -1)+V_{\text{cont\hspace{0.17em}}}^{\prime }\right)}^2}{{\left(F_{\text{L}/\text{C}}^{\prime }-\frac{\alpha +\beta +\gamma }{1+\beta +\gamma }\right)}^2}\\ {\frac{-2\left(V_{\text{tot\hspace{0.17em}}}{F}_{\text{L}/\text{C}}^{\prime }\right)\left(\frac{F_{\text{s}}^{\prime }}{F_{\text{cont\hspace{0.17em}}}^{\prime }}(\alpha -1)+V_{\text{cont\hspace{0.17em}}}^{\prime }\right)\cos \left({\varphi }_{\text{tot}}\right)}{{\left(F_{\text{L}/\text{C}}^{\prime }-\frac{\alpha +\beta +\gamma }{1+\beta +\gamma }\right)}^2})}^{1/2},\end{array}$$(5)

where the superscript (′) indicates that the quantities $V_{\text{cont\hspace{0.17em}}}^{\prime }$ and $F_{\text{cont\hspace{0.17em}}}^{\prime }$ are estimated outside the emission line spectral region, as opposed to Eq. (3). The line-to-continuum ratio $F_{\text{L}/\text{C}}^{\prime }$ is also normalized to the nearby continuum value. The parameters β and γ are the disk-to-star β=F_c/F_s and halo-to-star γ = F_h/F_s flux ratios.

The parameter α describes the star photospheric absorption, with 0 < α ⩽1, so that $$F_{\text{L}}=F_{\text{tot\hspace{0.17em}}}-F_{\text{cont\hspace{0.17em}}}^{\prime }+F_{\text{s}}(1-\alpha ),$$(6)

where α is equal to 1 when there is no absorption, and equal to zero when 100% of the star flux is absorbed (see Appendix C for a description of the photospheric absorption model).

Similarly, Eq. (4) is written in terms of known quantities as $$\sin {\varphi }_{\text{L}}=\sin {\varphi }_{\text{tot\hspace{0.17em}}}\frac{V_{\text{tot\hspace{0.17em}}}}{V_{\text{L}}}\frac{1}{\left(1-\frac{\alpha +\beta +\gamma }{1+\beta +\gamma }\frac{1}{F_{\text{L}/\text{C}}^{\prime }}\right)}.$$(7)

For a compact and marginally resolved gas emission component, we derive the photocenter displacement along each baseline from the pure-line differential phases as in Lachaume (2003): $$p=\frac{-{\varphi }_{\text{L}}}{2\pi }\cdot \frac{\lambda }{B},$$(8)

where p is the projection on the baseline B of the 2D photocen-ter vector with its origin on the continuum photocenter of the system.

4.2 Properties of the line-emitting region

We exploit the high-spectral resolution data of GRAVITY in the Brγ-line region to constrain the spatial scale and the kinematics of the hot-gas component. We used 10 out of 13 epochs for this analysis: in addition to epoch G3, for which no high-spectral resolution data could be acquired, epochs G6 and G10 are not considered, as no differential phase signal was detected using the small configuration.

From the pure-line visibilities, we estimated the characteristic size of the gas-emitting region at the peak wavelength of 2.1662 µm using a simple Gaussian disk model, as well as a ring model with 20% radial thickness for comparison. We obtain a Gaussian HWHM radius of 0.47 mas with a temporal standard deviation of 0.04 mas. Similarly, we obtain a radius of 0.50±0.05 mas for the thick ring. In both cases, the gas-emitting component model is found to be consistent with a face-on orientation, as we find an inclination of 5° ±5°. These values translate into a physical radius of ~0.3 au at a distance of 651 pc, in agreement with the results of Caratti o Garatti et al. (2015).

As the estimation of the characteristic size of the hot-gas region has been studied in previous works, we focus on the determination of the location of the compact hot-gas component in relation to the star, and further explore whether there is a correlation between the locations of the gas and the dust feature observed in the continuum. This can be investigated through an analysis of the interferometric differential phase, which provides precise information about the spatial location of the photocenter of the gas emission component on angular scales that surpass the nominal resolution of the interferometer.

The detection of a differential phase signal as reported in Fig. A.3 indicates that the photocenter positions of the continuum and gas emission components are not coincident. From the continuum-corrected (or pure-line) differential phases (Eq. (7)), we calculated the deprojected photocenter shifts of the hot-gas component with respect to the continuum photocenter for few spectral channels around the 0 km s^–1 velocity (Eq. (8)) following the formalism of GRAVITY Collaboration (2023). The location of the continuum photocenter may differ from the position of the central star if the dust emission is noncentrosymmetric. This is indeed the case for HD 98922, because our modeling of the continuum shows a strongly asymmetric time-variable brightness distribution of the inner circumstellar disk. Therefore, the location of the continuum photocenter is the parameter that most affects the position of the Brγ-line photocenter with respect to the central star. Figure G.1 shows the spatial location of the gas emission photocenter with respect to the continuum photocenter, and therefore with respect to the central star. For most of the epochs, the bulk of the Brγ emission appears offset – considering the retrieved characteristic size of the gas emission component – with respect to the stellar position by up to 0.5 mas. We also observe that the location of the compact gaseous component varies with time and that it is, on first order, consistently found in an area between the central star and the peak of the dusty feature, following its orbital motion. Looking closer at Fig. 5, we see that, for each epoch, the positions of the photocenter across the Brγ line are not spatially colocated but are distributed in a profile that qualitatively resembles a pattern of Keplerian motion.

5 Physical properties of the disk

In an earlier work, Hales et al. (2014) proposed a model for the disk of HD 98922 accounting for a distance of 507 pc and exploiting photometry data up to 160 µm. These authors proposed a disk model with an inner radius located at 1.5 au. We revisit this model in light of our new interferometric data.

5.1 Disk structure

We used the radiative transfer code MCMax (Min et al. 2009) to constrain the disk density and temperature structure based on archival broadband photometry data and on the spatial structure evidenced by our interferometric measurements. We started from a one-component model based on Hales et al. (2014). The inner disk radius is 1.5 au, the outer radius is 320 au, the flaring index γ=1, and the scale-height is 15 au at a radius of 100 au. We implement a grain population based on DIANA standard dust grains (Woitke et al. 2016) composed of 75% amorphous silicates (e.g., Mg₀₇Fe₀₃SiO₃), 25% porosity, by volume, and initially no amorphous carbon. The grain size a ranges from 1 µm to 2200 µm with a distribution of dn(a) ∝ a^–3.5, as large grains appear to dominate in this source (Bouwman et al. 2001; van Boekel et al. 2003). The surface density profile is based on a modified version from Li & Lunine (2003) with a power-law exponent of p_z = –1.5, a gas-to-dust ratio of 100, and a dust disk mass of 2 × 10^–5 M_⊙. The extinction is A_V = 0.5 mag (Guzmán-Díaz et al. 2021). The parameters of the central star are taken from Table 1.

With the revised stellar parameters, the stellar luminosity increases by a factor -3 in comparison to Hales et al. (2014) and the photospheric emission contribution results in overesti-mation of the K-band NIR flux in the SED by a factor ~1.2. On the other hand, our GRAVITY and PIONIER observations show that dust is present at separations of less than 1.5 au. As the SED fitting process is a degenerated problem, we attempt to better constrain the inner disk structure based on the consideration that a dust component at ≲1 au, producing ~75%–80% NIR excess and resolved with the VLTI has to be accounted for. Simply moving the inner radius to ~0.6 au in the one-component model described above is not possible because this leads to further over-estimation of the NIR flux in the SED. We therefore revise this disk model and propose a two-component disk model in which the inner component has a lower surface density in comparison to the outer component and an inner radius at 0.6 au.

The mass of the outer component is set to 2×10^–5 M_⊙ as in Hales et al. (2014). This value can be compared to the rough mass estimate obtained from the archival ALMA photometry point at 1.3 mm (ID: 2015.1.01600.S/PI Panic), which was not included in the work of these latter authors. We used the CARTA visualization tool and estimate a flux density at 1.3 mm of F_v=10.4 mJy for the unresolved HD 98922 source, with a typical conservative uncertainty of 10%. We compute the dust disk mass M_d following Beckwith et al. (1990). Assuming optically thin emission, we use the relation $$F_v\approx {\kappa }_v\frac{2k\overline{T}{v}^2}{c^2{d}^2}{M}_{\text{d}},$$(9)

where $\overline{T}$ is the dust mean temperature, k_v the dust absorption coefficient at 1.3 mm, d the distance, k the Boltzmann constant, and c is the speed of light. A large uncertainty subsists on the dust opacity value k_v and on the mean temperature, and therefore the mass value is only a crude approximation for first-order estimates. If we consider standard values of $\overline{T}=50\hspace{0.17em}\text{K}$, and k_v(1.3 mm) = 0.02 cm² g^–1 (Beckwith et al. 1990), we obtain a disk dust mass of M_d ~0.013 M_⊙. However, Woitke et al. (2016) underline the strong dependence of k_v on the fraction of amorphous carbon and grain size distribution, and report values as high as k_v(1.3 mm) ~5–10 cm² g^–1. The latter value results in a lower mass of M_d ~5 × 10^–5 M_⊙, which is comparable to that found by Hales et al. (2014). Because of the uncertain estimate of the disk dust mass, we choose to adopt the value reported by the latter authors for the outer disk component. For the disk composition, we considered two approaches: in the first case, we included a 25% fraction of carbon grains (model CS) to quench the strength of the silicate emission feature (van Boekel et al. 2003); in the second case, we consider an alternative case of an inner component only composed of quantum-heated particles (model Q) strongly coupled to the gas (e.g., Kluska et al. 2018; GRAVITY Collaboration 2021a).

We conducted a grid search for the inner disk dust mass M_d,in and for the transition radius between the low- and high-surface-density disk R_t (cf. Table 5) to match the SED as well as the NIR excess (F_c+ F_h) estimated by interferometry (cf. Table 3). The other disk parameters were kept identical for the inner and outer components. Table 6 shows our best result for the models CS and Q, respectively, for which our radiative transfer simulation produces the SED displayed in the top left panel of Fig. 6. For better readability, only the case of model CS is shown because the result for model Q is almost identical. Visually, these models provide the best overlap between the model and the data, except in the FIR region where the photometry data are more severely overestimated by the model. The FIR emission traces the disk at large radii not immediately relevant in our study. For the model CS, a transition radius Rt at -5.5 au and an inner disk dust mass of ~6×10^–9 M_⊙ are derived from the radiative transfer modeling (cf. Table 6). We note that our two best models CS and Q may still contain some level of degeneracy, because other parameters such as the scale height, the index of the surface density power law, or the flaring index may influence the NIR excess. However, we limited ourselves to the inner disk properties that can be directly constrained by the NIR flux ratios and spatial confinement of the emission derived from our interferometric measurements.

Fig. 5 Spatial and kinematic properties of the Bry-line-emitting gas region overlaid to the continuum model for two epochs. The complete sequence is shown in Appendix G. The star is centered on the origin. The black cross and the white circle around it show the position of the continuum photocenter and its uncertainty estimated from the continuum modeling. The red circle represents the extent of the gas-emitting region estimated at the peak of the line emission and centered at about ~0km s–1. The blue to red colored filled dots show the gas photocenter positions for the five spectral channels across the Bry line corresponding to velocities ranging from-100km s–1 to +100 km s–1, as color-coded in Fig. A.3.

5.2 Transition radius

For the adopted two-component disk model for HD 98922, the low- to high-density transition radiusR_t is found to be located at ~5.5 au, with a lower-limit at 2.5 au below which the NIR excess is systematically overestimated in the SED. Despite the degeneracy between R_t and M_d,in for the determination of the NIR excess flux in the SED fit, our interferometric measurements can set constrains on the inner radius, relative flux contribution, and compact spatial extent of the low-density component that dominates the NIR excess. Table 5 illustrates this point by comparing the relative flux contributions of the inner low- and outer high-density components inferred from radiative transfer to the estimates from the interferometric measurements, and this for two values of R_t, both giving a good fit of the SED. For smaller values of R_t, the NIR flux contribution of the outer high-density component increases, while that of the inner low-density component decreases. Qualitatively, this would translate for our geometrical models into wider rings, or into a higher contribution of the halo component, which we only find to be less than 10% in both bands. For instance, in the case R_t=3.5 au, the innermost low-density component probed by GRAVITY and PIONIER contributes only up to ~40%. However, it should be noted that the proposed argument suffers from some degree of uncertainty: as the modeling of our interferometric observables relies on single-ring geometrical models rather than on radiative transfer images, the accuracy on the determination of R_t remains limited. As the detailed interferometric modeling of radiative transfer images goes beyond the scope of this paper, we limit ourselves to propose that the structured disk of HD 98922 shows a dust density transition located no closer than ~4–6 au to the central star. Estimating an upper value for R_t is not feasible using only the presented data because the NIR excess (set on the first order by R_t, M_d,in and M_d,out) becomes a degenerated quantity for excessively large values of Rt.

5.3 Surface density profile

The dust surface density profile we implement has a sharp discontinuity at Rt, which does not correspond to a realistic case, given that a continuous transition profile would be more physically meaningful. Such sharp transitions are nonetheless proposed to differentiate between volatile-rich and volatile-free regions, as for instance in the case of the water snow line (Hayashi 1981; Lecar et al. 2006). More recent works have reported structured inner disks with regions of differentiated surface density (Tatulli et al. 2011; Matter et al. 2016). We derive a low- to high-density transition of two orders of magnitude at about 5 au (Fig. 6), with a dust surface density Σ_d ranging from ~2×10^–4 g cm^–2 to ~0.006 g cm^–2 in the low-density region. We emphasize that we cannot exclude that a gap is present in the disk of HD 98922 beyond ~5 au, but mid-infrared (MIR) interfero-metric data are necessary to explore this hypothesis. However, we also note that HD 98922 is classified as a group II source in the classification of Meeus et al. (2001) and based on the N-to-K size ratio (GRAVITY Collaboration 2019), suggesting a preferentially flat, gap-free disk.

5.4 Location and temperature of the dust

The dust emission in the H and K bands are expected to originate in a region close to the disk inner rim. In our case, the bulk of the emission is located at different radii for the two bands. We find the characteristic size of the K band is larger than that of the H band by a factor ~1.4 to 2 when considering either the ring annular radius a_r or the half-light radius a parameters. This trend is seen in previous works when comparing the results of Lazareff et al. (2017) and GRAVITY Collaboration (2019): out of 21 common sources, 17 show a larger K band by a factor ranging from 1.1 to ~3. At first glance, the difference in resolution between GRAVITY (1.75 mas) and PIONIER (1.22 mas) by a factor 1.4 seems capable of explaining the fact that PIONIER resolves a more compact emission, similarly to what observed for HD 163296 (GRAVITY Collaboration 2021b). A more physical, albeit simplistic, argument can be formulated based on the Wien’s temperature of the dust in a temperature-gradient disk (i.e., T(r) = T₀(r/r₀)^−q) emitting as a blackbody. For a typical value of q = 0.75 for a flat group II disk, we derive a size ratio of r_K/r_H ~(1350/1750)^1/0.75, or 1.41, in agreement with our findings.

One challenge posed by our results is the presence of H-band-emitting dust as close as ~0.6 au to the star. In the scenario of a passively irradiated disk with an optically thin inner cavity (Muzerolle et al. 2004; Monnier et al. 2005), the dust temperature at this location, assuming black-body grain emitters (ϵ = 1), is T_g ~ 2200 K, which is above the sublimation temperature of standard silicates. This is also seen in Fig. 6 for the temperature structure of our carbon-rich disk model. The problem of very hot dust grains inside the theoretical sublimation radius has been observed in other sources such as ZCMa, V1685Cyg, MWC297, and HD 190073 (Monnier et al. 2005; Hone et al. 2017; Setterholm et al. 2018), invoking the presence of refractory dust (e.g., corundum or iron), the sublimation temperature of which critically depends on the disk density (Kama et al. 2009). Alternatively, Benisty et al. (2010) proposed the possibility that a small fraction of highly refractory graphite grains with high sublimation temperature (T_s>2000 K) contribute to the NIR excess very close to HD 163296, and this might apply to HD 98922 as well. Clearly, further radiative transfer modeling is required to tackle this question by considering other values of a_min and p_Σ, as well as grain composition. We note that adding amorphous carbon to our initial pure-silicate model was essential to decrease the dust temperature at 0.6 au from T >3000 K to T ~ 2300 K.

Alternatively, for Herbig stars with masses similar to or larger than that of HD 98922, the undersized inner disk is explained by the presence of optically thick gas in the inner cavity, which is found to be favored in systems with accretion rates of Ṁ_acc ≳ 10⁻⁸ M_⊙ yr⁻¹ (Muzerolle et al. 2004), as for HD 98922.

Another possible explanation for the detection of excess emission very close to the star is the presence in the inner region of a fraction of quantum heated particles (QHPs), which can be stochastically heated by the strong UV radiation field from the central star, reaching temperatures higher than the equilibrium temperature and producing NIR continuum emission. This scenario was proposed for HD 100453, HD 179218, and HD 141569 (Klarmann et al. 2017; Kluska et al. 2018; GRAVITY Collaboration 2021a). Polycyclic aromatic hydrocarbons (PAHs), as an example of QHPs, are detected in HD 98922 (Geers et al. 2007; Acke et al. 2010). Interestingly, the 6.2-to-11.3 μm feature ratio is estimated from Seok & Li (2017)² to be I_6.2/I_11.3 ~3–4, with ratios larger than unity pointing at predominantly unshielded ionized PAH species in the disk inner regions that are directly exposed to the intense UV radiation field of the star (Maaskant et al. 2014). We tested a radiative transfer model with a PAH-based inner component (Model Q, Table 6), which allows a good fit of the SED with a very small mass of QHP grains (2×10⁻¹² M_⊙) and a transition radius of 3.5 au. This value of R_t is estimated as the balance between the increasing contribution of the outer high-density component (HDC) to the NIR SED (for R_t<3.5 au) and the increasing contribution of the inner low-density component (LDC) to the SED at λ ~1.5 μm (for R_t>3.5 au and due to the blue spectral index of the QHP emission). However, in the case of Model Q, the flux fractional contribution of the LDC is estimated to be only ~5% in the K-band and 10% in the H-band.

Ultimately, a mixture of refractory grains and QHPs could explain the presence of continuum emission as close as 0.6 au and detected in the H-band. The presence of predominantly ionized PAH species close to the star is therefore not to be discarded and requires further work.

Fig. 6 Radiative transfer modeling corresponding to Table 6. The top left panel shows the SED for Model CS, with the blue dashed line representing the stellar black-body function and the red line showing the modeled total emission. The black bars represent the photometric data. The top right panel shows the derived dust surface density profile as a function of the distance from the star. The bottom left plot shows the dust density structure, where the black dashed line represents the τ = 1 surface at 2.2 μm, and the red lines represent, from left to right and for both components, the density contours at 10−15, 10−16, 10−17, and 10−18 g cm−3, respectively. The bottom right plot shows the dust temperature structure, where the black lines represent, from left to right, the isothermal contours at 2300, 2000, 1700, 1500, and 1300 K, respectively.

6 Discussion

6.1 The origin of the time-variable inner disk asymmetry

Our analysis suggests a crescent-like asymmetric dust feature in the inner 1 au region of HD 98922 for which orbital motion is detected. Comparable azimuthal asymmetries are observed in other disks both at small and large scales: Varga et al. (2021) and GRAVITY Collaboration (2021b) characterized a time-variable crescent-like feature in the MIR and NIR in the innermost disk of HD 163296, whereas Ibrahim et al. (2023) directly imaged the complex asymmetric and time-variable inner rim of HD 190073 at similar wavelengths. Using ALMA, Casassus et al. (2015) and van der Marel et al. (2013) reveal a strong azimuthal asymmetry at separations of >50au in HD 142527 and Oph IRS 48, respectively. Infrared and submillimeter observations trace micronsized and millimeter-sized dust grains, respectively, and vortices or efficient dust traps are tentatively invoked for these sources.

Regarding HD 98922, it is remarkable that Caratti o Garatti et al. (2015) also detect large-scale disk asymmetry at ~40au from the central star. Imaged with SINFONI in the K-band, the emission is associated with scattered light, and the authors exclude the possibility that the arc-shaped structure is due to the presence of a possible stellar companion as close as ~20 au to the central source and with a mass of ≥ 0.5 M_⊙. Here, we explore a few hypotheses as to the nature and dynamical evolution of the asymmetric structure in the innermost disk of HD 98922.

6.1.1 Disk hydrodynamic instabilities

Inclination projection effects as in the case of FS CMa (Hofmann et al. 2022; Kluska et al. 2020) are very unlikely to be able to explain the azimuthal asymmetry in HD 98922. The strength of the closure phase signal – despite the low disk inclination (<30°) – and the detected orbital motion suggest a physical effect in the disk, possibly resulting from hydrodynamic instabilities.

Local perturbations in the inner disk may explain the variable azimuthal asymmetry under the assumption of strong dust–gas coupling. The Stokes parameter quantifies the degree of coupling and is given by Birnstiel et al. (2010): $$\text{St}=\frac{{\rho }_{\text{s}}\cdot a}{{\text{Σ}}_{\text{g}}}\frac{\pi }{2}.$$(10)

The condition for dust–gas coupling in the Epstein regime, namely St ≲ 1, is fulfilled for grain sizes of α ≲ 2∑_g/πρ_s. Considering a gas surface density at 1 au of 0.6 g cm⁻² (see Fig. 6) and a bulk density of silicate-rich grains of ρ_s ~3 g cm⁻³ (Pollack et al. 1994), we obtain a grain size of α ≲ 1300 μm, which suggests efficient dust–gas coupling for most grain sizes assumed in our model CS (see Table 6). Under this hypothesis, a vortex resulting from disk instabilities could result in the observed over-density. The vortex might be formed through Rossby wave instability (RWI, Meheut et al. 2010), which could be induced by the presence of a Jupiter-mass planetary companion. RWI-induced vortices occur preferentially in inviscid disks (typically α ≲ 10⁻⁴) associated to a low-turbulence environment, which allows the vortex to survive for thousands of orbits. It is therefore highly unlikely that we trace within the time baseline of our observations two unrelated over-densities that may have developed in different regions of the disk.

Regarding the morphology of the crescent-like feature, both Meheut et al. (2010) and Barge et al. (2017) show the formation of a vortex on one side of the disk, although the latter authors explore this aspect in the outer disk, at 60 au. The length/width aspect ratio of the simulated vortex in Meheut et al. (2010) and Varga et al. (2021) appears to be large, with the vortex covering at least one-third of the circumference. From our interferometric K-band observations, we model a crescent-like structure with an azimuthal span comparable to these authors (cf. Fig. 2). In Barge et al. (2017), the vortex appears azimuthally more compact. The latter authors include dust grains in their simulations and show that the contrast of the over-density is lowest for the smallest grains, with a value of ~ 10 for a size of 42 μm in the dust-density map. These latter authors report a fainter contrast of ~1.3 in the gas-density map for which no dust is included. Similarly, these levels of contrast compare well with the values derived from our geometrical modeling (see Sect. 3.3).

An azimuthally asymmetric feature can also arise from density waves due to magnetorotational turbulence (Flock et al. 2017). The observational signatures developed in this scenario also include a vortex visible in the form of an arc that develops at the edge of the dead zone, with a brightness contrast of comparable magnitude to that in the case of RWI.

It is noteworthy that, while the flux properties of the crescent-like feature are derived from a pure geometrical modeling of our interferometric data, the radiative transfer simulations predict an inner disk essentially optically thin at 2.2 μm due to the low surface density. It would therefore be interesting to include a nonaxisymmetric model in future radiative transfer simulations in order to better constrain the influence of the grain properties and size on the azimuthal NIR contrast of the crescent-like feature.

6.1.2 Disk warp induced by a close companion

An alternative scenario that can explain our observations is the presence of a close companion inside the cavity of what then becomes a circumbinary (CB) disk. Ragusa et al. (2017) show from 3D SPH gas and dust simulations that a crescent-like asymmetry develops at the CB disk inner edge only (and not further out in the disk) due to the dynamics of the central pair. The contrast of the nonaxisymmetric feature depends on the mass ratio of the components q, with q ⩾ 0.05 needed to produce a sufficiently contrasted over-density. This would translate in the case of HD 98922 into a companion mass of at least ~0.3 M_⊙ at a separation of ≲ 1 mas.

Baines et al. (2006) report a spectro-astrometric companion candidate with a separation of between 0.5″ and 3″ – that is, wider than 325 au –, which is not reported by Garufi et al. (2022) from the SPHERE total intensity images. Furthermore, the RUWE parameter reported from Gaia measurements is 0.94, which means that a single-star model provides, in principle, a good fit to the astrometric observations. Finally, a comparison of the photospheric lines obtained in 2005 and 2013 with the spectrograph/spectroplarimeter ESPaDOnS³ between 0.3 and 1 μm does not reveal any clear sign of variability in the profile or position of the line down to ~10 km s⁻¹ (Alecian et al., priv. comm.). With the realistic hypothesis that the period of such a companion does not coincide with the time lapse between the two ESPaDOnS epochs, we propose that no companion more massive than 0.3 M_⊙ at 0.1 au or 0.6 M_⊙ at 0.5 au is seen in the inner cavity. Of course, the low inclination of the system and the limited number of epochs mean that this preliminary conclusion requires further strengthening.

While the close binarity status of HD 98922 – in particular for companions lighter than 0.3 M_⊙ – cannot be totally excluded in the existing literature, we do not favor the idea that the presence of an inner stellar companion is inducing the crescent-like asymmetry.

6.1.3 Dynamics of the asymmetric feature

Thanks to our 11 yr temporal baseline, we suggest that we are tracing the same azimuthal asymmetric feature possibly interpreted as a vortex in the inner disk at ~1.0au. One hypothesis resulting from our work is that the bright over-density could show a sub-Keplerian motion with a period of ~ 1 yr. We emphasize however that this hypothesis requires further assessment in the future, since our disk data modeling only accounts for circular motion. Gas in protoplanetary disks is expected to orbit at sub-Keplerian rotation due to a radial pressure gradient with a deviation from Keplerian rotation of $\eta \propto c_{\text{s}}^2/v_{\text{K}}^2$. Even assuming gas temperatures of a few thousand Kelvin at ~1.0au, the deviation factor is η ~10⁻³, and would therefore not explain the difference we find between Keplerian and sub-Keplerian motion. Although the scenario of a close binary inside the disk cavity is not favored (Sect. 6.1.2), this is a case where more significant sub-Keplerian orbital velocity could arise. Ragusa et al. (2020) suggest that disk eccentricity due to the presence of a companion can result in the nesting of ellipses with a different eccentricity leading to a slowly precessing over-density through which the disk material passes. The caveat, in this case, is that the precession timescale of the over-density is of the order of ~100–1000 orbits, which should have resulted in a very steady feature considering our temporal baseline.

A more qualitative discussion concerns the structures formed in the disk by the presence of a planetary-mass companion. It is well known that the perturbing planet develops a pattern of wound spirals on either side of it (e.g., Kley 1999). The spiral arms are corotating with the planet, or, in other words, they are stationary in the reference frame of the planet. If the crescentlike feature in HD 98922 were associated with such a companion, this would lie at ~ 1.8 au to match the orbital period of 12 months, considering the literature value of 6 M_⊙ for the mass of the central star.

We underline that the question of the dynamics of asymmetric features in the innermost disk regions will become an important question for future interferometric observations. It is remarkable to note that in their recent imaging campaign of HD 190073, Ibrahim et al. (2023) present evidence that the detected sub-AU structure rotates two times slower than Keplerian, pointing at dynamics effects from the outer disk.

6.2 Kinematics of the hydrogen hot gas: Wind or companion?

The origin of the hydrogen Brγ-line in the system is debated: previous works suggest various scenarios such as a stellar wind, an X-wind, or magnetospheric accretion (Kraus et al. 2008), with the latter case being favored by these authors. On the other hand, strong variability over a timescale of a few days is instead detected in the blueshifted absorption lobe of the NaID and Balmer lines, which typically probes the inner dust-free cavity region (Aarnio et al. 2017). The blueshifted component can be modeled with a disk wind extending from 0.17 au to 0.42– 0.85 au and a wind mass-loss rate of 10⁻⁷ M_⊙ yr⁻¹. Similarly, a disk-wind model extending from ~0.1 au to ~1 au with a wind mass-loss rate of 2 × 10⁻⁷ M_⊙ along with an asymmetric continuum disk model can be successfully fitted to VLTI/AMBER Brγ-line interferometric data (Caratti o Garatti et al. 2015), which is in agreement with our continuum results.

Another interpretation involves the presence of an accreting very low-mass or planetary-mass companion. It is interesting to note that not only does the location of the compact Brγ-line emitting region change through the epochs, but so does the line luminosity. In the scenario of a young embedded planet, variability of the line luminosity resulting from changes in the accretion rates can be seen during the orbital motion of the planet (Szulágyi & Ercolano 2020).

By integrating the spectral line region through the observed wavelength-calibrated and continuum-normalized spectrum (Fig. A.3), we calculated the observed Brγ-line equivalent width $\left(W_{\mathrm{Br}\gamma }^{\text{obs}}\right)$ for each epoch, from which the line equivalent width corrected for photospheric absorption and veiling $\left(W_{\text{Br}\gamma }^{\text{corr}}\right)$ was derived (Table 7). Here, we discard the first two epochs G1 and G2 where the line is strongly affected by the continuum normalization and is therefore not highly reliable. We derive the flux density of the line F_Brγ in erg/sec/cm² by multiplying the spectral flux density of the nearby continuum, $$F_{\text{cont\hspace{0.17em}},\lambda }=F_{0,K}{10}^{-0.4\left(m_K-0.1A_V\right)}\text{,\hspace{0.17em}}$$(11)

with the corrected equivalent width $W_{\text{Br}\gamma }^{\text{corr}}$. The K band zero-point flux is F_0,K = 4.28×10⁻¹¹ erg/s/cm²Å⁻¹ (Rodrigo & Solano 2020). Finally, we compute the line luminosity as L_Brγ = 4πd² F_Brγ.

In Table 7, we see that in the four-month period from March 2019 (G4) to July 2019 (G8), the line luminosity L_Brγ is relatively constant at ~3.3 × 10⁻² L_⊙, while after 6 months (G9) it increases by ~10% and does not vary significantly for 11 months (G11), after which it increases again by ~5% (G12) and again by ~4% in the following week (G13), for a total increase of ~20%.

We note here that the values listed in Table 7 do not follow the known linear relationship between equivalent width and peak flux: W = A ⋅ F_peak + B. This could be due to the fact that our $W_{\text{Br}\gamma }^{\text{obs\hspace{0.17em}}}$ values are measured by integrating the spectral line region using the data themselves and not by integrating a Gaussian profile fitted to the data points. Additionally, the relative errors on W and F_peak may differ depending on the value of B. In any case, the key insight captured by the data is that the variation in both peak intensity and in $W_{\text{Br}\gamma }^{\text{obs\hspace{0.17em}}}$ is a positive one.

Under the hypothesis of an accreting planetary-mass companion, we use the empirical relation of Table 3 from Szulágyi & Ercolano (2020) that relates the Brγ-line luminosity to the mass of the accreting planet M_p: $$\log \left(\frac{L_{\text{Br}\gamma }}{L_{\odot }}\right)=a\cdot M_{\text{p}}+b,$$(12)

where the coefficients a and b depend on the opacity model. Considering the different opacity cases, we derive a mass M_p ranging between 10.3 M_J and 12.6 M_J.

A number of caveats must nevertheless be considered. First, the value of M_p should be treated as a loose estimate because the simulation from Szulágyi & Ercolano (2020) accounts for a model with a central stellar mass of 1 M_⊙, a circumstellar environment extending from 2 to 12.4 au, and a planet at 5.2 au. Furthermore, in this scenario, we make the implicit assumption that all the emission line flux is associated with the putative accreting companion, whereas a stellar component (e.g., in the form of a wind) may also contribute to the emission budget. In this scenario, the retrieved positions of the gas photocenter could result from the combined and weighted contributions from a stellar-driven component and a companion-driven component. Based on this argument, if it is a companion located further out at 1.8 au that causes the azimuthal asymmetry (see Sect. 6.1.3), a suitable weighting between the stellar-driven and companion-driven emission contributions could explain why the measured location of the combined photocenter is found between the star and the crescent-like feature. It is not the goal of this study to explore the case of a possible disk-embedded companion and its impact on the size, location, and dynamics of the Brγ-line-emitting region in detail; nonetheless, we speculate that a young accreting planetary-mass companion is a possible scenario to consider, with further modeling required.

7 Summary

We present new multi-epoch VLTI/GRAVITY observations of HD 98922 which, coupled with VLTI/PIONIER archival data, form an 11 yr observational-period data set for the system, accounting for a total of 33 different epochs between 2011 and 2022. This data set allows a unique interferometric study to test the potential time variability of the innermost circumstellar morphology. We can summarize the main conclusions of our work as follows:

• The system, which is spatially resolved by both instruments at all epochs, shows temporal variability dominated by changes in the asymmetric spatial distribution of brightness rather than in the flux ratio between the central star and its circumstellar environment. This is supported by the significant variability in the nonzero closure phase signals for observations obtained with comparable (u-v) coverage, whereas little time variability is observed in the continuum squared visibilities;

• Among the different modeling scenarios that we tested, the data are best explained by a crescent-like asymmetric dust feature radially extending from ~0.6 to ~2 au, which dominates the NIR excess. The revealed disk feature revolves around the central star with an estimated period of ~ 1 yr. This could point to sub-Keplerian orbital motion, which is also seen in the recently characterized system HD 190073; however, a shorter period cannot be excluded;

• The innermost location of the warm emitting dust traced by the interferometric data (~0.6 au) is coupled to revised radiative transfer models that suggest a radially structured disk formed by an inner low-density component (optically thin at 2.2 μm) and an outer high-density component with a transition between 3 to 5 au. The models favor either an inner region with a mass of 2× 10⁻⁹ M_⊙ made up of large, carbon-enriched silicate grains, or the presence of very small stochastically heated particles. The high dust temperature requires the presence of some level of refractory dust;

• The origin of the azimuthal asymmetry is preferentially connected to hydrodynamic effects, either induced by disk instabilities generating a vortex or by an undetected low-mass (substellar or planetary-mass) companion embedded in the low-density region of the disk and launching a wound spiral. An asymmetry resulting in disk fragmentation is unlikely considering the low surface density of <10 g cm⁻², and we do not find evidence for the presence of a close (0.1–0.5 au) stellar companion inside the inner dust rim that could drive eccentricity in the central cavity resulting in a crescentlike structure. However, we do believe that our data cannot completely rule out this hypothesis, in particular regarding a stellar companion less massive than ~0.3 M_⊙;

• The high-resolution GRAVITY data show a single-peaked Brγ emission line with a luminosity increasing by 20% in 3 yr. The location of the compact Brγ-line-emitting region is offset with respect to the central star and appears to be consistently located inside the dust cavity and between the dusty feature and the star;

• The interpretation of the interferometric observations of the hot, gaseous component in HD 98922 leaves some open questions: while a stellocentric magnetospheric accretion scenario is not favored, a wider asymmetric disk wind is a qualitatively plausible scenario, though more advanced numerical simulations would be needed to explain the time variability in its spatial distribution. We also speculate that a strongly accreting substellar or planetary companion with a mass of larger than ~10 M_J could explain our measurement of the Brγ emission line. However, this would trigger further questions regarding its dynamics and connection to the detected crescent-like structure in the inner disk.

Acknowledgements

The authors would like to thank Héloïse Meheut and Francois Ménard for the fruitful discussions on the inner disk of HD 98922. V.G. was supported for this research through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne, and from the Bonn-Cologne Graduate School of Physics and Astronomy (BCGS). A.C.G. has been supported by PRIN-INAF MAIN-STREAM 2017 and PRIN-INAF 2019 (STRADE), P. This work is based on observations made with ESO Telescopes at the La Silla Paranal Observatory under program IDs listed in Tables A.1 and A.2. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions in the Gaia Multilateral Agreement. We acknowledge the Gemini Observatory for the use of the IR spectrum model of the atmospheric transmission above Cerro Pachon. This research has made use of the model atmosphere grid NeMo, provided by the Department of Astronomy of the University of Vienna, Austria (Nendwich et al. 2004). NeMo was funded by the Austrian Science Fonds. This research has made use of the Jean-Marie Mariotti Center Aspro. This research has made use of the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIN/AEI/10.13039/501100011033/through grant PID2020-112949GB-I00. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.01600.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The authors thank the anonymous referee for assessing the quality of this work.

Appendix A Logs and observation data

Fig. A.1 HD98922 PIONIER data, squared visibilities, closure phases, and (u,v) plan coverage for each epoch. Colors refer to the different PIONIER spectral channels. Gray crosses represent the model described in Sect. 3.1 and Table 3. Gray circles in the bottom plots show the residuals.

Fig. A.2 HD 98922 GRAVITY FT data, squared visibilities, closure phases, and (u,v) plan coverage for each epoch. Colors refer to the different GRAVITY spectral channels. Gray crosses represent the model described in Sect. 3.1, and Table 3. Gray circles in the bottom plots show the residuals of the model fitting.

Fig. A.3 HD 98922 GRAVITY SC data for the different epochs. For each epoch, top plots show the wavelength-calibrated and continuum-normalized spectrum, left plots show the total squared visibilities, and right plots show the total differential phases. Circles represent the pure-line quantities. Colors refer to the different baselines.

Appendix B Global fit MCMC posterior distribution functions and azimuthal modulation parameters’ ${\chi }_{\text{r}}^2$ maps

Fig. B.1 Reduced chi-square ${\chi }_{\text{r}}^2$ curves of each nontime-variable parameter from the azimuthal modulation global fit model. The first two rows, where the curves are blue, refer to the PIONIER data model, while the last two rows, where the curves are black, refer to the GRAVITY data model.

Fig. B.2 Global fit MCMC posterior distribution functions of the fitted parameters for the PIONIER data set. In the one-dimensional histograms, the blue line identifies the median of the distribution, the dashed black lines identify the 16th and the 84th percentiles.

Fig. B.3 Global fit MCMC posterior distribution functions of the fitted parameters for the GRAVITY data set. In the one-dimensional histograms, the blue line identifies the median of the distribution, the dashed black lines identify the 16th and the 84th percentiles.

Fig. B.4 Global fit azimuthal modulation parameters ${\chi }_{\text{r}}^2$ maps. The red lines represent the 3 σ error bars of the parameters derived through the MCMC fitting procedure. Their intersections give the smallest ${\chi }_{\text{r}}^2$ for a given epoch.

Appendix C Spectrum wavelength calibration and star photospheric absorption model

The continuum normalization of the GRAVITY science object spectrum of each epoch (average of all the observations blocks files that come from the four different ATs for that epoch) was done by fitting the slope of the raw spectrum and dividing the latter by the resulting fit. The first step of the wavelength calibration was done by comparing the observed wavelength positions of the telluric lines in the HD 98922 spectrum for each epoch with respect to the positions of the telluric lines present in the IR spectrum of the atmospheric transmission above Cerro Pachon. We used pre-computed ATRAN models (Lord 1992) available from the Gemini Observatory website accounting for 4.3 mm water vapor column and an airmass of 1.5. The atmospheric transmission spectrum was convolved by a Gaussian with FHWM of 6 Å to have the same resolution as GRAVITY. The correction results in a blueshift between 0 and 5 Å depending on the epoch. The HD 98922 spectrum was then corrected for the star radial velocity (−4.9 km/s) and its proper motion with respect to the local standard of rest, with a value that depends on the epoch of observations (between ≈ − 35 and ≈2 km/s). The same corrections (telluric lines calibration, radial velocity correction, and local standard of rest correction) were also applied on the SC visibilities and differential phases for each epoch. The stellar atmospheric Brγ absorption was taken into account through a model selected from the Vienna New Model Grid of Stellar Atmospheres (Heiter et al. 2002) that best represents the star (see Table 1). The model accounts for a star effective temperature of 10500 K, a surface gravity logarithm of 3.4, a metallicity of [Fe/H] = −0.5, and a microturbulence of 2.0 km/s, a common value for a Herbig star. Finally, we included in the model the rotation broadening effect due to the rotation of the star (v sin i = 39.0 km/s) and we convolved the final model by a Gaussian with FHWM of 6 Å to have the same resolution as GRAVITY using SPECTRUM (Gray & Corbally 1994).

A systematic spectral shift of the line features was recently found in the GRAVITY SC data. The shift is seen when comparing data from an individual telescope with respect to the other telescopes, but also globally with regard to the known wavelength positions of the telluric lines in the K-band. It is thought that the shift is caused by the old grism of the science spectrometer which was upgraded in October 2019. GRAVITY data taken after this month should therefore no longer be affected. To check if our results are affected by the shift, we corrected the shift in the 19 March 2019 data following GRAVITY Collaboration 2023, and compared the results with the ones obtained without the correction. For a description of the correction we refer to Sect. 3.2 of GRAVITY Collaboration 2023. Figure C.1 shows the pure-line photocenter shift obtained with the original data (top panel) and that obtained with the corrected ones (bottom panel). We note a slight change in the shift when comparing photocenters of the same spectral channel. However, in this work we are not primarily interested in the relative position of the photocenters for each spectral channel, but rather in the overall location of the gas with respect to the star and the dusty feature. This choice is also based on the fact that our data were obtained with the ATs, meaning that we have only 4 or 5 spectral channels across the line. Data obtained with the UTs have instead around 12 spectral channels across the line, making these data more suitable for high-precision spectral analysis. From this point of view, the differences between the results obtained with the corrected and the uncorrected ATs data is negligible. Therefore, we decided to not apply this high-precision spectral calibration to our data.

Fig. C.1 Difference between pure-line photocenters shift obtained with data corrected with the GRAVITY data reduction software (Lapeyrere et al. 2014, top panel), and with the GRAVITY Collaboration (2023) high-precision spectral calibration (bottom panel). Same markers as in Fig. 5.

Appendix D Photometric data

Appendix E Interferometric variability and UV coverage

Appendix E.1 Variability in the visibilities and closure phases

Our data set on HD 98922 spans a 11 yr time period for the continuum emission, which gives us a unique opportunity to monitor intrinsic variability effects at a scale of 1au over an extended period of time. The time intervals between successive epochs range from 10 to 20 months for the PIONIER data set, and from 5 to 15 months for the GRAVITY one (see Table A.1 and Table A.2). Signatures of temporal variability in interfero-metric data are ideally investigated with configurations having comparable (u,v) plane coverage. For our data set, this is illustrated in Fig. E.1, which provides a visual estimate of comparable (u, v) coverage planes; although a more accurate assessment requires the comparison of the full configurations in Fig. A.1 and Fig. A.2.

In the top row of Fig. E.1, we report the length and PA of the longest baseline for a given configuration for PIONIER (P) and GRAVITY (G). The nomenclature can be followed in Table A.1 and A.2. Following this approach, we note qualitatively that a comparable (u, v) plane coverage is found for the following groups: (P1,15), (P2,6), (P3,4), (P7,8,9,10,17), (P11,12), (P13,14), and (P18, 19,20). Depending on the group, the difference in PA of the longest baseline varies between ~5° and 20°. For GRAVITY, the configurations with a similar (u, v) plane coverage are (G3,4,7,8,12), (G6,10), and (G9,13). The configurations of the remaining epochs differ more strongly. The central and bottom rows of Fig. E.1 are constructed as follows: first, we visually identify for a group of similar configurations the common spatial frequency for which a strong variation in the closure phase (CP) signal is observed and report the CP value (Fig. E.1, central panels). For this same spatial frequency, we report in the bottom panels the value of the squared visibility.

Fig. E.1 HD 98922 PIONIER and GRAVITY continuum data overview. Top row: Overview diagram summarizing Fig. A.1 and Fig. A.2 and showing array configurations with comparable (u,v) coverage. Each line shows the spatial frequency and PA probed by the longest baseline of a given four-telescope configuration. Two epochs have the same color and line style if they show a similar (u,v) coverage. The gray lines refer to epochs for which the (u,v) coverage cannot be reasonably compared to any other. Numbers identify epochs following Table A.1 and A.2; Middle row: Closure phase signal per epoch and instrument. The spatial frequency chosen is the one for which a strong variation in the CP signal between epochs with a comparable (u,v) plane coverage is observed. Bottom row: Squared visibilities signal for each epoch and instrument, for the same spatial frequency as in the middle row.

For those configurations with a similar (u,v) coverage, we are able to detect significant variations in the CP values across the epochs. For instance, a clear difference in the CP signal of ∆ϕ_max ≈ 25° is observed for G9 and G13, as well as in the medium configuration epochs G3, G4, G7, G8, G12 where the CP signals vary between ϕ_max ≈ −7° and ≈3°. Even though the error bars are larger, the CP variability in the PIONIER data set is also clearly observed. The epochs P1 and P15 show a significant closure phase variation (∆ϕ_max ≈ 10°), as does P17 with respect to P7, P8, P9, and P10 (∆ϕ_max ≈ 13°). On the contrary, when comparing observations taken with similar (u,v) plane coverage, we note that the variability of the corresponding squared visibilities is typically smaller than V² ~0.05.

Appendix E.2 Fit dependence on the (u,v) plane coverage

In the context of the study of spatial variability, we first verified using synthetic data whether the fit of different data sets related to the same nonvariable system but obtained with different baseline configurations would lead to the same or to a different solution. For this purpose, we generated with Aspro⁴ three synthetic data sets from an input image formed by a star and an azimuthally modulated ring observed with the GRAVITY small, medium, and large configurations, respectively, on the 4 February 2020, 11 July 2019, and 14 February 2022. We then fitted the different synthetic data sets with the same parametric model and then compared the results obtained for each parameter.

The input model image, shown in Fig. E.2, corresponds to an azimuthally modulated torus around a central star that contributes 15% of the total flux of the system, and with no spatially resolved halo contribution. The synthetic interferomet-ric visibilities and closure phases generated with Aspro are shown in Fig. E.3. The same error bars as for the real data are adopted. The fit of the synthetic data implements the nine free parameters given in Table E.1. The spectral indices of the star and the disk were not implemented and were set to zero. Four cases are explored, corresponding to the individual small, medium, and large configurations plus a case where all the three configurations are combined.

From the results in Table E.1, we find that the flux parameters (F_s, F_h, F_c) are the most affected by some spurious variability due to the varying (u,v) coverage, with differences of up to ~30 %. In addition, the stellar contribution is slightly (large configuration) to significantly (small configuration) overestimated with respect to the input model. The geometrical parameters (a_r, a_k) describing the characteristic size are found to be more in agreement (though not necessarily within the error bars), except for the small configuration that only loosely constrains the solution. The inclination value is consistent within the error bars between the four cases, although modestly constrained by the small and medium configurations. Also, the results are found to be generally consistent with a low-inclination disk. The same conclusion applies to the PA, with the small configuration expectedly failing to constrain this parameter. Finally, for the parameters c₁ and s₁ describing the azimuthal modulation, the fit converges towards very similar values for all four cases.

Simply accounting for the error bars on the fitted parameters, the modeling of the small-, medium-, and large-configuration data sets results in a consistent value for the ring inclination, the PA, c₁, and s₁, but clearly these parameters are more stringently constrained when going from the small to the large configuration. In Fig. E.2, we compare the input model image to the parametric models resulting from the fit of the different configurations. Besides the differences found in the relative flux contributions (see Table E.1), we can visually observe that the size properties are not in good agreement – in particular for the small configuration – when fitting the data of the different configurations separately, although the general structure is retrieved. In contrast, the azimuthal position of the disk asymmetry appears to be well constrained. This simple test reminds us that the intrinsic sparsity of infrared interferometric data implies that different spatial frequencies are probed with different configurations, which may result in a spurious variability of the fitted parameters that may not reflect a physical temporal variability of the system.

Based on this observation, we adopted here a different strategy to more robustly test the potential time-variability of the different parameters describing our model: we fitted the full PIO-NIER (respectively GRAVITY) continuum data set by forcing all the free parameters to be constant across the different epochs, except the parameter for which we wish to test the variability hypothesis. Throughout the paper we refer to this approach as the x global fit, where x is the selected time-variable parameter.

Fig. E.2 Aspro synthetic data test images for the small, medium, and large configurations of GRAVITY. The top-left model was used to generate the synthetic visibilities and closure phases with Aspro. The resulting fitted models are shown for the respective configurations. The intensity scale is linear and ranges from zero to the peak pixel value of the asymmetry. The central star is not displayed and its position is marked with the white cross. The images have a size of 20×20 mas. It is visible from this figure and Table E.1 that the small configuration properly constrains the absence of halo contribution in the system, but not the disk morphology, as opposed to the large configuration.

This is motivated by the fact that a rich (u,v) coverage remains desirable to constrain our model.

Fig. E.3 Aspro synthetic data (white markers), Vis2, closure phases, and (u,v) plane (left, center, and right panel, respectively) for the different baseline configurations. Blue markers represent the model described in Sect. E.2. Bottom panels show the residuals of the fitting process.

Appendix F Visualisation of the continuum geometrical models

Fig. F.1 Peak-normalized PIONIER continuum model images. The dashed white lines represent the ±3σ MCMC uncertainties on the PA of the azimuthal modulation. A description of the model is given in Sect. 3. The central object is not displayed in order to enhance the circumstellar emission. The contrast reported in Sect. 3.3 is calculated from the peak-normalized brightness distribution map as the ratio between the brightest pixel in the map and the corresponding centro-symmetric position in the disk. We report the standard deviation for the mean contrast value across the different PIONIER epochs.

Fig. F.2 Same as Fig. F.1 but for the GRAVITY data set.

Appendix G Pure-line photocenter displacements

Fig. G.1 GRAVITY continuum model images and Brγ-line-emitting gas region. The black filled dot shows the continuum (star + halo + disk) photocenter. The white circle around that point gives the error on the continuum photocenter position estimated from the continuum modeling. The red circle represents the gas region size estimated at the peak emission and is centered on the ~0 km s−1 point. The blue to red colored filled dots show the gas photocenters for different spectral channels, as in Fig. A.3.

References

All Tables

All Figures

	Fig. 1 HD 98922 PIONIER (first two rows) and GRAVITY (last two rows) data, squared visibilities (left panel), closure phases (central panel), and (u,v) plan coverage (right panel) for two epochs of the complete data set shown in Figs. A.1 and A.2. Colors refer to the different spectral channels.
In the text

	Fig. 2 Peak-normalized GRAVITY (top row) and PIONIER (bottom row) continuum model images. The dashed white lines represent the ±3σ uncertainty on the PA of the azimuthal modulation. The central object is not displayed but is marked with a star to enhance the circumstellar emission. North is up, east is to the left. See Appendix F for the full data set.
In the text

Fig. 3 Sine of the azimuthal modulation PAs as a function of time. The markers represent the sine of the variable dusty feature PAs depicted in Fig. 2, green for the PIONIER data and blue for the GRAVITY one. The full gray line represents the best fit (12.6 ± 0.1 months) of the uniform circular motion PA expressed as a sine function. The dashed gray line represents the fit solution corresponding to a period of 8.50 ± 0.25 months.

In the text

	Fig. 4 HD 98922 GRAVITY SC data for two different epochs (one per column) with similar (u,v) plan coverage. From top to bottom: Wavelength-calibrated and continuum-normalized spectrum, total squared visibilities (two panels), and total differential phases (two panels). Colors refer to the different baselines.
In the text

Fig. 5 Spatial and kinematic properties of the Bry-line-emitting gas region overlaid to the continuum model for two epochs. The complete sequence is shown in Appendix G. The star is centered on the origin. The black cross and the white circle around it show the position of the continuum photocenter and its uncertainty estimated from the continuum modeling. The red circle represents the extent of the gas-emitting region estimated at the peak of the line emission and centered at about ~0km s–1. The blue to red colored filled dots show the gas photocenter positions for the five spectral channels across the Bry line corresponding to velocities ranging from-100km s–1 to +100 km s–1, as color-coded in Fig. A.3.

In the text

Fig. 6 Radiative transfer modeling corresponding to Table 6. The top left panel shows the SED for Model CS, with the blue dashed line representing the stellar black-body function and the red line showing the modeled total emission. The black bars represent the photometric data. The top right panel shows the derived dust surface density profile as a function of the distance from the star. The bottom left plot shows the dust density structure, where the black dashed line represents the τ = 1 surface at 2.2 μm, and the red lines represent, from left to right and for both components, the density contours at 10−15, 10−16, 10−17, and 10−18 g cm−3, respectively. The bottom right plot shows the dust temperature structure, where the black lines represent, from left to right, the isothermal contours at 2300, 2000, 1700, 1500, and 1300 K, respectively.

In the text

	Fig. A.1 HD98922 PIONIER data, squared visibilities, closure phases, and (u,v) plan coverage for each epoch. Colors refer to the different PIONIER spectral channels. Gray crosses represent the model described in Sect. 3.1 and Table 3. Gray circles in the bottom plots show the residuals.
In the text

	Fig. A.2 HD 98922 GRAVITY FT data, squared visibilities, closure phases, and (u,v) plan coverage for each epoch. Colors refer to the different GRAVITY spectral channels. Gray crosses represent the model described in Sect. 3.1, and Table 3. Gray circles in the bottom plots show the residuals of the model fitting.
In the text

	Fig. A.3 HD 98922 GRAVITY SC data for the different epochs. For each epoch, top plots show the wavelength-calibrated and continuum-normalized spectrum, left plots show the total squared visibilities, and right plots show the total differential phases. Circles represent the pure-line quantities. Colors refer to the different baselines.
In the text

	Fig. B.1 Reduced chi-square ${\chi }_{\text{r}}^2$ curves of each nontime-variable parameter from the azimuthal modulation global fit model. The first two rows, where the curves are blue, refer to the PIONIER data model, while the last two rows, where the curves are black, refer to the GRAVITY data model.
In the text

	Fig. B.2 Global fit MCMC posterior distribution functions of the fitted parameters for the PIONIER data set. In the one-dimensional histograms, the blue line identifies the median of the distribution, the dashed black lines identify the 16th and the 84th percentiles.
In the text

	Fig. B.3 Global fit MCMC posterior distribution functions of the fitted parameters for the GRAVITY data set. In the one-dimensional histograms, the blue line identifies the median of the distribution, the dashed black lines identify the 16th and the 84th percentiles.
In the text

	Fig. B.4 Global fit azimuthal modulation parameters ${\chi }_{\text{r}}^2$ maps. The red lines represent the 3 σ error bars of the parameters derived through the MCMC fitting procedure. Their intersections give the smallest ${\chi }_{\text{r}}^2$ for a given epoch.
In the text

	Fig. C.1 Difference between pure-line photocenters shift obtained with data corrected with the GRAVITY data reduction software (Lapeyrere et al. 2014, top panel), and with the GRAVITY Collaboration (2023) high-precision spectral calibration (bottom panel). Same markers as in Fig. 5.
In the text

Fig. E.1 HD 98922 PIONIER and GRAVITY continuum data overview. Top row: Overview diagram summarizing Fig. A.1 and Fig. A.2 and showing array configurations with comparable (u,v) coverage. Each line shows the spatial frequency and PA probed by the longest baseline of a given four-telescope configuration. Two epochs have the same color and line style if they show a similar (u,v) coverage. The gray lines refer to epochs for which the (u,v) coverage cannot be reasonably compared to any other. Numbers identify epochs following Table A.1 and A.2; Middle row: Closure phase signal per epoch and instrument. The spatial frequency chosen is the one for which a strong variation in the CP signal between epochs with a comparable (u,v) plane coverage is observed. Bottom row: Squared visibilities signal for each epoch and instrument, for the same spatial frequency as in the middle row.

In the text

Fig. E.2 Aspro synthetic data test images for the small, medium, and large configurations of GRAVITY. The top-left model was used to generate the synthetic visibilities and closure phases with Aspro. The resulting fitted models are shown for the respective configurations. The intensity scale is linear and ranges from zero to the peak pixel value of the asymmetry. The central star is not displayed and its position is marked with the white cross. The images have a size of 20×20 mas. It is visible from this figure and Table E.1 that the small configuration properly constrains the absence of halo contribution in the system, but not the disk morphology, as opposed to the large configuration.

In the text

	Fig. E.3 Aspro synthetic data (white markers), Vis2, closure phases, and (u,v) plane (left, center, and right panel, respectively) for the different baseline configurations. Blue markers represent the model described in Sect. E.2. Bottom panels show the residuals of the fitting process.
In the text

Fig. F.1 Peak-normalized PIONIER continuum model images. The dashed white lines represent the ±3σ MCMC uncertainties on the PA of the azimuthal modulation. A description of the model is given in Sect. 3. The central object is not displayed in order to enhance the circumstellar emission. The contrast reported in Sect. 3.3 is calculated from the peak-normalized brightness distribution map as the ratio between the brightest pixel in the map and the corresponding centro-symmetric position in the disk. We report the standard deviation for the mean contrast value across the different PIONIER epochs.

In the text

	Fig. F.2 Same as Fig. F.1 but for the GRAVITY data set.
In the text

Fig. G.1 GRAVITY continuum model images and Brγ-line-emitting gas region. The black filled dot shows the continuum (star + halo + disk) photocenter. The white circle around that point gives the error on the continuum photocenter position estimated from the continuum modeling. The red circle represents the gas region size estimated at the peak emission and is centered on the ~0 km s−1 point. The blue to red colored filled dots show the gas photocenters for different spectral channels, as in Fig. A.3.

In the text