Terahertz continuum generation in the LCS lattice

Rabi oscillations in two-level Dirac systems have been shown to alter the frequency content of the system's nonlinear response. In particular, when considering Rabi oscillations in a quantum model beyond the semiclassical Boltzmann theory, even harmonics may be generated despite the centrosymmetric nature of these systems. This effect magnifies with increasing excitation intensity. In this work, we extend the Rabi theory to a three-level Dirac system arising from a line-centered-square optical lattice. In this case, the Dirac cones are bisected at the Dirac point by a flat band that persists throughout the Brillioun zone. Due to the presence of this flat band, we expect a significant enhancement of the coupling between Dirac states, resulting in a large increase of the Rabi effects and the associated nonlinearities, leading to continuum generation of terahertz radiation.


Introduction
Due to the novel electronic [1], thermal [2], mechanical [3] and optoelectronic [4] properties it possesses, graphene has drawn a considerable amount of attention in many fields. Its Dirac band structure contributes to the optical response of graphene, giving rise to novel effects including efficient light absorption and a large nonlinear susceptibility. The harmonic response of the current density in graphene, based on the semiclassical Boltzmann theory and the centrosymmetric structure of graphene itself, was believed to be restricted to only odd harmonic spectra. However, Lee et al [5] demonstrated that even in the presence of the centrosymmetry, the nonlinear optical response in graphene is not limited to odd harmonic spectra when the electron dynamics of Rabi oscillations are included in the current response [6]. Rabi oscillations can shift the current-induced harmonic spectra of graphene, and thus lead to a peak at the spectral position of the second harmonic of the incident light or even continuum generation (for a review of the Rabi effect in three-level atomic systems, see e.g. [7,8] and references therein). In the current paper, we extend the work of [5] to the three-band [9] LCS or Lieb [10,11] lattice where the additional band serves to resonantly enhance the nonlinearity giving rise to the Rabi spectrum.
The line-centered-square (LCS) lattice (figure 1(a)) formed by using ultracold atoms trapped in an optical lattice, was first described by Shen et al [11]. It is a two-dimensional counterpart of the face-centered-cubic lattice, and contains a single Dirac cone in the energy spectrum with an additional flat band energy state touching at the Dirac point. Thus, the LCS lattice exhibits three bands, i.e. the conduction band, the valence band and a central flat band ( figure 1(b)). The band structure of the LCS lattice satisfies a three-component quantum equation for pseudospin 1 Fermions.
In this paper, we calculate the Rabi contribution to the harmonic spectra of massless Dirac fermions in LCS lattice for various parameter sets. Following Lee et al [5], we begin with the time-dependent Dirac equation (TDDE) to model LCS electron dynamics. By transforming the TDDE to the dipole gauge (DG) [12], we obtain the Rabi frequency of the LCS lattice, and use it to calculate the Rabi contribution to the current response. We demonstrate that, in analogy with graphene, the nonlinear optical response of LCS lattice is not restricted to odd harmonic spectra when Rabi oscillation contributes strongly to the current response. We find the first nine orders of the harmonic spectra, ranging from 2 THz to 18 THz, fuse into a continuous spectrum in the presence Original content from this work may be used under the terms of the Creative Commons Attribution 3.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. of a 2 THz incident light field due to the nonlinearity inherent in the Rabi effect. Further, odd harmonics up to n=25 exhibit significant energy content for reasonable pump energy levels. The spectral broadening and continuum formation for the three-band system is enhanced relative to the two-band system due to the presence of a dispersionless (flat) band which separates the conduction and valence Dirac bands and provides a resonant enhancement of the Rabi effect.
The paper is organized as follows. In section 2, we begin with the TDDE in the DG to obtain the Rabi frequency and current response of the LCS lattice due to the terahertz (THz) light field. In section 3, we elucidate the Rabi coupling effect on harmonic spectra. Finally, we draw conclusions in section 4.

Eigenstates of the unperturbed Hamiltonian
In the LCS lattice, the Dirac Hamiltonian of a 2D massless electron can be described [11] as: where v 0 =±ta/ÿis the Fermi velocity, a is the LCS lattice constant, t is a hopping parameter, k x and k y are the perturbation components of k from Dirac point, and = | | k k . The eigenstates of H 0 are given by , and  + v k k x y 0 2 2 , and eigenvectors: The orthonormality of these eigenstates is given by ò , with ¢ = s s , 1, 2, 3 indexing the pseudospin. We also note here that the pseudospin indices 1, 2, and 3 correspond to the flat, valence, and conduction bands respectively.

TDDE in the dipole gauge with terahertz optical pump
We define the incident terahertz optical pump field in the Coulomb gauge by the vector potential: where E 0 is the peak of the incident electric-field strength, ζ is the polarization factor, ω 0 =2π f 0 is the central frequency and τ is the pulse width. With this incident light field, the Dirac Hamiltonian may be written: , S ± are the spin-1 ladder operators, k is the momentum operator, and J 3 is the 3×3 exchange matrix. Using this Hamiltonian, we obtain the TDDE in the Coulomb gauge: For a normally-incident THz optical field, the solution to equation (5) may be expressed as: , e x p t 6 D This represents a gauge transformation from Coulomb to dipole gauge where is the gauge generating function for the transformation, and Ψ D (r, t) is the spinor wavefunction in the dipole gauge, which we expand in the eigenstates of H 0 as: Substituting Ψ(r, t) into equation (5) yields the TDDE in the dipole gauge: is the incident electric field and I 3 is the 3×3 identity matrix. In our calculations, E(t) is assumed to be linearly polarized along x axis; thus and ζ=0. Substituting equation (7) into (8) , and using the real-space orthogonality property of the eigenstates, we obtain: Integrating equation (9) over the ¢ k space and using the transformations: e x p ,a n d e x p where ω p =v 0 k is the eigenfrequency. Similarly, following the procedure above with y ¢ Solving the coupled equations equations (10)-(12) numerically using a pulsed THz pump, we can observe the dynamical behavior of the system, including Rabi oscillations. To gain a detailed understanding of this behavior, we first consider resonant excitation with k=ω 0 /v 0 , directional angle of momentum f p =π/4, and light field intensity I=1 kW cm −2 . We assume that the electron population is initially in the valence band. Pump parameters f 0 =2 THz and τ=1 ps are used in our calculations. Figure 2 shows the temporal evolution of the normalized band populations | ( )| c t

Induced current density in LCS lattice
In order to simplify the notation used, at this point we transform the coefficients˜( ) . By introducing the continuity equation 2 , we obtain expressions for the single-particle current density , where ν=x, y indicates the current density component, as follows: 32 31 31 12 12 Finally, the net (total) current density is obtained by integrating the single-particle current density over the momentum space: where g s =2 is the spin-degeneracy factor.

Results and discussion
3.1. Single-particle harmonic spectra with Rabi oscillations In figure 3, we plot the single-particle current spectra excited by a Gaussian (figure 3(a)) or square ( figure 3(b)) pulse for resonant excitation at k=ω 0 /ν 0 and f k =π/4 with pump irradiance I=500 W cm −2 . The resultant single-particle current spectra are plotted in figures 3(c) and (d) respectively. These spectra exhibit energy at both even and odd harmonics. Further, the bandwidth over which significant energy content exists is quite large, ranging over at least the lowest five harmonics of the 2 THz fundamental excitation frequency. The detailed pulse shape is shown to impact the spectrum in a manner that is to be expected. The square pulse gives rise to sharper peaks in the spectrum than does the Gaussian pulse.

Total current spectral content
Following equation (14), we obtain the total current by integrating numerically over k. We have performed this calculation for pump irradiances ranging from I=100 W cm −2 to I=10 kW cm −2 . In what follows, we characterize the spectral content of the total current over this complete pump range, however in the interest of conserving space, we plot the results only for the I=10 kW cm −2 case in figure 4.
In figures 4(a)-(d), we illustrate the magnitude and direction of the single-particle current density over a range of times relative to the pump pulse. Figure 4(a) illustrates the thermal equilibrium current distribution prior to the arrival of the pump pulse. Figures 4(b)-(d) show the single-particle current density for times 800 fs, 6 fs, and 0 fs prior to the peak irradiance of the pump pulse respectively. These arrival times are noted as red dots on figure 4(e), which plots the temporal evolution of the total current density J. We note that in our model, the total current density does not decay exactly to zero as the current dynamics disappear due to the passing pump pulse. Such a result is a consequence the fact that our model does not include a relaxation term. The absence of a relaxation term, coupled with the fact that the area of the I=10 kW cm −2 pulse does not return the system exactly to its initial condition, results in a persistent offset in the plot of the temporal current density. Such an offset does not materially affect our conclusions.
Finally, in figure 4(f), we plot the spectrum of the total current density J. The spectrum exhibits a continuum that persists up through the ninth harmonic. This continuum rolls off by two orders of magnitude at the highfrequency cutoff. At higher frequencies, the energy is localized around the odd harmonics as would be expected from a perturbation analysis of the problem due to the symmetry inherent in the LCS lattice. Each of the odd harmonics exhibits a bifurcation into a higher and lower frequency lobe surrounding the odd harmonic. This large spectral continuum may be understood by considering the evolution equations, equations (10)- (12). In this set of equations, the Rabi frequencies are proportional to q ( ) E t k sin x k . As a result, the frequency content of the total current ranges from a component with infinite frequency at k=0 to components with frequency asymptotically approaching 0 as  ¥ k . These components are weighted proportionally to 2π k due to the increasing density of states as k increases.
Examining the resultant spectral content of the total current density for lower peak pump irradiances, we observe the following: for a peak pump irradiance of 1 kW cm −2 , the continuum spectrum is down by two orders of magnitude relative to the low-frequency component of the total current density. Harmonics 1 (fundamental) and 3 bifurcate, whereas harmonics 5-11 are visible but do not bifurcate. For a peak pump irradiance of 100 W cm −2 , the continuum is below the low frequency components by approximately the same two orders of magnitude, however only the fundamental frequency component bifurcates. Frequency components at harmonics 3-9 exist, but do not bifurcate at this irradiance.
The bifurcation and continuum formation of the spectral content of the total current may be understood by examining the form of the driving field in equations (10)- (12). The nonlinear response to the incident field gives rise to a Rabi frequency of the form where θ k is the angle of the k-vector relative to the x axis for each single-particle current component j(x, y). The Rabi frequency depends on k in two significant ways: (1) there is an angular q ( ) sin k dependence, and (2) there is a radial 1/k dependence. These two factors contribute Rabi frequency components ranging from 0 to ¥ to the total current spectrum. We further note that due to the mirror symmetry of the single-particle current density in the direction normal to the applied pump polarization, the total current in the direction normal to the polarization is zero.
When the incident THz pump field is circularly polarized, the current density and spectrum differ from the linearly-polarized case. Writing the total field in the form where E x and E y are harmonic fields in temporal quadrature with each other, equations (10)-(12) exhibit a k-dependent driving term of the form (k y E x (y)−k x E y (t))/k 2 . Therefore, the magnitude of the Rabi frequency does not change as a function of k-vector angle (the circularly-polarized E field driving the Rabi oscillation has a constant magnitude), and as a result, the angular dependence of the Rabi frequency is absent. This results in a continuum Figure 3. Illustration of the contribution from Rabi oscillations to the single-particle current spectra at k=ω 0 /v 0 with f p =π/4 in the LCS lattice. A Gaussian pulse (a) with peak irradiance (c) I=500 W cm −2 . A box-shaped pulse (b) with peak irradiance (d) I=500 W cm −2 . spectrum that persists only up through the second harmonic, for the circularly-polarized case (see figure 5). In addition, the presence of the two-component driving field results in current density components in both x and y directions. In the interest of economy of space, we plot only the x component in figure 5. The y component and its spectral content are similar. . Current density for a peak pump irradiance of I=10 kW cm −2 . Plots (a) through (d) illustrate the magnitude and direction of the single-particle current density j at increasing time samples throughout the pump pulse. Plot (e) shows the total current density J with the red dots indicating the temporal position of the single particle current density frames (a)-(d). Finally, plot (f) shows the spectral power density of the total current density J.

Conclusion
In conclusion, we have analyzed the three-level LCS lattice and obtained closed-form expressions for the carrier dynamics in this system under the influence of a picosecond THz pump pulse polarized in the plane of the lattice. The total current density J arising in this system exhibits a spectral content that evolves toward a continuum at relatively moderate pump irradiances. We provide a detailed analysis of that continuum for a pump irradiance of 10 kW cm −2 .