© ESO, 2016

1. Introduction

Stars showing the B[e] phenomenon are defined by the presence of permitted (Balmer and FeII) and forbidden emission lines along with a large infrared excess indicating circumstellar (CS) dusty material. Lamers et al. (1998) classified these stars in five categories including supergiants, pre-main-sequence stars, planetary nebulae, symbiotic stars, and unclassified stars. The diversity in the nature of stars showing the B[e] phenomenon is puzzling. It is still not clear if this phenomenon is caused by the same physical process because only a few objects have been studied in depth. Optical interferometry is a powerful tool to investigate the close environment of these stars (Weigelt et al. 2011; Wheelwright et al. 2012, 2013; Vural et al. 2014) and in some cases reveal their binarity (e.g. Millour et al. 2009; Kraus et al. 2012; Wang et al. 2012). Such a secondary component can disturb the environment of the primary star causing forbidden lines to appear. The most studied B[e] stars, however, remain unclassified because of their unknown and strongly debated evolutionary status. Some of these stars show hints of a pre-main-sequence nature but are not associated with any star-forming region. Studying the close environment of these massive and puzzling B[e] stars may bring new insights into our understanding of the baffling phenomena occurring around them.

MWC158 (also called HD 50138 and V743 Mon) is one of the brightest B[e] stars (see Table. 1) even though it remains unclassified (e.g. Lamers et al. 1998; Ellerbroek et al. 2015, hereafter E15). It was previously argued to be a young object (Morrison & Beaver 1995) and an evolved object in the main sequence or close to its ending (Borges Fernandes et al. 2009). Its lack of membership to a star-forming region (Maddalena et al. 1986) is the main obstacle to the classification of this object as a pre-main-sequence star.

MWC158 shows strong line profile variability (from years to hours, Hutsemekers 1985; Andrillat & Houziaux 1991), which has been linked to shell ejections or formation of a one-arm, ejection spiral in its environment (Borges Fernandes et al. 2009) that could indicate an evolved nature. Pogodin (1997) has reported variability of H_α, helium and sodium lines on timescales from days to months. Hutsemekers (1985) has reported variability in the UV that is associated with an outburst or a shell phase (Andrillat & Houziaux 1991). Recently, variability was reported in the visible spectra with timescales ranging from days (interpreted as stellar pulsations by Borges Fernandes et al. 2012) to months (linked with its environment, e.g. Borges Fernandes et al. 2009; E15). The discussion on the variability of this object has been based on spectrophotometry and until now there have been no spatially resolved constraints.

The disk-like morphology of the CS environment was first indicated by polarimetric and spectro-polarimetric observations (Vaidya et al. 1994; Bjorkman et al. 1998; Oudmaijer & Drew 1999; Harrington & Kuhn 2007) before they were confirmed by interferometric observations. These observations were performed in the near- and mid-infrared focusing on the size and orientation of the object (Monnier et al. 2009; Borges Fernandes et al. 2011). The inclination is found to be around 50°−55°degrees and the position angle (PA) around 70°. These studies confirmed the presence of a disk morphology but could not explore its shape in detail.

Here we report the analysis of a new and rich dataset taken with PIONIER/VLTI and spanning different epochs. It allows us to image the inner parts of the disk in more detail and to look for morphological variability. We present our observations in Sect. 2. In Sect. 3 we describe the image reconstructions on the interferometric data. In Sect. 4 we describe the parametric modelling of the interferometric data to take into account the observed morphological variability. We discuss our results in Sect. 5 and we conclude in Sect. 6.

2. Observations

We performed PIONIER/VLTI observations of MWC158 to reach high angular resolution in the near-infrared continuum. The PIONIER instrument (Le Bouquin et al. 2011) combines four telescopes in the H-band and provides six squared visibility (V²) measurements and four closure phases (CP). If the object is resolved well, the AT relocation capability can be used to pave the { u,v }-plane efficiently and enable image reconstruction. The preparation of observations was carried out with ASPRO2 (Bourgès et al. 2013). The calibrators were found using SearchCal (Bonneau et al. 2011, HD 46931, HD 47559, HD 50639, HD 51914, HD 52265). The data was reduced with pndrs (Le Bouquin et al. 2011).

The log of observations is shown in Table A.1. These observations were taken during the PIONIER Science Verification programme, the PIONIER Herbig Ae/Be Large Programme (190.C-0963, PI: Berger) and during 089.C-0211 (PI: Benisty).

Fig. 1 Temporal representation on data taken on MWC158 during 3 epochs.

The temporal distribution of all the observations is shown in Fig. 1. The target was observed in October and December 2010 (epoch A), April (epoch B) and December 2012, and January and February 2013 (epoch C).

Fig. 2 { u,v }-plane for the 3 epochs. Epoch A is indicated in red triangles, epoch B in orange circles, and epoch C in blue diamonds.

Fig. 3 Temporal representation on data taken on MWC158. The colours correspond to the wavelength. Dataset in V2 and CP. We notice strong differences in the CP between the epochs that may indicate variability.

Epoch A was observed with the large spectral dispersion mode (seven channels across the H-band), while data of epochs B and C were observed with the small spectral dispersion mode (three channels across the H-band). The { u,v }-planes of the three epochs correspond to a maximum resolution of 2.6, 9.6, and 2.3 mas, respectively, (with 72, 12, and 90 different baselines, respectively, and 504, 36, and 270 (u, v)-points, respectively) and are presented in Fig. 2. The associated datasets are presented in Fig. 3. A global V² decrease with spatial frequency indicates that the object is resolved. The flat asymptote in V² at the highest spatial frequencies indicates that the object is sufficiently resolved to discriminate the stellar (point-like) flux from the CS environment. In parallel, we see that at a given baseline there is a decrease in V² with increasing wavelengths. This is because of the chromatic effect that arises when there is a temperature difference between the star and its CS environment (Kluska et al. 2014). For the smallest baselines, we see a drop in V² that could be because of the presence of an extended emission.

The CP signal rises significantly above 15 degrees on epochs A and C, which have sufficiently long baselines to resolve any departure from point symmetry in the circumstellar environment. The CP signal seems to be different from epoch to epoch.

3. Image reconstruction

Fig. 4 Chromatic parameters determination for epoch A (red) and C (blue) at 1σ (solid lines) and 3σ (dashed lines). The red (epoch A) and blue (epoch C) points are the best values for the chromatic parameters. These values were used for the image reconstruction process. The vertical dashed line indicates the spectral index of the central star, which equals − 4 because the star is assumed to be in the Rayleigh-Jeans regime.

We first use an image reconstruction technique to retrieve a model-independent object morphology. Because image reconstruction needs a sufficiently well covered { u,v }-plane, we perform image reconstruction on the PIONIER datasets on epochs A and C only. The image reconstruction process requires a careful choice of the regularisation and its weight, which is described in Appendix D.

The images were made using the SPARCO approach (Kluska et al. 2014) together with the MiRA algorithm (Thiébaut 2008). In the SPARCO approach, the object is modelled as an unresolved central star plus an extended environment. Both components have different spectral behaviours and therefore different spectral indexes. The algorithm provides an image of the extended environment for a given stellar-to-total flux ratio at 1.65 μm, , and the spectral index of the environment, d_env (see Eq. (4) in Kluska et al. 2014). It is possible to derive these chromatic parameters from the SED of the object but the variability of the object and the lack of simultaneous photometry prevent us from doing so. We therefore use the interferometric data itself to estimate these chromatic parameters. Several maps of the environment are reconstructed for a large number of pairs of parameters, and d_env, sampled on a 100 × 100 grid ( from 0.05 to 0.25 for epoch A and 0.1 to 0.35 for epoch C, and d_env from 0 to 10 for epoch A and from −2 to 7 for epoch C). We compute the a posteriori image likelihood for each point in this grid (see Appendix D). We use different grids because we see that the ranges of chromatic parameters we need to sample are different. The grids must be large enough to have negligible likelihood values at their edges. We use bilinear interpolations to compare the grids. We thus choose the optimal parameters with the highest a posteriori image likelihood. The result is plotted in Fig. 4.

We use marginalised probability distributions to obtain the error bars on each parameter (see Fig. D.1). We notice that the result for the two epochs is different. Indeed, epoch A points to a value of the stellar-to-total ratio of = 11.9% ± 1.0% and a spectral index of the environment of d_env = 5.0 ± 0.8 ( K), while epoch C gives = 17.5% ± 2.8% and a spectral index of the environment of d_env = 1.6 ± 0.9 ( K). These results confirm the variable nature of the object and imply a change in the environment temperature of more than 200 K. The probability that the two epochs have the same chromatic parameters is less than 2%. Also, the stellar-to-total flux ratio is changing from 11.9% to 17.5% and agrees with the variability seen in photometric measurements (Lazareff, priv. comm.). We want to investigate the origin of this variability by performing image reconstruction on both epochs. We used the values with the highest probabilities = 11.8% and 16.7% and d_env = 5.1 and 1.8 for epoch A and C, respectively. To see the differences between the images at other flux ratios () and spectral indexes (d_env), we show image reconstruction with different pairs of chromatic parameters in Fig. D.2. No significant differences are seen in the image over the range of allowed chromatic parameters. More importantly, the departures from point symmetry in the image are not affected.

Fig. 5 Image reconstructions in H-band on PIONIER data, epochs A (left) and C (right). The colour bar on the right was redefined to have the maximum flux equal to one for the sake of readability. The beam for each epoch is plotted in blue in the bottom right corner of the images. The contours are at 1σ (dashed line), 3, and 5σ (solid lines) significance levels.

3.1. Image analysis

We apply a bootstrap method (e.g. Efron & Tibshirani 1994) on the baselines and telescope triplets used to compute the CP to determine the reliable features in the images. This consists of a random draw of the data points within the whole dataset. The generated dataset has the same number of points as the original dataset. Some data points can be drawn several times whereas some of them are not selected at all. This allows the effects of an imperfect { u,v }-plane to be revealed. This process was carried out 500 times and led to a per-pixel error estimation.

The mean images are presented in Fig. 5 with contours at 1 (dashed line), 3 and 5σ (solid lines). The of the images equal 2.0 for both epochs (see Fig. D.4).

The two images reveal an elongated morphology, detected at more than 5σ with a major axis of ≈6 mas and a minor axis of ≈3 mas. The PA of the major axes are similar for the two epochs and are in good agreement with previous interferometric measurements (Borges Fernandes et al. 2011). The differences mainly result from the relative positions between the star and the maximum brightness of the environment. While in epoch A the image seems to be symmetric with respect to the star, this is not the case in epoch C. At this epoch the environment is elongated in the north-east direction indicating a strong asymmetry detected at more than 20σ. Moreover, the epoch C image has an additional feature appearing at the south-west direction, which is detected at greater than 5σ.

3.2. Confirmation of the variability in Fourier space

Fig. 6 Image reconstruction fit on datasets of epoch C superimposed on dataset from epoch A (left) and inversely (right). The black squares represent the data and the blue triangles represent the images. The residuals normalised to the standard deviation are plotted in the bottom of each graph.

These two images appear very different from each other. The differences arise from the environment brightness distribution and from the stellar-to-total flux ratio. It is difficult to conclude the variability of an object if we do not measure precisely the same points in the { u,v }-plane. From the images we are able to extrapolate the values of the visibility through the entire Fourier plane. We can therefore compare the visibilities from the image at one epoch to the dataset on the other epoch. The superimposition of the data from the image of epoch C (A respectively) to the dataset of epoch A (C resp.) are in Fig. 6.

We can see that the differences in the squared visibilities are generally between 5 and 10σ (with peaks up to 60σ). The CP differences are smaller at baselines shorter than 80 m that correspond to an angular resolution of 4 mas. At longer baselines the differences increase significantly to 100σ. This is because of the elongations seen in the two epochs that appear at sizes smaller than 4 mas. The variability of the environment is therefore highly significative.

3.3. Asymmetry analysis

Fig. 7 Positive part of the asymmetric maps (Asym) of epochs A (left) and C (right). The orientation of the major axis (PA of 71 degrees as estimated in Borges Fernandes et al. 2011) is plotted in black in the image.

To emphasise the deviations from centro-symmetry in a model-independent way, we construct an asymmetric map (Asym). It is obtained taking the positive part of (1)where Img is the reconstructed image and Img₁₈₀ the same image rotated by 180 degrees. The asymmetry map is the flux in the image that produces the non-zero CP. It is the amount of flux that is not point-symmetric.

In Fig. 7 we plot the positive part Asym for both epochs.For epoch A, the positive parts of the asymmetry are on the southern side of the major axis (PA of 71 degrees as estimated in Borges Fernandes et al. 2011). This is what is expected for an inclined inner disk rim, for instance, where we see the further side of the rim because the closer part is obscured by the disk itself. This is not the case in epoch C anymore where the position of the asymmetric flux is close to the major axis. In both images the positions of the asymmetric flux seem to be at the same distance from the star. That could be an indication of rotation of the asymmetric feature around the central star. We use parametric models to quantify this spatial variability.

4. Modelling the disk with asymmetric emission

The analysis of the reconstructed images at different epochs leads us to the hypothesis that there could be a moving feature in the disk. We refer to this feature as a bright spot. We explore this hypothesis with a parametric model. Because parametric models do not require a { u,v }-plane that is as rich as for an image, we can divide the epochs into sub-epochs to better sample the trajectory of the asymmetry. We therefore divided the epochs A and C into the following sub-epochs (see Appendix A): A1 (1 night), A2 (1 night), A3 (2 nights), C1 (2 nights), C2 (3 nights), C3 (3 nights), and C4 (4 nights).

4.1. Model of a star, a disk, and an asymmetry

Our data supports the presence of a flattened morphology compatible with a disk in agreement with Borges Fernandes et al. (2011). We built a simple model of a star, disk, bright spot, and extended component. The extended component can be interpreted as stellar light scattered by the disk (see Pinte et al. 2008; Anthonioz et al. 2015). From the image reconstructions we do not see any evidence for a disk inner rim. We therefore modelled each component as a point source, inclined Gaussian disk, background flux, and second point source that can be shifted w.r.t. the primary star. This model is detailed in Appendix B. All the flux ratios are derived for λ₀ = 1.65 μm.

We used a Levenberg-Marquardt method to find the best parameters. We start from 2000 different initial solutions generated by a Monte Carlo method to avoid local minima. We estimate the errors on each parameter with the same bootstrap method already used in Sect. 3.1. We repeat it two thousand times in order to generate two thousand estimations of each parameter. We use the standard deviation on each parameter to estimate its error (see Fig. C.1).

4.2. Results

The best-fit parameters are presented in Tables 2 and 3. We find different values for the different epochs. Interestingly, the highest correspond to the observations with the highest spatial resolution (epochs A3 and C2). This indicates that our model is too simple and that we are still not catching all the complexity of the object on the smallest spatial scales.

We discuss the results of the fits for each parameter:

• Flux of the asymmetry: at epoch B, the flux is not significant enough (less than 2σ significance) because it is not resolved. For all the other epochs, the asymmetry has a flux that varies from 1.6% (epoch C4) to 20.5% (epoch C1).

• Flux of the star: the stellar flux ratio is oscillating around 30% with the highest value at epoch B (42.4%) and the lowest values at epochs A3 and C2. We notice that the epochs with the highest values have short baselines and the epochs with the lowest values have the longest baselines. The estimation of the stellar flux is dependent on the angular resolution achieved with respect to the size of the environment. If the angular resolution is sufficient to resolve the environment without resolving the star, then the stellar-to-total flux ratio is unbiased. In contrast, the flux of the central star can be considered the unresolved flux and include unresolved disk flux. Having simultaneous photometry would have helped to isolate the stellar contribution in the unresolved flux and claim a possible continuum variability.

• Flux of the Gaussian: the flux of the Gaussian that is representing the disk also varies. It seems to be correlated to the stellar flux variation for the epochs from B to C4. Interestingly the fluxes for epochs A1−A3 seem to be anti-correlated.

• Temperatures: the temperatures found for the disk are around 1400 K. We also notice that the temperatures are higher for epochs A2 and C1, which correspond to the observations with the less extended configurations, for which the environment is not entirely resolved. In these cases, the unresolved flux is also included in the stellar flux ratio (), thus the flux of the star cannot be considered anymore as in the Rayleigh-Jeans regime, as previously assumed, causing an overestimation of the disk temperature. These temperatures are quite different from those found by image reconstruction, especially for epoch A.

• Size of the Gaussian: the size of the Gaussian is stable (around 4.5 mas) but for epochs B, C1 and C2. Epochs B and C1 have the smallest angular resolution and it can affect the size determination. The best model for C2 epoch has the highest χ² and is not reproducing the data very well. We do not therefore consider the size variation as significative.

• Disk orientations: the inclination of the disk is consistent through the epochs with 60 degrees (a bit lower for epochs C3 and C4). Its position angle is consistent with a value between 55 and 60 degrees.

• Bright spot coordinates: the positions of the bright spot are different from epoch to epoch but they remain close to the disk. The bright spot seems to follow a trajectory around the central star that has the same orientation as the disk. It could therefore be a disk feature.

The variation of the stellar flux ratio, the disk temperature, and size seems to be dependent on the baseline length. This is due to a degeneracy. In our dataset the { u,v }-plane is not sufficient to resolve all of the circumstellar environment, the unresolved flux is therefore composed of the stellar flux and an unresolved circumstellar environment flux. In the fitting process the stellar-to-total flux ratio (), which is modelled by the unresolved stellar component, is therefore over estimated. Moreover the assumption that the model of the unresolved star is in the Rayleigh-Jeans regime is no longer valid. As the spectral index of the environment is defined relative to the index of the star, it is not properly estimated. As a result the temperature derived from it will be overestimated (for ex. 1987 K and 2417 K for epochs A2 and C1 respectively).

5. Discussion

A physical interpretation of the detected morphological variability is difficult from our interferometric dataset alone because of the insufficient temporal sampling of our observations, however, we present some possible interpretations.

5.1. A binary?

Fig. 8 Location of MWC158 primary star (black hexagon, see text for details) in the isochrones of Bressan et al. (2012) for two scenarios: the scenario in which MWC158 is young, from yellow (0.3 Myr) to red (0.5 Myr), and old, from blue (0.1 Gyr) to cyan (0.2 Gyr). The triangles indicate the position of 6 M⊙ as determined by Borges Fernandes et al. (2009) for the young track. For the old track, these positions are out of range of the plot. The dashed line indicates the upper limit on the absolute H magnitude of a hypothetical secondary of MWC158.

A binary can easily reproduce such a strong CP signal with rapid variability, however all of our binary fits alone (i.e. without a disk) are not able to reproduce the data (χ² of 15 to 250 for the different sub-epochs). Previous spectroscopic observational campaigns do not show any binary signal (Corporon & Lagrange 1999; Borges Fernandes et al. 2009, 2012) but they do not discard this possibility. From the parametric fit in Sect. 4.2, we saw that a fit with a star, disk, and bright spot is able to reproduce correctly most of the dataset, especially the closure phases.

With the flux ratios from the parametric fits, we can locate the primary star and the asymmetry (if related to a secondary) on a Hertzprung-Russel diagram. For the primary we used the effective temperature range derived by Borges Fernandes et al. (2009, 11 200 K < T_eff and the absolute H-magnitude range derived from the by our modelling. From isochrones of Bressan et al. (2012), we can see the presence of two scenarios for the primary star of MWC158 (see Fig. 8). A “young” scenario with ages from 0.3 to 0.5 Myr and an “old” scenario with ages from 0.1 to 0.2 Gyr. The masses associated with the “young” scenario range from 4.6 to 5.7 M_⊙ whereas the “old” scenario predicts masses from 3.7 to 4.8 M_⊙. These masses are in agreement with the literature (Borges Fernandes et al. 2009).Assuming that the flux of the asymmetry is composed of the flux of the secondary and from surrounding parts of the disk, we can estimate an upper limit for the secondary mass. If we take the minimum flux value for the asymmetry, which is 5% of the flux of the primary for epoch C4 (that corresponds to an absolute magnitude in H of 1.6), we can estimate the upper mass limit of the secondary using model isochrones. Assuming that the secondary would have the same age as the primary and that no mass exchange happened, the secondary would have a maximum mass of 0.6 M_⊙ in the “young” scenario and of 2.3 M_⊙ in the “old” scenario. The total mass of the binary system would be about 7 M_⊙. With such a mass and taking an orbital radius of 4.7 mas (the position of the spot of the C4 epoch), we expect a Keplerian period of almost one year (344 days). Looking at the sequence of C1−C4 epochs that lasts two months, however, we can see that the asymmetry changes its position on a shorter timescale than our prediction (we find the same conclusion for epochs A1−A3).

We notice that the coordinates of the asymmetry are distributed around the star and may trace an orbital trajectory around it. We tried to fit a circular orbit (see Appendix E) but it was unsuccessful. The trajectory of the spot could be more complex than a circular Keplerian orbit but we do not have good enough temporal sampling to characterise it more accurately. Another possibility is that our parametric model is too simple to recover all the complexity of the environment morphology.

5.2. An asymmetric environment?

From the parametric fit the positions of the bright spot are located close to the disk and this can be interpreted as an asymmetry lying in the disk. This is supported by a spectroscopic campaign over 20 yr on this object, which provides evidence for an asymmetric gaseous disk; this was carried out by Jeřábková et al. (2016). The authors could not find a clear periodical behaviour but they claimed several pseudo-periods ranging from years to months.

Several scenarios could explain such asymmetric features. Borges Fernandes et al. (2009) postulated the formation of a “one-armed spiral” due to matter ejection from a stellar hot spot. Asymmetric features were also observed in the circumstellar environment of the AGB star L₂ Puppis (Kervella et al. 2014; Lykou et al. 2015). Such a pattern is very likely generated by the interaction of a close companion (that we could not constrain around MWC158) with the strong stellar wind coming from the primary (van Winckel et al. 2009).

Another possible explanation for the observed asymmetry lies in shell ejection episodes. The photometric variability of this object has led several authors to postulate episodes of shell ejection for this star (Merrill 1931; Merrill & Burwell 1933; Doazan 1965; Hutsemekers 1985; Andrillat & Houziaux 1991; Bopp 1993; Pogodin 1997; Borges Fernandes et al. 2009, 2012; Jeřábková et al. 2016). The asymmetry could be created by these ejections that would carry a part of the gas and dust away. Such a phenomenon was already proposed by Bans & Königl (2012) with a disk wind model. The asymmetric shell ejection scenario could explain the fact that at epoch B the object was more extended that in epochs A and C. In order to reproduce our data the shell should be very asymmetric and its brightness cannot only evolve radially. In order to test such a possibility, the morphological variability of the environment of this object should have been observed more extensively to achieve a better temporal sampling.

6. Conclusions

The aim of these observations was to constrain the morphology of the environment of MWC158. The orientations and sizes that we derived are in clear agreement with what was found in previous studies (Borges Fernandes et al. 2011). We reconstructed two images at different epochs. For the first time we have detected a strongly asymmetric environment that is spatially evolving with time. We modelled the asymmetric environment with a simple model of a star, a disk and a bright spot. We found that the bright spot has different coordinates from epoch to epoch and seems to move around the central star. The result of fitting a circular orbit is inconclusive. It is either due to the lack of additional epochs or to the fact that the object is more complex than our modelling, as already suggested (Borges Fernandes et al. 2009; Jeřábková et al. 2016). Because it is located close to the star, it is likely that the bright spot is a disk asymmetry, however, the detected asymmetry cannot be linked directly with the presence of a secondary. It could be caused by spiral arms lying in the disk (Borges Fernandes et al. 2009; Jeřábková et al. 2016) or other disk features.

We do not exclude the scenario of an asymmetric shell ejection that could explain this data as well, but it needs to be further explored. The origin of the asymmetry could also be partly due to the circumstellar gas or an additional stellar companion that is not confirmed in our dataset. More observations with better temporal sampling are needed to assess the displacement of the asymmetry.

The environment of MWC158 cannot be considered as axisymmetric and static. Simultaneous interferometric observations of the continuum and circumstellar lines are necessary to constrain the dynamics in place in this object. This could be possible in the context of the Gravity instrument by observing the Br_γ line and the continuum simultaneously. Spatially resolving the kinematics will indicate whether the inner gaseous disk is disrupted in the same way as the dusty continuum. Nevertheless, to understand the variability of the environment, simultaneous spectrophotometric monitoring is needed to follow up the stellar-to-total flux ratio to correctly interpret the interferometric datasets. It would also be exciting to constrain the dust at larger scales, in order to detect the possible ejection shells or spirals caused by an inner companion, with direct imaging of the scattered light close to the star in the visible, for instance, with SPHERE/ZIMPOL.

In the context of B[e] stars, MWC158 is a key case because it is one of the brightest stars showing this phenomenon and it is relatively close. If the B[e] phenomenon is linked with ejections without being a binary, it will bring new elements to explain some stars showing the B[e] phenomenon.

Acknowledgments

We thank the anonymous referee for remarks that greatly improved the manuscript. We would like to thank Stefan Kraus for fruitful discussions. This work is supported by the French ANR POLCA project (Processing of pOLychromatic interferometriC data for Astrophysics, ANR-10-BLAN-0511). J.K. acknowledges support from an Marie Sklodowska-Curie CIG grant (Ref. 618910, PI: Stefan Kraus). F.S. was partly supported by the Sinergia project “Euclid” from the Swiss National Science Foundation. PIONIER is funded by the Université Joseph Fourier (UJF), the Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), the Agence Nationale pour la Recherche (ANR-06-BLAN-0421, ANR-10-BLAN-0505, ANR-10-LABX56), and the Institut National des Science de l’Univers (INSU PNP and PNPS). The integrated optics beam combiner is the result of a collaboration between IPAG and CEA-LETI based on CNES R&T funding. The authors want to warmly thank all the people involved in the VLTI project. This work is based on observations made with the ESO telescopes. It made use of the Smithsonian NASA Astrophysics Data System (ADS) and of the Centre de Données astronomiques de Strasbourg (CDS). All calculations and graphics were performed with the open source software Yorick. This research has made use of the Jean-Marie Mariotti Center ASPRO2 and SearchCal services co-developed by CRAL, IPAG, and FIZEAU.

References

Appendix A: Log of observations

Appendix B: Model of a star, disk with a “bright spot”, and an extended component

The star has a visibility of 1. The extended flux has a visibility of 0. The inclined Gaussian is geometrically described by its FWHM (w), its inclination (i), position angle (PA) as follows: (B.1)where V^g(ν′) is the visibility of the Gaussian and are the spatial frequencies from the { u,v }-plane points oriented with respect to the object such as with the spatial frequencies from the { u,v }-plane.

The bright spot is mainly described by its position x_bs and y_bs as follows: (B.4)where V^a(ν′′) is the visibility of the bright spot and ν′′ = u′′ + v′′ are the spatial frequencies from the { u,v }-plane points shifted with respect to the object such as The total visibilities are (B.7)with the stellar to total flux ratio, the extended-to-total flux ratio, the Gaussian-to-total flux ratio, and the bright spot-to-total flux ratios. To define these parameters, we have taken a reference wavelength λ₀ = 1.65 μm (centre of the H-band). We assumed the star () and the extended flux () being in the Rayleigh-Jeans regime ( and ) and that the Gaussian and the bright spot (, ) have their own temperatures so that their fluxes are scaled accordingly and . We assumed the temperature of the Gaussian and the bright spot to be equal. At λ₀, the sum of the stellar flux ratio (), extended flux ratio (), inclined Gaussian flux ratio (), and bright spot () equals unity.

Appendix C: Parameter error determination

In order to estimate the errors on the parameters from the parametric fit, we have generated 2000 new datasets from the original one adding a noise term. We assumed a Gaussian noise with a standard deviation equal to that from the data. We then have computed the mean value and the error on each parameter from the best fit to the generated dataset. The standard deviation for each parameter gave us the error at 1σ that we have used.

Fig. C.1 Histograms of best-fit values for all parameters of the parametric fit for each epoch. The red solid line represents the mean values and the red dashed lines represent the 1σ errors.

Fig. C.1 continued.

Appendix D: Image reconstruction

The incomplete sampling of the Fourier plane forces us to bring in additional information to reconstruct an image. The regularisation plays the role of the missing information by promoting a morphology in the image (i.e. a smooth image). In a Bayesian framework, these terms appear in the a posteriori likelihood where f = χ² + μf_rgl, where χ² is the mean square error, f_rgl is the regularisation function, and μ its weight (see Thiébaut 2008, for more details), which we maximise to find the image. Maximising the a posteriori likelihood of the image amounts to minimising the function f. We used the quadratic smoothness regularisation that promotes a smooth image and is well suited for this object morphology (Renard et al. 2011). We have set this weight to 10¹⁰ for both epochs, applying the L-curve method (see Fig. D.3). It consists of a grid of reconstructions with different weights of the regularisation. The results are plotted on a f_prior versus χ² plot. The resulting curve has an L shape. One asymptote is dominated by the data noise and the other by the regularisation. The location of the bend indicates the optimal regularisation weight.

Fig. D.1 Marginalised probabilities of the chromatic parameters (top) and denv (bottom) for epoch A (left) and C (right). The average value is indicated by solid vertical line and the standard deviation is indicated by dashed lines. The value chosen for the image reconstructions is represent by a vertical solid line, which is red for epoch A and blue for epoch C.

Fig. D.2 Images for different values of the chromatic parameters in the 3σ contours of the Fig. 4. The top line represents the image reconstructions for epoch A and the bottom line the image reconstructions for epoch C. We can see no strong morphological changes for both epochs. For epoch C, however, we can see that the asymmetry in the north-east direction appears with different brightness depending on the chromatic parameters.

Fig. D.3 Determination of the regularisation weight μ by the L-curve method. The reconstructions were performed on a grid of weights (logarithmic grid of 25 values from 107 to 1013) for the two epochs (A in red and C in blue). For each reconstruction, the values of χ2 and fprior were traced in this graph. The best value of the weight is that in the shoulder of the curves. We have chosen 1010 for both epochs.

Fig. D.4 Image reconstruction fit to the dataset on V2 (top) and on CP (bottom) for epoch A (left) and C (right). The black squares represent the data and the blue triangles represent the image reconstruction fit. The residuals normalised to the standard deviation are plotted in the bottom of each graph.

Appendix E: Trajectory of the asymmetry

We have fitted a circular trajectory to the position of the bright spot found in the parametric fit at all epochs but epoch B because the flux of the bright spot is not significant enough (see Fig. E.1). This trajectory can be oriented on the plane of the sky. We allow the centre of the trajectory to be shifted with respect to the central star along the minor axis (Δr). Such a shift can be produced when observing an inclined disk with an asymmetry orbiting the star above the equatorial plane. The best values are presented Table E.1.

Fig. E.1 Best-fit parameters of the circular trajectory. The red diamonds represent the positions of the asymmetry from the best fits. The blue hexagons represent the positions of the asymmetry for the best-fit trajectory. The black star represents the location of the central star. The green ellipse represents the fitted trajectory and the green cross shows its centre.

We obtained a non-satisfactory reduced χ² for the fit but it points out an 88-day period with a radius of 5.2 mas. The orientation is in agreement with the values derived from the model shown in Sect. 4.2 and from previous work (Borges Fernandes et al. 2011). There is a slight difference of 10 degrees on the inclination and the PA. The centre of the trajectory shows a large shift with respect to the star such that the trajectory is passing close to the location of the star.

The image reconstruction epochs A and C are 56 and 63 days long. This is of the order of the period of the bright spot. As a consequence the reconstructed images are probably “blurred”. The simple assumption of a displacement of a bright spot located in the disk is supported by the parametric analysis even if the models we are fitting are too simple to correctly reproduce the entire dataset, especially for the longest baselines. Moreover, the orbit fit is imperfect owing to the poor temporal sampling and its simplicity. It cannot satisfactorily describe the asymmetry evolution.

All Tables

All Figures

	Fig. 1 Temporal representation on data taken on MWC158 during 3 epochs.
In the text

	Fig. 2 { u,v }-plane for the 3 epochs. Epoch A is indicated in red triangles, epoch B in orange circles, and epoch C in blue diamonds.
In the text

	Fig. 3 Temporal representation on data taken on MWC158. The colours correspond to the wavelength. Dataset in V2 and CP. We notice strong differences in the CP between the epochs that may indicate variability.
In the text

Fig. 4 Chromatic parameters determination for epoch A (red) and C (blue) at 1σ (solid lines) and 3σ (dashed lines). The red (epoch A) and blue (epoch C) points are the best values for the chromatic parameters. These values were used for the image reconstruction process. The vertical dashed line indicates the spectral index of the central star, which equals − 4 because the star is assumed to be in the Rayleigh-Jeans regime.

In the text

	Fig. 5 Image reconstructions in H-band on PIONIER data, epochs A (left) and C (right). The colour bar on the right was redefined to have the maximum flux equal to one for the sake of readability. The beam for each epoch is plotted in blue in the bottom right corner of the images. The contours are at 1σ (dashed line), 3, and 5σ (solid lines) significance levels.
In the text

	Fig. 6 Image reconstruction fit on datasets of epoch C superimposed on dataset from epoch A (left) and inversely (right). The black squares represent the data and the blue triangles represent the images. The residuals normalised to the standard deviation are plotted in the bottom of each graph.
In the text

	Fig. 7 Positive part of the asymmetric maps (Asym) of epochs A (left) and C (right). The orientation of the major axis (PA of 71 degrees as estimated in Borges Fernandes et al. 2011) is plotted in black in the image.
In the text

Fig. 8 Location of MWC158 primary star (black hexagon, see text for details) in the isochrones of Bressan et al. (2012) for two scenarios: the scenario in which MWC158 is young, from yellow (0.3 Myr) to red (0.5 Myr), and old, from blue (0.1 Gyr) to cyan (0.2 Gyr). The triangles indicate the position of 6 M⊙ as determined by Borges Fernandes et al. (2009) for the young track. For the old track, these positions are out of range of the plot. The dashed line indicates the upper limit on the absolute H magnitude of a hypothetical secondary of MWC158.

In the text

	Fig. C.1 Histograms of best-fit values for all parameters of the parametric fit for each epoch. The red solid line represents the mean values and the red dashed lines represent the 1σ errors.
In the text

	Fig. C.1 continued.
In the text

	Fig. D.1 Marginalised probabilities of the chromatic parameters (top) and denv (bottom) for epoch A (left) and C (right). The average value is indicated by solid vertical line and the standard deviation is indicated by dashed lines. The value chosen for the image reconstructions is represent by a vertical solid line, which is red for epoch A and blue for epoch C.
In the text

Fig. D.2 Images for different values of the chromatic parameters in the 3σ contours of the Fig. 4. The top line represents the image reconstructions for epoch A and the bottom line the image reconstructions for epoch C. We can see no strong morphological changes for both epochs. For epoch C, however, we can see that the asymmetry in the north-east direction appears with different brightness depending on the chromatic parameters.

In the text

Fig. D.3 Determination of the regularisation weight μ by the L-curve method. The reconstructions were performed on a grid of weights (logarithmic grid of 25 values from 107 to 1013) for the two epochs (A in red and C in blue). For each reconstruction, the values of χ2 and fprior were traced in this graph. The best value of the weight is that in the shoulder of the curves. We have chosen 1010 for both epochs.

In the text

	Fig. D.4 Image reconstruction fit to the dataset on V2 (top) and on CP (bottom) for epoch A (left) and C (right). The black squares represent the data and the blue triangles represent the image reconstruction fit. The residuals normalised to the standard deviation are plotted in the bottom of each graph.
In the text

	Fig. E.1 Best-fit parameters of the circular trajectory. The red diamonds represent the positions of the asymmetry from the best fits. The blue hexagons represent the positions of the asymmetry for the best-fit trajectory. The black star represents the location of the central star. The green ellipse represents the fitted trajectory and the green cross shows its centre.
In the text