© GRAVITY Collaboration 2021

1 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 objectives of several current observing facilities. ALMA (ALMA Partnership 2015, Atacama Large Millimeter Arrray) and adaptive-optics assisted imagers like SPHERE (Spectro-Polarimetric High-contrast Exoplanet REsearch; 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 like gaps, rings, spirals, vortex, and shadows (Andrews et al. 2018; Long et al. 2018; Benisty et al. 2015, 2017; de Boer et al. 2016; Pohl 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 tounderstand 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 only made 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 (λ₀ ~ 1.65 μm) and K (λ₀ ~ 2.2 μm) bands (Lazareff et al. 2017; Anthonioz et al. 2015; GRAVITY Collaboration 2019). Recent reconstructed images of 15 Herbig stars in the H-band (Kluska et al. 2020) show that 60% exhibit disk asymmetries at a few au scale and 27% has a non-centrosymmetric morphology that could 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. ALMA DSHARP (Disk Substructures at High Angular Resolution 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, 2019; Pinte et al. 2018). Multi-epoch 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, 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 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 the instrument CLIMB (CLassic Interferometry with Multiple Baselines beam combiner) of the CHARA (Center for High Angular Resolution Astronomy) array. 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 (Precision Integrated-Optics Near-infrared Imaging ExpeRiment; 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 on 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 (Astronomical Multi-BEam combineR; λ₀ = 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. (2021) detected an asymmetric structure using MATISSE (the Multi AperTure mid-Infrared SpectroScopic Experiment) data in the L-band (λ₀ = 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 the K-band and archival PIONIER 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²) and closure phases (CPs). The manuscript is divided as follows: Sect. 2 presents the interferometric data used for this work; Sect. 3 describes our model and Sect. 4 provides our results; we present our discussion in Sect. 5 and our conclusion in Sect. 6.

2 Observations and data reduction

2.1 GRAVITY data

HD 163296 was observed with GRAVITY (Eisenhauer et al. 2008, 2011; Gravity Collaboration 2017) during 2018 and 2019 as part of the 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)¹. 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 function. Before analyzing the data, all squared visibilities (V²) with a signal-to-noise ratio (S/N) 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 re-binned from the native ~1700 channels of the high-resolution data sets down to eight 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 S/N statistics of the different data sets. We obtained average error bars of = 0.01 and σ_CPs = 0.3° for the 2019 data and = 0.01 and σ_CPs = 0.8° for the 2018data, respectively. Our analysis is based on the exploitation of data from the science channel of GRAVITY and not from the fringe-tracker channel of the instrument. Detailed analyses of emission lines like Brγ were discarded and are the subject of a forthcoming paper. GRAVITY u-v planes are included in Fig. A.1.

2.2 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 in the Optical Interferometry Data Base (OIDB²) 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 (λ₀ = 1.65 μm) in 2013, and 2017, respectively. The log of these observations are also listed in Table A.1.

3 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 Fourier transform (1)

where, to account for the possible elongation of the ring, the original u − v coordinates were rotated and inclined in the form (2) (3)

where, 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. The polar form of the spatial frequencies u_r -v_r is given by q exp (iψ), a_r is the ring angular radius, and ρ_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 was convolved with a Lorentzian kernel of the form (4)

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 (5)

where, F_c, F_s, and F_h are the flux contributions of the ring (denoted with the subscript “c”), the point-like object (denoted with the subscript “s”), and the over-resolved component (denoted with the subscript “h”). We explicitly constrain the model to have F_s + F_c + F_h = 1.

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₁ and s₁ 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 data sets were combined into a single epoch. Since this corresponds to the same data set 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 Tables 1 and 2. Figure 1 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 Figs. B.1 to B.10.

4 Results

4.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:

• 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 from 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 data set. 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 ( 1.78 mas; 1.19 mas and; 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, PA = 146°).

  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 IRDIS data with SPHERE (Spectro-Polarimetric High-Contrast Exoplanet REsearch) taken on 2016, 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 and 60%, and F_h contributes between 0 and 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. The values of the coefficients c₁ and s₁, which trace the modulation of the ring, change 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. We note 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 at least on temporal scales of one year for the GRAVITY data. These timescales are restricted by the temporal sampling of our data. However, we cannot exclude changes in the structure at shorter timescales and/or the possibility of having observed different bright structures instead of a persistent single one.

  The presence of an asymmetry and hints in the variability have been previously suggested by several authors in the literature. The two most recent are Kobus et al. (2020) and Varga et al. (2021). 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 closure phases and squared visibilities simultaneously. Therefore, they cannot prove the asymmetric nature of the source. The second work uses MATISSE data and parametric models to constrain the source in the 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 only traces the asymmetric structure of the extended emission at one position angle. In contrast, our findings go one step further in confirming the asymmetry of the source, its variable nature, and persistence (on a temporal baseline of seven years).

4.2 Limitations of the azimuthally modulated ring model

While the model presented in Sect. 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 PIONIER and GRAVITY. Similar ring models have been presented in the literature (Lazareff et al. 2017; Kluska et al. 2020; GRAVITY Collaboration 2019; Varga et al. 2021) 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) 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 a similar degree of accuracy (see the reduced χ² values reported in Tables 1 and 2 and Tables D.1–D.5). 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 a 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 continuous and monotonically asymptotic 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 versus 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 trace 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 included 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 better reproduce the visibility function of the target. In particular, those authors found problems in reproducing 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 in reproducing 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 that 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 as large as 0.1 for baselines above 100 m. This makes it difficult to ensure the monotonic decrement of the visibility.

  Furthermore, visibilities in the H-band appear to have a small rebound for the largest baselines above 200 m. Nevertheless, the H-band data lack intermediate baselines. This limitation does not provide us with more robust estimates of the visibility profile. The differences in the disk sizes reported in Sect. 4.1 support the idea that the inner structure of the disk does not have a sharp edge. In contrast, it appears that the inner disk is smooth. Therefore, the best-fit rings just correspond to the emission seen by the different interferometric arrays used for our observations (see Sect. 5.2).

Comparison with reconstructed images

our geometrical models reproduce the observed profiles in our data. However, those geometrical models can only explorea 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 data sets. The complete imaging process is described inAppendix 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 coverage, 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, it is possible to improve our knowledge of the asymmetric morphology of HD 163296.

From the best images obtained, we could not favor the Ring or 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 the idea 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.

Fig. 1 Best-fit azimuthally modulated ring model. Panels show a mean image of the best-fit azimuthally modulated ring model described in Sect. 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.

Fig. 2 Best-fit azimuthally modulated ring model. Upper panels: GRAVITY 2018 data with the observables extracted from the 2019 Ring model. Lower panels: opposite case. We note 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.

5 Discussion

5.1 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 that 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 expression (van Buren & McCray 1988) (6)

where L_*,38 is the UV luminosity of the star in units of 10³⁸ erg s⁻¹, 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. We computed T_d considering a power-law dust size distribution taken from Mathis et al. (1977) for dust sizes ranging from 10⁻² to 1 μm. To calculate the distribution of temperatures per epoch, we extracted 10⁴ different samples of the peak’s positions and dust grain sizes, and we set L_*,38 = 6.38×10³⁴ erg s⁻¹ (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 note that the reported range of temperatures only traces the material observed with the different interferometric data sets. Therefore, those temperature values do not necessarily 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 (that 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 considerably 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% versus 15 ± 10%). Similarly, the ring models reported by Varga et al. (2021) support the fact 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) that 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 being responsible for the thermal emission of the variable disk.

5.2 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. A.1). With larger baselines, the 2017 configurations filter out most of the extended emission of the source and, therefore, a more compact object is observed by the interferometer, the 2017 August epoch being the one with the smallest estimate of the ring size. These results indicate that our observations do not trace a sharp edge of the inner ring structure. In contrast, they support the existence of a smoother inner morphology of the disk.

From our best-fit models, the peak of the emission in the ring changes for each of the different epochs. Figure 3 displays the projected positions of the peak in the ring emission for PIONIER and GRAVITY, complemented with the peak’s position from the ring model applied to the MATISSE observations described in Varga et al. (2021). Unfortunately, the amount of data we have is too limited 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 seven years our observations span. 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 that we are 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 (millimeter) dust grains being trapped in an azimuthal direction on the disk (van der Marel et al. 2013, 2018; Pineda et al. 2019). These dust traps tend to create arc-like features similar to those observed in our ring-based model. Additionally, magneto-hydrodynamical 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. (2021) 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. (2021), the large scale vortex favors the production of small dust grains that 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 inhomogeneous distribution of dust above the mid-plane of the disk. This idea is supported by recent Hubble Space Telescope (HST) data of the outer disk structure. Rich et al. (2019, 2020) reported strong surface brightness variations at scales larger than 660 mas (66 au) on timescales lower than three 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) and Ellerbroek et al. (2014). Those authors also suggested that such a model is an important candidate for the origin of 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 to dedicated simulations of the object.

Fig. 3 Positions of the emission’s peak asymmetry in the Ring model. 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. (2021).

6 Conclusions

This work presents new near-infrared interferometric observations of HD 163296 taken with GRAVITY and 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. This 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 the belief 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 illuminationeffect. 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 non-centrosymmetric 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 timescales (less than a month) in order to properly combine them, and we need to be able to perform image reconstruction and more sophisticated parametric (radiative transfer) models to unveil the nature of the asymmetry.

Acknowledgements

The authors thank the anonymous referee for his/her valuable comments. J.S.B. acknowledges the full support from the UNAM PAPIIT project IA 101220 and from the CONACyT “Ciencia de Frontera” project 263975. This research has made use of the Jean-Marie Mariotti Center OiDB service available at http://oidb.jmmc.fr. A.C.G. acknowledges funding from the European Research Council under Advanced Grant No. 743029, Ejection, Accretion Structures in YSOs (EASY). A.A., M.F. and P.G. were supported by Fundação para a Ciência e a Tecnologia, with grants reference UIDB/00099/2020 and SFRH/BSAB/142940/2018. R.G.L. acknowledges support from Science Foundation Ireland under Grant No. 18/SIRG/5597.

Appendix A Interferometric observations (Log and u-v planes)

HD 163296: (u-v) planes

Fig. A.1 u-v sampling of the different instruments and epochs used for this study. Blue dots indicate the spatial frequencies sampled while violet ones indicate their complex conjugate.

Appendix B Interferometric observables and best-fit azimuthally modulated ring model

HD 163296 Ring Model (2013 June/July, V²)

Fig. B.1 Observations vs. Ring model for the 2013 run. Panels display the squared visibilities (black dots) from the data vs. spatial frequency. Each panel corresponds to a different wavelength (see labels on the panels). The synthetic observables obtained from the best-fit Ring model are over-plotted with colored dots in the different panels.

HD 163296 Ring Model (2013 June/July, CPs)

Fig. B.2 Observations vs. Ring model for the 2013 run. Panels display the closure phase (black dots) from the data vs. spatial frequency. Each panel corresponds to a different wavelength (see labels on the panels). The synthetic observables obtained from the best-fit Ring model are over-plotted with different colors in the panels.

HD 163296 Ring Model (2017 April, V²)

Fig. B.3 Observations vs. Ring model for the 2017 April run. Panels are plotted as in Fig. B.1.

HD 163296 Ring Model (2017 April, CPs)

Fig. B.4 Observations vs. Ring model for the 2017 April run. Panels are plotted as in Fig. B.2.

HD 163296 Ring Model (2017 Aug., V²)

Fig. B.5 Observations vs. Ring model for the 2017 August run. Panels are plotted as in Fig. B.1.

HD 163296 Ring Model (2017 Aug., CPs)

Fig. B.6 Observations vs. Ring model for the 2017 August run. Panels are plotted as in Fig. B.2.

HD 163296 Ring Model (2018 July, V²)

Fig. B.7 Observations vs. Ring model for the 2018 run. Panels are plotted as in Fig. B.1.

HD 163296 Ring Model (2018 July, CPs)

Fig. B.8 Observations vs. Ring model for the 2018 run. Panels are plotted as in Fig. B.2.

HD 163296 Ring Model (2019 July, V²)

Fig. B.9 Observations vs. Ring model for the 2019 run. Panels are plotted as in Fig. B.1.

HD 163296 Ring Model (2019 July, CPs)

Fig. B.10 Observations vs. Ring model for the 2019 run. Panels are plotted as in Fig. B.2.

Appendix C Off-centered Gaussian model

This model corresponds to a point-like object (denoted with the subscript "s"), an elongated Gaussian disk (denoted with the subscript "c"), and an over-resolved component (denoted with the subscript "h"), defined by the equation (C.1)

The point-like component is centered at the phase reference and the Gaussian could be displaced by a given amount Δx and Δy in Right Ascension (R. A.) and Declination (Dec.), respectively. The flux contributions between the different component are F_s, F_c and F_h, the spatial frequencies sampled with the interferometer are (u,v), and G_c(u_r, v_r) is the Fourier transform of a Gaussian which is equal to (C.2)

where Θ_FWHM is the full-width-at-half-maximum (FWHM) of the Gaussian. The u_r and v_r components are defined as in Eq. 2. Tables D.1 to D.5 display the best-fit parameters obtained for this model.

Appendix D Imaging

We used BSMEM (BiSpectrum Maximum Entropy Method; Buscher 1994; Lawson et al. 2004) to image our target. Each wavelength was imaged independently. BSMEM employs a maximum-entropy algorithm to recover the images. The code works simultaneously with CPs, V², and triple amplitudes for the reconstruction. In this case, only CPs and V² were extracted from our data, while triple amplitudes were extrapolated directly by BSMEM from the previous two observables. The code uses a gradient descent method to perform a regularized minimization of the general form (D.1)

where x is the most-probable image sought that reproduces our data, χ² is the negative log-likelihood between the synthetic observables retrieved from a model image and our data, R(x) is the prior term that includes the information that is known about the source, and α is an hyperparameter that weights between the likelihood and the prior. BSMEM uses the Gull-Skilling entropy as prior: (D.2)

where x_n is the value of the image at pixel n and is the initial image, which is recovered in the absence of data. This regularizer enforces positivity in the image and favors smooth extended structures over sparse-compact ones; BSMEM automatically selects the hyper-parameter value. Images recovered here uses a pixel scale of 0.1 mas and they are recovered over a pixel grid of 501 × 501 pixels. To not over-regularize the reconstruction with the initial conditions, the starting image was set to be a Gaussian (FWHM = 4 mas) centered in the middle of the pixel grid. Images are presented as recovered from BSMEM and no subsequent smoothing convolution was done. All the reconstructions converged, however, we noticed a significant increased in the χ² for the 2018 epoch, compared with the other four epochs. Still, the observables of the 2018 epoch are well reproduced by the images. Additionally to the recovered images obtained from the data, we also reconstructed images using observables generated from the best-fit Ring and Gaussian models, while keeping the S/N statistics of the data. These images serve to compare the morphology observed in the reconstructed images and the one obtained from the parametric models, after being recovered using the same setup and imaging algorithm. The images from the data per wavelength and the images from the parametric models can be consulted in Figs. D.1 to D.10. The χ² of the images, the hyper-parameter value, and the number of iterations of the reconstructions are listed in Tables D.1 to D.5

The reconstructed images from the observables extracted from the best-fit parametric models are quite similar to the reconstructed images from the data. Our imaging experiment does not show a clear preference over the two parametric models to discriminate which one could better reproduce the image from the data. While on average residuals remain at 10% of the peak in the image, we noticed that for some spectral channels we observed residuals as large as 20% when we compare the reconstructed images and the best-fit models. We suspect that this difference is caused by asymmetric structures in the morphology of the target not traced by our parametric models. To have more quantitativeestimates on these structures, a compelling u-v coverage is required to do a high-fidelity image of the source.

PIONIER (2013 June/July) reconstructed Images (Gaussian model).

Fig. D.1 BiSpectrum Maximum Entropy Method (BSMEM) reconstructed images for our 2013 PIONIER epoch and the best-fit Gaussian model. The first column displays maps of the best-fit parametric models and the second column corresponds to the images recovered from the data at the corresponding epoch and wavelength. The third column shows reconstructed images from synthetic data using the indicated best-fit model. The fourth column displays the residuals between thereconstructed images from the data and the ones from the best-fit models. The wavelength for each row is labeled on the panels of the first column. The white ellipse in the second and third column corresponds to the synthesized beam (at FWHM). The white contours in the images of the second and third columns correspond to 20%, 40%, 60%, 80%, and 95% of the image’s peak. The maximum and minimum values labeled in the fourth column show the relative percentage of the (positive and negative) residuals and the peak relative to the reconstructed image from the data.

PIONIER (2013 June/July) reconstructed images (Ring model)

Fig. D.2 BSMEM reconstructed images for our 2013 PIONIER epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.

PIONIER (2017 April) reconstructed images (Gaussian model)

Fig. D.3 BSMEM reconstructed images for our 2017 (April) PIONIER epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.

PIONIER (2017 April) reconstructed images (Ring model)

Fig. D.4 BSMEM reconstructed images for our 2017 (April) PIONIER epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.

PIONIER (2017 August) reconstructed images (Gaussian model)

Fig. D.5 BSMEM reconstructed images for our 2017 (Aug.) PIONIER epoch and the best-fit Gaussian model. Panels are plotted as in Fig. D.1.

PIONIER (2017 August) reconstructed images (Ring model)

Fig. D.6 BSMEM reconstructed images for our 2017 (Aug.) PIONIER epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.

GRAVITY (2018 July) reconstructed images (Gaussian model)

Fig. D.7 BSMEM reconstructed images for our 2018 GRAVITY epoch and the best-fit Gaussian model. Panels are plotted as in Fig. D.1.

GRAVITY (2018 July) reconstructed images (Ring model)

Fig. D.8 BSMEM reconstructed images for our 2018 GRAVITY epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.

GRAVITY (2019 July) reconstructed images (Gaussian model)

Fig. D.9 BSMEM reconstructed images for our 2019 GRAVITY epoch and the best-fit Gaussian model. Panels are plotted as in Fig. D.1.

GRAVITY (2019 July) reconstructed images (Ring model)

Fig. D.10 BSMEM reconstructed images for our 2019 GRAVITY epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.

References

All Tables

All Figures

Fig. 1 Best-fit azimuthally modulated ring model. Panels show a mean image of the best-fit azimuthally modulated ring model described in Sect. 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.

In the text

Fig. 2 Best-fit azimuthally modulated ring model. Upper panels: GRAVITY 2018 data with the observables extracted from the 2019 Ring model. Lower panels: opposite case. We note 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.

In the text

	Fig. 3 Positions of the emission’s peak asymmetry in the Ring model. 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. (2021).
In the text

	Fig. A.1 u-v sampling of the different instruments and epochs used for this study. Blue dots indicate the spatial frequencies sampled while violet ones indicate their complex conjugate.
In the text

	Fig. B.1 Observations vs. Ring model for the 2013 run. Panels display the squared visibilities (black dots) from the data vs. spatial frequency. Each panel corresponds to a different wavelength (see labels on the panels). The synthetic observables obtained from the best-fit Ring model are over-plotted with colored dots in the different panels.
In the text

	Fig. B.2 Observations vs. Ring model for the 2013 run. Panels display the closure phase (black dots) from the data vs. spatial frequency. Each panel corresponds to a different wavelength (see labels on the panels). The synthetic observables obtained from the best-fit Ring model are over-plotted with different colors in the panels.
In the text

	Fig. B.3 Observations vs. Ring model for the 2017 April run. Panels are plotted as in Fig. B.1.
In the text

	Fig. B.4 Observations vs. Ring model for the 2017 April run. Panels are plotted as in Fig. B.2.
In the text

	Fig. B.5 Observations vs. Ring model for the 2017 August run. Panels are plotted as in Fig. B.1.
In the text

	Fig. B.6 Observations vs. Ring model for the 2017 August run. Panels are plotted as in Fig. B.2.
In the text

	Fig. B.7 Observations vs. Ring model for the 2018 run. Panels are plotted as in Fig. B.1.
In the text

	Fig. B.8 Observations vs. Ring model for the 2018 run. Panels are plotted as in Fig. B.2.
In the text

	Fig. B.9 Observations vs. Ring model for the 2019 run. Panels are plotted as in Fig. B.1.
In the text

	Fig. B.10 Observations vs. Ring model for the 2019 run. Panels are plotted as in Fig. B.2.
In the text

Fig. D.1 BiSpectrum Maximum Entropy Method (BSMEM) reconstructed images for our 2013 PIONIER epoch and the best-fit Gaussian model. The first column displays maps of the best-fit parametric models and the second column corresponds to the images recovered from the data at the corresponding epoch and wavelength. The third column shows reconstructed images from synthetic data using the indicated best-fit model. The fourth column displays the residuals between thereconstructed images from the data and the ones from the best-fit models. The wavelength for each row is labeled on the panels of the first column. The white ellipse in the second and third column corresponds to the synthesized beam (at FWHM). The white contours in the images of the second and third columns correspond to 20%, 40%, 60%, 80%, and 95% of the image’s peak. The maximum and minimum values labeled in the fourth column show the relative percentage of the (positive and negative) residuals and the peak relative to the reconstructed image from the data.

In the text

	Fig. D.2 BSMEM reconstructed images for our 2013 PIONIER epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.
In the text

	Fig. D.3 BSMEM reconstructed images for our 2017 (April) PIONIER epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.
In the text

	Fig. D.4 BSMEM reconstructed images for our 2017 (April) PIONIER epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.
In the text

	Fig. D.5 BSMEM reconstructed images for our 2017 (Aug.) PIONIER epoch and the best-fit Gaussian model. Panels are plotted as in Fig. D.1.
In the text

	Fig. D.6 BSMEM reconstructed images for our 2017 (Aug.) PIONIER epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.
In the text

	Fig. D.7 BSMEM reconstructed images for our 2018 GRAVITY epoch and the best-fit Gaussian model. Panels are plotted as in Fig. D.1.
In the text

	Fig. D.8 BSMEM reconstructed images for our 2018 GRAVITY epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.
In the text

	Fig. D.9 BSMEM reconstructed images for our 2019 GRAVITY epoch and the best-fit Gaussian model. Panels are plotted as in Fig. D.1.
In the text

	Fig. D.10 BSMEM reconstructed images for our 2019 GRAVITY epoch and the best-fit Ring model. Panels are plotted as in Fig. D.1.
In the text