# Floquet engineering of binding in doped and photo-doped Mott insulators

Madhumita Sarkar  Jožef Stefan Institute, SI-1000, Ljubljana, Slovenia  Zala
Lenarčič  Jožef Stefan Institute, SI-1000, Ljubljana, Slovenia  Denis Golež
Jožef Stefan Institute, SI-1000, Ljubljana, Slovenia  Faculty of Mathematics
and Physics, University of Ljubljana, 1000 Ljubljana, Slovenia

###### Abstract

We investigate the emergence of bound states in chemically and photo-doped
Mott insulators, assisted by spin and $\eta$-pairing fluctuations within both
2-leg ladder and 2D systems. We demonstrate that the binding energies and
localization length in the chemically and photo-doped regimes are comparable.
To effectively describe the photo and chemically doped state on the same
footings, we employ the Schrieffer-Wolff transformation, resulting in a
generalized $t$-$J$ model. Furthermore, we show that manipulating the binding
is possible through external periodic driving, a technique known as Floquet
engineering, leading to significantly enhanced binding energies. We also
roughly estimate the lifetime of photo-doped states under periodic driving
conditions based on the Fermi golden rule. Lastly, we propose experimental
protocols for realizing Hubbard excitons in cold-atom experiments.

##  I Introduction

Doped Mott insulators are one of the most intensively studied solid-state
systems due to their significance for the high-temperature superconductors [1]
and a plead of intriguing exotic phases including pseudogap [2, 3, 4],stripes
[5, 6, 7], etc. One of the basic questions is how chemically doped charge
carriers interact with the antiferromagnetic background, and it was proposed
early on that carriers can form bound pairs whose binding energy originated
from the shared distortion of the spin background (string-based pairing) [8,
9, 10, 11, 12, 13, 14, 15]. For the binding of charge carriers, the
dimensionality of the system plays a particularly important role as in the
purely 1D systems spin-charge separation prevents pairing [16, 17, 18, 19],
making ladder [8, 20] and 2D systems [11, 15] minimal setups for studying of
bound pairs. Experimental progress in cold atoms simulators enabled a direct
observation of antiferromagnetic correlations [21], dressed spin polarons [22]
and pairing strings [23, 24]. Furthermore, combining cold atom experiments and
strong external fields enables a unique opportunity to manipulate microscopic
parameters, like hopping and superexchange amplitude [25, 26]. Recent examples
include enhanced binding of holes under a strong DC field in regimes where
superexchange becomes comparable to the hopping integrals [27, 28], a
situation not available in solid-state setups. These ideas stimulate exploring
whether one could employ periodic (Floquet) driving [29, 30, 31], whose
ability to manipulate microscopic parameters was already proved in cold-atoms
experiments [32], to stabilize bound pairs in doped Mott insulators even
further.

Photo doping of Mott insulators is an emerging direction to form exotic phases
of matter [33, 34]. While the photo-doped systems are metastable, their
lifetime can be exponentially long-lived for large gap Mott insulators [35,
36, 37, 38]. These states can exhibit novel non-thermal correlations, like the
formation of Hubbard excitons [39, 40, 41, 42, 43, 44] or $\eta$-pairing [45,
46, 47, 48, 49, 50, 51, 52], etc. In the high-dimensional systems, it was
shown that $\eta$-pairing fluctuations can lead to the condensation and form a
superfluid [53], while in the one-dimensional setup a spin-charge-$\eta$
separation takes place [54]. An important open question is to establish what
are the properties of photo-doped states beyond the two extreme case studies,
like the ladders and 2D situation, as these are the most relevant experimental
situations.

In this work, we compare the formation of a bound pair for chemically and
photo-doped systems in the ladder and 2D geometries. To describe the
chemically doped and metastable photo-doped systems on the same footing, we
resort to the canonical transformation, which perturbatively decoupled sectors
with different numbers of holons and doublons. We establish that the pairing
due to spin and $\eta$ fluctuations does not compete or actively cooperate but
rather leads to similar bindings between charge carriers in chemically and
photo-doped systems. We provide an approximate measure of the binding energy
for both chemically and photo-doped situations. Further, we show that periodic
driving with an electromagnetic field enables an efficient manipulation of
binding energies with its substantial increase for ladder and 2D systems and
more localized charge distribution in real space. For photo-doped systems, we
estimate the upper bound for the doublon-holon recombination rate in the
presence of Floquet driving and proposed driving regimes where the exciton
binding energies are larger than the lifetime. Finally, we comment on
preparing photo-doped states in cold atom experiments by either chirped
driving across the gap or adiabatic deformation of the lattice potential from
band to Mott insulator.

##  II Model and method

We consider the Fermi Hubbard model

$\displaystyle H$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
i,j\rangle,\sigma}(c^{\dagger}_{i\sigma}c_{j\sigma}+c^{\dagger}_{j\sigma}c_{i\sigma})+U\sum_{j}n_{j\uparrow}n_{j\downarrow},$
(1)

where sum $\langle i,j\rangle$ runs over nearest neighbour pairs of sites,
$\sigma$ is the spin index denoting either up or down spin, $t$ is the hopping
amplitude and $U$ the interaction strength. We will focus on the strong
interaction limit $U\gg t$, in which case the ground state at half-filling is
the Mott insulating state. We will consider two different geometries with $N$
sites: (i) quasi one-dimension two-leg ladder with periodic boundary
conditions along the elongated (chain) direction and (ii) two-dimensional
lattice with periodic boundary conditions along both directions. The aim of
our study is to understand and compare charge pairing and the role of spin and
$\eta$ fluctuations in photo and chemically doped correlated systems.

While chemically doped systems are stable, the photo-doped situation for the
large gap Mott insulator is only metastable with exponentially long lifetime
[36, 35, 38, 55]. To treat the two scenarios on the same footing, we will
perform a canonical Schrieffer-Wolf transformation which perturbatively
decouples sectors with different number of effectively charged holons and
doublons and obtain a generalized version of celebrated $t$-$J$ Hamiltonian
[56, 53, 47]. In the case of the photo-doped system, by neglecting higher-
order recombination terms, we approximate the metastable state with an
equilibrium state. Conveniently, such a description reduces the computational
complexity. We then estimate the Hubbard exciton lifetime by Fermi golden rule
in analogy with Ref. 38.

Figure 1: Schematic representation of terms entering the effective model
$H_{\text{eff}}$, see Eq. (II). Top: the spin exchanges process, which
exchanges the position of up and down spin; the $\eta$-exchange process, which
exchanges the position of a doublon and a holon. Bottom: The hopping process
of holon/doublon without changing the number of DH pairs.

We denote the onsite interaction term in Eq. (1) with
$H_{0}=U\sum_{j}n_{j\uparrow}n_{j\downarrow}$, and split the kinetic part into
terms that preserves the number of doublon-holon (DH) pairs $T_{0}$, that
increases the number of DH pairs by one $T_{1}$, and that decreases the number
of DH pairs by one $T_{-1}$. Following Ref. 56, the generator of the unitary
transformation $S$ is chosen such that it transforms out the processes that
change the number of DH pairs on the order $O(t)$ and retains them on the
higher $O(t^{2}/U)$ level of expansion. Hence, we choose $S=(T_{1}-T_{-1})/U$
such that $T_{1}$ and $T_{-1}$ are cancelled (in Eq. (II)),

$\displaystyle H^{\prime}$ $\displaystyle=$ $\displaystyle e^{iS}He^{-iS}$
$\displaystyle=$ $\displaystyle H_{0}+T_{0}+T_{1}+T_{-1}+[iS,H]+\dots$
$\displaystyle=$ $\displaystyle
H_{0}+T_{0}+\frac{1}{U}\big{(}[T_{1},T_{-1}]+[T_{0},T_{-1}]+[T_{1},T_{0}]\big{)}+\dots$

In the following, we distinguish the hopping amplitudes in x- and y-direction
and denote them as $t_{\parallel}$ and $t_{\perp}$ respectively. For 2D model,
$t=t_{\parallel}=t_{\perp}$, as can be easily imposed in the following
expressions. In the limit $U\gg t_{\parallel},t_{\perp}$, we retain terms of
the order $O(t_{\parallel}^{2}/U)$, $O(t_{\perp}^{2}/U$) and neglect all
higher orders, such that the effective Hamiltonian preserves the number of DH
pairs

$\displaystyle H_{\text{eff}}=$
$\displaystyle\sum\limits_{\alpha\in\\{\parallel,\perp\\}}\left[H_{\eta}^{\alpha}+H_{S}^{\alpha}+H_{H,hop}^{\alpha}+H_{D,hop}^{\alpha}\right]+H_{0}.$

The $\parallel$ terms act on the charge and spins along the x-direction (along
the chain) whereas the $\perp$ terms act in the y-direction (along the rung
for the two-leg model). When written out explicitly, the effective Hamiltonian
has the following terms[47]

$\displaystyle H_{S}^{\parallel(\perp)}$ $\displaystyle=$ $\displaystyle
J_{\parallel(\perp),S}\sum_{\langle
i,j\rangle,\sigma}\left(\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{j}-\frac{1}{4}\delta_{1,n_{i}n_{j}}\right),$
$\displaystyle H_{\eta}^{\parallel(\perp)}$ $\displaystyle=$ $\displaystyle-
J_{\parallel(\perp),\eta}\sum_{\langle
i,j\rangle,\sigma}\left(\boldsymbol{\eta}_{i}\cdot\boldsymbol{\eta}_{j}-\frac{1}{4}(1-\delta_{1,n_{i}})(1-\delta_{1,n_{j}})\right),$
$\displaystyle H_{D,hop}^{\parallel(\perp)}$ $\displaystyle=$ $\displaystyle-
t_{\parallel(\perp)}\sum_{\langle
i,j\rangle,\sigma}n_{i\bar{\sigma}}(c^{\dagger}_{i\sigma}c_{j\sigma}+h.c.)n_{j\bar{\sigma}},$
$\displaystyle H_{H,hop}^{\parallel(\perp)}$ $\displaystyle=$ $\displaystyle-
t_{\parallel(\perp)}\sum_{\langle
i,j\rangle,\sigma}\bar{n}_{i\bar{\sigma}}(c^{\dagger}_{i\sigma}c_{j\sigma}+h.c.)\bar{n}_{j\bar{\sigma}},$
(3)

where $J_{\parallel,S}=J_{\parallel,\eta}=\frac{4t_{\parallel}^{2}}{U}$ and
$J_{\perp,S}=J_{\perp,\eta}=\frac{4t_{\perp}^{2}}{U}$ are the spin/$\eta$
exchange coupling parameters along the chain and rung, respectively. When of
the same strength, we drop the subscript $S$ and $\eta$. The electron (hole)
density per spin is given by
$n_{i,\sigma}~{}(\bar{n}_{i,\sigma}=1-n_{i,\sigma})$, and $\bar{\sigma}$ is
the complement of $\sigma$, i.e., if $\sigma$ denotes up spin then
$\bar{\sigma}$ denotes down spin and vice versa. The $\eta$ operators in terms
of the fermion creation and annihilation operators can be expressed as
$\eta_{i}^{+}=\eta_{i}^{x}+i\eta_{i}^{y}=e^{i\boldsymbol{r}_{i}\cdot\boldsymbol{\pi}}c_{i,\uparrow}^{\dagger}c_{i,\downarrow}^{\dagger}$,
where $\boldsymbol{r}_{i}$ is the vector to the $i$th site and
$\boldsymbol{\pi}=\\{\pi,\pi,...\\}$ is a vector of the dimension of the
system with entries $\pi$. The $z$ component of the pseudospin is
$\eta_{i}^{z}=\frac{1}{2}(1-n_{i})$, $\eta_{i}^{-}=(\eta_{i}^{+})^{\dagger}$,
and the spin operators are defined as,
$\boldsymbol{S}_{i}=\frac{1}{2}c_{i\alpha}^{\dagger}\boldsymbol{\sigma}_{\alpha\beta}c_{i\beta}$,
where $\boldsymbol{\sigma}=\\{\sigma^{x},\sigma^{y},\sigma^{z}\\}$ are the
standard Pauli matrices.

Different terms in the Hamiltonian have the following physical meanings:
$H_{S}^{\parallel(\perp)}$ is the spin interaction term, familiar from the
standard t-J model; $H_{\eta}^{\parallel(\perp)}$ is the $\eta$-exchange term,
which exchanges the position of a doublon and holon and causes the interaction
between DH pair if they are on nearest neighbor sites;
$H_{H,hop}^{\parallel(\perp)}$ and $H_{D,hop}^{\parallel(\perp)}$ represent
hopping of holon and doublon conserving their number, respectively. These
processes are illustrated in Fig. 1. In the sectors with no doublon quasi-
particle, $H_{\text{eff}}$ reduces to the standard t-J model. $H_{0}$ is the
on-site Coulomb interaction and leads to unimportant energy shift between the
chemical and photo-doped system. Both chemically and photo-doped systems
reduce to the solution of the ground state problem in different sectors and we
employ the exact diagonalization using the Lanczos technique [57, 58].

We should note that in the full expansion there appear also three site terms
which are of two types: a) the holon/doublon correlated hopping terms of order
$J_{\parallel(\perp)}$ which conserve the number of holon and doublon pairs
and the analysis in Ref. 47 showed that their effect is small in the dilute
limit considered here, b) recombination terms of order $J_{\parallel(\perp)}$
which we will for a moment suppress to mimic the metastable state by an
equilibrium problem. Later on, we will estimate the lifetime of the metastable
phase (recombination time) by perturbative treatment of these 3-site
recombination terms.

##  III Results for two-leg ladder

Our aim is to study the properties of bound doublon-holon (Hubbard excitons)
or two holon states and compare the binding of charged carriers in photo and
chemically doped situation. An important goal is to understand the interplay
of spin and $\eta$ fluctuations on the doublon-holon binding, in systems
beyond one dimension, which is special due to the spin-charge and
spin-$\eta$-charge separation [47, 54].

We start our analysis with the simplest system which escapes the separation,
namely the two-leg ladder system, and later extend the study to the two-
dimensional situation.

###  III.1 Exciton binding energy

Figure 2: Binding energy versus the isotropic superexchange
$J=J_{\parallel,S/\eta}=J_{\perp,S/\eta}$ for (a) a holon-holon pair and for
(b) a doublon-holon pair at $L=12,16,20$ and $U=10$. Binding is present for a
large regime of $J$ and increases with $J$, since breaking of spin and $\eta$
fluctuations is more and more energetically penalized, leading to the
localization of charged pairs.

The tendency towards binding between charged particles can be quantified by
the binding energy [8, 59, 43, 38], which for the chemically doped state is
expresses as $E_{b}^{HH}=E_{0}^{HH}+E_{0}^{0}-2E_{0}^{H}$ and for the photo-
doped states by $E_{b}^{DH}=E_{0}^{DH}+E_{0}^{0}-E_{0}^{D}-E_{0}^{H}$, where
$E_{0}^{HH},E_{0}^{DH},E_{0}^{H},E_{0}^{D},E_{0}^{0},$ refer to the ground
state energy of $H_{\text{eff}}$ in the sector with two holons, one doublon-
holon pair, one holon, one doublon, and only singly occupied sites,
respectively. The criterion for the formation of a bound pair is a negative
binding energy, $E_{b}^{DH},E_{b}^{HH}<0$. In all of the results below, we
will measure the binding energy in the units of $t_{\parallel}$. We should
stress that the analysis below suffers from a particular finite size effect:
while $E_{0}^{0},E_{0}^{HH},E_{0}^{DH}$ are energies of ground states at zero
momentum, ground state for the sector with a single holon/doublon appear at a
finite momentum $k_{x}=\pi/2$, which might not be present for a given system
size $L$. In any case, we use the solution with lowest energy. To validate the
robustness of our conclusion despite finite size difficulties, we will further
compare the binding energy with real-space correlators to support the binding.

We start the analysis with the isotropic case
($J=J_{\parallel,S/\eta}=J_{\perp,S/\eta}$) and show that a finite binding is
present for a large parameter regime, $J\gtrsim 0.15$, and it increases as a
function of superexchange coupling, see Fig. 2. At large $J$, the magnitude of
binding between a DH pair and two holons become identical, while for small
values of $J$ the binding energies for the two kinds of exciton differ;
however the difference is getting smaller with increasing system size. The
increased binding in the limit of large superexchange $J$ comes from the fact
that the breaking of spin and $\eta$ fluctuations is more and more
energetically penalized, leading to localization of the charged pair.

Figure 3: Binding energy in the anisotropic system with varying perpendicular
superexchange $J_{\perp}$ and fixed parallel superexhange $J_{\perp}=0.4$ for
a doublon-holon pair (green) and for two holes (black lines). Initially,
binding energy grows with increasing $J_{\perp}=\frac{4t_{\perp}^{2}}{U}$, but
then the trend is reversed due to the competition of charge delocalization,
promoted by the increase in $t_{\perp}$, and localization aided by large
$J_{\perp}$. Parameters: $L=20$, $t_{\parallel}=1$, $U=10$.

Recent theoretical and experimental studies on binding of two holons showed
that by using external electic field [28, 27] one can create highly
anisotropic superexhange couplings, which can be favourable for binding. Here,
we would like to explore this perspective by tuning the anisotropy to
understand which parameter regime is favourable for both type of bindings.
Fig. 3 shows the binding energy as a function of $J_{\perp}$, keeping the
parallel exchange parameter fixed $J_{\parallel}=0.4$. Initially, binding
energy grows with increasing $J_{\perp}=\frac{4t_{\perp}^{2}}{U}$, but then
the trend is reversed due to the competition of charge delocalization,
promoted by the increase in $t_{\perp}$, and localization aided by large
$J_{\perp}$. This suggests that ideally, we would like to decouple $t_{\perp}$
and $J_{\perp}$ to increase only the latter. In Section III.3, we will
elucidate the methodology employed to accomplish this objective through
Floquet engineering. As a side note, we notice that the binding energy for the
two different charged pairs becomes equal for $J_{\perp}\gtrsim 0.25$.

Comparable binding energy for chemically and photo-doped systems means there
is no clear competition or cooperation for binding due to spin and $\eta$
fluctuations. We consider this as the first important message of this paper.
An useful starting point to understand the observation is to analyze the
ground state of a 2$\times$2 cluster for the anisotropic case,
$J_{\perp}>J_{\parallel}$, such that the DH or HH pairs are present in the
same rung to reduce the energy lost by the breaking of the spin bonds [8, 59].
The ground state manifold for two holons consists of two states,
$|\frac{0\uparrow}{0\downarrow}\rangle$ and
$|\frac{0\downarrow}{0\uparrow}\rangle$, where the top row represent the spin
configuration in the upper chain and bottom in the lower chain. The ground
state energy of this manifold is $\tilde{E}_{0}^{HH}=-J_{\perp,S}$. For the
case of one DH pair, we work in the manifold with four states,
$|\frac{0\uparrow}{\updownarrow\downarrow}\rangle$,
$\mid\frac{0\downarrow}{\updownarrow\uparrow}\rangle$,
$\mid\frac{\updownarrow\uparrow}{0\downarrow}\rangle$ and
$\mid\frac{\updownarrow\downarrow}{0\uparrow}\rangle$. The ground state energy
of this manifold is $\tilde{E}_{0}^{DH}=-(J_{\perp,S}+J_{\perp,\eta})/2$.
Notice that the binding energies $\tilde{E}_{0}^{DH}$ and $\tilde{E}_{0}^{HH}$
become identical only when we consider the $\eta$-exchange terms between the
doublon and holon. This outlines the importance of this term for the analysis
of Hubbard excitons [47, 53]. To obtain analytical estimations, we have
neglected the hopping terms, which crucially differ between holon and doublons
as well. While holon (as a quasiparticle) has a negative hopping sign, doublon
has a positive one since by moving a doublon we always hop over one fermion.
At least in 2D systems, this can crucially influence the symmetry of the bound
states, as pointed out in Refs. 43, 38, 60.

###  III.2 Charge and spin correlators

An alternative way to characterize the properties of bound states are real
space correlators. The charge correlators determine the relative spatial
distribution of charged carriers. For the DH and HH pair, they are defined as
$\chi_{c}^{DH}(\boldsymbol{\delta}{r}_{j})\equiv\frac{1}{N}\sum_{i}\langle
n_{i+j,\uparrow}n_{i+j,\downarrow}\bar{n}_{i,\uparrow}\bar{n}_{i,\downarrow}\rangle$
and
$\chi_{c}^{HH}(\boldsymbol{\delta}\textbf{r}_{j})\equiv\frac{1}{2N}\sum_{i}\langle\bar{n}_{i+j,\uparrow}\bar{n}_{i+j,\downarrow}\bar{n}_{i,\uparrow}\bar{n}_{i,\downarrow}\rangle$,
respectively. Furthermore, to capture the relative alignment of the up and
down spins we define the spin correlator
$\chi_{s}(\boldsymbol{\delta}\textbf{r}_{j})\equiv\frac{1}{N}\sum_{i}\langle
S_{i+r}^{z}S_{i}^{z}\rangle$.

Figure 4: (a) Inter-chain $\chi_{c}(\\{1,r\\})$ and (b) intra-chain charge
correlation $\chi_{c}(\\{0,r\\})$ reveal that charged particles prefer to sit
in different chains along the diagonal of the rung. Exponentially decaying
charge correlation is consistent with a bound nature of states considered. (c)
Inter-chain $\chi_{S}(\\{1,r\\})$ and (d) intra-chain spin correlation
$\chi_{S}(\\{0,r\\})$ for doublon-holon and holon-holon pairs reveal
antiferromagnetic local spin fluctuations which mediate the binding.
Parameters: $L=20$, $J_{\perp}=J_{\parallel}=0.4$.

We show the charge correlator (top) and spin correlator (bottom) between the
chains (left) and along the chain (right) in Fig. 4. For the isotropic case
$J_{\perp}=J_{\parallel}=0.4$ shown, not only the binding energies
$E_{b}^{HH}=E_{b}^{DH}$, c.f. Fig. 3, but also the correlators are
indistinguishable for the DH and HH pairs. Exponentially decaying charge
correlations for charges positioned in different legs, Fig. 4(a) are in
agreement with the state being bound. It is worth noting that the binding
along the same leg is considerably weaker, Fig. 4(b). The spin correlator,
Fig. 4(c,d), confirms anti-ferromagnetic local fluctuations, which mediate the
binding. While hopping terms tend to delocalise the charged pair, this causes
a disturbance in the antiferromagnetic spin pattern. The competition between
the two processes yields the existence of non-trivial bound states, with, for
example, the strongest occupation of charged particles along the diagonal and
not along the rung bond, as would be the case in the absence of hopping terms
[61]. For the DH pair, we have also evaluated the real-space superconducting
correlator (not shown), which shows a staggered long-range correlations
characteristic for $\eta$-pairing state and consistent with recent time-
dependent DMRG study [48].

###  III.3 Enhanced Binding by Floquet Drive

Increasing the exchange coupling
$J_{\parallel(\perp)}=\frac{4t_{\parallel(\perp)}^{2}}{U}$ by reducing the
interaction $U$ is one way to achieve the enhanced binding (Fig. 2), however,
this approach has limitations since we must remain in the perturbative regime
$4t_{\parallel(\perp)}\ll U$. Alternatively, tuning the rung hopping
$t_{\perp}$, resulting in the anisotropy $J_{\perp}\neq J_{\parallel}$, in
some cases leads toenhanced binding energy as well (Fig. 3), however, the
enhancement is limited by the increase in $t_{\perp}$ which tend to delocalize
the charges. Here we would like to explore whether manipulating the hopping
and the exchange coupling independently could even further increase the
stability of bounds states. Previously, this was achieved by introducing a
potential term in one of the rungs [27]. We further harness this idea by
introducing a time-periodic field and explore the Floquet engineering of the
effective model to enhance the stability of bounds states.

We drive the system with an external electric field applied along the rung and
introduced via Peierls substitution to the hopping term,
$t_{ij}(\tau)=t_{\perp}e^{i\boldsymbol{A}(\tau).\boldsymbol{r}_{ij}}$, between
nearest neighbor sites $i$ and $j$. Here, $t_{\perp}$ is the equilibrium rung
hopping parameter, $\boldsymbol{A}(\tau)$ is the vector potential and
$\boldsymbol{r}_{ij}$ is the vector from $j$ to the $i$ site. The charge of
the particles and distance between nearest neighbours is considered unity;
hence the electric field is related to the vector potential by
$\boldsymbol{E}(\tau)=-\partial\boldsymbol{A}(\tau)/\partial t$. For field in
the direction of the rung $\boldsymbol{A}(\tau)=\\{0,A\sin(\omega\tau)\\}$,
the hopping and exchange couplings along the chain ($t_{\parallel}$,
$J_{\parallel}$) remain the same as in the equilibrium. Let us write the
interaction strength as $U=l_{0}\omega+\Delta U$, where $U_{0}=l_{0}\omega$ is
the closest approximation to $U$ in multiples of $\omega$ and $\Delta U$ act
as the detuning. The Fourier components are then expressed as,

$\displaystyle t_{ij}(\tau)$
$\displaystyle=t_{\perp}\sum\limits_{l}A_{\boldsymbol{r_{ij}}}^{(l)}e^{il\omega\tau}$
$\displaystyle A_{\boldsymbol{r}_{ij}}^{(l)}$
$\displaystyle=\frac{1}{T}\int_{-T/2}^{T/2}e^{i\boldsymbol{r}_{ij}\cdot\boldsymbol{A}(\tau)}e^{-il\omega\tau}d\tau=\mathcal{J}_{l}(\boldsymbol{r}_{ij}\cdot\boldsymbol{A}),$
(4)

where $\mathcal{J}_{l}$ is the Bessel function of order $l$. Here, we have
used the relation
$\mathcal{J}_{l}(x)=\frac{1}{2T}\int_{-T}^{T}e^{i(x~{}\sin\omega\tau-l\omega\tau)}d\tau$.

Figure 5: Detuning parameter $\Delta U$ dependence on the driving frequency
$\omega$ with regions of validity for the off-resonant driving
$U,\omega,|\Delta U|\gg t_{\perp}$ shaded in green. The interaction parameter
is $U$=10.

Depending on the value of $\Delta U$ there are different regimes; far-from-
resonant driving with $U,\omega,|\Delta U|\gg t_{\perp}$ and near-resonant
driving with $U,\omega,\gg t_{\perp},|\Delta U|$, which can be treated with
appropriate high-frequency expansions [30, 29, 31, 62]. We will consider only
the off-resonant driving, in Fig. 5 denoted by the green regions, as heating
and doublon-holon recombination effects are suppressed in this regime, see Sec
V.

After performing the high-frequency expansion in the off-resonant regime and
neglecting higher order terms, see App. A, we arrive at the effective Floquet
Hamiltonian $H_{\text{eff}}^{(F)}$ of form (II), where parameters along the
chains equal to the equilibrium ones, while both hopping and exchange
parameters in the rung (field) direction depend on the drive parameters and
get modified as [30, 29, 31, 62]

$\displaystyle\mathfrak{t}_{\perp}$ $\displaystyle=$ $\displaystyle
t_{\perp}\mathcal{J}_{0}(A)$ (5) $\displaystyle\mathfrak{J}_{\perp,\eta}^{xy}$
$\displaystyle=$
$\displaystyle\sum_{l=-\infty}^{\infty}\frac{4t_{\perp}^{2}}{(U-l\omega)}A_{\boldsymbol{e}_{y}}^{(l)}A_{-\boldsymbol{e}_{y}}^{(l)}$
(6) $\displaystyle\mathfrak{J}_{\perp,\eta}^{z}$ $\displaystyle=$
$\displaystyle\mathfrak{J}_{\perp,S}=\sum_{l=-\infty}^{\infty}\frac{4t_{\perp}^{2}}{(U-l\omega)}\mid
A_{\boldsymbol{e}_{y}}^{(l)}\mid^{2}.$ (7)

Figure 6: Variation of the exchange parameters for off-resonant driving with
respect to the equilibrium value for (a) the spin channel
$J_{\perp}-\mathfrak{J}_{\perp,S}$ and (b) the $\eta$ exchange channel
$J_{\perp}-\mathfrak{J}_{\perp,\eta}^{xy}$ with driving frequency $\omega$ and
the driving strength $A$. These plots corresponds to $U=10$,
$t_{\parallel}=t_{\perp}=1$.

Applying a periodic drive thus gives the freedom to vary hopping and different
exchange parameters independently. In Fig. 6, we show the variation of (a)
$J_{\perp}-\mathfrak{J}_{\perp,S}$ and (b)
$J_{\perp}-\mathfrak{J}_{\perp,\eta}^{xy}$ as a function of the field strength
$A$ and frequency $\omega$. White regions in the plot correspond to areas that
do not obey the limit, $U,\omega,|\Delta U|\gg t_{\perp}$. Regions with
$J_{\perp}<\mathfrak{J}_{\perp,S}$ (stronger than equilibrium spin exchange
coupling) and $J_{\perp}<\mathfrak{J}_{\perp,\eta}^{xy}$ (stronger than
equilibrium $\eta$ pairing exchange coupling) are candidate regions in which
enhanced exchange couplings could mediate stronger binding between the charged
particles.

Figure 7: Binding energy in a Floquet driven setup for (a) holon-holon and (b)
doublon-holon pair as a function of driving amplitude $A$ at various driving
frequencies $\omega=4.01,6.01,8.01$. (c) Variation of the Floquet modified
superexchange in the spin channel $-\mathfrak{J}_{\perp,s}$ which is also a
2$\times$2 cluster estimate of the binding energy of holon-holon pair. (d)
2$\times$2 cluster estimation of the Floquet binding energy
$-(\mathfrak{J}_{\perp,S}+\mathfrak{J}_{\perp,\eta}^{xy})/2$ for the doublon-
holon pair. In all above cases, the horizontal lines mark the equilibrium
(non-driven) value. (e) Variation of Floquet modified hopping
$\mathfrak{t}_{\perp}$ as a function of field strenght $A$. Parameters:
$L=20$, $U=10$, $t_{\parallel}=t_{\perp}=1$.

Now, we will check how well such predictions work by analyzing the binding
energy of charge carriers within the effective Floquet Hamiltonian
$H_{\text{eff}}^{(F)}$. We will focus on the field amplitude $A$ dependence at
three frequencies, $\omega=4.01,6.01,8.01$. Figs. 7(a,b) shows the binding
energy for the HH and DH pair, respectively. In Figs. 7(b,c) we show for
comparison the binding energy estimates extracted from the 2$\times$2
clusters, for an HH pair given by $-\mathfrak{J}_{\perp,s}$ and for a DH pair
given by $-(\mathfrak{J}_{\perp,\eta}^{xy}+\mathfrak{J}_{\perp,S})/2$. Indeed,
the binding energy for HH pair is increased in the similar range $A\in(0-3.2]$
where also spin exchange coupling $\mathfrak{J}_{\perp,s}$ is increased
compared to the equilibrium value. Comparison of the DH binding energy (Fig
7(b)) with its cluster estimate (Fig 7(d)) shows that the DH binding is
enhanced in a broader range of the excitation strengths $A$ than anticipated
from the 2x2 cluster estimation. We attribute the difference to the fact that
the hopping strength $\mathfrak{t}_{\perp}$ is substantially reduced in this
regime of $A$ (Fig. 7(e)), which should be favourable for binding. Overall,
the conclusion is that the Floquet driving leads to a substantial modification
of the pair binding energies, which can be enhanced to more than three times
compared to the equilibrium values.

To further ensure that the enhanced binding energies indeed lead to a more
bound pair we consider real-space HH and DH correlators. In Fig. 8, we compare
the charge correlators for the system under the Floquet drive to its
equilibrium counterpart. The inter-chain correlators $\chi_{c}(\\{1,r\\})$
show that in the regions of the enhanced bindings (e.g. A=1.98), the charge
carriers are indeed residing at substantially closer distances both for HH and
DH pair, see Fig. 8(a) and (b). Furthermore, the intra-chain correlators
$\chi_{c}(\\{0,r\\})$ shows a strong reduction with respect to the equilibrium
value, which means that the two charge carriers will tend to form pairs in the
opposite chain along the rung. These results show that the Floquet driving can
substantially increase the tendency of charge carriers to form bound pairs.

Based on the above analysis, we find that Floquet engineering can result in
significantly enhanced binding energies due to the reduction of the hopping
integral and the enhancement of spin and $\eta$ fluctuations. At least in the
two-leg ladder system, enhancing HH binding appears in a broader regime of
drive parameters $(\omega,A)$ since it only requires increased $J_{\perp,S}$.
To increase DH binding, we preferentially want to increase both $J_{\perp,S}$
and $J_{\perp,\eta}^{xy}$, which seems to be exclusive in many cases. However,
the enhanced DH binding is still present in a broad parameter regime as the
Floquet driving decreases the effective hopping integral, leading to a higher
tendency for excitonic bindings.

Now, we analyze how these findings are translated to 2D systems.

Figure 8: Inter-chain change correlators $\chi_{c}(\\{1,r\\})$ for (a) holon-
holon and (b) doublon-holon pair for different driving frequencies
$\omega=4.01,6.01,8.01$ and in equilibrium (eq). (c,d) Same for the intra-
chain correlator. Parameters: $t_{\parallel}=t_{\perp}=1$, $U=10$, $A=1.98$
and $L=20$.

##  IV Results for two dimensional square lattice

We will consider 2D lattices with periodic boundary conditions for which the
$x$ and $y$ direction are equivalent; therefore, the $\perp$ and $\parallel$
terms in the effective model $H_{\text{eff}}$, Eq. (II), are equivalent. At
least in equilibrium, i.e., without any Floquet engineering, all the exchange
couplings are the same and equal
$J_{\parallel,S}=J_{\parallel,\eta}=J_{\perp,S}=J_{\perp,\eta}=\frac{4t^{2}}{U}\equiv
J$.

The question of HH and DH binding in 2D has already been addressed
theoretically [43, 63, 38, 64, 60, 11, 8], as of relevance for the high-
temperature superconductivity in effectively two dimensional doped cuprates
[8] and photo-doped charge transfer and Mott insulating materials [65, 66, 44,
67]. It has been recently experimentally confirmed that in photo-doped charge
transfer and Mott insulating materials, bound states of doublon-holon pairs,
i.e., Hubbard excitons, indeed appear as metastable states that form after the
pump pulse. They can be detected via a transfer of spectral weight from the
Drude peak to a finite frequency Lorenzian [67] and in the third order optical
response [44]. Another indirect clue for the presence of doublon-holon
excitons comes from an exponential functional form of the mid-infrared peak
decay in photo-doped materials [65, 66, 68], reflecting that the recombination
process of doublon-holon pair happens from a bound state. The recombination
rate is exponentially suppressed in the number of bosonic excitations that are
emitted in the process, in order to bridge the Mott gap [38, 64]. Furthermore,
it was shown that the non-local Coulomb interaction can cooperate with the
super-exchange coupling assisted excitons formation [38, 39, 67].

Here we revisit the equilibrium analysis and compare the HH and DH binding. As
opposed to the previous studies of DH Hubbard excitons, we include here the
effect of $\eta$ fluctuations, by keeping the $\eta$ term $H_{\eta}$ in the
effective Hamiltonian $H_{\text{eff}}$, Eq. (II), that we diagonalize using
Krylov technique in the sector with one holon-holon and one doublon-holon pair
on system of size up to $N=26$ sites.

###  IV.1 Equilibrium results

Figure 9: Binding energy for a doublon-holon pair (black) and holon-holon pair
(green) as a function of $J$ for $N=20,26$.

In Fig. 9, we show binding energies for holon-holon and doublon-holon pair as
a function of $J$ on systems of size $N=20,26$ sites. To minimize the finite
size effects on lattices that do not have the $\boldsymbol{k}=(\pi/2,\pi/2)$
wave vector (which is the case for $N=20,26$), we use for $H_{0}^{H}$ and
$H_{0}^{D}$ the fit obtained from $N=16,32$ in Ref. 69. As in the two-leg
ladder, increasing the exchange coupling $J$ enhances the binding energy,
since breaking the short-range antiferromagnetic spin fluctuations is
energetically penalized, leading to localization of the pair. Our results
suggest that in 2D, including the $\eta$ term $H_{\eta}$ slightly reduces the
binding energy of DH exciton, which suggests that $\eta$ fluctuations do not
help to stabilize the exciton. However, the $\eta$ contribution is not very
significant so the Hubbard excitons are stable even if this term is taken into
account. Unlike in the two-leg ladder setup, here, doublon-holon and holon-
holon binding energies are always distinct, and this is true with or without
including the $\eta$ fluctuations. We hypothesise that this can be traced back
to the difference in phases in the hopping term for holon and doublon, Eq.
(II), which is known to result in different wave function symmetry of the two
bound states: while the holon-holon bound states has $d$-wave symmetry, the
doublon-holon bound states has $s$-wave symmetry [43]. As some of us recently
discussed in relation to experimental observation of intra-excitonic
transitions [67], there actually exist several DH Hubbard excitons with
different symmetries. Here, however, we discuss only the lowest and most
stable one.

In Fig. 10(a), we plot the charge correlator $\chi_{c}^{DH}$, showing the
distribution of holon relative to the doublon (in the origin) for the Hubbard
exciton state at $J=0.4$. Results are in agreement with a bound state,
however, we can see that the excitonic state is rather extended: the largest
contribution to the excitonic wave function comes from holon and doublon being
only the third nearest neighbors. In Fig. 11(a), we plot an equivalent charge
correlator $\chi_{c}^{HH}$ for the holon-holon bound state, showing the
distribution of one holon (in the origin) relative to the other one for
$J=0.4$. Results are again in agreement with a bound state.

Figure 10: Charge correlator $\chi_{c}^{DH}$ showing the distribution of holon
relative to the doublon (in origin) for the doublon-holon exciton state (a) in
the equilibrium at $J=0.4$ ($A=0$) and (b-d) for the Floquet driving with
$A=1.0,1.5,2.0$, $\omega=6.01$ and $U=10$. The size of the dots is
proportional to the correlation on this site and suggests the presence of
binding. Periodic boundary condition considered on the $N=26$ sites lattice is
denoted by the square.

Figure 11: Charge correlator $\chi_{c}^{HH}$, showing the distribution of the
first holon relative to the second holon (in the origin) for the holon-holon
bound state (a) in the equilibrium at $J=0.4$ ($A=0$), (b-d) for the Floquet
driving with $A=1.0,1.5,2.0$, $\omega=6.01$ and $U=10$. Size of the dots is
proportional to the correlation on this site and suggests the presence of
binding. Periodic boundary condition on the considered $N=26$ sites lattice is
denoted by the square.

###  IV.2 Floquet engineering

Lastly, we revisit the Floquet engineering of exchange coupling parameters
from Sec. III.3. In the two-leg ladder, we considered the electric field along
the rung with $\boldsymbol{A}(\tau)=\\{0,A\sin(\omega t)\\}$ which modulated
only the exchange coupling along the rungs. Here, we consider the electric
field applied along the diagonal, with $\boldsymbol{A}(\tau)=\\{A\sin(\omega
t),A\sin(\omega t)\\}$. Then, the four-fold symmetry of the 2D lattice is
preserved and exchange coupling between all the nearest neighbors are the
same. We will again present results for the off-resonant regime with
$U,\omega,|\Delta U|\gg t$, where different exchange terms are given by Eq.
(6) and Eq. (7).

Figure 12: Binding energy for a (a) holon-holon pair and (b) doublon-holon
pair in the Floquet driven case as a function of driving amplitude $A$ at
$\omega=4.01,6.01,8.01$ and $U=10$ for $N=26$. The horizontal line corresponds
to the equilibrium value with $J=0.4$.

In Fig. 12, we show the variation of the binding energy for holon-holon and
doublon-holon pair as a function of $A$ for $\omega=4.01,6.01,8.01$ and
$N=26$. We choose the same frequencies as in the corresponding study in the
two-leg ladder. In this case, we use for $E_{0}^{H},E_{0}^{D}$ the energies at
the momentum with lowest g.s. energy for the given system size; the values of
effective exchange parameters vs. the hopping are outside of the validity of
fit from Ref. 69 used in the static case. Like in the ladder geometry, there
is an extended regime of field strengths for which the binding is enhanced via
the Floquet engineering of coupling parameters. This conclusion is supported
also by the DH and HH correlators, which are in Fig. 10(b,c,d) and Fig.
11(b,c,d) shown for field strengths $A=1.0,1.5,2.0$. It is clearly visible
that for a region of parameters $A$, the DH and HH pairs get more localized
than at equilibrium. In the extreme case of large driving, e.g. $A=2.0$ for DH
pair in Fig. 10, the DH pair resides at the nearest-neighbor sites in stark
contrast to the equilibrium value.

This result once again underlines that Floquet engineering is indeed a
relevant approach that can significantly increase the localization of the
bound states and should be a viable option that could be implemented in either
solid state or cold atoms setups. For the latter, we draft our suggestion in
Section VI.

##  V Recombination of Hubbard excitons

While Floquet driving can increase the binding of charged pairs by modifying
the effective hopping and exchange parameters, it also has an influence on the
doublon-holon recombination rate that determines the timescale on which
excitons are stable and observable. Namely, Floquet driving opens additional
recombination channels and therefore it is important to understand how stable
is the excitonic state in the driven case.

Due to the large energy of the Mott gap that must be instantaneously released
in the recombination process of a doublon-holon pair, the recombination
timescale is much slower than the timescale for the intrinsic
dynamics/relaxation of holons and doublons within each Hubbard band. Even
though recombination appears in the Hubbard model within the hopping term of
strength $t$, it is reasonable to treat it perturbatively: we have made this
formal in Sec. II by canonically transforming the original model, Eq. (II), so
that recombination term appears with a smaller prefactor (e.g., $t^{2}/U$) in
the transformed one. Transformation yields the effective Hamiltonian
$H_{\text{eff}}$, Eq. (II), that preserves the number of doublon-holon pairs,
using which we can elegantly estimate the metastable excitonic state.
Moreover, the transformation also singles out the dominant recombination, now
with a prefactor of the order $J$ [70, 38],

$H_{\text{rc}}=\frac{J}{2}\sum_{(ijk)ss^{\prime}}(1-n_{i,\bar{s}})c_{is}n_{ks^{\prime}}c_{k\bar{s}^{\prime}}\,\vec{\sigma}_{s\bar{s}^{\prime}}\cdot\boldsymbol{S}_{j},$
(8)

where $\bar{s}$ is spin opposite to $s$ and $i\neq k$ are nearest-neighbor
sites to $j$. In the static setup and at a low density of photoexcited DH
pairs, the only viable recombination channel is via emission of spin
excitations, which absorb the large energy of the Mott gap $\Delta$ as doublon
and holon recombine across the gap. The recombination rate can be numerically
estimated via Fermi’s golden rule expression for the transition from the
excitonic state $|\psi_{0}^{DH}\rangle$ into a highly excited state of spins
$|\psi_{m}^{0}\rangle$ with energy $E_{m}^{0}$,

$\Gamma=2\pi\sum_{m}|\langle\psi_{m}^{0}|H_{\text{rc}}|\psi_{0}^{DH}\rangle|^{2}\delta(E_{m}^{0}-E_{0}^{DH}).$
(9)

In Refs. 70, 64, some of us showed that the recombination rate $\Gamma$ is
roughy exponentially suppressed in the number of spin excitation $\Delta/J$
that are emitted in such a process. A more precise dependence on the model
parameters was derived in Ref. 71 by making a connection to the exciton-boson
model, where exciton recombines by emitting some general boson excitations. In
this case the decay rate can be approximated by

$\Gamma\approx g_{{\rm
rc}}^{2}e^{-\xi}\sqrt{\frac{2\pi}{\Delta\omega_{0}}}\exp\left[-\frac{\Delta}{\omega_{0}}\ln\left(\frac{\Delta}{e\xi\omega_{0}}\right)\right],$
(10)

where $g_{{\rm rc}}$ is the prefactor in the operator causing the
recombination, $\xi$ is the dimensionless coupling strength between the
exciton and bosons, and $\omega_{0}$ is the typical energy of bosons emitted.
For coupling to spins, $\omega_{0}=J$, $\xi$ was fitted to be roughly
$\xi\approx 3$, and $g_{{\rm rc}}$ was fitted to $g_{{\rm rc}}\approx J$,
corresponding to the actual prefactor in $H_{\text{rc}}$, Eq. (8).

In the case of Floquet driving with the original time-dependent Hubbard model,
a strategy to estimate the recombination rate is to perform a time-dependent
canonical transformation, which once again transforms out the recombination
term (from the time-dependent model) at order $t$. Inspired by Ref. 72, we
perform this transformation in App. B and obtain that in the lowest order of
the high-frequency expansion, the dominant recombination term has the same
structure as in the static case, but with a different prefactor,

$\displaystyle H_{\text{rc}}^{(F)}=\frac{1}{2}\sum_{(ijk)ss^{\prime}}\sum_{l}$
$\displaystyle\frac{t_{ij}t_{jk}\mathcal{J}_{-l}(\boldsymbol{A}\cdot\boldsymbol{r}_{ji})\mathcal{J}_{l}(\boldsymbol{A}\cdot\boldsymbol{r}_{jk})}{U-l\omega}$
$\displaystyle\times(1-n_{i,\bar{s}})c_{is}n_{ks^{\prime}}c_{k\bar{s}^{\prime}}\,\vec{\sigma}_{s\bar{s}^{\prime}}\cdot\boldsymbol{S}_{j},$
(11)

Figure 13: Upper bound on the recombination rate $\Gamma$ as a function of the
Floquet driving amplitude $A$, estimate from the expression (10) for
$\omega=6.01$ used in previous sections, c.f., Figs. 10, 11.

Most importantly, the recombination term $H_{\text{rc}}^{(F)}$ tells us that
Floquet driving opens additional recombination channels in which magnons and
photons can be emitted subsequently: emission of $l$ photons is suppressed by
$\sim\mathcal{J}_{l}(A)^{2}$ factor, but it implies that a considerably
smaller $\sim(U-l\omega)$ energy must be received by spin excitations.
Potentially, this can drastically shorten the recombination time and the
stability of the exciton. However, such formalism can not capture a process
where interference between the photon and the magnon emission would be
considered. Currently, we do not know how to treat such a process fully and we
will provide an upper bound for the recombination time assuming independent
photon/magnon emission. Since the form of the recombination operator is the
same in the static and Floquet case, we estimate the recombination rate by
expression (10). In the Floquet case, we use
$\Delta=U-\lfloor\frac{U}{\omega}\rfloor\omega$, where $\lfloor x\rfloor$ is
the floor function, as the minimal energy that must be emitted into spin
excitations,
$\omega_{0}=\sum_{l=-\infty}^{\infty}\frac{4t^{2}\mathcal{J}_{l}(A)^{2}}{U-l\omega}$
as a typical renormalized magnon energy, and $g_{{\rm
rc}}=\sum_{l=-\infty}^{\infty}\frac{4t^{2}\mathcal{J}_{l}(A)^{2}}{U-l\omega}$.
For expression (10) to be valid, $\omega_{0},g_{{\rm
rc}},\ln(\Delta/e\xi\omega_{0})>0$, which is not the case for every choice of
Floquet parameters. In addition, we should stress that this approach does not
correctly reproduce the recombination rate for the $A\to 0$ limit, when the
full Mott gap energy should be absorbed by the spins. Nonetheless, we show in
Fig. 13 this rather pessimistic upper bound estimate of the recombination rate
dependence on the field strength $A$ for $\omega=6.01$, in which case (10) is
in the valid regime. We anticipate that this result strongly overestimates the
recombination rate at small field strengths $A$. But it does also tell as that
for such a driving, the binding energy is always at least an order larger then
the width of the excitonic peak. While our study is by no means conclusive, it
open the question of better estimating the bound state’s lifetime in all
regimes, but most importantly in the regime where Floquet driving most
significantly enhances binding.

While holon-holon bound states do not suffer from recombination, the opposite
process of doublon-holon creation would plague the stability of HH binding in
the Floquet setup at the same rate $\Gamma$.

##  VI Photo-doped states in cold atoms experiments

While the existence of the doublon-holon excitons has been primarily discussed
in the context of photo-doped solid-state setups [73, 74, 75, 76, 77], cold-
atom experiments could serve as an excellent testbed to study properties of
photo-doped states. Recently, there has been tremendous progress in
understanding half-filled [21] and doped Fermi-Hubbard systems with cold atom
setups [78] exhibiting paradigmatic non-perturbative effects such as the
pseudogap [79], spin strings [80, 23], linear resistivity [81], etc. The
application of the Floquet engineering to the chemically doped systems
presented in this work can be straightforwardly extended to cold atom setups.
A new avenue is the extension to photo-doped systems and we propose simple
protocols enabling a systematic study of differences between the chemically
and photo-doped phase diagrams. We will consider two protocols to prepare the
photo-doped state: the adiabatic state preparation from band to Mott insulator
and (chirped) excitation across the Mott gap.

#### Adiabatic state preparation

In the first protocol, we start from a Mott-insulating state, where two sites
have large positive and negative on-site potential, leading to fully occupied
and empty states. Now, one can slowly reduce the on-site potential, allowing
for the virtual dressing of the empty and full state with charge and spin
fluctuations. The ramp velocity has to be slow enough so that the Landau-
Zenner tunnelling across the Mott gap is suppressed [82, 83, 70]. With this
respect, it might be convenient to first prepare the doublon-holon pair at
high values of the Mott repulsion $U$ and then reduce it to the desired value.
Such protocols can be easily extended to multiple holon-doublon pairs. This
would allow us to study a complete phase diagram of photo-doped Mott
insulators, including peculiar long-range phases such as the $\eta$-pairing
superconductivity or the charge-density wave [47, 53].

#### Excitations accross the gap

A standard approach to create holon-doublon pair is to either excite the
system with a resonant excitation or apply an electric field comparable to the
Mott gap. A naive implementation of the protocol would lead to a highly
excited state with broad distribution of holes and doublon within the Hubbard
band [84] which is detrimental for the formation of Hubbard excitons with
rather small binding energies and would require a substantial coupling with an
external heat bath. Instead, we propose using a chirped pulse (either
electromagnetic or parameter modulation) with the base frequency slightly
below the Mott gap, which is then slowly chirped across the Mott gap. It was
shown in Ref. [85, 86, 53, 86] that such a protocol creates a quasi-
equilibrium distribution with low effective temperatures and doublon-holon
excitons should emerge as a stable structure.

Upon the successful preparation of the doublon-holon state through either of
the prescribed protocols, one can initiate the application of the proposed
Floquet engineering technique (Eq. (III.3)) to manipulate the exciton binding
energy. Analogous superexchange manipulation can be obtained by directly
modifying hopping [32].

##  VII Conclusions

In this work, we compared the pairing of chemically and photo-doped charge
carriers mediated by spin and $\eta$ fluctuations in ladder and 2D systems. As
the photo-doped state is metastable, we employed the Schrieffer-Wolf
transformation to approximate the long-lived phase as a quasi-equilibrium
state described by an extended t-J model. We establish that the binding is
comparable in the two regimes by analysing binding energies and real-space
correlators. In equilibrium and within the single-band Hubbard model, the
superexchange $J$ and the hopping $t$ are connected and while the increase of
the superexchange would enhance charge pairing it naturally comes with the
increase of the hopping, leading to the tendency to delocalize charge
carriers. We used Floquet engineering to decouple the two and the analysis
yields an effective Floquet Hamiltonian, where both the superexchange
$\mathfrak{J}$ and the hopping $\mathfrak{t}$ are modified by the driving
frequency $\omega$ and the amplitude $A$. We demonstrate that a judicious
selection of drive parameters can substantially enhance the strength of
binding between exciton pairs both for ladders and 2D systems, although the
increase in the latter is limited in a narrower driving range. The tendency to
form bound pairs is most obvious in the analysis of the real-space
correlators, which show a dramatic reduction in the relative distance between
pairs. Finally, we provided an upper bound on the lifetime of the metastable
phase using a Fermi golden rule argument and commented on how our theoretical
predictions could be analyzed in cold-atom experiments.

Our estimation of the lifetime provides only an upper and most probably pretty
pessimistic limit for the lifetime of the holon-doublon pair. Extending the
description, which would consider both possible interferences between photon
and magnon-assisted recombination and properly capture the limit of weak
driving, is an important future problem. One direction would be a simulation
of the full Hubbard model prepared with a holon-doublon pair and the inclusion
of periodic driving. An alternative approach is to generalize the Floquet
Fermi golden rule [87] approach for photo-doped Mott insulators, where the
number of holons and doublons present an almost conserved quantity.

In the current work, we have focused on a dilute limit with a single pair;
however, the formation of long-ranged orders will depend on the interaction
between bound pairs. An exciting new perspective is to understand what is the
effective interaction between holon-doublon pairs and what would be the
symmetry of the $\eta$-paired condensate in the ladder and 2D systems. As the
Floquet driving can change the pairing of bound pairs, it will also have an
effect on the interaction between pairs and a systematic analysis of these
processes is needed to understand how to stabilize non-thermal states with
long-range orders.

###### Acknowledgements.

We thank Y. Murakami, T. Kaneko, M. Bukov, and M. Eckstein for several useful
discussions. We acknowledge the support by the projects J1-2463, N1-0318,
MN-0016-106 and P1-0044 program of the Slovenian Research Agency, the QuantERA
grant T-NiSQ by MVZI, QuantERA II JTC 2021, and ERC StG 2022 project DrumS,
Grant Agreement 101077265. ED calculations were performed at the cluster
‘spinon’ of JSI, Ljubljana.

## References

  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * Warren  _et al._ [1989] W. W. Warren, R. E. Walstedt, G. F. Brennert, R. J. Cava, R. Tycko, R. F. Bell, and G. Dabbagh, Phys. Rev. Lett. 62, 1193 (1989).
  * Norman  _et al._ [2005] M. R. Norman, D. Pines, and C. Kallin, Advances in Physics 54, 715 (2005).
  * Wu  _et al._ [2018] W. Wu, M. S. Scheurer, S. Chatterjee, S. Sachdev, A. Georges, and M. Ferrero, Phys. Rev. X 8, 021048 (2018).
  * Tranquada  _et al._ [1995] J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, and S. Uchida, Nature 375, 561 (1995).
  * Corboz  _et al._ [2014] P. Corboz, T. M. Rice, and M. Troyer, Phys. Rev. Lett. 113, 046402 (2014).
  * Wietek  _et al._ [2021] A. Wietek, Y.-Y. He, S. R. White, A. Georges, and E. M. Stoudenmire, Phys. Rev. X 11, 031007 (2021).
  * Dagotto [1994] E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
  * Maier  _et al._ [2006] T. A. Maier, M. S. Jarrell, and D. J. Scalapino, Phys. Rev. Lett. 96, 047005 (2006).
  * Maier  _et al._ [2008] T. A. Maier, D. Poilblanc, and D. J. Scalapino, Phys. Rev. Lett. 100, 237001 (2008).
  * Jaklič and Prelovšek [2000] J. Jaklič and P. Prelovšek, Advances in Physics 49, 1 (2000).
  * Trugman [1988] S. A. Trugman, Phys. Rev. B 37, 1597 (1988).
  * Bonča  _et al._ [1989] J. Bonča, P. Prelovšek, and I. Sega, Phys. Rev. B 39, 7074 (1989).
  * Chernyshev  _et al._ [1998] A. L. Chernyshev, P. W. Leung, and R. J. Gooding, Phys. Rev. B 58, 13594 (1998).
  * Bonča  _et al._ [2007] J. Bonča, S. Maekawa, and T. Tohyama, Phys. Rev. B 76, 035121 (2007).
  * Ogata and Shiba [1990] M. Ogata and H. Shiba, Phys. Rev. B 41, 2326 (1990).
  * SHIBA and OGATA [1991] H. SHIBA and M. OGATA, International Journal of Modern Physics B 05, 31 (1991).
  * Giamarchi [2003] T. Giamarchi,  _Quantum physics in one dimension_ , Vol. 121 (Clarendon press, 2003).
  * Essler  _et al._ [2005] F. H. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin,  _The one-dimensional Hubbard model_ (Cambridge University Press, 2005).
  * Dagotto  _et al._ [1992] E. Dagotto, J. Riera, and D. Scalapino, Phys. Rev. B 45, 5744 (1992).
  * Mazurenko  _et al._ [2017] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, Nature 545, 462 (2017).
  * Koepsell  _et al._ [2019] J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross, Nature 572, 358 (2019).
  * Chiu  _et al._ [2019] C. S. Chiu, G. Ji, A. Bohrdt, M. Xu, M. Knap, E. Demler, F. Grusdt, M. Greiner, and D. Greif, Science 365, 251 (2019), https://www.science.org/doi/pdf/10.1126/science.aav3587 .
  * Bohrdt  _et al._ [2019] A. Bohrdt, C. S. Chiu, G. Ji, M. Xu, D. Greif, M. Greiner, E. Demler, F. Grusdt, and M. Knap, Nature Physics 15, 921 (2019).
  * Trotzky  _et al._ [2008] S. Trotzky, P. Cheinet, S. Fölling, M. Feld, U. Schnorrberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, and I. Bloch, Science 319, 295 (2008), https://www.science.org/doi/pdf/10.1126/science.1150841 .
  * Dimitrova  _et al._ [2020] I. Dimitrova, N. Jepsen, A. Buyskikh, A. Venegas-Gomez, J. Amato-Grill, A. Daley, and W. Ketterle, Phys. Rev. Lett. 124, 043204 (2020).
  * Hirthe  _et al._ [2023] S. Hirthe, T. Chalopin, D. Bourgund, P. Bojović, A. Bohrdt, E. Demler, F. Grusdt, I. Bloch, and T. A. Hilker, Nature 613, 463 (2023).
  * Bohrdt  _et al._ [2022] A. Bohrdt, L. Homeier, I. Bloch, E. Demler, and F. Grusdt, Nature Physics 18, 651 (2022).
  * Bukov  _et al._ [2015] M. Bukov, L. D’Alessio, and A. Polkovnikov, Advances in Physics 64, 139 (2015), https://doi.org/10.1080/00018732.2015.1055918 .
  * Bukov  _et al._ [2016] M. Bukov, M. Kolodrubetz, and A. Polkovnikov, Phys. Rev. Lett. 116, 125301 (2016).
  * Mentink  _et al._ [2015] J. H. Mentink, K. Balzer, and M. Eckstein, Nature Communications 6, 6708 (2015).
  * Görg  _et al._ [2018] F. Görg, M. Messer, K. Sandholzer, G. Jotzu, R. Desbuquois, and T. Esslinger, Nature 553, 481 (2018).
  * Murakami  _et al._ [2023a] Y. Murakami, D. Golež, M. Eckstein, and P. Werner, Photo-induced nonequilibrium states in mott insulators (2023a), arXiv:2310.05201 [cond-mat.str-el] .
  * de la Torre  _et al._ [2021] A. de la Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef, Rev. Mod. Phys. 93, 041002 (2021).
  * Sensarma  _et al._ [2010] R. Sensarma, D. Pekker, E. Altman, E. Demler, N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz, and T. Esslinger, Phys. Rev. B 82, 224302 (2010).
  * Strohmaier  _et al._ [2010] N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz, T. Esslinger, R. Sensarma, D. Pekker, E. Altman, and E. Demler, Phys. Rev. Lett. 104, 080401 (2010).
  * Kollath  _et al._ [2007] C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007).
  * Lenarčič and Prelovšek [2013] Z. Lenarčič and P. Prelovšek, Phys. Rev. Lett. 111, 016401 (2013).
  * Bittner  _et al._ [2020] N. Bittner, D. Golež, M. Eckstein, and P. Werner, Phys. Rev. B 101, 085127 (2020).
  * Sugimoto and Ejima [2023] K. Sugimoto and S. Ejima, Pump-probe spectroscopy of the one-dimensional extended hubbard model at half filling (2023), arXiv:2305.09909 [cond-mat.str-el] .
  * Jeckelmann [2003] E. Jeckelmann, Phys. Rev. B 67, 075106 (2003).
  * Matsueda  _et al._ [2004] H. Matsueda, T. Tohyama, and S. Maekawa, Phys. Rev. B 70, 033102 (2004).
  * Tohyama [2006] T. Tohyama, Journal of the Physical Society of Japan 75, 034713 (2006).
  * Terashige  _et al._ [2019] T. Terashige, T. Ono, T. Miyamoto, T. Morimoto, H. Yamakawa, N. Kida, T. Ito, T. Sasagawa, T. Tohyama, and H. Okamoto, Science Advances 5, eaav2187 (2019), https://www.science.org/doi/pdf/10.1126/sciadv.aav2187 .
  * Rosch  _et al._ [2008] A. Rosch, D. Rasch, B. Binz, and M. Vojta, Phys. Rev. Lett. 101, 265301 (2008).
  * Kaneko  _et al._ [2019] T. Kaneko, T. Shirakawa, S. Sorella, and S. Yunoki, Phys. Rev. Lett. 122, 077002 (2019).
  * Murakami  _et al._ [2022] Y. Murakami, S. Takayoshi, T. Kaneko, Z. Sun, D. Golež, A. J. Millis, and P. Werner, Communications Physics 5, 23 (2022).
  * Ueda  _et al._ [2023] R. Ueda, K. Kuroki, and T. Kaneko, Photoinduced $\eta$-pairing correlation in the hubbard ladder (2023), arXiv:2310.10153 [cond-mat.str-el] .
  * Kaneko  _et al._ [2020] T. Kaneko, S. Yunoki, and A. J. Millis, Phys. Rev. Res. 2, 032027 (2020).
  * Ejima  _et al._ [2020] S. Ejima, T. Kaneko, F. Lange, S. Yunoki, and H. Fehske, Phys. Rev. Res. 2, 032008 (2020).
  * Tindall  _et al._ [2019] J. Tindall, B. Buča, J. R. Coulthard, and D. Jaksch, Phys. Rev. Lett. 123, 030603 (2019).
  * Shirakawa  _et al._ [2020] T. Shirakawa, S. Miyakoshi, and S. Yunoki, Phys. Rev. B 101, 174307 (2020).
  * Li  _et al._ [2020] J. Li, D. Golez, P. Werner, and M. Eckstein, Phys. Rev. B 102, 165136 (2020).
  * Murakami  _et al._ [2023b] Y. Murakami, S. Takayoshi, T. Kaneko, A. M. Läuchli, and P. Werner, Phys. Rev. Lett. 130, 106501 (2023b).
  * Eckstein and Werner [2011] M. Eckstein and P. Werner, Phys. Rev. B 84, 035122 (2011).
  * MacDonald  _et al._ [1988] A. H. MacDonald, S. M. Girvin, and D. Yoshioka, Phys. Rev. B 37, 9753 (1988).
  * Prelovšek and Bonča [2013] P. Prelovšek and J. Bonča, in  _Springer Series in Solid-State Sciences_ (Springer Berlin Heidelberg, 2013) pp. 1–30.
  * Lanczos [1950] C. Lanczos, Journal of Research of the National Bureau of Standards 45, 255 (1950).
  * Troyer  _et al._ [1996] M. Troyer, H. Tsunetsugu, and T. M. Rice, Phys. Rev. B 53, 251 (1996).
  * Shinjo  _et al._ [2021] K. Shinjo, S. Sota, and T. Tohyama, Phys. Rev. B 103, 035141 (2021).
  * White and Scalapino [1997] S. R. White and D. J. Scalapino, Phys. Rev. B 55, 6504 (1997).
  * Murakami  _et al._ [2023c] Y. Murakami, M. Schüler, R. Arita, and P. Werner, Phys. Rev. B 108, 035151 (2023c).
  * Takahashi  _et al._ [2002] M. Takahashi, T. Tohyama, and S. Maekawa, Phys. Rev. B 66, 125102 (2002).
  * Lenarčič and Prelovšek [2014] Z. Lenarčič and P. Prelovšek, Phys. Rev. B 90, 235136 (2014).
  * Okamoto  _et al._ [2010] H. Okamoto, T. Miyagoe, K. Kobayashi, H. Uemura, H. Nishioka, H. Matsuzaki, A. Sawa, and Y. Tokura, Phys. Rev. B 82, 060513 (2010).
  * Okamoto  _et al._ [2011] H. Okamoto, T. Miyagoe, K. Kobayashi, H. Uemura, H. Nishioka, H. Matsuzaki, A. Sawa, and Y. Tokura, Phys. Rev. B 83, 125102 (2011).
  * Mehio  _et al._ [2023] O. Mehio, X. Li, H. Ning, Z. Lenarčič, Y. Han, M. Buchhold, Z. Porter, N. J. Laurita, S. D. Wilson, and D. Hsieh, Nature Physics , 1 (2023).
  * Sahota  _et al._ [2019] D. G. Sahota, R. Liang, M. Dion, P. Fournier, H. A. Dąbkowska, G. M. Luke, and J. S. Dodge, Phys. Rev. Res. 1, 033214 (2019).
  * Leung and Gooding [1995] P. W. Leung and R. J. Gooding, Phys. Rev. B 52, R15711 (1995).
  * Lenarčič and Prelovšek [2012] Z. Lenarčič and P. Prelovšek, Phys. Rev. Lett. 108, 196401 (2012).
  * Karakonstantakis  _et al._ [2011] G. Karakonstantakis, E. Berg, S. R. White, and S. A. Kivelson, Phys. Rev. B 83, 054508 (2011).
  * Kitamura  _et al._ [2017] S. Kitamura, T. Oka, and H. Aoki, Phys. Rev. B 96, 014406 (2017).
  * Kishida  _et al._ [2001] H. Kishida, M. Ono, K. Miura, H. Okamoto, M. Izumi, T. Manako, M. Kawasaki, Y. Taguchi, Y. Tokura, T. Tohyama, K. Tsutsui, and S. Maekawa, Phys. Rev. Lett. 87, 177401 (2001).
  * Ono  _et al._ [2004] M. Ono, K. Miura, A. Maeda, H. Matsuzaki, H. Kishida, Y. Taguchi, Y. Tokura, M. Yamashita, and H. Okamoto, Phys. Rev. B 70, 085101 (2004).
  * Lu  _et al._ [2015] H. Lu, C. Shao, J. Bonča, D. Manske, and T. Tohyama, Phys. Rev. B 91, 245117 (2015).
  * Novelli  _et al._ [2012] F. Novelli, D. Fausti, J. Reul, F. Cilento, P. H. M. van Loosdrecht, A. A. Nugroho, T. T. M. Palstra, M. Grüninger, and F. Parmigiani, Phys. Rev. B 86, 165135 (2012).
  * Rincón and Feiguin [2021] J. Rincón and A. E. Feiguin, Phys. Rev. B 104, 085122 (2021).
  * Bohrdt  _et al._ [2021] A. Bohrdt, L. Homeier, C. Reinmoser, E. Demler, and F. Grusdt, Annals of Physics 435, 168651 (2021).
  * Brown  _et al._ [2020] P. T. Brown, E. Guardado-Sanchez, B. M. Spar, E. W. Huang, T. P. Devereaux, and W. S. Bakr, Nature Physics 16, 26 (2020).
  * Xu  _et al._ [2023] M. Xu, L. H. Kendrick, A. Kale, Y. Gang, G. Ji, R. T. Scalettar, M. Lebrat, and M. Greiner, Nature , 1 (2023).
  * Brown  _et al._ [2019] P. T. Brown, D. Mitra, E. Guardado-Sanchez, R. Nourafkan, A. Reymbaut, C.-D. Hébert, S. Bergeron, A.-M. S. Tremblay, J. Kokalj, D. A. Huse, P. Schauß, and W. S. Bakr, Science 363, 379 (2019), https://www.science.org/doi/pdf/10.1126/science.aat4134 .
  * Oka and Aoki [2005] T. Oka and H. Aoki, Phys. Rev. Lett. 95, 137601 (2005).
  * Oka and Aoki [2010] T. Oka and H. Aoki, Phys. Rev. B 81, 033103 (2010).
  * Oka [2012] T. Oka, Phys. Rev. B 86, 075148 (2012).
  * Werner  _et al._ [2019a] P. Werner, M. Eckstein, M. Müller, and G. Refael, Nature Communications 10, 5556 (2019a).
  * Werner  _et al._ [2019b] P. Werner, J. Li, D. Golež, and M. Eckstein, Phys. Rev. B 100, 155130 (2019b).
  * Ikeda and Polkovnikov [2021] T. N. Ikeda and A. Polkovnikov, Phys. Rev. B 104, 134308 (2021).

##  Appendix A Effective hamiltonian for the off-resonant regime

Here, we show the various steps involved in arriving at the effective
Hamiltonian $H_{\text{eff}}^{(F)}$ with the parameters shown in Eqs. (5-7). In
the off-resonant regime, we first apply a unitary transformation by
$\mathcal{U}(\tau)=\exp{(-iU\tau\sum_{j}n_{j,\uparrow}n_{j,\downarrow})}$. The
resultant Hamiltonian in the rotated frame is given by $H^{{\rm
rot}}(\tau)=\mathcal{U}^{\dagger}(\tau)H(\tau)\mathcal{U}(\tau)+i(\delta_{\tau}\mathcal{U}^{\dagger}(\tau))\mathcal{U}(\tau)$).
$H^{{\rm rot}}(\tau)$ can be written as a sum of terms that preserve the
number of DH pairs ($\tilde{T}_{0}(\tau)$), increase or decrease the DH pair
number by one ($\tilde{T}_{1}(\tau)$, $\tilde{T}_{-1}(\tau)$).

$\displaystyle H^{{\rm
rot}}(\tau)=-(\tilde{T}_{0}(\tau)+\tilde{T}_{1}(\tau)+\tilde{T}_{-1}(\tau))$
(12)

Introducing,
$g_{ij\sigma}=(1-n_{i\bar{\sigma}})c_{i\sigma}^{\dagger}c_{j\sigma}(1-n_{j\bar{\sigma}})+n_{i\bar{\sigma}}c_{i\sigma}^{\dagger}c_{j\sigma}n_{j\bar{\sigma}}$
and
$h_{ij\sigma}^{\dagger}=n_{i\bar{\sigma}}c_{i\sigma}^{\dagger}c_{j\sigma}(1-n_{j\bar{\sigma}})$,
the various terms in $H^{{\rm rot}}(\tau)$ can be written as

$\begin{split}\tilde{T}_{0}(\tau)&=\sum_{\langle
i,j\rangle,\sigma}t_{ij}(\tau)g_{ij\sigma}\\\
\tilde{T}_{1}(\tau)&=\sum_{\langle
i,j\rangle,\sigma}t_{ij}(\tau)e^{iU\tau}h_{ij\sigma}^{\dagger}\\\
\tilde{T}_{-1}(\tau)&=\sum_{\langle
i,j\rangle,\sigma}t_{ij}^{*}(\tau)e^{-iU\tau}h_{ij\sigma}\end{split}$ (13)

Using $t_{ij}(\tau)=t_{ij}\sum_{l}A_{\mathbf{r}_{ij}}^{(l)}e^{il\omega\tau}$
and Eq. (13) in Eq. (12),

$\displaystyle H^{{\rm rot}}(\tau)$ $\displaystyle=-t_{ij}\sum_{\langle
i,j\rangle,\sigma}\sum_{l}\Big{[}A_{\mathbf{r}_{ij}}^{(l)}e^{il\omega\tau}g_{ij\sigma}$
(14)
$\displaystyle+A_{\mathbf{r}_{ij}}^{(l)}e^{il\omega\tau}e^{iU\tau}h_{ij\sigma}^{\dagger}+A_{\mathbf{r}_{ij}}^{(l)*}e^{-il\omega\tau}e^{-iU\tau}h_{ij\sigma})\Big{]}$

Notice that $H^{{\rm rot}}(\tau)$ has two defining frequencies, U and
$\omega$. To make $H^{{\rm rot}}(\tau)$ periodic in both U and $\omega$ we
define a common frequency $\Omega_{0}$, such that $U=l_{0}\Omega_{0}$ and
$\omega=k_{0}\Omega_{0}$, where $l_{0}$ and $k_{0}$ are co-prime numbers. To
calculate the effective Floquet Hamiltonian we do a high-frequency expansion
of $H^{{\rm rot}}(\tau)$ by expanding it with respect to different orders of
the period $T=1/\Omega_{0}$. The first term in the expansion has the form
[30],

$\displaystyle H_{{\rm eff}}^{(0)}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\int_{0}^{T}H^{{\rm rot}}(\tau_{1})d\tau_{1}$
$\displaystyle=$ $\displaystyle H_{{\rm eff},1}^{(0)}+H_{{\rm
eff},2}^{(0)}+H_{{\rm eff},2}^{(0)\dagger}$

Now using Eq. (14) in Eq. (A), the terms of the effective Floquet Hamiltonian
at zeroth order are,

$\begin{split}H_{{\rm
eff},1}^{(0)}&=-\frac{1}{T}\int_{0}^{T}d\tau_{1}\tilde{T}_{0}(\tau_{1})\\\
&=-\sum_{\langle
i,j\rangle,\sigma}\frac{t_{ij}}{T}\int_{0}^{T}\sum_{l}A_{\mathbf{r}_{ij}}^{(l)}e^{il\omega\tau_{1}}g_{ij\sigma}d\tau_{1}\\\
&=-\sum_{\langle
i,j\rangle,\sigma}t_{ij}A_{\mathbf{r}_{ij}}^{(0)}g_{ij\sigma}\end{split}$ (16)

$\displaystyle H_{{\rm eff},2}^{(0)}$
$\displaystyle=-\frac{1}{T}\int_{0}^{T}d\tau_{1}\tilde{T}_{1}(\tau)$
$\displaystyle=-\sum_{\langle
i,j\rangle,\sigma}\frac{t_{ij}}{T}\int_{0}^{T}\sum_{l}A_{\mathbf{r}_{ij}}^{(l)}e^{il\omega\tau_{1}}e^{iU\tau_{1}}h_{ij\sigma}^{\dagger}d\tau_{1}$
$\displaystyle=-\sum_{\langle
i,j\rangle,\sigma}t_{ij}\sum_{l}A_{\mathbf{r}_{ij}}^{(l)}h_{ij\sigma}^{\dagger}\frac{1}{2\pi}\int_{0}^{2\pi}e^{i(lk_{0}+l_{0})\Omega_{0}\tau_{1}}d(\Omega_{0}\tau_{1})$
$\displaystyle=-\sum_{\langle
i,j\rangle,\sigma}t_{ij}\sum_{l}A_{\mathbf{r}_{ij}}^{(l)}h_{ij\sigma}^{\dagger}\delta_{lk_{0},-l_{0}}$
$\displaystyle=0$ (17)

Since $l_{0}$ and $k_{0}$ are co-prime numbers, $l_{0}\neq lk_{0}$, and the
above integral is zero. Similarly for the complex conjugate part $H_{{\rm
eff},2}^{(0)\dagger}=0$. Hence, at the zeroth order in the high-frequency
expansion the effective Hamiltonian is,

$\displaystyle H_{{\rm eff}}^{(0)}=\sum_{\langle
i,j\rangle,\sigma}t_{ij}A_{\mathbf{r}_{ij}}^{(0)}g_{ij\sigma}$ (18)

The next order term is given by

$\begin{split}H_{{\rm
eff}}^{(1)}&=\frac{1}{2iT}\int_{0}^{T}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}f(\tau_{1}-\tau_{2})\left[H^{{\rm
rot}}(\tau_{1}),H^{{\rm rot}}(\tau_{2})\right]\\\ &=H_{{\rm
eff},1}^{(1)}+H_{{\rm eff},2}^{(1)}+H_{{\rm eff},2}^{(1)~{}\dagger}+H_{{\rm
eff},3}^{(1)}\end{split},$ (19)

where $f(x)=1-2x/T$ is understood periodic with a period $T.$ The first term
is given by

$\begin{split}H_{{\rm
eff},1}^{(1)}&=\frac{1}{2iT}\int_{0}^{T}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}~{}f(\tau_{2}-\tau_{1})\left[\tilde{T}_{0}(\tau_{1}),\tilde{T}_{0}(\tau_{2})\right]\\\
&=\frac{1}{2iT}\int_{0}^{T}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}~{}f(\tau_{2}-\tau_{1})\left[\sum_{\langle
i,j\rangle,\sigma,l}t_{ij}A_{\mathbf{r}_{ij}}^{(l)}e^{il\omega\tau_{1}}g_{ij\sigma},\sum_{\langle
i^{\prime},j^{\prime}\rangle,\sigma^{\prime},l^{\prime}}t_{i^{\prime}j^{\prime}}A_{\mathbf{r}_{i^{\prime}j^{\prime}}}^{(l^{\prime})}e^{il^{\prime}\omega\tau_{2}}g_{i^{\prime}j^{\prime}\sigma^{\prime}}\right]\\\
&=\sum_{l>0}\frac{1}{l\Omega_{0}}\left[\sum_{\langle
i,j\rangle,\sigma}t_{ij}A_{\mathbf{r}_{ij}}^{(l)}g_{ij\sigma},\sum_{\langle
i^{\prime},j^{\prime}\rangle,\sigma^{\prime}}t_{i^{\prime}j^{\prime}}A_{\mathbf{r}_{i^{\prime}j^{\prime}}}^{(-l)}g_{i^{\prime}j^{\prime}\sigma^{\prime}}\right]=0\\\
\end{split}$ (20)

since all the commutators vanish. The second term is given by

$\begin{split}H_{{\rm
eff},2}^{(1)}&=\frac{1}{2iT^{2}}\int_{0}^{T}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}~{}f(\tau_{2}-\tau_{1})\left[\tilde{T}_{0}(\tau_{1}),\tilde{T}_{1}(\tau_{2})\right]\\\
&=\frac{1}{2iT^{2}}\int_{0}^{T}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}~{}f(\tau_{2}-\tau_{1})\left[\sum_{\langle
i,j\rangle,\sigma,l}t_{ij}A_{\mathbf{r}_{ij}}^{(l)}e^{il\omega\tau_{1}}g_{ij\sigma},\sum_{\langle
i^{\prime},j^{\prime}\rangle,\sigma^{\prime},l^{\prime}}t_{i^{\prime}j^{\prime}}A_{\mathbf{r}_{i^{\prime}j^{\prime}}}^{(l^{\prime})}e^{-i(l^{\prime}\omega+U)\tau_{2}}h_{i^{\prime}j^{\prime}\sigma^{\prime}}\right],\end{split}$
(21)

which rewritten in the frequency domain is,

$\begin{split}H_{{\rm eff},2}^{(1)}&=\sum_{\langle
i,j\rangle,\sigma}\sum_{\langle
i^{\prime},j^{\prime}\rangle,\sigma^{\prime}}\sum_{m}\frac{1}{m\Omega_{0}}A_{\mathbf{r}_{ij}}^{(m)}A_{\mathbf{r}_{i^{\prime}j^{\prime}}}^{(m-l_{0})}\left[t_{ij}g_{ij\sigma},t_{i^{\prime}j^{\prime}}h_{i^{\prime}j^{\prime}\sigma^{\prime}}\right]\end{split}$
(22)

$A_{\mathbf{r}_{ij}}^{(m)}$ is non-zero only when $m=lk_{0}$, and
$A_{\mathbf{r}_{ij}}^{(m-l_{0})}$ is non-zero only when
$m=l^{\prime}k_{0}+l_{0}$. This gives, $lk_{0}=l^{\prime}k_{0}+l_{0}$ or
$(l-l^{\prime})k_{0}=l_{0}$, which is never satisfied since $l_{0}$ and
$k_{0}$ are co-primes. As a result, $A_{\mathbf{r}_{ij}}^{(m)}$ and
$A_{\mathbf{r}_{ij}}^{(m-l_{0})}$ is never simultaneously non-zero and
$H_{{\rm eff},2}^{(1)}=0$ always. Similarly, the complex conjugate $H_{{\rm
eff},2}^{(1)\dagger}=0$.

The only non-zero integral at this order is,

$\begin{split}H_{{\rm
eff},3}^{(1)}&=\frac{1}{2iT}\int_{0}^{T}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}~{}f(\tau_{2}-\tau_{1})\left[\tilde{T}_{1}(\tau_{1}),\tilde{T}_{-1}(\tau_{2})\right]\\\
&=\frac{1}{2iT}\int_{0}^{T}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}~{}f(\tau_{2}-\tau_{1})\left[\sum_{\langle
i,j\rangle,\sigma,l}t_{ij}A_{\mathbf{r}_{ij}}^{(l)}e^{i(l\omega+U)\tau_{1}}h_{ij\sigma}^{\dagger},\sum_{\langle
i^{\prime},j^{\prime}\rangle,\sigma^{\prime},l^{\prime}}t_{i^{\prime}j^{\prime}}A_{\mathbf{r}_{i^{\prime}j^{\prime}}}^{(l^{\prime})*}e^{-i(l^{\prime}\omega+U)\tau_{2}}h_{i^{\prime}j^{\prime}\sigma^{\prime}}\right]\\\
&=\sum_{l}\frac{1}{(U+l\omega)}\left[\sum_{\langle
i,j\rangle,\sigma}t_{ij}A_{\mathbf{r}_{ij}}^{(l)}h_{ij\sigma}^{\dagger},\sum_{\langle
i^{\prime},j^{\prime}\rangle,\sigma^{\prime}}t_{i^{\prime}j^{\prime}}A_{\mathbf{r}_{i^{\prime}j^{\prime}}}^{(l)*}h_{i^{\prime}j^{\prime}\sigma^{\prime}}\right]\\\
&=(-1)\sum_{\langle
i,j\rangle,\sigma}\sum_{l}\frac{4t_{ij}^{2}}{U-l\omega}A_{\mathbf{r}_{ij}}^{-l}A_{-\mathbf{r}_{ij}}^{-l}\left(\mathbf{\eta}_{i}^{x}\mathbf{\eta}_{j}^{x}+\mathbf{\eta}_{i}^{y}\mathbf{\eta}_{j}^{y}\right)+\sum_{\langle
i,j\rangle,\sigma}\sum_{l}\frac{4t_{ij}^{2}}{U-l\omega}\mid
A_{\mathbf{r}_{ij}}^{(-l)}\mid^{2}\left(\mathbf{S}_{i}\mathbf{S}_{j}-\mathbf{\eta}_{i}^{z}\mathbf{\eta}_{j}^{z}\right)\end{split}$
(23)

So the effective Floquet Hamiltonian from Eq. (16) and Eq. (23) is,

$\begin{split}H_{\text{eff}}^{(F)}&=\sum_{\langle
i,j\rangle,\sigma}t_{ij}A_{\mathbf{r}_{ij}}^{(0)}g_{ij\sigma}-\sum_{\langle
i,j\rangle,\sigma}\sum_{l}\frac{4t_{ij}^{2}}{U-l\omega}A_{\mathbf{r}_{ij}}^{-l}A_{-\mathbf{r}_{ij}}^{-l}\left(\mathbf{\eta}_{i}^{x}\mathbf{\eta}_{j}^{x}+\mathbf{\eta}_{i}^{y}\mathbf{\eta}_{j}^{y}\right)+\sum_{\langle
i,j\rangle,\sigma}\sum_{l}\frac{4t_{ij}^{2}}{U-l\omega}\mid
A_{\mathbf{r}_{ij}}^{(-l)}\mid^{2}\left(\mathbf{S}_{i}\mathbf{S}_{j}-\mathbf{\eta}_{i}^{z}\mathbf{\eta}_{j}^{z}\right)\end{split}$
(24)

##  Appendix B Dominant recombination term in Floquet driven setup

To establish an approximate recombination rate in the Floquet driven setup
with the help of the static expression, Eq. (10), we perform a different time
dependent canonical transformation which, similar to the static case,
transforms out recombination term in the lowest order but retains it in the
next order.

Starting from the original Hamiltonian

$\begin{split}H(\tau)&=-\sum_{\langle
i,j\rangle,\sigma}t_{ij}(\tau)(c^{\dagger}_{i\sigma}c_{j\sigma}+c^{\dagger}_{j\sigma}c_{i\sigma})+U\sum_{j}n_{j\uparrow}n_{j\downarrow}\\\
&\equiv-T(\tau)+UD,\end{split}$ (25)

where we again split the hopping term into contributions that either retain,
increase or decrease the number of DH pairs,
$T(\tau)=\sum_{l}(T_{0,l}+T_{1,l}+T_{-1,l})e^{il\omega\tau}$ and define
$D=\sum_{j}n_{j\uparrow}n_{j\downarrow}$ as the double occupancy operator.

Following Ref. 72 we introduce a time-dependent canonical transformation
$e^{iS(\tau)}$, and obtain the transformed Hamiltonian $H_{{\rm r}}$ as,

$\displaystyle H_{{\rm
r}}(\tau)+UD=e^{iS(\tau)}\left[H(\tau)-i\partial_{\tau}\right]e^{-iS(\tau)}$
(26)

$S(\tau)$ is chosen to eliminate the terms off-diagonal in $D$. To determine
$H_{{\rm r}}(\tau)$ and $S(\tau)$ order by order we expand them in powers of
$\lambda$ as, $H_{{\rm r}}(\tau)=\sum_{n=1}^{\infty}\lambda^{n}H_{{\rm
r}}^{(n)}$ and $S(\tau)=\sum_{n=1}^{\infty}\lambda^{n}S^{(n)}$. We can re-
write Eq. (26) as

$\begin{split}H_{{\rm
r}}(\tau)&+U\left[D,iS(\tau)\right]+i\partial_{\tau}S(\tau)\\\ &=-\lambda
T(\tau)-\sum_{n=0}^{\infty}\frac{B_{n}}{n!}{\rm
ad}^{n}_{iS}\left[(-1)^{(n)}\lambda T(\tau)+H_{{\rm
r}}^{(n)}\right].\end{split}$ (27)

We truncate $S(\tau)=\lambda S^{(1)}$, since our aim is to keep recombination
term in the second order and compare with the equilibrium analysis to get an
upper bound for the recombination rate in the Floquet driven scenario. From
Eq. (27),

$\begin{split}H_{{\rm
r}}^{(1)}&=-T-U\left[D,iS^{(1)}\right]-\partial_{\tau}S^{(1)}\end{split}$ (28)

we find that
$iS^{(1)}=-\sum_{l}(\frac{T_{1,l}}{U-l\omega}+\frac{T_{-1,l}}{U+l\omega})e^{il\omega\tau}$,
which gives $H_{{\rm r}}^{(1)}=-\sum_{l}T_{0,l}e^{-il\omega\tau}$, i.e., it
eliminates the recombination in the lowest order.

From Eq. (27) using $S^{(2)}=0$, we get

$\begin{split}H_{{\rm r}}^{(2)}&=-\frac{1}{2}\left[iS^{(1)},T-H_{{\rm
r}}^{(1)}\right]\\\
&=-\frac{1}{2}\left[iS^{(1)},\sum_{l}(T_{-1,l}+2T_{0,l}+T_{1,l})e^{-il\omega\tau}\right]\end{split}$
(29)

The relevant term in the commutator causing the recombination is given by

$\begin{split}&\frac{1}{2}\left[\sum_{l}\frac{T_{-1,l}}{U+l\omega}e^{il\omega\tau},2\sum_{l^{\prime}}T_{0,l^{\prime}}e^{il^{\prime}\omega\tau}\right]\\\
&=\sum_{l,l^{\prime}}\frac{e^{i(l+l^{\prime})\omega\tau}}{U+l\omega}\left[T_{-1,l},T_{0,l^{\prime}}\right]\end{split}$
(30)

The effect of this term in the lowest order of the high-frequency expansion is
given by its time average over one time period,

$\begin{split}H_{{\rm
rc}}^{(F)}&=\frac{1}{T}\int_{0}^{T}d\tau\sum_{l,l^{\prime}}\frac{e^{i(l+l^{\prime})\omega\tau}}{U+l\omega}\left[T_{-1,l},T_{0,l^{\prime}}\right]\\\
&=\sum_{l}\frac{1}{U+l\omega}\left[T_{-1,l},T_{0,-l}\right]\end{split}$ (31)

We can explicitly write down $T_{0,l},T_{-1,l},T_{1,l}$ using $g_{ij\sigma}$,
$h_{ij\sigma}$ defined in Appendix A as

$\begin{split}T_{0,l}&=\sum_{\langle
i,j\rangle,\sigma}t_{ij}\mathcal{J}_{l}(\mathbf{r}_{ij}.\mathbf{A})g_{ij\sigma}\\\
T_{1,l}&=\sum_{\langle
i,j\rangle,\sigma}t_{ij}\mathcal{J}_{l}(\mathbf{r}_{ij}.\mathbf{A})h_{ij\sigma}^{\dagger}\\\
T_{-1,l}&=\sum_{\langle
i,j\rangle,\sigma}t_{ij}\mathcal{J}_{-l}(\mathbf{r}_{ij}.\mathbf{A})h_{ij\sigma}\end{split}$
(32)

Then by using Eq. (31), we find that the recombination Hamiltonian is given
by,

$\begin{split}H_{{\rm
rc}}^{(F)}&=\sum_{ijk}\sum_{l}\frac{1}{U-l\omega}\mathcal{J}_{-l}(\mathbf{r}_{ji}.\mathbf{A})\mathcal{J}_{l}(\mathbf{r}_{jk}.\mathbf{A})\left[h_{ij},g_{jk}\right]\\\
&=\frac{1}{2}\sum_{l}\sum_{(ijk)ss^{\prime}}\frac{t_{ij}t_{jk}\mathcal{J}_{-l}(\boldsymbol{r}_{ji}\cdot\boldsymbol{A})\mathcal{J}_{l}(\boldsymbol{r}_{jk}\cdot\boldsymbol{A})}{U-l\omega}\\\
&\times(1-n_{i,\bar{s}})c_{is}n_{ks^{\prime}}c_{k\bar{s}^{\prime}}\,\vec{\sigma}_{s\bar{s}^{\prime}}\cdot\boldsymbol{S}_{j}\end{split}$
(33)

which is Eq. (V) in Sec. V.