Floquet topological superconductivity induced by chiral many-body interaction

Non-equilibrium engineering is becoming a seminal way for realising novel quantum phases that are unimaginable in equilibrium. In particular, Floquet theory applied to quantum mechanics revealed that we can even control the band topology in semimetallic/insulating systems, while the straightforward application to topological superconductivity fails for typical superconductors because the supercondicting gap function does not couple to the electromagnetic field in a direct manner. Here we show that we can overcome this difficulty by taking account of correlation effects. Namely, we study how a d-wave superconductivity is changed when illuminated by circularly-polarised light (CPL) in the repulsive Hubbard model in the strong-coupling regime. We adopt the Floquet formalism for the Gutzwiller-projected effective Hamiltonian with the time-periodic Schrieffer-Wolff transformation. We find that CPL induces a topological superconductivity with a d + id pairing, which arises from the chiral spin coupling and the three-site term generated by the CPL. The latter term remains significant even for low frequencies and low intensities of the CPL. This is clearly reflected in the obtained phase diagram against the laser intensity and temperature for various frequencies red-detuned from the Hubbard U, with the transient dynamics also examined. The phenomenon revealed here can open a novel, dynamical way to induce a topological superconductivity. Driving quantum materials into non-equilibrium states using light-matter interactions is a way to induce novel quantum phases not attainable in equilibrium. Here, the authors theoretically demonstrate that circularly-polarised light can alter a d-wave superconductor in a strong-correlation regime into a topological superconducting state with broken time-reversal symmetry, as a combined effect of light and strong correlation.

I t 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 From being to becoming 1 , once said that he would have liked to entitle the book as 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 [2][3][4] . Namely, by applying a circularly-polarised light to honeycomb systems such as graphene [see Fig. 1a], 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 anom-alous Hall effect (i.e., quantum Hall effect in zero magnetic field) originally proposed for the static case by Duncan Haldane back in 1988 5 . Indeed, the effective Hamiltonian for the irradiated system in the leading order in the Floquet formalism exactly coincides with Haldane's model [ Fig. 1b], as shown by Takuya Kitagawa et al. 6 . The year 2019 witnessed an experimental detection of the Floquet topological insulator in graphene by James McIver et al. 7 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 in Table 1. 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 8 . Second, circularly-polarised laser lights can induce chiral spin interactions, ðŜ i Ŝ j Þ ÁŜ k , and consequent spin liquids for the strongly-coupled systems with the repulsion U much greater than the electron hopping t 0 [see Fig. 1c, d] [9][10][11][12] . There, we can realise that the Floquet-incuded chiral coupling becomes significant when the frequency ω of the laser is set close to the Hubbard U, with vastly different behaviours between the cases when ω 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 13 . Fourth, we can turn the usual pairing into an exotic η-pairing (condensation of pairs at Brillouin zone corners) by illuminating a superconductor with a (linearly-polarised) laser 14 . There, the Floquet effective pairhopping is drastically modulated when ω is close to U, again with Fig. 1 Various Floquet-induced quantum states. For Dirac fermions as in graphene, a illumination of a circularly-polarised light (CPL) induces the Floquet topological insulator, which arises from photon-assisted complex hopping (b), where colour shadings represent the piece-wise magnetic fluxes. c For magnetic phases in strongly-correlated systems, illumination of CPL induces a chiral spin coupling J χ (S i × S j ) ⋅ S k with S i being the spin at site i, for which d illustrates a resonant behaviour in the coupling strength at ω/U = 1, 1/2, 1/3, ⋯ with ω being the driving frequency, and U the strength of the onsite repulsive interaction. e For superconductors (present work), illumination of CPL on a d x 2 Ày 2 -wave superconductor produces pairing amplitudes hĉ i"ĉj# i across nearest neighbours (red: positive; blue: negative) along with imaginary diagonal bonds (magenta and green), leading to an emergent complexified gap function F(k) in k space, and we have a photo-induced chiral d x 2 Ày 2 þ id xy superconductivity (f), where arrows schematically indicate the phase of the complex gap function. different behaviours between ω 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 [15][16][17][18][19][20] , 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 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. 1f]. 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. 1e]. 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.

Results
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 whereĉ jσ annihilates an electron on the site j at R j with spin σ = ↑,↓, iσĉiσ is the density operator. For the hopping amplitude t ij , here we consider the second-neighbour (t 0 0 ) as well as the nearest-neighbour hopping (t 0 ), as necessitated for realising a CPL-induced topological superconductivity [as we shall see below Eq. (17)]. The onsite repulsion is denoted as U (>0).
The circularly-polarised laser field is introduced via the Peierls phase, which involves R ij ≡ R i − R j , and the vector potential, for the right-circularly-polarised case; replace E with E* for the left circulation. We set ℏ = e = 1 hereafter. Since the laser electric field is time-periodic, Hðt þ 2π=ωÞ ¼ĤðtÞ, we can employ the Floquet formalism 21 , 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ΛðtÞ withΛðt þ 2π=ωÞ ¼ΛðtÞ, we can introduce an effective static HamiltonianĤ F aŝ whereΛðtÞ is determined such thatĤ F becomes timeindependent. 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ΛðtÞ andĤ F in an exact manner, we can obtain their perturbative expansions in various situations. The high-frequency expansion 8, 22, 23 is a seminal example of such perturbative methods, where the effective HamiltonianĤ F is obtained as an asymptotic series in 1/ω, 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 ω, 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 ≪ U, we can employ the Floquet extension of the strongcoupling expansion (Schrieffer-Wolff transformation) 11,12,19,24,25 , when the photon energy ω exceeds the band width~t 0 of the doped holes (and chosen to be off-resonant with the Mott gap~U), i.e., t 0 ≪ ω < U. As we describe the detailed derivation in Methods, we can obtainΛðtÞ andĤ F in perturbative series in the hopping amplitudes t 0 ; t 0 0 . While the t 0 ≪ U assumption yields a low-energy effective model known in the undriven case as a t-J model 26 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 readŝ whereP G ¼ Q i ð1 Àn i"ni# Þ 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 0 ), 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, Γ i,j; k , are of second order, while the photon-generated chiral spin coupling, J χ 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 intot ij ¼ t ij J 0 ðA ij Þ due to the time-averaged Peierls phase 27,28 , where A ij = E|R ij |/ω and 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 24 which is a sum over m-photon processes, and appears in the second term on the first line of Eq. (5). 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 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 2, 6, 8 , 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. 2a, which is why a honeycomb lattice fits with the FTI.
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. 2b. This implies that such a time-reversal breaking term is present in (interacting) square lattice as well. The two-step correlated hopping Γ appears in the second-order perturbation in the presence of holes [the second line of Eq. (5)], which involves three sites and sometimes referred to as the "three-site term" 26 . While the three-site term also appears in the undriven case with the amplitude Γ 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,~g 2 δ in Eqs. (10), (11)]. In sharp contrast, when the system is driven by a CPL, Γ i,j; k acquires an important three-site imaginary part as with θ ij defined as R ij ¼ jR ij jðcos θ ij ; sin θ ij Þ. The term triggers a topological superconductivity even when it is small, as we shall see.
In addition to the two-step correlated hopping Γ, another timereversal breaking term arises if we go over to higher-order perturbations. In the previous studies of the half-filled case 11,12 , an emergent three-spin interaction (i.e., the scalar spin chirality term J χ ) is shown to appear in the fourth-order perturbation. This appears on the last line in Eq. (5) with a coefficient J χ ijk . While this term is basically much smaller than the second-order terms, it may become comparable with the contribution from Im Γ, as the spin chirality term does not accompany the dynamics of holes and has a larger Gutzwiller factor [~g 3 in Eqs. (10), (11). See also Fig. 3].
Mean-field decomposition under Gutzwiller ansatz. 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 26,29 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. (5) without the Gutzwiller factorP G if we renormalise the parameters. The mean-field solution is then formulated as a diagonalisation of the Bogoliubov-de Gennes Hamiltonian in the momentum space, Here we have where the expressions are independent of the site index i due to the spatial translational symmetry, and we have denoted the doping level as with N being the number of lattice sites. Here τ runs over τ = mx + ny with m; n 2 Z, and the label i + τ represents the site at R i+τ = R i + (m, n) in units of the lattice constant. The dependence on δ and a factor g = 2/(1 + δ) appears as a result of the Gutzwiller projection, which suppresses the contribution from charge dynamics as represented byt and Γ in the small δ regime. For the detailed derivation, see Methods. Throughout the present study, we take δ = 0.2. The mean-field Hamiltonian is self-consistently determined for the bond order parameter χ and the pairing amplitude Δ as 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 26 . The corresponding order parameter is for which the gap function F(k) becomes as derived from the first term in Eq. (11) with τ = ± x, ± y. Here we have dropped a small correction due to δ. Without a loss of generality, we assume that Δ is real. The d x 2 Ày 2 gap function has nodal lines along k y = ± k x , across which the gap function changes sign (see Fig. 4a-c below). Now let us look into the possibility for the laser field converting this ground state into topological superconductivity from the time-reversal breaking terms [the second line in Eq. (11)]. The leading term should be those for τ 0 ¼ ± x; ± y (with Δ τ 0 ¼ ± Δ), which is nonzero even for the original d x 2 Ày 2 -wave solution with no driving field and results in an imaginary gap function ∝ iΔ. In particular, gathering the terms with τ = ± (x + y), ± (x − y), we obtain the leading modulation to the gap function as where we have defined Thus we are indeed led to an emergence of a topological (chiral)  7) and ∑ lmn in Eq. (39)] are truncated such that the Taylor series up to E 20 is reproduced. This is accurate enough for ω > t 0 , while for the low-frequency regime ω < t 0 an intricate resonant structure due to the ω = U/integer resonance (seen also for t 0 < ω < 4t 0 ) is not captured.  1) has the second-neighbour hopping t 0 0 : Since the two-step correlated hopping Γ i,j; k is composed of two hoppings, k → i and j → k, we can see that the imaginary coefficient above has γ / t 0 t 0 0 . The chiral spin-coupling J χ / ðt 0 t 0 0 Þ 2 also necessitates the next-nearest-neighbour hopping, which can be deduced from Eq. (41) in "Methods".
Having revealed the essential couplings for the chiral superconductivity, let us now explore how we can optimise the field amplitude and frequency for realising larger topological gaps. For this, we set here t 0 0 ¼ À0:2t 0 and U = 12t 0 having cuprates with t 0 ≃ 0.4 eV in mind as a typical example.
We plot the essential γ and J χ against the driving amplitude E and frequency ω in Fig. 3a, 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 resonant fashion at ω/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 γ is given, in the leading order in the amplitude, as This expression, as a function of E, takes the maximal value around E ≃ 2.45ω/a, while diverges for ω → 0 or ω → 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 ≃ 2.4ω/a) can be small. If we turn to the spin-chirality term, J χ has a complicated form involving Bessel functions [see Eq. (41) below]. If we Taylorexpand it in E, we have which takes large values around ω = U due to the energy denominator in the above expression, while γ 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.
We can see in Fig. 3c-e how J χ , as well as the exchange interaction J and the renormalised nearest-neighbour hoppingt, vary with the CPL intensity for several representative frequencies for which the dynamical time-reversal breaking becomes prominent. The hopping amplitudet vanishes around the peak of γ (E ≃ 2.45ω/a), as seen in Fig. 3c, d. This involves a zero of the Bessel function (at E~2.40ω/a), and known as the dynamical localisation 27 . In Fig. 3d, e, we can see a strong enhancement of J. Panel d has ω = 0.45U (slightly red-detuned from U/2), while panel e has ω = 0.88U (slightly red-detuned from U). If we go back to Eq. (6), 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ω 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, γ attains a significantly large value ≃ 0.04t 0 , which yields F xy ðkÞ ' 0:3t 0 iΔ sin k x sin k y (for δ = 0.2). This is remarkably large and comparable with the undriven F x 2 Ày 2 ðkÞ ' 0:7t 0 Δð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.
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. 4. In the absence of the external field in Fig. 4a-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 ω = 10.5t 0 = 0.88U, E = 16t 0 /a for Fig. 4d-f, or with ω = 1.7t 0 = 0.14U, E = 3.5t 0 /a for Fig. 4g-i. The former represents the case where J χ is dominant (see Fig. 3e), while in the latter γ 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 Fig. 4f, 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 modifiedt and J.
To examine the robustness of the d + id superconductivity, we calculate |Δ| and Δ 0 ÀiΔ ± ðxþyÞ ¼ iΔ ± ðxÀyÞ at finite temperatures, which gives a phase diagram against T and field amplitude E in Fig. 5. We can see that the critical temperature T c , as delineated by the region for Δ ≠ 0, increases significantly as the  88U (a, b) or ω = 1.7t 0 = 0.14U (c, d), the superconducting order parameter Δ (a, c) and the time-reversal breaking order parameter Δ 0 ÀiΔ ± ðxþyÞ ¼ iΔ ± ðxÀyÞ (b, d) are colour-coded against the field amplitude E and temperature T. 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 0 ¼ À0:2t 0 , onsite interaction U = 12t 0 , and doping level δ = 0.2 as in Fig. 4. laser intensity E is increased. The enhanced T c originates from the fact that J [and Re Γ; see Eq. (38) below] are enhanced when the laser is applied, as we have seen in Fig. 3c-e. The time-reversal breaking order parameter Δ 0 emerges below T c down to T = 0, with the field dependence emerging from those of γ and J χ .
Transient dynamics. So far, we have investigated static properties of the effective HamiltonianĤ 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. (5) 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. (5) 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 H F (Eq. 5) with E ≠ 0. In such a treatment, hĤ F i is a conserved quantity for t > 0, and the driven state is expected to thermalise to the equilibrium state ofĤ F having a temperature that corresponds to the internal energy hĤ F i (which we call an effective temperature). Note that here we implicitly neglect the fact thatĤ F s for E = 0 and E ≠ 0 are on different frames [specified byΛð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. 5. 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 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~4t 0 /a in Fig. 5c, d is due to the band flipping (sign change of 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 13 . 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 |Δ| and jΔ 0 j in Fig. 6a, b, where we set the initial temperature at k B T initial = 0 with E = 7t 0 /a, ω = 10.5t 0 = 0.88U. Here, the unit of time ℏ/t 0 corresponds to~1 fs for t 0 ≃ 0.4 eV. We can see that |Δ| 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 30 . The time-reversal breaking order parameter jΔ 0 j 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Ĥ F , where the deviation arises because the pairbreaking scattering 30 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 Higgs modes in superconductors [31][32][33] , and thus the typical time scale for the oscillation (and the emergence of Δ 0 ) can be roughly estimated as the inverse of the superconducting gap (~1/2|F(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 34 if the change in the coupling constant is small, which reveals that the appearance of Δ 0 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. 5, we find an intriguing phenomenon: the superconducting state can appear even when we start from an initial temperature that is 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 "nonequilibriuminduced superconductivity" in Fig. 6c, d, where we set k B T initial = 0.2t 0 as an initial temperature, and choose E = 16t 0 /a, ω = 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 Δ = 10 −4 on the initial state. After the quench, both of the gap functions Δ and Δ 0 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 (|Δ| ≃ 0.18 and jΔ 0 j ' 0:03) expected for equilibrium with the effective temperature.

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 Fig. 6 Quench dynamics of the order parameters. Time evolution of |Δ| = |〈c i↑ c i+x,↓ 〉| and jΔ 0 j ¼ jhc i" c iþxþy;# ij after a sudden quench of the field amplitude from E = 0 to E = 7t 0 /a with an initial temperature k B T initial = 0 (a, b), or to E = 16t 0 /a with an initial k B T initial = 0.2t 0 (c, d). The latter corresponds to the case in which the superconductivity appears even when we start from an initial temperature above T c . |Δ| (a, c) represents the conventional d x 2 Ày 2 -wave component, while jΔ 0 j (b, d) represents the time-reversal breaking id xy -wave component. Horizontal lines indicate the equilibrium value at the effective temperature. The driving frequency is set to ω = 10.5t 0 = 0.88U with t 0 being the bare hopping amplitude.
proposal for the Floquet topological superconductivity. In the previous study of Floquet topological superconductivity in cuprates 17 , the topological transition is triggered by the modulation of the kinetic part of the Hamiltonian [diagonal components in Eq. (8)], 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 B T c (since the F x 2 Ày 2 ðkÞ and F xy (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~1 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Ĥ F in the thermodynamic limit at long times are believed to represent a featureless, infinitetemperature state (sometimes called the Floquet eigenstate thermalisation hypothesis [35][36][37][38][39][40] ). 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 [40][41][42][43], 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ω) with integers n, m, in which the enhancement of γ and J χ 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~10/t 0 for the emergence of Δ 0 (see Fig. 6b) lies within the validity of the Floquet t-J description, with more sophisticated methods including direct analyses of a timedependent 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 projectorP 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 44 along with related theories [45][46][47][48] , 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 ω = U/2 or ω = U, the required field intensity is E~10t 0 /a, which corresponds to E~10 2 MV/cm for the typical cuprates with t 0 ≃ 0.4 eV, a ≃ 3 Å. 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 timereversal breaking terms can be of second order 11,12 . This comes from the cancellation similar to Fig. 2a for the present square lattice [see Eq. (44) 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 γ, since the required intensity scales with the frequency ω. 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 49 ) 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.

Methods
Time-periodic Schrieffer-Wolff transformation. To obtain the effective lowenergy Hamiltonian in the presence of the laser electric field, we employ the timeperiodic Schrieffer-Wolff transformation (a canonical transformation). Following a previous study 11 (but extending it for the case where holes exist; see also ref. 19 ), we decompose the Hubbard Hamiltonian aŝ where λ is a bookkeeping parameter bridging the atomic limit λ = 0 to the system of interest at λ = 1, andD ∑ ini"ni# counts the number of doubly-occupied sites. The hopping operatorT d;m describes the process that increases the number of doubly-occupied sites by d, as given bŷ where σ Àσ. We also decompose the generatorΛðtÞ in Eq. (4) into e iΛðtÞ ¼ e iΛ h ðtÞ e iΛ c ðtÞ , whereΛ c eliminates the charge excitations (doubly-occupied sites) from the transformed Hamiltonian, whileΛ h makes the Hamiltonian static.
We first considerĤ c ðtÞ :¼P G e iΛ c ðtÞ ðĤðtÞ À i∂ t Þe ÀiΛ c ðtÞP G . Let us introduce a series solution,Λ While this Hamiltonian projected onto the spin subspace is evaluated in a previous study 11 , here we need to consider the expression for nonzero numbers of holes. We then obtain P GT þ1;n ;T À1;mÀn Note that, while we recover the Heisenberg spin-spin interaction for k = i, the socalled three-site terms with k ≠ i arise as well.
where we have used fĉ iσ ð1 Àn i σ Þ; ð1 Àn j σ 0 Þĉ y jσ 0 g ¼ δ ij ½δ σσ 0 ð1 Àn i =2Þ þ σ σσ 0 ÁŜ i , fð1 Àn i σ Þĉ y iσ ; ð1 Àn j σ 0 Þĉ y jσ 0 g ¼ 0. The effective Hamiltonian Eq. (4) is obtained up to the second order aŝ Note that the second term in the hoppingt ij describing the two-step hopping vanishes for the square lattice [see Fig. 2a], and the renormalised hopping is given ast ij ¼ t ð0Þ ij ¼ t ij J 0 ðA ij Þ if we put m = 0 in Eq. (26). For the Heisenberg term J ij we end up with Eq. (6) in the main text. The two-step correlated hopping is much more intricate, and reads In the main text, we have also considered the scalar spin-chirality term J χ ijk , which appears in the fourth-order perturbation. While the perturbative processes involving three sites are considered in previous studies 11,12 , we need to evaluate the processes involving four sites as well for the present case of the square lattice as we see in Eq. (44) below. By dropping the density-dependent part, where The first term on the second line (three-site term) coincides with the result in a previous study 11 . While the second term (four-site term) is absent for e.g. honeycomb and kagome lattices with no candidates for the fourth site h, square lattice accommodates this additional contribution. Then we can note that we have a combination of blue and red paths as depicted in Fig. 2a, which was invoked for the three-site term but also applies to the four-site term involving t kj t ji and t kh t hi . There, due to the factor R × R, the first and second terms cancel with each other for the square lattice (even if we consider many-body effects). The square lattice model does have a nonzero coefficient, however, if we go over to E 4 in the presence of the next-nearest-neighbour hopping, as we have seen Eq. (21).
Gutzwiller projection. Here we analyse the ground-state property of the system using the mean-field approximation with the Gutzwiller ansatz 26,29 . Namely, we consider an ansatz for the wavefunction, withP G ¼ Q i ð1 Àn i"ni# Þ being the Gutzwiller projection. Here,ĉ kσ ¼ N À1=2 ∑ iĉiσ e ÀikÁR i with N being the number of lattice sites. The ground state within this ansatz can be obtained by minimising the expectation value, To evaluate the Gutzwiller projection approximately, we replace the expectation value with the site-diagonal ones, where we decompose the projection operator aŝ withP ih ð1 Àn i" Þð1 Àn i# Þ;P iσ n iσ ð1 À n i σ Þ. Then we can evaluate the denominator as where f ¼ ð1 À δÞ=2; f ¼ ð1 þ δÞ=2, and we have used n!~(n/e) n on the last line. We can evaluate the numerators in the same manner. For the hopping, exchange, and spin-dependent three-site terms, we obtain, respectively, hP GŜi ÁŜ jPG i 0 ' ð f 2 δ À1 Þ Nδ f 2Nf À2 hŜ i ÁŜ j i 0 ð54Þ We also obtain hP GŜi Á ðŜ j Ŝ k ÞP G i 0 ' ½8=ð1 þ δÞ 3 hŜ i Á ðŜ j Ŝ k Þi 0 hP G i 0 in the same manner.
The remaining terms have an ambiguity in the Gutzwiller factor, because they are accompanied by the density operatorP GniPG ¼P G ½∑ σniσ ð1 Àn i σ ÞP G , which can be regarded either as an operator or a constant hn iσ ð1 Àn i σ Þi 0 ' f f , in the above scheme. Here, we discard the second-order fluctuation around the expectation value, ðn iσ À f Þð1 Àn i σ À f Þ, to approximate the projected density operator as ∑ σniσ ð1 Àn i σ Þ 'n i δ þ 2f 2 . Then the Gutzwiller factors for the remaining terms are evaluated as hP GninjPG i 0 ' 4 ð1 þ δÞ 2 hðn i δ þ 2f 2 Þðn j δ þ 2f 2 Þi 0 hP G i 0 ¼ 4δ 2 ð1 þ δÞ 2 hn inj i 0 hP G i 0 þ 2δð1 À δÞ 2 ð1 þ δÞ 2 hn i þn j i 0 hP G i 0 þ const:; ð58Þ hP Gĉ Note that the density-density term is known to have small effects in the variational Monte Carlo calculation in the low-doping regime, with which the present treatment is consistent. With these expressions, the minimisation of E 0 turns out to be reduced to that of hĤ F0 i 0 , whereĤ F0 is the Hamiltonian with the modified coupling constant but without the Gutzwiller projection, as given bŷ