© ESO, 2013

1. Introduction

This paper is a companion paper to the observational results published in Habibi et al. (2011), and it focusses on the simulation of the scintillation effects that were searched for. Cold transparent molecular clouds are one of the last possible candidates for the missing baryons of cosmic structures on different scales (Pfenniger & Combes 1994; Pfenniger & Revaz 2005 and McGaugh et al. 2010). Since these hypothesised clouds do not emit or absorb light, they are invisible for the terrestrial observer, so we have to investigate indirect detection techniques. Our proposal for detecting such transparent clouds is to search for the scintillation of the stars located behind the transparent medium, caused by the turbulence of the cloud (Moniez 2003 and Habibi et al. 2011). The objective of this technical paper is to describe the way we can connect observations to scintillation parameters through a realistic simulation. We used these connections in the companion paper (Habibi et al. 2011) to establish constraints both from null results (towards SMC) and from observations pointing to a possible scintillation effect (towards nebula B68). Similar studies of propagation through a stochastic medium followed by Fresnel diffraction have been made by 3 and for use in radio-astronomy by 9.

We first introduce the notations and the formalism in Sect. 2. Then we describe the different stages of the simulation pipeline up to the production of simulated light curves in Sect. 3. We study the observables that can be extracted from the light curve of a scintillating star, and in particular, we check the expected modulation amplitude properties in Sect. 4. The discussion is extended to non-Kolmogorv turbulence cases in Sect. 5. In Sect. 6 we use the results from the simulation pipeline to optimise the observational strategy for discovering scintillating stars, and indicate some perspectives in the conclusion.

Complementary information on observations made with the ESO-NTT telescope and on the analysis based on the present simulations are to be found in our companion paper (Habibi et al. 2011).

2. Basic definitions and formalism

Fig. 1 Geometric configuration. The source is located in the (x2,y2) plane, the screen contains the refractive structure, and the observer is located in the (x0,y0) plane. A1(x1,y1) and ${\mathit{A}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{y}}_{\mathrm{1}}\mathrm{\right)}$ are the amplitudes before and after screen crossing.

The formalism described in this section has been inspired and adapted from the radioastronomy studies (Narayan 1992). But at optical wavelength, the scintillation is primarily due to the refraction through dense clouds of H₂ + He instead of the interaction with the ionised interstellar medium. The origin of the stochastic phase fluctuation experienced by the electromagnetic wave when crossing the refractive medium is the phase excess induced by the stochastic fluctuation of the column density due to the turbulence (Moniez 2003): $\begin{array}{ccc}\mathit{\phi }\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{y}}_{\mathrm{1}}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\mathit{\alpha }}{\mathit{\lambda }}\mathit{Nl}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{y}}_{\mathrm{1}}\mathrm{\right)}\end{array}$(1)where x₁ and y₁ are the coordinates in the cloud’s plane, perpendicular to the sightline (see Fig. 1). Here φ(x₁,y₁) is the phase delay induced to the wavefront after crossing the cloud, Nl(x₁,y₁) is the cloud column density of H₂ molecules plus He atoms along the line of sight, α is the medium polarisability, and λ the wavelength. The phase delay here scales with λ^-1, in contrast to the radioastronomy where it scales with λ. Since other sources of phase delay, such as the geometrical delay induced by scattering from cloud inhomogeneities, are negligible, the thin screen approximation can be used, and the cloud can be considered as a 2D scattering screen whose optical properties are mapped by the phase screen φ(x₁,y₁). The statistical properties of the phase screen are described by the phase structure function D_φ(x₁,y₁). By assuming an isotropic turbulence (Narayan 1992):

$\begin{array}{ccc}& & {\mathit{D}}_{\mathit{\phi }}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{y}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}{\mathit{D}}_{\mathit{\phi }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\\ & & \end{array}$(2)where the first expression is averaged over the plane positions $\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}^{\hspace{0.17em}\mathrm{\prime }}\mathit{,}{\mathit{y}}_{\mathrm{1}}^{\hspace{0.17em}\mathrm{\prime }}\mathrm{\right)}$, $\mathit{r}\mathrm{=}\sqrt{{\mathit{x}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}_{\mathrm{1}}^{\mathrm{2}}}$, β is the turbulence exponent – equals 11/3 for Kolmogorov turbulence – and the diffraction radius, R_diff, is the transverse distance on the phase screen for which the root mean square of the phase variation in one radian. The diffraction radius can be expressed in terms of the cloud parameters (Habibi et al. 2011); assuming Kolmogorov turbulence it is given by $\begin{array}{ccc}{\mathit{R}}_{\mathrm{diff}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathrm{263}\hspace{0.17em}\mathrm{km}{\left[\frac{\mathit{\lambda }}{\mathrm{1}\hspace{0.17em}\mathit{\mu }\mathrm{m}}\right]}^{\frac{\mathrm{6}}{\mathrm{5}}}{\left[\frac{{\mathit{L}}_{\mathit{z}}}{\mathrm{10}\hspace{0.17em}\mathrm{AU}}\right]}^{\mathrm{-}\frac{\mathrm{3}}{\mathrm{5}}}{\left[\frac{{\mathit{L}}_{\mathrm{out}}}{\mathrm{10}\hspace{0.17em}\mathrm{AU}}\right]}^{\frac{\mathrm{2}}{\mathrm{5}}}{\left[\frac{{\mathit{\sigma }}_{\mathrm{3}\mathit{n}}}{{\mathrm{10}}^{\mathrm{9}}\hspace{0.17em}{\mathrm{cm}}^{-3}}\right]}^{\mathrm{-}\frac{\mathrm{6}}{\mathrm{5}}}\mathit{,}& & \end{array}$(3)where L_z is the cloud size along the sightline, L_out is the turbulence outer scale, and σ_3n is the dispersion of the volumic number density in the medium¹. Here we assume the gas to be a mixing of H₂/He with 24% He by mass (corresponding to the primordial abundances) and therefore ⟨ α ⟩ = 0.720 × 10^-30 m³. In this expression, the cloud parameters are scaled to the values given by the Pfenniger-Combes model for the clumpuscules (the tiniest cloudlets of the molecular cloud). In the NIR band, the diffraction radius of a typical clumpuscule is expected to be R_diff ~ 500 km.

The phase statistics of the screen can be equivalently described in Fourier space by the phase spectral density: ${\mathit{S}}_{\mathit{\phi }}\mathrm{\left(}{\mathit{q}}_{\mathit{x}}\mathit{,}{\mathit{q}}_{\mathit{y}}\mathrm{\right)}\mathrm{=}{\mathit{S}}_{\mathit{\phi }}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{diff}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{\beta }\mathrm{-}\mathrm{1}}\hspace{0.17em}\mathit{f}\mathrm{(}\mathit{\beta }\mathrm{)}}\hspace{0.17em}\mathrm{(}{\mathit{R}}_{\mathrm{diff}}\hspace{0.17em}\mathit{q}{\mathrm{)}}^{\mathrm{-}\mathit{\beta }}\mathit{,}$(4)where Fourier coordinates q_x and q_y have inverse length dimension, $\mathit{q}\mathrm{=}\sqrt{{\mathit{q}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{q}}_{\mathit{y}}^{\mathrm{2}}}$, and $\mathit{f}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\frac{{\mathrm{2}}^{\mathrm{-}\mathit{\beta }}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}\mathrm{-}\mathit{\beta }\mathit{/}\mathrm{2}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\mathit{\beta }\mathit{/}\mathrm{2}\mathrm{)}}$ is a constant.

After crossing the cloud, the distorted wavefront of a point source propagates toward the observer and produces an illumination pattern on the observer’s plane given by $\begin{array}{ccc}{\mathit{I}}_{\mathrm{0}}\mathrm{\left(}{\mathit{x}}_{\mathrm{0}}\mathit{,}{\mathit{y}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{s}}}{{\mathit{z}}_{\mathrm{1}}^{\mathrm{2}}}\hspace{0.17em}\mathit{h}\mathrm{\left(}{\mathit{x}}_{\mathrm{0}}\mathit{,}{\mathit{y}}_{\mathrm{0}}\mathrm{\right)}\mathit{,}& & \end{array}$(5)where I₀(x₀,y₀) is the light intensity on the observer’s plane, L_s is the source luminosity, z₁ is the source-screen distance (see Fig. 1), and h(x₀,y₀) is given by the Fresnel-Huygens diffraction integral after considering the Fresnel and the stationary phase approximations (Born & Wolf 2002; Moniez 2003): $\begin{array}{ccc}\mathit{h}\mathrm{\left(}{\mathit{x}}_{\mathrm{0}}\mathit{,}{\mathit{y}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}{\left|\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{F}}^{\mathrm{2}}}\hspace{0.17em}{”}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathrm{e}}^{\mathrm{i}\mathit{\phi }\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{y}}_{\mathrm{1}}\mathrm{\right)}}{\mathrm{e}}^{\mathrm{i}\frac{\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{-}{\mathit{x}}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{y}}_{\mathrm{0}}\mathrm{-}{\mathit{y}}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{R}}_{\mathrm{F}}^{\mathrm{2}}}}\mathrm{d}{\mathit{x}}_{\mathrm{1}}\mathrm{d}{\mathit{y}}_{\mathrm{1}}\right|}^{\mathrm{2}}\mathit{,}& & \end{array}$(6)where R_F = $\sqrt{\mathit{\lambda }\hspace{0.17em}{\mathit{z}}_{\mathrm{0}}\mathit{/}\mathrm{2}\mathit{\pi }}$ is the Fresnel radius² and z₀ is the screen-observer distance. The Fresnel radius can be expressed as $\begin{array}{ccc}{\mathit{R}}_{\mathrm{F}}\mathrm{=}\mathrm{2214}\hspace{0.17em}\mathrm{km}{\left[\frac{\mathit{\lambda }}{\mathrm{1}\mathit{\mu }\mathrm{m}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}\hspace{0.17em}{\left[\frac{{\mathit{z}}_{\mathrm{0}}}{\mathrm{1}\mathrm{kpc}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{·}& & \end{array}$(7)At the typical distance of a halo object (~10 kpc), R_F ~ 7000 km for NIR wavelengths. Because of the motion of the cloud with respect to the Earth-source line-of-sight, the illumination pattern sweeps the observer plane, so that a terrestrial observer receives fluctuating intensity light from the point source. This effect, the scintillation, has two different scattering regimes (Uscinski 1977; Tatarskii & Zavorotnyi 1980), a weak regime (R_diff > R_F) and a strong regime (R_diff < R_F). In the present studies, we concentrate on the strong regime, which clearly is easier to detect, but some realistic configurations may involve the intermediate regime studied by Goodman & Narayan (2006). For the strong regime, there are two different modes of flux variations (Narayan 1992; see also Rickett 1986; Rumsey 1975; Sieber 1982). The first one is the diffractive mode with length scale corresponding to the screen’s scale of phase variations R_diff given by Eq. (3). The resulting “speckles”, with typical size of the order of R_diff, are shown in Fig. 4. The corresponding time scale of the light fluctuations is t_diff = R_diff/V_T: $\begin{array}{ccc}& & {\mathit{t}}_{\mathrm{diff}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathrm{2.6}\mathrm{s}\\ & & \end{array}$(8)where V_T is the sightline relative transverse motion. Therefore fast flux variations are expected with a typical time scale of t_diff ~ few s. The second variation mode is the refractive mode associated to the longer length scale called refraction radius: $\begin{array}{ccc}{\mathit{R}}_{\mathrm{ref}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{\lambda }{\mathit{z}}_{\mathrm{0}}}{{\mathit{R}}_{\mathrm{diff}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}}\mathrm{~}\mathrm{30}\hspace{0.17em}\mathrm{860}\hspace{0.17em}\mathrm{km}\left[\frac{\mathit{\lambda }}{\mathrm{1}\mathit{\mu }\mathrm{m}}\right]\left[\frac{{\mathit{z}}_{\mathrm{0}}}{\mathrm{1}\hspace{0.17em}\mathrm{kpc}}\right]{\left[\frac{{\mathit{R}}_{\mathrm{diff}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}}{\mathrm{1000}\hspace{0.17em}\mathrm{km}}\right]}^{-1}\mathrm{·}& & \end{array}$(9)This natural length scale corresponds to the size, in the observer’s plane, of the diffraction spot from a patch of R_diff(λ) in the screen’s plane. This is also the size of the region in the screen where most of the scattered light seen at a given observer’s position originates. Our convention for R_ref differs from Narayan (1992) by a factor 2π since it emerges naturally from the Fourier transform we use for calculating the illumination pattern (see formula (14)), and it also matches the long distance-scale flux variations visible in Fig. 4. The corresponding time scale is given by t_ref = R_ref/V_T: ${\mathit{t}}_{\mathrm{ref}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{\simeq }\mathrm{5.2}\hspace{0.17em}\min \left[\frac{\mathit{\lambda }}{\mathrm{1}\mathit{\mu }\mathrm{m}}\right]\left[\frac{{\mathit{z}}_{\mathrm{0}}}{\mathrm{1}\hspace{0.17em}\mathrm{kpc}}\right]{\left[\frac{{\mathit{R}}_{\mathrm{diff}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}}{\mathrm{1000}\hspace{0.17em}\mathrm{km}}\right]}^{-1}{\left[\frac{{\mathit{V}}_{\mathrm{T}}}{\mathrm{100}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}}\right]}^{-1}\mathrm{·}$(10)

3. Simulation description

In this section, we describe the simulation pipeline with some numerical tricks, from the generation of the phase screen induced by turbulent gas, up to the light versus time curves expected from realistic stars seen through this gas. The steps in this pipeline are listed below:

• The simulation of the refractive medium: inSect. 3.1 we describe the generation of a phasescreen and examine the impact of the limitations caused by thesampling and by the finite size of the screen by comparing theinitial (theoretical) and reconstructed diffraction radii.

• The computation of the illumination pattern: in Sect. 3.2, we first describe the calculation of the illumination pattern produced on Earth by a point monochromatic source as seen through the refractive medium. We explain the technique for avoiding numerical (diffraction) artefacts caused by the borders of the simulated screen, and discuss the criterion on the maximum pixel size to avoid aliasing effects. The pattern computation is then generalised to extended polychromatic sources.

• The light curve simulation: we describe in Sect. 3.3 the simulation of the light fluctuations with time at a given position induced by the motion of the refractive medium with respect to the line of sight.

3.1. Simulation of the phase screen

Fig. 2 The phase-delay variations near the average for a simulated refractive screen with Nx × Ny = 20 000 × 20 000 pixels, Δ1 = 22.6 km, and Rdiff = 100 km. The grey scale ranges between ± 50 × 2π rad (clear regions correspond to an excess of phase with respect to the average). The zoom (inset in the lower-right corner) illustrates the self-similarity of the simulated screen (grey scale amplitude of 5 × 2π rad).

Numerical realisations of the 2D phase screen – made of N_x × N_y squared pixels of size Δ₁ – are randomly generated from the phase spectral density S_φ(q_x,q_y), determined by the choice of R_diff in relation (4). Such a phase screen – with the desired statistical properties – is obtained from the random realisation of a Fourier transform F_φ(q_x,q_y) in a way that makes the ensemble of such realisations satisfy the relation ${\mathrm{⟨}\mathrm{\left|}{\mathit{F}}_{\mathit{\phi }}\mathrm{\right(}{\mathit{q}}_{\mathit{x}}\mathit{,}{\mathit{q}}_{\mathit{y}}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathrm{⟩}}_{\mathrm{realisations}}\mathrm{=}{\mathit{L}}_{\mathit{x}}{\mathit{L}}_{\mathit{y}}{\mathit{S}}_{\mathit{\phi }}\left(\mathit{q}\mathrm{=}\sqrt{{\mathit{q}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{q}}_{\mathit{y}}^{\mathrm{2}}}\right)\mathit{,}$(11)where L_x and L_y are the screen physical size, and the average is the ensemble averaging over the different realisations. For each $\mathrm{\left(}{\mathit{q}}_{\mathit{x}}^{\mathit{i}}\mathit{,}{\mathit{q}}_{\mathit{y}}^{\mathit{j}}\mathrm{\right)}$ vector associated to a pixel (i,j), we generate a random complex number ${\mathit{F}}_{\mathit{\phi }}\mathrm{\left(}{\mathit{q}}_{\mathit{x}}^{\mathit{i}}\mathit{,}{\mathit{q}}_{\mathit{y}}^{\mathit{j}}\mathrm{\right)}\mathrm{=}\mathrm{(}{\mathit{f}}_{\mathrm{Re}}^{\mathit{i,j}}\mathrm{+}\mathrm{i}{\mathit{f}}_{\mathrm{Im}}^{\mathit{i,j}}\mathrm{)}\mathit{/}\sqrt{\mathrm{2}}$ where each component ${\mathit{f}}_{\mathrm{Re}}^{\mathit{i,j}}$ and ${\mathit{f}}_{\mathrm{Im}}^{\mathit{i,j}}$ spans the Gaussian distribution³ of zero mean, and ${\mathit{L}}_{\mathit{x}}{\mathit{L}}_{\mathit{y}}{\mathit{S}}_{\mathit{\phi }}\mathrm{(}\mathit{q}\mathrm{=}\sqrt{{\mathit{q}}_{\mathit{x}}^{\mathit{i}\hspace{0.17em}\mathrm{2}}\mathrm{+}{\mathit{q}}_{\mathit{y}}^{\mathit{j}\hspace{0.17em}\mathrm{2}}}\mathrm{)}$ width. In this way, relation 11 is automatically satisfied when averaging on a large number of such realisations, since the average of | F_φ(q_x,q_y)|² for the ensemble of (q_x,q_y) vectors with the same module q equals L_xL_yS_φ(q) by construction. The phase screen φ(x₁,y₁) is finally obtained by numerically computing the inverse Fourier transform of this phase spectrum F_φ(q_x,q_y). Figure 2 shows a phase screen generated by assuming Kolmogorov turbulence (β = 11/3).

3.1.1. Preliminary checks, limitations

To check the accuracy of the numerically generated phase screen (Fig. 2), we recomputed the phase structure function ${\mathit{D}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ (and consequently ${\mathit{R}}_{\mathrm{diff}}^{\mathrm{rec}}$) from the generated phase Fourier transform F_φ(q_x,q_y), and we compared it with the theoretical phase structure function (Eq. (2)). First, the spectral density is recomputed from the generated F_φ(q_x,q_y) using relation $\begin{array}{ccc}{\mathit{S}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\frac{\mathrm{⟨}\mathrm{\left|}{\mathit{F}}_{\mathit{\phi }}\mathrm{\right(}{\mathit{q}}_{\mathit{x}}\mathit{,}{\mathit{q}}_{\mathit{y}}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{q}\mathrm{=}\sqrt{{\mathit{q}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{q}}_{\mathit{y}}^{\mathrm{2}}}}}{{\mathit{L}}_{\mathit{x}}\hspace{0.17em}{\mathit{L}}_{\mathit{y}}}\mathit{,}& & \end{array}$(12)where the average is performed on the (q_x,q_y) coordinates⁴ spanning the circle of radius $\mathit{q}\mathrm{=}\sqrt{{\mathit{q}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{q}}_{\mathit{y}}^{\mathrm{2}}}$. The corresponding phase auto-correlation function is then given by Fourier transform:

$\begin{array}{ccc}{\mathit{\xi }}^{\mathrm{rec}}\mathrm{\left(}{r}\mathrm{\right)}& \mathrm{=}& ”{\mathit{S}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{\left(}{q}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}\hspace{0.17em}{q}\mathit{.}{r}}\mathrm{d}{q}\mathit{,}\\ {\mathit{\xi }}^{\mathrm{rec}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}{\mathit{q}}_{\max }\\ {\mathit{q}}_{\min }\end{array}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{q}\hspace{0.17em}{\mathit{S}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{\left(}\mathit{q}\mathrm{\right)}{\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathrm{i}\hspace{0.17em}\mathit{qr}\cos \mathit{\theta }}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{q}\\ & \mathrm{=}& {\mathrm{\int }}_{{\mathit{q}}_{\min }}^{{\mathit{q}}_{\max }}\mathrm{2}\mathit{\pi q}\hspace{0.17em}{\mathit{S}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\hspace{0.17em}{\mathit{J}}_{\mathrm{0}}\mathrm{(}\mathrm{2}\mathit{\pi }\mathit{qr}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{q,}\end{array}$(13)where J₀ is the Bessel function. The recomputed structure function is then given by ${\mathit{D}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ = 2(ξ^rec(0) − ξ^rec(r)), and the value of ${\mathit{R}}_{\mathrm{diff}}^{\mathrm{rec}}$ is deduced from its definition ${\mathit{D}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{\left(}{\mathit{R}}_{\mathrm{diff}}^{\mathrm{rec}}\mathrm{\right)}\mathrm{=}\mathrm{1}$.

Fig. 3 Phase structure functions Dφ(r) for a phase screen with Rdiff = 500 km. Blue line is the initial (theoretical) structure function. Red line is reconstructed from one of the realisations of the phase screen through simulation. The black curve is obtained from the numerical integration of the initial phase spectral density sampled as in the simulation.

In Fig. 3, we show the theoretical phase structure function of a turbulent medium with R_diff = 500 km – for which D_φ(r = 500 km) = 1 by definition –, and the recomputed (effective) structure function from one of the realisations of the screen. From this recomputed function, we find ${\mathit{R}}_{\mathrm{diff}}^{\mathrm{rec}}\mathrm{\approx }\mathrm{540}$ km since ${\mathit{D}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{(}\mathit{r}\mathrm{\approx }\mathrm{540}\hspace{0.17em}\mathrm{km}\mathrm{)}$ = 1. To find the origin of the difference with the input value R_diff = 500 km, we replaced ${\mathit{S}}_{\mathit{\phi }}^{\mathrm{rec}}\mathrm{\left(}\mathit{q}\mathrm{\right)}$ in Eq. (13) by the theoretical spectrum S_φ(q) sampled as in the simulation (number of pixels N_x × N_y ~ 14 000 × 14 000 with pixel size Δ₁ = 28.85 km). Then we computed the integral (13) numerically with the same integration limits (q_min,q_max)⁵ as the simulation. The integration result differs only by a few percent from the function recomputed from the simulated screen. We showed that the black curve approaches the blue curve when q_min → 0 and q_max → ∞. This means that the sampling is mainly responsible for the difference between D_φ and ${\mathit{D}}_{\mathit{\phi }}^{\mathrm{rec}}$. Since our simulation is limited by the number of pixels, we lose the contributions of the large and small scales in the recomputed R_diff. The only way to push back this limitation is to generate larger screens (larger N_x and N_y) with higher resolutions (smaller Δ₁) to cover wider interval of spatial frequencies which in return needs higher computational capacities (see also Sect. 3.4).

3.2. Illumination pattern

To obtain the illumination pattern on the observer plane, we numerically compute the integral (6) which can be written as a Fourier transform: $\begin{array}{ccc}\mathit{h}\mathrm{\left(}{\mathit{x}}_{\mathrm{0}}\mathit{,}{\mathit{y}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{F}}^{\hspace{0.17em}\mathrm{2}}}|\mathit{FT}\left[\mathit{G}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{y}}_{\mathrm{1}}\mathrm{\right)}\right]{|}^{\mathrm{2}}\hspace{0.17em}\left({\mathit{f}}_{\mathit{x}}\mathrm{=}\frac{{\mathit{x}}_{\mathrm{0}}}{\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{F}}^{\mathrm{2}}}\mathit{,}\hspace{0.17em}{\mathit{f}}_{\mathit{y}}\mathrm{=}\frac{{\mathit{y}}_{\mathrm{0}}}{\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{F}}^{\mathrm{2}}}\right)\mathit{,}& & \end{array}$(14)where $\begin{array}{ccc}\mathit{G}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{y}}_{\mathrm{1}}\mathrm{\right)}& \mathrm{=}& \end{array}$(15)Here, (x₁,y₁) and (x₀,y₀) are the screen and observer coordinates, respectively. Coordinates (f_x,f_y) are the conjugated variables in Fourier space. Before computing expression (14), a regularisation procedure for G(x₁,y₁) has been defined to avoid computational artefacts.

3.2.1. Screen regularisation

Since integral (14) is computed numerically, the coordinates (x₁,y₁) have discrete (integer) values describing pixel position centres on the screen, to allow simple combinations of illumination patterns with different pixel sizes (corresponding to different wavelengths). That the integration domain is limited is physically equivalent to computing the Fresnel integral within a diaphragm with the size of the screen. In this case, we face a parasitic effect: the light diffraction from the sharp edges of the diaphragm. This causes rapid intensity variations at the borders of the observer plane. To attenuate this effect and remove the resulting diffraction fringes, we multiply the screen intensity transmission by a 2D smoothing function. We define the following 1D smoothing function SF(x) (see Fig. 6): $\begin{array}{ccc}\mathit{SF}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \\ & \\ & \end{array}& & \end{array}$where L_m = 10R_F is the margin length from the borders of the screen with size L. We multiply the function G(x₁,y₁) by SF(x₁) × SF(y₁) in expression (14). Since the Fresnel integral is dominated by the contribution of the integrand within a disk that is within a few Fresnel radii, it is sufficient to smooth the discontinuity of G(x₁,y₁) within a distance of a few Fresnel radii (here L_m = 10 R_F). We tested the efficiency of this regularisation procedure by checking that the simulated illumination pattern from a point-source projected through a uniform phase screen was – as expected – also uniform beyond our precision requirements (<1%) (Habibi 2011). To define a reliable fiducial domain in the observer plane excluding the regions that are only partially illuminated owing to the diaphragm, we also delimited an R_ref/2 margin from the borders of the illumination pattern, corresponding to the typical radius (half-size) of the large-scale luminous spots (see also Coles et al. 1995). Figure 4 shows the pattern produced by a point source through a turbulent medium with R_diff = 100 km located at z₁ = 160 pc from the Earth at wavelength λ = 2.162 μm. Since the corresponding Fresnel radius R_F = 1300 km is larger than the diffraction radius, the regime is the strong scintillation regime. The hot speckles (of typical size R_diff ~ 100 km) can be distinguished from the larger dark/luminous structures that have a typical size of ${\mathit{R}}_{\mathrm{ref}}\mathrm{~}\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{F}}^{\mathrm{2}}\mathit{/}{\mathit{R}}_{\mathrm{diff}}\mathrm{\simeq }\mathrm{100}\hspace{0.17em}\mathrm{000}$ km.

Fig. 4 Typical illumination pattern from a point-like source. Here Rdiff = 100 km, the screen is at z0 = 160 pc, λ = 2.16 μm, then RF = 1300 km and Rref = 106 000 km. The typical length scale of the small-scale speckles is Rdiff, and the scale of the larger structures is Rref. The white square shows our fiducial zone with a margin of Lm = Rref/2 from the borders. Grey scale range from 0 to 4 times the mean intensity. The image has 20 000 × 20 000 pixels, each with a 22.6 km side.

Fig. 5 Typical illumination pattern from an extended source produced through the same screen as in Fig. 4, with source radius rs = 0.5 R⊙, located at z0 + z1 = 1 kpc + 160 pc (Rs ≃ 53 000 km). The small-scale speckles are smeared, and only the larger scale fluctuations survive. The white square shows the restricted fiducial zone with margin of Lm = Rref/2 + Rs from the borders. Grey scale ranges ± 20% around the mean intensity.

Fig. 6 The smoothing function SF(x). L is the screen size, Lm is the margin from the screen borders.

3.2.2. Effect of sampling

Here, we discuss some limitations caused by the pixellisation. The screen should be sampled often enough to avoid the aliasing effects. Aliasing happens when G(x₁,y₁) contains frequencies that are higher than the Nyquist frequency f_Nyq = 1/(2Δ₁), where Δ₁ is the pixel size. In relation (15), G contains two length scales, the diffraction and the Fresnel radii (R_diff and R_F). R_diff is the characteristic length of the phase screen variations φ(x₁,y₁); it is at least necessary that Δ₁ < R_diff/2 in order to sample phase variations up to 1/R_diff spatial frequency. R_F appears in the quadratic term $\exp \left[\mathrm{i}\hspace{0.17em}\frac{{\mathit{x}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{\mathit{R}}_{\mathrm{F}}^{\hspace{0.17em}\mathrm{2}}}\right]$; the oscillation of this term accelerates as x₁ and y₁ increase, and aliasing occurs as soon as the distance between two consecutive peaks is smaller than 2Δ₁. In one dimension we therefore expect aliasing if $\begin{array}{ccc}\frac{\mathrm{(}{\mathit{x}}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathrm{\Delta }}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}{\mathit{x}}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{\mathit{R}}_{\mathrm{F}}^{\mathrm{2}}}\mathit{>}\mathrm{2}\mathit{\pi ,}& & \\ \mathrm{i}\mathit{.}\mathrm{e}\mathit{.}\frac{{\mathit{x}}_{\mathrm{1}}}{{\mathrm{\Delta }}_{\mathrm{1}}}\mathit{>}\mathit{\pi }{\left[\frac{{\mathit{R}}_{\mathrm{F}}}{{\mathrm{\Delta }}_{\mathrm{1}}}\right]}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathit{,}& & \end{array}$(16)where x₁/Δ₁ is the distance to the optical axis, expressed in pixels. The condition Δ₁ = R_F/2 would be obviously insufficient here to avoid aliasing, since beyond ~11 pixels only from the optical axis (x₁/Δ₁ > 11 pixel), the quadratic term would be under-sampled. In practice, for the configurations considered in this paper, the Fresnel radius is >1000 km. When Δ₁ = 22 km in the simulation, the aliasing starts at x₁/Δ₁ > 6490 pixels from the centre of the image, corresponding to more than 140R_F, which is large enough to cover the sensitive domain related to the stationary phase approximation.

3.2.3. Extended source (spatial coherency)

The illumination pattern of a scintillating extended source is given by the convolution product of the illumination pattern of the point-like source with the projected source limb profile (Moniez 2003; & Habibi 2011), $\begin{array}{ccc}{\mathit{I}}_{\mathrm{ext}}& \mathrm{=}& \end{array}$(17)where L_s is the source luminosity, and the normalised limb profile is described as a uniform disk: $\begin{array}{ccc}{\mathit{P}}_{\mathit{r}}\mathrm{\left(}{\mathit{x}}_{\mathrm{0}}\mathit{,}{\mathit{y}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$where ${\mathit{R}}_{\mathrm{s}}\mathrm{=}\frac{{\mathit{z}}_{\mathrm{0}}}{{\mathit{z}}_{\mathrm{1}}}\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}$ is the projected source radius on the observer’s plane, and r_s is the source radius. Figure 5 shows the convolution of the pattern of Fig. 4 with the projected profile of a star with radius r_s = 0.5 R_⊙, located at z₀ (160 pc) +z₁ (1 kpc) = 1.16 kpc (R_s = 53 000 km). High-frequency fluctuations due to diffractive speckle disappear, and the pattern loses contrast, but the variations on R_ref scale remain visible. As the convolution involves a disk of radius R_s, we can perform the calculation only at a distance greater than R_s from the borders. We therefore define a new fiducial zone by excluding a margin of R_ref/2 + R_s from the initial borders. Any statistical analysis will be made within this zone to be safe from any border perturbation.

Fig. 7 Simulated illumination maps (20 000 × 20 000 pixel of 22.6 km side) produced on Earth by a source located at z0 + z1 = 1.18 kpc through a refracting cloud assumed to be at z0 = 160 pc with a turbulence parameter Rdiff(2.16 μm) = 100 km. Here Rref(2.16 μm) ≃ 100 000 km. Top-left and middle: illumination produced at λ = 2.16 μm from a point-source with a zoomed detail; the contrast is 100%. The grey scale ranges from 0 to 4 times the mean intensity. Top-right: the same from a K0V star (rs = 0.85 R⊙, MV = 5.9, at 1.18 kpc V = 16.3). The circle shows the projection of the stellar disk (RS = rs × z0/z1). Here the modulation index is only 3.3%, and the grey scale ranges from ±20% around the mean intensity. The bottom maps are the illuminations in Ks wide band (λcentral = 2.162 μm, Δλ = 0.275 μm), using the same grey scales as above. The modulation index is 55% for the point-source (left and centre) and 3.3% for the extended source (right). The two parallel straight lines show the sections sampled by two observers located about 10 000 km apart, when the screen moves with the transverse velocity VT.

3.2.4. Polychromatic source (time coherency)

The illumination patterns shown in Figs. 4 and 5 are computed for a monochromatic source (fixed λ), but observations are done through filters with finite-width passbands. To take the contributions of different wavelengths to the pattern into account, we superimpose the illumination patterns obtained with the same refractive structure (the same column density fluctuations) at different wavelengths⁶. We have considered the passband of the SOFI camera in K_s band and approximated it as a rectangular function over the transmitted wavelengths with a central value 2.162 μm and width 0.275 μm. Twenty-one illumination patterns were computed for 21 regularly spaced wavelengths within the [2.08, 2.28] μm interval, and were co-added to simulate the illumination pattern through the K_s passband. We checked that the spacing between successive wavelengths was small enough to produce a co-added image with a realistic residual modulation index, by studying this index as a function of the number of monochromatic components equally spaced within the full bandwidth. We found that an asymptotic value is reached with a co-added image made up of only about ten components.

Figure 7 (left and centre) shows a comparison between monochromatic (up) and K_s passband (down) illumination patterns of a point-like source. The speckle pattern is attenuated when the light is not monochromatic (or equivalently when time coherency is limited). This is caused by the small decorrelation bandwidth δλ_dec of the strong diffractive scintillation regime, which is given by $\mathit{\delta }{\mathit{\lambda }}_{\mathrm{dec}}\mathit{/}\mathit{\lambda }\mathrm{~}\mathrm{\left(}{\mathit{R}}_{\mathrm{diff}}\mathrm{\right(}\mathit{\lambda }\mathrm{)}\mathit{/}{\mathit{R}}_{\mathrm{F}}\mathrm{)}2$(18)(Narayan 1992; Gwinn et al. 1998). This chromatic effect results from the high sensitivity of the constructive interference condition with the wavelength⁷. In the case of Fig. 7, the decorrelation bandwidth is δλ_dec/λ ~ 1%. Since the K_s passband is δλ/λ ~ 0.1 ~ 10 × δλ_dec/λ, a first-order estimation of the modulation index for a K_s passband pattern is the value obtained when adding ten decorrelated patterns of speckle i.e. one multiplied by ~$\sqrt{\mathrm{1}\mathit{/}\mathrm{10}}$. The modulation index found with our simulation (~55%) has the correct order of magnitude.

By contrast, structures of size of R_ref are much less sensitive to the variations in λ, according to the combination of (3) and (9): $\frac{\mathrm{\Delta }{\mathit{R}}_{\mathrm{ref}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{ref}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{5}}\mathrm{\times }\frac{\mathrm{\Delta }\mathit{\lambda }}{\mathit{\lambda }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{5}}\mathrm{\times }\mathrm{10}\mathrm{\%}\mathrm{~}\mathrm{2}\mathrm{\%}\mathit{.}$(19)As a consequence of the large-scale smearing of the point-source pattern when considering an extended source, there is no significant difference between monochromatic and polychromatic patterns for such an extended source, as shown in Fig. 7 (right). Therefore, in the following we ignore the impact of the R_diff structures.

3.3. Simulation of light curves

What we observe with a single telescope is not the 2D illumination pattern but a light curve. Because of the relative motions, the telescope sweeps a 1D section of the pattern, at a constant velocity as long as parallax can be neglected. We therefore simulated light versus time curves by sampling the 2D pattern pixels along straight lines, with the relative speed of the telescope.

3.4. Computing limitations

We adopted the same number of pixels N × N and the same pixel scale Δ₁ to numerically describe both the screen’s and the observer’s planes; indeed, the light emerging from the screen is essentially contained in its shadow, and this choice optimises the screen filling. Since the two planes are conjugated, the following relation arises: $\mathit{N}{\mathrm{\Delta }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{F}}^{\mathrm{2}}\mathit{.}$(20)Because of this relation, N has to be as large as possible to minimize Δ₁ and also to get a wide useful area. This area is indeed restricted by the definition of the fiducial domain, which can be heavily reduced when simulating the illumination from an extended source. Also a large fiducial domain is essential when simulating long-duration light curves. As a consequence, memory limitations affect the maximum size of the screen and illumination 2D patterns. In the present paper, we were limited to patterns of 20 000 × 20 000 pixels.

4. Observables

The observable parameters of the scintillation process are the modulation index and the modulation’s characteristic time scales. The main observable we used in our data analysis (Habibi et al. 2011) is the modulation index. It is defined as the flux dispersion σ_I divided by the mean flux : . We first compare the effective modulation index of simulated screens with the theoretical expectations, then examine the precision we can reach on m from a light curve, and give a numerical example. We also show how we can connect the observed modulation index to the geometrical parameters (R_diff(λ), λ, R_s, etc.) through simulation.

Fig. 8 The effective intensity modulation index for simulated scintillating stellar illumination patterns as a function of the theoretical modulation index. The blue dots show the effective modulation index values for different realisations of the phase screens. The red squares represent the mean value of the effective modulation indices with the same theoretical value. This plot shows the expected agreement for m < 0.15, a domain where the simulated screens are large enough not to suffer from statistical biases (see text).

4.1. Modulation Index

For a point-like source in the strong scintillation regime the modulation index m ≈ 1. For an extended source (radius r_s, projected radius R_s = r_s.z₀/z₁) in the same regime, we always have m < 1, and Narayan (1992) showed that when the small-scale (diffractive) speckle is completely smeared, $\begin{array}{ccc}\mathit{m}\mathrm{\approx }{\left[\frac{{\mathit{R}}_{\mathrm{diff}}}{{\mathit{R}}_{\mathrm{F}}}\right]}^{\mathrm{1}\mathit{/}\mathrm{3}}{\left[\frac{{\mathit{\theta }}_{\mathrm{ref}}}{\mathrm{2}\mathit{\pi }{\mathit{\theta }}_{\mathrm{s}}}\right]}^{\mathrm{7}\mathit{/}\mathrm{6}}& & \end{array}$(21)in the case of the Kolmogorov turbulence. Here, θ_ref = R_ref/z₀ is the angular refraction radius⁸ and θ_s the source angular radius. In this expression, Narayan (1992) assumed that z₀ ≪ z₁, therefore z₀ + z₁ ≈ z₁. We use the following expression, which is formally identical to the previous one for z₀ ≪ z₁ but can also be used when z₀ is not negligible: $\begin{array}{ccc}\mathit{m}& \mathrm{\approx }& {\left[\frac{{\mathit{R}}_{\mathrm{diff}}}{{\mathit{R}}_{\mathrm{F}}}\right]}^{\mathrm{1}\mathit{/}\mathrm{3}}{\left[\frac{{\mathit{R}}_{\mathrm{ref}}}{\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{s}}}\right]}^{\mathrm{7}\mathit{/}\mathrm{6}}\\ & \mathrm{\approx }& \end{array}$This relation can be qualitatively justified by noticing firstly that R_diff quantity has to be considered relatively to R_F when considering its impact on diffraction, and secondly that the convolution of the point-source pattern (see Fig. 7-up) with the projected stellar profile (following expression (17)) makes the contrast decrease when the number of R_ref-size domains within R_S = r_Sz₀/z₁ increases. The “exotic” exponents are related to the Kolmogorov turbulence. It should be emphasised that this relation assumes the scintillation is strong and quenched (R_s/R_F > R_F/R_diff > 1).

We checked this relation by performing series of simulations with different phase screens and stellar radii. We generated series of screens with R_diff = 50 km to 500 km by steps of 50 km. For each screen, we considered different sources –at the same geometrical distances– with radii from 0.25 R_⊙ to 1.5 R_⊙ by steps of 0.25 R_⊙ and computed corresponding illumination pattern realisations. The modulation indices were estimated within the fiducial zone for each 2D illumination pattern. In Fig. 8, the thus estimated m for each generated pattern are plotted as a function of the theoretical value expected from expression (23). We note that the modulation has relatively large scatter. This is because the fiducial zone is not large enough, and it contains a limited number of regions with sizes of R_ref (the scale that dominates the light variations). In some cases there are very few distinct dark/luminous regions (as in Fig. 5). The number of such regions within the fiducial domain is ${\mathit{N}}_{\mathrm{R}}\mathrm{~}{\mathit{d}}^{\mathrm{2}}\mathit{/}{\mathit{R}}_{\mathrm{ref}}^{\mathrm{2}}$, where d is the size of the fiducial zone. After discarding the cases with N_R < 5, we computed the mean and root mean square of the effective m values, for each series of configurations with the same theoretical modulation index. Figure 8 shows that relation (23) is satisfied for a modulation index smaller than 0.15. The higher values of m are systematically underestimated in our simulation owing to the limitation of the screen size, which has a stronger impact on the number of large dark/luminous regions in this part of the figure. This number is indeed smaller (though it is larger than 5), and the chances for sampling the deepest valleys and the highest peaks are fewer, therefore biasing the modulation index towards low values.

Therefore, by measuring the modulation index from an observed light curve and using an estimate of star type (i.e. radius) and distance z₀ + z₁ ≈ z₁, we can constrain R_diff(λ) from this relation (23)⁹, even with poor knowledge of the screen’s distance z₀, considering the slow dependence with this parameter. This technique allowed us to infer constraints in Habibi et al. (2011) on the gas turbulence within galactic nebulae (z₀ ~ 80−190 pc and z₁ = 8 kpc), and upper limits on hidden turbulent gas within the galactic halo (assuming z₀ ~ 10 kpc and z₁ + z₀ = 62 kpc).

Fig. 9 Light curves extracted along the 3 horizontal white lines for two illumination patterns. Left column: 2D pattern from a point-like source in Ks band with Rdiff = 300 km and Rref ≈ 28 000 km. The modulation indices of the three light curves differ by less than 5% from the 2D pattern modulation index. Right column: the illumination pattern for an extended source with Rs ≈ 41 000 km, through the same refractive screen. The modulation indices fluctuate by more than 30% around the 2D pattern index, implying the necessity of longer light curves for a better statistical representativity. The distance scale is common to both patterns. The circle shows the projected star disk. The 3 light curves from the right column are not completely decorrelated, because of their common proximity to the same large positive fluctuation.

4.2. Information from the light curves

If a light curve is long enough (or equivalently if the observation time is long enough), its series of light measurements represents an unbiased subsample of the 2D pattern. For instance, the left hand panels of Fig. 9 represent the illumination pattern of a point-like source and three associated light curves. Here, R_ref ≈ 28 000 km and the modulation index m_point = 1.18. The light curves are extracted from three horizontal parallel lines with lengths of ~3.5 × 10⁵ km. The corresponding time scale depends on the relative transverse velocity. The lines are selected far from each other in order not to be affected by the fluctuations of the same regions. Modulation indices along the light curves differ from the 2D’s by less than 5%. When the modulation is characterised by both R_ref and R_diff, a light curve that spans a few R_ref is a sample of the 2D screen, which is large enough to provide a good approximation of the scintillation modulation index for a point-like source. The right hand panels of Fig. 9 show the 2D pattern for an extended source with R_s ≈ 41 000 km. Here, R_s > R_ref and the flux fluctuations are smoothed, characterised by the unique length scale R_ref, and have a much smaller modulation index m_extended = 0.04. The light curves are extracted in the same way as the point-like source within the corresponding restricted fiducial zone (see Fig. 5), therefore they span ~2.5 × 10⁵ km and statistically include a little bit less than 10 R_ref-scale variations. Because of this statistically short length, the light-curve-to-light-curve estimates of m_extended typically fluctuate by ~$\mathrm{1}\mathit{/}\sqrt{\mathrm{10}}\mathrm{\approx }\mathrm{30}\mathrm{\%}$. The fluctuations on m_extended estimates can only be reduced with longer light curves. To get a precision value of 5% on m_extended, the light curve should indeed be as long as ~400 × R_ref. As an example, if we observe through a turbulent gaseous core with R_diff = 200 km in B68 nebula located at z₀ = 80 pc at λ = 2.16 μm (R_ref ~ 27 000 km), and assuming V_T ~ 20 km s^-1, an observing time of ~150 h is needed to measure the modulation index with this precision of 5%. When searching for unseen turbulent media located at unknown distances from us, the diffraction and refraction radii are unknown, and we can only obtain a probability distribution of the observation time for a requested precision on the modulation index.

Fig. 10 Simulation of illumination scintillation patterns with associated light curves from an extended source assuming a refractive screen with non-Kolmogorov turbulences (left β = 3.1, right β = 3.9, see text). Here Rdiff = 100 km, Rref ≈ 8 × 104 km and the projected star radius is RS ~ 0.36Rref = 28 800 km.

Fig. 11 Top: two different phase screen power laws. Bottom: the corresponding modulation indices as a function of the expected modulation index for the Kolmogorov turbulence. The 3 indices plotted at a given abscissa correspond to screens with β = 3.1 (red), 3.67 (Kolmogorov, blue, along the diagonal), and 3.9 (black) with the same Rdiff. As the power law gets steeper, a larger modulation is expected. Blue dots show the Kolmogorov turbulence case.

The other important information carried by the light curves are the characteristic time scale between peaks t_ref = R_ref/V_T, associated to the refraction radius (Fig. 7 upper left), a correlation duration t_s = R_s/V_T, associated to the source radius (Fig. 7 upper right) and possibly t_diff = R_diff/V_T, associated to the small diffractive speckle structure (Fig. 7 lower left); the latter could be detected in exceptionally favourable cases (source with very small angular size, assuming the strong diffractive scintillation regime) with a powerful detection setup such as LSST¹⁰. The time power spectrum of the expected stochastic light curve should show a slope break around frequency ${\mathit{t}}_{\mathrm{ref}}^{-1}$, and possibly around ${\mathit{t}}_{\mathrm{diff}}^{-1}$ (Goodman & Narayan 1985), which should allow one to distinguish it from purely randomly fluctuating light curves due to photometric noise, and to extract constraints on the scintillation configuration. Estimating t_s would be challenging, since the extended source acts as a low passband filter on the point-source pattern; therefore, the imprint of the projected radius is essentially within the attenuation of the light curve autocorrelation (Goodman & Narayan 2006), and more specifically – as mentioned above – within the modulation index, that can be polluted by many observational artefacts. These potentialities of the time studies need more investigation and should be discussed in more detail in a forthcoming paper.

5. Probing various turbulence laws

Up to now, we have focussed on the standard Kolmogorov turbulence. In this section, we investigate a possible deviation from the Kolmogorov turbulence law¹¹. We chose two different phase spectra with β = 3.1 and β = 3.9 in relation (4). To study the corresponding scintillation modes, we generated two series of phase screens according to the spectra and computed the illumination patterns from an extended source with R_s ~ 0.36 R_ref through each of them. The patterns are represented in Fig. 10 with corresponding light curve samples. For both images, R_diff = 100 km, R_F = 1150 km, and R_ref ≈ 8 × 10⁴ km. The pattern with β = 3.1 (left) shows a small modulation index m(3.1) = 0.04, while the pattern with β = 3.9 (right) shows a much larger modulation index m(3.9) = 0.22. As can be seen from its visual aspect, the turbulence with a higher exponent produces stronger contrast on large scales compared to the other one (β = 3.1). To understand the origin of the difference, we compared the two phase spectra at the top of Fig. 11. The steeper spectrum (β = 3.9) has more power for fluctuations on large scales. Moreover, by integrating Eq. (4), we computed that the total power distributed from R_ref to R_diff is about an order of magnitude larger for β = 3.9 than for β = 3.1. Therefore, the larger the β exponent is the stronger flux fluctuations are produced. As a conclusion, the detection of scintillation should be easier for turbulences with steeper spectrum, because a larger modulation index is expected. This is illustrated at the bottom of Fig. 11 where we plot the modulation indices produced by three screens with different β values (including the Kolmogorov case), as a function of the “classical” Kolmogorov turbulence expected index. In this plot, each abscissa corresponds to a distinct R_diff value, which is common to the three screens. By increasing β, we increase the modulation, especially for high m values.

6. Discussion: guidelines provided by the simulation

The observations analysed in (Habibi et al. 2011) were intrepreted by using our simulation pipeline. But we have used our simulation not only to establish connections between the observed light curves and the scintillation configuration, but also to define observing strategies as follows.

Firstly, a correct sensitivity to the scintillation needs the ability to sample, with <1% photometric precision at a sub-minute rate, the light curves of small distant background stars (M ~ 20−21), which have a projected radius small enough to allow for a modulation index of a few percent (typically <R_⊙ at 10 kpc). Our study of the time coherency shows that the usual large passband filters can be used without significant loss of modulation index. Since the optical depth of the process is unknown, a large field of view seems necessary for the exploratory observations, either toward extragalactic stellar sources within LMC or SMC or through known gaseous nebulae. To summarise, an ideal setup for searching for scintillation with series of sub-minute exposures would be a ~4 m class telescope equipped with a fast readout wide field camera and a standard filter (optical passband to search for invisible gas towards extragalactic sources, infrared to observe stars through visible dusty nebulae).

Secondly, our work on the simulation provides us with a guideline to find an undisputable signature of scintillation. The first possibility consists in the search for chromatic effects. Subtle chromatic effects betwen the different regimes (associated to different time scales) have been shown, but will probably be hard to observe. Figure 7 (left and centre) shows the expected speckle image from a point source. The position/size of the small speckles (characterised by R_diff(λ)) are very sensitive to the wavelength, and it is clear that a desynchronisation of the maxima would be expected when observing such an idealised point source through different (narrow) passbands. But that real sources have a much larger projected radius than the speckle size completely screens this chromatic effect. The impact of the wavelength on the position/size of the wide (refractive) spots (Fig. 7-right) is much weaker (following λ^− 1/5 according to the combination of expressions (3) and (9)). Therefore, only a weak chromatic effect is expected from an extended source, even when observing with two very different passbands. For a clear signature of scintillation, it seems easier to take advantage of the rapidly varying luminosity with the observer’s position within the illumination pattern. As a consequence, simultaneous observations with two ~4 m class telescopes at a large separation (few 10³ km) would sample different regions of the illumination patterns (see Fig. 7 bottom right), and therefore measure (at least partially) decorrelated light curves as shown in Fig. 9 (right). The decorrelation will be complete if the distance between the two telescopes is greater than 2 × max(R_ref,R_s). A single observation of such a decorrelation will be sufficient to definitely confirm the discovery of a propagation effect that cannot be mimicked by an intrinsic variability.

7. Conclusions and perspectives

Through this work, we have simulated the phase delay induced by a turbulent refractive medium on the propagation of a wave front. We discussed the computational limitations of sampling the phase spectrum and to obtaining large enough illumination patterns. These limitations will be overriden in the near future with increasing computing capabilities. The illumination pattern on the observer’s plane has been computed for the promising strong regime of scintillation, and the effects of the source spatial and time coherencies have been included. We have established the connection between the modulation index (as an observable) of the illumination pattern with the geometrical parameters of the source and the strength of the turbulence (quantified by R_diff). Furthermore, we showed that when the spectral index of the turbulence increases (as in the case of supersonic turbulence), the detection of scintillating light curves should be easier.

These simulation studies and, more specifically, the modulation index topic were successfully used in our companion paper (Habibi et al. 2011) to interpret our light curve test observations of stars located behind known galactic nebulae and of stars from the Small Magellanic Cloud, in the search for hypothetical cold molecular halo clouds.

Time scales, such as t_ref = R_ref/V_T and t_s = R_s/V_T (where V_T is the relative velocity between the cloud and the line of sight), are observables that we plan to study further and in detail. Their extraction can be done through analysing the time power spectrum of the light curves, and it should give valuable information on the geometrical configuration, as well as on the turbulent medium.

The observing strategy has been refined through the use of the simulation, and we showed that observing desynchronised light curves simultaneously measured by two distant four-metre class telescopes would provide an unambiguous signature of scintillation as a propagation effect.

Acknowledgments

We thank J.-F. Lestrade, J.-F. Glicenstein, F. Cavalier, and P. Hello for their participation in preliminary discussions. We wish to thank our referee, Prof. B. J. Rickett, for his very careful review, which helped us to significantly improve the manuscript.

References

All Figures

	Fig. 1 Geometric configuration. The source is located in the (x2,y2) plane, the screen contains the refractive structure, and the observer is located in the (x0,y0) plane. A1(x1,y1) and ${\mathit{A}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{y}}_{\mathrm{1}}\mathrm{\right)}$ are the amplitudes before and after screen crossing.
In the text

Fig. 2 The phase-delay variations near the average for a simulated refractive screen with Nx × Ny = 20 000 × 20 000 pixels, Δ1 = 22.6 km, and Rdiff = 100 km. The grey scale ranges between ± 50 × 2π rad (clear regions correspond to an excess of phase with respect to the average). The zoom (inset in the lower-right corner) illustrates the self-similarity of the simulated screen (grey scale amplitude of 5 × 2π rad).

In the text

	Fig. 3 Phase structure functions Dφ(r) for a phase screen with Rdiff = 500 km. Blue line is the initial (theoretical) structure function. Red line is reconstructed from one of the realisations of the phase screen through simulation. The black curve is obtained from the numerical integration of the initial phase spectral density sampled as in the simulation.
In the text

Fig. 4 Typical illumination pattern from a point-like source. Here Rdiff = 100 km, the screen is at z0 = 160 pc, λ = 2.16 μm, then RF = 1300 km and Rref = 106 000 km. The typical length scale of the small-scale speckles is Rdiff, and the scale of the larger structures is Rref. The white square shows our fiducial zone with a margin of Lm = Rref/2 from the borders. Grey scale range from 0 to 4 times the mean intensity. The image has 20 000 × 20 000 pixels, each with a 22.6 km side.

In the text

Fig. 5 Typical illumination pattern from an extended source produced through the same screen as in Fig. 4, with source radius rs = 0.5 R⊙, located at z0 + z1 = 1 kpc + 160 pc (Rs ≃ 53 000 km). The small-scale speckles are smeared, and only the larger scale fluctuations survive. The white square shows the restricted fiducial zone with margin of Lm = Rref/2 + Rs from the borders. Grey scale ranges ± 20% around the mean intensity.

In the text

	Fig. 6 The smoothing function SF(x). L is the screen size, Lm is the margin from the screen borders.
In the text

Fig. 7 Simulated illumination maps (20 000 × 20 000 pixel of 22.6 km side) produced on Earth by a source located at z0 + z1 = 1.18 kpc through a refracting cloud assumed to be at z0 = 160 pc with a turbulence parameter Rdiff(2.16 μm) = 100 km. Here Rref(2.16 μm) ≃ 100 000 km. Top-left and middle: illumination produced at λ = 2.16 μm from a point-source with a zoomed detail; the contrast is 100%. The grey scale ranges from 0 to 4 times the mean intensity. Top-right: the same from a K0V star (rs = 0.85 R⊙, MV = 5.9, at 1.18 kpc V = 16.3). The circle shows the projection of the stellar disk (RS = rs × z0/z1). Here the modulation index is only 3.3%, and the grey scale ranges from ±20% around the mean intensity. The bottom maps are the illuminations in Ks wide band (λcentral = 2.162 μm, Δλ = 0.275 μm), using the same grey scales as above. The modulation index is 55% for the point-source (left and centre) and 3.3% for the extended source (right). The two parallel straight lines show the sections sampled by two observers located about 10 000 km apart, when the screen moves with the transverse velocity VT.

In the text

Fig. 8 The effective intensity modulation index for simulated scintillating stellar illumination patterns as a function of the theoretical modulation index. The blue dots show the effective modulation index values for different realisations of the phase screens. The red squares represent the mean value of the effective modulation indices with the same theoretical value. This plot shows the expected agreement for m < 0.15, a domain where the simulated screens are large enough not to suffer from statistical biases (see text).

In the text

Fig. 9 Light curves extracted along the 3 horizontal white lines for two illumination patterns. Left column: 2D pattern from a point-like source in Ks band with Rdiff = 300 km and Rref ≈ 28 000 km. The modulation indices of the three light curves differ by less than 5% from the 2D pattern modulation index. Right column: the illumination pattern for an extended source with Rs ≈ 41 000 km, through the same refractive screen. The modulation indices fluctuate by more than 30% around the 2D pattern index, implying the necessity of longer light curves for a better statistical representativity. The distance scale is common to both patterns. The circle shows the projected star disk. The 3 light curves from the right column are not completely decorrelated, because of their common proximity to the same large positive fluctuation.

In the text

	Fig. 10 Simulation of illumination scintillation patterns with associated light curves from an extended source assuming a refractive screen with non-Kolmogorov turbulences (left β = 3.1, right β = 3.9, see text). Here Rdiff = 100 km, Rref ≈ 8 × 104 km and the projected star radius is RS ~ 0.36Rref = 28 800 km.
In the text

Fig. 11 Top: two different phase screen power laws. Bottom: the corresponding modulation indices as a function of the expected modulation index for the Kolmogorov turbulence. The 3 indices plotted at a given abscissa correspond to screens with β = 3.1 (red), 3.67 (Kolmogorov, blue, along the diagonal), and 3.9 (black) with the same Rdiff. As the power law gets steeper, a larger modulation is expected. Blue dots show the Kolmogorov turbulence case.

In the text