© ESO, 2010

1. Introduction

When dealing with optical/IR interferometric data, one needs to invoke a model for the understanding of the astrophysical object under consideration. This is because of (1) the low coverage of the uv-plane and most of the time because of the lack of the visibility phase; and (2) because our aim is to extract physical parameters from the data. This is particularly true for the Mid-Infrared Interferometric Instrument (MIDI, Leinert et al. 2003) at the Very Large Telescope Interferometer (VLTI), on which our considerations will be focused. Some pure geometrical information can be recovered through a simple toy model such as Gaussians (see e.g.  Leinert et al. 2004; Domiciano de Souza et al. 2007).

However, this approach does not give any insights into the physical nature of the object. One would dream of having a fully consistent model to characterise the object under inspection. In many cases, if not all, a fully consistent model is out of reach and one uses at least a consistent treatment of the radiative transfer. Models based for instance on the Monte Carlo method are very popular (see e.g.  Ohnaka et al. 2006; Niccolini & Alcolea 2006; Wolf et al. 1999) for this purpose. Still, the medium density needs to be parameterised and it is not determined in a self-consistent way. For massive stars for instance, it would be necessary to take into account non-LTE effects including both gas and dust emission of the circumstellar material as well as a full treatment of radiation hydrodynamics. Fitting interferometric data this way is as yet impossible because of computing time limitations.

Of course, solving at least the radiative transfer in a self-consistent way is already very demanding for the computational resources. Consequently, model parameters cannot be determined in a fully automatic way and the model fitting process must be carried out mostly by hand, or automatised by systematically exploring the parameter-space , the “chi-by-eye” approach mentioned in Press et al. (1992). The followers of this approach consider the “best fitting” model as their best attempt: a model that is compatible with the data. It is admittedly not perfect, but it is in most of the cases the best that can be done given the difficulty of the task. It is remarkable that a thorough χ² analysis of VLTI/MIDI data of the Herbig Ae star AB Aurigae has been performed by di Folco et al. (2009) which remains to date one of the most achieved studies of this kind. From the χ² analysis, formal errors can be derived and at least the information concerning the constraints for the physical parameters can be quantified. Qualitative information about the correlation of parameters can be pointed out.

The next step after the toy models for the physical characterisation of the astrophysical objects can be made from the pure geometrical model towards the self-consistency by including and parameterising the object emissivity in the analysis. For instance, Lachaume et al. (2007) and Malbet et al. (2005) use optically thick (i.e. emitting as black bodies) and infinitely thin discs to model the circumstellar environment of pre-main-sequence and B[e] stars. Of course this approach has some restriction when modelling a disc: for instance it cannot handle nearly edge-on disc and an optically thin situation.

We propose an intermediate approach: between the use of simple geometrical models and sophisticated radiative transfer solvers. Indeed, it is a step backwards from the “self-consistent” radiative transfer treatment, which is in most cases too advanced with regard to the information provided by the interferometric data. For this intermediary approach, we assume a prescribed and parameterised emissivity for the medium. Our purpose is to derive the physical parameters that characterise this emissivity. In the process, we compute intensity maps and most particularly visibility curves from the knowledge of the medium emissivity with a fast ray-tracing technique (a few seconds depending on map resolution), taking into account the particular symmetries of a disc configuration. Then, the model fitting process can be undertaken in an automatic way with standard methods (see e.g. Levenberg 1944; Marquardt 1963). The techniques we present are designed to be quite general and not tailored to any particular emissivity except for the assumed axisymmetry of the problem under consideration.

Our purpose is twofold. On one hand – as already mentioned – we aim to estimate physical parameters and their errors characterising the circumstellar dusty medium under consideration with as few restrictive assumptions as possible; at least within the obvious limitations of the present model. On the other hand our purpose is to provide the user of a more detailed model, such as a Monte Carlo radiative transfer code, with a first characterisation of the circumstellar matter to start with.

In Sect. 2 we describe the general framework of the proposed ray tracing technique. In particular how to derive the observable from the astrophysical object emissivity. In Sect. 3 we describe the numerical aspects that are specific to the present ray-tracing technique. In particular, the use of a quadtree mesh and the symmetries that allow us to speed up the computation are detailed. In Sect. 4 we focus our attention on the circumstellar disc of B[e] stars and describe a parametric model of the circumstellar environment. In Sect. 5 we analyse artificial interferometric data generated both from the parametric model itself and from a Monte Carlo radiative transfer code (Niccolini & Alcolea 2006). Our purpose is not to fit any particular object, but to present our guideline to the following question: which physical information can we get from the data? A discussion of our results and the conclusions of our work are given in Sects. 6 and 7 respectively.

2. The ray-tracing technique

We describe here the FRACS algorithm, developed to study stars with CSEs from mid-IR interferometric observables (e.g. visibilities, fluxes, closure phases). Although FRACS could be extended to investigate any 3D CSE structures, we focus here on the particular case of axisymmetrical dusty CSEs. This study is motivated by the typical data one can obtain from disc-like CSE observed with MIDI, the mid-IR 2-telescope beam-combiner instrument of ESO’s VLTI (Leinert et al. 2003).

2.1. Intensity map

Intensity maps of the object are the primary outputs of the model that we need to compute the visibilities and fluxes that are directly compared to the observations. For this purpose, we integrate the radiation transfer equation along a set of rays (ray-tracing technique) making use of the symmetries of the problem (see Sect. 3 for details).

The unit vector along the line of sight is given by , i being the inclination between the z-axis and the line of sight and , ŷ the unit vector along the x et y-axis of a cartesian system of coordinates (see Fig. 1), referred to as the “model system” below. The problem is assumed to be invariant by rotation around the z-axis. We define a fictitious image plane by giving two unit vectors Ŷ =  −ŷ   cosi + ẑ   sini and . This particular choice is made making use of the axisymmetry of the problem. Note that for this particular coordinate system (X,Y) the disc position angle (whenever i ≠ 0) is always defined as 90°. The actual image plane, with the Y′ and X′ axis corresponding respectively to North and East, is obtained by rotating the axis of our fictitious image plane by an angle , where PA_d is the position angle of the disc with respect to North.

The dust thermal emissivity at wavelength λ and position vector r is given by (1)where is the absorption coefficient and B_λ(T(r) the Planck function at the medium temperature T(r) at r. is defined as , where is the absorption cross section and n(r) the number density of dust grains at r.

We neglect the scattering of the radiation by dust grains, optimising our approach to long wavelengths (from mid-IR to radio). This assumption simplifies the radiative transfer equation by removing the scattering term.

We obtain the intensity map at position (X,Y) in the image plane (inclined by i) and at wavelength λ by integrating the transfer equation along the particular ray that passes through the considered point of the image plane. Defining r_s(X,Y,i) (simply r_s for short) as the position vector along a ray, given in the model system of coordinates by (2)and by introducing the optical depth at wavelength λ and position s along the ray by (3)we obtain (4)where the extinction coefficient because scattering is neglected.

We assume that the CSE is confined within a sphere of radius R_out, s varies consequently from to (R² = X² + Y²) in Eq. (4) and in the definition of a ray Eq. (2). This hypothesis can be relaxed without altering the present considerations and the domain of integration of Eq. (4) suitably chosen.

If some radiation sources (e.g. black body spheres) are included in the analysis, an additional term must be added in Eq. (4) whenever a particular ray intersect a source. For a source with specific intensity this additional term is given by , s^(s) being the distance at which the ray given by X, Y and i (see Eq. (2)) intersects the outermost (along the ray) source boundary. In that case the lower integration limit in Eq. (4), that is , must also be replaced by s^(s).

2.2. Interferometric observables

From the monochromatic intensity maps at wavelength λ (Eq. (4)) we obtain both the observed fluxes F_λ and visibilities V_λ for an object at distance d, (5)and (6)where V_λ is obtained for a given baseline specified by its projected length B (on the sky, i.e. (X′,Y′) coordinates) and its polar angle PA from North to East (direction of the Y′ axis). Δ and j represent, respectively, PA_d − PA and .

3. Numerical considerations

We seek to produce intensity maps within seconds¹ and we aim for our numerical method to be sufficiently general in order to deal with a large range of density and temperature structures. Given these two relatively tight constraints, the numerical integration of Eq. (4) is not straightforward.

For example we have tested that the 5th order Runge-Kutta integrators of Press et al. (1992) with adaptive step-size (as discussed in Steinacker et al. 2006) doest not suit our constraints. Indeed, the step adaption leads to difficulties if sharp edges (e.g. inner cavities) are present in the medium emissivity.

3.1. Mesh generation

Regarding the above mentioned constraints and the different numerical approaches tested, we found that Eq. (4) is more efficiently computed with an adaptive mesh based on a tree data structure (quadtrees/octrees). The mesh purpose is twofold: first, it must guide the computation of Eq. (4) and distribute the integration points along the rays according to the variations of the medium emissivity; second, within the restriction of axis-symmetrical situations, the mesh must handle any kind of emissivity. Quad/octree meshes are extensively used in Monte Carlo radiative transfer codes (e.g. see Bianchi 2008; Niccolini & Alcolea 2006; Jonsson 2006; Wolf et al. 1999); the mesh generation algorithm is thoroughly described in Kurosawa & Hillier (2001).

The mesh we use is a cartesian quadtree. Cartesian refers here to the mesh type and not to the system of coordinates we use. Indeed, the mesh is implemented as a nested squared domain (cells) in the ρ − |z|  plane (). The whole mesh is enclosed by the largest cell (the root cell in the tree hierarchy) of size R_out in ρ and  |z| . The underlying physical coordinate system is cylindrical (with z > 0) and the mesh cells correspond to a set of two (for z > 0 and z < 0) tori, which are the actual physical volumes.

The mesh generation algorithm consists in recursively dividing each cell in four child cells until the following conditions are simultaneously fulfilled for each cell in the mesh (see Kurosawa & Hillier 2001, for more details): where V_ξ is the volume of cell ξ, V_tot is the volume of the root cell and α, β and η are parameters controlling the mesh refinement.

In the present work α and β have been fixed to 1, but higher values can be useful for some particular situations where the generated mesh must be tighter than the mesh generated directly from the and T variations. Typically, these situations show up for high optical depths (in this paper, optical depth values do not exceed  ≃1 at 10   μm along the rays). The practical choice of α, β and η is obtained from a compromise between execution speed and numerical accuracy of the Eq. (4) integration (e.g. typical values of η range from 10^-5 to 10^-4). When dealing with optically thick situations, a supplementary conditions can be added to Eqs. (7) and (8) in order to prescribe an upper limit to the cell optical depth. For instance, making use of the computation of the integral in Eq. (7), one can add the following criterion for cell ξ (whose centre is (ρ_ξ,z_ξ) and size Δ_ξ) (9)where Δτ_lim is the prescibed upper limit to the cell optical depth.

For the moderate optical depths reached in this work, with values of η down to 10^-5 and Δτ_lim set to 10^-2, the criteria of Eqs. (7) and (8) are the leading conditions to the mesh refinement.

Figure 2 shows the mesh obtained in the particular case of a B[e] circumstellar disc (see Sect. 4) for models whose parameters are given in Table 2 (see caption for more details).

The volume integrals in Eqs. (7) and (8) are estimated by Monte Carlo integration. For a quantity f(r) and for the cell ξ the integral d³r is approximated by (10)where we made explicit use of the mesh coordinates and where (ρ_k,z_k) with k = 1,···   ,N are chosen randomly and uniformly within the cell domain.

3.2. Symmetries

We can make use of the CSE symmetries to reduce the computation domain of an intensity map from Eq. (4) to only a fourth of it and consequently reduce the computation time.

Recalling the definition of a ray (Eq. (2)), we have two noticeable identities for any disc physical quantity Φ (e.g. , , T, n, ...) depending on rwhere Eq. (11) expresses the disc symmetry with respect to the y − z plane and Eq. (12) the point symmetry with respect to the origin of the model system of coordinates.

From the above identities it is straightforward to deduce their counterpart for the intensity map where . Note that the exponential factor in Eq. (14) has to be evaluated when computing I_λ(X,Y,i) anyway; no extra effort is required to derive I_λ(X,Y,i) from I_λ(X, − Y,i) except for the multiplication of I_λ(X, − Y,i) by this factor.

3.3. Intensity map

The fictitious image plane is split into a set of pixels whose positions X_j and Y_k are given by where Δ_X = Δ_Y is the pixel size in X and Y, and N is the number of pixels in X and Y, and where (17)where δ =  −1 for N even and δ = 0 otherwise and “ ÷ ” stands for the integer division. Taking into account the symmetries mentioned in Sect. 3.2 only a fourth of the pixels need to be considered.

The evaluation of the integral in Eq. (4) is carried out for each pixel (X_j,Y_k) and along the ray r_s(X_j,Y_k,i). The intersection points of the ray with the cell boundaries corresponds to a set of distances along the ray defined as where n_cells is the number of cells encountered along the ray, and Δs_l the distance crossed within the lth cell.

We estimate numerically the optical depth τ_λ(X,Y,i;s), defined in Eq. (3), via the midpoint rule quadrature by (20)where we defined for l = 0,···   ,n_cells − 1.

The numerical estimate of I_λ(X_j,Y_k,i) is obtained by (21)The results for all j and k can then be obtained from the discrete counterpart of the symmetry relations (13) and (14)

3.4. Interferometric observables

From the numerical estimate of I_λ(X_j,Y_k,i) given above we obtain (similarly to Eqs. (5) and (6)) the numerical fluxes and visibilities, which can be directly compared to the observed data. The numerical estimate of these quantities is again obtained through the mid-point rule.

The numerical flux F_λ(i) is computed by (24)The complex visibility is approximated numerically by (25)where B = (BcosΔ,BsinΔ) and R_kl = (X_k,Y_l).

3.5. Artificial data generation

The procedure described below aims to mimic the observables of the VLTI/MIDI instrument: the flux F_λ (Eq. (24)) and the modulus of the visibility  |V_λ|  (Eq. (25)). The wavelengths and baselines chosen for the artificial data generation correspond to accessible values to VLTI/MIDI with the Unit Telescopes (UTs): λ_j = 7,8,9,10,11,12, and 13   μm (j = 1,···   ,n_λ; n_λ = 7), and (B_k,PA_k) as shown in Table 1 (k = 1,···   ,n_B; n_B = 18). These values amount to 126 points covering the uv-plane.

For a given intensity map at λ_j, F_λj and  |V_λj|  are taken as the expectation values of the simulated data. The observed flux  is then generated assuming a Gaussian noise with an RMS (root mean square) corresponding to 10% relative error σ_F(j) = 0.1 × F_λj.

The artificial observed visibility amplitudes are obtained as (26)where ΔV_k is a wavelength independent shift that mimics the error in the observed visibilities, introduced by the calibration procedure commonly used in optical/IR interferometry. For each (B_k,PA_k), ΔV_k is computed assuming a Gaussian noise with an RMS corresponding to 10% relative error (typical for VLTI/MIDI) σ_V(k) = 0.1 ×  ⟨|V_λ(B_k,PA_k)|⟩ , where  ⟨|V_λ(B_k,PA_k)|⟩  is the wavelength mean visibility modulus.

3.6. Model fitting and error estimate

We describe here the procedure adopted in order to simultaneously fit observed fluxes and visibilities using FRACS models defined by a given set of input parameters. This procedure is applied to artificial data in the next sections.

In order to quantify the discrepancy between the artificial observations ( and ) and the visibilities and fluxes from a given model (|V_λj(B_k,PA_k)|  and F_λj) we use the χ² like quantities (27)and (28)To take into account both the mid-IR flux and the visibilities on the same level in the fitting process, we minimise the following sum (29)In the discussion below about the parameter and error determination we use the reduced χ² defined by (for n_free free parameters).

From a minimising algorithm the best-fit model parameters can then be found by determining the minimum : . The “error” estimate is obtained from a thorough exploration of the parameter space volume, defined by a contour level , where has been chosen equal to 1. This volume can be interpreted as a confidence region. The quantity defined in Eq. (29) is a weighted sum of χ² variables whose cumulative distribution function can be approximated by a gamma distribution (see Feiveson & Delaney 1968) with the same mean and variance. It is then possible to obtain a rough estimate of the confidence level associated with the confidence region given approximately by  ≃2σ.

The size of the confidence region is determined by considering all possible pairs of parameters for a given fitted model and computing maps for each. The procedure to estimate the errors can be summarised as follows:

• For a given map, i.e. for a given couple of parameters among the n_free × (n_free − 1)/2 possibilities, we identify the region bounded by around the minimum of this particular map.

• The boundaries of the projection of these regions on each of the two parameter axis considered are recorded for each map.

• The final errors on a given parameter are taken as the highest boundary values of the projected regions over all maps.

4. Astronomical test case: sgB[e] stars

In the following sections we apply FRACS to a theoretical interferometric study of dusty CSE of B[e] supergiants (sgB[e] in the nomenclature of Lamers et al. 1998). However, we emphasise that FRACS is in no way restricted to this particular class of objects.

sgB[e] stars reveal in particular a strong near- or mid-IR excess caused by hot dust emission. There is evidence (e.g. Zickgraf et al. 1985) that the stellar environment, and in particular dust, could be confined within a circumstellar disc. Our purpose is to characterise this class of objects and derive not only geometrical parameters (e.g. inner dust radius, disc position angle and inclination) but also physical parameters such as temperature gradients, dust formation region, material density, ...

The physical description of the CSE chosen for our study is the wind model with equatorial density enhancement. This is a classical CSE model commonly adopted for sgB[e] (e.g. Porter 2003).

In order to compute the model intensity maps we need to parameterise the emissivity of the disc. Consistently with FRACS assumptions, we consider only dust thermal emission without scattering by dust grains and the gas contribution to the medium emissivity. In the rest of this section we characterise the emissivity by describing the dust density law, the absorption cross section, and the temperature structure of the CSE.

4.1. Mass loss and dust density

Dust is confined between the inner and outer radius R_in and R_out respectively. We assume a stationary and radial mass loss; physical quantities will consequently depend only on the radial component r and the co-latitude θ. The disc symmetry axis coincides with the z axis of the model cartesian system of coordinates. The mass loss rate and velocity parametrisations are simplifications of the one adopted by Carciofi et al. (2010), and we refer the reader to their work for a complete description (see also Stee et al. 1995, for a similar description).

The mass loss rate per unit solid angle, at co-latitude θ, is parameterised as follows (30)with the help of two dimensionless parameters A₁ and m.

Even though our computations make no explicit use of the radial velocity field v_r(θ) (assumed to have reached the terminal velocity v_∞(θ) in the region under considerations, i.e. v_r(θ) ≈ v_∞(θ)), the dust density depends on v_r(θ) parameterised in a similar fashion (31)where we have introduced the supplementary dimensionless parameters A₂. From Eqs. (30) and (31) we see that A₁ and A₂ are the relative differences of the values of and v_r(θ) at the equator and the pole (relatively to the pole).

From the mass continuity equation one obtains the number density of dust grains (32)where n_in is the dust grain number density at R_in in the disc equatorial plane. In Eq. (32), the parameter m controls how fast the density drops from the equator to the pole, defining an equatorial density enhancement (disc-like structure).

Consistent with the accepted conditions for dust formation (Carciofi et al. 2010; Porter 2003) we assume that the dust can survive only in the denser parts of the disc. We thus define a dusty disc opening angle Δθ_d determined by the latitudes for which the mass loss rate has dropped to half of its equatorial value: (33)To summarise, the dust grains only exist (i.e., n(r,θ) ≠ 0) in the regions bounded by R_in ≤ r ≤ R_out and by .

4.2. Dust opacities

The absorption cross section for the dust grains is obtained from the Mie (1908) theory. The Mie absorption cross sections are computed from the optical indices of astronomical silicate (Draine & Lee 1984). Note that since scattering is neglected, , with being the extinction cross section.

For a power-law size distribution function according to Mathis et al. (1977) the mean cross sections (e.g. for ) are given by (34)where a_min and a_max are the minimum and maximum radii for the dust grains under consideration and β is the exponent of the power-law. The computation of the cross section in Eq. (34) was performed with the help of the Wiscombe (1980) algorithm.

4.3. Temperature structure

The dust temperature is assumed to be unique (i.e. independent of grain size) and described by a power-law (35)where T_in is the temperature at the disc inner radius R_in. We note that γ is not necessarily a free parameter because in the optically thin regime (large wavelength and radius) the temperature goes as with δ ≃ 1 (see Lamers & Cassinelli 1999).

4.4. Central continuum emission

The continuum emission from the central regions is composed by the emission from the star and from the close ionised gas (free-free and free-bound emission). This central source emission is confined to a small region of radius R_s (≪ R_in), which is unresolved (angular sizes of a few milliarcseconds) by mid-IR interferometers. Thus, in our modelling R_s is simply a scaling factor of the problem fixed to a typical radius value for massive stars. The specific intensity (in W   m^-2   μm^-1   str^-1) of this central source is parameterised as follows (36)where is the specific intensity at a reference wavelength λ₀ (= 10   μm in the following), and α gives the spectral dependence of the continuum radiation. In the mid-IR its value is expected to lie between α = 4 (pure black body) and α ≃ 2.6 (free-free emission) for an electron density proportional to r^-2 (Panagia & Felli 1975; Felli & Panagia 1981).

5. Study of the tested models

Following the description in the last sections we describe here the chosen sgB[e] model parameters used to simulate VLTI/MIDI observations (visibilities and fluxes) and the corresponding analysis, i.e. model fitting, using FRACS. The list of chosen parameters is summarised in Table 2. Two types of numerical tests are presented. Firstly, synthetic mid-IR interferometric data are generated from FRACS itself. In that way, it is possible to determine what information the mid-IR interferometric data contain under the optimistic assumption that we do have the true model. Secondly, this study is supplemented by the comparison of FRACS to a Monte Carlo radiative transfer computation. This confirms that FRACS can indeed mimic, under appropriate conditions, the results of a more sophisticated code as seen from the mid-IR interferometric eye.

5.1. Parameter description

The distance to the simulated object has been fixed to d = 1   pc, which is a typical distance for Galactic sgB[e].

The inner radius R_in = 30   R_s = 1800   R_⊙ value was chosen by considering the location of the hottest dust grains (see Lamers & Cassinelli 1999) with a condensation temperature of 1500   K assumed to be the T_in value. The value of R_out cannot be determined from the mid-IR data and has been fixed to 3000   R_s = 1.8 × 10⁵   R_⊙. The temperature gradient γ was fixed to 0.75 according to Porter (2003). PA_d was fixed to 125°.

The central source emission is supposed to have a radius R_s = 60   R_⊙. We recall that the central region is unresolved by the interferometer and that its radiation describes both the stellar and inner gas contribution to the continuum mid-IR emission. The specific intensity of this central source has been chosen to be 6500   W   m^-2   μm^-1   str^-1. If the central source was a pure blackbody this value would correspond to the 10   μm emission of a blackbody with an effective temperature around  ≃8000   K. However, this central emission is not a pure blackbody, and we adopt the spectral dependence of the central source emission to be α = 3, which is a compromise between α = 4 for a pure blackbody and a value of  ≃2.6 for free-free emission (Panagia & Felli 1975; Felli & Panagia 1981).

Spectroscopic observations of Hα and forbidden line emissions from B[e] CSE (Zickgraf 2003) reveal that typical values for A₂ are expected to range from  −0.95 to  −0.75. We adopt the value  −0.8 in our models. According to Lamers & Waters (1987), the values of A₁ range from 10² to 10⁴ in most cases (though values as low as 10 are not excluded). With this high value of A₁ the factor 1 + A₁(sinθ)^1/m in Eq. (32) of n(r,θ) is approximatively given by A₁(sinθ)^1/m for all pertinent values of θ, i.e. those close to π/2 within the disc. This leads to an evident degeneracy in n_in × A₁ in n(r,θ): we are only sensible to the product of the two parameters as a scaling factor for the density. Therefore, the value of A₁ is assumed to be fixed to 150.

To define the dust opacities the chosen value for β is that of Mathis et al. (1977), i.e. β =  −3.5. Because some sgB[e] show weak 9.7   μm silicate features in their spectrum (e.g. Porter 2003; Domiciano de Souza et al. 2007) we chose to use large grains in our test models: a_min = 0.5   μm and a_max = 50   μm. However, with this particular choice of large grains, the average albedo from 7 to 13   μm is 6.4%, with the highest value reached at 7   μm. We have checked with a Monte Carlo (MC) simulations (see Sect. 5.3) that the effect of scattering on our primary observables, visibilities, and fluxes is indeed negligeable by comparing the results obtained by switching the scattering process off and on². The mean relative differences are 3.5% and 3.0% for the visibilities and the fluxes respectively. These values must be compared to the effect of random noise in the MC simulation, estimated to be of the same order and to experimental errors, typically  ~10% for the visibilities and fluxes. We underline that whenever the albedo can be neglected, it is theoretically safe to compute visibilities and fluxes from the consideration presented in Sect. 2, in any other situations the effect of scattering on the observable must be carefully tested.

The parameters n_in, m and i were set to different values defining 12 test models to be analysed from their corresponding simulated data. Two n_in values (0.015 m^-3 and 0.15 m^-3) have been chosen in order to have an approximate disc-dust optical depth in the equatorial plane (from R_in to R_out) close to  ≃0.1 and  ≃1 in the wavelength range considered (from 7   μm to 13   μm). These values corresponds to a mass loss rate of Ṁ ≃ 2.5 × 10 ^{−7··· −6}   M_⊙   yr^-1. Two m values were chosen corresponding to a wide and a narrow opening angle, i.e. Δθ_d = 10° and Δθ_d = 60°. Three inclinations i were tested (20°, 50°, and 90°) corresponding to discs seen close to pole-on, intermediate inclination, and equator-on. These values of n_in, m, i, together with the parameters fixed above, define 12 test models that will be studied below.

From these 12 test models we have generated 12 sets of artificial VLTI/MIDI observations (visibilities and fluxes) following the procedure described in Sect. 3.5. We do not aim to present an exhaustive revue of all types of sgB[e] CSE. Rather, we focussed on the analysis of the parameter constraints one can hope to obtain from present and near-future mid-IR spectro-interferometry. The quantitative estimate of these constraints is derived from a systematic analysis of the variations with the parameters.

In our model fitting and analysis we concentrate on 10 free parameters (n_free = 10) that can be set into four different groups:

• the geometrical parameters: PA_d, i and R_in;

• the parameters related to the central source: and α;

• those describing the temperature structure: T_in and γ;

• and the number density of dust grains: A₂, Δθ_d (or equivalently m) and n_in.

The remaining parameters of the model are in general loosely constrained by mid-IR interferometric observations so that we kept them fixed to the values described above.

5.2. Model fitting and analysis of the 12 test models

We describe here the data analysis procedure adopted to study our 12 test models. The results of our analysis are summarised in Tables 3 and 4, and their physical interpretation is presented in Sect. 6.

As a first step we chose 2 of the 12 models, hereafter called models (a) and (b), to be exhaustively studied from a complete model fitting procedure. As an example we show the simulated observed mid-IR fluxes and visibilities for model (a) and (b) in Figs. 3−5. The parameters of models (a) and (b) are those of Table 2 with Δθ_d = 60°, i = 50° and n_in fixed to the values 0.015   m^-3 and 0.15   m^-3 respectively. These two models are those presenting some of the best constrained model parameters for the dust CSE. On the other hand, the contribution of the central regions to the total flux and visibilities is quite different in models (a) and (b) (see discussion in Sect. 6).

The study of models (a) and (b) have thus been performed as for real interferometric observations. The best-fit values of the parameters have been obtained by the Levenberg-Marquardt algorithm with a stopping criterion corresponding to a relative decrease in of 10^-3.

The errors on each model parameter have been obtained following the methods described in Sect. 3.6. The maps have been computed with a resolution of 21 × 21 around the best-fit values of the parameters. The map sizes have been adjusted in order to enclose the contour. This adjustment was performed until an upper limit for the map size of 100% of the best-fit parameter values was reached. This amounts to the computation of 3.969 × 10⁴ different models. The results, namely the mean relative error up to 100%, for these two particular models are summarised in Table 3.

The other ten models (numbered from 1 to 10 in Table 4) have been used in order to get some quantitative (but limited) information about how the uncertainties of the fitted parameters evolve as a function of three disc characteristics: its optical depth (τ by means of n_in parameter), its inclination (i) and its opening angle (Δθ_d, controlled by m). To perform this study we have decided to limit the exploration of the space parameter in a relative range of 25% on both sides from the model parameters. In order to reduce the computation time, the maps were not generated around the best-fit parameters which would have required to compute several thousands models more but around the true parameters themselves. This procedure has the supplementary advantage that we do not rely on any specific minimisation algorithm. We checked that estimating the best-fit parameters from the true ones is reasonable within a few percents using the Levenberg-Marquardt algorithm with a stopping criterion corresponding to a relative decrease in of 10^-3. The resolution of the maps have been reduced to 15 × 15. The total number of models to be computed is as large as 1.0125 × 10⁵.

5.3. Comparison with a Monte Carlo simulation

We generated synthetic data with the help of a Monte Carlo (MC) radiative transfer code (Niccolini & Alcolea 2006) for model b (see above) for the seven wavelengths considered in the problem and the baselines of Table 1. Again, the adopted procedure to generate the mid-IR interferometric data follows the considerations of Sect. 3.5. In the MC code, the source of photons is described by a blackbody sphere of radius R_s = 60   R_⊙ and an effective temperature of T_eff = 8000   K. The temperature of the CSE is not prescribed but computed from the Lucy (1999) mean intensity estimator. This choice of T_eff gives at the inner radius of model b a dust temperature of  ≃1150   K lower than the sublimation temperature. In this way, we can test if in the fitting process using FRACS, a spurious effect might not lead the minimisation algorithm to reach the upper limit for T_in of 1500   K, corresponding to the adopted dust sublimation temperature.

We obtained the best fitting parameters for the CSE model described in Sect. 4 with FRACS. For a comparison with the MC code, α has been set and fixed to 4 corresponding to the value of a blackbody. Depending of the disc optical depth, the temperature structure may show two separate regimes corresponding to (1) the inner regions with the strongest temperature gradient, optically thick to the stellar radiation and (2) the outer regions optically thin to the disc radiation with a flatter temperature gradient. In order to determine if mid-IR interferometric data are sensitive to two temperature regimes, we tested the effect of two parameterisations of the temperature structure: the unique power-law of Eq. (35) and a generalisation to two power-laws with a transition radius, R_T, and a second exponent γ′(37)for r ≥ R_T.

The best-fitting parameters for both parameterisations are shown in Table 5. The images of the disc at 10   μm generated with the MC code and their corresponding FRACS counterpart (best-fitting model) are shown in Fig. 6 for comparison.

6. Discussion

We first discuss the uncertainties in the parameters derived for the 12 models studied in Sect. 5.2. For each model we divided the parameters into three groups associated to a given level of constraints expressed by the relative errors: below 10%, between 10% and 25%, and above 25%. This information is summarised in Table 4. The exact relative errors for the two models studied in detail (models a and b) are shown in Table 3. Then, we analyse the results of Sect. 5.3 obtained from the best-fit of the data simulated with the MC radiative transfer code.

6.1. Central source

Table 4 shows that the central source parameters ( and α) can only be constrained with relative uncertainties  ≥10% for all test models. A deeper and more quantitative investigation of these parameters can be obtained from models (a) and (b). From Table 3 it can be seen that and α are much better constrained for model (a) (27% and 13%, resp.) than for model (b)(relative errors  ≃100%).

The key quantity for a good constraint for the central source parameters ( and α) is simply the relative contribution of the flux of the central source to the total flux of the object (source and disc). Indeed, the models in Table 3 only differ by this relative flux contribution of 5.3% in model (a) (τ = 0.1), while it is only 0.7% in model (b) (τ = 1).

Our analysis thus shows that interferometric data can constrain and α with a relative precision of  ≃15%−30% even when the central source contributes to (only) a few percent of the total mid-IR flux.

6.2. Geometrical parameters

The parameters PA_d, i, and R_in are those usually estimated from simple geometrical models (e.g. ellipses, Gaussians). However, their determination from geometrical models is quite limited, in particular for i, for which only an estimate can be derived from the axis-ratio of an ellipse, for example. In addition, the estimate of i from a simple analytical model such as a flat ellipse is only valid for configurations far from the equator (intermediate to low i). The use of a more physical and geometrically consistent model such as FRACS allows us to relax this constraint and makes the determination of i possible for all viewing configurations.

As expected, PA_d and i are better determined if the inclination of the disc with respect to the line of sight is away from pole-on (high i). In Fig. 8 we can clearly see this behaviour from the maps involving PA_d and i. Moreover, the uncertainties on PA_d and i do not seem to be strongly dependent on τ (equivalently n_in) and Δθ_d (equivalently m) for all models.

The inner dust radius R_in is not strongly dependent on any parameter (τ, Δθ_d or i), being very well constrained (better than 10%) for most tested models.

6.3. Temperature

The parameters related to the temperature structure of the CSE, T_in and γ, are well constrained in most models, with relative errors below 20% and 12% for both models (a) and (b). Indeed, γ has a strong impact on the IR emission across the disc, and consequently this parameter has a direct influence on the visibilities (see Figs. 4, 5). T_in has a lower influence, compared to γ, on the shape of the monochromatic image (radial dependence of intensity) and can be mainly considered as a scaling factor to it. On the other hand, the mid-IR flux imposes stronger constraints on T_in. From Table 4 we see that the CSE’s temperature structure is not highly dependent on τ (n_in) and Δθ_d for tested models.

6.4. Number density of dust grains

The parameters related to the density law, that is to say m, n_in and A₂, seems to be rather poorly constrained from the mid-IR data alone. From the results of model (a) and (b) corresponding to an intermediary inclination i = 50°, we found that only n_in is constrained somewhat moderately with a mean relative error of 46%. For m and A₂, according to the results of Table 3 it seems that nevertheless, upper limits to their values can be determined. Note that because A₂ is bounded (−1 ≤ A₂ ≤ 0) the mean relative errors, 53% and 59% for model (a) and (b) respectively, correspond approximatively to the limit values of A₂, which is consequently not constrained.

Table 4 confirms this trend for m, n_in and A₂ at least for the situations explored via the models presented here. From all maps computed within  ±25% of the true value, we always found that the mean relative error to these parameters is larger than 25% with no hint that it could be close to these limits.

From Fig. 9, comparing the maps for all pairs of n_in, m and A₂, we can see that the contours get sharper around the minimum value for model (a) (corresponding to lower optical depths along the line of sight) than for model (b). Indeed, the constraints on n_in and m are improved for lower optical depths, or equivalently for lower disc masses. Indeed, when the disc mass (or optical depth) decreases, the flux (mid-IR flux, intensity maps) emitted by the disc reflects the mass of the disc, while for high optical depths we only probe the regions of the disc very close to the projected surface revealed to the observer. A₂, however, is unaffected by the change in disc mass and remains undetermined anyway. From Fig. 9 it can be seen that n_in, m and A₂ are strongly correlated. This is expected from the expression of the density (see Eq. (32)) depending on these parameters. However, this dependence and the final correlation between these parameters are related through the computation of the visibilities and the mid-IR flux, as well as the comparison to the data and is, therefore, not straightforward.

To improve the situation concerning n_in, m and A₂, the mid-IR data can be supplemented by other types of observations such as for instance spectroscopic data, from which one can better determine A₂ (e.g. see Chesneau et al. 2005). We tested the effect of fixing the value of A₂, or equivalently of assuming that A₂ is fully determined, in the process of estimating the errors of the other parameters. For model (a), the relative errors on n_in and m go down to 33% and 71% respectively while for model (b), n_in and m are determined with an accuracy reaching 95% and 78% respectively. The precision to which other parameters are determined is not affected by the determination of A₂.

We also tested the influence of the determination of n_in, m and A₂ on other parameters by fixing their values and estimatng the relative errors on the remaining parameters for model (a) and (b). Only and T_in are more strongly affected by the determination of n_in, m and A₂: for model (a) (resp. model b) gets determined down to 19% (resp. 80%) and T_in down to 9% (resp. 18%). The influence is stronger with lower disc mass (model a compared to model b). This effect can be explained because if we have a good determination of the disc mass because we know n_in, m, and A₂, the determination of the parameters that scale the source and disc fluxes is improved accordingly for the visibility and the mid-IR flux.

n_in, m and A₂ shape the density structure of the circumstellar medium. Though they are not well constrained, they certainly have a strong influence on the temperature structure, which in turn is very well constrained. For the particular case of sgB[e] circumstellar discs, a natural evolution of FRACS is to include the direct heating of the medium by the central source of radiation assuming that the disc is optically thin to its own radiation. The temperature structure would not be parameterised, and its good determination would certainly put better constraints on n_in, m and A₂, while keeping an affordable computation time for the model-fitting procedure. This will be the purpose of a subsequent work.

Finally, one can derive a ranking of the parameter constraints according to two criteria: first the parameter must be constrained within the prescribed limits (100% for model a and b and 25% for model 1 to 10) and second the mean relative error must be as low as possible. The best-fitted parameters, most of the time according to these criteraria are by decreasing order of best determination: R_in, PA_d, γ, T_in, i, α, , n_in, m and A₂. This tendency can be seen in Table 4.

6.5. best-fit to the MC simulation

The values obtained for the two types of temperature parametrisations (one and two power-laws) are quite similar: 0.79 and 0.73 respectively. Regarding the data, both temperature parametrisations are indeed acceptable. In addition, these results show that we can actually obtain very good fits from data sets based on more physically consistent scenarios. A complete error analysis and study of the parameter determination has been presented in the previous sections for data generated from FRACS and will not be repeated here. In particular, parameter confidence intervals, from which errors were derived, have already been estimated. Here, we will instead focus on the true errors, i.e. the differences between the true model parameters and the best-fitting values for the parameters (see Table 5). The two types of errors must not be confounded. The true errors reflect the capability of FRACS to mimic the mid-IR interferometric data regarding the information it provides. Of course, with the sparse uv-plane coverage inherent to this kind of data as well as the experimental noise, one should not expect a full agreement of the fitted and the true parameters: they are indeed different.

From Table 5, we see that the geometrical parameters, PA_d, i and R_in, can be almost exactly recovered as expected. The source specific intensity , and the parameters related to the density, A₂, n_in and m can be recovered fairly well and have best-fitting values close to the true parameters.

The values of T_in, R_T, γ and γ′ reported as “true” in Table 5 are indeed the values of a fit to the average (over the co-latitude for a given r) computed temperature in the disc. The true relative differences for T_in do not exceed 7% independently of the adopted parameterisation of the temperature (one or two power-laws). The best-fitting values of γ, the inner temperature gradient, obtained with FRACS are very close to the true values with two and one power-law with true relative error of 1% and 7% respectively. This already suggests that the mid-IR data provide information on the inner and hottest region of the CSE, in particular on the inner temperature gradient γ.

Fitting the temperature computed with the MC code with a simple power-law, we obtain γ ≃ 0.64. This value is close to those of the best fitting models, especially with a unique power-law (6% relative difference). For comparison, the actual mean temperature gradient as derived from the MC simulation is  ≃0.60. For this particular data set, the values of γ′ and R_T recovered by FRACS differ by 28% and 45% respectively from the actual values. This again confirms the sensibility of the interferometric data to the temperature structure mostly in the inner (r ≲ R_T) regions of the disc. The best-fitting models (fourth and last columns in Table 5) as well as the MC results are shown in Fig. 7. Regarding the errors (estimated from the results given in Table 3) shown as a shaded area, we can see that both temperature parameterisations are essentially the same and show a better agreement with the MC results in the inner than in the outer regions of the disc.

We considered a “truncated” model with two power-laws (with parameter values listed in the third column of Table 5) in which the CSEs emission for r ≥ R_T, the “outer” regions, has been set to 0. We then compared the visibilities and the fluxes of this truncated model to the same model including the outer region emission. We obtained relative differences, averaged over all considered wavelengths and baselines (Table 1), of 18% and 17% respectively. These relative differences are larger, but are still close to the noise level. For this reason, one cannot expect to obtain much information on the outer temperature gradient γ′, at least for the particular configuration we considered.

7. Conclusion

We proposed and described here a new numerical tool to interpret mid-IR interferometric data. Even though we focussed on the special case of circumstellar disc observations, the numerical techniques have been developed with the aim to be as general as possible. The methods we employ rely on both parameterised physical models and the ray-tracing technique. The need for such a tool is evident because the nature of interferometric data imposes an interpretation through a model of the object to obtain any kind of information. On one hand, Monte-Carlo radiative transfer methods require too much computation time to associate the model-fitting to an automatic minimum search method. On the other hand, purely geometrical function fitting (such as ellipses or Gaussians) are too simple to envisage to obtain physical constraints on the observed disc. Hence, a tool like FRACS fills a blank in the model fitting approach for mid-IR interferometric data interpretation. The main advantages of FRACS are its speed and its flexibility, allowing us to test different physical models. Moreover, an exploration of the parameter space can be performed in different manners and can lead to an estimate of the sensitivity of the fit to the different model parameters, i.e. a realistic error estimate.

We applied these techniques to the special astrophysical case of B[e] star circumstellar environments by generating artificial data in order to analyse beforehand what constraints can be obtained on each parameter of the particular disc model in this work. The techniques will then be applied to real interferometric data of a sgB[e] CSE in a sequel to this paper.

We showed in our analysis that the “geometrical” parameters such as R_in, PA_d and i can be determined with an accuracy  ≲15%. Mid-IR interferometric data give access to a mean temperature gradient: the temperature structure (T_in and γ) can be very well determined (within  ≲20% and  ≲10% respectively). It is possible to have access to the central source emission (with an accuracy  ≳30%) when it has a significant contribution to the total flux of the object (a few % are sufficient). The remaining parameters of our disc model, namely n_in, m and A₂ are not very well constrained by MIDI data alone. n_in is at best determined with an accuracy of about  ≳50% in some cases. If A₂ can be estimated through spectroscopic observations, then the picture about the n_in and m determination improves somewhat.

FRACS can be used mainly for two purposes. First, it can be used by itself to try and determine physical quantities of the circumstellar matter. Admittedly, it is not a self-consistent model, i.e. the radiative transfer is not solved because the temperature structure is parameterised. From the usual habits in the interpretation of interferometric data it is nevertheless a step beyond the commonly use of toy models or very simple analytical models.

This approach has indeed been very successful in the millimetric wavelength range (e.g. see Guilloteau & Dutrey 1998). Second, it can be viewed as a mean to prepare the work of data fitting with a more elaborate model (such as a Monte Carlo radiative transfer code for instance) and to provide a good starting point.

FRACS is a tool that can help in the process of interpreting and/or preparing observations with second-generation VLTI instruments such as the Multi-AperTure mid-Infrared SpectroScopic Experiment (MATISSE) project (Lopez et al. 2006). In this respect, FRACS is not restricted to the mid-IR, and sub-millimeter interferometric data obtained with the Atacama Large Millimeter Array (ALMA) for instance can be tackled.

---

Acknowledgments

We thank the anonymous referee for her/his constructive comments that lead to significant improvements of the manuscript. We also would like to thank Alex Carciofi, Olga Suárez Fernández, Andrei and Ivan Belokogne for fruitful discussions and careful proof reading.

The Monte Carlo simulations have been carried out on a computer financed by the BQR grant of the Observatoire de la Côte d’Azur.

This work is dedicated to Lucien.

References

All Tables

All Figures