\section{Introduction}

\begin{figure*}[t]
\begin{centering}
\includegraphics[width=.8\linewidth]{overview.pdf}
\par\end{centering}
\caption{\textbf{Various Floquet-induced quantum states.}
For Dirac fermions as in graphene,
\textbf{a}
illumination of a circularly-polarised light (CPL)
induces the Floquet topological insulator,
which arises from photon-assisted complex
hopping (\textbf{b}), where colour shadings
represent the piece-wise magnetic fluxes.
\textbf{c} For magnetic phases in strongly-correlated
systems,
illumination of CPL induces
a chiral spin coupling, for which \textbf{d}
illustrates a resonant behaviour at
$\omega/U = 1, 1/2, 1/3, \cdots$.
\textbf{e} For superconductors (present work), illumination of CPL
on a $d_{x^2-y^2}$-wave superconductor produces
pairing amplitudes $\langle\hat{c}_{i\uparrow}\hat{c}_{j\downarrow}\rangle$ across nearest neighbours (red: positive;
blue: negative) along with
imaginary diagonal bonds (magenta and green),
leading to an emergent complexified
gap function $F(\bm{k})$,
and we have a photo-induced chiral $d_{x^2-y^2}+id_{xy}$ superconductivity (\textbf{f}), where arrows schematically indicate the phase of the complex gap function.}
\label{fig:overview}
\end{figure*}

It has become one of the key pursuits in the condensed-matter physics to design new quantum phases
as exemplified by unconventional superconductivity
and topological states.
Conventionally, materials
design is the accepted way, in which
we tailor the crystal structures and consituent
elements, as combined, if necessary, with
carrier doping, pressure applications, etc.
An entirely different avenue should be a
``non-equilibrium design", in which we envisage
to realise interesting quantum phases by putting
the systems out of equilibrium, typically by
illuminating intense laser lights.
This opens an in-situ
way to convert the system that would be unthinkable
in equilibrium, and is also of
fundamental interests as non-equilibrium physics.
Ilya Prigogine, in his book \textit{From being to becoming}~\cite{Prigogine1980},
once said that he would have liked to entitle the
book as \textit{Time, the forgotten dimension},
and we can indeed enlarge our horizon if we do not
forget the temporal dimension.

One of the most important pathways is the Floquet physics,
with which
we can seek various novel quantum states arising
from application of AC modulations to the system.
This is based on Floquet's theorem for time-periodic modulations, put forward by Gaston Floquet in 1883, which
is much older than 1928 theorem by Bloch for spatially-periodic
modulations.  A prime application of Floquet physics is
the
``Floquet topological insulator" proposed by Takashi
Oka and one of the present authors in 2009~\cite{Oka2009,Oka2010,Oka2011}.
Namely, by applying a circularly-polarised light to
honeycomb systems such as graphene [see Fig.~\ref{fig:overview}a], we can turn the
system into a topological insulator in a dynamical manner.  In other words, we are here talking about
matter-light combined states, since the electron is converted into a superposition
of the original, one-photon dressed, two-photon dressed,
..., electrons in the Floquet picture.
The Floquet topological insulator with a topological gap
exhibits a DC Hall effect despite the modulation being
AC.
The resulting state shows a kind of quantum anomalous
Hall effect (i.e., quantum Hall effect in zero
magnetic field) originally proposed for the static case
by Duncan Haldane back
in 1988~\cite{Haldane1988}.  Indeed, the effective Hamiltonian
for the irradiated system in the leading order
in the Floquet formalism exactly coincides with
Haldane's model [Fig.~\ref{fig:overview}b], as shown by Takuya Kitagawa \textit{et al.}~\cite{Kitagawa2011}
The year 2019 witnessed an
experimental detection of the Floquet topological
insulator in graphene by James McIver \textit{et al.}~\cite{McIver2019}

Thus we have a surge of interests in
Floquet topological phases, for many-body problems
as well as one-body cases.  Floquet physics has also
been extended to explore superconductivity in
AC-modulated situations.  It is out of the scope of
the present paper to review the whole field, but
let us briefly mention that the spectrum of
the many-body interests in Floquet physics
covers a range of quantum phases as

\begin{table}[h]
\begin{tabular}{c|c}
topology & superconductivity \\
\hline
FTI $\leftrightarrow$ Mott's insulator & Attraction-repulsion conversion \\
chiral spin states & $\eta$-pairing
\end{tabular}
\end{table}

\noindent First, when a repulsive Hubbard model
for correlated electrons is illuminated by laser,
phase transitions arise between
the Floquet topological insulator (FTI) and Mott's
insulator on a phase diagram
against the AC modulation intensity and the
repulsive electron-electron interaction.\cite{Mikami2016}
Second, circularly-polarised
laser lights can induce chiral spin interactions,
$(\hat{\bm{S}}_{i}\times\hat{\bm{S}}_{j})\cdot\hat{\bm{S}}_{k}$, and
consequent spin liquids for
the strongly-coupled systems with the repulsion $U$ much greater than the electron hopping $t_0$
[see Fig.~\ref{fig:overview}c.\ref{fig:overview}d]~\cite{Takayoshi2014,Takayoshi2014-2,Kitamura2017,Claassen2017}.
There, we can realise that the Floquet-incuded chiral coupling becomes significant
when the frequency $\omega$ of
the laser is set close to the Hubbard $U$,
with vastly different behaviours
between the cases when $\omega$
is slightly red-detuned from $U$ and
blue-detuned.

Third, if we turn to superconductivity, Floquet physics
can be used to even convert repulsive interactions into
attractive ones by applying (linearly-polarised) laser.\cite{Tsuji2011}
Fourth, we can turn the
usual pairing into an exotic $\eta$-pairing
(condensation of pairs at Brillouin
zone corners)
by illuminating a superconductor with a
(linearly-polarised) laser.\cite{Kitamura2016}
There, the Floquet effective
pair-hopping is drastically modulated
when $\omega$ is close to $U$, again
with different behaviours between $\omega$
red-detuned or blue-detuned from $U$.

Now the purpose of the present paper is to
encompass topological and superconducting
properties to seek whether a ``Floquet-induced
topological superconductivity" can arise.
While there have been various attempts at realising Floquet topological superconducting states~\cite{Ezawa2014,Zhang2015,Takasan2017,Chono2020,Kumar2021,Dehghani2021},
an obstacle is the fact that the pairing symmetry
is hard to be controlled in a direct manner, since the gap function does not couple to electromagnetic fields.
Thus, usually, one has to modulate the one-body part in a nontrivial manner (with the gap function kept intact), which e.g. necessitates the Rashba spin-orbit coupling.
Instead, the strategy of the present paper
is to exploit the peculiar \textit{laser-induced
interactions emergent in strongly-correlated systems}.
Namely, we illuminate a circularly-polarised light (CPL)
to the repulsive Hubbard model in the strong-coupling regime to
modulate the pairing interaction and the resultant gap function.
We shall show that a $d_{x^2-y^2}$-wave superconductor is indeed changed into a topological $d_{x^2-y^2}+id_{xy}$ wave [Fig.~\ref{fig:overview}f].
This will be shown in the Floquet formalism for the Gutzwiller-projected
effective Hamiltonian with the time-periodic Schrieffer-Wolff transformation.
The $d+id$ pairing is shown to arise from the
chiral spin coupling along with the three-site term
caused by the CPL [Fig.~\ref{fig:overview}e].  The latter
term turns out to remain significant
even for low frequencies and low intensities of the CPL.
This will be reflected in a phase diagram against the laser intensity
and temperature obtained here for various frequencies red-detuned from
the Hubbard $U$, along with transient dynamics.

\section{Results}

\subsection{Low-energy effective Hamiltonian}

We take in the present study the hole-doped Mott insulator on a square lattice, which is irradiated by an intense circularly-polarised light (CPL).
This can be minimally modelled by the repulsive
Hubbard model, with a Hamiltonian
\begin{align}
\hat{H}(t)&=-\sum_{ij\sigma}t_{ij}e^{-i\bm{A}(t)\cdot\bm{R}_{ij}}\hat{c}_{i\sigma}^\dagger\hat{c}_{j\sigma}+\frac{U}{2}\sum_i \hat{n}_{i}(\hat{n}_{i}-1),\label{eq:Hubbard}
\end{align}
where $\hat{c}_{j\sigma}$ annihilates an electron on the site $j$ at $\bm{R}_j$ with spin $\sigma=\uparrow,\downarrow$, and $\hat{n}_i=\sum_{\sigma}\hat{n}_{i\sigma}=\sum_{\sigma}\hat{c}_{i\sigma}^\dagger\hat{c}_{i\sigma}$ is the density operator.
For the hopping amplitude $t_{ij}$, here
we consider the second-neighbour ($t_0^\prime$) as well
as the nearest-neighbour hopping ($t_0$), as
necessitated for realising a CPL-induced
topological superconductivity
[as we shall see below Eq.~(\ref{eq:delta-f})].
The onsite repulsion is denoted as $U\;(>0)$.

The circularly-polarised laser field is introduced via the Peierls phase, which involves
$\bm{R}_{ij}\equiv\bm{R}_{i}-\bm{R}_{j}$,
and the vector potential,
\begin{align}
\bm{A}(t)&=\frac{\bm{E}}{2i\omega}e^{-i\omega t}-\frac{\bm{E}^\ast}{2i\omega}e^{i\omega t},\\
\bm{E}&=E(1,i),
\end{align}
for the right-circularly polarised case; replace $\bm{E}$ with $\bm{E}^\ast$ for the left circulation.   We set $\hbar=e=1$ hereafter.

Since the laser electric field
is time-periodic, $\hat{H}(t+2\pi/\omega)=\hat{H}(t)$,
we can employ the Floquet formalism~\cite{Oka2019}, where the eigenvalue, called quasienergy, of the discrete time translation plays a role of energy.
Because the quasienergy spectrum is invariant under the time-periodic unitary transformation, $e^{-i\hat{\Lambda}(t)}$ with $\hat{\Lambda}(t+2\pi/\omega)=\hat{\Lambda}(t)$,
we can introduce an effective static Hamiltonian $\hat{H}_{\text{F}}$ as
\begin{align}
\hat{H}_{\text{F}} = e^{i\hat{\Lambda}(t)}[\hat{H}(t)-i\partial_t]e^{-i\hat{\Lambda}(t)},\label{eq:canonical}
\end{align}
where $\hat{\Lambda}(t)$ is determined such that $\hat{H}_{\text{F}}$ becomes time-independent.
With such a transformation, we can analyse the time-dependent original problem using various approaches that are designed for static Hamiltonians.

While it is difficult to determine $\hat{\Lambda}(t)$ and $\hat{H}_{\text{F}}$
in an exact manner,
we can obtain their perturbative expansions in various situations.
The high-frequency expansion~\cite{Bukov2015,Eckardt2015,Mikami2016} is a seminal example of such perturbative methods, where the effective Hamiltonian $\hat{H}_{\text{F}}$ is obtained as an asymptotic series in $1/\omega$, and widely used for describing the Floquet topological phase transition.
However, since the onsite interaction $U$ in Mott insulators is typically much greater than the photon energy $\omega$,
the high-frequency expansion is expected to be invalidated for the present case.

Still, in the present system in a strongly-correlated regime $t_0\ll U$, we can employ the Floquet extension of the strong-coupling expansion (Schrieffer-Wolff transformation)~\cite{Mentink2014,Bukov2016,Kitamura2017,Claassen2017,Kumar2021}, when the photon energy $\omega$ exceeds the band width $\sim t_0$ of the doped holes (and chosen to be off-resonant with the Mott gap $\sim U$), i.e., $t_0\ll\omega< U$. As we describe the detailed derivation in Methods, we can obtain $\hat{\Lambda}(t)$ and $\hat{H}_{\text{F}}$ in perturbative series in the hopping amplitudes $t_0,t_0^\prime$.
While the $t_0\ll U$ assumption yields a low-energy effective model known in the undriven case as a $t$-$J$ model~\cite{Ogata2008} which describes the dynamics of holes in the background of localised spins, an essential difference in the Floquet states under CPL is that extra terms emerge.
Namely, the Hamiltonian of the irradiated case reads
\begin{align}
\hat{H}_{\text{F}} & =-\sum_{ij\sigma}\tilde{t}_{ij}\hat{P}_{G}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}\hat{P}_{G}
+\frac{1}{2}\sum_{ij}J_{ij}\left(\hat{\bm{S}}_{i}\cdot\hat{\bm{S}}_{j}-\frac{1}{4}\hat{n}_{i}\hat{n}_{j}\right)
\nonumber \\
 & +\sum_{ijk\sigma\sigma^{\prime}}\Gamma_{i,j;\,k}\hat{P}_{G}\left[(\hat{c}_{i\sigma}^{\dagger}\bm{\sigma}_{\sigma\sigma^{\prime}}\hat{c}_{j\sigma^{\prime}})\cdot\hat{\bm{S}}_{k}-\frac{1}{2}\delta_{\sigma\sigma^{\prime}}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}\hat{n}_{k}\right]\hat{P}_{G}
 \nonumber \\
& +\frac{1}{6}\sum_{ijk}J^\chi_{ijk}(\hat{\bm{S}}_{i}\times\hat{\bm{S}}_{j})\cdot\hat{\bm{S}}_{k}
 ,\label{eq:tjjchi}
\end{align}
where $\hat{P}_G=\prod_i(1-\hat{n}_{i\uparrow}\hat{n}_{i\downarrow})$ projects out doubly-occupied configurations. Here, we have retained all the processes up to the second order (in the hopping amplitudes $t_0,t_0^\prime$), and additionally taken account of the fourth-order processes for the last line of the above equation.  Thus
the photon-modified exchange interaction, $J_{ij}$,
and the photon-induced correlated hopping, $\Gamma_{i,j;\,k}$,
are of second order, while the photon-generated
chiral spin coupling,
$J^{\chi}_{ijk}$, is of fourth order.

Let us look at the effective model in the laser field term by term.
First, the hopping amplitude of the holes is renormalised from $t_{ij}$ into $\tilde{t}_{ij}=t_{ij}\mathcal{J}_0(A_{ij})$
due to the time-averaged Peierls phase~\cite{Dunlap1986,Eckardt2009}, where $A_{ij}=E|\bm{R}_{ij}|/\omega$ and $\mathcal{J}_m$ is the $m$-th Bessel function.
Then we come to the interactions.  The spin-spin interaction
has a coupling strength dramatically affected
by the intense electric field, since it is mediated by kinetic motion of electrons. Namely, the static kinetic-exchange interaction $J$ is modulated as~\cite{Mentink2014}
\begin{align}
J_{ij}= \sum _{m=-\infty}^\infty \frac{4t_{ij}^2\mathcal{J}_m(A_{ij})^2}{U-m \omega},\label{eq:heisenberg}
\end{align}
which is a sum over $m$-photon processes, and
appears in the second term on the first line of Eq.~(\ref{eq:tjjchi}).

While these modulations of $t$ and $J$ also occur for the case of linearly-polarised lasers (i.e., for time-reversal symmetric modulations),
a notable feature of the CPL is the emergence of the \textit{time-reversal breaking many-body interactions}
specific to the strong correlation.
In the one-body Floquet physics, the two-step hopping (i.e., a perturbative process with hopping twice) in CPL
becomes imaginary (and thus breaks the time-reversal symmetry)
already on the noninteracting level if we take a honeycomb lattice~\cite{Oka2009,Kitagawa2011,Mikami2016},
while such an imaginary hopping does not arise for the (noninteracting) square lattice due to a cancellation of contributions from different paths as depicted in Fig.~\ref{fig:two-step}a, which is why a honeycomb
lattice fits with the FTI.

\begin{figure}[t]
\begin{centering}
\includegraphics[width=0.9\linewidth]{two-step.pdf}
\par\end{centering}
\caption{\label{fig:two-step}\textbf{Schematic picture for two-step hopping.}
\textbf{a} Photon-induced hopping processes in the noninteracting square lattice, where contributions from two different paths (red and blue dotted arrows) cancel with each other. \textbf{b} For a strongly-correlated model,
the photon-induced hopping processes are sensitive to the charge and spin configurations on the path, which lifts the cancellation in general.}
\end{figure}

If we move on to correlated systems,
on the other hand, the hopping process involves interactions for spinful
electrons on the path,
which lifts such cancellations, as exemplified in Fig.~\ref{fig:two-step}b.
This implies that such a time-reversal breaking term is present
in (interacting) square lattice as well.
The two-step \textit{correlated} hopping $\Gamma$ appears in the second-order perturbation in the presence of holes [the second line of Eq.~(\ref{eq:tjjchi})], which
involves three sites and sometimes referred to as the ``three-site term"~\cite{Ogata2008}.
While the three-site term also appears in the undriven case with the amplitude $\Gamma_{i,j;\,k}=J_{ij}/4$,
this term involving the dynamics of holes (i.e., the change of holes' positions)
is usually neglected because its contribution is small in the low-doping regime [see the Gutzwiller factor in the next subsection, $\sim g^2\delta$ in Eqs.~(\ref{eq:dispersion}), (\ref{eq:gapfunction})].
In sharp contrast, when the system is driven by a CPL, $\Gamma_{i,j;\,k}$ acquires an important
\textit{three-site imaginary part} as
\begin{align}
\text{Im}\,\Gamma_{i,j;\,k} & =\sum_{m=1}^\infty\frac{2t_{ik}t_{kj}\mathcal{J}_{m}(A_{ik})\mathcal{J}_{m}(A_{kj})}{m\omega(1-m^2\omega^2/U^2)}\sin m(\theta_{jk}-\theta_{ik})\label{eq:imgamma}
\end{align}
with $\theta_{ij}$ defined as $\bm{R}_{ij}=|\bm{R}_{ij}|(\cos\theta_{ij},\sin\theta_{ij})$.
The term triggers a topological superconductivity even
when it is small, as we shall see.

In addition to the two-step correlated hopping $\Gamma$, another
time-reversal breaking term arises
if we go over to higher-order perturbations.
In the previous studies of the half-filled case~\cite{Kitamura2017,Claassen2017},
an emergent three-spin interaction (i.e., the scalar \textit{spin chirality term} $J^\chi$) is shown
to appear in the fourth-order perturbation.
This appears on the last line
in Eq.~(\ref{eq:tjjchi}) with a coefficient $J^\chi_{ijk}$.
While this term is basically much smaller than the second-order terms,
it may become comparable with the contribution from $\text{Im}\,\Gamma$, as the spin chirality term does not accompany the dynamics of holes
and has a larger Gutzwiller factor [$\sim g^3$ in Eqs.~(\ref{eq:dispersion}), (\ref{eq:gapfunction})].

\subsection{Mean-field decomposition under Gutzwiller ansatz}

\begin{figure*}[t]
\begin{centering}
\includegraphics[width=0.8\linewidth]{params.pdf}
\par\end{centering}
\caption{\label{fig:params}\textbf{Dependence of the coupling constants on the laser intensity for the driven model.}
\textbf{a}, \textbf{b} Two-step correlated hopping $\gamma$ (\textbf{a}), and scalar spin-chirality term $J_\chi$ (\textbf{b}) against circularly-polarised-light amplitude $E$ and driving frequency $\omega$.
Fig.~\ref{fig:overview}d corresponds to a cross section of \textbf{b} at $E=6t_0/a$.
\textbf{c}-\textbf{e} The renormalised hopping amplitude $\tilde{t}$, exchange coupling $J$, along with
the two-step correlated hopping $\gamma$ and the chiral spin-coupling $J_\chi$, against $E$. The driving frequency is chosen as
$\omega=1.7t_0=0.14U$ (\textbf{c}), $\omega=5.4t_0=0.45U$ (\textbf{d}), or $\omega=10.5t_0=0.88U$ (\textbf{e}) with
$t_0$: bare hopping amplitude, $a$: lattice constant.
Note a difference in the horizontal scale between \textbf{c} and \textbf{d,e}.
We take here the bare second-neighbour hopping $t_0^\prime=-0.2t_0$, the onsite interaction $U=12t_0$,
and a doping level $\delta=0.2$.
Here, the infinite summations of Fourier components [$\sum_m$ in Eq.~(\ref{eq:imgamma}) and $\sum_{lmn}$ in Eq.~(\ref{eq:chirality})] are truncated such that the Taylor series up to $E^{20}$ is reproduced. This is accurate enough for $\omega> t_0$, while for the low-frequency regime $\omega< t_0$ an intricate resonant structure due to the
$\omega=U/$integer resonance (seen also for $t_0<\omega<4t_0$) is not captured.
}
\end{figure*}

Let us investigate the fate of the $d$-wave superconductivity in the presence of the time-reversal breaking terms.
To this end, we here adopt the Gutzwiller ansatz~\cite{Ogata2003,Ogata2008} for the ground-state wavefunction and perform a mean-field analysis, as elaborated in Methods.
After the approximate evaluation of the Gutzwiller projection, the problem of determining the ground state is turned into an energy minimisation of Eq.~(\ref{eq:tjjchi}) without the Gutzwiller factor $\hat{P}_G$ if we renormalise the parameters.
The mean-field solution is then formulated as a diagonalisation of the the Bogoliubov-de Gennes Hamiltonian in the momentum space,
\begin{align}
\hat{H}_{\text{F}} & =\sum_{\bm{k}}\begin{pmatrix}\hat{c}_{\bm{k}\uparrow}\\
\hat{c}_{-\bm{k}\downarrow}^{\dagger}
\end{pmatrix}^{\dagger}\begin{pmatrix}\varepsilon(\bm{k}) & F(\bm{k})\\
F(\bm{k})^{\ast} & -\varepsilon(-\bm{k})
\end{pmatrix}\begin{pmatrix}\hat{c}_{\bm{k}\uparrow}\\
\hat{c}_{-\bm{k}\downarrow}^{\dagger}
\end{pmatrix}\label{eq:BdG}\\
 & =\sum_{\bm{k}}\begin{pmatrix}\hat{c}_{\bm{k}\uparrow}\\
\hat{c}_{-\bm{k}\downarrow}^{\dagger}
\end{pmatrix}^{\dagger}\left[\sum_{\tau}\begin{pmatrix}\varepsilon_\tau & F_\tau\\
F_\tau^{\ast} & -\varepsilon_\tau
\end{pmatrix}\cos\bm{k}\cdot\bm{R}_{i,i+\tau}\right]
\begin{pmatrix}\hat{c}_{\bm{k}\uparrow}\\
\hat{c}_{-\bm{k}\downarrow}^{\dagger}
\end{pmatrix}.
\end{align}
Here we have
\begin{align}
\varepsilon_\tau & =-g\delta\tilde{t}_{i,i+\tau}-g^2\frac{3-\delta^{2}}{8}J_{i,i+\tau}\,\chi_{\tau}
\nonumber\\
& -g^2\delta\sum_{\tau^{\prime}}\text{Re}\left(\frac{1-\delta^{2}}{4}\Gamma_{i,i+\tau;\,i+\tau^{\prime}}+\frac{3-\delta}{2}\Gamma_{i,i+\tau+\tau^{\prime};\,i+\tau}\,\chi_{\tau^{\prime}}\right) \nonumber \\
&-\frac{3}{16}g^3\sum_{\tau^{\prime}}J^{\chi}_{i,i+\tau,i+\tau+\tau^{\prime}}\text{Im}(\Delta_{\tau+\tau^{\prime}}^{\ast}\Delta_{\tau^{\prime}}),\label{eq:dispersion}
\\
F_\tau  &=g^{2}\frac{3+\delta^{2}}{8}J_{i,i+\tau}\Delta_{\tau}+g^{2}\delta\frac{3+\delta}{2}\sum_{\tau^{\prime}}\text{Re}\Gamma_{i,i+\tau+\tau^{\prime};\,i+\tau}\Delta_{\tau^{\prime}}\nonumber\\
 & +ig^{2}\sum_{\tau^{\prime}}\left(\delta\frac{3+\delta}{2}\text{Im}\Gamma_{i,i+\tau+\tau^{\prime};\,i+\tau}+\frac{3}{8}gJ_{\chi,i,i+\tau,i+\tau+\tau^{\prime}}\,\chi_{\tau+\tau^{\prime}}\right)\Delta_{\tau^{\prime}},
\label{eq:gapfunction}
\end{align}
where the expressions are independent of the site index $i$ due to the spatial translational symmetry, and we have denoted the doping level as
\begin{equation}
\delta=1-\frac{1}{N}\sum_i\langle\hat{n}_i\rangle
\end{equation}
with $N$ being the number of lattice sites.
Here $\tau$ runs over $\tau=mx+ny$ with $m,n\in\mathbb{Z}$, and the label $i+\tau$ represents the site at $\bm{R}_{i+\tau}=\bm{R}_i+(m,n)$ in units of the lattice constant.
The dependence on $\delta$ and a factor $g=2/(1+\delta)$ appears as a result of the Gutzwiller projection, which suppresses the contribution from charge dynamics as represented by $\tilde{t}$ and $\Gamma$ in the small $\delta$ regime.
For the detailed derivation, see Methods.
Throughout the present study, we take $\delta=0.2$.
The mean-field Hamiltonian is self-consistently determined
for the bond order parameter $\chi$ and
the pairing amplitude $\Delta$ as
\begin{align}
\chi_{\tau}&=\frac{1}{N}\sum_i\langle\hat{c}_{i\uparrow}^{\dagger}\hat{c}_{i+\tau\uparrow}+\hat{c}_{i\downarrow}^{\dagger}\hat{c}_{i+\tau\downarrow}\rangle,\\
\Delta_{\tau}&=\frac{1}{N}\sum_i\langle\hat{c}_{i\uparrow}\hat{c}_{i+\tau\downarrow}-\hat{c}_{i\downarrow}\hat{c}_{i+\tau\uparrow}\rangle.
\end{align}

In the absence of the external field, the system undergoes a phase transition from normal to the $d_{x^2-y^2}$-wave superconductivity when the temperature is sufficiently low~\cite{Ogata2008}. The corresponding order parameter is
\begin{equation}
\Delta_{\pm x}=-\Delta_{\pm y}\equiv \Delta,
\end{equation}
for which the gap function $F(\bm{k})$ becomes
\begin{equation}
F^{x^2-y^2}(\bm{k})=\frac{3}{4}g^2J\Delta(\cos k_x-\cos k_y),
\end{equation}
as derived from the first term in Eq.~(\ref{eq:gapfunction}) with $\tau=\pm x,\pm y$. Here we have dropped a small correction due to $\delta$.
Without a loss of generality, we assume that $\Delta$ is real.
The $d_{x^2-y^2}$ gap function has nodal lines along $k_y=\pm k_x$, across which the gap function changes sign (See Figs.~\ref{fig:gap}a-\ref{fig:gap}c below).

Now let us look into the possibility for the laser field
converting this ground state into \textit{topological superconductivity} from the time-reversal breaking terms [the second line in Eq.~(\ref{eq:gapfunction})].
The leading term should be those for $\tau^\prime=\pm x,\pm y$ (with $\Delta_{\tau^\prime}=\pm\Delta$), which is nonzero even for the original $d_{x^2-y^2}$-wave ansatz and results in an imaginary gap function $\propto i\Delta$.
In particular, gathering the terms with $\tau=\pm(x+y),\pm(x-y)$, we obtain the leading modulation to the gap function as
\begin{align}
F^{xy}(\bm{k}) \simeq 3i
g^{2}\left[4\delta\gamma+g(J_{\chi}\,\chi_{x}+J_{\chi}^{\prime}\,\chi_{2x+y})\right]\Delta\sin k_{x}\sin k_{y}\label{eq:delta-f}
\end{align}
where we have defined
\begin{gather}
\gamma = \text{Im}\,(\Gamma_{i-x,i;\,i+y}
-\Gamma_{i-x-y,i+x;\,i}),\\
J_{\chi} = J_{i,i+y,i+x}^{\chi}\;, \; \;
J_{\chi}^{\prime} = J_{i-x,i,i+x+y}^{\chi}.
\end{gather}
Thus we are indeed led to
an emergence of a topological (chiral) $d_{x^2-y^2}+id_{xy}$ superconductivity with
$F(\bm{k})\simeq F^{x^2-y^2} +iF^{xy}$
with a full gap.  For the full expression of the modulated gap function, see Eq.~(\ref{eq:gap-full}) in Methods.
The key interactions, $\gamma$ and $J_\chi$, take nonzero values
when the original Hubbard Hamiltonian Eq.~(\ref{eq:Hubbard}) has the second-neighbour hopping $t_0^\prime$:
Since the two-step correlated hopping $\Gamma_{i,j;\,k}$ is composed of two hoppings,  $k\to i$ and $j\to k$,
we can see that the imaginary coefficient above has
$\gamma\propto t_0t_0^\prime$.
The chiral spin-coupling $J_\chi\propto (t_0t_0^\prime)^2$ also necessitates the next-nearest-neighbour hopping,
which can be deduced from Eq.~(\ref{eq:chirality-details}) in Methods.

Having revealed the essential couplings for the chiral superconductivity,
let us now explore how we can
\textit{optimise} the field amplitude and frequency for realising larger topological gaps.  For this,
we set here $t_0^\prime=-0.2t_0$ and $U=12t_0$
having cuprates with $t_0\simeq0.4$~eV in mind as a typical example.

We plot the essential $\gamma$ and $J_\chi$ against the driving amplitude $E$ and frequency $\omega$ in Figs.~\ref{fig:params}a and \ref{fig:params}b.
If we look at the overall picture,
we can see that the dynamical time-reversal breaking is strongly enhanced along some characteristic frequencies in
a \textit{resonant} fashion at $\omega/U = 1/$integer, which we can capture from the expressions for the coupling constants in the small-amplitude
regime as follows.

The two-step correlated hopping $\gamma$ is given,
in the leading order in the amplitude, as
\begin{equation}
\gamma=\frac{2t_{0}t_{0}^{\prime}\mathcal{J}_{2}\left(\frac{Ea}{\omega}\right)\mathcal{J}_{2}\left(\sqrt{2}\frac{Ea}{\omega}\right)}{\omega(1-4\omega^{2}/U^{2})}+O(E^{12}).
\end{equation}
This expression, as a function of $E$,
takes the maximal value around $E \simeq 2.45\omega/a$,
while diverges for $\omega\to0$ or $\omega\to U/2$ for each
value of $E$,
as seen from the energy denominator.
The enhancement occurs even
in the low-frequency regime, which should be advantageous
for experimental feasibility,
since the required field amplitude ($E\simeq 2.4\omega/a$)
can be small.

If we turn to the spin-chirality term, $J_\chi$ has a complicated form involving Bessel functions [see Eq.~(\ref{eq:chirality-details}) below]. If we
Taylor-expand it in $E$, we have
\begin{equation}
J_{\chi}\sim\frac{2(Ea)^4 t_{0}^{2}t_{0}^{\prime2}(2U^6+75\omega^{2}U^4-399\omega^{4}U^2-164\omega^{6})}{\omega(U^2-\omega^{2})^{3}(U^2-4\omega^{2})^{3}},
\label{eq:Jchi-fourth}
\end{equation}
which takes large values around $\omega=U$ due to the
energy denominator in the above expression, while $\gamma$ is small in this regime.
The expression reveals that a dynamical time-reversal breaking also occurs
as a fourth-order nonlinear effect with respect to the field strength $E$.

\begin{figure*}[t]
\begin{centering}
\includegraphics[width=\linewidth]{gap.pdf}
\par\end{centering}
\caption{\textbf{Gap functions.} \textbf{a},\textbf{d},\textbf{g} Absolute values of the gap function, $|F(\bm{k})|$. \textbf{b},\textbf{e},\textbf{h} Phases, $\arg F(\bm{k})$.
 The dashed lines represent the Fermi surface at the critical temperature $T=T_c$.
\textbf{c},\textbf{f},\textbf{i} Density of states against
energy in units of the bare hopping amplitude $t_0$. These are shown for \textbf{a}-\textbf{c} $d_{x^2-y^2}$-wave superconductivity in the absence of driving field, and \textbf{d}-\textbf{i} $d_{x^2-y^2}+id_{xy}$-wave superconductivity in the circularly-polarised laser.  The driving frequency $\omega$ and the field amplitude $E$ are chosen as $\omega=10.5t_0=0.88U, E=16t_0/a$ (\textbf{d}-\textbf{f}), or $\omega=1.7t_0=0.14U, E=3.5t_0/a$ (\textbf{g}-\textbf{i}).
We take a parameter set: the bare second-neighbour hopping $t_0^\prime=-0.2t_0$, the onsite interaction $U=12t_0$,
and a doping level $\delta=0.2$.}
\label{fig:gap}
\end{figure*}
\begin{figure*}[t]
\begin{centering}
\includegraphics[width=\linewidth]{phase.pdf}
\par\end{centering}
\caption{\label{fig:phase}\textbf{Finite-temperature phase diagram.}
The superconducting order parameter $|\Delta|$ (\textbf{a}, \textbf{c}), and the time-reversal breaking order parameter $\Delta^\prime \equiv -i\Delta_{\pm (x+y)}=i\Delta_{\pm (x-y)}$
 (\textbf{b}, \textbf{d}) plotted against the field amplitude $E$ and temperature $T$.
The driving frequency is set to
 $\omega=10.5t_0=0.88U$ (\textbf{a}, \textbf{b}) or $\omega=1.7t_0=0.14U$ (\textbf{c}, \textbf{d}). Dotted lines represent the effective temperature of the Floquet Hamiltonian when the system is quenched from $E=0$ at a
temperature that is the leftmost starting point
of each dotted line. We take the same
bare second-neighbour hopping $t_0^\prime=-0.2t_0$, onsite interaction $U=12t_0$,
and doping level $\delta=0.2$ as in Fig.4.}
\end{figure*}

We can see in Figs.~\ref{fig:params}c-e
how $J_{\chi}$, as well as the exchange interaction
$J$ and
the renormalised nearest-neighbour hopping $\tilde{t}$,
vary with the CPL
intensity for several representative frequencies
for which the dynamical time-reversal breaking becomes prominent.
The hopping amplitude $\tilde{t}$
vanishes around the peak of $\gamma$ ($E \simeq 2.45\omega/a$), as seen in Figs.~\ref{fig:params}c,d.
This involves a zero of the Bessel
function (at $E\sim2.40\omega/a$), and known as the dynamical localisation~\cite{Dunlap1986}.
In Figs.~\ref{fig:params}d,e, we can see a strong enhancement of $J$.  Panel~d has $\omega = 0.45U$ (slightly red-detuned from $U/2$), while panel~e
has $\omega = 0.88U$ (slightly red-detuned from $U$).
If we go back to
Eq.~(\ref{eq:heisenberg}), the former has to do with
the energy denominator for $m=2$, while the
latter for $m=1$. While the driving frequency is set
to be red-detuned from the resonance (which we can call
``$U-m\omega$ resonance") in these results, the blue-detuned cases would give negative contributions from these terms,
which will lead to a ferromagnetic exchange interaction, unfavouring the spin-singlet $d$-wave.
With an optimal choice of the driving field, $\gamma$
attains a significantly large value $\simeq0.04 t_0$,
which yields $F^{xy}(\bm{k})\simeq 0.3t_0\times i\Delta\sin k_{x}\sin k_{y}$ (for $\delta=0.2$).
This is remarkably large and \textit{comparable} with the
undriven $F^{x^2-y^2}(\bm{k})\simeq 0.7t_0\times \Delta (\cos k_{x}-\cos k_{y})$,
even though the coupling constant itself
is much smaller than $J$.  This is the first key result
of the present work.

\subsection{Phase diagram}

Now we investigate the ground state of the effective static Bogoliubov-de Gennes Hamiltonian for several choices of the CPL parameters.
Let us show the gap function and the density of states in Fig.~\ref{fig:gap}.
In the absence of the external field in Figs.~\ref{fig:gap}a-\ref{fig:gap}c,
the gap function has the $d_{x^2-y^2}$ symmetry with nodal lines that give the zero gap at the Fermi energy.
We now switch on the CPL, with $\omega=10.5t_0=0.88U$, $E=16t_0/a$ for Figs.~\ref{fig:gap}d-\ref{fig:gap}f, or with $\omega=1.7t_0=0.14U$, $E=3.5t_0/a$ for Figs.~\ref{fig:gap}g-\ref{fig:gap}i.
The former represents the case where
$J_\chi$ is dominant (see Fig.~\ref{fig:params}e), while in the latter $\gamma$ plays the central role.
In both cases we have the $d_{x^2-y^2}+id_{xy}$ pairing, where
the nodal lines are gapped out due to the complex gap function.
We can indeed see in Figs.~\ref{fig:gap}f,i clear energy gaps in the density of states with a gap size comparable with the original $d_{x^2-y^2}$ superconducting gap (as measured by
the energy spacing between the two peaks in the
density of states).  Note that
the band width and the $d_{x^2-y^2}$ superconducting gap
are renormalised there, due to the modified $\tilde{t}$ and $J$.

To examine the robustness of the $d+id$ superconductivity, we calculate $|\Delta|$ and $\Delta^\prime \equiv -i\Delta_{\pm (x+y)}=i\Delta_{\pm (x-y)}$ at finite temperatures,
which gives a phase diagram against $T$ and field amplitude $E$ in Fig.~\ref{fig:phase}.
We can see that the critical temperature $T_c$, as delineated by the region for $\Delta\neq0$, increases significantly
as the laser intensity $E$ is increased. The enhanced $T_c$ originates from the fact that $J$ [and $\text{Re}\,\Gamma$; see
Eq.~(\ref{eq:Gamma}) below] are enhanced when the laser is applied, as we have seen in Figs.~\ref{fig:params}c-\ref{fig:params}e. The time-reversal breaking order parameter $\Delta^\prime$ emerges below $T_c$ down to $T=0$,
with the field dependence emerging from those of $\gamma$ and $J_\chi$.

\begin{figure*}[t]
\begin{centering}
\includegraphics[width=\linewidth]{quench.pdf}
\par\end{centering}
\caption{\label{fig:dynamics}\textbf{Quench dynamics of the order parameters.} Time evolution of
$|\Delta|=|\langle c_{i\uparrow} c_{i+x,\downarrow} \rangle|$ (\textbf{a}, \textbf{c}) and $|\Delta^\prime|=|\langle c_{i\uparrow} c_{i+x+y,\downarrow} \rangle|$ (\textbf{b}, \textbf{d})
after a sudden quench of the field amplitude from $E=0$ to $E=7t_0/a$ with an initial temperature $k_BT_{\text{initial}}=0$ (\textbf{a}, \textbf{b}), or to $E=16t_0/a$ with an initial
$k_BT_{\text{initial}}=0.2t_0$ (\textbf{c}, \textbf{d}).
The latter corresponds to the case in which
the superconductivity appears even when
we start from an initial temperature above $T_c$.
Horizontal lines indicate the equilibrium value at the effective temperature.
The driving frequency is set to $\omega=10.5t_0=0.88U$ with $t_0$ being the bare hopping amplitude.
}
\end{figure*}

\subsection{Transient dynamics}
So far, we have investigated static properties of the effective Hamiltonian $\hat{H}_{\text{F}}$. While the results clearly show the presence of Floquet topological superconductivity in a wide parameter region,
an important question from a dynamical viewpoint is whether we can achieve the Floquet topological superconductivity within
a short enough time scale, over which the description with the effective Hamiltonian Eq.~(\ref{eq:tjjchi}) is valid (i.e., over which we can neglect
the coupling to phonons, the heating process due to higher-order perturbations, etc).

Thus let us study the dynamics of the superconducting gap, by solving a quench problem formulated as follows.
We first prepare the initial state as the mean-field ground state of the equilibrium $t$-$J$ model [i.e., Eq.~(\ref{eq:tjjchi}) with $E=0$],
and then look at the transient dynamics when the
CPL electric field is
suddenly switched on, by changing
the Hamiltonian at $t=0$ to the effective Hamiltonian $\hat{H}_{\text{F}}$ (\ref{eq:tjjchi}) with $E\neq0$.
In such a treatment, $\langle\hat{H}_{\text{F}}\rangle$ is a conserved quantity for $t>0$, and
the driven state is expected to thermalise to the equilibrium state of $\hat{H}_{\text{F}}$ having a temperature that
corresponds to the internal energy $\langle\hat{H}_{\text{F}}\rangle$ (which we call an effective temperature).
Note that here we implicitly neglect the fact that $\hat{H}_{\text{F}}$s for $E=0$ and $E\neq0$ are on different frames [specified by $\hat{\Lambda}(t)$].

Before exploring the dynamics, it is important to check
whether the expected steady state of the quench problem remains superconducting.  We can do this by looking at
the above-defined effective temperature on the phase diagram as indicated by dotted lines in Fig.~\ref{fig:phase}.
While the effective temperature rises as we increase the field strength $E$ when we start from temperatures below $T_c$, we can see that
the effective temperature does \textit{not} exceed $T_c$ in a wide parameter region, which implies that the Floquet topological superconductivity should indeed appear as a steady state of the quench dynamics.
A sudden increase of temperature around $E\sim 4t_0/a$ in Figs.~\ref{fig:phase}c, \ref{fig:phase}d is
due to the band flipping (sign change of $\tilde{t}$), which makes the kinetic energy of the initial state quite large. As we further increase the field strength, the effective temperature eventually reaches infinity and will even become negative~\cite{Tsuji2011}.

Now, the question is whether the equilibration (known as the Floquet prethermalisation) occurs
fast enough (e.g., faster than the neglected heating processes,
within the experimentally accessible pulse duration, etc.).
We compute the time evolution of the superconducting order parameter in the time-dependent mean-field approximation (with the Gutzwiller ansatz),
and plot the time evolution of $|\Delta|$ and $|\Delta^\prime|$ in Figs.~\ref{fig:dynamics}a, \ref{fig:dynamics}b, where we set the initial temperature at $k_BT_{\text{initial}}=0$ with  $E=7t_0/a$, $\omega=10.5t_0=0.88U$.
Here, the unit of time $\hbar/t_0$ corresponds to $\sim1$~fs for $t_0\simeq0.4$~eV.
We can see that $|\Delta|$ rapidly evolves with an overshooting behaviour, and converges to a certain value with a damped oscillation.
Similar behaviour can be found in a previous study of the quench problem for the $d$-wave (but within the $d_{x^2-y^2}$ pairing) superconductor~\cite{Peronaci2015}.
The time-reversal breaking order parameter $|\Delta^\prime|$ also quickly converges to a nonzero value. This is the second key result in the present work.
We note that the obtained steady state slightly deviates from the equilibrium state with respect to $\hat{H}_{\text{F}}$, where the deviation arises because the pair-breaking scattering~\cite{Peronaci2015} is lacking in the mean-field treatment of the dynamics.

The curious oscillation in the amplitude of the gap function can be interpreted as an excitation of the \textit{Higgs modes} in superconductors~\cite{Matsunaga2014,Shimano2020,Katsumi2018},
and thus the typical time scale for the oscillation (and the emergence of $\Delta^\prime$) can be roughly estimated
as the inverse of the superconducting gap ($\sim1/2|F(\bm{{k}})|$ with an appropriate $k$-average).
This should be much faster than the time scale for heating, although its nonempirical evaluation would be difficult.
We can analyse the present quench dynamics in a linearised form~\cite{Schwarz2020} if the change in the coupling constant is small,
which reveals that the appearance of $\Delta^\prime$ is described as the $A_{2g}$ Higgs mode in the $d$-wave sector with an amplitude proportional to $E^4$.

If we have a closer look at the effective temperature in Fig.~\ref{fig:phase}, we find an intriguing
phenomenon: the superconducting state can appear even when
we start from an initial temperature that is \textit{above} the $T_c$ before laser illumination.
Namely, some dotted lines that start from $T$ above $T_c$
at $E=0$ do plunge into the superconducting
region as $E$ is increased. We show
the time evolution of the order parameter for this
``nonequilibrium-induced superconductivity" in Figs.~\ref{fig:dynamics}c, \ref{fig:dynamics}d, where we set $k_BT_{\text{initial}}=0.2t_0$ as an initial temperature, and choose $E=16t_0/a$, $\omega=10.5t_0$.
Since the homogeneous mean-field ansatz without fluctuations cannot describe
the spontaneous symmetry breaking,
we instead inspect the growth of a tiny (homogeneous) perturbation $\Delta=10^{-4}$ on the initial state.
After the quench, both of the gap functions
$\Delta$ and $\Delta'$ grow exponentially and converge
respectively to nonzero values,
although the damped oscillation is slower than the previous case, and the converged values are far below those ($|\Delta|\simeq0.18$ and $|\Delta^\prime|\simeq0.03$) expected for equilibrium with the effective temperature.

\section{Discussion}
In the present paper we reveal that Floquet-induced
interactions do indeed give rise to a
novel way for creating topological superconductivity.
Let us recapitulate advantages of the present proposal for the Floquet topological superconductivity.
In the previous study of Floquet topological superconductivity in cuprates~\cite{Takasan2017},
the topological transition is triggered by the modulation of the kinetic part of the Hamiltonian [diagonal components in Eq.~(\ref{eq:BdG})], with the pairing symmetry of the gap function remaining the same.
The nontrivial structure of the kinetic part necessitates the presence of a strong Rashba spin-orbit coupling,
which would limit the applicable class of materials.
In the present approach, by contrast, we exploit the correlation effects themselves,
with which we can directly modulate the pairing symmetry to obtain the topological superconductivity.
Our approach does not require tailored structures in the one-body part, either,
and thus a simple square-lattice Hubbard model suffices for inducing the topological transition.
The present approach is also advantageous in achieving a large topological gap.  Namely,
the topological gap we conceive is comparable with $k_BT_c$ (since the $F^{x^2-y^2}(\bm{k})$ and $F^{xy}(\bm{k})$ components have the same order of magnitude), while in the previous studies the size of the topological gap is bounded e.g. by the Rashba spin-orbit coupling and estimated to be $\sim1$~K.

While we have evaluated the increase of the effective temperature upon a sudden change in the field amplitude (which might be evaded by adiabatic ramping of the field),
there is also a many-body heating process that is dropped in the present formulation.
In isolated, nonintegrable many-body systems in general, the exact eigenstates of $\hat{H}_{\text{F}}$ in the thermodynamic limit  at long times are believed to represent a featureless, infinite-temperature state
(sometimes called the Floquet eigenstate thermalisation hypothesis~\cite{DAlessio2014,Lazarides2014,Lazarides2014-2,Bukov2016-2,Seetharam2018,Mori2018}).
Thus, the present strong-coupling expansion for the nontrivial structure should be interpreted as an asymptotic expansion (with a vanishing radius of convergence) of this featureless Hamiltonian.
The expanded Hamiltonian truncated at an optimal order
accurately describes dynamics towards a long-lived state (dubbed as Floquet prethermalisation~\cite{Abanin2017,Mori2016,Kuwahara2016,Mori2018}), while the inevitable small truncation error separately describes a slow heating process towards the infinite temperature.

Specifically, the present expansion is characterised by the energy denominators of the form $(nU-m\omega)$ with integers $n,m$, in which
the enhancement of $\gamma$ and $J_\chi$ grows as
the resonance is approached.  We have to note, however, that the
accuracy of the asymptotic expansion becomes degraded,
i.e., the time scale over which the effective Hamiltonian
remains valid shrinks (with heating becoming faster) as we come closer to the resonance with vanishing denominators.
Thus, we will have to
examine whether the time scale $t\sim10/t_0$ for the emergence of $\Delta^\prime$ (see Fig.~\ref{fig:dynamics}b)
lies
within the validity of the Floquet $t$-$J$ description,
with more sophisticated methods including direct analyses of a time-dependent problem, in future works.
While the emergence of the topological superconductivity is verified here with a numerically economical approach,
we do obtain in the present paper
the effective Hamiltonian with the crucial
time-reversal breaking by
fully taking account of the noncommutative nature of the Gutzwiller projector $\hat{P}_\text{G}$
(whereas the usual hopping terms are commutative).
It should be thus important to retain such correlation effects in performing sophisticated calculations.

As for ``nonequilibrium-induced superconductivity"
(i.e., laser illumination making
a system superconducting even when we start from $T$ above $T_c$) discussed in the present paper,
there is existing literature that reports laser-induced
phenomena~\cite{Budden2021} along with related theories~\cite{Knap2016,Babadi2017,Murakami2017,Kennes2019},
although the time scale and proposed mechanism are
quite different from the present paper.
Whether the present theory has some
possible relevance will be a future problem.

Let us turn to the required field amplitude and frequency for experimental feasibility.
Specifically, if we want to employ the enhancement of time-reversal breaking around $\omega=U/2$ or $\omega=U$,
the required field intensity is $E\sim 10t_0/a$, which
corresponds to $E\sim10^2$~MV/cm for the typical
cuprates with $t_0\simeq0.4$~eV, $a\simeq3\AA$.
One reason why the strong intensities are required
derives from
the fact that the time-reversal breaking terms evoked
here are of fourth order in $E$,
while usually the time-reversal breaking terms can be of second order~\cite{Claassen2017,Kitamura2017}.
This comes from
the cancellation similar to Fig.~\ref{fig:two-step}a
 for the present square lattice [See Eq.~(\ref{eq:cancel}) for details],
and can be evaded in e.g. honeycomb and kagome lattices.
So the application of the present mechanism
for the Floquet topological superconductivity to
a wider class of materials with various
lattice structures and/or multi-orbitals may
be an interesting strategy.
More trivially, going to low frequencies is another practical route for the enhanced $\gamma$, since the required intensity scales with the frequency $\omega$.
A possibility of chiral superconductivity triggered by the same interaction but with a different mechanism
(such as proposed for the doped spin liquid on triangular lattice~\cite{Jiang2020}) is also of interest.

We can further raise an intriguing possibility that, since the emergence of the topological superconductivity is intimately related to the excitation of Higgs modes as stressed above, the topological signature
might be enhanced by resonantly exciting the Higgs modes,
which provides another future problem.


\section{Introduction}

The nature of dark matter (DM) represents a fundamental question for the understanding of our Universe.
Observational hints at cosmological and galaxy scales such as the discrepancy between the measured rotation curves of galaxies and their theoretical predictions, the formation of large structures, and the anisotropies of the Comic Microwave Background show that DM makes up about 85\% of the total matter.

The search for DM has therefore become a priority in the scientific community where a collective effort has been made in indirect, direct, and collider searches in order to unravel its mystery. Its detection would also be a milestone in searches for Physics beyond the Standard Model (SM).
In this talk, we focus on the indirect detection using the observations made by five gamma-ray experiments towards twenty dwarf spheroidal galaxies (dSphs).
These dSphs represent one of the most promising targets for DM indirect searches due to their high DM content and their negligible astrophysical background~\cite{PhysRevD.69.123501}.
They are all located at high Galactic latitude and no sign of very high energy emission has been detected so far in the dSphs' directions.
Gamma rays have the advantage of being neutral and do not get deflected by magnetic fields. Thus, the regions of $\gamma$-ray production can be traced back from the incident direction. The observations are therefore performed based on this property where the telescopes are directly pointing to the sources.

This work represents a collective effort between three imaging atmospheric Cherenkov telescope (IACT) arrays  H.E.S.S., MAGIC, and VERITAS, the water Cherenkov array HAWC, and the space-borne telescope {\it Fermi}-LAT, which agreed on sharing their data previously published individually.
The goal of our study is to combine the individual upper limits published by each collaboration in order to optimize the statistics and increase the sensitivity to potential DM signals. The combination brings the novelty of extending the upper and lower boundaries of the energy range and the derivation of the upper limits on the DM annihilation cross-section over the widest DM particle mass range ever.
In this work, each of the five collaborations performed the analysis of their own data sets using a common DM model to optimize the data handling at different energy, angular resolutions, and sensitive energy ranges of the various instruments.
By following this procedure, we also avoid the need for sharing raw data and instrument response functions (IRFs) outside the collaborations.
As no significant excess was detected from the selected sources, nor in their combination, we derive upper limits on the DM annihilation cross-section as a function of the DM particle mass by combining the likelihood functions of all dSphs and all experiments.

\section{Experiments} \label{Sec:Experiments}

\subsection{Fermi-LAT}

The {\it Fermi}-Large Area Telescope ({\it Fermi-LAT}) is the collaboration which operates the pair conversion Large Area Telescope (LAT) carried by the Fermi satellite orbiting the Earth at an altitude of 565 km since 2008. The telescope has a wide field of view covering about 20\% of the sky and scans the whole sky every 3 hours in the energy range between $\sim 20$ MeV and 1~TeV. {\it Fermi}-LAT thus covers the lowest energy region of this study. Detailed descriptions of the detector and its performance can be found in \cite{2012ApJS..203....4A}.

\subsection{HAWC}

The High-Altitude Water Cherenkov (HAWC) detector is a high-energy $\gamma$-ray telescope located at Sierra Negra, Mexico at 4100~m altitude and consists of an array consisting of 300 water Cherenkov detectors (WCD) covering an area of 22,000~$\rm{m}^2$. The WCD are sensitive to $\gamma$-ray events of energies ranging from 300 GeV to a couple hundred TeV \cite{abeysekara}. The experiment covers a field of view of 15\% of the sky at all times.

\subsection{H.E.S.S.}

The High Energy Stereosocpic System (H.E.S.S.) experiment is an array consisting of five IACTs designed to detect brief and faint flashes of Cherenkov radiation generated by very high energy $\gamma$ rays between $\sim$30 GeV and $\sim 100$ TeV. The telescope array is located in central Namibia in the Khomas Highland region at 1,800 m above sea level~\cite{Aharonian:2006pe} at 110 km south west of Windhoek. The four small telescopes are equipped with a 12~m reflector while the central one is 28~m. The array collects the $\gamma$ rays within a field of view of 5$^\circ$.

\subsection{MAGIC}
The Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescope array consists of two telescopes of 17 m diameter reflector situated at the Roque de los Muchachos Observatory on the Canary Island of La Palma, Spain, 2,200~m above sea level. MAGIC is sensitive to very high energy $\gamma$-ray events above $\sim 50$ GeV~\cite{Sitarek} and is equipped with fast imaging cameras with a field of view of 3.5\degree.

\subsection{VERITAS}
The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is an array of four telescopes of 12~m reflector located at the Fred Lawrence Whipple Observatory in Southern Arizona. The telescope array is sensitive to a very high energetic band from $\sim$ 85~GeV up to $\sim 30$~TeV whose events are recorded within a $3.5\degree$ field of view~\cite{Park}.

\section{DM signal} \label{sec:DM}

The differential flux of $\gamma$ rays from the self-annihilation of Majorana DM particles is given by:
\begin{equation}
    \label{eq:DM_flux}
    \frac{\text{d}^2\Phi \left(\langle \sigma v \rangle, J \right)}{\text{d}E \text{d}\Omega} = \frac{1}{4\pi} \frac{\langle \sigma v \rangle}{2m^{2}_{\chi}} \sum_{f} \text{BR}_{f}\frac{\text{d}N_{f}}{\text{d}E} \times \frac{ \text{d}J}{\text{d}\Omega}.
\end{equation}

The first term contains the mass $m_{\chi}$ of the DM particles in GeV and their annihilation cross-section averaged over the velocity distribution $\langle \sigma v \rangle$ in $\rm{cm}^3\rm{s}^{-1}$. It also carries the differential spectrum $dN_f /dE$ for a given annihilation channel $f$.
Since we do not assume any specific particle physics model, each channel is treated individually where the branching ratio $\text{BR}_{f} = 100\%$.
The second term, known as the astrophysical $J$ factor, describes the amount of DM annihilations within a source or a region of the sky.

The differential $J$-factor is defined as the integral of the square of the DM density distribution $\rho_{\rm{DM}}$ along the line-of-sight (l.o.s.):
\begin{equation}
  \frac{\text{d}J}{\text{d}\Omega} = \displaystyle{\int_\mathrm{l.o.s.} \: \rho_{\textrm{DM}}^2(r(s,\theta)) \:ds },
  \label{dwarf_J}
\end{equation}
where $\rho_{\textrm{DM}}$ is assumed to be spherically symmetric for all considered dSphs and depends on the distance to the center of the source $r$. This distance can also be expressed in terms of the distance $s$ from Earth along the line of sight, and the angular distance $\theta$ with respect to the center of the source, as $r(s, \theta) = (s^2 +d^2 -2sd\cos\theta)^{1/2}$, where $d$ is the distance between the Earth and the source.
The $J$ factor computation is usually performed through a Jeans analysis based on the spherical Jeans equations~\cite{Bonnivard:2014kza, Geringer-Sameth:2014yza, binney2011galactic, Bonnivard:2015xpq}.
This technique relies on the spectroscopic data and assumes that the galaxies are in steady-state hydrodynamic equilibrium, have a spherical symmetry, and are non-rotating systems to reconstruct the galactic dynamics.
In this work, we use the $J$ factors produced by Geringer-Sameth \emph{et al.}~\cite{Geringer-Sameth:2014yza}.

\section{Joint likelihood analysis}\label{Sec:Likelihood}

\subsection{Dataset}

Twenty classical and ultrafaint dwarf spheroidal galaxies are selected for the combination. All were observed by one or more instruments and previously published by individual collaborations. Table~\ref{tab:dsphs_list} presents the list of dwarf galaxies used this project and the experiments with which they were observed.

\begin{table}[h!]
\small{
\centering{
\begin{tabular}{cccc}
\hline
\hline
Source name          & Experiments  & Distance & $\log_{10}J$           \\
 & & \scriptsize{(kpc)} &  \scriptsize{$\log_{10}(\rm{GeV}^2 \rm{cm}^{-5}\rm{sr})$}   \\
\hline
\CellTopTwo{}
Bootes I             & {\it Fermi}-LAT, HAWC, VERITAS &$66$ & $18.24^{+0.40}_{-0.37}$  \\
\CellTopTwo{}
Canes Venatici I     & {\it Fermi}-LAT   &$218$ & $17.44^{+0.37}_{-0.28}$           \\
\CellTopTwo{}
Canes Venatici II     & {\it Fermi}-LAT, HAWC & $160$& $17.65^{+0.45}_{-0.43}$    \\
\CellTopTwo{}
Carina     & {\it Fermi}-LAT, H.E.S.S. & $105$&  $17.92^{+0.19}_{-0.11}$  \\
\CellTopTwo{}
Coma Berenices     & {\it Fermi}-LAT, HAWC, H.E.S.S., MAGIC & $44$&  $19.02^{+0.37}_{-0.41}$   \\
\CellTopTwo{}
Draco     & {\it Fermi}-LAT, HAWC, MAGIC, VERITAS & $76$ &  $19.05^{+0.22}_{-0.21}$  \\
\CellTopTwo{}
Fornax    & {\it Fermi}-LAT, H.E.S.S.   & $147$ & $17.84^{+0.11}_{-0.06}$  \\
\CellTopTwo{}
Hercules   & {\it Fermi}-LAT, HAWC  & $132$ & $16.86^{+0.74}_{-0.68}$  \\
\CellTopTwo{}
Leo I     & {\it Fermi}-LAT, HAWC  & $254$ & $17.84^{+0.20}_{-0.16}$  \\
\CellTopTwo{}
Leo II     & {\it Fermi}-LAT, HAWC  & $233$ & $17.97^{+0.20}_{-0.18}$   \\
\CellTopTwo{}
Leo IV   & {\it Fermi}-LAT, HAWC  & $154$&  $16.32^{+1.06}_{-1.70}$  \\
\CellTopTwo{}
Leo T     & {\it Fermi}-LAT  &$417$& $17.11^{+0.44}_{-0.39}$    \\
\CellTopTwo{}
Leo V     & {\it Fermi}-LAT & $178$&  $16.37^{+0.94}_{-0.87}$  \\
\CellTopTwo{}
Sculptor     & {\it Fermi}-LAT, H.E.S.S.  & $86$& $18.57^{+0.07}_{-0.05}$  \\
\CellTopTwo{}
Segue I     & {\it Fermi}-LAT, HAWC, MAGIC, VERITAS   & $23$ &  $19.36^{+0.32}_{-0.35}$  \\
\CellTopTwo{}
Segue II     & {\it Fermi}-LAT   & $35$ & $16.21^{+1.06}_{-0.98}$ \\
\CellTopTwo{}
Sextans    & {\it Fermi}-LAT, HAWC &$86$ & $17.92^{+0.35}_{-0.29}$  \\
\CellTopTwo{}
Ursa Major I     & {\it Fermi}-LAT, HAWC  & $97$ & $17.87^{+0.56}_{-0.33}$  \\
\CellTopTwo{}
Ursa Major II     & {\it Fermi}-LAT, HAWC, MAGIC  &$32$& $19.42^{+0.44}_{-0.42}$   \\
\CellTopTwo{}
Ursa Minor     & {\it Fermi}-LAT, VERITAS  & $76$ & $18.95^{+0.26}_{-0.18}$  \\
\hline
\hline
\end{tabular}
}
\caption{Summary of the relevant properties of the dSphs included in the combination of {\it Fermi}-LAT, HAWC, H.E.S.S., MAGIC, and VERITAS likelihood functions. The list of the observed dwarf galaxies is presented in column 1 with the instruments that performed the observations in column 2. Their heliocentric distance and  $J$ factor with their estimated $\pm 1\sigma$ uncertainties are listed in columns 3 and 4 respectively. The $J$ factors are given for a source extension truncated at the outermost observed star with their estimated $\pm 1\sigma$ uncertainties. }
\label{tab:dsphs_list}}
\end{table}

\subsection{Combination principle}

Our search for DM is carried out using a technique of maximum likelihood in which the profile likelihood ratio $\lambda$ is a function of the annihilation cross-section $\langle \sigma v \rangle$, \emph{i.e.} the parameter of interest, and reads as:

\begin{equation}
    \label{Eq:likelihood_profile}
\lambda \left( \langle \sigma \mathit{v} \rangle \mid \boldsymbol{\mathcal{D}_{\text{dSphs}}} \right) = \frac{\mathcal{L} \left( \langle \sigma \mathit{v} \rangle ; \boldsymbol{\hat{\hat{\nu}}} \mid \boldsymbol{\mathcal{D}_{\text{dSphs}}} \right)}{\mathcal{L} \left( \widehat{\langle\sigma v \rangle} ; \boldsymbol{\hat{\nu}} \mid \boldsymbol{\mathcal{D}_{\text{dSphs}}} \right)} ,
\end{equation}
where $\boldsymbol{\mathcal{D}_{\text{dSphs}}}$ is the dataset, $\boldsymbol{\nu}$ represents the nuisance parameters, $\widehat{\langle\sigma v \rangle}$ and $\boldsymbol{\hat{\nu}}$ are the values that maximize $\mathcal{L}$ globally, and $\boldsymbol{\hat{\hat{\nu}}}$ the values that maximize $\mathcal{L}$ for a given value of $\langle \sigma \mathit{v} \rangle$.
The joint likelihood function $\mathcal{L}$ describing all measurements is the product of the individual likelihood functions of all instruments and all dSphs and is given by:
\begin{equation}
    \label{Eq:dSph_combination}
    \mathcal{L} \left( \langle \sigma \mathit{v} \rangle ; \boldsymbol{\nu} \mid \boldsymbol{\mathcal{D}_{\text{dSphs}}} \right) = \prod_{l=1}^{N_{\text{dSphs}}} \mathcal{L}_{\text{dSph},l} \left( \langle \sigma \mathit{v} \rangle ; J_{l}, \boldsymbol{\nu_{l}} \mid \boldsymbol{\mathcal{D}_{l,\text{measured}}} \right) \times \mathcal{J}_l \left( J_{l} \mid J_{l,\text{obs}}, \sigma_{\log{J_l}} \right).
\end{equation}
The quantity $ N_{\text{dSphs}} = 20 $ is the total number of dSphs; $ \boldsymbol{\mathcal{D}_{l,\text{measured}}} $ is the dataset from gamma-ray observations for the $l$-th dSph; $ \boldsymbol{\nu_{l}} $ is the set of nuisance parameters associated to the $l$-th dSph, excluding $ J_{l} $; and $ J_{l} $ is the total $J$ factor of the $l$-th dSph, whose value can be found in Tab.~\ref{tab:dsphs_list}
; $\log_{10} J_{l,\text{obs}}$ and $ \sigma_{\log{J_l}} $ are obtained from the fit (see Jeans analysis in Sec.~\ref{sec:DM}) of a log-normal function of $J_{l,\text{obs}}$ to the posterior distribution of $ J_{l} $ \cite{2015PhRvL.115w1301A}.

\subsection{Shared data format}
In order to perform the combination of the observations, a table of test statistic (TS) values is provided by each experiment for the annihilation channels, $b\bar{b}$ and $\tau^+\tau^-$, for each set of $m_{\chi}$ and $\langle \sigma v \rangle$. All collaborations agreed on 63 DM masses ranging from 10~GeV to 100~TeV for all continuum channels following the mass spacing of~\cite{Cirelli} to avoid an interpolation. The $\langle \sigma v \rangle$ range is defined between $10^{-28}~\rm{cm}^3\rm{s}^{-1}$ and $10^{-18}~\rm{cm}^3\rm{s}^{-1}$ and is logarithmically spaced in 1001 values.

\subsection{Statistical uncertainty bands}
The $68\%$ (1$\sigma$) and $95\%$ (2$\sigma$) containment bands are derived by individual experiments by performing 300 Poisson realizations of the background events. Each collaboration provides the results of their statistical uncertainties in the same format as for the nominal values which are then combined following the same procedure as the combination of the nominal upper limits.

\section{Results and discussion}\label{Sec:Results}
No significant DM signal has been observed by any of the five instruments.
We therefore present the results of the combined upper limits at 95\% C.L. on the DM annihilation cross-section $\langle \sigma v \rangle$ in the case of two annihilation channels, $b\bar{b}$ and $\tau^+ \tau^-$, using all the data collected towards the twenty dSphs. We note that we selected these hadronic and leptonic channels as the follow up of our previous results presented at ICRC 2019~\cite{Oakes:2019ywx}.
We set our upper limits by solving $\mathrm{TS} = -2 \ln{ \lambda(\langle \sigma \mathit{v} \rangle)}$ for $\langle \sigma \mathit{v} \rangle$, with $\mathrm{TS}= 2.71 $. The value 2.71 represents the 95\% confidence level of a one-sided distribution assuming the test statistics behaves like a $\chi^2$ distribution with one degree of freedom. The combination is performed using two independent public analysis software packages, \texttt{gLike} \cite{javier_rico_2021_4601451} and \texttt{LklCombiner} \cite{tjark_miener_2021_4597500}, that provide compatible results.
The combined upper limits are presented in Fig.~\ref{fig:bb_tautau_UL} and are given with their 68\% ($1\sigma$) and 95\% ($2\sigma$) containment bands. These limits (solid black lines) are expected to be close to the median limit (dashed black lines) as no signal is present.
We obtain upper limits within the 2 $\sigma$ expected bands for the two annihilation channels $b\bar{b}$ and $\tau^+ \tau^-$.
The individual limits produced by each experiment are also indicated in the figures as a comparison to our new combined results.
Below \textasciitilde500~GeV, the DM limits are largely dominated by the {\it Fermi}-LAT experiment.
Between \textasciitilde500 GeV to \textasciitilde10 TeV, {\it Fermi}-LAT continues to dominate for the hadronic DM channel then above \textasciitilde10 TeV, the IACTs (H.E.S.S., MAGIC, and VERITAS) and HAWC take over.
In the case of the leptonic channel, both the IACTs and HAWC contribute significantly to the DM limit from \textasciitilde1 TeV to \textasciitilde100 TeV.

\vspace{-0.2cm}

\begin{figure}[h!]
\includegraphics[width=0.93\textwidth]{gd_ICRC2021_v3_crop.png}
\vspace{-0.1cm}

\caption{Upper limits at 95\% confidence level on $\langle \sigma \mathit{v} \rangle$ as a function of the DM mass for the annihilation channels $b\bar{b}$ (left) and $\tau^+ \tau^-$ (right), using the set of $J$ factors from Ref.~\cite{Geringer-Sameth:2014yza}. The black solid line represents the observed combined limit, the black dashed line is the median of the null hypothesis corresponding to the expected limit, while the green and yellow bands show the 68\% and 95\% containment bands. Combined upper limits for each individual detector are also indicated as solid, colored lines. }

\label{fig:bb_tautau_UL}
\end{figure}

We observe that the combined DM constraints from all five telescopes are 2 to 3 times stronger than any individual telescope for multi-TeV DM.
The selection of multiple targets increases statistics used to probe these sources and allows us to derive upper limits spanning the largest mass range of any WIMP DM search. We note that these limits depend on the choice of the annihilation channels and are driven by the objects with the highest $J$ factors that can be observed. The ultrafaint dSphs, containing a few tens of bright stars only, can be subject to large systematic uncertainties for the determination of their $J$-factors such as Segue I.
The derivation of upper limits through 6 additional annihilation channels is currently in progress, with 5 other continuum channels and the monoenergetic channel $\gamma \gamma$. A further analysis using a second $J$ factor set derived by~\cite{Bonnivard:2014kza, Bonnivard:2015xpq} is also yet to come in order to study the systematics induced by the choice of $J$ factor.

{\small
\bibliographystyle{ieeetr}