Simulation of optical interstellar scintillation
^{1}
Laboratoire de l’Accélérateur Linéaire, IN2P3 − CNRS, Université de
ParisSud,
BP 34,
91898
Orsay Cedex,
France
email:
moniez@lal.in2p3.fr
^{2}
School of Astronomy, Institute for Research in Fundamental
Sciences (IPM), PO Box
193955531, Tehran,
Iran
^{3}
Department of Physics, Sharif University of
Technology, PO Box
113659161, Tehran,
Iran
^{4}
Perimeter Institute for Theoretical Physics, 31 Caroline Street North,
Waterloo, Ontario
N2L 2y5,
Canada
Received:
6
August
2012
Accepted:
6
February
2013
Aims. Stars twinkle because their light propagates through the atmosphere. The same phenomenon is expected on a longer time scale when the light of remote stars crosses an interstellar turbulent molecular cloud, but it has never been observed at optical wavelengths. The aim of the study described in this paper is to fully simulate the scintillation process, starting from the molecular cloud description as a fractal object, ending with the simulations of fluctuating stellar light curves.
Methods. Fast Fourier transforms are first used to simulate fractal clouds. Then, the illumination pattern resulting from the crossing of background star light through these refractive clouds is calculated from a Fresnel integral that also uses fast Fourier transform techniques. Regularisation procedure and computing limitations are discussed, along with the effect of spatial and temporal coherency (source size and wavelength passband).
Results. We quantify the expected modulation index of stellar light curves as a function of the turbulence strength – characterised by the diffraction radius R_{diff} – and the projected source size, introduce the timing aspects, and establish connections between the light curve observables and the refractive cloud. We extend our discussion to clouds with different structure functions from Kolmogorovtype turbulence.
Conclusions. Our study confirms that current telescopes of ~4 m with fastreadout, widefield detectors have the capability of discovering the first interstellar optical scintillation effects. We also show that this effect should be unambiguously distinguished from any other type of variability through the observation of desynchronised light curves, simultaneously measured by two distant telescopes.
Key words: dark matter / Galaxy: disk / Galaxy: halo / Galaxy: structure / local insterstellar matter / ISM: molecules
© 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 radioastronomy 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 nonKolmogorv 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 ESONTT 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 (x_{2},y_{2}) plane, the screen contains the refractive structure, and the observer is located in the (x_{0},y_{0}) plane. A_{1}(x_{1},y_{1}) and are the amplitudes before and after screen crossing. 

Open with DEXTER 
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_{2} + 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): (1)where x_{1} and y_{1} are the coordinates in the cloud’s plane, perpendicular to the sightline (see Fig. 1). Here φ(x_{1},y_{1}) is the phase delay induced to the wavefront after crossing the cloud, Nl(x_{1},y_{1}) is the cloud column density of H_{2} 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_{1},y_{1}). The statistical properties of the phase screen are described by the phase structure function D_{φ}(x_{1},y_{1}). By assuming an isotropic turbulence (Narayan 1992):
(2)where the first expression is averaged over the plane positions , , β 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 (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^{1}. Here we assume the gas to be a mixing of H_{2}/He with 24% He by mass (corresponding to the primordial abundances) and therefore ⟨ α ⟩ = 0.720 × 10^{30} m^{3}. In this expression, the cloud parameters are scaled to the values given by the PfennigerCombes 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: (4)where Fourier coordinates q_{x} and q_{y} have inverse length dimension, , and 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 (5)where I_{0}(x_{0},y_{0}) is the light intensity on the observer’s plane, L_{s} is the source luminosity, z_{1} is the sourcescreen distance (see Fig. 1), and h(x_{0},y_{0}) is given by the FresnelHuygens diffraction integral after considering the Fresnel and the stationary phase approximations (Born & Wolf 2002; Moniez 2003): (6)where R_{F} = is the Fresnel radius^{2} and z_{0} is the screenobserver distance. The Fresnel radius can be expressed as (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 Earthsource lineofsight, 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}: (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: (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 distancescale flux variations visible in Fig. 4. The corresponding time scale is given by t_{ref} = R_{ref}/V_{T}: (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 phasedelay variations near the average for a simulated refractive screen with N_{x} × N_{y} = 20 000 × 20 000 pixels, Δ_{1} = 22.6 km, and R_{diff} = 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 lowerright corner) illustrates the selfsimilarity of the simulated screen (grey scale amplitude of 5 × 2π rad). 

Open with DEXTER 
Numerical realisations of the 2D phase screen – made of N_{x} × N_{y} squared pixels of size Δ_{1} – 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 (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 vector associated to a pixel (i,j), we generate a random complex number where each component and spans the Gaussian distribution^{3} of zero mean, and 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})^{2} for the ensemble of (q_{x},q_{y}) vectors with the same module q equals L_{x}L_{y}S_{φ}(q) by construction. The phase screen φ(x_{1},y_{1}) 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 (and consequently ) 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 (12)where the average is performed on the (q_{x},q_{y}) coordinates^{4} spanning the circle of radius . The corresponding phase autocorrelation function is then given by Fourier transform:
(13)where J_{0} is the Bessel function. The recomputed structure function is then given by = 2(ξ^{rec}(0) − ξ^{rec}(r)), and the value of is deduced from its definition .
Fig. 3 Phase structure functions D_{φ}(r) for a phase screen with R_{diff} = 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. 

Open with DEXTER 
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 km since = 1. To find the origin of the difference with the input value R_{diff} = 500 km, we replaced 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 Δ_{1} = 28.85 km). Then we computed the integral (13) numerically with the same integration limits (q_{min},q_{max})^{5} 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 . 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 Δ_{1}) 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: (14)where (15)Here, (x_{1},y_{1}) and (x_{0},y_{0}) 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_{1},y_{1}) has been defined to avoid computational artefacts.
3.2.1. Screen regularisation
Since integral (14) is computed numerically, the coordinates (x_{1},y_{1}) 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): where L_{m} = 10R_{F} is the margin length from the borders of the screen with size L. We multiply the function G(x_{1},y_{1}) by SF(x_{1}) × SF(y_{1}) 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_{1},y_{1}) 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 pointsource 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 (halfsize) of the largescale 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_{1} = 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 km.
Fig. 4 Typical illumination pattern from a pointlike source. Here R_{diff} = 100 km, the screen is at z_{0} = 160 pc, λ = 2.16 μm, then R_{F} = 1300 km and R_{ref} = 106 000 km. The typical length scale of the smallscale speckles is R_{diff}, and the scale of the larger structures is R_{ref}. The white square shows our fiducial zone with a margin of L_{m} = R_{ref}/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. 

Open with DEXTER 
Fig. 5 Typical illumination pattern from an extended source produced through the same screen as in Fig. 4, with source radius r_{s} = 0.5 R_{⊙}, located at z_{0} + z_{1} = 1 kpc + 160 pc (R_{s} ≃ 53 000 km). The smallscale speckles are smeared, and only the larger scale fluctuations survive. The white square shows the restricted fiducial zone with margin of L_{m} = R_{ref}/2 + R_{s} from the borders. Grey scale ranges ± 20% around the mean intensity. 

Open with DEXTER 
Fig. 6 The smoothing function SF(x). L is the screen size, L_{m} is the margin from the screen borders. 

Open with DEXTER 
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_{1},y_{1}) contains frequencies that are higher than the Nyquist frequency f_{Nyq} = 1/(2Δ_{1}), where Δ_{1} 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_{1},y_{1}); it is at least necessary that Δ_{1} < R_{diff}/2 in order to sample phase variations up to 1/R_{diff} spatial frequency. R_{F} appears in the quadratic term ; the oscillation of this term accelerates as x_{1} and y_{1} increase, and aliasing occurs as soon as the distance between two consecutive peaks is smaller than 2Δ_{1}. In one dimension we therefore expect aliasing if (16)where x_{1}/Δ_{1} is the distance to the optical axis, expressed in pixels. The condition Δ_{1} = R_{F}/2 would be obviously insufficient here to avoid aliasing, since beyond ~11 pixels only from the optical axis (x_{1}/Δ_{1} > 11 pixel), the quadratic term would be undersampled. In practice, for the configurations considered in this paper, the Fresnel radius is >1000 km. When Δ_{1} = 22 km in the simulation, the aliasing starts at x_{1}/Δ_{1} > 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 pointlike source with the projected source limb profile (Moniez 2003; & Habibi 2011), (17)where L_{s} is the source luminosity, and the normalised limb profile is described as a uniform disk: where 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_{0} (160 pc) +z_{1} (1 kpc) = 1.16 kpc (R_{s} = 53 000 km). Highfrequency 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 z_{0} + z_{1} = 1.18 kpc through a refracting cloud assumed to be at z_{0} = 160 pc with a turbulence parameter R_{diff}(2.16 μm) = 100 km. Here R_{ref}(2.16 μm) ≃ 100 000 km. Topleft and middle: illumination produced at λ = 2.16 μm from a pointsource with a zoomed detail; the contrast is 100%. The grey scale ranges from 0 to 4 times the mean intensity. Topright: the same from a K0V star (r_{s} = 0.85 R_{⊙}, M_{V} = 5.9, at 1.18 kpc V = 16.3). The circle shows the projection of the stellar disk (R_{S} = r_{s} × z_{0}/z_{1}). 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 K_{s} wide band (λ_{central} = 2.162 μm, Δλ = 0.275 μm), using the same grey scales as above. The modulation index is 55% for the pointsource (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 V_{T}. 

Open with DEXTER 
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 finitewidth 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^{6}. 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. Twentyone illumination patterns were computed for 21 regularly spaced wavelengths within the [2.08, 2.28] μm interval, and were coadded to simulate the illumination pattern through the K_{s} passband. We checked that the spacing between successive wavelengths was small enough to produce a coadded 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 coadded 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 pointlike 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 (18)(Narayan 1992; Gwinn et al. 1998). This chromatic effect results from the high sensitivity of the constructive interference condition with the wavelength^{7}. In the case of Fig. 7, the decorrelation bandwidth is δλ_{dec}/λ ~ 1%. Since the K_{s} passband is δλ/λ ~ 0.1 ~ 10 × δλ_{dec}/λ, a firstorder 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 ~. 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): (19)As a consequence of the largescale smearing of the pointsource 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 Δ_{1} 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: (20)Because of this relation, N has to be as large as possible to minimize Δ_{1} 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 longduration 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). 

Open with DEXTER 
4.1. Modulation Index
For a pointlike 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_{0}/z_{1}) in the same regime, we always have m < 1, and Narayan (1992) showed that when the smallscale (diffractive) speckle is completely smeared, (21)in the case of the Kolmogorov turbulence. Here, θ_{ref} = R_{ref}/z_{0} is the angular refraction radius^{8} and θ_{s} the source angular radius. In this expression, Narayan (1992) assumed that z_{0} ≪ z_{1}, therefore z_{0} + z_{1} ≈ z_{1}. We use the following expression, which is formally identical to the previous one for z_{0} ≪ z_{1} but can also be used when z_{0} is not negligible: 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 pointsource pattern (see Fig. 7up) with the projected stellar profile (following expression (17)) makes the contrast decrease when the number of R_{ref}size domains within R_{S} = r_{S}z_{0}/z_{1} 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 , 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_{0} + z_{1} ≈ z_{1}, we can constrain R_{diff}(λ) from this relation (23)^{9}, even with poor knowledge of the screen’s distance z_{0}, 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_{0} ~ 80−190 pc and z_{1} = 8 kpc), and upper limits on hidden turbulent gas within the galactic halo (assuming z_{0} ~ 10 kpc and z_{1} + z_{0} = 62 kpc).
Fig. 9 Light curves extracted along the 3 horizontal white lines for two illumination patterns. Left column: 2D pattern from a pointlike source in K_{s} band with R_{diff} = 300 km and R_{ref} ≈ 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 R_{s} ≈ 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. 

Open with DEXTER 
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 pointlike 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^{5} 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 pointlike 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 pointlike source within the corresponding restricted fiducial zone (see Fig. 5), therefore they span ~2.5 × 10^{5} km and statistically include a little bit less than 10 R_{ref}scale variations. Because of this statistically short length, the lightcurvetolightcurve estimates of m_{extended} typically fluctuate by ~. 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_{0} = 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 nonKolmogorov turbulences (left β = 3.1, right β = 3.9, see text). Here R_{diff} = 100 km, R_{ref} ≈ 8 × 10^{4} km and the projected star radius is R_{S} ~ 0.36R_{ref} = 28 800 km. 

Open with DEXTER 
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 R_{diff}. As the power law gets steeper, a larger modulation is expected. Blue dots show the Kolmogorov turbulence case. 

Open with DEXTER 
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^{10}. The time power spectrum of the expected stochastic light curve should show a slope break around frequency , and possibly around (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 pointsource 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^{11}. 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^{4} 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 subminute 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 subminute 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. 7right) 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^{3} 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 fourmetre class telescopes would provide an unambiguous signature of scintillation as a propagation effect.
A cloud column of width L_{z} can include several turbulent structures with outer scale L_{out}. The direct relation of L_{out} with the turbulence strengh explains why R_{diff} increases in conjunction with this parameter. By contrast, since the column density increases with L_{z}, then the refraction also increases, thus decreasing R_{diff}.
For a given physical screen characterised by the column density Nl(x_{1},y_{1}), R_{diff} varies with λ^{6/5}, as shown in Eq. (3).
Nevertheless, one should consider that the gaseous screen is a nondispersive medium for optical wavelengths (dielectric medium with an index independent of λ), unlike the radioastronomy case (plasma), which may result in weaker chromaticity sensitivity. This specificity deserves more detailed studies, which are beyond the scope of this paper, considering the minor impact of the bandwidth compared to the smearing produced by the source size.
In Habibi et al. (2011), we used a simplified relation, with no significant impact on the resulting constraints.
The detection condition would be r_{s}/z_{1} ≲ 10 × R_{diff}/z_{0} (for which the projected stellar disk includes less than ~100 speckle spots), assuming a setup able to sample the target star with ≲1% photometric precision every few seconds. The scintillation of an A0 type star (R = 2.4 R_{⊙}) in the LMC (magnitude V = 19.4) seen through a screen with R_{diff} = 100 km at distance 30 pc is an example of such a favourable configuration, which could be discovered by the LSST (Abell 2009).
For the observed supersonic turbulence (with β > 11/3) see 10.
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
 Abell P. A., et al. (LSST Science Collaboration) 2009, LSST Science Book, submitted [arXiv:0912.0201] [Google Scholar]
 Born, M., & Wolf, E. 2002, Principles of Optics (Cambridge: University Press), 7th edn. [Google Scholar]
 Coles, Wm. A., Filice, J. P., Frehlich, R. G., & Yadlowsky, M. 1995, Appl. Opt., 34, 2089 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J., & Narayan, R. 1985, MNRAS, 214, 519 [NASA ADS] [Google Scholar]
 Goodman, J., & Narayan, R. 2006, ApJ, 636, 510 [NASA ADS] [CrossRef] [Google Scholar]
 Gwinn, C. R., Britton, M. C., Reynolds, J. E. et al. 1998, ApJ, 505, 928 [NASA ADS] [CrossRef] [Google Scholar]
 Habibi, F. 2011, Ph.D. Thesis, ParisSud and Sharif Universities [Google Scholar]
 Habibi, F., Ansari, R., Moniez, M., & Rahvar, S. 2011, A&A, 525, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamidouche, M., & Lestrade, J.F. 2007, A&A, 468, 193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Larson, R. B. 1981, MNRAS, 194, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Lovelace, R. V. E. 1970, Ph.D. Thesis, Cornell University [Google Scholar]
 Lyne, A. G., & GrahamSmith, F. 1998, Pulsar Astronomy (Cambridge University Press) [Google Scholar]
 McGaugh, S. S., Schombert De Block, W. J. G., & Zagursky, M. J. 2010, ApJ, 708, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Moniez, M. 2003, A&A, 412, 105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Narayan, R. 1992, Phil. Trans. Roy. Soc. Lond. A, 341, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Pfenniger, D., & Combes, F. 1994, A&A, 285, 94 [NASA ADS] [Google Scholar]
 Pfenniger, D., & Revaz, Y. 2005, A&A, 431, 511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rickett, B. J. 1986, ApJ, 307, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Rumsey, V. H. 1975, Radio Sci., 10, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Sieber, W. 1982, A&A, 113, 311 [NASA ADS] [Google Scholar]
 Tatarskii, V. I., & Zavorotnyi, V. U. 1980, Progr. Opt., 18, 207 [Google Scholar]
 Uscinski, B. J. 1977, The elements of wave propagation in random media (New York: McGrawHill) [Google Scholar]
All Figures
Fig. 1 Geometric configuration. The source is located in the (x_{2},y_{2}) plane, the screen contains the refractive structure, and the observer is located in the (x_{0},y_{0}) plane. A_{1}(x_{1},y_{1}) and are the amplitudes before and after screen crossing. 

Open with DEXTER  
In the text 
Fig. 2 The phasedelay variations near the average for a simulated refractive screen with N_{x} × N_{y} = 20 000 × 20 000 pixels, Δ_{1} = 22.6 km, and R_{diff} = 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 lowerright corner) illustrates the selfsimilarity of the simulated screen (grey scale amplitude of 5 × 2π rad). 

Open with DEXTER  
In the text 
Fig. 3 Phase structure functions D_{φ}(r) for a phase screen with R_{diff} = 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. 

Open with DEXTER  
In the text 
Fig. 4 Typical illumination pattern from a pointlike source. Here R_{diff} = 100 km, the screen is at z_{0} = 160 pc, λ = 2.16 μm, then R_{F} = 1300 km and R_{ref} = 106 000 km. The typical length scale of the smallscale speckles is R_{diff}, and the scale of the larger structures is R_{ref}. The white square shows our fiducial zone with a margin of L_{m} = R_{ref}/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. 

Open with DEXTER  
In the text 
Fig. 5 Typical illumination pattern from an extended source produced through the same screen as in Fig. 4, with source radius r_{s} = 0.5 R_{⊙}, located at z_{0} + z_{1} = 1 kpc + 160 pc (R_{s} ≃ 53 000 km). The smallscale speckles are smeared, and only the larger scale fluctuations survive. The white square shows the restricted fiducial zone with margin of L_{m} = R_{ref}/2 + R_{s} from the borders. Grey scale ranges ± 20% around the mean intensity. 

Open with DEXTER  
In the text 
Fig. 6 The smoothing function SF(x). L is the screen size, L_{m} is the margin from the screen borders. 

Open with DEXTER  
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 z_{0} + z_{1} = 1.18 kpc through a refracting cloud assumed to be at z_{0} = 160 pc with a turbulence parameter R_{diff}(2.16 μm) = 100 km. Here R_{ref}(2.16 μm) ≃ 100 000 km. Topleft and middle: illumination produced at λ = 2.16 μm from a pointsource with a zoomed detail; the contrast is 100%. The grey scale ranges from 0 to 4 times the mean intensity. Topright: the same from a K0V star (r_{s} = 0.85 R_{⊙}, M_{V} = 5.9, at 1.18 kpc V = 16.3). The circle shows the projection of the stellar disk (R_{S} = r_{s} × z_{0}/z_{1}). 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 K_{s} wide band (λ_{central} = 2.162 μm, Δλ = 0.275 μm), using the same grey scales as above. The modulation index is 55% for the pointsource (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 V_{T}. 

Open with DEXTER  
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). 

Open with DEXTER  
In the text 
Fig. 9 Light curves extracted along the 3 horizontal white lines for two illumination patterns. Left column: 2D pattern from a pointlike source in K_{s} band with R_{diff} = 300 km and R_{ref} ≈ 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 R_{s} ≈ 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. 

Open with DEXTER  
In the text 
Fig. 10 Simulation of illumination scintillation patterns with associated light curves from an extended source assuming a refractive screen with nonKolmogorov turbulences (left β = 3.1, right β = 3.9, see text). Here R_{diff} = 100 km, R_{ref} ≈ 8 × 10^{4} km and the projected star radius is R_{S} ~ 0.36R_{ref} = 28 800 km. 

Open with DEXTER  
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 R_{diff}. As the power law gets steeper, a larger modulation is expected. Blue dots show the Kolmogorov turbulence case. 

Open with DEXTER  
In the text 