The GRAVITY young stellar object survey XIII. Tracing the time-variable asymmetric disk structure in the inner AU of the Herbig star HD98922

Temporal variability in the photometric and spectroscopic properties of protoplanetary disks is common in YSO. However, evidence pointing toward changes in their morphology over short timescales has only been found for a few sources, mainly due to a lack of high cadence observations at mas resolution. We combine GRAVITY multi-epoch observations of HD98922 at mas resolution with PIONIER archival data covering a total time span of 11 years. We interpret the interferometric visibilities and spectral energy distribution with geometrical models and through radiative transfer techniques. We investigated high-spectral-resolution quantities to obtain information on the properties of the HI BrG-line-emitting region. The observations are best fitted by a model of a crescent-like asymmetric dust feature located at 1 au and accounting for 70% of the NIR emission. The feature has an almost constant magnitude and orbits the central star with a possible sub-Keplerian period of 12 months, although a 9 month period is another, albeit less probable, solution. The radiative transfer models show that the emission originates from a small amount of carbon-rich (25%) silicates, or quantum-heated particles located in a low-density region. Among different possible scenarios, we favor hydrodynamical instabilities in the inner disk that can create a large vortex. The high spectral resolution differential phases in the BrG-line show that the hot-gas component is offset from the star and in some cases is located between the star and the crescent feature. The scale of the emission does not favor magnetospheric accretion as a driving mechanism. The scenario of an asymmetric disk wind or a massive accreting substellar or planetary companion is discussed. With this unique observational data set for HD98922, we reveal morphological variability in the innermost 2 au of its disk region.


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 ⋆ GRAVITY is developed in a collaboration by the Max Planck Institute for Extraterrestrial Physics, LESIA of the Paris Observatory, and IPAG of the Université Grenoble Alpes/CNRS, the Max Planck Institute for Astronomy, the University of Cologne, the Centro de Astrofísica e Gravitação, and the European Southern Observatory.
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 et al. 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) Article number, page 1 of 45  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 et al. 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 et al. 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 ≲ 500 au (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;GRAV-ITY Collaboration et al. 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.

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 et al. 2019) and the H band emission inside ∼1.2 au (Lazareff et al. 2017).

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 Article number, page 2 of 45 (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 minutes 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 GRAV-ITY dataset spans an 11 year period from 2011 to 2022.

Data reduction
The PIONIER data are archival reduced and calibrated data retrieved from the JMMC Optical interferometry DataBase 1 .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 et al. ( 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 wavelengthcalibrated 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.

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, ) 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, ) 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 2 ∼ 0.05.

Modeling methodology
The observational results show that, for the epochs using VLTI configurations with a similar (u, ) 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.

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 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 infinitesimally 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 f Lor; and the azimuthal modulation parameters c 1 and s 1 .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, ) and at wavelength λ is therefore described by a linear combination of the three components as 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 circumstellar environment, respectively, at λ 0 = λ K,0 = 2.15 µm (the wavelength of the central spectral channel of the GRAVITY FT), or at λ 0 = λ 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.645at λ 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 (cf.Table 3) 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.minimizeusing 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 χ 2 r by computing the error on the mean of the stochastic variable T i described in Eq. 8 in GRAVITY Collaboration et al. (2021a).

Global fit modeling
As the epochs are constituted by observations using different array configurations, different epochs probe different spatial fre- quencies 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.

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 1 ,s 1 )) 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 f Lor 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 χ 2 r analysis that the time-variable azimuthal modulation model #6 leads to the smallest χ 2 r value.The nonvariable system (#1) gives the poorest fit.This is less marked for the PIONIER data set for which the χ 2 r values are larger.However, the trend in terms of decreasing χ 2 r suggests an analogous behavior for both the PI-ONIER 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

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 1 ,s 1 diagrams.Figure 2 displays the azimuthal uncertainties in the form a white-line cone as derived from the 3σ error bars on c 1 and s 1 .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 GRAV-ITY).We also generally observe that the continuum emission appears azimuthally more compact in the H-band than in the Kband.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.

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 year 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 es-   timate.The best sine fit gives a period of 12.6±0.1 months.The χ 2 r 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 χ 2 r 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 = (a 2 r + a 2 k ) 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 highspectral-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.4months, 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 χ 2 r 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.

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 (2) 12.6 ± 0.1 38 ± 17 3.0±0.8 (3) Notes: (1) period estimated from the fit of the time-variable PA of the continuum azimuthal feature; (2) same as ( 1) but with the group of epochs at 60 months removed.The mass of the central star is derived by assuming that the feature has a separation ranging from 0.6 to 1.3 au as given by the ring radius a r ; (3) same as ( 1) but with the mass of the central star derived by assuming that the feature has a separation ranging from 1.2 to 1.7 au as given by the half-light radius a.

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 Å) 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, singlepeaked, double-peaked, or S-shape) for the different epochs and baselines.

Modeling methodology
To estimate the Brγ gas region size, kinematics, and displacement with respect to the continuum emission from the GRAV-ITY 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 gasemitting region, the visibility V L and differential phase ϕ L , are related as follows: The quantities reported in Eq. 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 Article number, page 7 of 45 A&A proofs: manuscript no.aa46926-23_corr emission is required.One can show that where the superscript ( ′ ) indicates that the quantities V ′ cont and F ′ cont are estimated outside the emission line spectral region, as opposed to Eq. 3. The line-to-continuum ratio F ′ L/C 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 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 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): where p is the projection on the baseline B of the 2D photocenter vector with its origin on the continuum photocenter of the system.

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.05mas 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 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  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.

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.

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 0.7 Fe 0.3 SiO 3 ), 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 Σ = −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 overestimation 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 overestimation 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 ν =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 where T is the dust mean temperature, κ ν 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 κ ν 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 T =50 K, and κ ν (1.3 mm) = 0.02 cm 2 /g (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 κ ν on the fraction of amorphous carbon and grain size distribution, and report values as high as κ ν (1.3 mm) ∼ 5-10 cm 2 /g.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 et al. 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-surfacedensity disk R t (cf.Table 6) 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 5 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 R t at ∼ 5.5 au and an inner disk dust mass of ∼ 6×10 −9 M ⊙ are derived from the radiative transfer modeling (cf.Table 5).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 Article number, page 9 of 45 disk properties that can be directly constrained by the NIR flux ratios and spatial confinement of the emission derived from our interferometric measurements.

Transition radius
For the adopted two-component disk model for HD 98922, the low-to high-density transition radius R 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 6 illustrates this point by comparing the relative flux contributions of the inner lowand outer high-density components inferred from radiative trans- fer 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 highdensity 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 PIO-NIER 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 R t .

Surface density profile
The dust surface density profile we implement has a sharp discontinuity at R t , 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) interferometric 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 et al. 2019), suggesting a preferentially flat, gap-free disk.

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)   temperature-gradient disk (i.e., T (r) = T 0 (r/r 0 ) −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 Hband-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 Z CMa, V1685 Cyg, MWC 297, 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 −8 M ⊙ yr −1 (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 et al. 2021a).Polycyclic aromatic hydrocarbons (PAHs), as an example of QHPs, are detected in Article number, page 11 of 45 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) 2 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 5), which allows a good fit of the SED with a very small mass of QHP grains (2×10 −12 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.

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: Regarding HD 98922, it is remarkable that Caratti o Garatti et al. ( 2015) also detect large-scale disk asymmetry at ∼40 au 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.

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 2 The ratio is uncertain due to incomplete wavelength coverage 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): The condition for dust-gas coupling in the Epstein regime, namely St≲1, is fulfilled for grain sizes of a ≲ 2Σ g /πρ s .Considering a gas surface density at 1 au of 0.6 g cm −2 (see Fig. 6) and a bulk density of silicate-rich grains of ρ s ∼ 3 g.cm−3 (Pollack et al. 1994), we obtain a grain size of a ≲ 1300 µm, which suggests efficient dust-gas coupling for most grain sizes assumed in our model CS (see Table 5).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 −4 ) 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 crescentlike 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.

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 Article number, page 12 of 45 GRAVITY Collaboration: V. Ganci et al.: The GRAVITY young stellar object survey 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 ESPaDOnS3 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 −1 (Alecian et al., private communication).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.

Dynamics of the asymmetric feature
Thanks to our 11 year temporal baseline, we suggest that we are tracing the same azimuthal asymmetric feature possibly interpreted as a vortex in the inner disk at ∼1.0 au.One hypothesis resulting from our work is that the bright over-density could show a sub-Keplerian motion with a period of ∼1 year.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 η ∝ c 2 s /v 2 K .Even assuming gas temperatures of a few thousand Kelvin at ∼1.0 au, the deviation factor is η ∼ 10 −3 , 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 Na I D 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.85au and a wind mass-loss rate of 10 −7 M ⊙ yr −1 .
Similarly, a disk-wind model extending from ∼ 0.1 au to ∼ 1 au with a wind mass-loss rate of 2× 10 −7 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 (W obs Brγ ) for each epoch, from which the line equivalent width corrected for photospheric absorption and veiling (W corr Brγ ) 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 2 by multiplying the spectral flux density of the nearby continuum, with the corrected equivalent width W corr Brγ .The K band zero-point flux is F 0,K = 4.28×10 −11 erg/s/cm 2 /Å (Rodrigo & Solano 2020).Finally, we compute the line luminosity as L Brγ = 4πd 2 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 −2 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 This could be due to the fact that our W obs Brγ 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 obs Brγ 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 : 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 stellardriven 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.

Summary
We present new multi-epoch VLTI/GRAVITY observations of HD 98922 which, coupled with VLTI/PIONIER archival data, form an 11 year 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-) 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 year.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 −9 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 −2 , 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 crescent-like 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 years.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.
Article number, page 14 of 45 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. 1 .
Fig. 1.HD 98922 PIONIER (first two rows) and GRAVITY (last two rows) data, squared visibilities (left panel), closure phases (central panel), and (u, ) 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.
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.

Fig. 2 .
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 .
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 GRAV-ITY 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.

Fig. 4 .
Fig.4.HD 98922 GRAVITY SC data for two different epochs (one per column) with similar (u, ) 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.

Fig. 5 .
Fig. 5. Spatial and kinematic properties of the Brγ-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 ∼ 0 km s −1 .The blue to red colored filled dots show the gas photocenter positions for the five spectral channels across the Brγ line corresponding to velocities ranging from -100 km s −1 to +100 km s −1 , as color-coded in Fig. A.3.
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 Article number, page 8 of 45 GRAVITY Collaboration: V. Ganci et al.: The GRAVITY young stellar object survey component with respect to the continuum photocenter for few spectral channels around the 0 km/s velocity (Eq.8) following the formalism of Gravity Collaboration et al. (

Fig. 6 .
Fig. 6.Radiative transfer modeling corresponding to Table 5.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.
Varga et al. (2021) and GRAVITY Collaboration et al. (2021b) characterized a timevariable 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 > 50 au in HD 142527 and Oph IRS 48, respectively.Infrared and submillimeter observations trace micron-sized and millimeter-sized dust grains, respectively, and vortices or efficient dust traps are tentatively invoked for these sources.

Fig. B. 4 .
Fig. B.4. Global fit azimuthal modulation parameters χ 2 r maps.The red lines represent the 3 σ error bars of the parameters derived through the MCMC fitting procedure.Their intersections give the smallest χ 2 r for a given epoch.

Fig. C. 1 .
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 et al. (2023) high-precision spectral calibration (bottom panel).Same markers as in Fig. 5.

Table 2 .
Results of the global fit for different parameters tested for variability (Models #2 to #6).Model #1 corresponds to a nonvariable system.The uncertainty on the χ 2 r is reported.

Table 3 .
Nonvariable parameters of the azimuthal modulation global fit continuum model.F c is not a free parameter, but is obtained followingF c =1-F h -F s .The σ MCMC error estimates derived through the MCMC fitting procedure are given by the 16 th and 84 th percentiles of the samples in the MCMC marginalized distributions.
Article number, page 5 of 45

Table 4 .
; GRAV-ITY Collaboration et al. 2019) to ∼23 % (Caratti o Garatti et al. 2015).However, when comparing these values, we must keep in Article number, page 6 of 45 GRAVITY Collaboration: V. Ganci et al.: The GRAVITY young stellar object survey Orbital period estimates of the dust azimuthal asymmetry ± 0.1 38 ± 17 0.9±0.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).

Table 6 .
RT results for Model CS.The sign( * )indicates that the relative flux contributions of the disk and halo components are added.LDC and HDC stand for low-density and high-density component, respectively.

Table 7 .
Table 7 do not follow the known linear relationship between equivalent width and peak Properties of the Brγ-line emission.The relation between W obs

Table D .
1. HD 98922 photometry data Article number, page 35 of 45