The GRAVITY Young Stellar Object Survey. VI. Mapping the variable inner disk of HD 163296 at sub-au scales

Protoplanetary disks drive some of the formation process (e.g., accretion, gas dissipation, formation of structures, etc.) of stars and planets. Understanding such physical processes is one of the main astrophysical questions. HD 163296 is an interesting young stellar object for which infrared and sub-millimeter observations have shown a prominent circumstellar disk with gaps plausibly created by forming planets. This study aims at characterizing the morphology of the inner disk in HD 163296 with multi-epoch near-infrared interferometric observations performed with GRAVITY at the Very Large Telescope Interferometer (VLTI). Our goal is to depict the K-band (lambda_0 ~ 2.2 um) structure of the inner rim with milliarcsecond (sub-au) angular resolution. Our data is complemented with archival PIONIER (H-band; lambda_0 ~ 1.65 um) data of the source. We performed a Gradient Descent parametric model fitting to recover the sub-au morphology of our source. Our analysis shows the existence of an asymmetry in the disk surrounding the central star of HD 163296. We confirm variability of the disk structure in the inner ~2 mas (0.2 au). While variability of the inner disk structure in this source has been suggested by previous interferometric studies, this is the first time that it is confirmed in the H- and K-bands by using a complete analysis of the closure phases and squared visibilities over several epochs. Because of the separation from the star, position changes, and persistence of this asymmetric structure on timescales of several years, we argue that it is a dusty feature (e.g., a vortex or dust clouds), probably, made by a mixing of sillicate and carbon dust and/or refractory grains, inhomogeneously distributed above the mid-plane of the disk.


Introduction
The formation and evolution of protoplanetary disks is crucial in the formation process of stars and planets.They are key laboratories for magneto-hydrodynamic, radiative and astro-chemical processes.Understanding the processes at play in these dust-and gas-rich circumstellar disks is one of the main science cases of several current observing facilities.ALMA (ALMA Partnership et al. 2015) and adaptive-optics assisted imagers, like SPHERE (Beuzit et al. 2019), have revealed the complexity and the diversity of the outer parts of these disks (e. g., ∼ 20-500 au).Most of them exhibit diverse features as gaps, rings, spirals, vortex, and shadows (Andrews et al. 2018;Long et al. 2018;Benisty et al. 2015;de Boer et al. 2016;Pohl et al. 2017;Benisty et al. 2017;Avenhaus et al. 2018) whose origin is still a matter of debate.Knowing the disk properties at different spatial scales and, in particular, in the innermost regions close to, or within the dust sublimation front, is critical to understand the conditions for planet formation and migration in proto-stellar disks around young stars (Flock et al. 2019).
Zooming into these innermost regions requires angular resolution down to a few milliarcseconds (mas) or less, which is made only possible with optical long-baseline interferometry.
Several statistical studies of the dust inner rims of young stellar objects (YSOs) have been conducted using instruments at the Very Large Telescope Interferometer (VLTI) in H (λ 0 ∼ 1.65 µm) and K (λ 0 ∼ 2.2 µm) bands (Lazareff et al. 2017;Anthonioz et al. 2015;GRAVITY Collaboration et al. 2019).Recent reconstructed images of 15 Herbig stars in the H−band (Kluska et al. 2020) show that 60% exhibits disk asymmetries at a few au scale, and 27% has a non-centrosymmetric morphology that can be due to a non-axisymmetric and/or variable environment (Chen et al. 2019).
The Herbig Ae star HD 163296 has a mass of 2.23 ± 0.22 M (Alecian et al. 2013) and is located at a distance of 101 ± 1.2 pc (Vioque et al. 2018).This makes it one of the closest Herbig stars to us.Its protoplanetary disk has been well studied in a wide spectral range, from the optical to the millimetric part.The ALMA DSHARP project (Andrews et al. 2018) reveals several rings around this star and the authors report detailed azimuthal asymmetries in the outer disk (Huang et al. 2018;Isella et al. 2018).Such features may be produced by planets, whose presence is also suggested by deviations from Keplerian rotational motions (Teague et al. 2018(Teague et al. , 2019;;Pinte et al. 2018).Multiepoch study of HD 163296 in the optical with coronagraphic imagers revealed a temporally variable, non-azimuthally symmetric illumination of the outer disk (Rich et al. 2019(Rich et al. , 2020)), while J-band polarimetric observations detected an off-center ring around HD 163296, which appears to be caused by scattering of the upper layers (Monnier et al. 2017).Also, it is one of the few Herbig stars to drive a prominent jet (Ellerbroek et al. 2014).
Optical interferometric studies of HD 163296 show the presence of an elongated dusty disk with a diameter of ∼5 mas with an inclination of ∼45 • and a position angle of ∼130 • (Lazareff et al. 2017;Setterholm et al. 2018;GRAVITY Collaboration et al. 2019).In these works, the closure phases show a deviation from zero on the long baselines, which hints at an asymmetric structure or clumpy emission in the innermost disk region.In particular, Setterholm et al. (2018) show strong closure phase variations for the longest baselines sampled with CLIMB-CHARA.These authors suggest the presence of an asymmetric emission at scales smaller than ∼2 mas (∼0.2 au).
Recently, Kluska et al. (2020) reconstructed an image of HD 163296 using PIONIER-VLTI (H−band).These authors discussed how a possible asymmetry located at the inner edge of the disk would move during the observing time spanning the interferometric observations, and how this could be affecting the interpretation of the observed morphology.For example, a body in Keplerian rotation at a radius of 3 mas (or 0.3 au) would complete a full orbit in ∼40 days, assuming a central mass of 2.23 M .This puts strong constraints in the way interferometric data could be combined for imaging.In contrast, it also motivates monitoring campaigns to depict the nature of this asymmetry.Additionally, Kobus et al. (2020) suggested symmetric variability at sub-au scales based on the relative changes of the squared visibilities observed in PIONIER and AMBER (λ 0 = 2.2 µm) data.However, these analyses do not include the characterization of the closure phases.Therefore, the asymmetry of the source could not be fully depicted.More recently, Varga et al. (2020) detected an asymmetric structure using MATISSE data in the L-band (λ 0 = 3.8 µm), which produces a brighter side on the extended emission in the inner 4 mas (0.4 au) of the source.This asymmetric brightness distribution has been identified before by Lazareff et al. (2017) and Kluska et al. (2020).However, the interesting aspect is that the position angle of the asymmetry detected with MATISSE appears to be significantly different from the one reported in the PIONIER observations.Thus, this suggests a varying morphology of the inner disk.However, the limited MATISSE data on large baselines did not allow a proper characterization of the asymmetry and its supposedly variable nature.
This work is a further step in the analysis of the inner ≤ 1 au (≤ 10 mas) of HD 163296.We aim at depicting the morphology of the inner disk by using multi-epoch interferometric observations obtained with GRAVITY-VLTI in K-band and archival PI-ONIER H−band data.Our specific goal is to characterize the predicted asymmetric structure in a systematic way across the H and K bands.For this purpose, we present a geometrical model to the squared visibilities (V 2 ) and closure phases (CPs).The manuscript is divided as follows: Sec.2 presents the interferometric data used for this work; Sec. 3 describes our model and Sec.4 our results; in Sec. 5, we present our discussion and; in Sec.6 our conclusions are presented.

Observations and Data Reduction
2.1.GRAVITY data HD 163296 was observed with GRAVITY (Eisenhauer et al. 2008(Eisenhauer et al. , 2011;;Gravity Collaboration et al. 2017) during 2018 and 2019 as part of the Young Stellar Object (YSO) Guaranteed Time Observations (GTO).The log of the observations is reported in Table A.1.The observations were done using the highest spectral resolution (R ∼ 4000) of the instrument, deploying the science source as Fringe Tracker (R ∼ 22, sampling ∼1 kHz; Gillessen et al. 2010).This allows us to stabilize the fringes of the science beam combiner for up to several tens of seconds.All the data were recorded using the 1.8 m Auxiliary Telescopes (ATs) with the intermediate array (D0-G2-J3-K0) 1 .The interferometric observables (squared visibilities and closure phases) were obtained using the instrument's data reduction software provided by ESO (Lapeyrere et al. 2014).
Absolute calibration of the science data is implemented through the observations of point-like sources interleaved with the science sequences.The data reduction software estimates the instrumental transfer function by correcting the observed calibrator visibilities with the theoretical ones according to the estimated angular size of the calibrator.Finally, the algorithm corrects the science raw observables by the estimated transfer func-HD 163296 -Best-fit azimuthally modulated ring model tion.Before analyzing the data, all squared visibilities (V2 ) with a signal-to-noise ratio (SNR) lower than five and closure phases (CPs) with uncertainties above 40 • were discarded.
Post-processing of the data was done with custom-made python routines.In this step, science data were wavelength rebinned from the native ∼1700 channels of the high-resolution data sets down to 8 spectral channels across the K−band, which increased the signal-to-noise for our analysis of the continuum.We account for calibration errors, while keeping the SNR statistics of the different data sets.We obtained average error bars of σ V 2 = 0.01 and σ CPs = 0.3 • for the 2019 data and σ V 2 = 0.01 and σ CPs = 0.8 • for the 2018 data, respectively.Our analysis is based on the exploitation of data from the science channel of GRAV-ITY and not from the fringe-tracker channel of the instrument.Detailed analysis of emission lines like Brγ were discarded and are the subject of a forth-coming paper.GRAVITY u-v planes are included in Figure A.1 in the Appendix.

Ancillary Data
To complement the analysis and results obtained from our GRAVITY data, we used the available PIONIER-VLTI data sets of HD 163296 taken during 2013, and 2017, which are included into the Optical Interferometry Data Base (OIDB 2 ) hosted by the Jean-Marie Mariotti Center (JMMC).These data are reduced and calibrated in a consistent way, which allows us to compare and merge them directly.PIONIER data are of low-spectral resolution, sampling three and six spectral channels across the H−band (λ 0 = 1.65 µm) in 2013, and 2017, respectively.The log of these observations are also listed in Table A.1 in the Appendix.2.12 2.9 2.17 3.0 2.9

Azimuthally Modulated Ring Model
To analyze the flux distribution at the inner disk of HD 163296, the PIONIER and GRAVITY observables were fitted with an azimuthally modulated ring (hereafter called Ring model) based on the prescription proposed by Lazareff et al. (2017).It is composed of an infinitesimal wire-frame ring with the following Fourier transform: here, to account for the possible elongation of the ring, the original u-v coordinates were rotated and inclined in the following form: here, PA is the position angle (measured from East to North) of the semi-major axis of the ring, and cos(i) is the elongation factor.q exp (iψ) is the polar form of the spatial frequencies u r -v r .a r is the ring angular radius.ρ j exp(iθ j ) is the polar representation of the modulation amplitude, c j + is j , applied to the ring profile.The index j corresponds to the order of the modulation, in this case, we only tested models with j = 1.
To provide a width to the ring wire-frame, its profile is convolved with a Lorentzian kernel of the following form: where a k is the Kernel angular radius.To avoid degeneracy in the fitting process, we fixed the quotient a k / a r = 0.3.The complete visibility expression is the following one: here, F c , F s and F h are the flux contributions of the ring (denoted with the subscript "c"), point-like object (denoted with the subscript "s") and over-resolved component (denoted with the subscript "h").We explicitly constrain the model to have For the minimization, we used the Gradient Descent Least-Squares algorithm implemented in lmfit.The starting points of the parameters use values close to the position angle, inclination, size and flux ratios reported in Lazareff et al. (2017).The c 1 and s 1 variables, which define the degree of asymmetry in the ring, were initially set to zero.Initial values were the same for all epochs and instruments.Each spectral channel was fitted independently.The 2013 datasets were combined into a single epoch.Since this corresponds to the same dataset used by Lazareff et al. (2017) and Kluska et al. (2020), this serves us to compare our best-fit models directly with those previous estimates.For the PIONIER data sets taken in 2017, we only combined data taken with a maximum separation of three days and with a minimum of four snapshots.The 2018 and 2019 GRAVITY data are used as separate epochs.The best-fit parameters obtained for our model are reported in Table 1 and Table 2. Figure 1

Constraints on the geometry of the target
Our parametric model reproduces the observables for all the epochs.The morphological characteristics of the source obtained with our best-fit model are the following ones: 1. Geometry: The radii of the ring vary from 1.75 -2.7 mas (0.175 -0.27 au) for the 2013, 2018 and 2019 epochs, and of 0.65 -1.2 mas (0.065 -0.12 au) for both 2017 epochs.The 2013 PIONIER data helped us to compare our results with previous estimates in the literature.In this regard, the estimated radius of the ring, a r , is in agreement with the 1.81 mas reported by Lazareff et al. (2017) using the same dataset.
We expect to observe a change in the radius of the ring between the H and K bands due to a radial gradient in the temperature profile of the dust.However, we observed that there is an important difference in the obtained values, even within the PIONIER epochs (a 2013,June/July r ∼ 1.78 mas; a 2017,April r ∼ 1.19 mas and; a 2017,Aug.r ∼ 0.8 mas).Our model predicts a projected disk semi-major axis position angle (North to East) of 130 • -140 • .This range of parameters is consistent with previous infrared interferometric findings, as well as with ALMA observations of the disk structure at large scale (see e.g., Andrews et al. 2018 The inclination angle of our model varies from 35 • to 55 • depending on the epoch.It is interesting to notice that, despite covering similar u-v frequencies, the GRAVITY data show a strong difference in the derived inclination angle between the 2018 (i = 53.1 ± 1.5) and the 2019 (i = 34.8± 0.7) epochs.We consider that this apparent difference in the inclination angle derived by our model is caused by the asymmetric and variable structure of the ring obtained from the closure phase information.Muro-Arena et al. (2018) suggest a small misalignment from +1 • to +3 • between the spatially unresolved inner disk (with a radius between 4 and 10 mas) and the outer structure of the disk (which has an inclination of i = 45 • ).The models presented by those authors support this claim in order to explain the outer shadow casts observed in the disk.It is interesting to notice that the inner inclination angle derived by Muro-Arena et al. (2018), based on SPHERE/IRDIS 2016 data, is consistent with the value derived with our 2017 April data (i = 48.3± 0.4).2. Flux Distribution: F * contributes between 33% and 45% of the total flux, depending on the epoch, while F c contributes between 40% to 60%, and F h contributes between 0% to 20%.The observed changes in the flux contributions of the PIONIER data could be, partly, explained by the sampling of different u-v frequencies between the different epochs.We notice that these large variations (more than 3σ) are also present in the GRAVITY data.Since our GRAVITY epochs have coincident baselines, the observed variations in the GRAVITY data cannot be related to flux filtering due to the different u-v spacing.In contrast, our findings suggest that there is an intrinsic flux variability associated with changes in the morphology of the source.This is in direct line with the results from Kobus et al. (2020) who interpret the variability of the squared visibilities for similar baselines at different epochs as changes in the structure of the target.3. Asymmetry: We can appreciate that the loci of the brightest side of the Ring model is changing, which strongly supports the variability of the inner ring structure.Notice that the values of the coefficients c 1 and s 1 , which traces the modulation of the ring, are changing between epochs.Observing these changes is of particular interest for the GRAVITY data, since the 2018 and 2019 epochs sample quasi-coincident spatial frequencies (in length and position angle).Figure 2 displays the GRAVITY observables for the 2018 epoch over-plotted with the observables extracted from the best-fit models obtained from the 2019 epoch, together with the opposite case.
It can be noticed that observables extracted from the model in one epoch are not able to reproduce the observables of the other one, despite having quasi-coincident u-v frequencies.This demonstrates changes in the morphology of the object in, at least, temporal scales of 1 year for the GRAV-ITY data.These time-scales are restricted by the temporal sampling of our data.However, we cannot exclude changes in the structure at shorter time-scales and/or the possibility of being observing different bright structures instead of a persistent single one.
The presence of an asymmetry and hints in the variability were suggested before by several authors in the literature.
The two most recent ones are Kobus et al. (2020) and Varga et al. (2020).The first work uses archival PIONIER data and parametric models applied to the squared visibilities in order to derive morphological changes in the structure of the source.These authors infer that the source must have symmetrical variations.Nevertheless, they do not include models to fit simultaneously closure phases and squared visibilities.Therefore, they cannot prove the asymmetric nature of the source.The second work uses MATISSE data and parametric models to constrain the source in L-band.However, this study includes a limited number of large baselines, where the asymmetry is detected.Therefore, the authors cannot confirm the variability of the target solely with the MATISSE data.Moreover their work just traces the asymmetric structure of the extended emission at one position angle.In contrast, our findings go one step forward to confirm the asymmetry of the source, its variable nature, and persistence (on a temporal baseline of 7 years).

Limitations of the azimuthally modulated ring model
While the model presented in Sec. 3 reproduces the different trends observed in the data, we have identified the limitations discussed below: Comparison with other parametric models: Our azimuthally modulated ring model uses a valid prescription to describe the morphology of the target at the spatial scales traced with PIO-NIER and GRAVITY.Similar ring models have been presented in the literature (Lazareff et al. 2017;Kluska et al. 2020;GRAV-ITY Collaboration et al. 2019;Varga et al. 2020) with good results to reproduce the morphology of the target.Nevertheless, we consider that our model presents several limitations that need to be taken into account when interpreting the morphology of the target.These limitations are related to: (i) a limited u-v coverage, particularly for the GRAVITY data; (ii) the lack of baselines larger than 100 m to clearly sample the visibility trend after the first rebound in the Fourier transform of the ring and; (iii) the lack of similar baselines and wavelength sampling on all epochs.
We complemented our analysis by exploring another parametric model based on an off-centered Gaussian plus a central point-like source and an over-resolved component (hereafter called Gaussian model).A detailed description of this Gaussian model is included in Appendix C. By comparing the two models, we confirmed the following points: -Both models are able to reproduce the trend observed in the squared visibilities and closure phases with similar degree of accuracy (see the reduced χ 2 values reported in Tables 1 -2 and Tables D.1 -D.5 in Appendix C).Therefore, both are equally valid solutions to describe the morphology of HD 163296.Similar inclination and position angles are found between the Ring model and the Gaussian one.The Gaussian model also shows the difference in the inclination angle of 50 • and 30 • between the 2018 and 2019 GRAVITY data.However, we find a larger discrepancy in the contribution of F s to the total flux compared with the Ring model.While the Ring model predicts values of F s between 33% and 45% of the total flux, the Gaussian model predicts values between 20% and 30%.This discrepancy can be explained by the fact that the Fourier transform of a Gaussian is a continuous and monotonically asymptotes to the F s value as spatial frequencies increase.Therefore, an additional contribution from the Gaussian is always present on top of the flux contribution of the central source to reproduce the visibility trends.This fact compensates the lower percentage of F s obtained with the best-fit Gaussian model compared with the Ring one.
The off-centered Gaussian model also reproduces the asymmetries of the source structure.In this model, the direction of the peak's displacement (asymmetry vs central source) is in agreement with the brightest side of the Ring model for each different epoch.Nevertheless, the amplitude of the displacements of the Gaussian component is one order of magnitude smaller than the position of the brightest part of the Ring.This is because the displacements of the Gaussian component are tracing the flux-centroid position of the extended morphology and not the position of the asymmetry in the structure.Therefore, caution must be taken when comparing those values with the position of the brightest side of the Ring model.-Given the limitations of our Ring model, it is interesting to compare it with previous formulations in the literature, which use a Gaussian to describe the extended morphology of the target.For example, Setterholm et al. (2018) use radially symmetric models: a Gaussian, a uniform disk and an infinitesimal ring.All of them include a centered point-like source The data used by those authors consisted in a combination of different instruments from the VLTI and CHARA.Their u-v sampling include baselines up to 350 meters and their models were applied to the H and K bands.Their results suggest that an on-axis elongated Gaussian and a point-like source reproduce better the visibility function of the target.
In particular, those authors found problems to reproduce the short angular scales (the ones traced with baselines < 50 m.) with the uniform disk and infinitesimal ring models.However, we do not find a similar problem to reproduce the visibilities at short spatial scales using our Ring model.One important difference between the ring model presented by Setterholm et al. (2018) and ours is that those authors use an infinitesimal ring, which produces more pronounced rebounds in the visibilities after the first minimum.In our case, the wire-frame of our azimuthally modulated ring is convolved with a Lorentzian kernel which smooths the profile of the ring and produces flatter rebounds after the first minimum in the visibility trend.Additionally, the uncertainties in the K−band data presented by those authors have values Comparison with reconstructed images: Our geometrical models reproduce the observed profiles in our data.However, those geometrical models can only explore a limited degree of asymmetries in the data.In order to better explore the asymmetries traced by the closure phases, we performed image reconstruction on our datasets.The complete imaging process is described in Appendix D. The regularized minimization used for image reconstruction is able to trace more complex asymmetric structures in the data.It is true that imaging is more suitable for rich u-v coverages, therefore, our limited data constrain the quality of the reconstructions.Still, by comparing independent images from the data and reconstructed ones from our parametric models is possible to improve our knowledge about the asymmetric morphology of HD 163296.
From the best images obtained, we could not favor the Ring nor the Gaussian model as the best one to describe the target.This is because reconstructed images obtained from their simulated observables are quite similar to the recovered images obtained using our data.However, we observed that, for some spectral channels, there are residuals as large as 20% of the peak in the images.This is similar for both the Gaussian and Ring models.Therefore, this supports that the degree of asymmetry of the source traced by our parametric-models is underestimated and the morphology of the target is more complex than what we map with the geometrical models.

Temperature of the asymmetry:
The main result of our analysis is the confirmation of the asymmetric and variable structure of the inner disk morphology of HD 163296.We expect to have an optically thick inner disk.Therefore, estimating the temperature of the dust (T d ) at the positions of the asymmetry cannot be obtained directly with its surface brightness.However, we can compute a rough estimate of T d by assuming dust grains of given sizes directly heated by the UV radiation of the central source, which are located at a distance r from the star at the position of the brightest point in our ring model.For this purpose, we use the following expression (van Buren & McCray 1988): here, L * ,38 is the UV luminosity of the star in units of 10 38 erg s −1 ; r pc is the distance to the dust from the star in units of parsecs; and ,a, is the size of the dust particle in microns.T d was computed considering a power-law dust size distribution Mathis et al. (1977), for dust sizes ranging from 10 −2 to 1 µm.To calculate the distribution of temperatures per epoch, we extracted 10 4 different samples of the peak's positions, dust grain sizes, and we set L * ,38 = 6.38×10 34 erg s −1 (Acke & van den Ancker 2004).The deprojected peak's positions were obtained from the Ring model, assuming an inclination i = 40 • .This produces the temperatures reported in Table 3.These values ( 1240 K < T d < 1600 K) are in agreement with temperatures between the sublimation point of the silicate (T s = 1500 K) and carbon (T c = 1800 K) dust grains.However, we should notice that the reported range of temperatures only traces the material observed with the different interferometric datasets.Therefore, those temperature values do not necessary correspond to the upper temperature limit of the most heated dust in the disk.
Previous near-IR interferometric studies (Benisty et al. 2010;Setterholm et al. 2018) suggested the possibility of having refractory dust grains (which survive temperatures above T > 2000 K).At mid-infrared wavelengths, measurements obtained with MIDI-VLTI and reported by van Boekel et al. (2004) found a considerable larger fraction of cristallinity within the central 20 mas (2 au) in HD 163296, compared with the outer 20-200 mas (2-20 au) of the disk (40% ± 20% vs 15% ± 10%).Similarly, the ring models reported by Varga et al. (2020) support that around 20% of the surface brightness near the star comes from a region where small micron-sized standard dust grains cannot survive.Hence, those authors, also, suggest the presence of refractory grains with small cooling efficiencies ( ∼ 0.1 = 0.18) which survive temperatures above 2300 K.These constraints and our derived T d values support the presence of a mixing dust species or refractory dust grains as responsible for the thermal emission of the variable disk.

On the nature of the asymmetry
The derived changes in the size of the fitted ring cannot be explained solely by the temperature gradient since the different epochs are sampled at the same wavelength.Therefore, they can only be explained by the changes in the effective resolution between the different epochs (see Fig.  ture.In contrast, they support the existence of a more smooth inner morphology of the disk. From our best-fit models, the peak of the emission in the ring changes for each one of the different epochs.Figure 3 displays the projected positions of the peak in the ring emission for PI-ONIER and GRAVITY, complemented with the peak's position from the ring model applied to the MATISSE observations described in Varga et al. (2020).Unfortunately, we have too limited amount of data to clearly trace the orbital motion of the material and to test the presence of a single persistent structure, instead of several different ones, over the 7 years spanning our observations.Furthermore, due to the change in angular resolution, the apparent size of the ring changes by a factor of two between the PIONIER and the GRAVITY models (see Fig. 1).Our limited resolution, does not allow us to clearly resolve the asymmetry, therefore, we cannot determine whether its forming material is distributed on a well-localized structure or, if it is more extended over several angular scales.In this section, we discuss several physical scenarios to explain the origin of the observed asymmetry.
Due to the change in distance and position angle of the asymmetry across the different epochs and instruments, we discard the possibility of being observing an illumination effect due to a fixed inclination of the disk.A possible cause of such an asymmetry is that the inner disk is warped and, therefore, casts different shadows on the outer disk.Very recently, Kraus et al. (2020) discovered a highly misaligned and warped disk around GW Orionis.These authors propose disk tearing (Facchini et al. 2013) as the hydrodynamic effect that causes the inner disk to change its orientation and precession.However, this mechanism is only possible if there is at least a stellar binary as central engine for the system.In the case of HD 163296, there is no evidence of a secondary stellar companion.Furthermore, the inner disk shows small or no precession, as indirectly seen from the jet/counter jet opening angle (∼2 • ; Wassell et al. 2006).
Another possibility to explain the nature of the asymmetry is the presence of a local perturbation on the disk material.An interesting hypothesis is the presence of a pressure bump produced by a vortex originated from an unseen planetary or dwarf companion.Until now, efficient dust traps produced by an anticyclonic vortex have been presented as plausible explanations for large (mm) dust grains to be trapped in an azimuthal direction on the disk (van der Marel et al. 2013Marel et al. , 2018;;Pineda et al. 2019).These dust traps tend to create arc-like features similar to the observed ones in our ring-based model.Additionally, magnetohydrodynamical simulations conducted by Flock et al. (2017) show that a local pressure maximum inside the disk's dead zone favors the creation of vortex, which can cast non-axisymmetric shadows on the outer disk.
More recently, the simulations performed by Varga et al. (2020) suggest that a large scale vortex produced by a Rossby wave instability could be the cause of the asymmetry in HD 163296.It is important to mention that density enhancement is not enough to create a change in the brightness distribution of the ring profile.This is because the emission of the ring is expected to be optically thick.In the scenario proposed by Varga et al. (2020), the large scale vortex favors the production of small dust grains, which modify the local temperature profile of the disk and, therefore, produce an increment in the emission at the position of the vortex.
Finally, the observed asymmetry could be explained by an in-homogeneous distribution of dust above the mid-plane of the disk.This idea is supported by recent HST data of the outer disk structure.Rich et al. (2019Rich et al. ( , 2020) ) reported strong surface brightness variations at scales larger than 660 mas (66 au) on time scales lower than 3 months.These results suggest that the origin of the moving shadows is material located at distances smaller than 5 mas (0.5 au) from the central star.To produce shadows in the disk at large scales, the material must reside at 0.8 mas (0.08 au) above the mid-plane of the disk, assuming coplanarity of the inner and outer disk.These authors also reported the presence of two dipper events in 2018, probably caused by variations of the scale height of the inner rim.The material must reside at the inner 4.1 mas (0.41 au) and at a scale height above 3.7 mas (0.37 au), suggesting the presence of a dusty wind.
A theoretical dusty-wind model that lifts material above the mid plane has been proposed by Bans & Königl (2012); Ellerbroek et al. (2014).Those authors also suggested that such a model is an important candidate for the origin of the strong outflows, like the one present in HD 163296.This would support the existence of material ejected above the mid-plane of the disk, not homogeneously distributed, which might be linked to the variable structure that we observe.To conclude which of the aforementioned scenarios is more plausible to explain the asymmetry in HD 163296, more observations (to improve considerably our u-v plane) are required in addition with dedicated simulations of the object.

Conclusions
This work presents new near-infrared interferometric observations of HD 163296 taken with GRAVITY, complemented with archival PIONIER data.Our multi-epoch campaign allows us to characterize the asymmetric and variable inner structure of the target.For this purpose, we used a parametric model of an azimuthally modulated ring.This model reproduces the squared visibilities and closure phases of each one of the epochs analyzed.To test the limitations of our model, we also fitted the data with an off-centered Gaussian model and, we conducted image reconstruction.These additional model and images confirmed the asymmetry of the inner disk and its variability.The inclination and position angle of the disk found with our parametric models are in agreement with previous estimates.However, the changes in the size of the ring across the different epochs make us support that the disk does not have a sharp inner edge but a smooth brightness profile.
Due to the variability of the disk morphology, we hypothesize that the nature of the asymmetry is not caused by an illumination effect.More plausible explanations include the presence of a local perturbation, like a vortex, or the presence of ejected dust above the mid-plane of the disk.Our estimation of the temperature of the asymmetry favors the existence of a mixed population of carbon and silicate dust grains or, as previously suggested, the presence of refractory dust grains.New data taken with MATISSE add further evidence for the presence of a noncentrosymmetric structure over different angular scales across the H, K and L bands.To fully determine the nature of a such structure, it is necessary to combine several interferometric observations with different baselines and wavelengths.Due to the high variability of the source, it is critical to obtain data over short time-scales (less than a month) to properly combine them and being able to perform image reconstruction and more sophisticated parametric (radiative transfer) models to unveil the nature of the asymmetry.

Fig. 1 .
Fig. 1.Panels show a mean image of the best-fit azimuthally modulated ring model described in Sec. 3. Instrument and epochs are displayed on each panel.These mean images were created from the best-fit models of the individual wavelengths per epoch.All panels are plotted with the same colormap and normalized flux scale.The white contours in the images represent the 20%, 40%, 60%, 80% and 95% of the peak's emission.The white star in the center of the images traces the position of the star in the model.
displays a mean image of the best-fit model for each epoch.A comparison between the observables from the data and the ones recovered from our model can be consulted in Figures B.1 to B.10 included in Appendix B.
Fig. 2. The upper panels display the GRAVITY 2018 data with the observables extracted from the 2019 Ring model.The lower panels display the opposite case.Notice how the observables from one epoch are not reproduced by the model of the other one.This is more evident for the closure phases.Therefore, this test supports the presence of a variable-asymmetric structure at the inner disk in HD 163296.

Fig. 3 .
Fig. 3.The plot shows the positions of the emission's peak extracted from our Ring model and complemented with the position of the peak obtained from the MATISSE data and the ring model presented by Varga et al. (2020).

Table 1 .
PIONIER -Best-fit parameters of the azimuthally modulated ring model

Table 2 .
GRAVITY -Best-fit parameters of the azimuthally modulated ring model

Table 3 .
HD 163296 -Temperature of the asymmetry V 2 ∼ 0.1 for baselines above 100 m.This makes difficult to ensure the monotonic decrement of the visibility.Furthermore, visibilities in H−band appear to have a small rebound for the largest baselines above 200 m.Nevertheless, the H−band data lack of intermediate baselines.This limitation does not provide us with more robust estimates of the visibility profile.The differences in the disk' sizes reported in Sec.4.1 support that the inner structure of the disk does not have a sharp edge.In contrast, it appears that the inner disk is smooth, being the best-fit rings of our model the emission traced by the different interferometric arrays used (see also Sec. 5.2).

Table D .
1. PIONIER (2013 June/July) -Best-fit parameters of the parametric models Table D.4.GRAVITY (2018 June) -Best-fit parameters of the parametric models Article number, page 22 of 33 GRAVITY Collaboration (): J. Sanchez-Bermudez et al.: The GRAVITY Young Stellar Object Survey Table D.5.GRAVITY (2019 June) -Best-fit parameters of the parametric models