Conventionally, nonlinear media are employed to transform the fundamental frequency of the input laser to those of high harmonics to obtain extreme ultraviolet and even X-ray beams using table-top lasers^[1^,2^]. High harmonic generation (HHG), as an ideal tool, can be applied in probing electronic dynamics on the atomic or molecular scale, microscale imaging, etc. Gases or metals are utilized to generate high-order harmonics for lasers with intensities around 10^14–16 W/cm^2[3^–7^], and a plasma target is ideal for lasers with much higher intensities, i.e., approximately or over 10¹⁸ W/cm^2[8^–15^]. In general, all harmonics couple with the input laser(s), and they are not separated outside the interaction domain^[16^-20^]; consequently, additional spectrometers are usually utilized to separate a specific-order harmonic for application^[21^-23^]. Angular separation of high harmonics from solid plasma targets has been achieved experimentally in a single-beam scheme by using grating targets^[24^]. Enhancement of harmonics in the grating scheme via surface plasmon excitation has also been demonstrated^[25^]. On the other hand, angularly isolated high harmonics from a gas jet were first reportedly generated by two overlapping 10¹⁴ W/cm² non-collinear laser beams^[22^]. This non-collinear HHG geometry presents numerous advantages^[22^,26^-29^], including the separation of the harmonics from the pump beams, angular separation of different order harmonics, generation of high-energy photons, and production of isolated attosecond bursts. Evidently, this HHG process meets the phase-matching requirement, which stipulates that the output momentum of an nth-order harmonic photon (k_n) is only related to the momentum of the input photons (k₁ and k₂) by the momentum conservation principle, k_n = n₁k₁ + n₂k₂, where n₁ and n₂ are the numbers of photons absorbed from each of the input beams, and n = n₁ + n₂^[22^,23^,30^]. As far as we know, all this interesting research is based on the gas ionization process, where the laser intensity is generally below 10¹⁶ W/cm²_, and the crossing angle is small enough to ensure adequate interaction.

In this study, based on a different HHG mechanism, we propose an approach for generating intense isolated harmonics through the interaction between relativistic lasers and a plasma target, but with a quite different phase-matching selection rule. In our approach, two relativistic laser pulses impinge on a solid thin foil at an angle. The generated intense harmonics are naturally separated from each other and from the input beams at considerably large angles, which can be predicted. The proposed approach offers several exclusive benefits. First and most importantly, both our theoretical analysis and three-dimensional particle-in-cell (3D PIC) simulations prove that the nth harmonic carrying the information of both input lasers is converted by its adjacent (n – 1)th harmonic which has the smallest crossing angle with the nth harmonic owing to the highest occurrence possibility and one fundamental input beam, rather than directly by the two input beams in the gas HHG case. That is, considering phase matching, the output momentum (emission direction) of an nth-order harmonic photon depends on the momentum of the (n – 1)th-order harmonic photon and that of one input laser photon, according to the momentum conservation, k_n = k_n–1 + k₁₍₂₎, for the input beams with the same wavelength. Second, plasma, as a nonlinear medium, can endure lasers nearly without the intensity limitation. Therefore, the intensity of isolated harmonics can be in the relativistic region^[31^,32^]. Finally, because the crossing angle between the input lasers can be sufficiently large, the individual harmonic can be emitted at an angle nearly two orders of magnitude higher than that (usually in the milliradian level) in gas HHG. This point is very important for the measurement and application of this method in experiments. It should be noted that a large number of papers on the selection rule for the plasma HHG in terms of the oscillating mirror model, first proposed by Bulanov et al.^[9^] and then studied extensively by Lichters et al.^[10^], have been published. All of them are about a single incident laser or collinear incident lasers, and all of harmonics are coupled spatially. The problems such as proper phase matching do not exist. However, in our case, the harmonic emission is determined by the phase-matching relation and hence isolated angularly.

The proposed scheme is shown in Figure 1. Two circularly polarized (CP) laser pulses a₁ and a₂ impinge on a solid thin foil symmetrically at a large crossing angle 2θ. For simplification, two identical input pulses with the same frequency (ω) are considered, and the angle θ = π/4 is the incident angle of each laser pulse. The electrons in the target experience the overlapped fields in the thin foil and oscillate. The harmonics radiated from the oscillating electrons are isolated angularly and propagate in the reflected side and also in the transmitted side if the target is sufficiently thin. The total normalized laser amplitude acting on the target is a = a₁+a₂ ($a= eA/{m}_e{c}^2$, where A is the vector potential, $c$ is the speed of light in vacuum, ${m}_e$ is the electron mass, and $e$ is the electron charge), where a₁ and a₂ are the driving laser beams used as

$$\begin{align}{a}_1&={a}_0[\sin (\omega t- kx\cos \theta - ky\sin \theta )\notag\\{}&\times (\overset{\frown }{y}\cos \theta -\overset{\frown }{x}\sin \theta )+\cos (\omega t- kx\cos \theta - ky\sin \theta )\overset{\frown }{z}],\notag\\{}{a}_2&={a}_0[\sin (\omega t- kx\cos \theta + ky\sin \theta )\notag\\{}&\times (\overset{\frown }{y}\cos \theta +\overset{\frown }{x}\sin \theta )\pm \cos (\omega t- kx\cos \theta + ky\sin \theta )\overset{\frown }{z}],\end{align}$$ (1)

$$\begin{align}{a}_{\mathrm{CR}\perp}&\sim \sin (\omega t- kx\cos \theta )\notag\\{}&\times [\cos \theta \cos ( ky\sin \theta )\overset{\frown }{y}+\sin ( ky\sin \theta )\overset{\frown }{z}],\notag\\{}{a}_{\mathrm{SR}\perp}&\sim \cos ( ky\sin \theta )\notag\\{} &\times [\cos \theta \sin (\omega t- kx\cos \theta )\overset{\frown }{y}+\cos (\omega t- kx\cos \theta )\overset{\frown }{z}],\end{align}$$ (2)

As an example, we adopt the z component of the input fields, ${a}_z={a}_{1z}+{a}_{2z}$, which is related to the harmonics, to analyze the HHG process. Here ${a}_{1z}=\cos \left(\omega t- kx\cos \theta - ky\sin \theta \right)\overset{\frown }{z}$ and ${a}_{2z}=-\cos (\omega t- kx \cos \theta + ky\sin \theta )\overset{\frown }{z}$. The nonlinear part, ${a}_z/\gamma$ (where $\gamma =\sqrt{1+{a}^2}$), in Maxwell’s equation, ${\nabla}^2a-{\partial}^2a/\partial {t}^2={na}_z/\gamma$, contributes to the harmonic generation (where n is the particle density). After the Fourier expansion of this source term, we can derive the expressions for the harmonics generated, as follows:

$$\begin{align}\left\{\cos \left[2m\left(\omega t- kx\cos \theta \right)+\omega t- kx\cos \theta - ky\sin \theta \right]\right.\notag\\{}\left.+\cos \left[2m\left(\omega t- kx\cos \theta \right)+\omega t- kx\cos \theta + ky\sin \theta \right]\right\}\overset{\frown }{z},\end{align}$$ (3)

$$\begin{align}\tan {\alpha}_n=\frac{\sin \theta +\left(n-1\right)\sin {\alpha}_{n-1}}{\cos \theta +\left(n-1\right)\cos {\alpha}_{n-1}},\notag\\{}\left({k}_{nx},{k}_{ny}\right)=\left(n\cos {\alpha}_n,n\sin {\alpha}_n\right).\end{align}$$ (4)

The PIC simulation configuration is shown in Figure 2(a). The simulations are carried out using the EPOCH code^[35^]. Two CP laser pulses, a₁ and a₂, with the same frequency ω (the corresponding wavelength λ=1 μm) reach the target at the same angle θ=π/4 symmetrically. A peak normalized amplitude (a₀) of 3 and a focal spot size of 5 μm are used. The laser profile is ${\sin}^2\left[\pi t/\left(2{t}_0\right)\right]$, where t₀ = 7T and T is the period of the input laser. The simulation box is $20\;\unicode{x3bc} \mathrm{m}\;(x)\times 40\;\unicode{x3bc} \mathrm{m}\;(y)\times 40\;\unicode{x3bc} \mathrm{m}\;(z)$, corresponding to a window with $1000\times 1000\times 1000$ cells and one particle per cell. A thin foil with a thickness of $1\;\unicode{x3bc} \mathrm{m}$ and a density of ${n}_0=20{n}_c$, is employed, where ${n}_c=1.1\times {10}^{21}\ {\mathrm{cm}}^{-3}$ is the critical density for the input laser pulse. Here we note this thin foil is plasma which is initially solid and ionized by the laser pedestal. At $t=0$, the laser pulses enter the simulation box. Considering that the harmonics in the reflected direction are more intense, we focus our analysis on the harmonics in the reflected side from the simulations. As shown in Figure 2(b), where the fundamental frequency components are filtered out, harmonics containing the information of two input pulses are emitted in the reflected side, and they are much more intense than those in the same-rotation case (see Figure s1 in the Supplementary Material).

Figure 3 shows the spectrum distribution of the harmonics in k-space after the Fourier transformation of the reflected electric field, E_z, in Figure 2(b). As expected, the harmonics include odd and even orders in the reflected directions. The harmonic locations marked with the small white circles in Figure 3 are obtained from Equation (4), showing the good agreement with the simulation results. Harmonics such as 2ω₂₀, 3ω₃₀, 4ω₄₀, 2ω₀₂, 3ω₀₃, and 4ω₀₄ only contain the information of one input laser pulse and propagate along the incident directions. Here we focus on the harmonics containing the information of both input laser pulses. They are generated and angularly isolated, including 2ω₁₁ of the second harmonic, 3ω₂₁ and 3ω₁₂ of the third harmonics, and 4ω₂₂, 4ω₃₁, and 4ω₁₃ of the fourth harmonics, as shown in Figure 3. The propagating directions of these harmonics are determined by vector addition according to their generation way, in which the phase-matching rule should be satisfied.

For example, there are four ways to generate the third harmonics in the photon picture according to energy conservation, i.e., 3ω=3×1ω (a₁), 3ω=3×1ω (a₂), 3ω=2ω+1ω (a₁), 3ω=1ω (a₂) +2ω. However, it implies there are four different emission directions, as shown with 3ω₃₀, 3ω₀₃, 3ω₂₁, and 3ω₁₂ in Figure 3, because of the different phase-matching relation. Harmonics 3ω₃₀ or 3ω₀₃, are undoubtedly emitted along the incident directions because only one linear momentum is related. However, for the other two ways, it is questionable that the photon with an energy of 2ω originates from 2ω₂₀, 2ω₀₂, or 2ω₁₁. Taking 3ω₂₁ as an example, each photon of 3ω₂₁ may be transformed by one photon of a₁ and one photon of 2ω₁₁, or by one photon of a₂ and one photon of 2ω₀₂ (this selection rule is the same as that in non-collinear gas HHG). For the former, according to Equation (4), the coordinate of 3ω₂₁ in the k-spectrum distribution should be (2.9, 0.76), and the emission angle (α) should be 0.25 rad. Conversely, for the latter, the coordinate of 3ω₂₁ in the k-spectrum distribution should be (2.85, 0.95), and the emission angle should be 0.32 rad. The former selection rule is evidenced by the PIC simulations, as shown in Figure 3. In addition, according to our chain selection rule, each photon of 4ω₃₁ should be transformed by one photon of a₁ and one photon of 3ω₂₁, and its coordinate in the k-spectrum distribution is expected to be (3.71, 1.5) and the emission angle 0.37 rad, instead of the (3.58, 1.79) coordinate and the 0.464 rad emission angle in the gas HHG way (one photon of a₂ and one photon of 3ω₃₀). This expectation is also confirmed by the simulation result in Figure 3.

To further verify this conclusion, we show the second harmonic (2ω₁₁), third harmonic (3ω₂₁), and the fourth harmonic (4ω₃₁) in Figure 4. Figures 4(d)–4(f) demonstrate the electric field distributions of the constant phase planes, where the emission angle is selected according to Equation (4). The distributions meet our expectations. We also note that the intensities of the relatively low-order harmonics are in the relativistic range in the present a₀ = 3 case, with the normalized amplitude around 1, basically following the high energy conversion efficiency rule in plasma HHG. The emission angle is about two orders of magnitude higher than that (usually, in the milliradian level) in the non-collinear gas HHG.

Here, two points are elucidated. According to the angular momentum conservation, the n = 2m order harmonics in the normal direction should be linearly polarized, because one photon of such harmonics stems from m photons of a₁ and m photons of a₂, and the total spin angular momentum is zero. The other harmonics are circularly (or elliptically) polarized, because one photo of such harmonics is transformed by one photon of its adjacent (n – 1)th-order harmonic and one photon of the input laser (CP laser). In fact, the PIC simulation result confirms this point; see the k-space spectrum distribution of E_y in the Supplementary Material where the even harmonics in the normal direction are missing. We also compared these results with those obtained in the cases of purely linearly polarized incident beams. In the case of p-polarized laser case, harmonics include both the odd and even orders in y component. In the case of s-polarized case, there are only odd order harmonics in z component and even order harmonics in y component. These results are in agreement with the conclusion in Ref. [10]. The emission angle for the same order harmonic is identical with that in the CP incident lasers case because it is determined by the same selection rule. Second, phase matching ensures that there is no wave vector mismatch between the emitted harmonics and the input beams which contributes to its emission, usually in the case of small enough crossing angle between the input beams, and these harmonics at different locations can hardly be distinguished spatially. However, in the present large crossing angle geometry, the projection of the input beam wave vectors along the harmonic emission direction should be considered, which results in a phase mismatch, expressed as $\Delta k={k}_n-{k}_{n-1}\cos \left({\alpha}_n-{\alpha}_{n-1}\right)-{k}_{1(2)}\cos \left(\theta -{\alpha}_n\right)$. As demonstrated in Figures 1 and 3, the harmonic wave vector is slightly larger than the sum of the wave vector projections of the beams which contributes to the harmonic. In this case, the phase matching must be accomplished by considering the whole interaction system including the plasma.

For the HHG from the laser–plasma interaction, one of the most significant findings is that the roll-off of the harmonic spectrum exhibits a power-law dependence of I_ωn ~ ω^-8/3, and the cutoff frequency increases with the laser intensity by ω_cut-off ~ I_L^2.5[11^,34^]. Our results basically follow these findings for the same plasma HHG mechanism. Our selection rule is still valid for the different intensities and crossing angles, confirmed by additional simulations in the cases of a₀ = 0.2, 0.5, 3, 10, 20, 50 when the half crossing angle is π/4 and a₀ = 3 when the half crossing angle is π/8 and π/3. That is, for the same crossing angle, the intensity and cut-off frequency of the harmonics rise with the input laser intensity. On the other hand, for the same input lasers, the superposed electric field acting on electrons in the surface layer will be strengthened with the decrease of the crossing angle, which will also enhance the harmonics and extend the cut-off frequency. According to the simulation results, we found the conversion efficiency follows the power-law dependence of ~ω^-4.1 as shown in Figure 5(a) for the harmonics in the reflected directions. For the same order harmonics emitted at different angles, the conversion efficiency (intensity) decreases with the emission angle, as shown in Figure 5(b).

Although the emission angles of the harmonics obtained from the phase-matching equation agree well with the simulation results, we find a slight angle mismatch, as shown in Figures 3 and 4. The main reason is the deformation of the surface plane and Doppler effect of the reflected laser pulse. If the highly relativistic laser pulse is used, the radiation pressure dominates in the laser–plasma interaction, and there will be a curve for the surface electron layer, which may deflect the reflected pulses and impact the emission angle slightly. At the same time, red shift effect will occur for the lasers (harmonics) reflected from the election layer. We know that the motion velocity of the surface electron layer increases with the laser intensity and decreases with the target density, i.e., ${v}_p\propto a/\sqrt{n}$^[36^–39^]. For the same target density, the velocity of the surface electron layer increases with the laser intensity, which broadens the spectrum of the harmonics. However, for the given high-density target, our theoretical model is still predictable for the highly relativistic laser case, though the harmonic spectrum would be slightly modulated, as shown by the comparison of the results of different laser intensity cases in Figure 6.

Harmonics can also be modulated angularly directly from laser irradiated gratings, which has been demonstrated theoretically and experimentally. The approach using a single-beam scheme with grating targets results from the target (grating) modulation directly. In our approach, by overlapping two non-collinear lasers on the plasma target, the spatial structure of the driving field is imprinted on the electron dynamics and, hence, on the harmonic generation. The generated harmonics contain the information of the two incident laser pulses, and thus one can adjust the incident lasers for the different harmonic emission. Most importantly, we found the new phase-selection rule, which is the feature for the non-collinear plasma HHG mechanism, and the harmonic emission direction can be obtained. In addition, the non-collinear plasma HHG process is not limited to mixing two input beams of the same wavelength. When two-color input pulses are used, the selection rule is more complicated because there are additional electron oscillation frequencies except 1ω and 2ω which determines the HHG process. The selection rule of the phasing matching to generate harmonics is similar to what we obtained.

The computational precision in the PIC simulation may also affect the simulation results, particularly for those of order higher than the fourth. We check this by conducting another two-dimensional (2D) PIC simulation with 10,000×10,000 grids instead of 1000×1000 for the x×y window in the aforementioned 3D simulation. Apparently, there is nearly no difference between the harmonics in this case and those analyzed earlier (see the k-space spectrum distribution in the 2D case in the Supplementary Material).

In summary, we have proposed an approach, based on the non-collinear plasma HHG process, for generating angularly isolated harmonics at large emission angles. When two counter-rotation CP lasers hit a plasma target at the same side in V-shape, the harmonics carrying the information of the two input lasers are isolated. Using the momentum vector addition of photons, the propagating direction has been obtained. With the present scheme, the origin of these harmonics radiated in the relativistic laser and plasma interaction was elucidated for the first time. Both theoretical analysis and 3D PIC simulations confirmed the newly derived phase-matching selection rule for generating high harmonics. Thus, it is meaningful for practical applications and for understanding the plasma HHG process.

Multi-angle imaging is a potential application for these angularly isolated harmonics. Conventionally, to obtain more detailed imaging information, more than one laser pulses with different wavelengths are required to probe the sample simultaneously from different directions. The application of multiple pump lasers in an optical band will challenge the experimental layout and temporal synchronization. Therefore, angularly isolated multi-wavelength lasers generated simultaneously using the non-collinear HHG approach present significant application promise, particularly in biomedical imaging and material detection. In addition to such application as the gas harmonic source, one can extend its application in laser–plasma interaction, taking advantage of the high intensity of the plasma harmonic source. Although the plasma harmonic source has not efficiently been used so far owing to the strict requirement for the laser condition (intensity, contrast, etc.), there are some exclusive advantages for this harmonic source, such as high intensity/flux and large crossing angle, and this will certainly broaden the application field. Moreover, with the development of the laser technique and the experimental skill, the plasma harmonic source application is expected to flourish in the near future.