Theory of circularly polarized harmonic generation using bi-colour lasers in underdense plasmas

Circularly polarized (CP) extreme ultraviolet- and x-ray radiation is an essential tool for analyzing the magnetic properties of materials. Elliptically polarized high harmonic generation (HHG) has been demonstrated by focusing bi-chromatic (800 + 400 nm wavelengths), counter-rotating CP laser pulses into gas targets (Fleischer et al 2014 Nat. Photonics 8 543). More recent theoretical studies indicate that a bi-circular laser driver can also work in both under- and overdense plasmas with analogous selection rules to those in gases: for example, every third harmonic is suppressed and adjacent harmonics have opposite helicity for counter-polarized CP ω 0 and 2ω 0 pumps. In this work, an analytical theory of bi-circular HHG from underdense plasmas is formulated which provides quantitative predictions of harmonic efficiency scaling, selectivity and helicity for both co- and counter-polarized drivers of arbitrary frequency ratio. This is compared to a fully non-linear, one-dimensional fluid model and particle-in-cell simulations, showing good agreement with both.


Introduction
Conversion of laser light to multiples of its fundamental frequency, or high harmonic generation (HHG), has been Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. studied and applied as a means of producing coherent extreme-ultraviolet (EUV) and attosecond pulses for some time now. HHG can occur both in gas and plasma targets through non-linear interaction between a laser field and bound (gas) or free (plasma) electrons. In a gas target this process occurs via recombination of ionized electrons with the parent ions, which is commonly described by a three-step model [1,2,3]. In a fully ionized plasma target, HHG arises through acceleration and deceleration of free electrons in the field of the laser. Interest in producing short-wavelength light by this means began soon after invention of the first lasers in 1960 [4], when a number of theoretical studies appeared devoted to laser interaction with matter and feasibility of harmonic generation by non-linear Thomson scattering [5].
The first experiments using gas targets date back to the late 1980s [6,7], when picosecond laser pulses were focused into a jet of noble gas. To date, this method has been refined considerably to become one of the choice techniques for producing EUV light with a tabletop laser [8]. There are many ways of improving HHG efficiency, such as phase-matching [9], polychromatic drive [10], or extending the cutoff frequency [11], but the single-shot photon yield-or harmonic intensity at a particular wavelength-is ultimately limited by ionization of the gas at high laser intensity.
HHG can also occur in a fully ionized gas target (underdense plasma) in the relativistic regime, a fact that prompted a series of early analytical studies [12][13][14][15] using onedimensional (1D) perturbational models which predicted rather low efficiencies. For this reason most theoretical and experimental attention has since been devoted to studying plasma harmonics from overdense foil targets [16], where different mechanisms prevail depending on the laser pulse contrast ratio [17][18][19]. Underdense plasmas have mainly been examined for low-order harmonics [20] but have the advantage of being experimentally easier to handle with the same optical setup as for gas-generated HHG, thus benefiting from continuous target medium replenishment more suited to kHz laser operation than solid targets.
Normally, harmonics are completely suppressed when using a CP driver pulse both in gas and plasma media. In the former case, bound electrons ionized by a CP electric field will not return to the atom and so recombination events with associated high-energy photon emission cannot take place. In an ionized medium, electrons remain on circular orbits with constant speed, suppressing the high-frequency density fluctuations and relativistic nonlinearities driven by linearly polarized (LP) light which act as the primary harmonic source.
Generation of efficient CP harmonics in gas was first demonstrated by Fleischer et al in a breakthrough experiment [31] using two bichromatic drivers. This work has since prompted a flurry of activity, including the demonstration of bright, phase-matched CP high harmonics [32], non-collinear [33] generation of CP high harmonics, harmonic generation driven by time-delayed, few-cycle ω-2ω pulses [34] and with controlled helicity for ω-2ω and ω-3ω fields [35]. The selectivity of gas harmonics has been recently studied by Venzke et al [36]. Obliquely incident, counter-rotating CP laser pulses have been used by Hernandez-Garcia et al [37] as a means of generating isolated attosecond pulses. Bicircular harmonic generation can be used to probe dynamical symmetry [38,39].
Polarization control of harmonic generation on a plasma surface with the aim of obtaining a single attosecond x-ray pulse was first proposed by Baeva et al [40]. Subsequently, CP harmonic generation from a relativistic plasma mirror was studied by Chen et al [41] and Ma et al [42] using an obliquely incident CP pulse. The idea of mixing bi-chromatic pulses as proposed in gas harmonic generation has also been explored numerically for overdense plasmas at ultrarelativistic laser intensities [43] and analytically for underdense plasmas and more moderate laser amplitudes [44].
In this paper we examine the harmonic generation mechanism from free electrons in a plasma target, and in particular, explore options for achieving CP EUV light with high efficiency and selectivity of the harmonic number. Due to the small wavelength of higher harmonics in the nanometer (nm) range, compared to optical laser wavelengths (400 + 800 nm), numerical study of this phenomenon normally demands very high resolution particle-in-cell (PIC) simulation and therefore, computation time. A new, complementary non-linear plasma fluid approach is introduced here which offers a faster, noisefree, hi-fidelity alternative to PIC while retaining the essential harmonic generation physics. This model is able to guide both experiments and more complete multi-dimensional PIC simulations.

Analytical fluid model: harmonic selection rules, helicity and efficiency
Starting from the standard EM wave equation for the vector potential, we observe that harmonics are generated by the source term in the right hand-side of the equation, therefore, to determine the harmonic content of the radiation, the fully non-linear current J = en e p/γm must be evaluated, which in turn amounts to solving for the longitudinal density perturbation and the relativistic factor. In order to fully determine the density we need three additional equations: the Lorentz force equation together with Poisson's equation and the continuity equation, which in one dimension for arbitrary amplitudes n e , ϕ and p x can be written as: where p x represents the dimensionless longitudinal momentum (normalized to mc), u = p x /γ the corresponding fluid velocity (c), ϕ the scalar potential (mc 2 /e), A ⊥ the laser wave vector amplitude (mc/e) and n e is the electron density normalized to the background ion density (n 0 ). Henceforth we will assume that the laser driving force ∂A 2 ⊥ /∂x in equation (2) consists of two bi-colour CP pulses travelling in the x-direction with arbitrary amplitudes which can be represented by the transverse wave vectors: where θ = ω 0 t − k 0 x, θ ′ = qω 0 t − k q x and q is an integer which is the frequency ratio of the two co-(+) or counterpolarized (−) laser pulses. The present work will mainly consider drivers with a ratio q = 2 corresponding to the commonly available laboratory laser configuration ω 0 and 2ω 0 . Keeping this ratio arbitrary for now, we first observe that: where Θ ± = θ ± θ ′ . For counter-polarized pulses, the driving term is cosΘ + ; in the case of co-polarized pulses, this is cosΘ − . This implies that the fluid momentum p x , density n e and scalar potential all oscillate with harmonics of the driver, so that we can expand each fluid quantity accordingly: Likewise, the relativistic parameter γ −1 expanded up to O(ϵ 2 ) order is: where ε is an ordering parameter such that the driving term is considered to be O(ϵ) and γ 0 = 1 + A 2 0 /2 + A ′2 0 /2. Substituting these expansions into equation (2) leads to a set of algeb-raic equations which can be solved by iteration. In the first iteration, the O(ϵ) equations yield: Poisson's equation: Continuity equation: where K ± q = (k 0 ± k q ), Ω ± q = (q ± 1)ω 0 and k p = ω p /c the plasma wave number. Substituting p 1 and ϕ 1 into (7) gives an explicit expression for the O(ε) plasma density amplitude: for counter-(+) and co-polarized (−) pulses. The corresponding first order currents are then: for counter polarized pulses and, for co-polarized pulses, where α ± = c 2 K ± q 2 /(Ω ± q 2 − ω ′2 p ) and ω ′2 p = ω 2 p /γ 0 . Note that in the low density limit (ω p ≪ ω 0 ), the coefficient α ± − 1 term reduces to (ω ′ p /Ω ± q ) 2 , meaning that the first order non-linear current will scale as O(ω ′ p ) 4 , or the square of the electron density, just as in the original analyses of third harmonic generation via a single-frequency, LP laser [13,14,45].
In the second iteration, equation (2) is solved to O(ϵ 2 ) for n e , ϕ and p x by substituting the expansions from equation (5) and then gathering all second order terms. This gives a new set of algebraic equations analogous to equations (7)- (9): We see immediately that these expressions also contain n 1 and p 1 , for which we can substitute the solutions already found from equations (9) and (10). In this way, equations (13)- (15) can be solved for the second order density amplitude: whereα ± = (Ω ± q /cK ± q )α. Inspection of these coefficients reveals that there appear to be contributions O(ω ′ p ) 2 , but it turns out that these end up cancelling other terms of the same order in the non-linear current, finally leaving a leading contribution scaling as (ω ′ p ) 6 to this order. Putting together the results from equations (6), (10) and (16) we can now calculate complete expressions for J y,z = ω 2 p (1 + n)γ −1 A y,z up to O(ϵ 2 ) in the low density limit, or ω p ≪ ω 0 : Inspection of the phases of the non-linear current components (J y,z ≡ nA y,z ) in equations (17)- (20) immediately reveals that the combination of bi-colour pumps, whether counter-or copolarized, always results in CP harmonics. Moreover, the helicity of the harmonics is determined automatically through the relative sign of each current component: for counter-polarized pulses, the terms θ, θ ′ + 2θ, 2θ ′ + 3θ, have the same helicity (LCP) and are counter-polarized relative to θ ′ , 2θ ′ + θ, 3θ ′ + 2θ (RCP). In other words, for the q = 2 case (ω 0 + 2ω 0 pumps), modes ω 0 , 4ω 0 , 7ω 0 are LCP, whereas 2ω 0 , 5ω 0 , 8ω 0 are RCP. For co-polarized drivers, 2θ ′ − θ, 3θ ′ − 2θ, has the same helicity (LCP) as the main pumps, θ and θ ′ , which are counter-polarized relative to θ ′ − 2θ, 2θ ′ − 3θ, (RCP). For q = 2 this results in same helicity (LCP) for ω 0 , 2ω 0 , 3ω 0 , 4ω 0 . Selection rules for this and higher values of the pump frequency ratio q are summarized in table 1. It is worth taking a moment to consider the origin of the selection rules which arise using bi-circular pulses. To do this, we extend the expansion of the γ −1 term in the current (equation (6)) to arbitrary order: Note that cos m Θ ± term in γ −1 expands to cosmΘ harmonics. Multiplying this term by the density oscillating at the  (17)- (20).
harmonics. If we now expand A y = a 1 cos θ + a 2 cos 2θ + · · · + a m cos mθ, the product of (1 + n)γ −1 with A y automatically excludes (q ± 1)m harmonic terms. This generalizes the selection rules: only harmonics with mode number (q + 1)m ± 1 can be generated in case of counter-rotating pulses or (q − 1)m ± 1 when the pumps are co-polarized. The (−1) m factor contributes the sign of each term and therefore the helicity of harmonics. Expanding γ −1 one can show that for counter-polarized pulse all (q + 1)m + 1 are left circularly polarized (LCP) and (q + 1)m − 1 harmonics are right circularly polarized (RCP). In case of co-polarized pulses the (q − 1)m + 1 harmonics have the same helicity as the main pumps (LCP), opposite to the helicity of the (q − 1)m − 1 harmonics. We also note that these findings are consistent with the work of Sharma et al [44], who examined the non-linear current up to O(ϵ) in our notation, or harmonics (q + 1)m ± 1 ⩽ 5.
We now proceed to calculate the harmonic amplitudes. Substituting the A y = a 1 cos θ + a 2 cos 2θ + · · · + a m cos mθ in equation (1) gives the dispersion relation: where m ⩾ 2. The harmonic non-linear (NL) current for counter-and co-polarized pulses according to equations (17)-(20) will be: with cos terms indicating the y-and sin terms the z-component of the non-linear current density. To obtain the harmonic efficiency we equate equation (23) with the matching resonant cosmΘ term on the left-hand side of the expanded wave equation (1). This will lead to corresponding steady-state harmonic amplitudes: which by extrapolating the series, suggests a general rule for the leading contribution (see also the appendix) to the amplitude of each harmonic to arbitrary order of: for m ∈ [0, 1, 2, …] when pulses are counter-polarized. Whereas for co-polarized pulses the leading contribution (see also the appendix) in the amplitude of each arbitrary harmonic number will be: for m ∈ [0, 1, 2, …] and a −1 ≡ a 1 . From now on we limit our calculations to counter-polarized case. The power scaling is obtained using the electric field strength of the harmonic P n ∝ I n ∝ E 2 n ∝ ω 2 n a 2 n :

Nonlinear fluid model
In order to confirm that the findings deduced from the preceding analytical study are equally applicable to higher harmonic orders, and to gain deeper insight into the propagation dynamics of CP harmonics from underdense plasmas, we have developed a fully non-linear 1D3V plasma fluid model. This is constructed by combining a dispersion-free EM wave solver with the well-known 1D solution for the longitudinal plasma wakefield. In the following we describe the algorithm in detail before applying it to verify some of the results of the preceding theoretical analysis.

The directional splitting field solver
Unlike the FDTD scheme where the following normalized Maxwell's equations are solved explicitly in each grid point at each time step ∆t: in the directional splitting (DS) scheme, fields are integrated along their propagation direction on a single computational grid [46]. Equations (28) and (29) can then be rewritten as: where ± F y = (B z ± E y )/2 and ± F z = (B y ∓ E z )/2. Therefore, the electromagnetic field is updated based on the combination of westward ( + F) and eastward ( − F) propagating waves: and current densities are space-centred with: E-and B-fields are updated at each grid point and time-step using directional fields: Making use of the identity p y,z = eA y,z in this geometry, we observe that we can apply the DS integration scheme a second time to obtain the vector potential and thus the transverse electron fluid momentum: ( and finally current density will be updated using the fluid momentum: where γ = (1 + p 2 ⊥ ) 1/2 /(1 − β x ) 1/2 with β x the normalized longitudinal plasma velocity and p 2 ⊥ = p 2 y + p 2 z .

Plasma oscillation: quasi-static approximation
The field solver introduced in section 3.1 updates the EM-field via the transverse source terms J y and J z which contain contributions from the transverse momentum of the electrons in the laser contribution and the longitudinal plasma oscillation n e . As shown in section 2, in this context the density oscillation acts to nearly cancel the contribution from the relativistic γ factor, reducing the final harmonic amplitude, so it is essential to treat this component accurately. In 1D geometry, the density can be found by solving Poisson's equation at each iteration of the EM field. To do this we adopt the quasi-static approximation (QSA) of Esarey et al [45] to calculate the scalar potential in a coordinate system moving with the laser pulse with an arbitrary group velocity v g , so that Poisson's equation then becomes [45,47]: where ψ = β g = √ 1 − n e /n c the normalized group velocity, γ g its corresponding Lorenz factor; ϕ and a the normalized scalar and vector potential, respectively. Within the QSA, the normalized density and the velocity can be found from algebraic expressions involving the intermediate ψ parameter thus: Therefore, the current densities J y,z = en e p y,z /γ in each EM field-solver iteration of equation (30) will be found after updating not only the p y and p z , but also β x , which in turn contributes to the new γ: To maintain consistency with the earlier theoretical analysis, our fluid model solves the Poisson's equation in the low density limit where n e ≪ n c , β g → 1, so that equations (37) and (38) reduce to:

Fluid model results
Here the results of the fluid code introduced in the previous section are presented. The simulations are made for an underdense plasma slab of ∼49 µm with 10 µm distance from left boundary in a box of total length 120 µm. For the pump frequency ratio of q = 2, based on the selection rules deduced in section 2, every third harmonic for counter-polarized pumps should be suppressed, whereas all harmonics should be present for co-polarized pumps. Figure 1 shows the numerical result for an underdense plasma with n e = 0.1n c and a combined laser intensity of I total = I ω + I 2ω = 1.1 × 10 18 W cm −2 with duration of cτ = 9 µm FWHM where I 2ω0 /I ω0 = 2.0 and q = 2. As can be seen the charge density oscillates predominantly at 3ω 0 for counter and ω 0 for co-polarized case. This is the O(ϵ) density perturbation in equation (5). Signatures of higher order contributions are visible in the current density spectra, with every third harmonic suppressed for counter-polarized pumps, but generation of all harmonics for co-polarized pumps-figures 1(a), (b) and (d). Figures 1(c) and (e) show the total electric field at the time where the pulses reach the middle of the target for counter-and co-polarized pulses, respectively. Figure 2 shows examples of the corresponding harmonic spectra computed using the fluid code for pump frequency ratios q = 2, 3, 4, 5, with same total laser intensity and target density. As expected from section 2, for counter-polarized (solid red line) only (q + 1)m ± 1 and for co-polarized (blue dashed line) (q − 1)m ± 1 appear; modes not adhering to these rules are forbidden. Note that the efficiency (power scaling P n /P 0 as defined in section 2) in case of co-polarized pulses falls off more rapidly with harmonic number than in the counter-polarized case.
We have also compared the harmonic efficiencies from the fluid model and the analytical model equation (27) for counter-polarized pulses, figure 3. For a total laser intensity I = 1.1 × 10 18 W cm −2 and five target densities n e /n c = 0.0005, 0.005, 0.01, 0.03, 0.05, we confirm that increasing the density also increases the efficiency with the expected scaling. Furthermore, the decrease in harmonic efficiency from fourth and fifth harmonic to seventh and eighth is the same for both numerical and analytical model, figure 3(a). For both models, harmonic efficiency increases with increasing laser intensity from 1.1 × 10 17 to 1.1 × 10 19 W cm −2 for the same target density (n = 0.05n c ), although it saturates for strongly relativistic laser intensities, figure 3(b). Although the analytical model has been formulated for weakly relativistic laser intensities, the numerical fluid model is valid for any arbitrary laser intensity, which allows us to test the validity and applicability of the former in this regime.

Free electron trajectory
Gas target harmonic generation is highly dependant on recombination of electrons and the atom. In the field of a CP laser pulse, electrons do not have the chance to recombine with the atom. Using bi-colour counter-rotating pulses, electron recombination is only possible because of the special return trajectory [28,32]. By contrast, a free electron can generate harmonics in a bi-circular field without necessarily returning back to its origin. In this case, harmonic generation occurs by electrons being accelerated and decelerated during one cycle of their trajectory.
In order to illustrate this fundamental difference between harmonic generation in gas and plasma media, it is helpful to examine single electron trajectories. We start from Lorentz equation dp/dt = −e(E + v × B) and associated energy equation d(γmc 2 )/dt = −e(v · E), where γ = (1 + p 2 x + p 2 ⊥ ) 1/2 . Integrating the transverse and longitudinal components of the Lorentz force for an electromagnetic wave described by equation (3) leads to the well-known general solution for the normalized momenta [48] p ⊥ = A; p x = p 2 ⊥ /2, if we assume that the electron is at rest (p x = p ⊥ = 0, γ = 1) (c) q = 4, counter-polarized: 5m±1, co-polarized 3m±1; (d) q = 5, counter-polarized 6m±1, co-polarized 4m±1. The spectra are on a logarithmic scale. Note that harmonic fall-off, in particular for the co-polarized case, is quite rapid. This causes higher order harmonics (m or ε higher than 3 for this specific laser and target configuration) to fall below the background level. before the laser arrives. For the bi-circular pulses defined in equation (3), this results in the following explicit expressions for the momentum components: Integration of equations (41) gives the electron orbit in the field of two counter-and co-polarized pulses of different frequencies and arbitrary amplitude: . This is easy to show from the solution in equation (42), where for q = 2 the electron orbit has y = z = 0 at phases θ = 0, 2π/3 and 4π/3 only if A 0 = 2A ′ 0 . Of course, the electron also experiences a longitudinal drift just as it does for single-pulse irradiation.

Comparison to PIC simulation
Additional numerical confirmation of our findings comes from 1D and 2D PIC simulations performed using the EPOCH code [49]. All 1D (2D) simulations were made using a 120 µm (120 × 40 µm 2 ) box discretized by a computational grid with dimensions n x = 60 000 (n x × n y = 60 000 × 2000). An underdense preionized helium plasma of ∼49 µm length is placed at 10 µm from the left boundary with 100 particles per cell and immobile ions. The laser pulse has 0.8 µm wavelength and duration of τ = 30 fs. The target densities are n = 0.005n c , 0.05n c and the laser has normalized intensities a 0 = 0.2, 0.7. To get the harmonic spectra, the transmitted pulse is allowed to propagate beyond the right edge of the plasma, where it is then Fourier transformed (x to k) in order to reduce the transient noise within the plasma. This is the equivalent of taking the frequency spectrum of a time-series collected from a probe placed in the vacuum region, which is the procedure used for numerical fluid model.
The chosen frequency ratio is q = 2, where we expect generation of 3m±1 harmonics. According to the figure 5, 1D PIC simulation confirms the selection rules predicted in sections 2 and 3. Increasing the density increases the efficiency of the harmonics as expected from equation (27), see red and green spectra in figure 5. Higher laser intensities result in higher har-  monic efficiencies as well, green and blue spectra. However, for the higher laser intensities, increasing the target density increases the noise level and harmonics are not distinguished. This is because parametric instabilities such as Raman scattering and relativistic self-modulation will play an increasingly dominant role for high intensity/density combinations. In fact sidebands at ω 0 ± ω p and 2ω 0 ± ω p are clearly visible around the pumps for the n e /n c = 0.05 case, which is a classic signal of forward Raman scattering. Note that Raman scattering is in principle also included in the fluid model, but grows from a much lower level than in the PIC simulation. Despite these  [44] for target density of ne = 0.05nc and total laser intensity of I total = 1.1 × 10 18 W cm −2 .
quantitative differences, the selection rules are as predicted; with every third harmonic absent in this case.
Finally, the numerical results from fluid model and PIC code are compared to the analytical results from section 2 and Sharma et al [44] and the results are shown in figure 6. For fourth and fifth harmonics, the efficiency computed by PIC simulation and fluid model are in a fair agreement with our analytical model, but discrepancies start to appear at higher frequencies, which might be due to neglect or suppression of parametric instabilities in the fluid model compared to PIC. Similar differences between analytical and PIC results for higher-order LP-generated harmonics were previously noted by Mori [13]. We note that Sharma et al predicts higher values Table 2. Comparison of harmonic efficiency (Pn/P 0 ) for analytical and numerical models for target density of ne = 0.05nc and total laser intensity of I total = 1.1 × 10 18 W cm −2 . The values here correspond to the harmonic efficiencies in figure 6 and not the maximum possible values of each harmonic.  for 4ω 0 and 5ω 0 -table 2, which may be partly due to their omission of the γ 0 relativistic correction.

Dephasing length
According to the dispersion relation in section 2 the phase velocity of harmonics v ϕm = ω m /k m exhibit differences inside the plasma depending on the harmonic number. After a certain dephasing length, harmonics and their pump waves will be π out of phase which ultimately results in an oscillatory rather than constantly growing solution [13,44,45]. For the case with q = 2, as discussed in section 2, every third harmonic is suppressed. This behaviour is demonstrated in figure 7, where harmonics up to fifth from the analytical model, numerical fluid model, and the PIC code are plotted. Dephasing happens after a propagation distance π/∆k m , where for the fourth, fifth, seventh and eighth counterpolarized driven harmonics gives rise to the following relations: where k m is computed from equation (22). In fact the ∆k m is the difference between the harmonic wave number in the plasma and vacuum. Therefore, the higher order harmonics dephasing length can be calculated, as shown in table 3, which compares the dephasing lengths calculated from equation (43) with the values measured from the numerical fluid model and PIC code, respectively.
The fluid code has a number of advantages over the PIC code in this context: first, the fluid code utilizes a electromagnetic field-solver based on the DS scheme (see section 3.1), which has much less numerical dispersion than the FDTD scheme used here for the EPOCH code [46]. The biggest factor however is the avoidance of particle-grid mapping needed in PIC to obtain the currents, effectively eliminating the usual cause of noise. High noise levels are transmitted to the harmonics, making it difficult to determine the dephasing length accurately, unless much higher resolution or computationally expensive, higher-order schemes are used.
Not surprisingly, the short dephasing length is consistent with that found in the early work on HHG from underdense plasmas [13,14,45]. Various schemes have been proposed to overcome this limitation such as density ripples to allow pumps and harmonic to 'coast' back into phase [15,50,51], or non-collinear driver pulses with the help of higher-order modes in a focused beam [20].

Discussion and conclusion
One of the main outcomes of this study is the construction of a non-linear analytical model for the properties of harmonics driven by bi-circular pumps. The model predicts a general set of selection rules in which harmonic orders (q + 1)m ± 1 appear for counter-polarized CP laser drivers and orders (q − 1)m ± 1 in case of co-polarized drivers, which is a natural consequence of a density oscillation driven at (q ± 1)mω 0 . These results are confirmed both by fluid model results and PIC simulations. For the case when the frequency of the second pulse is twice as high as the first pulse, q = 2, these selection rules lead to suppression of every third harmonic when the pulses are counter-polarized but permit the generation of all harmonics for co-polarized pulses. These findings are consistent with the selection rules given by Sharma et al [44] for low-order harmonics 3, 4 and 5.
We note that although the harmonic generation mechanism in plasmas is physically distinct from that responsible for gas harmonics, the selection rules appear to be the same [32]. Further simulations show that these selection rules also apply to the case of an overdense plasma slab at normal laser incidence, confirming the results obtained previously [43]. However, according to our analytical model, harmonic selectivity in plasma is a natural consequence of the non-linear current composition, and is not something that requires symmetry effects and/or corresponding conservation laws as argued by Chen [43].
Due to strong dependence of the harmonic power on the plasma density, P m /P 0 ∝ (n e /n 0 ) 2m , the harmonic efficiency falls off rapidly for higher order harmonics which makes their experimental detection quite challenging. For a femtosecond laser with total power of 1 TW and a plasma density of n e = 0.05n c , the number of photons decreases from~10 9 for the fourth and fifth harmonics to~10 4 for seventh and eighth, and next higher orders are accordingly at least four orders of magnitude lower (n e /n 0 ) 2 ∼ 10 −4 . Therefore, in a single-shot mode, a UV spectrometer may only be able to pick up lower order harmonics. One the other hand, the CP harmonic efficiency could potentially be increased by making use of one or more of the phase-matching schemes noted above in section 5.
Although a rigorous extension of our model to overdense plasmas is a topic for future work, we note that the selection rules summarized in table 1 should still apply to a thin overdense slab or near-critical plasma targets, the main difference being that there will also be reflected as well as transmitted CP harmonics.
The non-linear fluid model introduced here already offers a fast, hi-fidelity tool for studying harmonic generation in a regime difficult to access with conventional PIC codes, although we still observe some quantitative discrepancies with the latter. For a resolution of 2 nm in a simulation box of 120 µm the 1D non-parallel fluid model currently implemented in Python takes one order of magnitude less computation time comparing to the 1D PIC code. The same 2D simulation with 2 nm resolution in laser propagation direction and 20 nm in transverse direction consume three orders of magnitude more computation time than the 1D PIC simulation. Good agreement with analytical theory in the 1D limit paves the way for extending this approach to two or even three spatial dimensions, which would enable phase-matching schemes to be investigated in more realistic optical geometry. ( for m ∈ [1, 2, 3, …], where (q + 1)m − q is the correction to (q + 1)m + 1 and (q + 1)m − 1 is the correction to (q + 1)m + q. In fact in equations (17)- (20), we have only included the main contributions for the higher order harmonics except the driver harmonics, i.e., θ and θ ′ . The same considerations apply for the co-polarized pulses, whose corrections will be, respectively: for m ∈ [1, 2, 3, …] and a −1 ≡ a 1 . For an underdense plasma, these terms are smaller than the leading terms by a factor of ∼ (ω p /ω 0 ) 2 for a chosen laser amplitude. Therefore, these corrections have been omitted when comparing the analytical results with numerical results.