Electro-optic shocks from blowout laser wakefields

Laser wakefield accelerators operating in the blowout regime exhibit unique forward scattering characteristics. When a cavitation region is formed, the drive pulse excites currents in the electron sheath that lead to conical emission of harmonic radiation in the near-forward direction. Theoretical expressions for the spectral-angular distribution in the far field are derived. These expressions allow for arbitrary temporal profiles of both the laser pulse and the sheath density. The characteristics of the emission can be related to key characteristics of the accelerating structure. The theoretical results are compared with the results of the experiment.


Introduction
Early laser wakefield accelerators were driven in gas jets by laser pulses approximately half a picosecond in duration. For typical gas jet densities (∼10 19 cm −3 ), many plasma periods fit within the time envelope of such a pulse. In this regime, forward scattering is characterized by the generation of well-defined Stokes and anti-Stokes sidebands. The generation of such sidebands is evidence of the generation of plasma waves, which serve as the accelerating structure. Furthermore, the broadening of the sidebands is evidence of wavebreaking and the onset of electron trapping [1].
In modern laser wakefield accelerators, the driving laser pulse is typically about 50 fs while the plasma density remains nearly 10 19 cm −3 . There are typically about two plasma periods in the time envelope of the laser pulse. This leads to the generation of inherently broad sidebands, or even pulse broadening that resembles self-phase modulation more than a Ramantype process [2,3]. In this case, the broadening of the forward scattered sidebands cannot be easily related to electron trapping since the sidebands are broad from the outset.
In the case where the intensity of the drive pulse is high enough so that a ≈ 2 or higher, where a = e A/mc −2 is the normalized vector potential, a cavitation region is formed due to ponderomotive blowout. This is accompanied by the formation of an electron sheath that surrounds the ion-rich 'bubble'. Currents driven in the sheath by the drive pulse lead to the generation of electro-optic shocks at the harmonics of the laser frequency. The electro-optic shock radiation is emitted conically in the near-forward direction. Because the emission process is connected with the electron sheath, the characteristics of the radiation are related to the blowout wakefield and its evolution. These blowout wakefields, or 'plasma bubbles', have been associated with the production of quasi-monoenergetic electrons in both the forward direction [4][5][6] and at large angles [7]. They are therefore of great importance for laser wakefield accelerators.
In a previous publication [8], we presented particle-in-cell simulations of electro-optic shock generation, along with a simplified theory. In this paper, the theory is extended to allow for arbitrary laser pulse shapes and is compared with experimental results. The focus is on the far-field angular-spectral distribution of the electro-optic shock, which is a useful experimental observable. Other aspects of the experimental results are discussed elsewhere [9].

Electro-optic shocks
The concept of an electro-optic shock was introduced by Auston in the context of subpicosecond lasers propagating in nonlinear crystals [10,11]. In this context, a slowly varying nonlinear polarization is driven by optical rectification. The frequency of the emitted radiation is of the order of the inverse of the laser pulse length, which is typically in the THz range. If the THz frequencies propagate more slowly than the laser pulse that causes their emission, the radiation pattern resembles that of Cherenkov radiation (i.e. an electromagnetic shock). The radiation is emitted conically with half-angle Figure 1. Wavenumber matching diagram for electro-optic shock generation in a plasma. If the plasmon wavevector, k p , is perpendicular to the pump photon wavevector, k(ω 0 ), then the electro-optic shock wavevector, k(2ω 0 ), must lie on the Cherenkov cone. Here, it is assumed that plasmon frequency is negligible.
optical rectification also leads to second harmonic generation. However, second harmonic radiation is not considered in traditional electro-optic shock theory. There are two possible reasons for this. Firstly, depending on the material, the dispersion could be such that (ω 0 ) > (2ω 0 ), in which case the Cherenkov angle would not be real. Secondly, since the scale of the drive pulse is much larger than the second harmonic wavelength, shocks emitted from different regions can destructively interfere. A strongly driven laser wakefield has (ω 0 ) < (2ω 0 ) and also produces structures that are smaller than the second harmonic wavelength. In particular, particle-in-cell simulations predict that the electron sheath that envelops the cavitation region has a thickness that is less than an optical wavelength. Hence, if second harmonic currents are driven in the sheath, an electro-optic shock is emitted.
Harmonic generation in a plasma is generally thought to be highly inefficient due to the impossibility of phase matching [12]. For example, a frequency-doubled photon has more momentum than two fundamental photons. Furthermore, attention is often restricted to odd harmonics because of the isotropy of the plasma [13,14]. In the presence of a strongly driven wake, however, the situation is different. Such a wake has a broad wavenumber spectrum, which can be regarded as a reservoir of plasmons with a broad range of momenta. The appropriate plasmon momentum can be added to the momentum of two pump photons to give the momentum of a frequency-doubled photon, as illustrated in figure 1. The choice of k p ⊥ k(ω 0 ) is motivated by the fact that the sheath is thinnest in the radial direction. It should also be noted that a nonlinear wake consists of both real and virtual plasmons, so that plasmons with ω = ω p may be present in the plasmon reservoir. In addition to providing wavenumber matching, the nonlinear wake provides an anisotropy in the form of density gradients. The laser pulse itself can also be considered anisotropic because of the ponderomotive force. Second harmonic generation in a plasma is discussed in [15][16][17][18].
The interpretation of figure 1 is as follows. Firstly, the plasmon wavenumber k p does not represent a simple plasma wave propagating perpendicular to the pump photons. Rather, it is a Fourier component selected from a broad spectrum. Secondly, although the interaction conserves momentum, it is not a traditional phase-matched interaction. Phase matching implies that a wave and a source function stay in phase for a finite time, resulting in nonlinear growth of the wave energy. In the situation considered here, the source region is infinitesimally thin, and the emitted wave immediately leaves this region due to the cone angle θ c . This results in secular growth of the wave energy. The situation is somewhat analogous to non-collinear phase matching with an infinitesimally narrow pump beam.

Perturbation theory for strongly driven wakes
Analyses of laser-plasma interactions often begin with a perturbation expansion in powers of the normalized vector potential, a. This approach is not directly applicable in the blowout regime due to the limited radius of convergence of the binomial expansion. In particular, the expansion of the relativistic Lorentz factor, γ ≈ (1 + a 2 ) 1/2 , only converges for a < 1. This difficulty is overcome by making use of two observations. Firstly, the general character of the solution is known beforehand from particle-in-cell simulations [19,20]. According to these simulations, the ponderomotive force expels electrons radially and forms a 'plasma bubble' that propagates quasi-statically with the laser pulse. Secondly, there may be interesting source currents in localized regions where a < 1. In the case of a plasma bubble, this localized region is the electron sheath. In the cavitation region, where the fields are strongest, there are very few electrons to contribute to the current. Instead, most electrons stream around the high-field region. If these sheath electrons follow orbits that take them outside the region where a > 1, a perturbation analysis becomes possible.
In order to perturbatively calculate the currents in the sheath, it is necessary to first specify the prescribed 'equilibrium' that is to be perturbed. The first part of this equilibrium is a quasistatically propagating plasma bubble. It is associated with an electron density n (0) and a fluid velocity v (0) . In this paper, it is assumed that the equilibrium velocity v (0) does not contribute to the radiation signal of interest. This is motivated by the expectation that the sheath density is highest where the radial component of v (0) is smallest. The plasma bubble is also associated with electromagnetic potentials A (0) and φ (0) , which are important for acceleration, but are not needed here.
The second part of the prescribed equilibrium is the laser field. This is represented by the normalized vector potential a (1) . The largest magnitude of a (1) that interacts with the sheath is the expansion parameter in the perturbation analysis. In deriving the perturbed current density, the Coulomb gauge, ∇ · a (1) = 0, is used. Note that the roles of the wake fields and the laser fields are somewhat reversed when compared with the usual perturbation theory. In the usual theory, the first-order laser fields perturb a zero-order uniform plasma to produce a second-order wake. Here, a high-order process produces a plasma bubble that is regarded as the new zeroorder equilibrium. This is then perturbed by the first-order laser field to produce a second-order scattered field.
The second-order perturbed current density can be written as j (2) = −n (0) ev (2) − n (1) ev (1) . To put this in a more useful form requires expressions for n (1) , v (1) and v (2) in terms of the sheath density n (0) and the laser potential a (1) . The momentum equation for a cold, unmagnetized plasma with immobile ions is This is written as a series by taking γ ≈ (1 + a 2 ) 1/2 and making a binomial expansion. To first order, To second order,

5
In both these equations, φ can be neglected provided the laser frequency is much higher than the plasma frequency, and provided it is the high-frequency component of v (2) that is of interest. The first-order velocity is then simply v (1) = ca (1) , consistent with the assumption γ ≈ (1 + a 2 ) 1/2 . The perturbed density n (1) is determined from the continuity equation where the first-order momentum equation and the gauge condition were used. While the Coulomb gauge is convenient for the calculation of j (2) , the scattered vector potential A (2) is most conveniently calculated in the Lorentz gauge. With this choice, the scattered wave is described by In the far field, the electric and magnetic fields are determined entirely by A (2) . In particular, (2) × n/c and E (2) = B (2) × n, where n is the unit vector pointing from the source to the observation point [21]. The above equations are sufficient to describe electro-optic shocks from blowout laser wakefields.

Second-order current density
In experiments with ultrashort pulse lasers, measuring the spectral properties of scattered radiation is often simpler than measuring temporal properties. Moreover, the measurement apparatus is often far from the interaction region. Hence, connection with experiment is conveniently made by deriving expressions for the angular-spectral distribution in the far field. The Fourier transform convention used here is The axial Fourier transformation F z is defined similarly. The circumflex is used to denote a Fourier-transformed variable in either t or z. The forward and reverse convolution theorems take the forms The scattered radiation is described in terms of ω and the spherical coordinates r = (x 2 + y 2 + z 2 ) 1/2 , θ = cos −1 z/r and ϕ = tan −1 y/x. However, the current density is described using the cylindrical coordinates ρ = (x 2 + y 2 ) 1/2 , ϕ and z. The Cartesian unit vectors are denoted by e x , e y and e z . The spherical unit vectors are denoted by e r , e θ and e ϕ . It should be noted that all expressions herein assume 0 θ π .
To express the nonlinear current density in terms of the sheath density and the incident wave, consider a thin cylindrical sheath propagating with velocity v in the z-direction: Here, n 0 is the sheath density, ρ is the sheath thickness, δ is the Dirac delta function, is the temporal profile of the sheath and h(z) measures the change in sheath density as the bubble propagates. Performing Fourier transformation giveŝ Let the incident wave have the form where a 0 is assumed real, but α may be complex. It should be noted that this form cannot account for group velocity dispersion. Performing Fourier transformation givesâ (1) = e x â + +â − , whereâ andâ In the following, attention is restricted to positive frequencies, which are contained inâ + . Negative frequencies can be accounted for using the symmetry properties of F. Cross correlations such asâ + * â − are ignored because they generate only low frequencies that are not of interest. With this in mind, Fourier analysis of the perturbed equations for n (1) and v (2) givesĵ The first term is due to electrons quivering in the direction of a density gradient. The second is due to ponderomotive effects. Carrying out the convolutions in equation (16) gives for the ponderomotive contribution where and δω = ω − 2ω 0 . The ∇n contribution iŝ where The total second harmonic current density for positive frequencies is thenĵ (2) + =ĵ +, p +ĵ +,∇n .

Second-order scattered field
In frequency space, the wave equation for the scattered vector potential becomes the Helmholtz equation where c 2 k 2 = ω 2 − ω 2 p . The solution to the Helmholtz equation iŝ In the far field, the positive frequency components arê Carrying out the integrations giveŝ where and κ = kρ 0 sinθ . Here, J 0 and J 1 are the Bessel functions of the first kind. The Cherenkov angle emerges from the factorĥ(K ), which appears in all terms. As the bubble propagation distance increases,ĥ(K ) approaches a delta function. It then holds that K = 0, which leads to equation (1) in the limit δω → 0. Calculation of the scattered electric field proceeds fromÊ (2) = −(iω/c)(Â (2) × e r ) × e r . Note that this expression implicitly accounts for the scalar potential (which satisfies a wave equation in the chosen gauge), so that any component ofÊ can effectively depend on any other component ofÂ. Expressing the components ofÊ in spherical coordinates giveŝ where E 0 = J 0 (κ) sin θ e θ k F n cos θ cos 2 ϕ − k 0 F z − e ϕ k F n sin 2ϕ 2 (31) and E 1 = 1 2 k ⊥ J 1 (κ) e θ cos θ 2F n cos 2 ϕ − F p − e ϕ F n sin 2ϕ .
Note that E 0 was dropped in [8] based on the assumption that θ 1. However, E 0 /E 1 is actually of order (k 0 /k ⊥ )tanθ , which is not small in general. The expression for the electric field can still be simplified by assuming θ 1 and δω ω 0 , but with the additional constraint k 0 /k ⊥ 1. When this is done, and The scattered spectral intensity,Î + = (c/2)Ê (2) + ·Ê (2) * + , is then given bŷ Finally, the total energy radiated is given by Useful expressions for U can only be obtained by making additional approximations.

Time domain solution
In considering the solution for the scattered field in the time domain, the fact that the diagnostic apparatus is far from the plasma has to be considered. The error introduced by the plasma-vacuum transition takes two forms. Firstly, there may be refraction as the scattered rays leave the plasma, altering the predicted angular distribution. This affects even the frequency domain solution, but is typically a small error. Secondly, group velocity dispersion is arrested as the scattered rays leave the plasma, altering the predicted temporal distribution.
To account for the lack of group velocity dispersion in the vacuum region it is assumed that the plasma region is small enough for the group velocity dispersion to be wholly unimportant. For typical parameters, this limits the extent of the plasma to the centimeter scale. In order to deliberately neglect group velocity dispersion, the overall phase factor e ikr is simply replaced by e iωr/c .
To carry out the inverse Fourier transformation ofÊ (2) + , the approximation ω ≈ 2ω 0 has to be used wherever the frequency dependence is weak. This results in where andŜ (r, θ, ω) =ĥ(K )F p (ω)e iωr/c (41) 9 is the only factor that still depends on ω. This is simplified further by taking the laser pulse to be long compared with the bubble length so thatα(ω − ω 0 ) → δ(ω − ω 0 ), in which case F p (ω) =ĝ(δω)/2ω 0 . Then, the inverse transform ofŜ is and the scattered field in the time domain is The expression for S(r, θ, t) can be put in a more useful form as follows. The factor involvingĝ can be written as In order to expand the factor involvingĥ, the wavenumber is expanded about 2ω 0 , giving where η 2 = (4ω 2 0 − ω 2 p ) 1/2 /2ω 0 is the refractive index at the second harmonic frequency. The refractive index at the fundamental frequency is η 1 = (ω 2 0 − ω 2 p ) 1/2 /ω 0 . Taking the bubble propagation speed, v, to be the group velocity, cη 1 , and using cos θ c = η 1 /η 2 , we obtain This leads to The small angle approximation is made by first expanding in δθ = θ − θ c , then expanding in θ c and then expanding in δθ/θ c . This gives which leads to where τ = t − r/c. An important limiting case is when the bubble propagates for a long enough distance so that This corresponds to L /θ 2 c , where L is the characteristic propagation distance and is the characteristic length of the bubble. Then, the factor involving h may be pulled out of the convolution integral. Furthermore, it is safe at this point to take η 1 ≈ 1, as is the case in a tenuous plasma. Then, In this limit, frequency is perfectly correlated with angle, the angular distribution relative to θ c is the bubble spectrumĝ, and the length of the scattered pulse is Lθ 2 c /c. The total energy emitted is proportional to L Q 2 , where Q is the charge in the sheath.
The finite length of the scattered pulse is a manifestation of the non-ideal nature of the shock. In a plasma, the group velocity and the phase velocity are always unequal, even to lowest order. In particular, when the refractive index η is less than unity, as is the case in a plasma, then the lowest order phase velocity is c/η and the lowest order group velocity is cη. To see the effect of this, consider the case of infinite propagation distance so that the frequency-dependent Cherenkov angle is found by setting K = 0. Defining the phase velocity v(ω) = ω/k and taking the bubble velocity v to be the group velocity at the fundamental, v g (ω 0 ), we obtaiñ If ω = 2ω 0 , the expression in brackets reduces to the ratio of phase velocities as in equation (1).
On the other hand, if it were possible to take v(ω 0 ) = v g (ω 0 ), then one would have the more familiar case where the expression in brackets is the ratio of the emitted phase velocity to the source group velocity. In that case, it would be possible to obtain perfect shocks by keeping only the lowest order of dispersion (here, a perfect shock is one where the wave energy is localized on the surface of a cone at every instant). In the present case, there seems to be no logical way to consider perfect shocks, even in principle. The shock is fundamentally non-ideal.

Angular-spectral distributions
In order to plot the angular-spectral distribution given in equation (35), the temporal profiles of the wake and laser pulse must be given. The basis for choosing these profiles is twofold. Firstly, the density of the ambient plasma and the vacuum laser pulse parameters are chosen to be consistent with experiments in which high-energy electrons are observed. Secondly, these are used as initial conditions in particle-in-cell simulations, and the simulated characteristics of the laser pulse and wake, after propagating some distance in the plasma, are used to construct a simple model. The assumed form of the laser envelope is where τ L is the pulse width, which may be shorter than the vacuum pulse width. The assumed form of the temporal profile of the sheath is where H is the Heaviside step function, and τ d is the wake decay time. The assumed form of the propagation history (i.e. the change in sheath density with propagation distance) is  The dashed line is α(t) (laser profile) and the solid line is g(t) (wake profile). The function g(t) measures the temporal variation in the sheath density. The laser pulse is presumed to be shorter than the vacuum pulse length due to effects such as self-focusing and self-modulation.
where L is the propagation distance. The particular parameters used here are given in table 1. Four cases are considered: two different ambient densities and two different propagation lengths. The sheath density is held in fixed proportion to the ambient density. The functions α and g are plotted in figure 2.
Using the parameters of table 1, four angular-spectral distributions are plotted in figure 3. The distributions are evaluated at ϕ = 0, i.e. in the plane of polarization. The units are converted into energy per unit frequency per unit solid angle so that the distance to the measurement apparatus cancels out. The low-density case, n e = 10 19 cm −3 , is plotted in figures 3(a) and (b) for L = 100 and 500 µm, respectively. The high-density case, n e = 4 × 10 19 cm −3 , is plotted in figures 3(c) and (d), for L = 100 and 500 µm, respectively. In the low-density case, there is only one strong bubble within the laser pulse envelope, as shown in figure 2(a). This leads to a single-peaked second harmonic spectrum. As the propagation distance is increased, the angular spread at a given frequency narrows, consistent with the prediction that angle and frequency are perfectly correlated as L → ∞. The central angle in either case is close to the Cherenkov angle, θ c = 0.066. In the higher density case, two bubbles are inside the laser pulse envelope, as shown in figure 2(b). This leads to a modulated second harmonic spectrum. The correlation between angle and frequency is again strengthened as the propagation distance is increased. The central angle is again close to the Cherenkov angle, θ c = 0.13. The total energy emitted was found by numerically integrating equation (37) in the four cases corresponding to figures 3(a)-(d). In the low-density case, the total energy was U = 14 µJ for L = 100 µm and U = 92 µJ for L = 500 µm. In the high-density case, the total energy was U = 190 µJ for L = 100 µm and U = 1250 µJ for L = 500 µm. These results depart from the simple scalings implied by equation (54) by up to 30%. Some of this error is due to taking k → 2k 0 in certain factors when forming the time domain solution and some is due to finite propagation length.
An experimental angular-spectral distribution is shown in figure 4. This is obtained by firing a 10 TW, 50 fs, 0.8 µm Ti : sapphire laser pulse, focused with an f /8 off-axis parabolic mirror, into a 1 mm long gas jet. An f = 15 cm spherical mirror is used to collect the light reflected from a pellicle and focus it back through the pellice onto a diffuser placed a distance f from the spherical mirror. This maps the polar angle θ to the radial coordinate on the diffuser. The diffuser is then imaged onto the slit of an imaging spectrometer to produce an angular-spectral distribution. For the case considered here, interferometric measurements give the ambient electron density as n e = 3 × 10 19 cm −3 . Figure 4 clearly indicates conical emission for ω/ω 0 > 1.7. The measured angle is close to the Cherenkov angle, θ c = 0.11. The image also includes a feature corresponding to direct forward scattering near ω/ω 0 = 1.6. This is the extremity of an anti-Stokes-like sideband. Modulations in the scattered spectrum are also evident, suggesting the presence of multiple bubbles within the laser envelope, as in figures 3(c) and (d). The angular spread in the conical emission suggests that the bubble propagation distance is much shorter than the 1 mm length of the gas jet, although the resolution of the optical system must be taken into account. Finally, the correlation between angle and frequency is similar to the predicted correlation, i.e. higher frequencies are scattered into smaller angles.
The experimental result differs from the theoretical in the overall shift of the scattered light. While the theoretical distribution has side lobes on either side of ω/ω 0 = 2, the experimental one has side lobes only on the red side. This might be due to nonlinear frequency shifts in the pump radiation itself, or it might be related to the frequency response of the measurement apparatus. The manufacturer of the spectrometer quotes the high-frequency cutoff as λ = 0.365 µm, or ω/ω 0 ≈ 2.2.

Conclusion
Electro-optic shocks emitted from blowout laser wakefields provide new diagnostic opportunities for laser wakefield accelerator research. Theoretical expressions for the angularspectral distribution reproduce key features of experimentally observed electro-optic shocks. These features include conical emission at the Cherenkov angle, inverse correlation between angle and frequency and the development of a modulated second harmonic spectrum at high densities. Because the emission process is closely tied to the electron sheath associated with a plasma bubble, the electro-optic shock characteristics can be related to the characteristics of strongly driven wakes. For example, if the electro-optic shock spectrum is modulated, multiple bubbles reside within the laser pulse. If the angular and spectral distributions are highly correlated, the bubble propagates quasi-statically over a long distance. These characteristics have important implications for electron trapping and acceleration and are therefore of great interest in laser wakefield accelerator research.