Modifying angular and polarization selection rules of high-order harmonics by controlling electron trajectories in k-space

Recent advances in generation of strong laser pulses have enabled the acceleration of electrons in solids into regions far away from the band edge. Because nonlinear currents can be generated by laser-driven carriers in the non-parabolic region, tailored laser fields may allow control of optical properties of high-order harmonics (HHs). So far, investigations on laser-induced nonlinear optical phenomena have focused on the simple electron motion induced by linearly or elliptically polarized fields. However, more complex trajectories can be important for development of novel optoelectronic devices. Here, we show that a weak laser field (optical frequency is ω2) applied in a direction orthogonal to a strong main field (ω1) enhances certain HH intensity components of bulk GaSe by a factor of 100. Good agreement between the experiments and calculations shows that manipulation of the electron trajectory allows breaking inversion symmetry of the electronic states felt by the accelerated electrons and leading to a modification of selection rules for frequency-mixing processes of HHs. Owing to our usage of non-integer multiples for ω1 and ω2, it is found that the generation of HHs constitutes a novel way of ultrafast control of light polarization and optical switching.


I. Crystal angle dependences of generation of odd-order HHs for weak excitation and generation of even-order HHs
The angle dependences of both harmonic orders with m + n = even that originate from interband polarization (Figs. 1a-d) and orders with m + n = odd for weak excitation (Figs. 1e,f), can be explained by perturbative nonlinear optics [1]. Figures 1a-d show the φ dependences of the 1 HH and 2 HH components for (1,1) and (2,2), which is representative of the harmonic orders with m + n = even. Regarding the even order (1,1), 1 HH vanishes at φ = 30 and 90° (Fig. 1a), and 2 HH vanishes at φ = 0, 60, and 120° ( Fig.   1b). To explain the periodicity for m + n = even in terms of conventional nonlinear optics, the nonlinear polarization ( + ) with a nonlinear susceptibility χ can be used. In general, ( + ) reads [2] ( + ) = 1⋯1 ⏟ 2⋯2 ⏟ 1 2 , ( = 1,2). The behaviors of components with m + n = odd observed under weak excitation conditions can be explained by perturbative nonlinear optics. Figures 1e and 1f show the measured  dependences of 1 HH and 2 HH for the (1,2) order as a representative for the group of (odd,even) orders. It can be seen that under a weak applied field (E1 = 5 MV/cm),

II. Non-perturbative region observed in the intensity dependence
One of the characteristics of HHa is that, above a certain applied field (E1), the intensity dependence does not follow the scaling law of nonlinear polarization expected from perturbative nonlinear optics. This threshold is often regarded as a characteristic value that separates perturbative and non-perturbative regimes. However, this interpretation is not generally valid. For example, the intensity dependences in Fig. 1c (E1 up to 10 MV/cm) in the main text follow the scaling law. Nevertheless, we observe a nonperturbative behavior in the angle dependence of order (1,2) already at E1 = 10 MV/cm. Furthermore, the angle dependence does not change near this threshold value. Figure 2a shows an extended view of the E1 dependences in Fig. 1c  suggests that the HHa mechanism does not change within this range of electric field strengths. Because a non-perturbative behavior can be observed in the angle dependence even under weaker applied electric fields, it can be said that the angle dependence ( Fig. 3 in the main text) is more sensitive to the band anharmonicity than the intensity dependence ( Fig. 1c in the main text). In addition, we performed the experiment with the higher field of E2, and confirmed that the sample was destroyed already at E2 ≈ 3 MV/cm.
Immediately before that, the confirmed angle dependences were similar to those shown in Fig. 3 in the main text.

III. Derivation of analytical formula
To calculate the electric field generated by the 2D electron motion in the band, we start from the semi-classical model [4,5]. The temporal evolution of the electron wave packet driven by the input fields 1 ( ) and where is the elementary charge, ℏ is the Dirac constant, and the initial state is The group velocity is defined as g ( ) = (1/ℏ)( ( )/ ) using the energy dispersion ( ) of the band in which the wave packet moves. The intraband current intra ( ) can be obtained from the product of g ( ) and the electron density distribution ( , ) , i.e., intra ( ) = − ∫ g ( ) ( , ) . By integrating over the whole Brillouin zone under the assumption of a delta-like electron density ( , ) = ( − ( )) for the single electron, a simple expression for the intraband current, intra ( ) = − g ( ( )), can be deduced. The electric field of the HH, HH ( ), which is generated by the intraband current, can be derived from the time derivative of intra ( ): The component of 2 ( )/ 2 is a rank-two tensor representing the curvature of the band, and Eq. (3) can be expressed using each tensor component as shown below.
( 1 cos( 1 ) 2 cos ( 2 ) ) , In order to calculate this curvature tensor analytically, the band dispersion of the first conduction band of GaSe was modeled in the plane using a cosine function: By fitting the band structure calculated by DFT to this cosine function, we obtained Eg = 1.9 eV, a = 13 Å, and b = 0.5 eV. When 1 ( ) and 2 ( ) are applied in the directions that are tilted by the rotation angle φ with respect to the and axes, respectively, the relationship between , and 1 , 2 becomes ( , ) = ( 1 cos − 2 ( ) It can be seen that all three equations contain the same expression for the cosine function with the square brackets. By substituting Eq. (2) into this expression and expanding it in a series of Bessel functions of the first kind, we can rewrite this common expression as follows [7]: Here, = /ℏ defines the Bloch frequency, and i,j is the Kronecker delta. Eq.
(10) describes a sum of discrete oscillation terms and contains only terms for which the sum of the orders of the two frequency components [(2M+1)±(2N+1)] is even. Using this expansion and Eqs. (4) and (7-9), the generated harmonic electric field can be expanded as follows.  Figure 3a shows the crystal rotation angle dependences obtained using two-color laser fields with parallel polarization ( 1 HH is shown on the left, and 2 HH on the right side).

IV. HHG using two parallel polarized laser fields
The employed maximum electric field amplitudes inside the sample are E1 = 10 MV/cm and E2 = 0.5 MV/cm. All harmonic orders that satisfy m + n = odd have a six-fold rotation symmetry for 1 HH and a twelve-fold rotation symmetry for 2 HH . These angle and polarization dependences are same as those obtained under single-color excitation, e.g.
the same dependence as that of (5,0). In order to calculate these results analytically as in the case of orthogonal polarization, the term 1 cos( 1 ) + 2 cos ( 2 ) is substituted for the term 1 cos( 1 ) and the value zero is substituted for the term 2 cos( 2 ) in Eq. (4). The left and right panels of Figure 3b plot the 2D maps of the prefactors that Supplementary Figure 3. Crystal angle dependences of high-order harmonics with parallel electric fields E1 and E2. a, Crystal rotation angle dependence of the 1 HH and 2 HH components of the HHG spectrum obtained using two-color laser fields with parallel polarization (E1=10 MV/cm, E2=0.5 MV/cm), which are experimentally obtained. b, Calculated angle dependence of the prefactors for the harmonic orders (1,2) and (2,1). Each plot is normalized.
reflect the electric field amplitude in Eqs. (14-17), that is, , for the (1,2) and (2,1) orders, respectively. The theoretical results accurately reproduce the experiment. Similar to the case of single linearly polarized excitation, we find that, when electrons oscillate along a one-dimensional trajectory, the sum of the orders m and n determines the angle dependence.  [8]. It should be noted that the HHG spectrum of the total current includes the cross term from interband and intraband currents, 2 [̃i nter

VI. Calculation for the other orders and timing jitter
In this section, we briefly discuss the possible effects of coherence/timing jitter between the two excitation pulses described by E1(t) and E2(t). In our experiment, the about 12 cycles of E1(t) and 23 cycles of E2(t) are included within the pulse duration of 100 fs. Our laser pulses have random time delays (a small timing jitter is present between the two pulses due to, for example, thermal drift, although a common pump laser was employed) and thus the experimental results are the values averaged over the timing jitter.
To confirm that this timing jitter does not significantly affect our reported results for the angular dependence, we calculated the angular dependences for different phase offsets as shown in Fig. 4 (phase delay 0 is the timing of the field peaks of the two pulses overlap).
The results show that the symmetry of the angular dependence does not change for a different timing between the pulses and we were able reproduce the behaviors of the experimentally observed angular dependences. This data reproduces the experimentally observed angular selection rules. We can interpret this important result as follows: because the optical frequencies of the two pulses are incommensurate, even if the phase of one pulse shifts, the peak-to-peak situation occurs at another timing (within the 100 fs). symmetry with respect to (t). Therefore, the selection rule resulting in the 12-fold and 6-fold rotational symmetries for orders like (4,1) for the interband current is the same as that for the intraband current. Here, the energy ratio of E1 = 10 MV/cm to E2 = 0.5 MV/cm is 400:1. Note that the peaks of (4,0), (5,0), and (6,0) (that is, HHs that only include the frequency of E1) are generated only in the component 1 HH , and the sum frequency peaks of (2,1), (3,1) and (4,1), which are at energies close to those of the (4,0), (5,0), and (6,0) orders, respectively, are generated only in the component 2 HH . In particular, the peak intensity of (2,1), which is the strongest order in the E2 direction, is 1.5 times stronger than that of (4,0). This means that, if the frequency of E2 is twice the frequency of E1 and both fields are mutually orthogonal, the HH generated at frequency 4ω1 in the E2 direction can be stronger than that in the E1 input pulse, it can be converted into linearly polarized light by choosing an appropriate phase delay, and the polarization direction of the electric field can be rotated by about 50°

VIII. Efficiency of conversion via the orthogonal current
with respect to the E1 direction. This constitutes one approach to realize polarization control using the HHa process.