Floquet DMFT Analysis of High Harmonic Generation in Mott Insulators

We study the high harmonic generation (HHG) in Mott insulators using Floquet dynamical mean-field theory (DMFT). We show that the main origin of the HHG in Mott insulators is the doublon-holon recombination, and that the character of the HHG spectrum differs depending on the field strength. In the weaker-field regime, the HHG spectrum shows a single plateau as in the HHG from gases, and its cut-off energy $\epsilon_{\rm cut}$ scales linearly with the field strength $E_0$ as $\epsilon_{\rm cut}=\Delta_{\rm gap} + \alpha E_0$, where $\Delta_{\rm gap}$ is the Mott gap. On the other hand, in the stronger-field regime, multiple plateaus emerge and the $m$-th cut-off scales as $\epsilon_{\rm cut,m}=U + m E_0$. We show that this difference originates from the different dynamics of the doublons and holons in the weak- and strong-field regimes. We also comment on the similarities and differences between HHG from Mott insulators and from semiconductors. This proceedings paper complements our recent work, Phys. Rev. Lett. 121, 057405 (2018), with additional results and analyses.


Introduction
High harmonic generation (HHG) is a nonlinear phenomenon originating from strong light-matter coupling. While HHG has been initially studied in gases and is technologically important [2,3], recent observations of HHG in semiconductors open the possibility of HHG in condensed matter systems (CMs) [4][5][6][7][8][9][10][11]. HHG signals from condensed matter can be more intensive than those from gases, which is useful for developing stronger lasers. The understanding of the HHG mechanism may also lead to spectroscopic techniques which yield information on the material and its electron dynamics.
Recently, HHG from CMs beyond semiconductors were observed in amorphous systems [12,13] and liquids [14,15]. These experiments demonstrate the possibility of HHG from correlated systems where the single-particle picture is not sufficient anymore and raise questions about the role of correlations in HHG in CMs. In order to understand the effects of correlations on HHG and explore new candidates for HHG sources, in this work, we theoretically study the HHG from Mott insulators (MIs), which is a prototypical class of insulating states in CMs. By means of the single-band Hubbard model and the nonequilibrium extension of dynamical mean-field theory (DMFT), we reveal the characteristic features of the HHG in MIs, the main process for HHG, and the similarities and differences between MIs and semiconductors.

Model and Method
In this work, we focus on the Hubbard model at half-filling driven by an AC electric field, where c † i,σ is the creation operator of an electron at site i with spin σ, J i j indicates the hopping parameter, and U is the on-site interaction. The effect of the external electric field is treated by the Peierls substitution with a vector potential Here q is the electron charge. We note that A(t) is related to the electric field by E(t) = −∂ t A(t) and the setup is equivalent to a pure scalar potential in the Hamiltonian in a different gauge. Furthermore, we take account of the effects of the environment by attaching a thermal bath H bath to the Hubbard model. Here we use a thermal bath of noninteracting electrons (the Büttiker model) [1,16,17]. Due to the coupling to a thermal bath, when the system is continuously excited by an external field with frequency Ω, it reaches a time-periodic nonequilibrium steady state (NESS) with a period T ≡ 2π Ω . In the following, we calculate the NESS using the Floquet DMFT [1,[16][17][18][19]. As in the usual DMFT, we map the original lattice problem to an effective impurity problem so that the impurity selfenergy and the Green's function are identical to the local self-energy and the local Green's function of the lattice problem. The Green's functions are formulated on the Keldysh contour since the NESS is determined by the balance between the excitation and the energy dissipation to the bath [18,19]. To use DMFT, we focus on a hyper-cubic lattice with lattice spacing a in the limit of infinite spatial . We apply the field along the body diagonal, A(t) = A(t)e 0 with e 0 = (1, 1, · · · , 1) and A(t) = A 0 sin Ωt. Hence the field strength along a given axis is E(t) = −A 0 Ω cos Ωt ≡ −E 0 cos Ωt. The Büttiker model with a finite band width W bath is used, −ImΣ R bath (ω) = Γ 1 − (ω/W bath ) 2 , and we set q, a = 1 and use J * as the unit of energy.
We consider parameters for which the Mott gap is larger than the Hubbard bands, typically U = 8, β = 2.0, Γ = 0.06, W bath = 5, Ω = 0.5. Here β is the inverse temperature of the bath. This parameter choice justifies the use of the non-crossing approximatoin (NCA) and allows us to focus on the fundamental aspects of HHG in MIs. The excitation frequency is chosen such that it is much smaller than the gap, as in the semiconductor experiments. The HHG intensity is evaluated from the square of the Fourier transform of the dipole acceleration d dt j(t) as I hh (nΩ) = |nΩ j(nΩ)| 2 with n ∈ N.
Since we consider time-periodic NESSs, we only have a nonvanishing current at nΩ with n ∈ N.

Results
In Fig. 1, we show the behavior of HHG spectra from the MI for various field strengths. When the field is strong compared to the hopping (J * ), one can observe multiple plateaus, see Fig. 1(a). On the other hand, when it is comparable to or weaker than J * , there is only one plateau as in atomic gases, see Fig. 1(b). The HHG spectrum is plotted in the plane of the harmonic frequency (nΩ) and the field strength (E 0 ) in Fig. 1(c). As the field is increased, the cut-off energies of the plateaus in the HHG spectrum increase. Both in the weaker-field regime and the stronger-field regime, the cutoff energy (ǫ cut ) scales linearly with the field strength. More specifically, in the weaker-field regime ǫ cut ≃ ∆ + αE 0 with α > 1, while in the stronger-field regime we find ǫ cut ≃ U + mE 0 with m ∈ N.
The first question we address is the main process responsible for the HHG in the MI. In semiconductors, strong fields induce photo carriers (electrons in the conduction band and holes in the  Comparison of the HHG intensity extracted from the current originating from different processes in the stronger-field (a) and weaker-field (b) regimes. "Full" corresponds to the full FDMFT calculation, j rc is the current from the recombination processes, and j hop is the current from the hopping processes of doublons and holons, which are evaluated from the generalized tunneling formula presented in Ref. [19]. (c) Time-averaged local spectral functionĀ(ω) of the nonequilibrium steady state as a function of E 0 . valence band). The current in the semiconductors originates from two different processes, namely the intraband and the interband current [20,21]. The former corresponds to the dynamics of electrons and holes within the corresponding bands, while the latter is associated with the recombination of an electron and a hole. In the MI, strong fields induce doubly occupied sites (doublons) and empty sites (holons), which act as photo carriers. As in the semiconductor case, the current can be separated into two contributions, the doublon/holon hopping process and the doublon-holon recombination process. The former process does not alter the number of doublons and holons (the carrier number) and is analogous to the intraband current in semiconductors. On the other hand, the latter process changes the number of doublons and holons and is analogous to the interband current in semiconductors. Now the question is which process is dominant. By means of the generalized tunneling formula derived in Ref. [19], one can evaluate the currents from these two different processes. The results are shown in Fig. 2(a,b). In all regimes the recombination current dominates over the hopping current. Therefore, we conclude that both in the weak-and strong-field regimes recombination is the dominant process for HHG. We also note that, in some regimes, the contribution from the hopping process roughly follows that of the recombination process. This is because, rigorously speaking, these processes cannot be fully decoupled since the states with different doulbon-holon numbers is slightly mixed via the finite transfer integral. (A similar effect has been reported in a semiconductor study under electric fields [20].) In Fig. 2(a) for nΩ > 20, we can see the values of these contributions become closer and cancel each other (destructive interference).
The next question is the origin of the different behavior depending on the field strength. To find out the origin, we first show the (time-averaged) local spectrum of the NESS,Ā(ω) ≡ − 1 π ImḠ R loc (ω), see Fig. 2(c). Here G R loc indicates the retarded part of the local Green's function, andḠ R loc (ω) = 1 T T 0 dt av dt r G R loc (t r ; t av )e iωt r , where t r denotes the relative time and t av the average time. One can see that for E 0 J * , there is only a small renormalization of the Hubbard bands, while for stronger E 0 , they are strongly renormalized and one can observe Wannier-Stark side peaks. Note that the width of the bands indicates the mobility of injected test charges. Hence these features indicate that under the strong fields (E 0 J * ) carriers tend to be localized. These observations suggest the following picture. In the stronger-field regime, after the doublons and holons have been created, they cannot move freely because of the field-induced localization. At some point, they recombine by emitting light, which leads to the U + mE(t) cutoff for the localized pairs of doublons and a holons separated by m sites. On the other hand, in the weaker-field regime, after the doublons and holons have been created, they can move more freely because of the not so strong field and acquire kinetic energy. In the recombination, they emit the minimum energy necessary to create a doublon and holon plus the kinetic energy.
In the following, we provide additional data and analyses to support these scenarios. For the stronger-field regime, we perform a sub-cycle analysis to obtain a characteristic frequency emitted at each time (t probe ) within one cycle [21]. More specifically, we calculate a windowed Fourier transformation of j(t), j(ω; t probe ) = dte itω j(t)W(t; t probe ), and evaluate I hh (ω; t probe ) ≡ |ω j(ω; t probe )| 2 . Here W(t; t probe ) is the Blackman window function with a half-window of length 2 centered at t = t probe . As is expected from the scenario of recombination of localized doublon-holon pairs, the signal follows U ± E(t) in the normal-scale plot Fig. 3(a). This corresponds to the recombination of the nearestneighbor pairs. In the log-scale plot, Fig. 3(b), the signal suddenly becomes weak for |ω| |U ±2E(t)|, and in some time domain a clear strong signal is seen on the |U ± 2E(t)| curve. The latter correspond to the recombination of next-nearest-neighbor pairs. These results support the scenario of localized pair recombination in the stronger-field regime. We also note that this scenario can naturally explain the multiple plateaus and their scaling. The m-th neighbor pairs yield a signal at U ± mE(t), which oscillates from U − mE 0 to U + mE 0 . Hence, the m-th plateau, whose edge scales as U + mE 0 , originates from the recombination of m-th neighbor pairs.
To analyze the weaker-field regime, we studied the cutoff energies for different values of U. According to the scenario of the recombination of itinerant doublons and holons mentioned above, we expect that the emitted energy is the sum of the minimum energy necessary to create the doublonholon pair (the Mott gap energy in the present case) and the kinetic energy which the carriers acquired under the influence of the periodic field (E kin ). In gas systems, the former corresponds to the ionization energy and the latter corresponds to the ponderomotive energy [2,3]. Indeed, the offset ∆, determined from extrapolations E 0 → 0, essentially coincides with the Mott gap (∆ gap ), which is determined from the local spectrum independently. Hence, it is consistent with the above scenario and our numerical results indicate that the acquired kinetic energy scales linearly against the field, E kin = αE 0 . Now we compare the HHG from MIs and semiconductors. We consider a semiconductor model with a valence band and a conduction band, which corresponds to the upper and lower Hubbard bands, respectively: D α denotes the band center for band α = {v, c}, and we choose D v = −U/2 and D c = U/2. We introduce a transfer integral between the different semiconductor orbitals at neighboring sites, in order to mimic the fact that the hopping of electrons in MI leads to the creation of a doublon/holon pair on neighboring sites. The effect of the electric field is included via the Peierls substitution and we consider the NESS by attaching a Büttiker-type thermal bath. We consider two sets of hopping parameters, (i) J c = J v = J cv = 0.5J (we call this "type 1" model) and (ii) J c = −J v = J cv = 0.5J (we call this "type 2" model). The type 1 model corresponds to the indirect gap semiconductors, whose dispersion is the same as that from the Hubbard 1 approximation (H1A), which is based on the atomic-limit self-energy Σ R (ω) = U 2 4ω [22]. Namely, for a momentum k, the energy of the upper band and lower band is expressed as ǫ k,± = (ǫ k ± ǫ 2 k + U 2 )/2 with the free particle energy ǫ k = −2J d a=1 cos(k a ). Indeed, the dispersion of the single-particle spectrum obtained from DMFT follows nicely the dispersion of the type 1 model (or the H1A), see Fig. 4(a). On the other hand, the type 2 model corresponds to a direct gap semiconductor. The HHG spectrum of the type 1 semiconductor is shown in Fig. 4(b). Intriguingly, the structure of the HHG spectrum turns out to be very similar to that from the H1A, and one can observe cutoff energies that scale with U + E 0 and U + 3E 0 . However, the model strongly underestimates the HHG spectrum in the weak to intermediate field regimes, because electron-hole pairs are not efficiently created. On the other hand, interestingly, the type 2 model shows qualitatively similar features (scaling and intensity distribution) as the MIs, although its dispersion is very different from that of MIs, see Fig. 4(a). These comparisons demonstrate that the different nature of MIs and semiconductors is reflected in a very different relation between the HHG spectrum and the single-particle spectrum, and intuition from the semiconductor HHG does not work when one just looks at the single-particle dispersion. We note that we have also studied the type1 semiconductor with disorders. The correlation from impurities produces a finite width in the spectrum at each k and makes the spectrum even closer to that of MIs, but such correlation effects do not result in a HHG spectrum which more closely resembles that of MIs.  Fig. 1 and Ω = 0.5, U = 8, β = 2.0, Γ = 0.06.

Conclusion
In this proceedings, we analyzed the HHG from Mott insulating states described by the single band Hubbard model using the Floquet dynamical mean-field theory. We showed that the main origin of the HHG is the doublon-holon recombination and that the character of the HHG spectrum changes due to the different dynamics of the doublons and holons in the weaker-and stronger-field regimes. We also discussed the similarities and differences between the HHG in Mott insulators and semiconductors. Recent experiments have motivated numerous theoretical efforts to understand the HHG in correlated systems such as disordered systems [23,24], strongly correlated systems [25][26][27] as well as spin systems [28]. We hope that our present work and these studies will stimulate further explorations of new HHG sources in consented matters.