E Siminos

Stefan Skupin

A Sävert

Jason M Cole

S P D Mangles

M C Kaluza

Modeling ultrafast shadowgraphy in laser-plasma interaction experiments

Keywords: numbers: 42.30.Kq, 42.30.Va, 52.38.Kd, 41.75.Jv, 52.65.Rr, 52.70.-m

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The interaction of intense laser pulses with plasma can lead to the excitation of nonlinear structures, such as wakes, solitons, shocks and filaments [1]. The time evolution of such structures is highly non-trivial and a lot of effort has been invested in developing theoretical tools and experimental diagnostics to investigate it. When such nonlinear structures move close to the speed of light c, as is the case for a wake driven by an intense laser pulse in an underdense plasma, additional difficulties arise due to the requirement of short interaction times and micrometer length-scales in order to resolve the relevant dynamics.

In a typical laser-wakefield acceleration (LWFA) experiment, a relativistically intense laser pulse (the 'pump' pulse) excites an electron plasma wave (the 'wake') which propagates with a phase velocity v p equal to the pump pulse's group velocity v g c through the plasma. Ambient electrons can be trapped and accelerated in the strong electromagnetic fields of the wake, producing quasi-monoenergetic beams [2][3][4] up to an energy of 4 GeV [5]. Stable operation of such accelerators requires a thorough understanding of nonlinear laser-plasma interaction processes motivating the development of novel diagnostics. The characteristic length scale of the wake is the plasma wavelength λ p . Therefore, optical prob- * evangelos.siminos@gmail.com ing of such a wake with a transversely propagating probe pulse requires a probe pulse duration τ pr λ p /c. Typically this corresponds to a duration τ pr of a few optical cycles τ L = 2π/ω L , where ω L is the center frequency of the probe pulse. Recent pump-probe experiments utilized probe pulses of few femtosecond duration to obtain the first direct shadowgraphic images of a laserdriven wake propagating in an underdense plasma [6,7]. These measurements offer far more information about the plasma wave than previous ones using longitudinal [8][9][10] or longer transverse probe pulses [11,12].

Conventional shadowgraphy can be analyzed using geometrical optics and be formulated as an inverse problem: for a given shadowgraphic image local gradients in index of refraction can be retrieved, see, e.g., Ref. [START_REF] Panigrahi | Schlieren and Shadowgraph Methods in Heat and Mass Transfer[END_REF]. In ultrafast shadowgraphy of laser induced wakes such retrieval is in general not straightforward. The effect of longitudinal motion of the wake as the probe pulse traverses, the presence of strong magnetic fields as well as relativistic effects are not negligible and need to be taken into account. At the same time diffraction effects can be significant when the probe pulse wavelength λ pr is comparable to the characteristic length scale of density gradients in the plasma, i.e., the skin depth λ s = c/ω pe , where ω pe = q 2 e n 0 / 0 m e 1/2 is the electron plasma frequency for a plasma of density n 0 . For parameters used in recent experiments [7] this is the case and the effect of diffraction has to be accounted for.

In this paper we analyse how ultra fast probe pulses form shadowgraphic images by using 3D particle-in-cell arXiv:1509.08846v1 [physics.plasm-ph] 29 Sep 2015 (PIC) simulations of the full pump-probe interaction. Propagation of the pump pulse through the plasma is fully simulated and at different delays a probe pulse propagating transversely to the pump is launched. It is worth noticing that the computational overhead due to the probe pulse is very small, because only (transverse) propagation through the interaction zone (plasma wake) needs to be simulated by the PIC code. After the probe traverses the wake, post-processing in Fourier space allows the effect of a typical imaging setup in shadowgraphic image formation to be taken into account. We show that synthetic shadowgrams generated using this methodology turn out to be crucial for the correct interpretation of experimental measurements.

This paper is structured as follows. In Sec. II we present our PIC simulations of pump and probe propagation. The methodology we use to post-process the results in order to obtain artificial shadowgrams is explained in Sec. III. The direct confrontation with experimental shadowgrams in Sec. IV illustrates the importance of synthetic shadowgrams in the interpretation of experimental results. In Sec. V we present several case studies illustrating the effect of varying object plane position, duration and wavelength of probe pulse on the shadowgrams and show that time-dependent information can be recovered from a single chirped probe pulse by spectral filtering. Finally, in Sec. VI we discuss our findings and present our conclusions.

II. PIC SIMULATIONS

A. Pump propagation

We perform simulations of LWFA using the PIC code EPOCH [14] with parameters similar to the ones used in recent experiments that utilized ultrafast shadowgraphy to study injection in LWFA [7]. We consider a plasma with super-Gaussian electron density profile n e = n 0 e -(x-x0) 4 /w 4 p in propagation direction of the pump pulse, where n 0 = 1.7 × 10 19 cm -3 , x 0 = 1.17 mm and w p = 1.198 mm are taken from a fit of a typical experimental density profile of Ref. [7]. The computational domain is a 'sliding window' of size 150 × 70 × 70 µm 3 moving at c. We use 2700 × 525 × 525 cells and two particles per cell, with fifth order particle weighting to reduce noise. Ions are considered as a stationary backgound. A Gaussian pump laser pulse propagating along the x direction, linearly polarized along the y direction, with an intensity full width at half maximum (FWHM) duration of 36 fs and a central wavelength λ L = 810 nm is focused to a spot size of 18.84 µm (intensity FWHM) at x f = 300 µm. The pump pulse's maximum intensity (in vacuum) is I 0 = 2.5 × 10 18 W/cm 2 . To ensure better dispersion properties of both the pump and probe pulses, while keeping the resolution requirements low, the sixth order (in space) finite-difference-time-domain (FDTD) Maxwell solver available in EPOCH is used.

B. Probe propagation

Probe propagation in the plasma is also fully simulated in 3D with EPOCH. While simulating pump propagation, so called 'restart dump' files, which contain the complete data representing the present state of the PIC simulation are saved at regular intervals. This allows us to start several probe simulations at different 'delays' using a single simulation of pump propagation. For each probe simulation the moving window is stopped and a probe pulse is injected from the side of the box, propagating along the negative y-direction, i.e. perpendicularly to the pump propagation direction, see Fig. 1. Because the pump propagates during the probing stage (without a moving window) additional space along the x-direction is required. Also, depending on the probe pulse duration τ pr , enlarging the computational box in y-direction may be necessary. Moreover, it is important to take into account diffraction of the probe pulse from the computational boundaries along the x and z directions. To keep it to a minimum we use a flat-top spatial profile of the probe pulse with sinusoidal ramps of 10 µm length. In our reconstruction of the shadowgrams we only take into account the part of the computational domain that is not affected by numerical diffraction from the boundaries. All this results in a larger computational domain compared to a typical LWFA simulation for the same parameters, however, this larger domain is only necessary during the probe simulations. Moreover, each probe simulation is still much smaller than the full pump simulation because the pump propagates for millimetres whereas the probe only propagates for tens of microns. Although we investigate the influence of various probe pulse parameters, see Sec. V, in all other runs the pulse parameters are similar to the ones in recent experiments [7]: central wavelength λ pr = 0.75 µm, bandwidth limited duration of τ BL = 4.4 fs (intensity FWHM), a negative linear chirp that increased its duration to τ pr = 12 fs, and maximum intensity I pr = 8.6 × 10 14 W/cm 2 . The transverse size of the plasma is limited to few tens of microns in both the experiments and our simulations, and for such short propagation distances the use of sixth order (FDTD) solver guarantees probe propagation with-out significant lengthening of the pulse due to numerical dispersion for relatively low spatial resolution (see, e.g., Ref. [START_REF] Shlager | [END_REF] for a study of dispersion characteristics of high order FDTD solvers). We note that our probe pulse intensity is four orders of magnitude higher than typical experimental values [7] in order to have a better signal-tonoise ratio in the simulated shadowgrams (in PIC simulations the noise level in the electromagnetic fields is much higher than in experiments due to the discretization by macroparticles). However, we verified that varying the probe intensity by one order of magnitude does not alter the shadowgrams, i.e., the results can be considered linear in the probe intensity (as one would indeed expect for non-relativistic probe pulses). We note that one has to be careful when using increased probe intensity if ionization effects are included in the simulations (neglected in the present study).

In the simulations we allow the probe pulse to propagate past the wake structure, here approximately until its center reaches y S = -20µm. In Fig. 2(a-c) we track the evolution of the envelope of the probe electric field (absolute of the complex electric field, see Sec. III A) as it crosses the wake. For presentation purposes only, we subtract the background fields using a simulation with identical initial conditions but no probe pulse. Figure 2(b-c) indicates that modulations in the envelope of the probe pulse occur as it interacts with density gradients in the wake. The highest envelope modulations occur at the front of the bubble where the pump drives the wake.

Once the probe pulse reaches y S there are no further density perturbations in the plasma, and therefore all local phase differences induced by the wake have been imprinted to the probe pulse. However, if we simply try to reconstruct a shadowgram by recording the timeintegrated intensity passing through the plane y S , as shown Fig. 2(d), there are two problems with the image we obtain. In the front of the wake there is strong scattering of pump light which is of much higher intensity than the probe intensity. This scattered light is not present in the experimental shadowgrams [6,7]. As we will see the reason for this is that the aperture of the imaging system eliminates most of the side-scattered light. Moreover, the wake structure appears blurred in Fig. 2(d). The reason for this is strong diffraction, since the length scale for density gradients (λ s = 1.3 µm here) is comparable to the wavelength of the probe pulse (λ pr = 0.75µm).

In the experiments diffraction is compensated by optical imaging and thus, in order to be able to compare our PIC simulation results to experimental shadowgrams, we have to take into account the role of a typical imaging system in shadowgram formation. This could be done by solving a vacuum diffraction integration through the entire optical system. However, as we will describe in Sec. III, connecting the near field results of our PIC simulations to the Fourier description of optical systems offers a more straightforward and easy-to-implement method.

III. POST-PROCESSING AND RECONSTRUCTION OF SHADOWGRAMS

In order to reconstruct shadowgrams from our PIC simulation results, we use two major simplifications.

Firstly, the transverse size of the simulation box in our PIC simulations is limited to a few of tens of microns, and so the PIC simulation cannot directly account for the optical imaging system. However, this is not really necessary, because once the probe pulse has traversed the plasma wake, all information related to the laserplasma interaction process are already imprinted in its wavefronts. At the same time, the wavefronts are distorted by diffraction, but this will be compensated by the imaging system. In order to take this into account, we assume that the probe pulse propagates further on in vacuum, and neglect the influence of any low density plasma and/or gas which may be present between the main interaction region and the imaging system. By doing so, we can use much more efficient spectral pulse propagation methods. Secondly, the probe pulse, after it has passed the plasma wake, propagates mainly in one direction, and the imaging system has a certain numerical aperture (NA), i.e., acceptance angle. Here, we assume that the NA of our imaging system is small enough to justify a paraxial description of the imaging process.

A. Post-processing: Theory An optical imaging system in the paraxial approximation can be conveniently described by means of Fourier optics. Here, we will restrict ourselves to a generic 4f imaging system, which can be found in any textbook on optics, see, e.g., Ref. [START_REF] Saleh | Fundamentals of Photonics[END_REF]. However the results are broadly equivalent to a more general optical system. Figure 3 shows a schematic setup. When applying the thin element approximation for the two lenses, one can relate any transverse electric (or magnetic) field component û in frequency space (see Appendix A for definitions) in the object plane (y = y o ) and the image plane (y = y i ) via

û(x, y i , z, ω) = -e i4f ω/c × k ⊥ < ωD 2cf ū(k x , y o , k z , ω)e -ikxx-ikzz dk x dk z . (1)
Here, D is the diameter of the aperture ( = diameter of the mask in the Fourier plane), f the focal length of the two identical lenses, and ū(k x , y o , k z , ω) represents the field component in the object plane. In our configuration (propagation of the probe pulse in negative y direction, û = Êx , Bz , Êz , or Bx ), we have y i = y o -4f , and the aperture blocks all Fourier components with transverse wave vector 1) is derived by assuming paraxial propagation in vacuum, and the two (perfect) lenses and the mask in the Fourier plane are treated in the thin element approximation.

k ⊥ = k 2 x + k 2 z > ωD 2cf . Equation (
From our PIC simulations, we have access to a temporal snapshot of the fields E( r, t S ), B( r, t S ), at some time t = t S when the probe pulse has already passed through the interaction region and its center is at y S , i.e., at a distance ∆y = y Sy o from the object plane. So in general, at time t = t S the probe pulse has already passed the object plane, which usually lies somewhere inside the interaction region, but has not yet reached the first lens. By symmetry of the 4f -system, the electromagnetic fields at y S are the same as the fields at y S = y i + ∆y at a distance ∆y from the image plane [up to a constant phase shift and image reversal, cf. Eq. ( 1)]. The fields at the image plane could then be obtained by propagating the fields from point y S to the image plane or, by propagating the fields from position y S backwards to the object plane. We choose to do the latter, thus circumventing the need to describe the full propagation of the pulse through the imaging system. In our argument we ignored the effect of the aperture. However, since its role is to merely eliminate directions for which k ⊥ < ωD 2cf , it can be taken into account independently of the position of the pulse, a fact that we will use in the following.

The (backwards) propagation from y S to the object plane is achieved by using the temporal snapshot of the electromagnetic field as initial datum to Maxwell's equations in vacuum. To this end, it is useful to remember that any solution to Maxwell's equations in vacuum E( r, t), B( r, t) can be decomposed into plane waves

E k ( r, t) = E k,0 e i k• r-iωt (2) B k ( r, t) = k × E k,0 ω e i k• r-iωt , ( 3
)
where k is the wave vector, E k,0 is the amplitude vector, k • E k,0 = 0, and the vacuum dispersion relation dictates

ω = kc with k = | k| = k 2 x + k 2 y + k 2 z .
For technical convenience, in the following we resort to complex field notation

E( r, t) = E( r, t) + E * ( r, t) = 2Re E( r, t) (4) B( r, t) = B( r, t) + B * ( r, t) = 2Re B( r, t) (5) E( r, t) = Θ(ω) ˆ E( r, ω)e -iωt dω (6) B( r, t) = Θ(ω) ˆ B( r, ω)e -iωt dω, (7)
and Θ(x) represents the usual Heaviside step function.

For given E( r, t S ), B( r, t S ) we can thus compute the fields at any time t and at any position (in particular in the object plane) via

E( r, t) = ¯ E( k, t S )e -ikc(t-t S ) e i k• r d 3 k (8) B( r, t) = ¯ B( k, t S )e -ikc(t-t S ) e i k• r d 3 k, (9)
where ¯ E( k, t S ), ¯ B( k, t S ) denote the respective complex fields in spatial Fourier space.

Even though pump and probe pulse may overlap in position space at time t S of our snapshot, they are well separated in Fourier space due to their different propagation directions. Our 4f system selects plane wave components inside a cone defined by k ⊥ < kD 2f , and therefore only the probe pulse gets imaged (see Fig. 4). Moreover, it is safe to assume that no light propagates in positive y direction, i.e., in opposite direction to the probe pulse. This additional assumption allows us to write (see Appendix B for the general case)

¯ E( k, t S ) = Θ(-k y ) ¯ E( k, t S ) (10) ¯ B( k, t S ) = Θ(-k y ) ¯ B( k, t S ). (11)
Hence, by using Eq. ( 1), the complex fields of the probe pulse in the image plane at any time t are given by

E(x, y i , z, t) = - k ⊥ < kD 2f ¯ E( k, t S )e -ik[c(t-t S )-4f ] e -ikxx+ikyyo-ikzz d 3 k. ( 12
) B(x, y i , z, t) = k ⊥ < kD 2f ¯ B( k, t S )e -ik[c(t-t S )-4f ] e -ikxx+ikyyo-ikzz d 3 k. (13)
Finally, the shadowgram can be constructed by computing the y-component of the time integrated Poynting vector in the image plane

S y (x, z) ∝ e y • +∞ -∞ E(x, y i , z, t)× B(x, y i , z, t)dt. ( 14
)
We note that the propagation in time is transformed to propagation in space in the practical implementation of our post-processing methodology, as explained in Sec. III B.

FIG. 4. The 4f system selects plane wave components in a cone defined by k ⊥ < kD 2f only. Because the probe pulse propagates in negative y direction, its complex electric field in Fourier space Ē( k, t) is nonzero for ky < 0 only. The same is true for the complex magnetic field B( k, t).

B. Post-processing: Implementation

We have seen in the previous section that in order to simulate shadowgrams, we have to 1. filter the relevant Fourier components of the snapshot E( r, t S ), B( r, t S ) from the PIC simulation in order to extract the probe electromagnetic field components which will form the image, 2. propagate the probe fields (in vacuum) to the object plane (resp. image plane), and compute the time integrated Poynting vector to obtain the shadowgram.

Let us now develop a post-processing algorithm doing exactly that. From our PIC code, once the probe pulse has passed the interaction region, we record the real field components

E (i,j,m) α , B (i,j,m) α
, α = x, z for each point (i, j, k) on our computational grid, i = 0, . . . , N x -1, j = 0, . . . , N y -1, m = 0, . . . , N z -1. For ease of implementation N x , N y , N z are assumed to be even (if not, we truncate the matrices accordingly). We cell-center the data which originally reside on a Yee grid. We perform a three dimensional spatial discrete Fourier transform of the field components in order to obtain Ē(kx,ky,kz)

α , B(kx,ky,kz) α , where k x ∈ 2π
Nx ∆x [-N x /2, -N x /2 + 1, . . . , N x /2 -1], etc. In the following, wavenumbers will be used as indices of matrices storing the spatial Fourier transform of fields.

The first step is to set all matrix elements of matrices Ēα , Bα for which k ⊥ > kD/2f holds equal to zero. This models the effect of the aperture of the imaging system on our data, and also eliminates the field of the pump pulse (see Fig. 5). Note that because we computed the spatial Fourier transform of real fields, both cones (positive and negative k y ) feature nonzero Fourier amplitudes (cf. Fig. 4). Because typical values of D/f are small, most entries of the filtered matrices are zero. Thus, copying only the relevant matrix elements to new matrices Ē α , B α of dimensions N x , N y , N z significantly increases the efficiency of all further post-processing. The second step is to propagate (in vacuum) the components of the probe electromagnetic field obtained above to the position of the object plane at y = y o and compute the time integrated Poynting vector, i.e., the shadowgram. According to Eqs. ( 8)-( 9), propagation by distance ∆y = y o -y is achieved by multiplying components Ē (kx,ky,kz) α , B (kx,ky,kz) α with P (∆y) = e -i sign(ky) k ∆y .

(

) 15
Note that here we translated propagation time into propagation distance via c∆t = ∆y, and the sign factor sign(k y ) in the exponent accounts for the fact that we deal with Fourier transforms of real fields, and not complex fields as in Eqs. ( 8)- (9). Since in experiments the position of the object plane can vary from shot to shot by a few microns with respect to the position of the center of the wake, it is useful to compute shadowgrams for different positions of the object plane y o using the following approach: We introduce a matrix S y (i,j,m)

of the same dimensions as E α to store time-integrated

Poynting flux along y. The fields are then propagated for ∆y ∈ [∆y (min) , ∆y (min) + δy, . . . , ∆y (max) ] where δy λ pr . Typically, one may use the spatial grid vector elements y s = y j for j < j max , where j max should be chosen such that the probe pulse does not reach the boundaries of the numerical box. Then, for each propagation distance y s , we transform the propagated fields back to position space, form e y • E × B and sum up all contributions in S y . By construction, slices of the matrix S y along the direction j (corresponding to y) are proportional to the time-integrated Poynting flux for an object plane at this position. Of course, out of the N y possible object plane positions, we can only use those through which the full length of the probe pulse was propagated during the above procedure.

A typical simulated shadowgram obtained by this procedure is shown in Fig. 2(e) for y o = 0. We can see that the scattered light from the pump which overwhelmed the image in Fig. 2(d), has been eliminated through the introduction of an aperture and subsequent suppression of large transverse wavevector components, and thus at the expense of spatial resolution. At the same time, as demonstrated in Sec. IV, correcting for the object plane position has allowed to obtain an image that matches the typical experimental shadowgrams [7], where special care is taken to image approximately the center of the interaction region.

C. Post-processing: Frequency Filtering

In some cases it is desirable to account for frequency filters or polarizers in the experimental setup, or even the sensitivity curve of the CCD sensor itself. Such filters can be simulated easily after step one of the above postprocessing algorithm, i.e., in spatial Fourier space. We just have to translate the frequency filter function given in ω to wavenumbers k via the vacuum dispersion relation k = ω/c, and multiply all field components. However, special care must be taken when the frequency filter is such that it would lead to an effective increase of the duration of the probe pulse, for example a narrow bandpass filter. In such case, we may have to increase the y range of the computational box to accommodate the filtered pulse, i.e., add some zero elements to the respective field matrices E α , B α . Without these additional measures, we can have the situation that the filtered pulse overlaps with itself by virtue of the periodic boundary conditions implied by the use of discrete Fourier transformations.

IV. COMPARISON WITH EXPERIMENTS: AN EXAMPLE

The ability to track the probe propagation in our simulations, as in Fig. 2(a-c), allows to interpret shadowgraphic images, such as those of Fig. 2(e), and associate them with transverse (to the direction of probe propagation) density gradients in the plasma. As we will show in this section, this is important in order to interpret experimental shadowgraphic images and to deduce quantitative information from them.

In Fig. 6(a-b) we plot cross-sections of the density distribution from a snapshot of our PIC simulation, while the corresponding shadowgram is shown in Fig. 6(d). As Fig. 2 suggests, the main contribution is due to the strong density gradients near the symmetry axis of the wake. Moreover, the strongest modulations originate at the front of the bubble, where the laser-pulse drives the wake.

In Fig. 7 experimental shadowgrams (bottom row) taken with the few-cycle probing setup available at the JETI-laser facility at the Institute of Optics and Quantum Electronics in Jena, Germany, are shown. In the middle row we show computed shadowgrams obtained using the same probe characteristics as in Fig. 6, but also taking into account the sensitivity versus wavelength of the CCD sensor used in the experiments. Comparison of the experimental shadowgrams with the computed shadowgrams reveals that essential quantities like the period of the plasma oscillation, which can be easily retrieved in the electron density profile, are retained in the shadowgram. In Fig. 7 the length of the second oscillation period is defined by the distance between the dark regions in the intensity profile. In the PIC simulations the front of the bubble is readily identified as a peak in electron density. However the experimental shadowgrams exhibit strong intensity modulations in this region which are not readily recognisable as this density peak. The calculated shadowgrams also exhibit these modulations and hence allow us to confirm that the centre of these modulations directly corresponds to the front of the bubble. Crucially this then allows direct and quantitative measurements of the bubble shape and size from the experimental shadowgrams. Exemplary shots showing different stages during the evolution of the laser driven wake are plotted in Fig. 8. In the first stage, in front of the position of the vacuum focus, only a small modulation in the plasma is visible. The phase fronts of the wave are slightly tilted due to a spatio-temporal asymmetry in the pump pulse. Later the pump pulse drives a plasma wave with a high amplitude. The transverse extent is of the order of the focal spot size and the phase fronts show the characteristic curvature associated to nonlinear wakes [START_REF] Bulanov | [END_REF]. After wavebreaking and electron injection, the shape of the plasma wave, especially the back of the bubble, changes drastically. At the same time a strong modulation at the bubble front appears, which is also visible in the simulations.

We observe that the largest density gradient in Fig. 6(a-b) is associated with the injected electron bunch inside the bubble. However, this bunch is not visible in the shadowgrams. The reason is that the index of refraction for a relativistic plasma is: where ω pr is the probe frequency, γ = 1 + p 2 /m 2 c 2 1/2 and p is the cycle-averaged electron fluid momentum. The electrons of the injected bunch are therefore energetic enough to be invisible to the probe pulse, as can be seen in Fig. 6(d). This can be seen as another facet of relativistic transparency [18][19][20][21][22], where relativistic motion of electrons is induced by the wake field, while the modified index of refraction is experienced by the probe pulse. We note that no signature of injected electrons could be detected in the experimental shadowgrams of Ref. [7], even when injection occurred before the shadowgram was taken, as could be seen by the emission of wavebreaking radiation [23]. Such an example experimental shadowgram is shown in Fig. 9. The simultaneously measured electron spectrum (not shown here) exhibits a broad energy range from 60 MeV to 100 MeV. At this time of the snapshot the laser pulse has propagated 1130 µm into the gas jet. The bright spot on the left is attributed to wavebreaking radiation [23]. This is a direct signature for the onset of self-injection in the experiment. Relative to this point the laser pulse has traveled 230 µm and electrons should be accelerated to an energy of ≈ 40 MeV (γ = 80) assuming an effective electrical field of 160 GV/m. The typical charge measured by a spectrometer for the conditions of this experiment is 40 pC. Assuming that the bunch is ellipsoidal with volume 2 µm 3 we find a number density n b = 1.25 × 10 20 cm -3 . This implies no significant difference in the refractive index compared to the electron cavity surrounding the electron bunch. Thus the bunch remains invisible in the shadowgram.

η = 1 -ω 2 p /(γ ω 2 pr ) 1/2 ,
v g t = 1450µm v g t λ b λ p λ b λ p a) b) λ b λ p λ b λ p λ b λ p λ b λ p ξ = x-ct (µm) ξ = x-ct (µm) x-vgt (µm) ξ = x-ct (µm) ξ = x-ct (µm) x-vgt (µm)

V. EFFECT OF PROBE PULSE AND IMAGING SYSTEM PARAMETERS

In this section we investigate changes in the simulated shadowgrams when parameters of the probe pulse or the imaging systems are varied, which could be taken into account in the design of future experiments. As an important application we show that chirped pulses together with spectral filtering can, in principle, be used to obtain time-dependent information from a single shot.

A. Effect of object plane position

The post-processing methodology described in Sec. III allows us to adjust the object plane position for the same probe propagation simulation. The effect of object plane position on the produced shadowgram is shown in Fig. 10, where two shadowgrams have been reconstructed for different object plane positions y o = -10µm and y o = 10µm from the same pump-probe simulation used in Fig. 6. Changing object plane position not only brings different parts of the wake into focus, but also results in enhanced blur due to diffraction compared to the case y o = 0 shown in Fig. 6(d).

B. Effect of probe pulse wavelength

Using a shorter wavelength probe pulse would be expected to reduce diffraction and improve resolution. In Fig. 11(a) and (b) we see that this is indeed the case when λ pr is reduced to 0.6µm and 0.45µm, respectively. At the same time, since the refraction index of the plasma depends on the ratio of density to the critical density (for the probe frequency) the variation in the index of refraction becomes smaller for shorter wavelength, leading to somewhat lower contrast in the obtained images.

C. Effect of aperture

The presence of an aperture in the imaging system prevents light scattered at large angles from reaching the image plane, while at the same time it limits the spatial resolution of the shadowgrams. We illustrate the effect of varying the aperture acceptance angle in Fig. 12, where N A = 0.35. Compared to Fig. 6(d) for which N A = 0.26 (matching the experimental one) we see that in the present case the resolution is somewhat higher at the expense of increased noise due to scattered light.

D. Effect of probe pulse chirp

Introducing a linear chirp in the probe pulse increases its duration (for a given spectral width) and this is expected to result in a reduction in resolution of the shadowgrams due to motion blur. In our simulations a linear chirp is introduced through a quadratic phase in the probe pulse temporal profile,

φ = C(t -t 0 ) 2 /w 2 t , (16)
where C = ± (τ pr /τ BL ) 2 -1 1/2 is the chirp parameter, τ BL is the bandwidth limited probe pulse duration, t 0 corresponds to the pulse maximum, and w t is connected to the intensity FWHM duration of the probe pulse by w t = τ pr / √ 2 ln 2 and the positive (negative) sign of C corresponds to a positive (negative) chirp. This is illustrated in Fig. 13. We see that the probe pulse of Fig. 13(a) with no chirp, i.e., having the bandwidth-limited duration of τ pr = τ BL = 4.4 fs results in a shadowgram of higher resolution than the one of Fig. 6(d), which was obtained with a pulse that had the same bandwidth but a duration τ pr = 12 fs (due to a negative chirp). Introducing a negative chirp such that τ pr = 30 fs results in the shadowgraphic image of Fig. 13(b) in which the structure of the wake is no longer visible. Since the length of such a pulse is c τ pr 9 µm, which is longer than λ p = 8.1 µm, the wake structure cannot be resolved in this case.

E. Single shot movie generation through spectral filtering

The achievable temporal resolution for a linearly chirped pulse is given by the bandwidth limited pulse duration [24], and the information contained in the chirped 30 fs pulse of Fig. 13(b) can be recovered by spectral filtering. Introducing a step-like bandpass filter with central wavelength λ c in the range 0.6 to 0.9 µm and width of 0.02 µm allows us to recover the shadowgrams of Fig. 14. In these shadowgrams the signal to noise ratio is increased compared to the ones taken with shorter pulses, since the fraction of the pulse energy contributing to the image formation is reduced. Since the probe pulse has a negative chirp, shorter wavelengths interact with the wake earlier than longer wavelengths and each frequency-filtered shadowgram captures a different snapshot of the motion of the wake in time. This indicates that 1D time evolution of the wake could be obtained from a single pump-probe shot using a spectrometer after the interaction. More complicated setups involving splitting and filtering the chirped probe pulse after the interaction could in principle allow the recovery of 2D timedependent information. For example, we can estimate the wake phase velocity v p by measuring the position of a certain intensity maximum in the different images of Fig. 14. In general one should take into account the plasma dispersion relation when analysing probe propagation, but here we will assume that the probe pulse chirp does not change significantly over the course of propagation in the plasma and that we can use the vaccuum dispersion relation to connect frequency to wavelength. This is justified for our simulations since the phase velocity of light in a tenuous plasma is given by v p c(1 + n 0 /n c ), where n c = 0 m e ω 2 /e 2 is the critical density for frequency ω. For the range of frequencies that we consider here, 5 × 10 -3 n 0 /n c 10 -2 and the probe propagates in plasma for only a few tens of microns before reaching the interaction region in our simulations. Then, from Eq. ( 16), we find that the frequency of the probe pulse varies in time as

ω = 2C(t -t 0 )/w 2 t , (17)
and using the vacuum dispersion relation ω = c k to associate with wavelength λ c a time of interaction t c , we find

t c = t 0 + πcw 2 t C 1 λ c - 1 λ pr . ( 18
)
Measuring the position of a certain intensity peak in the different images of Fig. 14 and using Eq. ( 18) to estimate the corresponding time of interaction of the specific wavelength, cf. Table I, allows us to estimate the phase velocity of the wake to be v p = 2.8 × 10 8 m/s, which is pretty close to the value we find by tracking the corresponding density maximum directly in the PIC simulations, v p = 2.9 × 10 8 m/s. 14 for different central wavelength λ1 of the applied filter and corresponding relative time of interaction ∆t = t1 -t0 computed using Eq. ( 18).

The length of the probe pulse in our computations is limited by the transverse extent of the computational domain to about 30fs. During this time the length of the bubble does not change significantly. However, longer pulses could capture the lengthening and change of shape of the bubble in a single shot. We estimate that a pulse with transform limited duration τ pr < 20 fs chirped to 300 -600 fs would be sufficient to capture the lengthening of the bubble induced by pump pulse self-focusing and self-compression, according to the data presented in Ref. [7]. We note that a related technique which uses a train of short pulses of different wavelengths has recently been published [25].

VI. DISCUSSION AND CONCLUSIONS

We presented a computational study of few cycle shadowgraphy of LWFA. Our study relies on a combination of 3D PIC simulations of both pump and probe propagation in the plasma and subsequent post-processing that takes into account the effect of a typical imaging system. Both steps are crucial in order to analyze this technique. We have compared the simulated shadowgrams with experimental ones and have shown that the former reproduce all features found in the experiments. The ability to track probe propagation in the PIC simulations allowed us to identify the signature of the front of the bubble in the shadowgrams. This in turn facilitated the interpretation of the experimental shadowgrams and the measurement of important quantities such as the length of the bubble [7].

In order to facilitate the choice of probe pulse and imaging system parameters for future experiments, we studied the effect of probe pulse wavelength and chirp, numerical aperture of the imaging optics and object plane position. Moreover, we have shown that time-dependent information is retained in a pulse with short transform limited duration (adequate to resolve the plasma wave) that is linearly chirped to a longer duration. We envisage that if such a pulse is split after the interaction and filtered at different wavelengths, snapshots revealing the motion and evolution of the wake can be obtained.

In summary, we have proposed a method that allows us to construct synthetic ultrafast shadowgrams of plasma wakes by relating a near-field snapshot of the probe pulse to the integrated intensity at the image plane of a typical imaging system. The technique presented here could be used to analyse similar experimental methods, e.g., sequentially timed femtophotography [25], or have applications in other fields which involve the simulation of interaction of short probe pulses with matter, e.g., nanophotonics [26].

FIG. 1 .

 1 FIG.1. PIC simulation setup: The pump propagates along the x-direction and is polarized along y, while the probe propagates along the (negative) y-direction and is polarized along x.

FIG. 2 .

 2 FIG. 2. Probe pulse propagation across the wake. Panels (a)-(c): Three snapshots from the PIC simulation showing the probe electric field envelope squared |E probe x | 2 and contours of the plasma density corresponding to n = 2 n0. (d) Image obtained by recording the time-integrated Poynting flux Sy passing trough the plane yS = -20 µm (e) Time-integrated Poynting flux Sy after adjusting for focusing optics, assuming the object plane is at yo = 0.

FIG. 3 .

 3 FIG.3. Schematic setup of a so-called 4f -system. In paraxial approximation and by using thin-element approximation for the two lenses, we find the spatial Fourier transform of the object plane in the Fourier plane, and a (mirrored) image in the image plane. The aperture of the system can be represented by a (thin) mask in the Fourier plane.

18 |FIG. 5 .

 185 FIG. 5. Effect of taking into account the aperture in spatialFourier space, visualized in the kx, ky resp. kx, kz plane. The value of the numerical aperture used here is N A = 0.26. Because we are using Fourier transforms of the real fields, both cones (positive and negative ky) are present (cf. Fig.4).

FIG. 6 .

 6 FIG. 6. Snapshots of density distributions n/n0 from the PIC simulation for (a) z = 0 and (b) y = 0 at the time when the probe passes through y = 0. We introduced the coordinate ξ = x -ct. (c) Map of the y = 0 cross-section of n/(γ n0), (d) artificial shadowgram obtained with probe pulse with duration τpr = 12 fs, central wavelength λpr = 750 nm and object plane at y f = 0 (see Sec. II B). For the remainder of this paper we use the same pump-probe delay as in this figure and the same probe pulse characteristics unless otherwise specified. The color scale in all the following figures is the same as in panel (d).

FIG. 7 .

 7 FIG. 7. Measuring the plasma wave length at different propagation distances a) vgt = 527 µm and b) vgt = 1214 µm. Top: electron density, center: computed shadowgram, bottom: experimental shadowgram.

FIG. 8 .FIG. 9 .

 89 FIG. 8. (a) Experimental shadowgrams (b) computed shadowgrams, (c) density cross-sections from the simulation.

FIG. 10 .FIG. 11 .

 1011 FIG. 10. Effect of varying the object plane position on the shadowgram (a) yo = -10µm, (b) yo = 10µm. All images have been obtained from the same pump-probe simulation used in Fig. 6. The object plane position has been varied in post-processing. To be compared to Fig. 6(d) where yo = 0.

FIG. 12 .

 12 FIG. 12. Same as Fig. 6(d) but with N A = 0.35.

FIG. 13 .

 13 FIG. 13. Effect of varying the probe pulse duration through the introduction of a linear chirp. (a) Fourier limited probe pulse duration τpr = 4.4 fs, (b) τpr = 30 fs. To be compared to Fig. 6(d) where τpr = 12 fs. All images have been obtained for the same delay and the object plane position is fixed to yo = 0.

  FIG. 14.Effect of introducing a bandpass filter after the τpr = 30 fs chirped pulse of Fig.13(b) has propagated through the wake. The filters allow transmission of wavelengths (a) 0.6 ± 0.01 µm, (b) 0.7 ± 0.01 µm, (c) 0.8 ± 0.01 µm, (d) 0.9 ± 0.01 µm. The red arrow indicates the position of the intensity maximum that we track in order to deduce the phase velocity of the wake.

TABLE I .

 I Measurement of position x of intensity maximum indicated by a red arrow in Fig.

	λ1 (µm) x (µm) ∆t (fs)
	0.6	-25.0	-30.2
	0.7	-19.6	-8.6
	0.8	-14.7	7.5
	0.9	-11.1	20.0

ACKNOWLEDGMENTS

For the computations, the supercomputer at the Max Planck Computing and Data Facility at Garching was used. EPOCH was developed under UK EPSRC grants EP/G054940/1, EP/G055165/1 and EP/G056803/1. This study was supported by DFG (Grants No. TR18 B9, and No. KA 2869/2-1), BMBF (Contracts No. 05K10SJ2 and No. 03ZIK052), and European Regional Development Fund (EFRE). The collaboration was funded by LASERLAB-EUROPE (Grant No. 284464, ECs Seventh Framework Programme).

Appendix A: Definition of Fourier transforms

According to usual notations used in the context of optics, we define Fourier transforms with respect to space and time as follows:

Appendix B: How to compute complex electromagnetic fields from a temporal snapshot of real fields

In the general case, without any prior knowledge on the propagation direction, one can compute the complex fields ¯ E( k, t S ), ¯ B( k, t S ) from the real fields as follows: In spatial Fourier space we have according to Eqs. ( 4)-( 5)

we can rewrite Eq. (B2) as