A&A 377, 721-734 (2001)
DOI: 10.1051/0004-6361:20011047

Interferometric-Doppler imaging of stellar surface structure

S. Jankov 1,4 - F. Vakili 1 - A. Domiciano de Souza Jr. 1 - E. Janot-Pacheco 2,3


1 - Observatoire de la Côte d'Azur, Département FRESNEL, CNRS UMR 6528, 2130 route de l'Observatoire, Caussols, 06460 St Vallier de Thiey, France
2 - Observatoire de Meudon, DASGAL, 92195 Meudon, France
3 - Instituto Astronômico e Geofísico, Universidade de São Paulo, CP 9638, 01065-970 São Paulo, SP, Brasil
4 - Astronomical Observatory Beograd, Volgina 7, Yugoslavia

Received 7 February 2001 / Accepted 16 July 2001

Abstract
In this paper we discuss the combination of two basic approaches which should potentially generate images of spatially unresolved stars: differential interferometry and classical spectroscopy. Doppler Imaging provides indirect observational information on stellar surface structures by modeling the rotational modulation of the observed flux distribution across spectral lines. Similarly, differential interferometry makes it possible to measure the disturbances of photocentroid location of an unresolved star as a function of wavelength and to deduce the corresponding stellar map. We present the general formalism to reconstruct images from spectroscopy and differential interferometry data for sources with spatially unresolved structures, and we discuss how their combination improves the image reconstructions. This technique, that we call Interferometric-Doppler Imaging (IDI), leads to significant progress in solving some long-standing problems of Doppler Imaging, such as latitude smearing and bias as well as the non-uniqueness of the solution in the special case of an equator-on star. We treat explicitly the most delicate case of non-radial stellar pulsations, for which the cancellation of opposite sign temperature or velocity fields introduces additional difficulties. The performance of the method is demonstrated, using the indirect imaging code built on the basis of the developed approach to reconstruct an input image from a series of generated noisy spectra. The problem of image reconstruction from two-aperture interferometry data has been particularly addressed since it represents the case of most presently operating interferometers.

Key words: techniques: image processing - techniques: interferometric - techniques: spectroscopic - stars: activity - stars: imaging - stars: oscillations


1 Introduction

1.1 Differential interferometry method

Over the last three decades, significant progress has been achieved in the development of interferometry at optical wavelengths. The introduction of single aperture speckle interferometry (Labeyrie 1970) and long baseline interferometry (Labeyrie 1975) demonstrated the potential of the technique which uses the instantaneous fine structure in object images ("speckles'') to attain angular resolution close to the diffraction limit of the instrument, despite atmospheric turbulence.

Beckers (1982) described a technique of narrow spectral band speckle interferometry which derives spatial information with a resolution higher than the instrument diffraction limit, using the variations in the spectrum across the astronomical object under study. This technique, called "Differential Speckle Interferometry'' (DSI), measures the relative shift of an object in different spectral bands, and its application is based on the assumption that the shift between two speckle images can be related to the spatial structures of the object under study. Although the phase of the spatial coherence function is corrupted by instrumental and atmospheric phase errors, different methods of data processing can be used to overcome this difficulty. If the two sets of speckle interferograms of the same object are recorded simultaneously at two different wavelengths, but close enough for the Point Spread Function to be the same for both channels, then their ensemble average cross-spectrum provides the information about the relative shift as a fraction of the object size.

Aime et al. (1984), estimated the sensitivity of the method using a cross-spectrum analysis technique which makes it possible to measure sub-milliarcsecond shifts between two speckle patterns at two close wavelengths, demonstrating the feasibility of the DSI technique. The theoretical estimation of expected signal-to-noise ratios in differential speckle interferometry (Petrov et al. 1986; Chelli 1989) demonstrated the practical applicability of the technique to a wide number of sources. Chelli & Petrov (1995b) described the application of the DSI method to the measurement of angular diameters, rotational velocities and position angles of the rotation axes of rotating stars and to the measurement of the angular vectorial separation and radial velocity differences in binary systems. They also presented a method to compute the uncertainties in these parameters using a numerical computation of the limiting performances of DSI for these applications. Chelli & Petrov (1995a) presented a formalism for the estimation of the parameters, modeling the brightness and velocity fields of unresolved or partially resolved astronomical candidates using the average cross spectrum of speckle patterns recorded in different spectral channels to define the seeing independent estimator. They applied it to unresolved sources, when all the information is contained in the object spectrum and in the angular vectorial function  $\vec{\varepsilon}({\lambda})$ representing the variation of the object photometric barycenter (photocenter) with wavelength, giving a simple formal expression of the parameters errors.

1.2 Doppler imaging method

In parallel to the development of interferometric techniques, the method of Doppler Imaging allowed access to the spatial structure of non-resolved stars through observation and interpretation of temporal spectroscopic variability along the stellar rotational phase. The fact that the wavelength axis of a star is spatially resolved, due to the rotational Doppler shift, is the basis for the method invented by Deutsch (1970), who used the variation of line equivalent widths to adjust parameters describing the development in spherical harmonics of the stellar surface inhomogeneities. Methods making full use of the profile have been further developed (e.g. Khokhlova & Ryabchikova 1975; Vogt & Penrod 1983). However, these methods suffered from a rather arbitrary procedure and a considerable uncertainty in deciding when a correct solution has been obtained.

The inverse problem in Doppler Imaging is fundamentally ill-conditioned; the solution is very unstable against small perturbations in the data and can be stabilized only by progressively adding more weight to the regularization constraint (Jankov & Foing 1987). In fact, the current status of Doppler Imaging has been achieved by using sophisticated regularization methods able to stabilize the inversion.

Different formulations of the Doppler imaging method have been proposed or applied to various observations: Goncharskij et al. (1982) formulated the problem of finding local abundances for Ap stars, in terms of an integral equation; a Lagrange multiplier method is used to minimize the error between calculated and observed profiles, with a stabilization Tikhonov functional.

Vogt et al. (1987) described the method of Doppler imaging for spotted late-type stars, expressing the relation between local surface intensities and the observed spectral profile in a matrix form, by approximating the projection matrix as the marginal response of data pixel to changes in image pixel. Assuming the shape independence of the local line profile, Jankov (1987) gave a formulation for the indirect imaging method, in terms of a matrix formalism, treating explicitly the problem of nonlinearity of the image data transformation due to variable continuum flux level of spotted late-type stars. Both methods used the Maximum Entropy approach in order to stabilize the inversion.

However, in the case of an equator-on star the solution is intrinsically non-unique and only additional information (for example, eclipse modulation of spectral lines) can help to resolve the north-south ambiguity. The problem of the non-uniqueness of the solution in that special case was particularly addressed by Petrov (1988) who showed that the knowledge of $\vec{\varepsilon}({\lambda})$ constrains the model of the star, securing the uniqueness of the reconstructed image.

Jankov et al. (1992) discussed the common situation of regularization of the inverse problem for the interferometric and tomographic methods, while the full mathematical formulation of the problem of stellar Tomographic Imaging from projections provided by spectroscopic and photometric data is given by Jankov & Foing (1992). Although the method was developed for imaging of late-type stars, the basic principles are generally applicable and their formalism is used throughout the present paper. Jankov et al. (1999) presented an approach for imaging of non-radial stellar pulsations that is fully applicable to one of the most promising applications of DSI: the study of the surface of rapidly rotating stars with strong non-radial oscillations (Vakili 1990).

1.3 Stellar atmospheric structures resolved by differential interferometry

High angular resolution differential interferometry imaging has already provided particularly important new information on the atmospheric structure, wind and non-radial pulsations in extended envelopes of Be stars. Vakili et al. (1994) discussed the optical resolution of Be star envelopes in the context of data from the Grand Interféromètre à 2 Télescopes (GI2T) on Plateau de Calern (Mourard et al. 1994).

In recent years, one of the most studied stars with long baseline interferometry has been the Be star $\gamma$ Cassiopeiae. The H$\alpha$emitting region was resolved for the first time by Thom et al. (1986), using optical interferometric measurements from the Interféromètre à 2 Télescopes (I2T), while details in the rotating envelope of the star have been resolved by Mourard et al. (1989). The model presented by Stee et al. (1995), based on spectroscopic and interferometric data collected with the GI2T, provided the detailed physics and geometry of the star's envelope, constraining both the density and velocity relationships present in the equatorial plane. Stee et al. (1998) concluded that the size of the emitting region (as observed in He I $\lambda 6678$ Å, continuum $\lambda 4800$ Å, continuum $\lambda 6500$ Å, H$\beta$ and H$\alpha$) increases from 2.3 to 17 stellar radii. Using differential interferometry observations, Sanchez et al. (1997) reported the non-axisymmetric envelope of $\gamma$ Cas. Berio et al. (1999) revealed azimuthally asymmetric variations of the H$\alpha$ profile due to a prograde one-armed oscillation precessing in the equatorial disk of the star.

A few other stars have been observed using spectrally resolved interferometry: Harmanec et al. (1996) investigated jet-like structures in the Be star $\beta$ Lyrae from an extensive collection of interferometric, spectroscopic and photometric observations. Using the high spatial resolution provided by the GI2T, Vakili et al. (1997) detected subtle structures in the wind of P Cygni by analyzing spectrally resolved fringes in ${\rm H}_\alpha$ and He I $\lambda 6678$ Å. Vakili et al. (1998) reported the first interferometric detection of a prograde one-armed oscillation in the equatorial disk of $\zeta$ Tauri.

2 Basic equations

2.1 Differential Interferometry

Let $\tilde{I'}(u,v,\lambda)$ denote the spatial Fourier transform (complex visibility function) of sky brightness  $I'(s,p,\lambda)$ at the wavelength $\lambda$, where the Fourier frequency pair u,v is the conjugate of the spatial Cartesian coordinates s,p, respectively. Applying the derivative and definite integral theorem to the Maclaurin expansion of  $\tilde{I'}(u,v,\lambda)$ one obtains (see Appendix):
 
$\displaystyle \tilde{I'}(u,v,\lambda)$ =$\displaystyle \int\int_\Omega {I'}(s,p,\lambda) {\rm d}s{\rm d}p
- i 2 \pi u \i...
... d}s{\rm d}p - i 2 \pi v \int\int_\Omega p {I'}(s,p,\lambda) {\rm d}s{\rm d}p ,$ (1)

where the domain of integration in Eq. (7) has been replaced with the visible stellar surface $\Omega$, and we have retained only the zero and first order terms. This is justified by the fact that, when applied to an unresolved source, the integrands of Eq. (7) are determined by the powers of us and vp and since $us \ll 1$ and $vp \ll 1$, they are vanishing in high-order terms of this equation. For the same reason the imaginary ($\Im$) part of the Eq. (1) can be considered as much smaller than the corresponding real part ($\Re$), consequently the argument of the complex visibility  $\tilde{I'}(u,v,\lambda)$: ${\rm Arg} [\tilde{I'}(u,v,\lambda)]
=\arctan [\Im \tilde{I'}/\Re \tilde{I'}]
\approx [\Im \tilde{I'}/\Re \tilde{I'}]
$can be expressed as the sum of two terms:

\begin{displaymath}-{\rm Arg} [\tilde{I'}(u,v,\lambda)]
= 2 \pi u \varepsilon^{s}(\lambda) + 2 \pi v \varepsilon^{p}(\lambda) ,
\end{displaymath}

where:

 \begin{displaymath}\varepsilon^{\rm s}(\lambda) =
{{\int\int_\Omega s {I'}(s,p,\...
...p}
\over{\int\int_\Omega {I'}(s,p,\lambda) {\rm d}s {\rm d}p}}
\end{displaymath} (2)

and

 \begin{displaymath}\varepsilon^{\rm p}(\lambda) =
{{\int\int_\Omega p {I'}(s,p,\...
...p}
\over{\int\int_\Omega {I'}(s,p,\lambda) {\rm d}s {\rm d}p}}
\end{displaymath} (3)

represent two orthogonal shifts in the space domain. The vectorial function $\vec{\varepsilon}(\lambda),$ with the projections $\varepsilon^{s}(\lambda), \varepsilon^{p}(\lambda)$represents the angular shift of the stellar photocenter at different wavelengths and can be evaluated measuring the relative phase of the interferometric signal along a spectral line with respect to the continuum. It provides the first order moment of the spatial brightness distribution, in addition to the zero order moment spectroscopic information.

2.2 Rotationally broadened stellar spectra

Let us consider a star rotating with projected rotational velocity  $V_{\rm e} \sin i$, observed in a spectral line centered at $\lambda_0$. The intensity in the line $I'(s,p,\lambda)$is related to the normalized spectrum $H(M,\lambda)$, as it would be observed in the spatially resolved nonrotating star at a point M(s,p), in the stellar atmosphere by the relation:

 \begin{displaymath}I'(M,\lambda-\lambda_0-\lambda_r)=H(M,\lambda-\lambda_0-\lambda_r) I'_{\rm c}(M,\lambda_0)
\end{displaymath} (4)

with $\lambda_r= y {\lambda_0} (V_{\rm e} \sin i) / c$. The observer's coordinate system is defined by the axis z, chosen to be oriented in the direction of the projection of the stellar rotational axis on the tangent plane and by the axis y orthogonal to it in the same plane (see Fig. 1). From this figure the longitude L and latitude $B=\pi / 2 - \theta$ can be associated with y,z as well as with the angle $\Theta $ between the line of sight x and the normal to the surface at point M:
 
                                            x =$\displaystyle \cos \Theta = \sin B \cos i + \cos B \sin i \cos \phi,$  
y =$\displaystyle \cos B \sin \phi,$  
                              z =$\displaystyle \sin B \sin i - \cos B \cos \phi \cos i,$ (5)

where we have assumed spherical geometry (radius equal to unity), $\phi=L+2\pi\Phi$ is the azimuth related to the longitude L and the particular rotational phase $\Phi $. The intensity in the continuum as observed on the apparent disc of the star $I'_{\rm c}(M,\lambda_0)$ is related to the surface intensity $I_{\rm c}(M,\lambda_0)$ by

 \begin{displaymath}I'_{\rm c}(M,\lambda_0) =
I_{\rm c}(M,\lambda_0)
D(M, \lambda_0),
\end{displaymath} (6)

where the limb darkening law for the continuum intensity distribution on the apparent disc of the star $D(M, \lambda_0)$is usually described by the linear combination of the limb darkening coefficient $\varepsilon_{\lambda_0}$ and the angle $\Theta $:

\begin{displaymath}D(M, \lambda_0) = (1-\varepsilon_{\lambda_0} + \varepsilon_{\lambda_0} \cos \Theta) .
\end{displaymath}


  \begin{figure}
\par\includegraphics[width=6.1cm,clip]{spn.ps}\end{figure} Figure 1: Connection between stellar and observer's coordinate systems. The position of a point M on the stellar surface is determined by its longitude L and colatitude $\theta $. The position of the zeroth meridian is determined by rotational phase $\Phi $. The observer's coordinate system in the tangent plane is defined by the axis z chosen to be oriented in the direction of the projection of stellar rotational axis on to the tangent plane and orthogonal axis y. Axis x, oriented toward the observer, is tilted to the stellar rotational axis by the angle i and to the surface normal at M by the angle $\Theta $. The arc i as well as the arc $\Theta $ are projected as straight lines. With respect to y,z the instrumental coordinate system s, p is rotated in the tangent plane by the angle $\eta $. The shadowed region represents the constant radial velocity strip (its width $\Delta y$ is determined by spectral resolution) corresponding to Doppler shift $\lambda _r$.
Open with DEXTER

In the spectroscopic case we cannot observe directly the specific intensities at the surface, but only their integrals over the apparent stellar disc, i.e. the fluxes, in the line $F(\lambda-\lambda_0)$ and in the continuum $F_{\rm c}(\lambda_0)$. Consequently the observed normalized flux spectrum in the spectral line is given by the relation

\begin{displaymath}R(\lambda - \lambda_0)
= {F(\lambda-\lambda_0)\over F_{\rm c...
...da_r)
\cos \Theta {\rm d}S\over
F_{\rm c}(\lambda_0)
} \cdot
\end{displaymath}

Using Eqs. (4) and (6) and introducing

 \begin{displaymath}A(M,\lambda_0) =
D(M, \lambda_0)
\cos \Theta ,
\end{displaymath} (7)

the observed normalized profile can be written in the form

 \begin{displaymath}R(\lambda - \lambda_0)
={\int\int_\Omega
H(M,\lambda - \lambd...
...bda_0)
A(M,\lambda_0) {\rm d}S\over
F_{\rm c}(\lambda_0)
} ,
\end{displaymath} (8)

where the observed continuum flux is

 \begin{displaymath}F_{\rm c}(\lambda_0)=\int\int_\Omega
I_{\rm c}(M,\lambda_0) A(M,\lambda_0) {\rm d} S.
\end{displaymath} (9)

Dividing both the numerator and denominator of the Eqs. (2) and (3) by the continuum flux $F_{\rm c}(\lambda_0)$and substituting Eqs. (4) and (6), with ${\rm d}s {\rm d}p = \cos \Theta \,{\rm d}S,$ one obtains (using Eq. (8))

 \begin{displaymath}E^s(\lambda) =
{\int\int_\Omega
s H(M,\lambda - \lambda_0 - ...
...ambda_0)
A(M,\lambda_0) {\rm d}S\over
F_{\rm c}(\lambda_0)
},
\end{displaymath} (10)

and

 \begin{displaymath}E^p(\lambda) =
{\int\int_\Omega
p H(M,\lambda - \lambda_0 - ...
...ambda_0)
A(M,\lambda_0) {\rm d}S\over
F_{\rm c}(\lambda_0)
},
\end{displaymath} (11)

where $E^s(\lambda)=R(\lambda)\varepsilon^s(\lambda)$ and $E^p(\lambda)=R(\lambda)\varepsilon^p(\lambda)$, are the spectra of normalized first order moments with respect to the coordinate pair $s,\,p$ which can be associated to an arbitrary instrumental coordinate system. For example, it can be related to the slit and perpendicular-to-slit directions in the case of one aperture speckle interferometry or to two orthogonal baselines of a multi aperture interferometer. But in our case, it is more appropriate to evaluate the moments in the coordinate system related to stellar rotation:

\begin{displaymath}E^{\perp}(\lambda)=R(\lambda)\varepsilon^{\perp}(\lambda), \q...
...parallel}(\lambda)=R(\lambda)\varepsilon^{\parallel}(\lambda),
\end{displaymath}

where $\varepsilon^{\perp}(\lambda)$ is the spectrum of photocenter shift orthogonal to z and $\varepsilon^{\parallel}(\lambda)$ is the spectrum of photocenter shift parallel to z. The coordinates y,z can be related to s and p in terms of a transformation under rotation of the coordinate system:

\begin{eqnarray*}s &=& y \sin \eta + z \cos \eta, \\
p &=& - y \cos \eta + z \sin \eta,
\end{eqnarray*}


where the time evolution of the angle $\eta $ between s and the projection of the stellar rotational axis zis provided by book astrometry, assuming that it is known for a reference position $\eta_0$. Substituting s and p into the Eqs. (10) and (11) it follows:
 
$\displaystyle E^{\rm s}(\lambda)$ =$\displaystyle E^{\perp}(\lambda) \sin \eta + E^{\parallel}(\lambda) \cos \eta,$  
$\displaystyle E^{\rm p}(\lambda)$ =$\displaystyle - E^{\perp}(\lambda) \cos \eta + E^{\parallel}(\lambda) \sin \eta,$ (12)

with

 \begin{displaymath}E^{\perp}(\lambda) =
{\int\int_\Omega
H(M,\lambda - \lambda_0...
...mbda_0)
B(M,\lambda_0) {\rm d}S\over
F_{\rm c}(\lambda_0)
},
\end{displaymath} (13)


 \begin{displaymath}E^{\parallel}(\lambda) =
{\int\int_\Omega
H(M,\lambda - \lamb...
..._0)
C(M,\lambda_0) {\rm d}S\over
F_{\rm c}(\lambda_0)
}\cdot
\end{displaymath} (14)

Here we have introduced

 \begin{displaymath}B(M,\lambda_0) =
y D(M, \lambda_0)
\cos \Theta =
\end{displaymath} \begin{displaymath}(\cos B \sin \phi)
D(M, \lambda_0)
\cos \Theta ,
\end{displaymath} (15)
and

 \begin{displaymath}C(M,\lambda_0)=
z D(M, \lambda_0)
\cos \Theta =
\end{displaymath} \begin{displaymath}(\sin B \sin i - \cos B \cos \phi \cos i )
D(M, \lambda_0)
\cos \Theta .
\end{displaymath} (16)

Note that the implementation to stars with non-spherical surfaces will affect not only Eq. (5), but also the limb darkening law $D(M, \lambda_0)$ in Eqs. (7), (15) and (16). In order to obtain the appropriate model this should be taken into account for rapidly rotating and close binary stars by using the Roche geometry and gravity darkening law (e.g. Townsend 1997).

3 Application to non-radial stellar pulsations

In order to illustrate the application of the preceeding basic equations, we consider the complex brightness distribution on the surface of a non-radially pulsating star. The non-radial pulsator model describes the surface intensity distribution and velocity field in terms of the associated Legendre functions $P_{l}^{m}(\theta , \phi_{\rm R})$ and time t as:

 \begin{displaymath}I_{\rm c}(M,\lambda_0,t) = I_0(\lambda_0) +
A_{{\rm I}}
P_{ l}^{ m}(\cos \theta) {\rm e}^{i({ m}\phi_{\rm R} + \omega t)},
\end{displaymath} (17)

where $I_0(\lambda_0)$ is the immaculate intensity of the star, $A_{{\rm I}}$ is the oscillation amplitude, $\omega $ is the angular time frequency in the observer's frame, the quantum numbers l, m denote the pulsation degree and azimuthal order, respectively, and  $\phi_{\rm R}$ is the azimuth related to the rotational frequency of the star  $\Omega_{\rm R}$. In order to establish the stationarity (which is necessary for Doppler Imaging), the rotational phase should be related to the frequency  $\Omega_{\rm W}$ defined as:

\begin{displaymath}m \Omega_{\rm W} t = m \phi = m \phi_{\rm R} + \omega t ,
\end{displaymath}

where $\phi$ is the corresponding azimuth defined by the same relation.

Figure 2 displays the stellar image $I_{\rm c}(M,\lambda_0)$ corresponding to the surface intensity perturbation due to the non-radial pulsation m=4, l=5 mode in the stationary reference frame. In our example, the artificial star has a rotation axis tilted by $i=45^\circ $ to the line of sight and rotates with a projected equatorial velocity $V_{\rm e} \sin i$ of 150 kms-1. In a rapidly rotating star the pulsation velocity field acts as a small perturbation to the dominant rotational velocity field. This perturbation is mapped onto a wavelength corresponding to the rotationally induced Doppler shift. Using Eqs. (8), (13) and (14) we calculated the profiles of spectrum normalized flux $R(\lambda )$, shifts orthogonal to the rotational axis $\varepsilon^{\perp}(\lambda)$ and parallel to the rotation $\varepsilon^{\parallel}(\lambda)$, corresponding to the surface intensity $I_{\rm c}(M,\lambda_0)$ distribution as shown in Fig. 2. The synthetic intensity spectra and photocenter shifts are calculated using the He I $\lambda 6678$ Å intrinsic profile corresponding to a typical OB star, with Hubeny & Lanz's (1995) spectral synthesis. The oscillation amplitude $A_{{\rm I}}$ was set to reproduce the strength of bumps ( $\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil
$\displaystyle ...1%) observed across this line in some Be stars (e.g. Jankov et al. 2000), and plots associated with the rotational phase $\Phi = 0$ are presented on top of Figs. 3, 4 and 5. The synthetic data set is then computed, consisting of normalized flux and photocenter shift profiles evenly spread throughout the rotational cycle. Artificial zero-average Gaussian noise ( $1/\sigma =1000$), is added to produce the dynamic residual spectrum (the immaculate star line profiles are subtracted from perturbed star line profiles) of normalized intensity, orthogonal and parallel photocenter shifts as shown in bottom of Figs. 34 and 5 respectively.

  \begin{figure}
\par\hspace*{5mm}\includegraphics[width=5.7cm,clip]{et45.ps}\par\includegraphics[angle=-90,width=6.5cm,clip]{map_45.ps}\end{figure} Figure 2: Stellar surface brightness perturbation due to the non-radial pulsation m=4, l=5 mode on a star tilted at $i=45^\circ $. Top: Spherical projection. Bottom: Mercator projection of the visible surface.
Open with DEXTER


  \begin{figure}
\par\includegraphics[angle=-90,width=8.2cm,clip]{fs.ps}\par\includegraphics[angle=-90,width=8.3cm,clip]{d45s.ps}\end{figure} Figure 3: Synthetic normalized flux in the He I spectral line. Top: the spectrum calculated for $\Phi = 0$, using the model from Fig. 2 (solid line). The dashed line represents the spectrum corresponding to the homogeneous stellar surface. Bottom: dynamic spectrum of residuals (immaculate star profiles are subtracted from perturbed profiles) evenly spread throughout the rotational cycle. Artificial Gaussian noise ( $1/\sigma =1000$) is added.
Open with DEXTER


  \begin{figure}
\par\includegraphics[angle=-90,width=8.2cm,clip]{phy.ps}\par\includegraphics[angle=-90,width=8.3cm,clip]{d45y.ps}\end{figure} Figure 4: Same as in Fig. 3, but for photocenter shifts orthogonal to rotation.
Open with DEXTER


  \begin{figure}
\par\includegraphics[angle=-90,width=8.2cm,clip]{phz.ps}\par\includegraphics[angle=-90,width=8.3cm,clip]{d45z.ps}\end{figure} Figure 5: Same as in Fig. 3, but for photocenter shifts parallel to rotation. The shifts corresponding to the homogeneous star are equal to zero across the entire profile.
Open with DEXTER

4 Matrix formulation of image-data transformation

Let us consider the stellar surface divided by $N_{\rm L}$ meridians and $N_{\rm B}$ parallels. Let the $N_{\rm L}$ meridians and $N_{\rm B}$ parallels determine the image resolution pixels. We number the resolution pixels from 1 to J, each pixel representing a part of the stellar surface with constant intensity Ij where

\begin{displaymath}j = m + N_{\rm B}(l-1) \qquad \quad
m=\overline{1,N_{\rm B}}\qquad l=\overline{1,N_{\rm L}}
\end{displaymath}

are the coordinates of the resolution pixels. In order to calculate the projection with a given precision (determined by the signal-to-noise ratio of the observations) we divide each resolution pixel by $N_{\rm P}$ projection pixels (for more details see Jankov & Foing 1992). The projection pixel coordinate frame is chosen so that elementary surface of the image (surface of the projection pixel) is equal to the surface element corresponding to adopted geometry.

With regard to this discretization, we define the spectrum of the normalized flux (Eq. (8)) measured in N discrete wavelengths across the spectral line:

\begin{displaymath}R_n={{\sum_{r=1}^N \sum_{j=1}^J \sum_{k=1}^{N_{\rm P}}
H_{n-r...
...{rk}} }\over {F_{\rm c}(\lambda_0)}},
\qquad n=\overline{1,N},
\end{displaymath}

where the subscripts r and n represent the wavelength (measurement in the nth position of the detector pixel), the summation  $\sum_{r=1}^N$ concerns the broadening with the  $H(M,\lambda)$ profile denoted by Hn-r,k, the summation  $\sum_{j=1}^{J}$concerns all the resolution pixels, and the summation $\sum_{k=1}^{N_{\rm P}}$concerns all the projection pixels inside the resolution pixel j.

The profile Hn-r,k describes the dependence of the intensity of the local spectrum as a function of wavelength (index r corresponding to the strip $y_r,y_r + \Delta y$) and of position on the stellar surface (index k corresponding to the projection pixel inside the resolution element j). The profiles  $H(M,\lambda)$ can be obtained from observations of non-rotating template stars (Jankov & Foing 1992); otherwise they can be calculated using spectral synthesis, but in that case they should be convolved with the instrumental function. This function is independent of M, in contrast with the broadening induced by stellar rotation in conjuction with finite exposure time. The latter broadening depends on the position at the stellar surface M and should be convolved with the local profile $H(M,\lambda)$ in each resolution pixel j. The quantity Ark is defined by the Eq. (7), and the domain of integration associated with the equal radial velocity strip r is introduced by:

\begin{displaymath}b_{rk} = \ \ \ \cases{1 & $y \in
(y_r,y_r+\Delta y)$\space \cr
0 & $y \not\in (y_r,y_r+\Delta y) $\space \cr}.
\end{displaymath}

Here brk is unity if the particular projection pixel k is covered by the strip r, and zero otherwise. Note that this definition of the surface of the projection pixel brk imposes the dimensional equivalence between intensities and fluxes in our presentation. Since we assume invariability of the physical conditions inside the resolution pixel j:

\begin{displaymath}I_1 = I_2 =...= I_k =...= I_{N_{\rm P}} = I_j
\end{displaymath}

and

\begin{displaymath}H_{n-r,1} = H_{n-r,2} =...= H_{n-r,k} =...= H_{n-r,N_{\rm P}} = H_{n-r,j}
\end{displaymath}

it follows that
 
$\displaystyle R_n={{\sum_{j=1}^J I_j \sum_{r=1}^N H_{n-r,j}\sum_{k=1}^{N_{\rm P...
...m_{r=1}^N H_{n-r,j}\sum_{k=1}^{N_{\rm P}}
{b_{rk}A_{rk}} }\over {^{(q)}C}}\cdot$     (18)

For the qth spectrum this equation can be written in matrix notation:

\begin{displaymath}^{(q)}C\cdot^{(q)}R_N={{^{(q)}S_{NJ} \star X_J}},
\end{displaymath}

where $X_j=I_j / I_0(\lambda_0)$ and $^{(q)}C = F_{\rm c}(\lambda_0) / I_0(\lambda_0)$; $I_0(\lambda_0)$ being the immaculate intensity of the star. Here we introduced the matrix (q)SNJ with the components

\begin{displaymath}S_{nj}=\sum_{r=1}^N H_{n-r,j} D_{rj}
\end{displaymath}

where Drj is defined by:

\begin{displaymath}D_{rj} = \sum_{k=1}^{N_{\rm P}} b_{rk}A_{rk},
\end{displaymath}

brk and Ark being dependent on j.

Similarly we can deduce from Eqs. (13) and (14) that

 
$\displaystyle ^{(q)}C\cdot^{(q)}E^{\perp}_N$ = $\displaystyle {{^{(q)}O_{NJ} \star X_J}}$  
$\displaystyle ^{(q)}C\cdot^{(q)}E^{\parallel}_N$ = $\displaystyle {{^{(q)}P_{NJ} \star X_J}},$ (19)

where

\begin{displaymath}O_{nj}=\sum_{r=1}^N H_{n-r,j} F_{rj}, \qquad
F_{rj} = \sum_{k=1}^{N_{\rm P}} b_{rk}B_{rk}
\end{displaymath}


\begin{displaymath}P_{nj}=\sum_{r=1}^N H_{n-r,j} G_{rj}, \qquad
G_{rj} = \sum_{k=1}^{N_{\rm P}} b_{rk}C_{rk}
\end{displaymath}

with Brk and Crk defined by Eqs. (15) and (16), respectively. The quantities $E^{\perp}$ and  $E^{\parallel}$ are not directly observable, but can be related to the observables $E^{\rm s}$ and $E^{\rm p}$ through Eq. (12).

The problem of non-linearity due to the variable continuum level (q)C has been explicitly treated by Jankov & Foing (1992). Simultaneous photometry can be used to measure (q)C, but if it is not available, the transformation matrices (q)SNJ, (q)ONJ and (q)PNJ should be replaced by matrices (q)SNJ - (q)VNJ, (q)ONJ - (q)YNJ, (q)PNJ - (q)ZNJ, with the components: $V_{nj}=R_n \sum_{r=1}^N D_{rj}$, $Y_{nj}=E^{\perp}_n \sum_{r=1}^N D_{rj}$ and $Z_{nj}=E^{\parallel}_n \sum_{r=1}^N D_{rj}$. In that case the left-hand side of the equations describing the image-data transformation should be set to zero. Nevertheless, the image-data transformation remains non-linear when the local line profile dependence on temperature is taken into account. In that case the transformation matrix depends on XJ and should be recalculated in each iteration, with respect to the temperature distribution over the stellar surface from the previous iteration.

The vectors (q)RN, $^{(q)}{E}^{\perp}_N$, and $^{(q)}{E}^{\parallel}_N$, $(q=\overline{1,Q}),$have the components Rn, $E^{\perp}_n$, $E^{\parallel}_n$, $(n=\overline{1,N}),$ representing the measured intensity in the nth position of the detector pixel for the observed frame q. In order to obtain the image-data transformation in the form ${Y}_I=R_{IJ} \star X_J$, with Q observations consisting of normalized flux spectra:

(1)RN,(2)RN,...,(q)RN,...,(Q)RN

first-order moment orthogonal to rotation:

\begin{displaymath}^{(1)}{E}^{\perp}_N,^{(2)}{E}^{\perp}_N,...,^{(q)}{E}^{\perp}_N,...,^{(Q)}{E}^{\perp}_N
\end{displaymath}

and first-order moment parallel to rotation:

\begin{displaymath}^{(1)}{E}^{\parallel}_N,^{(2)}{E}^{\parallel}_N,...,^{(q)}{E}^{\parallel}_N,...,^{(Q)}{E}^{\parallel}_N,
\end{displaymath}

we arrange these vectors (and corresponding transformation matrices ) into the data vector YI and block matrix RIJ with components Yi and Rij, where $i=n+M \cdot (q-1)$, with M equal to Q or 2Q, depending on possible combinations. By reason of redundancy, we do not use together the normalized flux spectra (q)RN and moments orthogonal to rotation $^{(q)}{E}^{\perp}_N$.

5 Numerical experiments of image reconstruction

5.1 Multi-baseline interferometric observations

Here we consider the case when the quantities $E^{\perp}$ and  $E^{\parallel}$ can be related to the observables Es and Ep through Eq. (12). This is likely to be met in practice for measurements obtained with multi-baseline interferometers as well as for single aperture speckle differential interferometry where the source intensity distribution can projected on different directions (e.g. Lagarde 1994).

In principle, an image can be reconstructed using only the stellar flux spectra, or any of the photocenter shift projections. For example, for the photocenter projection orthogonal to rotation, the corresponding data vector YI and the projection matrix RIJ have $I=Q\cdot N$ components and the image-data transformation can be written in the form of matrix product ${Y}_I=R_{IJ} \star X_J$:

\begin{displaymath}{{Y}_I}
= {C}_Q \cdot {\pmatrix{ ^{(1)}{E}^{\perp}_{N}\cr
^...
...
\vdots \cr
^{(q)}O_{NJ}\cr
\vdots \cr
^{(Q)}O_{NJ}\cr}},
\end{displaymath}

where the row vector

\begin{displaymath}{C}_Q = {\pmatrix{^{(1)}C, ^{(2)}C,\ldots, ^{(q)}C, \ldots, ^{(Q)}C}}
\end{displaymath}

represents the normalized continuum flux corresponding to each observation. Due to the cancellation effect, the medium-degree NRP mode (as presented in Fig. 2) shows a very low-amplitude light variability. Consequently in the following the entries of the vector CQ were considered to be constant without affecting the image reconstruction.

We have applied our indirect imaging code to reconstruct the input stellar image using 30 spectra equally spaced in time. These spectra have been extracted from the dynamic residual spectrum of photocenter shifts orthogonal to rotation as presented in Fig. 4, with a wavelength step corresponding to spectral resolution $\lambda /\Delta \lambda =15\,000$. Complete sampling of the stellar surface corresponding to this spectral resolution requires our relatively large number of spectra. Detailed discussion of the appropriate observing strategy, in order to obtain optimal image resolution, is given by Jankov & Foing (1992). However, the maps reconstructed with half of this number have approximately the same quality of reconstruction but lower image resolution. Note that in the limit, when spectral resolution and number of spectra approaches infinity (as presented in Figs. 3, 4 and 5), and in the absence of noise, the dynamic difference spectrum represents the projection from image space into a data space. In practice, the dynamic difference spectrum, being an approximation of the projection from the image space into a data space, can be used as an input to perform the deprojection, but the image reconstruction procedure should be regularized. The reconstructed image shown in Fig. 6 (top), as well as all reconstructions throughout this paper, has been obtained using Maximum Entropy regularization, where the solution is chosen by minimizing:

\begin{displaymath}-\Lambda \sum_{j=1}^J X_j {\rm ln}(X_j/^0X_j) +
\sum_{i=1}^I {\Bigl ( (Y_i-\sum_{j=1}^J R_{ij} X_j)/
{\sigma_i}^2\Bigr )}^2,
\end{displaymath}

where the standard deviation $\sigma_i$ represents the measurement error and $\Lambda$can be determined subject to the data statistics constraint (right-hand term of the expression): $\chi^2 = \chi^2_0$, where $\chi^2_0$ is chosen to construct a confidence domain enclosing (say) $99\%$of the probability that the true object contains more information than the maximum entropy estimate. In common with the regularization by that structure condition, the maximum entropy approach imposes a positivity constraint ( $X_j~>~0, \quad j=\overline{1,J}$) on the unknown image function, if the background level $^0X_j\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil
$\displaysty...
...offinterlineskip\halign{\hfil$\scriptscriptstyle ... is used. The benefits of the positivity constraint are emphasized in the prior knowledge that only dark spots (or bright plages) are present at the stellar surface. Unfortunately it cannot be used in the case of non-radial stellar pulsations since both positive and negative disturbances are affecting the surface. However, comparing with Fig. 2, one notices that the global structure of the input image is reasonably well reconstructed (spurious features being due to low effective signal-to-noise ratio of the data), but the reconstruction shows an appreciable latitude loss of contrast, the quality of the reconstruction being latitude dependent.

Similar image-data transformations can be formulated for the photocenter projection parallel to rotation, for which the Maximum Entropy reconstruction is shown in Fig. 6 (bottom) and for the normalized flux spectra, for which the reconstruction is shown in Fig. 7 (top).

The simulations presented here indicate that stellar surface imaging is feasible not only from flux spectra but also from any component of photocenter shift, but with different reliability of latitude information. In order to achieve better image reconstructions, the image data transformation can be defined to incorporate both photocenter projections or flux spectra and photocenter projection parallel to rotation. The corresponding data vector YI and the projection matrix RIJ have $I=2Q\cdot N$ components and, for the later case, the image-data transformation can be written in the form ${Y}_I=R_{IJ} \star X_J,$where

 \begin{displaymath}{{Y}_I}
= {C}_{2Q} \cdot {\pmatrix{ ^{(1)}{E}^{\parallel}_{N...
...r
\vdots \cr
^{(q)}S_{NJ}\cr
\vdots \cr
^{(Q)}S_{NJ}\cr}},
\end{displaymath} (20)

with

\begin{displaymath}{C}_{2Q} = {\pmatrix{^{(1)}C, \ldots, ^{(q)}C, \ldots, ^{(Q)}C,
^{(1)}C, \ldots, ^{(q)}C, \ldots, ^{(Q)}C}}.
\end{displaymath}

The derived image-data transformations have been employed to reconstruct the input image by inverting a set of time-resolved normalized flux spectral profiles and photocenter shifts parallel to rotation as shown in Figs. 3 and  5. A map of the brightness distributions has been produced by means of Maximum Entropy regularization and the resulting image is shown in Fig. 7 (bottom).
  \begin{figure}
\par\includegraphics[angle=-90,width=8.2cm,clip]{map1_45y.ps}\par\includegraphics[angle=-90,width=8.2cm,clip]{map1_45z.ps}\end{figure} Figure 6: Maximum Entropy reconstructions from two orthogonal photocenter shift spectra. Only 30 equidistantly spaced spectra from Figs. 4 and 5 were used as input, with a wavelength step corresponding to a spectral resolution of $\lambda /\Delta \lambda =15\,000$. Top: Reconstruction from $E^{\perp }(\lambda )$, the component orthogonal to rotation. Bottom: Reconstruction from $E^{\parallel }(\lambda )$, the component parallel to rotation. In both figures, one notices a latitude smearing and the presence of spurious features induced by noise.
Open with DEXTER


  \begin{figure}
\par\includegraphics[angle=-90,width=8.2cm,clip]{map1_45s.ps}\par\includegraphics[angle=-90,width=8.2cm,clip]{map1_45sz.ps}\end{figure} Figure 7: Maximum Entropy reconstructions from normalized flux spectra alone, and from flux spectra together with photocenter shifts parallel to rotation. The sampling of input dynamic spectra (Figs. 3 and 5) was identical to that of the reconstructions in Fig. 6. Top: Maximum Entropy reconstruction from flux spectra  $R(\lambda )$ alone. Note the intrinsic limitations on the reconstruction of features in the hemisphere in which the rotational pole is hidden. Bottom: Reconstruction from flux spectra and the component parallel to rotation $E^{\parallel }(\lambda )$, using the image-data transformation from Eq. (20). Note the significant improvement with respect to previous reconstructions.
Open with DEXTER

A map recovered from this reconstruction should be compared to the original one (Fig. 2) as well as to the maps obtained using single photocenter projections (Fig. 6) or spectra alone (Fig. 7, top). One can notice a significant improvement with respect to previous reconstructions. In particular, the recovery of information about the distribution in latitude is improved.

5.2 Single-baseline interferometer

To illustrate the capabilities of the technique under more challenging conditions we consider another test case, that of a single-baseline interferometer for which only one component (say Es) can be measured with the required spatial resolution. In that case the dispersed interference pattern contains only the information in the direction of axis s representing the projection of the interferometer baseline onto the sky. In order to calculate the input moments  $E^{s}(\lambda)$, we used the dynamic residual spectrum of a corresponding photocenter shift  $\varepsilon^{s}(\lambda)$, over one rotational cycle of the star displayed in Fig. 2. This dynamic spectrum cannot be visually distinguished from the one presented in Fig. 4, because the residual component of $\varepsilon^{\parallel}(\lambda)$ is much smaller than $\varepsilon^{\perp}(\lambda)$.

The corresponding image-data transformation can be obtained using Eqs. (12) and (19). For example, using only the photocenter shift measurements in the direction of the axis s, the image-data transformation can be written in the form

\begin{displaymath}{{{Y}_I}
= {C}_Q \cdot {\pmatrix
{ ^{(1)}{E}^{s}_{N} \cr
^{...
..._{N} \cr
\vdots \cr
^{(Q)}{E}^{s}_{N} \cr
}}}, \nonumber \\
\end{displaymath}


 \begin{displaymath}{R_{IJ} = \pmatrix{
^{(1)}O_{NJ}\sin \eta_1 + ^{(1)}P_{NJ}\...
...\cr
^{(Q)}O_{NJ}\sin \eta_Q + ^{(Q)}P_{NJ}\cos \eta_Q \cr
}}.
\end{displaymath} (21)

Including spectroscopy, the image-data transformation can be arranged as has been done for Eq. (20), but  (q)PNJ, figuring there, should be replaced by

\begin{displaymath}^{(q)}O_{NJ}\sin \eta_q + ^{(q)}P_{NJ}\cos \eta_q .
\end{displaymath}

For $\eta_q =0^\circ,$ this reduces to Eq. (20), but $\eta_q = \pm 90^\circ$implies the redundant combination of normalized flux spectra and moments orthogonal to rotation. A possible way to treat the problem of redundancy is weighting the q-th moment equation with $\cos \eta_q$. With this weighting, for $\eta_q = \pm 90^\circ$, the reconstruction concerns only the normalized flux spectra. However, an optimal strategy of observation would consist in observing the star when the interferometer basis is aligned with the projection of stellar rotational axis ( $\eta \approx 0^\circ$), in order to achieve the best reconstruction, as presented in Fig. 7 (bottom).

The derived image-data transformation (Eq. (21)) has been employed to reconstruct the input image from the calculated spectra  $E^{s}(\lambda)$. In order to show clearly that the reconstruction is not due to the super-synthesis effect (observing the source at different hour angles so that the brightness distribution is projected onto the baseline with different orientations), we assumed the angle $\eta =45^\circ $ to be constant. The resulting reconstruction (Fig. 8 top) shows that reasonably good image reconstructions can be obtained for the case of single baseline interferometric measurements, even under such conditions.

Of course, the visibility data taken throughout a night represent more information and better reconstructions are expected allowing $\eta $ to vary.

5.2.1 Uncertainty in the position angle of rotation axis

A number of parameters enter the image-data transformation in Doppler Imaging. The uncertainty in classical ones, such as inclination, projected rotational velocity, limb darkening law etc. have been discussed in details in literature (e.g. Vogt et al. 1987). Note that the stellar inclination is traditionally the most uncertain value, but it can be determined by using IDI observations. High resolution spectroscopy provides the projected rotational velocities $V_{\rm e} \sin i$, while interferometry (together with precise parallaxes from Hipparcos) provides the absolute stellar radii which give access (supposing the stellar rotational frequency to be known by time-series analysis) to the rotational velocity $V_{\rm e}$.

In order to estimate the influence of the uncertainty of the new parameter $\eta $, the same time series of $E^{s}(\lambda)$ have been used again to perform the reconstruction, but with $15^\circ$ error in the position angle of the stellar rotation axis entering the image-data transformation (Eq. (21)). The reconstructed image is presented in Fig. 8 (bottom). Comparison of two maps shows that the corresponding errors in $\eta $ have not affected significantly the global structure of the image, but only the contrast of reconstructed structures.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.2cm,clip]{map1_45e.ps}\par\includegraphics[angle=-90,width=8.2cm,clip]{map1_45er.ps}\end{figure} Figure 8: Maximum Entropy reconstructions from photocenter shifts measured along a single baseline direction. The reconstructions were performed using the image data transformation (Eq. (21)). Top: Using the correct value for the angle $\eta =45^\circ $ between the projection of interferometer baseline and stellar rotation axis. Bottom: Using a wrong value $\eta =30^\circ $. The global structure of the image is not significantly affected, but only the contrast of reconstructed structures.
Open with DEXTER

6 Discussion

6.1 Latitude smearing and bias

Previous examples show clearly that the interferometric constraint makes the mapping of stellar surfaces much more reliable and informative than Doppler imaging alone. For the interpretation of maps of stellar surfaces obtained from intensity spectra it is important to recognize basic limitations of the technique. Particularly striking is the loss of contrast of features below the equator. For lower signal-to-noise ratios the reconstructions from stellar flux spectra alone show the loss of such structures (Fig. 9 top).

  \begin{figure}
\par\includegraphics[angle=-90,width=8.2cm,clip]{noi1_45s.ps}\par\includegraphics[angle=-90,width=8.2cm,clip]{noi1_45sz.ps}\end{figure} Figure 9: Same as in Fig. 7, but for the spectra affected by noise $1/\sigma =200$, instead of $1/\sigma =1000$. Top: Maximum Entropy reconstruction from flux spectra $R(\lambda )$ alone. Note the vanishing of structures in the hemisphere bellow the stellar equator. Bottom: Reconstruction from flux spectra and the component parallel to rotation $E^{\parallel }(\lambda )$. Despite the spurious features introduced by high noise level, the global structure of the original image (Fig. 2) is conserved.
Open with DEXTER

The reasons for this are clear: the contribution to the line profiles in such areas is significantly reduced due to foreshortening and limb darkening effects (see Fig. 2 top). Moreover, features are visible only briefly, so they contribute to the observed profiles for only a few rotational phases.

This degradation of restored maps is much less present in the reconstructions performed using the moments parallel to rotation in conjunction with normalized flux spectra or moments orthogonal to rotation (Fig. 9 bottom). The corresponding image reconstructions demonstrate that the features in the hemisphere in which the rotational pole is hidden are much better reproduced. When using the moments parallel to rotation, the corresponding stellar regions are reinforced by weighting with z-coordinate (as can be seen from Eqs. (14) and (16)), and consequently better reproduced in the reconstructed map.

6.2 Case of an equator-on star

The particularly challenging case is that of a star with $i \approx 90^\circ$, for which the spectral inversions are not unique. In both cases (normalized flux  $R(\lambda )$ and moment orthogonal to rotation  $E^{\perp }(\lambda )$), when $i=90^\circ$ the entries of the matrix RIJ do not contain the term $\sin B$ (Eqs. (5), (7), (15), (18), (19)), and being dependent only on $\cos B$, produce the situation where the latitude sign is undetermined.

To illustrate this situation we constructed an artificial star non-radially pulsating in m=4, l=5 mode, and using $i=85^\circ $, we computed synthetic dynamic spectra (Fig. 11) as has been done for the previous numerical experiments. Comparing the input image (Fig. 10) and the reconstructed map (Fig. 12 top), it is obvious that the inversion from spectra alone fails to reproduce the correct image.

  \begin{figure}
\par\hspace*{5mm}\includegraphics[width=5.8cm,clip]{et85.ps}\par\includegraphics[angle=-90,width=6.6cm,clip]{map_85.ps}\end{figure} Figure 10: Stellar surface brightness perturbation due to the non-radial pulsation m=4, l=5 mode on a star tilted at $i=85^\circ $. Top: Spherical projection. Bottom: Mercator projection of the visible surface.
Open with DEXTER


  \begin{figure}
\par\includegraphics[angle=-90,width=8.3cm,clip]{d85s.ps}\par\includegraphics[angle=-90,width=8.3cm,clip]{d85z.ps}\end{figure} Figure 11: Dynamic spectrum of residuals (evenly spread throughout the rotational cycle) for a star tilted at $i=85^\circ $. Top: Normalized stellar flux in the He I spectral line with artificial Gaussian noise $1/\sigma =1000$ added. Bottom: Photocenter shift parallel to rotation, in the same line with artificial Gaussian noise $1/\sigma =200$ added.
Open with DEXTER

In fact, the reconstructed mirror symmetry with respect to stellar equator is consistent with zero first-order moments parallel to rotation  $E^{\parallel }(\lambda )$. But the corresponding moments in that case are about one order of magnitude higher than that presented in Fig. 5. Introducing this constraint in the reconstruction, the entries of the matrix RIJ contain also the term $\cos B$ in Eqs. (16), (19) and in that case the latitude sign is determined even for $i=90^\circ$. This is illustrated in Fig. 12 (bottom) where the reconstruction is dramatically improved.

  \begin{figure}
\par\includegraphics[angle=-90,width=8.2cm,clip]{map1_85s.ps}\par\includegraphics[angle=-90,width=8.2cm,clip]{map1_85z.ps}\end{figure} Figure 12: Top: Maximum Entropy reconstruction from normalized flux spectra. One notices a characteristic mirroring due to the north-south ambiguity. Bottom: Maximum Entropy reconstruction from photocenter shifts parallel to rotation. The north-south ambiguity is removed.
Open with DEXTER

6.3 Potential of the method from present and future instruments

The high signal-to-noise ratios ( $1/\sigma =1000$) considered throughout this paper have been used in order to show the performances of the method with respect to data quality from the instruments that will operate in the near future (VLTI). But our examples of image reconstruction from data with lower signal-to noise ratio ( $1/\sigma =200$) show clearly that the technique can be successfully applied even in that case, providing the information that is not accessible from spectroscopy alone (Figs. 9 and 12).

The requirement for a high signal-to-noise ratio could seem somehow unrealistic at present (particularly concerning the photocenter shift measurements), but observational results approaching this limit have already been obtained with the present configuration of the GI2T, using a single spectral line (Vakili et al. 1997). With the simultaneous use of several spectral lines and the least-square deconvolution (Donati et al. 1997), such signal-to-noise ratios could be attained from present instruments modified in order to obtain wide spectral regions (échelle spectrograph) and adaptive optics (Vèrinaud 2000).

In order to become a routinely usable technique, IDI needs to be compared to images of a number of well-known stars for which Doppler-Imaging or Zeeman-Doppler-Imaging techniques have already produced reliable maps. IDI, as well as classical Doppler Imaging, can be applied only to the structures which do not change appreciably during the time needed to cover one rotational cycle of the star. The class of such stars is extended by using multi-site observations and this is the main advantage of existing spectrographs. However, with the advent of future interferometric arrays the multi-site strategy should become available to IDI observations, as well.

7 Conclusions

We have shown the imaging potential of the IDI technique, which combines time-resolved spectroscopy and long baseline interferometry. It improves dramatically the reliability and details of reconstructed stellar surface maps, providing information that cannot be otherwise obtained with each of these techniques taken alone. Success in synthesizing images obtained from this method is achieved by bringing together continuous spatial resolution in the direction of spectral dispersion, provided by Doppler shifts in a rapidly rotating star, and first order moments of brightness distribution provided by interferometric measurements of photometric barycenter shifts.

We have carried out a number of numerical experiments with realistic spectral and/or spatial resolutions expected for operating (the GI2T) or close-to-operating long baseline interferometers (the VLTI, Keck). We conclude that at the reasonable spectral resolution of a few thousands and a desired signal-to-noise ratio of less than a thousand, accurate maps of stellar non-radial pulsations can be obtained by using regularized inversion Maximum Entropy methods. Interestingly enough IDI can solve a number of intrinsically ambiguous cases of stellar configurations like an equator-on rapidly rotating/pulsating star where high resolution and signal-to-noise spectroscopy fails.

In fact, IDI is evidently applicable to other classes of stellar surface structure imaging. For instance, magnetic activity of Ap, Bp or RS CVn stars, surface temperature and/or chemical abundance inhomogeneities.

Acknowledgements
A. D. de S. Jr. acknowledges CAPES (Brazil) (contract BEX 1661/98-1) for financial support. EJP acknowledges support from the FAPESP (Brazil) through grant No. 99/02506 and from the CNRS (France). We thank P. Stee and F. Wilkin for helpfull comments.

Appendix

Expanding a function $\tilde{I'}(u,v,\lambda)$ about the origin u=v=0, we write symbolically the Maclaurin series as:

 \begin{displaymath}\tilde{I'}(u,v,\lambda)=\sum_{n=0}^\infty {{1}\over{n!}}
\lef...
...\right)^n
\left. \tilde{I'}(u,v,\lambda) \right\vert _{u=v=0}.
\end{displaymath} (22)

The derivatives can be expressed using the derivative theorem: If $I'(s,p,\lambda)$ is the inverse Fourier transform of $\tilde{I'}(u,v,\lambda)$ then $-i 2 \pi s {I'}(s,p,\lambda)$ is the inverse Fourier transform of $\partial / \partial u \ \tilde{I'}(u,v,\lambda)$ and $-i 2 \pi p {I'}(s,p,\lambda)$ is the inverse Fourier transform of $\partial / \partial v \ \tilde{I'}(u,v,\lambda)$. Applying the definite integral:

\begin{displaymath}\tilde{I'}(0,0,\lambda) =
\int\int_{-\infty}^\infty {I'}(s,p,\lambda) {\rm d}s {\rm d}p
\end{displaymath}

to the previous theorem one has:

\begin{displaymath}{{\partial}\over{\partial u}} \left. \tilde{I'}(u,v,\lambda) ...
...t\int_{-\infty}^\infty s {I'}(s,p,\lambda) {\rm d}s {\rm d}p ,
\end{displaymath}


\begin{displaymath}{{\partial}\over{\partial v}} \left. \tilde{I'}(u,v,\lambda) ...
...t\int_{-\infty}^\infty p {I'}(s,p,\lambda) {\rm d}s {\rm d}p .
\end{displaymath}

The higher order derivatives can be easily deduced from these expressions. Here we give the second order ones:

\begin{displaymath}{{\partial^2}\over{\partial u^2}} \left. \tilde{I'}(u,v,\lamb...
...int_{-\infty}^\infty s^2 {I'}(s,p,\lambda) {\rm d}s {\rm d}p ,
\end{displaymath}


\begin{displaymath}{{\partial^2}\over{\partial v^2}} \left. \tilde{I'}(u,v,\lamb...
...int_{-\infty}^\infty p^2 {I'}(s,p,\lambda) {\rm d}s {\rm d}p ,
\end{displaymath}


\begin{displaymath}{{\partial^2}\over{\partial u \partial v}}
\left. \tilde{I'}...
...int_{-\infty}^\infty p s {I'}(s,p,\lambda) {\rm d}s {\rm d}p .
\end{displaymath}

Introducing these derivatives into the Eq. (22) we get:

\begin{displaymath}\tilde{I'}(u,v,\lambda)=
\int\int_{-\infty}^\infty {I'}(s,p,...
...int\int_{-\infty}^\infty s {I'}(s,p,\lambda) {\rm d}s{\rm d}p
\end{displaymath}


\begin{displaymath}- i 2 \pi v \int\int_{-\infty}^\infty p {I'}(s,p,\lambda) {\r...
...nt\int_{-\infty}^\infty sp {I'}(s,p,\lambda) {\rm d}s {\rm d}p
\end{displaymath}


 
 \begin{displaymath}- 2 \pi^2 u^2 \int\int_{-\infty}^\infty
s^2 {I'}(s,p,\lambda...
...\infty}^\infty
p^2 {I'}(s,p,\lambda){\rm d}s{\rm d}p + \cdots
\end{displaymath} (23)

References

 
Copyright ESO 2001