# Anomalous temperature dependence of high-harmonic generation in Mott
insulators

Yuta Murakami  Department of Physics, Tokyo Institute of Technology, Meguro,
Tokyo 152-8551, Japan  Center for Emergent Matter Science, RIKEN, Wako,
Saitama 351-0198, Japan  Kento Uchida  Department of Physics, Graduate School
of Science, Kyoto University, Sakyo-ku, Kyoto 606-8502, Japan  Akihisa Koga
Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551,
Japan  Koichiro Tanaka  Department of Physics, Graduate School of Science,
Kyoto University, Sakyo-ku, Kyoto 606-8502, Japan  Institute for Integrated
Cell-Material Sciences, Kyoto University, Sakyo-ku, Kyoto, 606-8501, Japan
Philipp Werner  Department of Physics, University of Fribourg, 1700 Fribourg,
Switzerland

###### Abstract

We reveal the crucial effect of strong spin-charge coupling on high-harmonic
generation (HHG) in Mott insulators. In a system with antiferromagnetic
correlations, the HHG signal is drastically enhanced with decreasing
temperature, even though the gap increases and the production of charge
carriers is suppressed. This anomalous behavior, which has also been observed
in recent HHG experiments on Ca2RuO4, originates from a cooperative effect
between the spin-charge coupling and the thermal ensemble, and the strongly
temperature-dependent coherence between charge carriers. We argue that the
peculiar temperature dependence of HHG is a generic feature of Mott
insulators, which can be controlled via the Coulomb interaction and
dimensionality of the system. Our results demonstrate that correlations
between different degrees of freedom, which are a characteristic feature of
strongly correlated solids, have significant and nontrivial effects on
nonlinear optical responses.

High-harmonic generation (HHG) is a fundamental nonlinear optical phenomenon
with potentially important technological applications. It was first reported
in atomic gases Ferray  _et al._ (1988) and is utilized in attosecond laser
sources as well as spectroscopiesKrausz and Ivanov (2009). Recently its scope
is extended to condensed matters because of the observation of HHG in solids,
in particular semiconductors and semimetals Ghimire  _et al._ (2011); Schubert
_et al._ (2014); Luu  _et al._ (2015); Vampa  _et al._ (2015a); Langer  _et
al._ (2016); Hohenleutner  _et al._ (2015); Ndabashimiye  _et al._ (2016); Liu
_et al._ (2017); You  _et al._ (2017); Yoshikawa  _et al._ (2017); Hafez  _et
al._ (2018); Kaneshima  _et al._ (2018); Yoshikawa  _et al._ (2019); Cheng
_et al._ (2020); Schmid  _et al._ (2021). HHG in semiconductors and semimetals
can be well described by the dynamics of independent electrons (independent-
particle picture) Golde  _et al._ (2008); Vampa  _et al._ (2014, 2015b); Wu
_et al._ (2015); Otobe (2016); Ikemachi  _et al._ (2017); Tancogne-Dejean  _et
al._ (2017); Luu and Wörner (2016); Hansen  _et al._ (2017); Osika  _et al._
(2017); Ikeda  _et al._ (2018); Tamaya  _et al._ (2016); Floss  _et al._
(2018); Lysne  _et al._ (2020a); Chacón  _et al._ (2020); Wilhelm  _et al._
(2021); Taya  _et al._ (2021), which enables the HHG spectroscopy of band
information such as dispersion relations Vampa  _et al._ (2015c); Luu  _et
al._ (2015); Li  _et al._ (2020); Luu and Wörner (2018); Uchida  _et al._
(2021). On the other hand, the effects of electronic correlations are often
taken into account phenomenologically and a detailed understanding of their
role in solid-state HHG is lacking Kemper  _et al._ (2013); Kruchinin  _et
al._ (2018); Ghimire and Reis (2019); Du and Ma (2022). This understanding is
however essential for the exploration of HHG and the application of HHG
spectroscopy in correlated materials.

Figure 1: (a) Experimental HHG intensity at the indicated HHG peaks as a
function of the optical gap for Ca2RuO4 (Mott insulator) and InAs
(semiconductor), reproduced from Ref. Uchida  _et al._ , 2022. The temperature
$T$ is modified in the range $T\in[290$ K, $50$ K]. (b) DMFT results for the
intensity at the indicated HHG peaks as a function of the Mott gap
($\Delta_{\rm Mott}$) for the single-band Hubbard model in the Mott insulating
phase.

The new research frontier of HHG in strongly correlated systems (SCSs) has
attracted considerable interest both on the theoretical Silva  _et al._
(2018); Murakami  _et al._ (2018); Murakami and Werner (2018); Lysne  _et al._
(2020b); Tancogne-Dejean  _et al._ (2018); Imai  _et al._ (2020); Chinzei and
Ikeda (2020); Orthodoxou  _et al._ (2021); Murakami  _et al._ (2021); Shao
_et al._ (2022); Hansen  _et al._ (2022); Masur  _et al._ (2022); Udono  _et
al._ (2022) and experimental Grånäs  _et al._ (2020); Bionta  _et al._ (2021);
Uchida  _et al._ (2022) sides. In contrast to semiconductors, which can be
described in terms of electrons and holes, the driven state of SCSs involves
various types of many-body elemental excitations. This makes the mechanism and
features of HHG in SCSs nontrivial. Previous studies revealed the direct
connection between many-body excitations and HHG in SCSs Murakami  _et al._
(2018); Imai  _et al._ (2020); Murakami  _et al._ (2021), suggesting possible
spectroscopic applications of HHG to detect many-body states Murakami  _et
al._ (2021) as well as photoinduced phase transitions Silva  _et al._ (2018).
On the other hand, very recently, an unexpected exponential enhancement of the
HHG signal with increasing gap size is reported in the Mott insulator Ca2RuO4
Uchida  _et al._ (2022), see Fig. 1 (a). Since a larger gap should suppress
the excitation of charge carriers, this increase is opposite to the behavior
expected in semiconductor HHG. Such a counter-intuitive result calls for a
deeper theoretical understanding of HHG in SCSs. A hallmark of SCSs is the
coupling between different degrees of freedom, such as charges, orbitals and
spins. These correlations are at the origin of rich physical properties
observed in equilibrium SCSs Imada  _et al._ (1998); Dagotto (1994). However,
their role in highly nonlinear optical phenomena such as HHG is hardly known.

In this letter, we reveal the crucial role of spin-charge coupling on HHG in
Mott insulators analyzing the Hubbard model. Previous works showed that HHG in
Mott insulators originates from the coherent dynamics of a pair of local many-
body states – a doublon (doubly occupied state) and holon (empty state) –
generated by strong fields, where the three-step model picture is applicable
Murakami  _et al._ (2018, 2021). The kinematics of doublons and holons is
strongly correlated with spins, since their hopping disturbs the spin
background. We demonstrate that this spin-charge coupling and its cooperation
with thermal fluctuations produces a drastic enhancement of the HHG intensity,
accompanied with an increasing Mott gap, as observed in Ca2RuO4 (Fig. 1(b)).
These results demonstrate that strong correlations between active degrees of
freedom in SCSs can result in counter-intuitive behaviors of highly nonlinear
optical phenomena such as HHG.

We focus on the single-band Hubbard model, which is a standard model for SCSs.
The Hamiltonian is

$\displaystyle\hat{H}(t)=-t_{\rm hop}\sum_{\langle
ij\rangle}e^{i\phi_{ij}(t)}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$
(1)

where $\hat{c}^{\dagger}_{i\sigma}$ is the creation operator for an electron
with spin $\sigma$ at site $i$, $\langle ij\rangle$ indicates a pair of
neighboring sites, and
$\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}$. $t_{\rm
hop}$ is the hopping parameter and $U$ the onsite interaction. The electric
field is included via a Peierls phase $\phi_{ij}$, see Supplemental Material
(SM) 111Supplemental Material [url], which includes Ref. Werner  _et al._
(2018).. We mainly use the nonequilibrium dynamical mean-field theory (DMFT)
Georges  _et al._ (1996); Eckstein and Werner (2010); Aoki  _et al._ (2014);
Sandholzer  _et al._ (2019); Schüler  _et al._ (2020) to solve this problem,
and focus on the Bethe lattice for simplicity Werner  _et al._ (2017). The
qualitatively same results are obtained for the two-dimensional square
lattice, see Note (1). In the following, we use the quarter of the bandwidth
at $U=0$ as the energy unit, and mainly consider $U=6$. If our energy unit
corresponds to $0.5$ eV, the Mott gap ($\Delta_{\rm Mott}\simeq 3$, see below)
corresponds to 1.5 eV. This is a typical gap size of cuprates, which are often
described by the Hubbard model.

We consider the half-filled system, which becomes a Mott insulator for large
enough $U$ in equilibrium. While the Mott insulator can be realized in the
paramagnetic (PM) phase, the system on the bipartite lattice exhibits an
antiferromagnetic (AF) phase below the Néel temperature $T_{c}$ $(\simeq
0.15)$. The corresponding evolution of the single-particle spectra is shown in
Fig. 2(a). With decreasing temperature $T$, the Mott gap $\Delta_{\rm Mott}$
increases. In the PM phase, the upper and lower Hubbard bands are featureless.
On the other hand, in the AF phase, peak structures develop within the bands,
indicating the formation of spin-polarons Martinez and Horsch (1991); Dagotto
(1994); Sangiovanni  _et al._ (2006). When an electron is added to (removed
from) the system, a doublon (holon) is created, see Note (1) for schamatics.
When this doublon (holon) moves around, it can disturb the spin background at
the cost of multiples of the exchange energy $J_{\rm ex}$ $(=\frac{4t_{\rm
hop}^{2}}{U})$. This results in strong spin-charge coupling, of which the
spin-polaron is one manifestation.

Figure 2: (a) Local spectral functions, $A_{\rm loc}(\omega)$, in equilibrium.
(b) HHG spectra of the Mott insulator computed with DMFT for various $T$. (c)
The intensity at the peaks of the HHG spectra as a function of $T$. The peak
intensity is normalized by the value at $T=0.2$ (PM phase). For (a-c), we use
$U=6$. (d) $U$-dependence of the increase ratio of the HHG peaks. In order to
take into account the change of the Mott gap, we compare the $(2U+n)$th HHG
peaks. We use $T=0.2$ for the PM phase, while we use $T=0.3/U$ for the AF
state to take account of the change of $J_{\rm ex}\propto\frac{1}{U}$. The
excitation parameters are $E_{0}=0.8$, $\Omega=0.5$, $t_{0}=75$ and
$\sigma=15$.

Now we discuss the kinematics of doublons and holons accompanied by a
disturbance of the spin configurations, and its effect on highly nonlinear
optical phenomena. We study the $T$-dependence of HHG in Mott insulators
excited with frequency $\Omega$ smaller than the Mott gap $\Delta_{\rm Mott}$.
We mainly use $\Omega=0.5$ in the following. If our energy unit corresponds to
$0.5$ eV, this is a mid-infrared excitation with $0.25$ eV, whose period
$T_{p}$ is about $16$ fs. From the $T$-dependence of the spectral functions,
one would naively speculate that the HHG intensity is suppressed by lowering
temperature, since the enhancement of the gap reduces the tunneling
probability (see Note (1)) and the formation of the spin-polarons suggests a
reduced mobility of the charge carriers. However, the $T$-dependence turns out
to be opposite to this naive expectation.

Applying a Gaussian electric field pulse $E(t)$, we evaluate the HHG intensity
$I_{\rm HHG}(\omega)$ from the Fourier transformation of the current $J(t)$ as
$I_{\rm HHG}(\omega)=|\omega J(\omega)|^{2}$. The pulse is characterized by
the standard deviation $\sigma$, the center $t_{0}$ and the maximum field
strength $E_{0}$. We show the resulting HHG spectra for various temperatures
in Fig. 2(b) and plot the $T$-dependence of the relative intensity of the HHG
peaks in Fig. 2(c). $I_{\rm HHG}(\omega)$ is strongly enhanced above
$\Delta_{\rm Mott}$ and the width of the HHG plateau is enhanced with
decreasing temperature. The increase in the ratio of HHG signals is larger for
the higher harmonic peaks. Above $T_{c}$, the $T$-dependence becomes very
weak. As a function of the gap, the intensity increases almost exponentially,
as illustrated in Fig. 1(b). Importantly, the DMFT results of the simple
Hubbard model reproduce the qualitative features of the HHG spectrum and the
empirical scaling law observed in Ca2RuO4 Uchida  _et al._ (2022) (see Fig.
1(a) and Note (1)).

To reveal the origin of this $T$-dependence, we consider the $U$-dependence of
the relative increase of the HHG signal (Fig. 2(d)). For large $U$ the
bandwidth of the upper and lower Hubbard bands is insensitive to $U$ and the
Mott gap scales almost linearly with $U$. Therefore, in order to focus on the
contribution from the kinetic energy of the doublon-holon pair, we compare
$I_{\rm HHG}(\omega)$ for the same $\omega-U$. It turns out that the increase
ratio monotonically decreases with increasing $U$. Since $J_{\rm ex}$ is
reduced with increasing $U$, the disturbance of the spin background costs less
energy, and the spin-charge coupling becomes weaker. Hence, the $U$-dependence
of the HHG increase ratio suggests that the anomalous $T$-dependence of HHG is
related to the spin-charge coupling.

Next we perform a subcycle analysis considering a windowed Fourier transform
$J(\omega,t_{p})=\int dte^{i\omega t}F_{\rm window}(t-t_{p})J(t)$ and
evaluating $I_{\rm HHG}(\omega,t_{p})\equiv|\omega J(\omega,t_{p})|^{2}$.
$I_{\rm HHG}(\omega,t_{p})$ provides the time-resolved spectral features of
the emitted light around $t_{p}$. Since HHG in Mott insulators mainly
originates from the recombination of doublon-holon pairs Murakami  _et al._
(2018, 2021), the subcycle spectra reveal the recombination time of the pairs
and their energy at that time. In Figs. 3(a),(b), we show $I_{\rm
HHG}(\omega,t_{p})$ in the PM and AF phases. In both cases, the dominant
intensity appears at early times within one period, suggesting that only short
trajectories of the doublon-holon pairs contribute to the HHG signal. In other
words, the coherence time of the doublon-holon pair is very short ($<T_{p}/4$)
compared to one cycle of the pulse field and to the coherence times typically
considered in the analysis of semiconductors, see e. g. Fig. 6 in Ref. Vampa
_et al._ , 2015b. The kinematics estimated from the peak position of $I_{\rm
HHG}(\omega,t_{p})$ at each $\omega$ as a function of $t_{p}$ is represented
with red dashed (blue dot-dashed) lines for the AF (PM) phase in Figs.
3(a),(b). These lines define the function $f_{\omega}(t_{p})$. The difference
between the blue and red lines is mostly explained by the difference in the
gap size ($\simeq 1.1$), indicating that the trajectory of the doublon-holon
pair is almost the same in the AF and PM phases. The main difference is the
coherence time of the pair.

Figure 3: (a),(b) Subcycle spectra $I_{\rm HHG}(\omega,t_{p})$ for $U=6$ at
(a) $T=0.2$ (PM phase) and at (b) $T=0.05$ (AF phase). A Gaussian window with
standard deviation $\sigma^{\prime}=0.9$ is used. The red dashed (blue dot-
dashed) lines indicate the maxima of $I_{\rm HHG}(\omega,t_{p})$ at $T=0.05$
($T=0.2$) at a given $\omega$ as a function of $t_{p}$ around $t_{p}=80$,
which define the function $f_{\omega}(t_{p})$. The vertical dashed lines
indicate the times when the electric field $E(t)=0$. (c) Intensity $I_{\rm
HHG}(\omega,t_{p})$ along the lines $f_{\omega}(t_{p})$ for $U=6$. (d)
Normalized intensity $I_{\rm HHG}(\omega,t_{p})$ along the lines
$f_{\omega}(t_{p})$ for the indicated values of $U$. We use $T=0.2$ for the PM
phase and $T=0.3/U$ for the AF states to take account of the change of $J_{\rm
ex}\propto\frac{1}{U}$. $I_{\rm HHG}(f_{\omega}(t_{p}),t_{p})$ is renormalized
by the value at $t_{p}=79.5$ in each case. The excitation parameters are the
same as in Fig. 2.

To quantify this, we show in Fig. 3(c) the intensity along the peaks, $I_{\rm
HHG}(f_{\omega}(t_{p}),t_{p})$. The results indeed show that for the higher
$T$ the intensity decays faster, suggesting that the dephasing time of the
doublon-holon pair is shorter. This is in a stark contrast with the behavior
of the charge distribution, where the absence of the AF spin background at
high $T$ leads to a slower relaxation Lenarčič and Prelovšek (2013); Golež
_et al._ (2014); Eckstein and Werner (2016). On the other hand, with
increasing $U$, the behavior of $I_{\rm HHG}(f_{\omega}(t_{p}),t_{p})$ in the
AF and PM phases becomes more similar, see Fig. 3(d). Furthermore, the peak in
$I_{\rm HHG}(f_{\omega}(t_{p}),t_{p})$ becomes clearer, which indicates that
the intensity coming from longer-time trajectories of the doublon-holon pairs
and hence the coherence time are increased. This feature appears counter-
intuitive, because the single-particle spectrum becomes highly incoherent for
large $U$ Dagotto (1994); Sangiovanni  _et al._ (2006), and demonstrates that
HHG in SCSs is not directly related to the single-particle spectra, in
contrast to semiconductors Murakami  _et al._ (2018); Imai  _et al._ (2020);
Murakami  _et al._ (2021).

Figure 4: (a) $I_{\text{HHG}}$ for different $B^{z}_{\text{stagg}}$ and
subcycle analysis for $B^{z}_{\rm stagg}=0.001$ (b) and $B^{z}_{\rm
stagg}=0.2$ (c). A Gaussian window with $\sigma^{\prime}=0.9$ is used. The
colored markers indicate the energy emitted at $t_{p}$ by the recombination of
a doublon-holon pair, which is predicted from the three-step model using the
doublon and holon dispersions from the Bethe ansatz Murakami  _et al._ (2021).
The color indicates the time interval between the recombination and the
creation of the doublon-holon pair $t_{\rm pair}$, and
$T_{p}=\frac{2\pi}{\Omega}$. (d) Phase of the Fourier component of $J(\omega)$
at $\omega=n\Omega$ ($n$ is an integer). In all panels, we set $U=8$, and the
excitation parameters are $\Omega=0.5,E_{0}=0.8$, $t_{0}=60$ and $\sigma=15$.

These behaviors can be consistently explained in terms of the spin-charge
coupling. To directly compare cases with and without spin-charge coupling, we
switch to the one-dimensional (1D) Hubbard model with a staggered magnetic
field $B^{z}_{\rm stagg}$. In one dimension, without $B^{z}_{\rm stagg}$, the
kinematics of the doublons and holons is independent of the spin-degrees of
freedom (spin-charge separation), while for $B^{z}_{\rm stagg}\neq 0$, the
hopping of a doublon (holon) creates a mismatch between the staggered field
and the spin configuration, as it happens in higher-dimensional systems
without field, see Note (1) for schematics. With this set-up, the 1D model can
mimic the spin-charge coupling in higher dimensions. The infinite time-
evolving block decimation (iTEBD) Vidal (2003) allows to compute accurate
results for this model at $T=0$ in the thermodynamic limit.

We show the HHG spectra for various $B^{z}_{\rm stagg}$ in Fig. 4(a) and the
corresponding subcycle analysis in Figs. 4(b),(c). For small $B^{z}_{\rm
stagg}$, the expected HHG peaks at $(2n+1)\Omega$ in $I_{\rm HHG}(\omega)$ are
not clear, suggesting that the system is not fully time periodic during the
pulse. This is attributed to the long coherence time of the doublon-holon
pair, which leads to the interference of many quasi-classical trajectories
within the three-step model Vampa  _et al._ (2015b). Indeed, the subcycle
spectra for small $B^{z}_{\rm stagg}$ suggest that long trajectories of
doublon-holon pairs strongly contribute to the HHG signal, see Fig. 4(b)
Murakami  _et al._ (2021). With increasing $B^{z}_{\rm stagg}$, the HHG
intensity becomes weaker but the HHG peaks become clearer at $(2n+1)\Omega$.
Here, $B^{z}_{\rm stagg}$ is chosen to be comparable to $J_{\rm ex}$. In the
subcycle spectrum, the weight is shifted to earlier times in one period, see
Fig. 4(c), as it is the case in the DMFT results in Fig. 3(b) at low $T$.
These results show that the coherence time of the doublon-holon pair is
efficiently suppressed by the spin-charge coupling, which consistently
explains the behavior of the DMFT results. The short coherence time reduces
the interference between different quasi-classical trajectories and results in
clear HHG peaks both in the DMFT data and the iTEBD data for nonzero
$B^{z}_{\rm stagg}$.

The reduction of the coherence time of the doublon-holon pair with increasing
$T$ can be understood as a cooperative effect of the spin-charge coupling and
the thermal ensemble. At nonzero temperatures, the initial equilibrium state
is described by an ensemble of eigenstates, represented by the density matrix
$\hat{\rho}\propto e^{-\beta\hat{H}}$. In such a system, the total current
induced by the field can be calculated as the ensemble average over the
individual currents evaluated for these eigenstates. With increasing $T$, the
weight of the high-energy states increases. In our case, at higher
temperatures, spin configurations different from the AF ground state are
activated, see Note (1). The dynamics of the doublon or holon is different for
each configuration, since the energy transfer to the spin background during an
excursion depends on the spin configuration. This should produce emitted light
with different phases for different spin configurations, resulting in phase
cancellations after the ensemble average, and thus reduce the coherence
between the doublon-holon pairs with increasing $T$. Note that this effect
does not rely on long-range magnetic ordering and is also relevant in the PM
phase, but is absent without spin-charge coupling. Namely, for small $J_{\rm
ex}$, weaker cancellations between different spin configurations are expected,
which explains the results in Fig. 3(d) and the reduction of the enhancement
of the HHG signal with larger $U$ in Fig. 2(d). To exemplify that the spin-
charge coupling can indeed provide such phase shifts, in Fig. 4(d), we show
the $B^{z}_{\rm stagg}$-dependence of the phase of $J(\omega)$ for
$\omega=n\Omega$ (with $n$ some integer). The result suggests that the phase
is sensitive to $B^{z}_{\rm stagg}$, which supports the above argument. Hence,
the modification of the coherence time due to the spin-charge coupling and
thermal fluctuations dominates over the reduction of the tunneling rate by the
gap opening, leading to an enhancement of $I_{\text{HHG}}$ at lower
temperatures.

The strong $T$-dependence of the HHG spectrum observed in Mott insulators is
not expected in typical semiconductors. In the theoretical analysis of HHG in
semiconductors, a short dephasing time $T_{2}$ of a few fs for an electron-
hole pair is often used. The main origin of the fast dephasing is the
experimental setup, i.e. the dephasing by the propagation of light and the
inhomogeneity of the field strength Floss  _et al._ (2018); Kilen  _et al._
(2020), which is insensitive to temperature. Another relevant factor is the
electron-electron scattering among excited carriers in semiconductors Becker
_et al._ (1988); Nagai  _et al._ (2021). Still, this is also expected to be
insensitive to temperature, since thermal fluctuations cannot efficiently
excite carriers across the gap. These considerations are supported by the
experimental HHG spectrum for the semiconductor InAs shown in Fig. 1(a).

In summary, our theoretical study revealed important effects of strong spin-
charge coupling on the coherent carrier dynamics in Mott insulators, which
lead to the counter-intuitive enhancement of HHG accompanied by a gap
enhancement. Spin-charge coupling is inevitable in Mott insulators in
dimensions larger than one, so that this peculiar behavior should be a generic
feature of HHG in SCSs (see Note (1)). In addition, in multiorbital systems
like Ca2RuO4, the orbital-charge coupling should have a similar effect as the
spin-charge coupling (see Note (1)) Strand  _et al._ (2017). These insights
demonstrate the important role of correlations in highly nonlinear optical
responses and provide useful guidance for the future exploration of HHG in
SCSs. On the one hand, our results suggest that the $T$-dependence of
$I_{\text{HHG}}$ can be controlled by changing the ratio $\frac{U}{t_{\rm
hop}}$, which is feasible with the application of chemical or physical
pressure. On the other hand, to realize a strong HHG signal, 1D Mott systems
are more favorable than higher-dimensional ones due to the absence of spin-
charge coupling. The recovery of coherence and the possible increase of the
HHG intensity due to the reduction of the dimensionally could be
systematically analyzed by exploiting the dimensional crossover in ladder-type
compounds such as Srn-1Cun+1O2n Imada  _et al._ (1998). Furthermore, the
sensitivity of HHG to the temperature and spin-charge coupling suggests
possible HHG-based techniques for detecting and characterizing thermal and
non-thermal phases, and for measuring the strength of the spin-charge
coupling. In the future, it will also be interesting to study HHG with more
sophisticated methods such as cluster DMFT to reveal the role of magnetic
fluctuations.

###### Acknowledgements.

The calculations have been performed on the Beo05 cluster at the University of
Fribourg. This work is supported by Grant-in-Aid for Scientific Research from
JSPS, KAKENHI Grant Nos. JP20K14412(Y.M.), JP21H05017(Y.M.,K.U.,K.T.),
JP19H05821, JP18K04678, JP17K05536 (A.K.), JST CREST Grant No. JPMJCR1901
(Y.M.), and ERC Consolidator Grant No. 724103 (P.W.). The nonequilibrium DMFT
calculations have been implemented using the open source library Nessi Schüler
_et al._ (2020).

##  I Non-equilibrium dynamical mean-field theory

Dynamical mean-field theory (DMFT) is a powerful theoretical framework that
can deal with strongly correlated systems Georges  _et al._ (1996), and in
particular Mott physics. Nonequilibrium DMFT is the extension of DMFT to
nonequilibrium problems Aoki  _et al._ (2014). The DMFT approach gives
reliable results for high-dimensional systems. This has been confirmed also
for the nonequilibrium dynamics by a recent ab-initio comparison between
nonequilibrium DMFT and cold atom quantum simulators for a three dimensional
system.Sandholzer  _et al._ (2019) DMFT is based on the Green’s function
formalism, and in nonequilibrium DMFT the Green’s functions $G(t,t^{\prime})$
are defined on the so-called L-shape contour, which includes the Matsubara
(imaginary time) branch and a real-time contour.Aoki  _et al._ (2014) The
information on the temperature enters through the Matsubara branch, which
describes the initial equilibrium system (at $t=0$) at a given temperature.

To evaluate the Green’s functions, in DMFT, we map the lattice system to an
effective impurity model with a self-consistently determined time-dependent
noninteracting bath Aoki  _et al._ (2014). In the present case, we consider
the single-band Hubbard model with a possible antiferromagnetic (AF) order,
and assume that the spins can be polarized along the $z$ axis. In the AF
phase, the two sublattices A and B show opposite magnetizations and effective
impurity models are introduced for each sublattice. The action of the impurity
model for the sublattice $\alpha$ can be expressed as

$\displaystyle\mathcal{S}^{\alpha}_{\rm imp}=-i\int
dtdt^{\prime}\sum_{\sigma}\hat{d}^{\dagger}_{\sigma}(t)\Delta^{\alpha}_{\sigma}(t,t^{\prime})\hat{d}_{\sigma}(t^{\prime})-i\int
dt\hat{H}_{\rm loc}(t),$ (2)

where $\Delta$ is the hybridization function and $\hat{H}_{\rm
loc}(t)=-\mu\sum_{\sigma}\hat{d}^{\dagger}_{\sigma}(t)\hat{d}_{\sigma}(t)+U\hat{d}^{\dagger}_{\uparrow}(t)\hat{d}_{\uparrow}(t)\hat{d}^{\dagger}_{\downarrow}(t)\hat{d}_{\downarrow}(t)$.
We note that here the times $t,t^{\prime}$ are defined on the L-shaped
contour. In DMFT, $\Delta$ is self-consistently determined such that the local
Green’s functions and the self-energies of the lattice are the same as those
of the impurity model. We solve the impurity model using the non-crossing
approximation (NCA) Eckstein and Werner (2010), which yields reliable results
in the strong coupling regime.

The lattice self-consistency condition for the hybridization function used in
the main text is
$\Delta^{\alpha}_{\sigma}(t,t^{\prime})=\sum_{\xi=\pm}\Delta^{\alpha}_{\sigma,\xi}(t,t^{\prime})$,
where $\xi=\pm$ corresponds to the positive/negative bond direction (relative
to the polarization of the field) and
$\Delta^{\alpha}_{\sigma,\pm}(t,t^{\prime})=\frac{t_{\rm hop}^{2}}{2}e^{\pm
iA(t)}G_{\rm imp,\sigma}^{\bar{\alpha}}(t,t^{\prime})e^{\mp iA(t^{\prime})}$,
with $A$ the vector potential of the field pulse Werner  _et al._ (2017).
Here, the bond length $a$ and the electron charge are set to unity. The vector
potential is related to the electric field $E(t)$ by $E(t)=-\partial_{t}A(t)$.
We choose the vector potential $A(t)=\frac{E_{0}}{\Omega}F_{\rm
G}(t,t_{0},\sigma)\sin(\Omega(t-t_{0}))$ with $F_{\rm
G}(t,t_{0},\sigma)=\exp[-\frac{(t-t_{0})^{2}}{2\sigma^{2}}]$. $E_{0}$
indicates the maximum value of the electric field. The corresponding current
(per site) can be computed as $J(t)={\rm
Im}[\sum_{\sigma,\xi=\pm}\xi\Gamma^{\alpha}_{\sigma,\xi}(t)]$, where
$\Gamma^{\alpha}_{\sigma,\pm}(t)\equiv-i[G_{\rm
imp,\sigma}^{\alpha}*\Delta_{\sigma,\pm}^{\alpha}]^{<}(t,t)$. The self-
consistency condition represents a Bethe lattice with $d$ bonds connected to
each lattice site, where we take the limit of $d\rightarrow\infty$ with a
rescaled hopping parameter $t_{\rm hop}/\sqrt{d}$. In the free system ($U=0$),
the full bandwidth becomes $W=4t_{\rm hop}$. For the Bethe lattice, the self-
consistency condition is simplified, compared to other lattices, which reduces
the numerical cost and enables a systematic analysis. The qualitative features
of the HHG spectrum are expected to be insensitive to the choice of the
lattice. Below we confirm this point with simulations for the two-dimensional
square lattice, although the scope of this analysis is limited.

In equilibrium, we define the momentum-averaged single-particle spectral
function as

$\displaystyle A_{\rm loc}(\omega)=-\frac{1}{\pi}{\rm Im}\;G^{R}_{\rm
loc}(\omega).$ (3)

Here, we introduced the retarded Green’s function
$G^{R}_{ij,\sigma}(t-t^{\prime})=-i\theta(t-t^{\prime})\langle\\{\hat{c}_{i\sigma}(t),\hat{c}^{\dagger}_{j\sigma}(t^{\prime})\\}\rangle$,
and $G_{\rm loc}=\frac{1}{N}\sum_{i}G_{ii,\sigma}$. $G^{R}(\omega)$ is defined
as $\int dte^{i\omega t}G^{R}(t)$. In practice, since the simulation is
limited to finite $t$, we use a Gaussian window function with standard
deviation $\sigma^{\prime}=8$ in the Fourier transformation. The Mott gap
$\Delta_{\rm Mott}$ is determined by the criterion $A_{\rm loc}(\Delta_{\rm
Mott}/2)=\delta$ with $\delta=0.005$. We checked that the choice of $\delta$
has no qualitative effect on the results.

Figure 5: (a) Temperature dependence of the magnetization,
$m=|\langle\hat{n}_{i\uparrow}\rangle-\langle\hat{n}_{i\downarrow}\rangle|/2$.
The vertical dashed line indicates the transition temperature $T_{c}$. We use
$U=6$ and consider the Bethe lattice. (b) Schematic pictures of the spin-
charge coupling accompanying the kinematics of a doublon (circle). Panel i)
shows the spin configuration of an antiferromagnetic state, ii) shows a
doublon added to this state, and iii) shows the dynamics of the doublon, which
disturbs the spin configuration (zigzag lines) at the cost of multiples of the
exchange energy $J_{\rm ex}$.

##  II Infinite time-evolving block decimation

The Hamiltonian of the one-dimensional Hubbard model considered here is

$\displaystyle\hat{H}(t)=$ $\displaystyle-t_{\rm
hop}\sum_{i,\sigma}[e^{-iA(t)}\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{i+1,\sigma}+h.c.]$
(4)
$\displaystyle+U\sum_{i}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}+B^{z}_{\rm
stagg}\sum_{i}(-1)^{i}\hat{S}_{z,i},$

where $\hat{S}_{z,i}=\frac{1}{2}(\hat{n}_{i,\uparrow}-\hat{n}_{i,\downarrow})$
and $t_{\rm hop}$ is set to unity. We analyze this model with the infinite
time-evolving block decimation (iTEBD) method Vidal (2003). In iTEBD, assuming
translational invariance, we express the wave function of the system as a
matrix product state (MPS). iTEBD directly treats the thermodynamic limit and
we use the cut-off dimension $D=2000$ for the MPS to obtain converged results.
In the implementation, we use the conservation laws for the numbers of spin-up
and spin-down electrons to improve the numerical efficiency.

Figure 6: Single-particle spectral functions $A_{\rm loc}(\omega,t_{\rm av})$
and $A^{<}_{\rm loc}(\omega,t_{\rm av})$ after the application of the electric
field. $A^{<}_{\rm loc}(\omega,t_{\rm av})$ measures the density of carriers
excited by the field. We set $U=6$ and the pulse parameters are (a)
$\Omega=0.5,E_{0}=0.8,t_{0}=75$ and $\sigma=15$, and (b) $E_{0}=0.8,t_{0}=75$
and $\sigma=15$. In all cases, we take $t_{\rm av}=150$.

##  III Supplementary results for the Bethe lattice

In this section, we present supplementary results obtained with DMFT for the
Hubbard model on the Bethe lattice. First, we show the magnetization in
equilibrium in Fig. 5(a). As explained in the main text, it shows the
appearance of an AF phase below the transition temperature
$T_{\text{N\'{e}el}}\approx 0.15$. Figure 5(b) illustrates the kinematics of a
doublon added to the AF state. It explains how the doublon dynamics couples to
the modification of the spin background. Panel i) shows the spin configuration
of an antiferromagnetic state, ii) shows a doublon added to this state, and
iii) shows the dynamics of the doublon, which disturbs the spin configuration
(zigzag lines) at the cost of multiples of the exchange energy $J_{\rm ex}$.

Figure 7: (a) HHG spectra $I_{\rm HHG}(\omega)$ of the Mott insulator
evaluated with DMFT for various temperatures. (b) The intensity at the peaks
of the HHG spectra $I_{\rm peak}$ as a function of the temperature. (c) The
intensity at the peaks of the HHG spectra as a function of the Mott gap size
$\Delta_{\text{Mott}}$. In panels (b) and (c), the peak intensity is
normalized by the value at $T=0.2$. Here, we set $U=6$ and consider the Bethe
lattice. The pulse parameters are $\Omega=1,E_{0}=0.8,t_{0}=75$ and
$\sigma=15$.

Next, we show the temperature dependence of the number of charge carriers
created by the field. To this end, we compute the time-dependent single-
particle spectral functions Aoki  _et al._ (2014)

$\displaystyle A_{\rm loc}(\omega,t_{\rm av})=-\frac{1}{\pi}{\rm
Im}\;G^{R}_{\rm loc}(\omega,t_{\rm av}),$ (5) $\displaystyle A^{<}_{\rm
loc}(\omega,t_{\rm av})=\frac{1}{2\pi}{\rm Im}\;G^{<}_{\rm loc}(\omega,t_{\rm
av}).$

Here, we introduced the Green’s functions
$G_{ij,\sigma}(t,t^{\prime})=-i\langle
T_{\mathcal{C}}\hat{c}_{i\sigma}(t)\hat{c}^{\dagger}_{j\sigma}(t^{\prime})\rangle$
with $T_{\mathcal{C}}$ being the contour ordering operator, $G_{\rm
loc}=\frac{1}{N}\sum_{i}G_{ii,\sigma}$, and $G^{R}$ and $G^{<}$ are the
retarded and lesser parts of the Green’s function. $G^{R/<}(\omega,t_{\rm
av})$ is defined as $\int dt_{\rm rel}e^{i\omega t_{\rm rel}}G^{R/<}(t_{\rm
rel},t_{\rm av})$, where $G^{R/<}(t_{\rm rel},t_{\rm
av})=G^{R/<}(t,t^{\prime})$, $t_{\rm rel}=t-t^{\prime}$ and $t_{\rm
av}=\frac{t+t^{\prime}}{2}$. $A_{\rm loc}(\omega,t_{\rm av})$ yields the
energy spectrum at time $t_{\rm av}$ and $A^{<}_{\rm loc}(\omega,t_{\rm av})$
the occupied states. Figure 6 shows how many charge carriers are created after
the pulse field. Panel (a) plots the results for the same condition as in the
main text, while panel (b) shows the results for a pulse without oscillations
($E(t)=E_{0}F_{\rm G}(t,t_{0},\sigma)$), i.e. close to a DC excitation. Both
panels demonstrate that a much smaller amount of charge carriers is produced
when the temperature is lowered and the Mott gap is increased, as expected.

Next, we show the supplemental data for the HHG spectra for $\Omega=1$. This
excitation frequency is twice larger than what is used in the main text, and,
in terms of the value of $\Delta_{\rm Mott}/\Omega$, this choice is closer to
the situation in the experiments on Ca2RuO4 Uchida  _et al._ (2022). In Fig. 7
(a), we show the HHG spectra for $U=6$ for different temperatures. The
corresponding temperature dependence of the intensity of the HHG peaks is
shown in Fig. 7(b), while the HHG peak intensity is shown as a function of the
Mott gap in Fig. 7(c). The intensity of the $n$th HHG peak $I_{\rm peak}$ is
defined as the maximum value of $I_{\rm HHG}(\omega)$ for
$\omega\in[(n-1)\Omega,(n+1)\Omega]$. Since now the gap size is two to three
times the excitation frequency, one may expect that multi-photon excitation
processes play an important role. Still, the results shown in Fig. 7 indicate
that the qualitative temperature dependence of the HHG spectrum is
qualitatively the same as that for $\Omega=0.5$ and that the characteristic
increase of the HHG intensity with decreasing temperature is insensitive to
the excitation condition. Furthermore, the match between the experiment for
Ca2RuO4 and the theory is better for $\Omega=1$ than for $\Omega=0.5$, compare
Fig. 7(c) with Fig. 1 in the main text.

Figure 8: The increase ratio of the HHG peak intensity for the $n$th order
harmonics as a function of the emitted energy
$\omega_{\textrm{emit}}=n\Omega$. Note that both axes are plotted on a log
scale. The vertical dashed line indicates the gap at low temperatures
($\Delta(T\rightarrow 0)$). The dot-dashed line shows the fit to the data for
$\omega_{\textrm{emit}}>\Delta(T\rightarrow 0)$, which supports the feature
(iii). Here, we set $U=6$ and consider the Bethe lattice. The pulse parameters
are $E_{0}=0.8$, $t_{0}=75$ and $\sigma=15$.

Now we discuss to what extent the empirical scaling relation for the HHG
intensity introduced in Ref. Uchida  _et al._ , 2022 for Ca2RuO4 is reproduced
by our theoretical results. The empirical relation has the form

$\displaystyle I_{\rm peak}(n,T)/I_{\rm
peak}(n,T_{0})\simeq\Bigl{(}\frac{n\Omega}{\Omega_{\text{th}}}\Bigl{)}^{\frac{\Delta_{\rm
Mott}(T)-\Delta_{\rm Mott}(T_{0})}{\Delta_{\text{th}}}},$ (6)

where $I_{\rm peak}(n,T)$ is the peak intensity of the $n$th HHG peak at
temperature $T$, $T_{0}$ is a reference temperature, and $\Omega_{\rm th}$ and
$\Delta_{\rm th}$ are fitting parameters. This equation implies the following
three features of HHG:

  1. (i)

The HHG intensity for a given order of harmonics exponentially increases with
$\Delta(T)-\Delta(T_{0})$.

  2. (ii)

The increase ratio of the HHG peak intensity is determined only by the
emission energy $\omega_{\text{emit}}(=n\Omega)$.

  3. (iii)

At fixed temperature, the increase ratio exhibits a power law as a function of
$\omega_{\text{emit}}$ with positive exponent.

These three features are the necessary and sufficient conditions for Eq. (6).
The present theory does not perfectly, but well reproduce these features in
the following sense. Property (i) is almost satisfied as pointed out in the
main text. A close inspection reveals that the theoretical results show a
slight saturation at large $\Delta$, which makes the fitting by Eq. (6) not
perfect, but the main trend is consistent. Property (ii) is nicely reproduced
by the present analysis as shown in Fig. 8, where we compare the increase
ratio for different excitation frequencies. Property (iii) is also
consistently reproduced as long as we focus on
$\omega_{\text{emit}}\geq\Delta_{\rm Mott}(T\rightarrow 0)$, see Fig. 8. Note
that in Ref. Uchida  _et al._ , 2022 the empirical law has been proposed for
the harmonics above the maximum gap of the system. In the theoretical results,
the ratio for the harmonics below the gap does not follow the scaling law (6).
This qualitatively different behavior below the gap likely originates from a
different HHG mechanism in this regime. Above the gap, HHG mainly originates
from the recombination of doublon-holon pairs. On the other hand, around or
below the gap, the contribution from the hopping of doublons or holons becomes
relevant, as in the case of the intraband current in semiconductors. The
resulting change in the HHG mechanism is the likely origin of the deviation
from the formula.

Figure 9: Schematic pictures of the interplay between thermal ensembles and
spin-charge coupling. Panels i)-iii) show cases with different spin
configurations activated by thermal fluctuations. The difference between case
i) and case ii) is indicated by the blue dashed rectangle, while that between
case i) and case iii) is indicated by the green dashed rectangle. The zigzag
lines in each panel show the sites where the mismatch in the spin
configuration occurs as the doublon moves around. The horizontal arrows with
different colors indicate the difference in the kinematics of the doublon due
to the spin-charge coupling.

We note that a systematic derivation and explanation of the empirical formula
is difficult and beyond the scope of this paper. Here we simply propose
plausible scenarios for some of the above features and leave a detailed
analysis to some future work. As for (i), the strong increase of the HHG
intensity originates from the large change in the coherence time of the
doublon-holon pairs. This is indeed an important message of this paper. The
exponential enhancement of the HHG intensity as a function of the gap size can
be explained in the following manner. It is natural to assume that the factor
in the intensity originating from the dephasing can be expressed as
$\exp(-t_{p}/T_{2})$. Here $t_{p}$ is the time interval between the creation
and recombination of a relevant doublon-holon pair, and $T_{2}$ is the
coherence time. This $T_{2}$ decreases with increasing temperature, which is
accompanies by a decrease of the gap. In this situation, if the inverse of the
coherence time behaves as $\frac{1}{T_{2}}=C_{1}-C_{2}\Delta_{\text{Mott}}$,
with $C_{1}$ and $C_{2}$ some positive constants, an exponential enhancement
of the HHG signal as a function of the gap size is expected. As for (ii), we
need to remember that the strength of the electric field ($E_{0}$) is fixed.
Therefore, when the excitation frequency ($\Omega$) is increased, the
amplitude of the vector potential ($A_{0}=E_{0}/\Omega$) is reduced, and vice
versa. The increase of the excitation frequency tends to reduce the time
interval $t_{p}$, while the decrease of the vector potential tends to increase
it. These two effects may compensate each other such that the interval $t_{p}$
for a given emission energy ($\omega_{\rm emit}$) becomes insensitive to a
change in the excitation frequency. The ratio of the interval $t_{p}$ and the
coherence time determines how much the pairs contribute to the emission. Since
this ratio is almost unchanged, the effects of the dephasing for a given
$\omega_{\text{emit}}$ should be almost independent of the excitation
frequency. Thus, given that the change of the HHG intensity mainly originates
from the change of the dephasing time, the increase ratio should also be
insensitive to the excitation frequency. As for (iii), we can only explain why
the exponent should be positive, i.e. that the increase ratio is larger for a
larger emission energy. For a larger emission energy, the corresponding
interval $t_{p}$ increases, and therefore the effect of the change in the
coherence time becomes larger. Further systematic analyses are needed to
reveal the origin of the power law. Finally, we note that, for a detailed
discussion, the efficiency of the creation of the doublon-holon pair should
also be taken into account. However, this effect should be subdominant for the
analysis of the temperature dependence of a given harmonic, and the above
discussion should cover the main effects. We note that the creation becomes
more efficient as we increase the temperature and suppress the gap, as shown
in Fig. 2 of the Supplemental Material. However, the HHG intensity is enhanced
if temperature is decreased, which suggests that the dephasing effect
dominates.

Figure 10: Local single-particle spectral functions $A_{\rm loc}(\omega)$ and
momentum-resolved single-particle spectral functions $A_{\bf k}(\omega)$ for
the Hubbard model on the two-dimensional square lattice. Panels (a) is for
$T=0.2$ in the PM phase and (b) is for $T=0.05$ in the AF phase. Here, $U=6$,
$t_{\rm hop}=0.5$ and $(L_{x},L_{y})=(16,16)$. In the momentum-resolved
spectra, $\Gamma,X$ and $M$ follow the canonical notation for the Brillouin
zone of the PM phase.  Figure 11: (a)(b) HHG spectra of the Mott insulator
evaluated with DMFT for the indicated temperatures. Here, we consider the
Hubbard model on the two-dimensional square lattice and set $U=6$ and $t_{\rm
hop}=0.5$. The pulse parameters are (a)
$\Omega=1,E_{0}=0.8,t_{0}=30,\sigma=7.5$ and $\theta=0$, and (b)
$\Omega=1,E_{0}=0.8,t_{0}=30,\sigma=7.5$ and $\theta=\pi/4$. (c-f) The
corresponding polarization dependence of the HHG peak intensity. In all cases,
we evaluate the HHG intensity along the direction of the applied electric
field.

In Fig. 9, to support the explanations about the effect of thermal
fluctuations in the main text, we show schematic pictures of the dynamics of a
doublon-holon pair in different spin backgrounds. Panels i)-iii) show cases
with different spin configurations activated by thermal fluctuations. The
difference between case i) and case ii) is indicated by the blue dashed
rectangle, while that between case i) and case iii) is indicated by the green
dashed rectangle. The zigzag lines in each panel show the sites where the
mismatch in the spin configuration occurs as the doublon moves around. The
curved arrows with different colors indicate the difference in the kinematics
of the doublon due to the spin-charge coupling.

Finally, we comment on the relevance of the physics of the single-band Hubbard
model for Ca2RuO4. A more accurate modeling of Ca2RuO4 should involve a three-
band Hubbard model. In multi-orbital strongly correlated systems, spin as well
as orbital degrees of freedom become important. Still, we believe that similar
physics as discussed for the single-band Hubbard model is relevant in this
case. Firstly, the spin degrees of freedom should play a similar role also in
the three band case, since the dynamics of excited local multiplets disturbs
the spin background. Secondly, in multi-orbital systems, the orbital-charge
coupling can play a similar role as the spin-charge coupling, as has been
discussed for the charge relaxation in a previous study.Strand  _et al._
(2017) More specifically, the dynamics of the excited multiples disturbs the
orbital configurations at the cost of the Hund energy. We expect that, in
Ca2RuO4, both the spin-charge coupling and the orbital-charge coupling should
contribute to the peculiar behavior of HHG, while the physical mechanism is
essentially the one discussed for the single-band Hubbard model.

##  IV Supplementary results for the square lattice

In this section, we present supplementary DMFT results for the Hubbard model
on the two-dimensional square lattice at half filling. We show that the
equilibrium features and the characteristic temperature dependence of the HHG
are essentially the same as in the case of the Bethe lattice, although a
systematic analysis is difficult since the DMFT calculation becomes
numerically more demanding. We assume that the system only has nearest
neighbor hopping, $t_{\rm hop}$. We set $t_{\rm hop}=0.5$ so that the
bandwidth of the free system becomes $4$ as in the case of the Bethe lattice.
The lattice constant $a$ is set unity. In the following the system size is
$(L_{x},L_{y})=(16,16)$, where $L_{x}$ ($L_{y}$) is the number of sites along
the $x$ ($y$) axis, and we use periodic boundary conditions.

In Fig. 10, we show the local single-particle spectral function $A_{\rm
loc}(\omega)$ and the momentum-resolved single-particle spectral function
$A_{\bf k}(\omega)$ in equilibrium. The latter is defined as $A_{\bf
k}(\omega)=-\frac{1}{\pi}{\rm Im}G^{R}_{\bf k,\sigma}(\omega)$, where
$G^{R}_{\bf k,\sigma}(\omega)$ is the Fourier component of the retarded part
of the Green’s function $G_{\bf k,\sigma}(t)$. Here, $G_{\bf
k,\sigma}(t)=-i\langle T_{\mathcal{C}}\hat{c}_{\bf
k\sigma}(t)\hat{c}^{\dagger}_{\bf k\sigma}(0)\rangle$ and
$\hat{c}^{\dagger}_{\bf k\sigma}=\frac{1}{\sqrt{N}}\sum_{l}e^{i{\bf
k}\cdot{\bf r}_{l}}\hat{c}_{l\sigma}^{\dagger}$. As in the case of the Bethe
lattice, the local spectral function $A_{\rm loc}(\omega)$ is featureless
above the transition temperature $T_{c}$. Below $T_{c}$ the Mott gap is
enhanced and peaks corresponding to spin polarons emerge. In the PM phase, the
dispersions of the upper and lower Hubbard bands are parallel to each other,
similar to the prediction of the Hubbard I approximation, see Fig. 10(a). On
the other hand, in the AF phase, the momentum-resolved spectral function
consists of many almost-flat spin-polaron bands, see Fig. 10(b). Also, at each
momentum, the spectral function covers a wider frequency range than in the PM
phase, and as a whole the dispersion relation is less well-defined than in the
PM phase, i.e. the spectra look more incoherent. Based on the experience from
HHG in semiconductors, these observations lead to the natural expectation that
the HHG intensity should decrease with decreasing temperature.

Figure 12: (a)(b) Subcycle analysis of the HHG signal, $I_{\rm
HHG}(\omega,t_{p})$, for the field along the $(1,1)$ direction
($\theta=\frac{\pi}{4}$). Panel (a) is for $T=0.2$ (PM phase) and (b) is for
$T=0.05$ (AF phase). The pulse parameters are $\Omega=0.5$,
$E_{0}=0.9\sqrt{2}$, $t_{0}=60$, $\sigma=15$. The vertical dashed lines
indicate the times when $E(t)=0$.

Next, we discuss the HHG spectrum of this system. We apply linearly polarized
light along the ${\bf e}_{\theta}=[\cos(\theta),\sin(\theta)]$ direction:

$\displaystyle A_{x}(t)$ $\displaystyle=\frac{E_{0}}{\Omega}\cos(\theta)F_{\rm
G}(t,t_{0},\sigma)\sin(\Omega(t-t_{0})),$ (7) $\displaystyle A_{y}(t)$
$\displaystyle=\frac{E_{0}}{\Omega}\sin(\theta)F_{\rm
G}(t,t_{0},\sigma)\sin(\Omega(t-t_{0})).$ (8)

We extract the HHG signal along the ${\bf e}_{\theta}$ direction from the
Fourier transformation of $J_{\theta}(t)$, where $J_{\theta}(t)={\bf
e}_{\theta}\cdot{\bf J}(t)$ and ${\bf J}(t)=[J_{x}(t),J_{y}(t)]$. The
resulting HHG spectra are shown in Figs. 11(a)(b) for $\Omega=1$ and $U=6$.
These data show that the temperature dependence of the HHG intensity is
qualitatively the same as that obtained for the Bethe lattice, and the
enhancement of the HHG signal with decreasing temperature is independent of
the polarization. In Figs. 11(c-f), we show the polarization dependence of the
intensity of the HHG peaks for a given order. In the AF phase, the intensity
becomes largest when the field is polarized along the bond direction, which is
consistent with the experiment on Ca2RuO4 Uchida  _et al._ (2022). However,
interestingly, in the PM phase, the polarization dependence changes
qualitatively for the higher harmonic components ($n\geq 5$).

For lower frequencies, like $\Omega=0.5$ used in the main text, the DMFT self-
consistency becomes too expensive for the full simulation of the time
evolution. Still, one can simulate halfway, and analyze the subcycle features
around the peak of the pulse, see Fig. 12. Again, we find qualitatively the
same behavior as in the Bethe lattice, which supports the generality of the
physics discussed in the main tex.

##  V Supplemental results for other models

In the main text, we revealed the important effect of the correlations between
different degrees of freedom on HHG in strongly correlated systems, focusing
of the single-band Hubbard model, where the spin-charge coupling is strong. We
showed that this simple model can explain many peculiar HHG features reported
in an experimental study of Ca2RuO4.Uchida  _et al._ (2022) The underlying
physics should be applicable to a wide range of Mott insulators. In
particular, orbital-charge coupling in addition to spin-charge coupling can
yield similar phenomena in multi-orbital Mott insulators, since a previous
study Strand  _et al._ (2017) showed that orbital-charge coupling can have a
similar effect on nonequilibrium charge carriers as spin-charge coupling. More
specifically, in multi-orbital systems with Hund coupling, the dynamics of
local multiplets disturbs the orbital configurations in a similar way as the
doublon dynamics disturbs the spin background in the single-band Hubbard
model. Here, we briefly show some results for two other models to support the
above statements.

Figure 13:  (a) $I_{\text{HHG}}$ for different $B^{z}_{\text{stagg}}$ for the
charge transfer insulator described by the one-dimensional $d$-$p$ model.
(b)(c) The corresponding subcycle analysis for $B^{z}_{\text{stagg}}=0.01$ (b)
and $B^{z}_{\text{stagg}}=0.3$ (c). A Gaussian window with
$\sigma^{\prime}=0.9$ is used. In all panels, we set
$v_{dp}=3,\epsilon_{d}-\epsilon_{p}=-10,U_{dd}=20,U_{pp}=0,V_{dp}=0$ and
consider a three-quarter filled system. The excitation parameters are
$\Omega=0.5,E_{0}=0.8$, $t_{0}=60$ and $\sigma=15$.

###  V.1 $d$-$p$ model for charge transfer insulators

In this section we analyze the $d$-$p$ model in one dimension using iTEBD. The
Hamiltonian consists of the two parts

$\displaystyle\hat{H}_{dp}=\hat{H}_{dp,\text{kin}}+\hat{H}_{dp,\text{int}}.$
(9)

If we assume that the $d$ orbital corresponds to $d_{x^{2}-y^{2}}$ and the $p$
orbital to $p_{x}$, the explicit for of these terms is

$\displaystyle\hat{H}_{dp,\text{kin}}=$
$\displaystyle-\sum_{i,\sigma}[v_{dp}(t)\hat{p}^{\dagger}_{i,\sigma}\hat{d}_{i,\sigma}+h.c.]$
$\displaystyle+\sum_{i,\sigma}[v_{dp}(t)\hat{d}^{\dagger}_{i,\sigma}\hat{p}_{i-1,\sigma}+h.c.]$
(10)
$\displaystyle+\epsilon_{d}\sum_{i}n^{d}_{i}+\epsilon_{p}\sum_{i}n^{p}_{i}+B^{z}_{\text{stagg}}\sum_{i}(-1)^{i}\hat{S}^{d}_{z,i},$
$\displaystyle\hat{H}_{dp,\text{int}}$
$\displaystyle=U_{dd}\sum_{i}n^{d}_{i\uparrow}n^{d}_{i\downarrow}+U_{pp}\sum_{i}n^{p}_{i\uparrow}n^{p}_{i\downarrow}$
(11) $\displaystyle\;\;\;\;+V_{dp}\sum_{i}n^{d}_{i}(n^{p}_{i}+n^{p}_{i-1}).$

Here, $\hat{d}^{\dagger}$ represents the creation operator of the $d$-orbital
electron, while $\hat{p}^{\dagger}$ represents the creation operator of the
$p$-orbital electron.
$\hat{n}_{i\sigma}^{d}=\hat{d}^{\dagger}_{i\sigma}\hat{d}_{i\sigma}$,
$\hat{n}_{i\sigma}^{p}=\hat{p}^{\dagger}_{i\sigma}\hat{p}_{i\sigma}$,
$\hat{n}_{i}^{d}=\sum_{\sigma}\hat{n}_{i\sigma}^{d}$,
$\hat{n}_{i}^{p}=\sum_{\sigma}\hat{n}_{i\sigma}^{p}$ and
$\hat{S}^{d}_{z,i}=\frac{1}{2}[\hat{d}^{\dagger}_{i\uparrow}\hat{d}_{i\uparrow}-\hat{d}^{\dagger}_{i\downarrow}\hat{d}_{i\downarrow}]$.
$v_{dp}$ is the hopping parameter between the neighboring $d$ and $p$
orbitals, $\epsilon_{a}$ indicates the energy level of the orbital $a$, and
$B^{z}_{\text{stagg}}$ is the staggered magnetic field applied to the $d$
orbitals. $U_{dd}$ and $U_{pp}$ are the on-site interactions for the $d$ and
$p$ orbitals, respectively, and $V_{dp}$ is the nearest neighbor interaction.
We assume that the $d$ and $p$ orbitals are equally distanced, and set the
lattice constant and the charge to unity. The effect of the electric field is
considered through the Peierls substitution as in the Hubbard model, i.e.,
$v_{dp}(t)=v_{dp}e^{iA(t)/2}$. The current is defined as
$\hat{J}(t)=-\frac{\delta\hat{H}_{dp}}{\delta A(t)}$. We focus on the three-
quarter filled system, which corresponds to a Mott or charge transfer (CT)
insulator, and simulate the time evolution by iTEBD. The HHG intensity is
evaluated by the field induced current as in the case of the Hubbard model. We
note again that the staggered magnetic field mimics the effect of the spin-
charge coupling expected in higher dimensions.

In Fig. 13, we show the results for a parameter set corresponding to a CT
insulator. We choose the parameters such that the bandwidth of the upper
Hubbard band and the Zhang-Rice singlet band, as well as the band gap, become
similar to those in the Hubbard model studied in the main text. One can
identify essentially the same HHG features as in the results for the single-
band Hubbard model discussed in the main text. Namely, when
$B^{z}_{\text{stagg}}$ is small, the HHG spectrum does not show clear HHG
peaks, while with larger $B^{z}_{\text{stagg}}$ the HHG peaks become clearer,
see Fig. 13(a). Also the subcycle analysis shows that the coherence time of
the pairs of excited local multiplets is efficiently reduced by
$B^{z}_{\text{stagg}}$, see Fig. 13(b)(c). These results suggest that also in
CT insulators, the spin-charge coupling plays a crucial role in the HHG
mechanism and that the physics discussed for the single-band Hubbard model is
relevant.

Figure 14: Results for paramagnetic Mott insulators described by the two band
Hubbard model. (a) Equilibrium local spectral functions
$A_{\text{loc}}(\omega)$ computed with DMFT for various $T$. (b) HHG spectra
of the Mott insulator computed with DMFT. In all panels, we set $U=6,J=0.5$
and consider half filling. The excitation parameters are
$\Omega=0.5,E_{0}=0.8$, $t_{0}=75$ and $\sigma=15$.

###  V.2 Two-band Hubbard model

As a second example, we consider a two-orbital Hubbard model with Hamiltonian

$\displaystyle\hat{H}_{\text{2bH}}=\sum_{\langle
ij\rangle}\sum_{\alpha=1,2}\sum_{\sigma}v_{ij}^{\alpha}(t)\hat{c}^{\dagger}_{i,\alpha\sigma}\hat{c}_{j,\alpha\sigma}$
$\displaystyle+\sum_{i}\sum_{\alpha=1,2}\Bigl{[}U\hat{n}_{i,\alpha\uparrow}\hat{n}_{i,\alpha\downarrow}-\mu(\hat{n}_{i,\alpha\uparrow}+n_{i,\alpha\downarrow})\Bigl{]}$
(12)
$\displaystyle+\sum_{i}\sum_{\sigma}\Bigl{[}(U-2J)\hat{n}_{i,\alpha\sigma}\hat{n}_{i,\alpha\bar{\sigma}}+(U-3J)\hat{n}_{i,\alpha\sigma}\hat{n}_{i,\alpha\sigma}\Bigl{]}.$

Here $\hat{c}^{\dagger}_{\alpha\sigma}$ is the creation operator of the
electron in orbital $\alpha$ with spin $\sigma$ and
$\hat{n}_{i,\alpha\sigma}=\hat{c}^{\dagger}_{i,\alpha\sigma}\hat{c}_{i,\alpha\sigma}$.
$v_{ij}^{\alpha}$ is the hopping parameter between site $i$ and $j$ for
orbital $\alpha$, $\mu$ the chemical potential, $U$ the intra-orbital
interaction and $J$ the Hund coupling. For simplicity, we neglect the spin-
flip and pair-hopping terms. In the following, we put $v_{ij}^{1}=v_{ij}^{2}$,
and consider the Bethe lattice as in the Hubbard model case, following Refs.
Werner  _et al._ , 2018; Lysne  _et al._ , 2020b. We set the quarter of the
free electron bandwidth to unity. Focusing on half filling, we simulate the
time evolution of the system under an electric field pulse with DMFT.

Figure 15: Schematic pictures of the effects of the staggered magnetic field
$B^{z}_{\rm stagg}$ applied to the one-dimensional Hubbard chain, which can be
directly compared to Fig. 5 (b). i) In equilibrium, each site is occupied by
one electron, and the spin is aligned anti-parallel to the staggered field.
ii), iii) As the doublon moves, it produces a mismatch between the spin
configuration and the staggered magnetic field at the cost of the Zeeman
energy, about $B^{z}_{\rm stagg}$ (zigzag lines). Figure 16:  Momentum-
integrated spectrum $A_{\rm loc}(\omega)$ and momentum-resolved spectrum
$A_{k}(\omega)$ for the one-dimensional Hubbard model with staggered magnetic
field at half filling. Panel (a) is for $U=8$, $B_{\rm stagg}^{z}=0.001$, (b)
is for $U=8$, $B_{\rm stagg}^{z}=0.2$ and (c) is for $U=8$, $B_{\rm
stagg}^{z}=0.3$. Here, we use the iTEBD method.

In Fig. 14, we show the local spectrum averaged over spins and orbitals
($A_{\text{loc}}(\omega)$) and the HHG spectra for different temperatures.
Here, we simulate the system in the paramagnetic (PM) phase. As we decrease
the temperature, the band gap increases and peak structures develop within the
bands [see Fig. 14(a)]. This is in contrast to the single-band case, where the
local spectrum in the PM phase is featureless. The difference originates from
the Hund coupling and the resultant orbital-charge coupling. When a test
charge is injected into the half-filled system, a triplon is created and its
dynamics disturbs the orbital background at the cost of the Hund energy. This
results in a strong orbital-charge coupling, which plays a similar role as the
spin-charge coupling in the single-band Hubbard model. Figure 14(b) shows that
the HHG intensity increases as the temperature decreases (the gap increases),
as in the single-band Hubbard model. These results suggest that indeed the
orbital-charge coupling in multi-orbital systems has a similar effect on HHG
as the spin-charge coupling in the single-band Hubbard model, and thus the
physics discussed in the paper should also be applicable for this case. We
note that, if we go into the antiferromagnetic phase in the two-orbital model,
the spin-charge coupling also starts to play a role and a further increase of
the HHG signal is expected.

##  VI Supplementary data for the one-dimensional system

In this section, we present supplementary results obtained with iTEBD for the
one-dimensional Hubbard model with staggered magnetic field. First, we show a
schematic picture to illustrate the analogy between this model and the Hubbard
model in higher dimensions, which is discussed in the main text, see Fig. 15.
i) In equilibrium, each site is occupied by one electron, and the spin is
aligned anti-parallel to the staggered field. ii), iii) As the doublon moves,
it produces a mismatch between the spin configuration and the staggered
magnetic field at the cost of the Zeeman energy, about $B^{z}_{\rm stagg}$
(zigzag lines).

In Fig. 16, we show how the single-particle spectra change with the magnetic
field $B_{\rm stagg}^{z}$. One can see that, with increasing field strength,
there emerges a separated band at the bottom (top) of the upper (lower)
Hubbard band in the momentum-resolved spectrum $A_{k}(\omega)$. Furthermore,
replicas of these separated bands can be identified. In the momentum-
integrated spectrum $A_{\rm loc}(\omega)$, these sub-bands correspond to the
peaks that emerge with increasing field strength. These features are
qualitatively the same as those of the DMFT spectra, and originate from the
effective spin-charge coupling and the resulting spin-polarons that are
induced by the staggered magnetic field, as mentioned in the main text.

## References

  * Ferray  _et al._ (1988) M. Ferray, A. L’Huillier, X. F. Li, L. A. Lompre, G. Mainfray, and C. Manus, Journal of Physics B: Atomic, Molecular and Optical Physics 21, L31 (1988).
  * Krausz and Ivanov (2009) F. Krausz and M. Ivanov, Rev. Mod. Phys. 81, 163 (2009).
  * Ghimire  _et al._ (2011) S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and D. A. Reis, Nat. Phys. 7, 138 (2011).
  * Schubert  _et al._ (2014) O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, Nat. Photon. 8, 119 (2014).
  * Luu  _et al._ (2015) T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T. Hassan, and E. Goulielmakis, Nature (London) 521, 498 (2015).
  * Vampa  _et al._ (2015a) G. Vampa, T. J. Hammond, N. Thire, B. E. Schmidt, F. Legare, C. R. McDonald, T. Brabec, and P. B. Corkum, Nature (London) 522, 462 (2015a).
  * Langer  _et al._ (2016) F. Langer, M. Hohenleutner, C. P. Schmid, C. Pöllmann, P. Nagler, T. Korn, C. Schüller, M. Sherwin, U. Huttner, J. Steiner, S. Koch, M. Kira, and R. Huber, Nature (London) 533, 225 (2016).
  * Hohenleutner  _et al._ (2015) M. Hohenleutner, F. Langer, O. Schubert, M. Knorr, U. Huttner, S. Koch, M. Kira, and R. Huber, Nature (London) 523, 572 (2015).
  * Ndabashimiye  _et al._ (2016) G. Ndabashimiye, S. Ghimire, M. Wu, D. A. Browne, K. J. Schafer, M. B. Gaarde, and D. A. Reis, Nature (London) 534, 520 (2016).
  * Liu  _et al._ (2017) H. Liu, Y. Li, Y. S. You, S. Ghimire, T. F. Heinz, and D. A. Reis, Nat. Phys. 13, 262 (2017).
  * You  _et al._ (2017) Y. S. You, D. A. Reis, and S. Ghimire, Nat. Phys. 13, 345 (2017).
  * Yoshikawa  _et al._ (2017) N. Yoshikawa, T. Tamaya, and K. Tanaka, Science 356, 736 (2017).
  * Hafez  _et al._ (2018) H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green, N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Teichert, Z. Wang, K.-J. Tielrooij, Z. Liu, Z. Chen, A. Narita, K. Mullen, M. Bonn, M. Gensch, and D. Turchinovich, Nature (London) 561, 507 (2018).
  * Kaneshima  _et al._ (2018) K. Kaneshima, Y. Shinohara, K. Takeuchi, N. Ishii, K. Imasaka, T. Kaji, S. Ashihara, K. L. Ishikawa, and J. Itatani, Phys. Rev. Lett. 120, 243903 (2018).
  * Yoshikawa  _et al._ (2019) N. Yoshikawa, K. Nagai, K. Uchida, Y. Takaguchi, S. Sasaki, Y. Miyata, and K. Tanaka, Nat. Comm. 10, 1 (2019).
  * Cheng  _et al._ (2020) B. Cheng, N. Kanda, T. N. Ikeda, T. Matsuda, P. Xia, T. Schumann, S. Stemmer, J. Itatani, N. P. Armitage, and R. Matsunaga, Phys. Rev. Lett. 124, 117402 (2020).
  * Schmid  _et al._ (2021) C. P. Schmid, L. Weigl, P. Grössing, V. Junk, C. Gorini, S. Schlauderer, S. Ito, M. Meierhofer, N. Hofmann, D. Afanasiev, J. Crewse, K. A. Kokh, O. E. Tereshchenko, J. Güdde, F. Evers, J. Wilhelm, K. Richter, U. Höfer, and R. Huber, Nature 593, 385 (2021).
  * Golde  _et al._ (2008) D. Golde, T. Meier, and S. W. Koch, Phys. Rev. B 77, 075330 (2008).
  * Vampa  _et al._ (2014) G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B. Corkum, and T. Brabec, Phys. Rev. Lett. 113, 073901 (2014).
  * Vampa  _et al._ (2015b) G. Vampa, C. R. McDonald, G. Orlando, P. B. Corkum, and T. Brabec, Phys. Rev. B 91, 064302 (2015b).
  * Wu  _et al._ (2015) M. Wu, S. Ghimire, D. A. Reis, K. J. Schafer, and M. B. Gaarde, Phys. Rev. A 91, 043839 (2015).
  * Otobe (2016) T. Otobe, Phys. Rev. B 94, 235152 (2016).
  * Ikemachi  _et al._ (2017) T. Ikemachi, Y. Shinohara, T. Sato, J. Yumoto, M. Kuwata-Gonokami, and K. L. Ishikawa, Phys. Rev. A 95, 043416 (2017).
  * Tancogne-Dejean  _et al._ (2017) N. Tancogne-Dejean, O. D. Mücke, F. X. Kärtner, and A. Rubio, Nat. Comm. 8, 745 (2017).
  * Luu and Wörner (2016) T. T. Luu and H. J. Wörner, Phys. Rev. B 94, 115164 (2016).
  * Hansen  _et al._ (2017) K. K. Hansen, T. Deffge, and D. Bauer, Phys. Rev. A 96, 053418 (2017).
  * Osika  _et al._ (2017) E. N. Osika, A. Chacón, L. Ortmann, N. Suárez, J. A. Pérez-Hernández, B. Szafran, M. F. Ciappina, F. Sols, A. S. Landsman, and M. Lewenstein, Phys. Rev. X 7, 021017 (2017).
  * Ikeda  _et al._ (2018) T. N. Ikeda, K. Chinzei, and H. Tsunetsugu, Phys. Rev. A 98, 063426 (2018).
  * Tamaya  _et al._ (2016) T. Tamaya, A. Ishikawa, T. Ogawa, and K. Tanaka, Phys. Rev. Lett. 116, 016601 (2016).
  * Floss  _et al._ (2018) I. Floss, C. Lemell, G. Wachter, V. Smejkal, S. A. Sato, X.-M. Tong, K. Yabana, and J. Burgdörfer, Phys. Rev. A 97, 011401 (2018).
  * Lysne  _et al._ (2020a) M. Lysne, Y. Murakami, M. Schüler, and P. Werner, Phys. Rev. B 102, 081121 (2020a).
  * Chacón  _et al._ (2020) A. Chacón, D. Kim, W. Zhu, S. P. Kelly, A. Dauphin, E. Pisanty, A. S. Maxwell, A. Picón, M. F. Ciappina, D. E. Kim, C. Ticknor, A. Saxena, and M. Lewenstein, Phys. Rev. B 102, 134115 (2020).
  * Wilhelm  _et al._ (2021) J. Wilhelm, P. Grössing, A. Seith, J. Crewse, M. Nitsch, L. Weigl, C. Schmid, and F. Evers, Phys. Rev. B 103, 125419 (2021).
  * Taya  _et al._ (2021) H. Taya, M. Hongo, and T. N. Ikeda, Phys. Rev. B 104, L140305 (2021).
  * Vampa  _et al._ (2015c) G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, D. D. Klug, and P. B. Corkum, Phys. Rev. Lett. 115, 193603 (2015c).
  * Li  _et al._ (2020) L. Li, P. Lan, L. He, W. Cao, Q. Zhang, and P. Lu, Phys. Rev. Lett. 124, 157403 (2020).
  * Luu and Wörner (2018) T. T. Luu and H. J. Wörner, Nat. Comm. 9, 916 (2018).
  * Uchida  _et al._ (2021) K. Uchida, V. Pareek, K. Nagai, K. M. Dani, and K. Tanaka, Phys. Rev. B 103, L161406 (2021).
  * Kemper  _et al._ (2013) A. F. Kemper, B. Moritz, J. K. Freericks, and T. P. Devereaux, New J. Phys. 15, 023003 (2013).
  * Kruchinin  _et al._ (2018) S. Y. Kruchinin, F. Krausz, and V. S. Yakovlev, Rev. Mod. Phys. 90, 021002 (2018).
  * Ghimire and Reis (2019) S. Ghimire and D. A. Reis, Nat. Phys. 15, 10 (2019).
  * Du and Ma (2022) T.-Y. Du and C. Ma, Phys. Rev. A 105, 053125 (2022).
  * Uchida  _et al._ (2022) K. Uchida, G. Mattoni, S. Yonezawa, F. Nakamura, Y. Maeno, and K. Tanaka, Phys. Rev. Lett. 128, 127401 (2022).
  * Silva  _et al._ (2018) R. E. F. Silva, I. V. Blinov, A. N. Rubtsov, O. Smirnova, and M. Ivanov, Nat. Photon. 12, 266 (2018).
  * Murakami  _et al._ (2018) Y. Murakami, M. Eckstein, and P. Werner, Phys. Rev. Lett. 121, 057405 (2018).
  * Murakami and Werner (2018) Y. Murakami and P. Werner, Phys. Rev. B 98, 075102 (2018).
  * Lysne  _et al._ (2020b) M. Lysne, Y. Murakami, and P. Werner, Phys. Rev. B 101, 195139 (2020b).
  * Tancogne-Dejean  _et al._ (2018) N. Tancogne-Dejean, M. A. Sentef, and A. Rubio, Phys. Rev. Lett. 121, 097402 (2018).
  * Imai  _et al._ (2020) S. Imai, A. Ono, and S. Ishihara, Phys. Rev. Lett. 124, 157404 (2020).
  * Chinzei and Ikeda (2020) K. Chinzei and T. N. Ikeda, Phys. Rev. Research 2, 013033 (2020).
  * Orthodoxou  _et al._ (2021) C. Orthodoxou, A. Zaïr, and G. H. Booth, npj Quantum Materials 6, 76 (2021).
  * Murakami  _et al._ (2021) Y. Murakami, S. Takayoshi, A. Koga, and P. Werner, Phys. Rev. B 103, 035110 (2021).
  * Shao  _et al._ (2022) C. Shao, H. Lu, X. Zhang, C. Yu, T. Tohyama, and R. Lu, Phys. Rev. Lett. 128, 047401 (2022).
  * Hansen  _et al._ (2022) T. Hansen, S. V. B. Jensen, and L. B. Madsen, arXiv:2201.06879 (2022).
  * Masur  _et al._ (2022) J. Masur, D. I. Bondar, and G. McCaul, arXiv:2202.06895 (2022).
  * Udono  _et al._ (2022) M. Udono, K. Sugimoto, T. Kaneko, and Y. Ohta, Phys. Rev. B 105, L241108 (2022).
  * Grånäs  _et al._ (2020) O. Grånäs, I. Vaskivskyi, P. Thunström, S. Ghimire, R. Knut, J. Söderström, L. Kjellsson, D. Turenne, R. Y. Engel, M. Beye, J. Lu, A. H. Reid, W. Schlotter, G. Coslovich, M. Hoffmann, G. Kolesov, C. Schüßler-Langeheine, A. Styervoyedov, N. Tancogne-Dejean, M. A. Sentef, D. A. Reis, A. Rubio, S. S. P. Parkin, O. Karis, J. Nordgren, J. E. Rubensson, O. Eriksson, and H. A. Dürr, arXiv:2008.11115 (2020).
  * Bionta  _et al._ (2021) M. R. Bionta, E. Haddad, A. Leblanc, V. Gruson, P. Lassonde, H. Ibrahim, J. Chaillou, N. Émond, M. R. Otto, A. Jiménez-Galán, R. E. F. Silva, M. Ivanov, B. J. Siwick, M. Chaker, and F. m. c. Légaré, Phys. Rev. Research 3, 023250 (2021).
  * Imada  _et al._ (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * Dagotto (1994) E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
  * Note (1) Supplemental Material [url], which includes Ref. Werner  _et al._ (2018).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Eckstein and Werner (2010) M. Eckstein and P. Werner, Phys. Rev. B 82, 115115 (2010).
  * Aoki  _et al._ (2014) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
  * Sandholzer  _et al._ (2019) K. Sandholzer, Y. Murakami, F. Görg, J. Minguzzi, M. Messer, R. Desbuquois, M. Eckstein, P. Werner, and T. Esslinger, Phys. Rev. Lett. 123, 193602 (2019).
  * Schüler  _et al._ (2020) M. Schüler, D. Golež, Y. Murakami, N. Bittner, A. Herrmann, H. U. Strand, P. Werner, and M. Eckstein, Computer Physics Communications 257, 107484 (2020).
  * Werner  _et al._ (2017) P. Werner, H. U. R. Strand, S. Hoshino, and M. Eckstein, Phys. Rev. B 95, 195405 (2017).
  * Martinez and Horsch (1991) G. Martinez and P. Horsch, Phys. Rev. B 44, 317 (1991).
  * Sangiovanni  _et al._ (2006) G. Sangiovanni, A. Toschi, E. Koch, K. Held, M. Capone, C. Castellani, O. Gunnarsson, S.-K. Mo, J. W. Allen, H.-D. Kim, A. Sekiyama, A. Yamasaki, S. Suga, and P. Metcalf, Phys. Rev. B 73, 205121 (2006).
  * Lenarčič and Prelovšek (2013) Z. Lenarčič and P. Prelovšek, Phys. Rev. Lett. 111, 016401 (2013).
  * Golež  _et al._ (2014) D. Golež, J. Bonča, M. Mierzejewski, and L. Vidmar, Phys. Rev. B 89, 165118 (2014).
  * Eckstein and Werner (2016) M. Eckstein and P. Werner, Scientific Reports 6, 21235 (2016).
  * Vidal (2003) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
  * Kilen  _et al._ (2020) I. Kilen, M. Kolesik, J. Hader, J. V. Moloney, U. Huttner, M. K. Hagen, and S. W. Koch, Phys. Rev. Lett. 125, 083901 (2020).
  * Becker  _et al._ (1988) P. C. Becker, H. L. Fragnito, C. H. B. Cruz, R. L. Fork, J. E. Cunningham, J. E. Henry, and C. V. Shank, Phys. Rev. Lett. 61, 1647 (1988).
  * Nagai  _et al._ (2021) K. Nagai, K. Uchida, S. Kusaba, T. Endo, Y. Miyata, and K. Tanaka, “Effect of incoherent electron-hole pairs on high harmonic generation in atomically thin semiconductors,” (2021), arXiv:2112.12951 [physics.optics] .
  * Strand  _et al._ (2017) H. U. R. Strand, D. Golež, M. Eckstein, and P. Werner, Phys. Rev. B 96, 165104 (2017).
  * Werner  _et al._ (2018) P. Werner, H. U. R. Strand, S. Hoshino, Y. Murakami, and M. Eckstein, Phys. Rev. B 97, 165119 (2018).