Issue 
A&A
Volume 564, April 2014



Article Number  A80  
Number of page(s)  11  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201322926  
Published online  10 April 2014 
SPARCO : a semiparametric approach for image reconstruction of chromatic objects
Application to young stellar objects
^{1}
Institut de Planétologie et d’Astrophysique de Grenoble, UJF,
CNRS
414 rue de la piscine, 38400 Saint Martin
d’
Hères
France
email:
jacques.kluska@obs.ujfgrenoble.fr
^{2}
European Southern Observatory, Alonso de Cordóva 3107, Casilla 19001 Vitacura,
Santiago,
Chile
^{3}
Center for High Angular Resolution Astronomy, Georgia State
University, PO Box
3969, Atlanta,
GA
30302,
USA
^{4}
University of Michigan Astronomy Department,
941 Dennison Bldg,
Ann Arbor, MI
481091090,
USA
^{5}
CRAL, Observatoire de Lyon, CNRS, Univ. Lyon 1, École Normale Supérieure de Lyon,
69364
Lyon,
France
Received:
28
October
2013
Accepted:
11
March
2014
Context. The emergence of optical interferometers with three and more telescopes allows image reconstruction of astronomical objects at the milliarcsecond scale. However, some objects contain components with very different spectral energy distributions (SED; i.e. different temperatures), which produces strong chromatic effects on the interferograms that have to be managed with care by image reconstruction algorithms. For example, the gray approximation for the image reconstruction process results in a degraded image if the total (u,v)coverage given by the spectral supersynthesis is used.
Aims. The relative flux contribution of the central object and an extended structure changes with wavelength for different temperatures. For young stellar objects, the known characteristics of the central object (i.e., stellar SED), or even the fit of the spectral index and the relative flux ratio, can be used to model the central star while reconstructing the image of the extended structure separately.
Methods. We present a new method, called SPARCO (semiparametric algorithm for the image reconstruction of chromatic objects), which describes the spectral characteristics of both the central object and the extended structure to consider them properly when reconstructing the image of the surrounding environment. We adapted two imagereconstruction codes ( Macim , Squeeze , and MiRA ) to implement this new prescription.
Results. SPARCO is applied using Macim , Squeeze , and MiRA on a young stellar object model and also on literature data on HR 5999 in the nearinfrared with the VLTI. We obtain smoother images of the modeled circumstellar emission and improve the χ^{2} by a factor 9.
Conclusions. This method paves the way to improved aperturesynthesis imaging of several young stellar objects with existing datasets. More generally, the approach can be used on astrophysical sources with similar features, such as active galactic nuclei, planetary nebulae, and asymptotic giant branch stars.
Key words: methods: numerical / techniques: high angular resolution / techniques: interferometric / protoplanetary disks / stars: variables: T Tauri, Herbig Ae/Be
© ESO, 2014
1. Introduction
The number of aperturesynthesis images based on optical longbaseline interferometry measurements has recently increased thanks to easier access to visible and infrared interferometers. The interferometry technique has now reached a technical maturity level that opens new avenues for numerous astrophysical topics that require milliarcsecond modelindependent imaging (Berger et al. 2012). Image reconstruction (see review Thiébaut 2013) is the key to achieve the most probable, noncommittal images following some global constraints (image positivity, size of the support, regularization, etc.). A thorough study (Renard et al. 2011) has shown the limitations of image reconstruction, but never challenged the type of regularizations in use. The first images of the inner regions of the environment have been obtained around the young stellar objects HD 163296 (Renard et al. 2010) and HR 5999 (Benisty et al. 2011), revealing the structure around the central objects. However, one main caveat of the images reconstructed from spectrally dispersed instruments is that the visibilities measured at different wavelengths have been assumed to come from a gray object. In addition, the central star has a much higher surface brightness than the surrounding emission, a contrast problem that limits the reliability and the quality of image reconstruction and complicates serious analysis of the morphology of the circumstellar material. These concerns motivate one to reconstruct an image of the envelope alone without the star in the image, an approach considered first in selfaperture masking techniques (Monnier et al. 2004).
Some of the interferometric visibilities obtained with spectral resolution on young stellar objects (YSO), mainly in the nearinfrared (NIR), have been seriously affected by strong spectral dependence. For example, the visibilities of the object MWC 158, as measured with VLTI/AMBER in the K band and in the H band are spread over a broad range (see Fig. 5 of Borges Fernandes et al. 2011). First believed to be quality (...) clearly lower compared to the Kband data (Borges Fernandes et al. 2011), the increase of the visibilities of the Hband data has been observed not only on VLTI/AMBER, but also with VLTI/PIONIER. This indication caused us to consider an astrophysical interpretation: the change of visibility is generated by the different chromatic behavior between the central star that peaks in the visible and the circumstellar material that radiates out mainly in the NIR (Kluska et al. 2012).
The star contribution is highest at the shortest wavelengths and becomes moderate at longer wavelengths. Classic gray algorithms assume the same brightness distribution for all wavelengths. Consequently, they cannot satisfactorily reconstruct an image of a chromatic object with spectral supersynthesis. We need to include the chromatic effect induced by the physical properties of the target in the image processing to retrieve a good intensity map of the observed target.
In this paper, we present a semiparametric approach that includes the knowledge of the relative stellar and environment spectral properties in optical interferometry image reconstruction. As previously demonstrated in parametric modeling of optical interferometric data (e.g. Kraus et al. 2012a), we can directly take the object chromaticity into account in the process of image reconstruction, which improves the final fit.
This approach is called SPARCO (semiparametric approach for image reconstruction of chromatic objects) and consists of separating a wellknown object (e.g. the central star of a YSO) and its complex, unknown environment (e.g. its dusty disk). The star is modeled by a parametric model (that can include hydrostatic models or binaries) and the environment by the reconstructed image. The chromatism is reproduced by changing the flux ratio between the two components across the observed bandwidth.
To present the methods, we focus on YSOs because the star can be modeled at first order in the NIR by an unresolved component and the flux ratio can be represented by a power law with a good approximation. Moreover, the environment is poorly known and is complex. The application of this technique to this type of objects is therefore important.
However, this method does not apply only to YSO, but can also be used in any system where a known source is present that displays a spectral behavior different from the rest of the emitting material in the optical. For instance, the accretion disk of an active galactic nuclei is considered as unresolved in order to retrieve its environment (Kishimoto et al. 2013). A method of separation of the star from its environment was invoked for asymptotic giant branch stars (Hillen et al. 2013), or planetary nebulae (Lagadec et al. 2006).
In Sect. 2 we demonstrate that the SPARCO approach allows proper modeling of the interferometric observables of YSOs, especially their chromatic content. We show how standard image reconstruction algorithms can be modified accordingly. In Sect. 3 we validate the method on the model of a realistic YSO. In Sect. 4, we further discuss important aspects of the method in detail. We finish in Sect. 5 by applying SPARCO on actual data used by Benisty et al. (2011) to reveal the circumstellar environment of HR 5999, and we compare it with a previous analysis.
2. SPARCO method
2.1. Fluxes from the star and the environment
Fig. 1 Left: SED of a young stellar object. In blue: the stellar photosphere (represented as a black body) at 10 000 K; in red: the environment black body at 1500 K; and in black: the sum of the two contributions. The vertical color lines are the spectral channels of PIONIER. The two contributions cross in the H band. In the top right corner as a dashed line: the flux ratios of every component from the SED, as a full line: the result of modeling the fluxes with power laws. Right: in black, expected visibilities from the extended structure alone. In color: the star contribution has been added and colors follow the standard color code (blue: shortest wavelength, red: longest wavelength). The visibilities increase with shorter wavelengths because the stellar contribution is higher. 
A YSO, consisting of a star and a dusty environment, has two main components in the spectral energy distribution^{1}. The photosphere dominates from the UV domain to the visible. The contribution of the environment, which is mainly dusty and occurs at lower temperature (T< 1500 K), prevails from the infrared to radio wavelengths (see Fig. 1 left).
For example, the effective temperature of Herbig Ae/Be stars is T ≈ 10 000 K. The spectral domain where the total flux is no longer dominated by the stellar flux is in the NIR where the emission of the environment increases steeply. Typically, the stellar SED is in the RayleighJeans regime in the H band (1.65 μm), while the environment is in the Wien regime. In these spectral bands (J, H and K bands for young stellar objects), the shortest wavelengths are dominated by the photosphere flux, the longest ones by the environment flux. These two components (the star and the dusty environment) have different spatial extents, which can be resolved by current optical interferometers.
Some interferometric instruments (AMBER, PIONIER on the VLTI, and MIRC on the CHARA interferometer) cover the H band with several spectral channels. The change of flux ratio (star/environement) for each channel following the physical laws described above implies a strong chromatic effect in the visibilities called star/environment chromatism (see Fig. 1 left).
When considering low spectral resolution within a single NIR band, this effect can be described by the following two parameters:

: the stellartototalflux ratio at wavelength λ_{0}. This parameter sets the flux balance between the two components at the reference wavelength λ_{0}. λ_{0} is arbitrary chosen (for example, the central wavelength of observations).

d_{env}: the spectral index for the circumstellar environment. For interferometric data, only the difference between the spectral indexes of the two components matters to account for the chromatism. In the NIR, the star emission occurs in the RayleighJeans regime and can be approximated by . The parameter that sets the index difference is the absolute spectral index of the dust, .
The total flux f_{tot} at a wavelength λ normalized by the total flux at λ_{0} can therefore be written as (1)This model is sufficient to describe the continuum emission of the object in a NIR band (especially the H band). With additional effort, complex flux distribution (e.g. adding lines) can be used to model observations at higher spectral resolution, following Eq. (2).
2.2. Visibilities from the star and the environment
The total complex visibilities are the sum of the stellar visibilities and visibilities from the environmental components weighted by their fluxes at the corresponding wavelengths (λ), (2)with f_{tot}(λ) given by Eq. (1), f_{∗}(λ) the stellar flux, f_{env}(λ) the flux of the environment, the stellar visibility, the visibility of the environment. Finally, b is the interferometric baseline vector projected onsky. The visibility depends on b/λ, which is the spatial frequency. In this paper all the quantities with a tilde (as ) are complex numbers.
In YSO, the apparent size of the central star can be considered as unresolved. For instance, a Herbig AeBe star at a typical distance of 140 pc (distance of the Taurus starforming region) and with a radius of 5 R_{⊙} has an angular radius of ≈0.17 mas. For a 100 m baseline and for the H band, this corresponds to a visibility V ≈ 0.997. Hereafter we model this star with an unresolved component with V_{∗} ≈ 1. By developing the total flux term, we have (3)where is the visibility of the environment alone, derived from its brightness distribution by a Fourier transform. With the two parameters described in Sect. 2.1 we can rewrite Eq. (3) as follows: (4)
2.3. Image reconstruction
is retrieved by making the Fourier transform of its image. This image is obtained by a wellknown imagereconstruction process. The goal is to retrieve the most probable image given the dataset and some assumption called regularizers (e.g. image positivity). To reconstruct the image we have to solve an illposed inverse problem by minimizing the function defined as (Thiebaut & Giovannelli 2010) (5) being the global distance to minimize, the distance to the data (reduced χ^{2}), the regularization distance, x the image pixel values, and μ the regularization weight. The choice of the hyperparameter μ is discussed in Appendix B. For more information about the regularization see Renard et al. (2011).
Various image reconstruction algorithms exist (e.g. MiRA Thiébaut 2008; Macim Ireland et al. 2006; and Squeeze Baron et al. 2010). They mainly differ by the way they minimize the function (gradient descent or Monte Carlo Markov chain (MCMC) minimization).
The imagereconstruction process included in SPARCO is monochromatic. In other words, the object intensity distribution is wavelength independent. This image changes its flux ratio across the observed band only thanks to Eq. (3).
Because most algorithms compute the complex visibilities of the image, our method can be implemented in any of them. They have to be modified to include Eq. (4) in their iterative computation of visibilities ( being the complex visibilities of the image). If the imagereconstruction algorithm is based on a gradient descend, the gradient needs to be multiplied by the factor (which is the environmenttototalflux ratio). Current MCMC reconstruction algorithms (i.e. Macim and Squeeze ) use stochastic steps whose modifications to the χ^{2} are computed by finite difference, and thus do not require any gradient evaluation.
The proposed method couples (1) the fitting of a parametric model for one part of the object with (2) a simultaneous image reconstruction of the remaining part. The details of the algorithm are described in Appendix C.
3. Numerical validation
To validate our method and estimate its capability to retrieve object features, we built a synthetic model. This geometrical model of an unresolved star and its surrounding environment includes the chromatic effect: the star is hotter than its environment. This object was used to simulate interferometric observations in a realistic configuration.
3.1. Model description
Fig. 2 Top: image of the ring with the star represented in red (model), the (u,v)plan used in the simulation, and the simulated dataset. Bottom: image reconstructions with the four cases discussed in Sect. 3.2 : 1) classical gray image reconstruction; 2) and d_{env} = d_{star}; 3) and d_{star} = − 4 and d_{env} = 1 with MiRA/SPARCO; 4) idem with MACIM/SPARCO. 
We used an analytical model consisting of an unresolved star and its environment. The star was assumed to be in the RayleighJeans regime so that its spectral dependence is: . Its spectral index is then − 4. The star contributes 40% of the flux at 1.65 μm ().
The environment was modeled by a Gaussian ring with a radius of 6 mas inclined by 60 degrees. An azimuthal modulation set the ring flux to vary as the cosine of the azimuthal angle. This reproduces the asymmetry generated by radiative transfer effects (Monnier et al. 2006). The Gaussian thickness of the ring (the FWHM of the Gaussian with which we convolved the infinitesimal ring) is 2.4 mas. The center of the ring is shifted by 3 mas to the south with respect to the star. This shift can reproduce a perturbation induced by a companion or an inclination effect. The environment spectral index is d_{env} = 1, which is the logarithmic black bodycurve derivative at λ_{0} and at a temperature of ≈1400 K. At 1.5 μm the star carries ~55% of the flux and ~38% at 1.8 μm. This chromatism is strong enough to significantly affect the data. All the parameters of this model are summarized in Table 1.
From this model we simulated a realistic dataset as it would have been obtained by the PIONIER instrument (Le Bouquin et al. 2011). We used a actual (u,v)plan made with PIONIER (see Fig. 2) consisting of 14 pointings of four telescopes on three different configurations available at the VLTI (two pointings on small configuration, nine on the medium one, and three on the large one). One pointing represents 40 min of observation (calibrators included). Each point on each baseline is spectrally dispersed onto seven channels across the H band. PIONIER only provides V^{2} and closure phases (CP; no complex visibilities, no complex differential visibilities). We added realistic noise to the data by selecting three regimes: high, intermediate, and low flux (see Table 2).
The resulting dataset is shown in Fig. 2, top. This artificial dataset is qualitatively similar to real observations obtained with existing optical interferometers (e.g. Kluska et al. 2012). The overall circumstellar structure is resolved by the longest baselines, and the V^{2} data display strong chromaticity effects.
Model parameters.
Model noise estimated from Le Bouquin et al. (2011).
3.2. Validation of the chromatic imagereconstruction method
In this subsection we assume a perfect knowledge of the object’s chromatic parameters. We compare the images obtained with the classical “gray” approach to the SPARCO approach and check the validity in different signaltonoise regimes.
Except when explicitly stated, the images were reconstructed with the MiRA algorithm which computes the image following a gradient method. The regularization used was “total variation”, which was described in Renard et al. (2011) as the most successful regularization for astrophysical objects. The process of selecting the optical regularization weight (μ, see Eq. (5)) is described in Appendix B. In all cases, we performed the reconstruction with μ = 1500. We emphasize that the regularization and its hyperparameter has to be adapted to every reconstructed object. One has to be careful when choosing the value of the hyperparameter μ although we show (in Appendix B) the low sensitivity of the method to the choice of regularization type. Appendix B also presents several methods that help to chose the regularization. We reconstructed images (see Fig. 2) with 256 × 256 pixels of 0.2 mas each with three different methods:

1.
A classic, gray image reconstruction, that is, with (in which case the d_{env} parameter has no meaning).

2.
An image reconstruction considering a central point source (), but without taking into account the difference of spectral indicies (d_{env} = d_{star} = −4).

3.
An image reconstruction with the full SPARCO approach, that is, considering the difference of spectral index between the star and the environment (). To show the effect of the SPARCO approach we implemented it in three different imagereconstruction algorithms ( MACIM , Squeeze , and MiRA ).
The reconstructed images are shown in Fig. 2 bottom. In the classical gray case the algorithm poorly fits the data (χ^{2} = 9.3). It is possible to lower the χ^{2} by significantly increasing the field of view of the image and the number of pixels. With these additional degrees of freedom, the algorithm is able to reproduce the chromatic effect by adding ripples in the Fourier space. This creates strong artifacts at both small and large separations in the image. Somehow it replaces chromatic effects by spatial artifacts linked with the Fourier sampling. For instance, the ring PA is related to the (u,v)plan orientation. These effects are still observed when subtracting the star monochromatically.
The full SPARCO approach allows one to reach a good χ^{2} (≈1.2). The reconstructed image does not show strong artifacts. The PA of the ring, its offset, and its azimuthal modulation are correctly reproduced. We conclude that subtracting the star and taking into account the difference of spectral index differences is mandatory to reconstruct a reliable image of the environment.
Additionally, we verified that the conclusions remain unchanged when using other imagereconstruction algorithms than MiRA . We performed the same exercise using Macim and Squeeze with a regularization based on the Laplacian on the image. These algorithms rely on MCMC methods that were modified to handle the SPARCO approach. We obtained the same χ^{2}(1.2) and the same conclusions (see Fig. 2 bottomright panel).
3.3. Determination of chromatic properties
Fig. 3 MiRA image reconstructions as a function of the assumed chromatic parameters: the stellartototalflux ratio and the dust spectral index d_{env}. The star is represented in red at the center of each image. The true values are and d_{env} = 1. Top right: χ^{2} map of the reconstructions as a function of the chromatic parameters for reconstructions with the total variation regularization. All χ^{2}sup3 are represented in black. The blue crosses represent the location of the images in the χ^{2} map. We can clearly see that this map is degenerated. 
The previous section presented a validation of the method in the ideal case where the chromatic parameters ( and d_{env}) are known. In this section, we explore the effect of varying and d_{env} (in other words: the values used for the image reconstruction do not correspond to those used to simulate the dataset). This allows us to conclude whether it is possible to recover these parameters from the dataset.
The range of values for the environment spectral index d_{env} was chosen to cover the temperatures of dust sublimation (from 2200 K to 1100 K, that is, d_{env} = − 1 to 3 respectively). The range of values of the flux ratio was 0.2 to 0.6. The results are presented Fig. 3. We found that the image morphology weakly depends on d_{env} in the considered range of values. On the other hand, significantly influences the morphology of the reconstructed image. If the stellar flux ratio is too high, the algorithm will make an image with a large inner hole that affects the ring. In the opposite case, if the startototal flux ratio () is too low, the algorithm will add flux at the star position. But the added flux is at the dust temperature. To compensate for this, the algorithm has to create artifacts in the image of the environment.
The topright inset of Fig. 3 displays the χ^{2} value obtained at the end of imagereconstruction processes for the grid of and d_{env}. The χ^{2} is an indicator of the distance of the model (parameterized by the image pixel values, and d_{env}) to the dataset. The map shows a single global minimum. This validates the capability of reconstructing the image of the environment while simultaneously fitting the chromatic parameters ( and d_{env}). Appendix D demonstrates that this joint minimization is robust: it does not depend on the starting point or on the choice of the regularization.
However, the χ^{2}map shows a correlation between the parameters and d_{env}. We verified that this degeneracy is not caused by regularization or the S/R ratio but is intrinsically linked with the interferometric data. As pointed out in Sect. 3.2, artifacts in the reconstructed image compensate for an error on the chromatic parameters. This degeneracy seems to be carried out by these artifacts. Solutions to that degeneracy problem are discussed in the next sections.
4. Discussion
4.1. Need for spectrophotometry
We showed in Sect. 3.3 that one can fit and d_{env} when reconstructing the image, but that these parameters are strongly correlated. Nevertheless, one can at least derive an upper limit to the star flux ratio () and a lower limit to the environment spectral index (d_{env}).
Improved image reconstruction is possible using additional constraints on and/or d_{env} from other observations, for instance, spectrophotometry. This is possible when the flux emitted by the central object can be safely extrapolated in the NIR, for instance when the central object is a star whose spectral type and extinction are known. In this case, spectrophotometric observations in the NIR allow one to recover and/or d_{env}. However, this simplistic interpretation of the SED should be applied with caution, as illustrated in Sect. 5.
4.2. Effect of the gradient temperature in YSO’s disks.
Disks around young stellar objects display a temperature gradient in the radial direction. It has the same signature in the visibility as the star/environment chromatism (e.g., Kraus et al. 2012b). So far, this effect is not taken into account in SPARCO: we assume that the spectral index of the environment is not changing with radius.
Fig. 4 Reconstructions with an intensity gradient in the disk (left) and a temperature gradient in the disk ∝ (right). The reconstructions have a χ^{2} ≈ 1. Bottom: cut of the reconstructed images (black) and of the disk models. The colors are for the different wavelengths in the right panel (blue: 1.55 μm, red: 1.8 μm). The flux level is given on a logarithmic scale. 
To test this, we built a disk model with a temperature gradient corresponding to a flared accretion disk with (Kenyon & Hartmann 1987). We also built a reference model with an intensity gradient matching the average intensity distribution of the first model. The flux distribution is 50% for the star and 50% for the disk in the middle of the H band. The disks extend from 5 mas to 10 mas. The temperature of the inner rim is 1500 K.
Figure 4 shows the reconstructed images with the SPARCO approach setting and letting d_{env} free to vary. They have acceptable χ^{2} values (≈1.0). The recovered d_{env} corresponds to 1400 K for the temperature gradient case and 1500 K for the intensity gradient case. Interestingly, 1400 K corresponds to the average temperature of the accretion disk. This is therefore a satisfactory value. The overall shape (size and radial extension) of the disk is qualitatively recovered in both cases. The inner part of the disk dominates the emission. The differences between the two models are significantly smaller than the dynamic range of the image reconstructions.
We also made the test with a flat accretion disk () and arrived at the same conclusions. We conclude that the effect of a temperature gradient is negligible for image reconstruction of YSOs over the H or the K band.
However, for multiband datasets we recommend splitting the data into several spectral bands, which better adapts the flux description. In the mid or farinfrared, where the star is negligible, we would need to modify the method for image reconstruction.
4.3. Evolution of the parametric model
The parametric model of the star as an unresolved object is sufficient to detect interesting features in an object, but it may be inefficient for imaging objects that are more complex for observations with an interferometer with a sufficient resolution to resolve the central star. We can modify the algorithm to subtract other shapes than just an unresolved star. In Eq. (2) we can replace by any model with an analytical formula or even an image. For a binary we can upgrade this method to subtract both of the components and even find their positions.
In addition, the adaptation to other spectral bands for YSOs or to instruments with better spectral resolutions will pass through the modification of the flux laws for the star and the environment. We chose power laws because in the H and K bands (which are quite narrow) they reproduce the spectral behaviors of the two components well. Another parametric flux description can replace the current one in the Eq. (3).
5. Application on actual datasets: HR 5999
In this section we apply the SPARCO method to the actual dataset from Benisty et al. (2011) on HR 5999, a Herbig Ae star (Tjin A Djie et al. 1989). This dataset was obtained by the VLTI/AMBER instrument, which is a threebeam combiner in the NIR. The object shows an excess in near and midinfrared spectral energy distribution. It was first imaged by Benisty et al. (2011) (see Fig. 6). The photosphere contributes 22% of the flux in K band. Moreover, there is an inner disk that is marginally resolved (between 0.43 mas and 2.8 mas) and contributes 38% of the flux. It might be interpreted as a gaseous disk or as refractory dust grains. The remaining flux comes from the inner dusty rim, which is believed to be at the dustsublimation radius.
Fig. 5 SED of HR 5999 from Benisty et al. (2011). We can see that in the Kband (at log(2.2 μm) = 0.35), the emission is dominated by the environment which has two components. The star is weaker but has to be taken into account for the image reconstruction. 
This object is quite complex because there is a third component that contributes to the NIR flux: the inner disk. Our method is not yet designed to answer this complexity. Nevertheless, we adopted a strategy in three steps to apply our method to this complex case:

Reconstruction by subtracting the photosphere only from theimage: using the SED (Fig. 5), we set the stellar tototal flux ratio to the value of and the temperature of theenvironment to 1500 K. The image (Fig. 6 center) showsan unresolved component in the center and two patterns at bothsides of the unresolved component. The subtracted flux is locatedat the position of the red star. As expected, the algorithm sees theunresolved flux, which is the inner disk, but has difficulties toreconstruct the environment correctly. We note that the positionangle is similar to that found in Benistyet al. (2011).

Reconstruction by subtracting the photosphere and the inner disk: the flux in the unresolved component is now set to 60% (as show in the SED). The environment is still set to be at 1500 K (d_{env} ≈ 1). There is no feature at the center of the image (see Fig. 6). We can see the inner dusty disk rim. The image looks cleaned in the center because the fluxes (star + inner disk) that we are subtracting are not necessarily unresolved (i.e., the inner disk was marginally resolved).

Fit of the chromatic parameters in the reconstruction process: if we did not have the information on the photometry on these objects, we would try the extended method described in Sect. 3.3. We converge to subtracted flux of 40% and a relative spectral index d_{env} of −1 (which translates into a temperature of 2100 K if we assume that we only subtract the photosphere). We subtract the photosphere and half of the inner disk. The rest of the inner disk is considered as resolved.
The first conclusion of the application of SPARCO on this dataset is that we retrieve the inner dust rim if we subtract the photosphere and the inner disk. We have to be careful to correct the spectral index of the environment d_{env}, which depends on the temperature of the unresolved flux. If we only subtract the photosphere, SPARCO reconstructs the unresolved flux from the inner gaseous disk.
Fig. 6 Image reconstructions on HR 5999. Left: Benisty et al. (2011). Center: SPARCO with 22% of flux in the central red star. This corresponds to the amount of flux from the photosphere. Right: SPARCO with 60% of flux (corresponds to the photosphere and inner disk contribution) in the central star. 
Second, if we let the method fit the chromatic parameters in the reconstruction process, it finds an that is larger than the photosphere contribution and a spectral index d_{env} that indicates a temperature of more than 2000 K. This is higher than the dust sublimation temperature of 1500 K. This higher temperature is explained by the fact that the unresolved component around the star has a spectral index in K band “cooler” than the star itself, so the difference in spectral indices favored from the chromatic effect in the interferometric data is lower. This will automatically increase the derived temperature of the extended component if we still assume that the unresolved component is in the RayleighJeans regime.
However, Benisty et al. (2011) suggested that the origin of the emission might be refractory dust grains or a gaseous disk. This indicates that the hypothesis of the photosphere alone as the unresolved component can be excluded.
The method is sensitive to an unresolved flux that does not come from the photosphere alone. The SED clearly indicates an inner component, as suggested in Benisty et al. (2011).
6. Conclusion
For image reconstruction based on interferometric data with spectral dispersion, visibilities are determined not only by the geometry of the object, but also by the (differential) spectral slope of its components.We developed a method that includes the knowledge that we have on the relative spectral behaviors of two components to reconstruct the intensity distribution of the extended one. This method allowed us to improve the χ^{2} by one order of magnitude in our validation.
The SPARCO method includes an analytical description of the stellar contribution and the chromatic ratio between the star and its environment in the imagereconstruction algorithms used in optical interferometry. The first component is modeled parametrically and an image is reconstructed for the second one.
SPARCO was used correctly on a young stellar object with the following hypothesis:

The environment spatial distribution is assumed to be wavelength independent. The chromatic dependence of visibilities only arises because the star has a different spectral dependence from the environment. This also implies that the environment has the same spectral dependence in the whole image.

The star is very close to the parametric description given in the method.

The fluxes (F_{λ}) of the components are approximated by power laws as a function of the wavelength. Since the chromatic parameter space can be degenerate when considering interferometric constraints alone, obtaining independent spectrophotometric observations helps to retrieve an image of the observed environment. Despite the low sensitivity of the method to the choice of regularization, one has to be careful when choosing the value of the hyperparameter μ .
The method is simple but can be adapted to more complex models. As demonstrated in Sect. 5, if the object is more complex some more information needs to be added to the model (add a parametric disk for the inner part of the object). In Sect. 4.3 we showed that we can easily model the parametric part by a uniform disk, a binary, and even an image.
This tool is complementary to other observations of the target, especially with spectrophotometric observations. Since the chromatic parameter space can be degenerate, these types of observations are very important to correctly retrieve an image of the observed target. It is difficult to retrieve them by simultaneous fitting.
The SPARCO method will be intensively used on a Large Program dataset on Herbig Ae/Be stars gathered by PIONIER at the VLTI. Unraveling the image of the close environment of young stars will help us to constrain the effects of the inclination of the inner regions of YSO and therefore to detect early signs of planet formation very close to the star. This method cannot be applied only to young stars, but can also be used in any system whith a pointlike source and that displays a different spectral behavior from the rest of the emitting material in the optical such as active galactic nuclei.
Algorithms that implement a fully polychromatic approach are currently being developped, e.g. MiRA (Thiébaut et al. 2013) and Squeeze (Baron et al., in prep.). Combining the SPARCO approach with polychromatic reconstruction will allow imaging any stellar environment with limited perturbation from its star, and greatly enhance our capability to study variations in the environment morphology (e.g. temperature gradient in the disk).
The presence of an additional hot inner component as suggested by Tannirkulam et al. (2008) and Benisty et al. (2011) is discussed Sect. 5.
Acknowledgments
This work is supported by the French ANR POLCA project (Processing of pOLychromatic interferometriC data for Astrophysics, ANR10BLAN0511). We acknowledge Myriam Benisty for fruitful discussions and for providing us the dataset on HR 5999. We acknowledge Christophe Pinte, Michel Tallon and Isabelle TallonBosc for interesting discussions that led to this paper.
References
 Baron, F., Monnier, J. D., & Kloppenborg, B. 2010, in SPIE Conf. Ser., 7734, 478 [Google Scholar]
 Benisty, M., Renard, S., Natta, A., et al. 2011, A&A, 531, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berger, J.P., Malbet, F., Baron, F., et al. 2012, A&ARv, 20, 53 [NASA ADS] [CrossRef] [Google Scholar]
 BorgesFernandes, M., Meilland, A., Bendjoya, P., et al. 2011, A&A, 528, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillen, M., Verhoelst, T., Van Winckel, H., et al. 2013, A&A, 559, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ireland, M. J., Monnier, J. D., & Thureau, N. 2006, SPIE Conf. Ser., 6268, 458 [NASA ADS] [Google Scholar]
 Kenyon, S. J. & Hartmann, L. 1987, ApJ, 323, 714 [NASA ADS] [CrossRef] [Google Scholar]
 Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2013, ApJ, 775, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Kluska, J., Malbet, F., Berger, J.P., et al. 2012, in SPIE Conf. Ser., 8445 [Google Scholar]
 Kraus, S., Calvet, N., Hartmann, L., et al. 2012a, ApJ, 752, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Kraus, S., Monnier, J. D., Che, X., et al. 2012b, ApJ, 744, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Lagadec, E., Chesneau, O., Matsuura, M., et al. 2006, A&A, 448, 203 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LeBouquin, J.B., Berger, J.P., Lazareff, B., et al. 2011, A&A, 535, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Monnier, J. D., MillanGabet, R., Tuthill, P. G., et al. 2004, ApJ, 605, 436 [NASA ADS] [CrossRef] [Google Scholar]
 Monnier, J. D., Berger, J.P., MillanGabet, R., et al. 2006, ApJ, 647, 444 [NASA ADS] [CrossRef] [Google Scholar]
 Renard, S., Malbet, F., Benisty, M., Thiébaut, E., & Berger, J.P. 2010, A&A, 519, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Renard, S., Thiébaut, E., & Malbet, F. 2011, A&A, 533, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tannirkulam, A., Monnier, J. D., MillanGabet, R., et al. 2008, ApJ, 677, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Thiébaut, E. 2008, in SPIE Conf. Ser., 7013 [Google Scholar]
 Thiébaut, É. 2013, in EAS Pub. Ser. 59, eds. D. Mary, C. Theys, & C. Aime, 157 [Google Scholar]
 Thiebaut, E., & Giovannelli, J.F. 2010, IEEE Signal Processing Magazine, 27, 97 [CrossRef] [Google Scholar]
 Thiébaut, É., Soulez, F., & Denis, L. 2013, J. Opt. Soc. Am. A, 30, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Tjin A Djie, H. R. E., The, P. S., Andersen, J., et al. 1989, A&AS, 78, 1 [NASA ADS] [Google Scholar]
Appendix A: Reconstructions with different signaltonoise levels
Fig. A.1 Intermediate and high photon noise regimes, left: squared visibilities, right: closure phases. 
Fig. A.2 MiRA image reconstructions with SPARCO . The topleft image is the original model image. In the top right we show the reconstruction with the lowphoton noise level. In the bottom left corner there is an image reconstruction with the average photon noise level. In the bottomright corner the reconstruction was made with the high photon noise level. The red star represents the central star in all the images. With this method, the star is not represented in the image as flux in the pixels. In the bottomright corner of each image the dirtybeam is represented. The dirtybeam is the equivalent of the point spread function in interferometry. It is computed from the (u,v)plan. 
We show in Fig. A.1 the simulated datasets with intermediate and high S/N.
We compare the reconstructions with different signal to noise ratios (S/N). The starting point of the reconstructions is an image with flux in one central pixel alone. The results with the three different noise levels are presented Fig. A.2. We used the chromatic parameters corresponding to those we set in our model.
The ring is well reproduced in the sense that its size, width, and orientation are good. Furthermore, as expected, the algorithm is able to shift the environment with respect to the star.
Moreover, the azimuthal modulations are well reproduced even if there is a quality variation with the noise. It does not interpret the closure phases by asymmetric shapes. For the noise contribution, the global shape of the object and the details are well represented even with the high photon noise level. We conclude that there is no strong limitation on retrieving the global shapes and asymmetries of the objects in the PIONIER operating regime (see Fig. A.2).
Appendix B: Choice of the regularization
Appendix B.1: Choice of the hyperparameter μ
We reconstructed the images using the total variation regularization, because this regularization has been shown to be welladapted to many types of objects (Renard et al. 2011). As mentioned in Sect. 2.3, we need to tune the hyperparameter μ in Eq. (5) to set the weight of the regularization. This needs to be done carefully to obtain a strong regularization and to still stick to the data. We present in Fig. B.1 three graphs that helped us to find a good value of the hyperparameter μ. Ultimately, the best way to find the best value for μ (which is noted μ^{+}) is to compare the reconstructed image with the real one (when the model is available). The images were normalized to a total flux of unity. We summed the differences of pixel values throughout the image. This is called the quality criterium (QC): (B.1)where is the reconstructed image pixel with the regularization weight μ, andx^{ref} are the pixels of the model image. The disadvantage of this method is that we must have the “true” image, which is not possible with real data. That is why this QC method indicates whether the definition of μ^{+} found by another method is accurate. We made the QC analysis for our model with reconstruction using SPARCO/MiRA with totalvariation regularization (see Fig. B.1 right). Each reconstruction starts with all the flux in one central pixel. We found μ^{+} ≈ 1500. The right graph in Fig. B.1 shows a wide range of μ that corresponds to a good image reconstruction (typically 500 <μ< 2000).
For real datasets, we have no access to the real image. We aim for an image that is well regularized but still sticks to the data. The χ^{2}) curve (Fig. B.1 on the left) increases sharply with μ>μ^{+}.
This is confirmed by the L curve μ^{+} determination (see Fig. B.1 center). We plot vs. in a log log graph. The μ^{+} location should be in the corner of the curve. It actually is just in front of a strong deviation from the minimum of .
Two regimes are present for μ>μ^{+}. We can see them on all the graphs of Fig. B.1. This is the effect of the gradientdescent method that remains in local minima for some μ values and moves to a better minima for others. To correct that we could start with a prior image that is closer to what we expect to see. This effect should not be noticeable with Macim and Squeeze since the Monte Carlo method is less likely to remain in local minima.
Fig. B.1 Tuning the hyperparameter μ for the totalvariation regularization with MiRA/SPARCO . Left: cost functions vs. the hyperparameter μ. Red line: χ^{2} (=). Black dashed line: (multiplied by 1000 for better graph visibility). Black solid line: . Blue line: . The vertical dotted line corresponds to the best regularization weight μ^{+} found by the quality criterium. Center: The L curve is one of the criteria to choose the regularization weight μ. Black line: the values of vs. for a range of μ values from 1 to 10^{6}. Blue triangles: points where μ = {1,10,10^{2},10^{3},10^{4},10^{5},10^{6} }. Red triangle: the value where μ = μ^{+}. Right: the quality criterium graph is the distance between the pixels of the reconstructed image to the real (here the model) one. Blue line: the value of the quality criterium vs. the hyperparameter μ. Dotted line: the position of the minimum of the curve and the corresponding value of μ = μ^{+}. 
In all the following reconstructions we took μ to 1500. To allow comparisons we did not change its value between the different methods.
Appendix B.2: Effect of regularization on chromaticparameter χ^{2}map
We can compare the χ^{2} map of Fig. B.2 realized with the smoothness regularization with Fig. 3, which was realized with the totalvariation regularization. These two maps are very similar. The degeneracy is present for both regularizations. Even when we change the hyperparameter μ, the area where the χ^{2} values are acceptable is still large.
Fig. B.2 Choice of the μ parameter and its influence on the fit of the chromatic parameters. χ^{2} map of the reconstructions as a function of the chromatic parameters for reconstructions with the smoothness regularization. Black represents χ^{2}> 3. The blue crosses represent the location of the images of Fig. 3 in the χ^{2} map. The four figures represent μ = { 10^{7},10^{8},10^{9},10^{10} } respectively. The χ^{2} does not change significantly with the μ hyperparameter, but the degeneracy is always strong. 
Appendix C: Method of chromatic parameter fitting
We propose the following algorithm iterative sequence (where i is the iteration index) for a gradientdescentbased algorithm for image reconstruction (e.g. MiRA ):

1.
Computing the total visibilities (): to compute the total complexvisibilities, we use Eq. (4) using (which is the Fourier transform of x^{i}) and the chromatic parameters and .

2.
Computing the χ^{2}gradient: if the modified algorithm is based on gradient minimization, it computes the flux gradient on every pixel. The first step is to compute the errors and the gradients on the total visibilities. Then the environment visibility gradient is obtained by multiplying the total visibility gradient by a factor (which is the environmenttototalflux ratio). The last step is a Fourier transform of these gradients to obtain the flux gradients on the pixels.

3.
Making a step in the pixel fluxes: the algorithm will add to or subtract some flux from on every pixel depending on the gradient and the boundaries. We then have x^{i + 1}.

4.
Computing the visibilities from the image: this step is made with and . We then have .

5.
Fitting the chromatic parameters: the chromatic parameters are fit between two steps of the image reconstruction given the current image, in this way we obtain . This step can be made using a LevenbergMarquardt minimization. Then we return to the first step.
where at the ith iteration we have an image x^{i} and chromatic parameters and .
Appendix D: Convergence of the chromatic parameter fitting
To test the results of a fit of the parameters without an SED constrain we produced a series of image reconstructions with fitting processes using a hundred of different starting points that formed a regular grid on the chromatic parameters (f_{*}^{0} from 10% to 90% and d_{env} from −4 to 4). The results are given in Fig D.1. The results of all the fitting processes end in the χ^{2} valley, but not always close to the good values. This underlines that it is important to know the spectrophotometric information of the observed target and that the degeneracy is deeply embedded in the interferometric data.
To conclude, the fitting process finds the global minimum valley. The starting points do not influence the convergence (see Fig. D.1).
Fig. D.1 χ^{2} map with the fit departures and arrivals for the model. The color represent the χ^{2} for every pair of chromatic parameters and d_{env}. The white dots are the starting pair for image reconstructions. The red points represent the found parameters for the final image. The fit always ends in the χ^{2} minimum but far from the correct values ( and d_{env} = 1). The black solid lines link the start and the final values of one reconstruction. 
All Tables
All Figures
Fig. 1 Left: SED of a young stellar object. In blue: the stellar photosphere (represented as a black body) at 10 000 K; in red: the environment black body at 1500 K; and in black: the sum of the two contributions. The vertical color lines are the spectral channels of PIONIER. The two contributions cross in the H band. In the top right corner as a dashed line: the flux ratios of every component from the SED, as a full line: the result of modeling the fluxes with power laws. Right: in black, expected visibilities from the extended structure alone. In color: the star contribution has been added and colors follow the standard color code (blue: shortest wavelength, red: longest wavelength). The visibilities increase with shorter wavelengths because the stellar contribution is higher. 

In the text 
Fig. 2 Top: image of the ring with the star represented in red (model), the (u,v)plan used in the simulation, and the simulated dataset. Bottom: image reconstructions with the four cases discussed in Sect. 3.2 : 1) classical gray image reconstruction; 2) and d_{env} = d_{star}; 3) and d_{star} = − 4 and d_{env} = 1 with MiRA/SPARCO; 4) idem with MACIM/SPARCO. 

In the text 
Fig. 3 MiRA image reconstructions as a function of the assumed chromatic parameters: the stellartototalflux ratio and the dust spectral index d_{env}. The star is represented in red at the center of each image. The true values are and d_{env} = 1. Top right: χ^{2} map of the reconstructions as a function of the chromatic parameters for reconstructions with the total variation regularization. All χ^{2}sup3 are represented in black. The blue crosses represent the location of the images in the χ^{2} map. We can clearly see that this map is degenerated. 

In the text 
Fig. 4 Reconstructions with an intensity gradient in the disk (left) and a temperature gradient in the disk ∝ (right). The reconstructions have a χ^{2} ≈ 1. Bottom: cut of the reconstructed images (black) and of the disk models. The colors are for the different wavelengths in the right panel (blue: 1.55 μm, red: 1.8 μm). The flux level is given on a logarithmic scale. 

In the text 
Fig. 5 SED of HR 5999 from Benisty et al. (2011). We can see that in the Kband (at log(2.2 μm) = 0.35), the emission is dominated by the environment which has two components. The star is weaker but has to be taken into account for the image reconstruction. 

In the text 
Fig. 6 Image reconstructions on HR 5999. Left: Benisty et al. (2011). Center: SPARCO with 22% of flux in the central red star. This corresponds to the amount of flux from the photosphere. Right: SPARCO with 60% of flux (corresponds to the photosphere and inner disk contribution) in the central star. 

In the text 
Fig. A.1 Intermediate and high photon noise regimes, left: squared visibilities, right: closure phases. 

In the text 
Fig. A.2 MiRA image reconstructions with SPARCO . The topleft image is the original model image. In the top right we show the reconstruction with the lowphoton noise level. In the bottom left corner there is an image reconstruction with the average photon noise level. In the bottomright corner the reconstruction was made with the high photon noise level. The red star represents the central star in all the images. With this method, the star is not represented in the image as flux in the pixels. In the bottomright corner of each image the dirtybeam is represented. The dirtybeam is the equivalent of the point spread function in interferometry. It is computed from the (u,v)plan. 

In the text 
Fig. B.1 Tuning the hyperparameter μ for the totalvariation regularization with MiRA/SPARCO . Left: cost functions vs. the hyperparameter μ. Red line: χ^{2} (=). Black dashed line: (multiplied by 1000 for better graph visibility). Black solid line: . Blue line: . The vertical dotted line corresponds to the best regularization weight μ^{+} found by the quality criterium. Center: The L curve is one of the criteria to choose the regularization weight μ. Black line: the values of vs. for a range of μ values from 1 to 10^{6}. Blue triangles: points where μ = {1,10,10^{2},10^{3},10^{4},10^{5},10^{6} }. Red triangle: the value where μ = μ^{+}. Right: the quality criterium graph is the distance between the pixels of the reconstructed image to the real (here the model) one. Blue line: the value of the quality criterium vs. the hyperparameter μ. Dotted line: the position of the minimum of the curve and the corresponding value of μ = μ^{+}. 

In the text 
Fig. B.2 Choice of the μ parameter and its influence on the fit of the chromatic parameters. χ^{2} map of the reconstructions as a function of the chromatic parameters for reconstructions with the smoothness regularization. Black represents χ^{2}> 3. The blue crosses represent the location of the images of Fig. 3 in the χ^{2} map. The four figures represent μ = { 10^{7},10^{8},10^{9},10^{10} } respectively. The χ^{2} does not change significantly with the μ hyperparameter, but the degeneracy is always strong. 

In the text 
Fig. D.1 χ^{2} map with the fit departures and arrivals for the model. The color represent the χ^{2} for every pair of chromatic parameters and d_{env}. The white dots are the starting pair for image reconstructions. The red points represent the found parameters for the final image. The fit always ends in the χ^{2} minimum but far from the correct values ( and d_{env} = 1). The black solid lines link the start and the final values of one reconstruction. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.