Time-dependent Gutzwiller simulation of Floquet topological superconductivity

Periodically driven systems provide a novel route to control the topology of quantum materials. In particular, Floquet theory allows an effective band description of periodically-driven systems through the Floquet Hamiltonian. Here, we study the time evolution of $d$-wave superconductors irradiated with intense circularly-polarized laser light. We consider the Floquet $t$-$J$ model with time-periodic interactions, and investigate its mean-field dynamics by formulating the time-dependent Gutzwiller approximation. We observe the development of the $id_{xy}$-wave pairing amplitude along with the original $d_{x^2-y^2}$-wave order upon gradual increasing of the field amplitude. We further numerically construct the Floquet Hamiltonian for the steady state, with which we identify the system as the fully-gapped $d+id$ superconducting phase with a nonzero Chern number. We explore the low-frequency regime where the perturbative approaches in the previous studies break down, and find that the topological gap of an experimentally-accessible size can be achieved at much lower laser intensities.


INTRODUCTION
Topological superconductors host robust gapless excitations at the boundary or vortex cores due to the topological structure of the superconducting gap function [1][2][3][4].In particular, Majorana fermions that emerge in topological superconductors provide a platform for fault-tolerant quantum computation [5].Theoretical proposals for creating a topological superconductor include topological insulators in the proximity of s-wave superconductors [6] and semiconductors with spinorbit coupling in the proximity of s-wave superconductors [7].Despite intense experimental efforts to confirm topological superconductivity and Majorana fermions in those setups, their existence is still elusive [8].Hence, seeking an alternative platform for topological superconductivity remains an important issue.
Periodically driven systems provide a novel route to control the topology of quantum materials.In particular, Floquet theory allows an effective band description of periodicallydriven systems through the Floquet Hamiltonian.Thus the dynamical control of quantum phases has recently been studied actively, called "Floquet engineering" [9][10][11][12][13].A canonical example of Flouqet engineering of a topological phase is the quantum anomalous Hall state that emerges in graphene irradiated by circularly-polarized light (CPL) [14][15][16][17].In graphene, CPL induces an effective complex hopping for the next-nearest neighbor in the Floquet Hamiltonian, which takes the same form as in the Haldane model for the quantum anomalous Hall state [18].
Applying the concept of Floquet engineering to topological superconductors, Floquet topological superconductivity has been explored [19][20][21][22][23][24][25].For example, a honeycomb lattice with s-wave pairing interaction is predicted to exhibit topological superconductivity under the irradiation of CPL [19,23], where CPL induces a mass term to gapless excitations around point nodes, leading to topological superconductivity.A similar strategy for Floquet topological superconductivity was also pursued for cuprate superconductors.Specifically, a dwave superconductor on a square lattice was shown to support CPL-induced topological superconductivity when strong spinorbit coupling (SOC) is present [21].All these approaches to Floquet topological superconductivity essentially rely on the presence of additional internal degrees of freedom supporting nontrivial geometry (i.e.sublattice in the honeycomb lattice or spins with SOC) for turning the systems into topological phases.
Recently, it was revealed that d-wave superconductors exhibit Floquet topological superconductivity purely from the many-body effect without invoking the internal degrees of freedom [26].To incorporate strong correlation effects, they derive the Floquet t-J Hamiltonian using the Schrieffer-Wolff transformation [27][28][29][30] and the high-frequency expansion (HFE) [10,31,32].Time-reversal symmetry breaking terms appear from the interaction terms in the Floquet t-J model and induce topological d x 2 −y 2 + id xy pairing upon a mean field treatment (See Fig. 1).This approach enables us to broaden the class of candidate materials for the Floquet topological superconductivity.
However, experimental implementation of the above theoretical proposals is still challenging because of the required field intensity of ∼ 100 MV/cm [21,26].This stringent requirement essentially stems from the fact that the HFE has been employed in these previous studies.The field-induced effective coupling in the HFE typically scales with the am-arXiv:2309.06069v2[cond-mat.str-el]7 Dec 2023 plitude of the vector potential, A ∝ E/ω, which implies that a strong electric field E is necessary for a large driving frequency ω.Also, driving in the low frequency regime below the electronic band gap is desirable to avoid heating of the system and achieve a coherent control.Thus it is essential to develop a theoretical framework that is applicable to the low frequency regime.
In this paper, we study Floquet topological superconductivity at low frequencies by performing a time-dependent Gutzwiller simulation.We formulate the time-dependent Gutzwiller approximation based on the action principle combined with the mean-field approximation.Specifically, to simulate the time evolution of the many-body state in a tractable way, we apply the Gutzwiller approximation to the Lagrangian formalism to deduce a time-dependent Schrödinger equation with a mean field approximation for pairing amplitudes and bond orders.The time-dependent simulation shows that the d x 2 −y 2 -wave superconductor evolves into the topological d x 2 −y 2 + id xy -wave superconductor under the CPL driving, both in the high and low frequency regimes.We further analyze the obtained time-periodic superconducting state in terms of the Floquet Hamiltonian, revealing the full-gap nature and the nontrivial winding of the gap function.We find that the topological gap of the order of 3 K emerges for an electric field of ∼ 10 MV/cm, which will be feasible for experimental measurements.

Formalism
In this section, we derive an effective low-energy Bogoliubov-de Gennes (BdG) Hamiltonian in the presence of the CPL, employing time-periodic Schrieffer-Wolff transformation (a canonical transformation) with Gutzwiller ansatz [26].
We consider a periodically-driven Hubbard model defined on a square lattice, having in mind a cuprate superconductor.The time-dependent Hamiltonian is given by ĤHub where ĉiσ is an electron annihilation operator at site i with spin σ =↑, ↓, and niσ = ĉ † iσ ĉiσ is the spin-density operator.Here we set ℏ = e = 1 for simplicity.The latter term is the on-site Coulomb repulsion terms with the Hubbard interaction U, while the former term is the hopping terms with modulated hopping amplitude t i j e −iA(t)•R i j , where t i j is the hopping amplitude between site i and site j.We introduce R i j = R i − R j , where R i is the location of the i-th site.Here we consider CPL, for which the vector potential A(t) is given by To deduce low-energy dynamics of the driven Hubbard model, we consider the Lagrangian of this system, where |Ψ(t)⟩ is a state vector of the many-body system.First, we perform the time-periodic Schrieffer-Wolff transformation [28,29], where the transformed state vector is represented as PG e i Ŝ (t) |Ψ(t)⟩, with the unitary transformation e i Ŝ (t) and the Gutzwiller projection PG = i (1 − ni↑ ni↓ ).Here Ŝ (t) should be chosen such that the transformed Hamiltonian, ĤSW (t) ≡ PG (e i Ŝ (t) ĤHub e −i Ŝ (t) − e i Ŝ (t) (i∂ t e −i Ŝ (t) )) PG , becomes diagonal in the charge sector (eliminating charge excitations), and thus commutes with PG (For a similar method for the Hubbard model not based on the Schrieffer-Wolff transformation but with a generalized projection operator, see Refs.[33,34]).Then the transformed state vector only takes the configurations that have no doubly-occupied site.By adopting the Gutzwiller ansatz where e i Ŝ (t) |Ψ(t)⟩ is chosen to be the BCS wave function |0⟩, here we approximate the original Lagrangian L by L G , as We conduct the Schrieffer-Wolff transformation up to the second-order of the hopping and obtain [26] ĤSW where Ŝ j is a spin operator.The coupling constants are given by Ji Γi Here, t (m) i j is the m-th Fourier component of the modulated hopping amplitude t i j e −iA(t)•R i j , whose explicit form is ob-FIG.2. Schematic picture of the id xy -wave component of the pairing amplitude ∆ idxy that arises from three-site terms.Through the virtual hopping process from the site k to the site i via the site j, the id xywave component ∆ idxy (magenta) is induced from the original d x 2 −y 2wave order (orange).
tained by the Jacobi-Anger expansion as with the m-th Bessel function J m (x) and Θ i j defined as the polar angle of R i j , i.e., R i j = |R i j |(cos Θ i j , sin Θ i j ).In the absence of the external field, Eq. ( 6) is known as the t-J model [35].The first term is the hole hopping term in the configurations that have no doubly-occupied site.The second term is composed of Heisenberg interaction Ŝi • Ŝ j and the density-density interaction ni n j .The third term, representing an interaction of the form σσ is the so-called three-site term.This Hamiltonian is known to yield singlet (d-wave) Cooper pairing in equilibrium, so that here we assume the singlet pairing as well.A remarkable point here is that the three-site terms in the present case break time-reversal symmetry, due to the circular polarization of the light field.When we consider the spin part of the three-site terms σσ ′ ĉ † iσ σ σσ ′ ĉkσ ′ • Ŝ j where ij is a next-nearest neighbor bond and j-k is a nearest neighbor bond, the three-site terms generate pairing amplitude on each bond as j↑ and ĉ j↑ ĉk↓ correspond to d xy -wave components and d x 2 −y 2 -wave components of pairing amplitudes respectively.When the coefficients of the three-site terms Γ(t) in Eq. ( 6) become complex due to the CPL, the next-nearest neighbor pairing amplitude ĉ † i↓ ĉ † j↑ becomes complex, embodying the id xy -wave component of the pairing amplitude (as shown in Fig. 2).The same applies to the density part of the three-site term σ ĉ † iσ ĉkσ n j .These terms have already been shown to induce topological superconductivity by performing the HFE of Eq. ( 6) in Ref. [26].The time-periodic Schrieffer-Wolff expansion is a perturbative expansion in condition t 0 ≪ U meanwhile there is no resonance between ω and U (i.e.t (m)  i j ≪ U − nω).For the detailed derivation, see Schrieffer-Wolff transformation subsection in Methods.
Second, we perform the Gutzwiller approximation [33,34,36,37], which replaces the Gutzwiller projection PG with c-numbers renormalizing each term in the Hamiltonian and effectively incorporates the reduction of double occupancies.
We perform the Gutzwiller approximation to the second term of Eq. ( 5) and derive an effective Hamiltonian ĤG as Hereafter, ⟨• • •⟩ represents the expectation value in terms of BCS wavefunction |Ψ BCS (t)⟩.As described in Gutzwiller approximation subsection in the Methods, we obtain the effective Hamiltonian ĤG (t), where , and f = 1 − f .Similarly, we perform Gutzwiller approximation to the first term of Eq. ( 5) and obtain For details of derivation, see Gutzwiller approximation subsection in Methods.Third, we derive the BdG Hamiltonian from the timedependent Gutzwiller Hamiltonian in Eq. (12).For order parameters, we consider two SU(2)-symmetric orders, i.e., the bond order amplitude χ τ (t) and the superconducting pairing amplitude ∆ τ (t), where τ = mx + ny represents a bond connecting two sites that are distant by ma in the x direction and na in the y direction, with a being the lattice constant.By applying the meanfield approximation to each term in Eq. ( 12) and performing a Fourier transformation, we arrive at the effective BdG Hamiltonian in the momentum-space representation as (g) Notations for the superconducting pairing amplitudes on the bonds and their color codes for (b-e,h,i).We adopted the parameter set: the driving frequency ω = 5.8t 0 , the next nearest neighbor hopping t ′ 0 = −0.2t0 , the onsite interaction U = 12t 0 and the hole doping level δ = 0.2. with Here, matrix elements are given by Finally, we obtain the effective Lagrangian as L G ≃ ⟨Ψ BCS (t)| (i∂ t − ĤBdG (t)) |Ψ BCS (t)⟩.By using the action principle with δu * k (t), δv * k (t) as variants, we end up with the timedependent Schrödinger equation, where ⃗ ψ k (t) = (v k (t), u k (t)) T .We solve this Schrödinger equation (20), by computing the order parameters χ τ and ∆ τ in Eqs. ( 14) and ( 15) consecutively.

Time evolution of superconducting pairing amplitudes
In this section, we show the results of the time evolution generated by Eq. (20).We consider hopping amplitudes, bond order amplitudes, and superconducting pairing amplitudes up to next-nearest neighbors.
Throughout this paper, we set the next-nearest-neighbor hopping t ′ 0 = −0.2t0 , the onsite interaction U = 12t 0 , and the hole doping level δ = 0.2, with t 0 being the nearest-neighbor hopping amplitude.As for the other doping levels and onsite interaction, see Figs.S1-S3 in Supplementary Note 1.
We show the time evolution of superconducting pairing amplitudes ∆ τ for ω = 5.8t 0 in Figs.3(b,c).The time profile of the driving field amplitude is shown in Fig. 3(a), and the convention for the superconducting pairing amplitudes is depicted in Fig. 3(g).As shown in Fig. 3(a), we increase the electric field amplitude linearly until ωt/2π = 200 and keep it constant after that.We adopted this time profile of E rather than just quenching the field, to connect the initial state to the dynamically-stabilized topological phase while mitigat- ing damping oscillations in the order parameters.We show four superconducting pairing amplitudes: 1] in Fig. 3 (Other components can be obtained from the relationship ∆ [x,y] = ∆ [−x,−y] for singlet pairing).We perform simulations with the initial state in the d x 2 −y 2 -wave superconducting state which is the ground state of Eq. 16 in the absence of the field.Namely, as shown in Figs.3(b,c), we set positive ∆ [1,0] and negative ∆ [0,1] with the same magnitude while the others are zero, with which the d x 2 −y 2 -wave component given by (Re∆ [1,0] − Re∆ [0,1] )/2 has a nonzero value.
Once we turn on the electric field, the phases of the superconducting pairing amplitudes start to rotate, as shown in Figs.3(b,c).To make it easier to see the relative phases between different ∆ τ 's, we perform a phase rotation with the phase θ(t) of the pairing amplitude ).In Figs.3(d,e ])/2, which indeed belongs to the same irreducible reprenentation of C 4v as the d xy pairing.Figure 3(f) shows the time evolution of ∆ id xy , and Fig. 3(h) is its blowup of the shaded area.While ∆ id xy is zero in the initial state without the CPL, ∆ id xy becomes finite once the CPL is applied, clearly indicating that the CPL irradiation induces id xy -wave component of the pairing amplitude leading to topological superconductivity.
Next, to study the magnitude of the id xy -wave pairing amplitude in the steady state, we compute the time average of ∆ id xy over ωt/2π ∈ [240, 400], denoted as ∆ id xy .We show a color plot of ∆ id xy as a function of the driving amplitude E and the frequency ω in Fig. 4(a).Note that for ω = 4t 0 , 6t 0 , Eqs. ( 18) and ( 19) contain divergent terms with the denomi-nator U − nω for U = 12t 0 , with which the effective Hamiltonian becomes ill-defined because of breaking down of the Schrieffer-Wolff transformation.We avoided those parameters in plotting Fig. 4(a).For comparison, we show the id xy -wave component ∆ HFE id xy obtained from the high-frequency expansion of ĤSW in Fig. 4 is the superconducting pairing amplitude in the ground state of the Floquet Hamiltonian obtained by applying the HFE to Eq. 6 [26].The present time-dependent Gutzwiller approach does not rely on the HFE and is also applicable to the low-frequency region in contrast to the HFE approach.Here we fix Ea/ω = 1 and compare ∆ id xy with ∆ HFE id xy as functions of ω in Fig. 4(c).∆ id xy and ∆ HFE id xy show similar behaviors in the high-frequency region ω > 5t 0 , while they are significantly different in the low-frequency region ω < 2t 0 .∆ id xy even shows a sign change around ω ∼ t 0 as opposed to ∆ HFE id xy .This topological phase transition can also be seen in Fig. 4(d), a color plot of ∆ id xy as a function of E and ω in the low-frequency region.Specifically, the topological phase transition clearly occurs around E = 1.05t 0 /a.In this way, the present method can capture detailed behaviors of the Floquet topological superconductivity even in the low-frequency region where the HFE is not applicable.In particular, topological phase transitions can be induced by changing the field strength of CPL in the lowfrequency region.

Floquet theory analysis
The Floquet theory is a useful tool to describe periodically driven systems via an effective static Hamiltonian (Floquet Hamiltonian).The topological properties of the steady state in the CPL driven system are characterized by the Chern number of the Floquet Hamiltonian.In this section, we show the Floquet bands and Chern numbers obtained in the present method.

Floquet Hamiltonian
The Floquet Hamiltonian is defined as H(k, t)dt for a time-periodic Hamiltonian H(k, t) of a period T .As can be seen from Figs. 3(b,c,h,i), ĤBdG (t) is not time-periodic due to the phase rotation of ∆ τ over a period, while ∆ τ e −iθ becomes (approximately) time-periodic.In the following, we derive the Floquet Hamiltonian using this phase rotation for ∆ τ .
In the present system, once the system arrives at the steady state, the superconducting pairing amplitudes satisfy as can be seen from Figs. 3(b,c).We can remove this phase rotation by a gauge transformation ⃗ ψ ′ k ≡ e iϵtτ z ⃗ ψ k with −ϵT ≡ α/2, where ⃗ ψ ′ k becomes time-periodic.Here, τ z is a Pauli z matrix acting on the Nambu space.Correspondingly, the Hamiltonian becomes time-periodic within this gauge and enables us to define the Floquet Hamiltonian as where U(T, 0) is the original time evolution operator.For detailed derivation, see Time-periodic Hamiltonian and Floquet Hamiltonian subsection in Methods.

Floquet band and Chern number
We compute the Floquet Hamiltonian for H(k, t) in Eq. ( 20) with the above method and obtain the associated Floquet bands.We then compute the Berry curvature and the Chern number for the Floquet bands by using the Fukui-Hatsugai-Suzuki method [38].Note that the results below do not take into account the occupation of the Floquet band.
Figure 5 shows the obtained Floquet bands, the phases of gap functions F(k) of the Floquet Hamiltonian, and the Berry curvature.The results for the system without driving fields are shown in Figs.5(a,d,g).Figure 5(d) indicates that the gap function has the conventional d x 2 −y 2 -wave symmetry with nodal lines.The black curve in Fig. 5(d) indicates the Fermi surface in the normal state.The point nodes appear at the crossing points of nodal lines and the Fermi surface as seen in Fig. 5(a).
Next, let us look at the results for the systems with the CPL drivings.We study two cases of CPL driving.One is in the high-frequency regime with the driving frequency ω = 5.8t 0 and the field amplitude E = 5.4t 0 /a, shown in Figs.5(b,e,h).The other is in the low-frequency regime with ω = 1.03t 0 and E = 0.96t 0 /a, shown in Figs.5(c,f,i).With the CPL drivings, we find that the point nodes are gapped out in the Floquet band structures in both cases in Figs.5(b,c).We show the phase of gap functions obtained from the Floquet Hamiltonians with the CPL drivings in Figs.5(e,f).We remark that the Floquet Hamiltonian H(k) and the gap function depend on the initial time t 0 of the time evolution operator U(t 0 + T, t 0 ) (while the eigenvalues are independent).Although the gap function shown in Figs.5(e,f) breaks the C 4 symmetry, the C 4 -breaking component rotates with t 0 , so that the Floquet state preserves C 4 symmetry in time average.
In these cases, the black curves correspond to the k points where the diagonal (normal) component of the Floquet Hamiltonian is zero (Fermi energy in the equilibrium cases) or ω/2 (Floquet Brillouin zone boundary).In the high-frequency case (ω = 5.8t 0 , E = 5.4t 0 /a), the black curve in Fig. 5(e) almost coincides with the Fermi surface in the equilibrium case [Fig.5(d)].Due to the emergence of the id xy -wave component of the pairing amplitude, the phase of the gap function rotates twice in the counterclockwise direction along the black curve [Fig.5(e)].As a consequence, a negative Berry curvature appears around the area where the point node was orig- inally located [Fig.5(h)].In the low-frequency case (ω = 1.03t 0 , E = 0.96t 0 /a), Fig. 5(f) shows additional black curves in addition to the one corresponding to the Fermi surface in the equilibrium.These additional black curves correspond to band crossings arising from the folding of copies of Floquet bands (photon-dressed states) onto the Floquet Brillouin zone.If we focus on the black curve from the original Fermi surface, the phase of F(k) rotates twice in the clockwise direction along the black curve in Fig. 5(f).Correspondingly, a positive Berry curvature appears around the gapped node in Fig. 5(i).As can be seen from Figs. 5(h,i), the Berry curvature B(k) satisfies the relation B(k) = B(−k), which indicates that the time-reversal symmetry is broken while the inversion symmetry is preserved, as is consistent with the d x 2 −y 2 + id xy pairing.
Furthermore, the sign of the Berry curvature coincides with the sign of the id xy -wave component of the pairing amplitude ∆ id xy in Fig. 4.This clearly indicates that the CPL induced Berry curvature originates from the gap opening at the point node with the id xy -wave component ∆ id xy .Specifically, the electronic structure of the original point node is described by a gapless Dirac fermion.Once the mass gap is introduced by ∆ id xy to the Dirac fermion, the Berry curvature of the sign of the mass ∆ id xy emerges.
We show the phase diagram of Floquet topological superconductivity by computing Chern numbers of the Floquet bands in Fig. 6.We show the phase diagrams for the highfrequency region in Fig. 6(a) and the low-frequency region in Fig. 6(b).The location of topological phases in Fig. 6 is generally consistent with the sign of the id xy -wave component of the pairing amplitude ∆ id xy in Fig. 4. Again, we avoid ω = 4t 0 , 6t 0 in plotting Fig. 6 because of the ill-defined effective Hamiltonian with the divergent denominator U − nω in Eqs.(18) and (19).We note that there are some regions with C ±2 showing jumps of the Chern number in Fig. 6.They are mainly due to the Berry curvature around the additional black curves in Fig. 5(i).Since they appear from the folding of copies of Floquet bands, usually the occupation of energy states does not change abruptly across the energy gap (See Occupation of Floquet band subsection in Methods for detail).Namely, although those Floquet bands have large Chern numbers, the states in the both lower and upper bands around those minigaps are equally occupied and their contributions to the expectation value of the (thermal) Hall response in state |Ψ(t)⟩ cancel out.Therefore, Berry curvature around the original point node gives a dominant contribution to the (thermal) Hall response when the electron occupation is not far from the equilibrium one.

DISCUSSION
We have demonstrated the emergence of topological superconductivity even at low frequencies.To discuss the experimental feasibility, let us turn our attention to evaluating the size of topological gaps and the required magnitude of electric fields and frequencies.The low-frequency region shown in Fig. 6(b) corresponds to ℏω ≃ 0.4-0.5 eV, E ≃ 8-17 MV/cm for the typical cuprates with t 0 ≃ 0.4 eV, a ≃ 3 Å.One of the advantages of the present time-dependent Gutzwiller approach in the low-frequency regime is that it does not require the large electric fields as large as 100 MV/cm as mentioned in previous research [21,26], due to the effective coupling through the vector potential A ∝ E/ω.The size of the topological gap appearing in Fig. 5(c) is approximately 0.03 times that of the original d x 2 −y 2 -wave superconducting gap (the gap function at the Fermi surface in the antinodal direction).Therefore, assuming that T c of cuprate is 89 K [39], experimental observation of the topological gap can be achieved at around 3 K.While the magnitude of the topological gap is smaller than that indicated in Ref. [26], the gap described in that reference scales as ∝ O(E 4 ).At E ∼ 10 MV/cm, the magnitude of the topological gap obtained by the present method is overwhelmingly larger.In addition, as for the impurity effect, since the present d + id-wave superconductivity is based on the d x 2 −y 2 -wave superconductivity in the undriven cuprates, we consider that the suppression of the topological superconductivity by the impurities will be qualitatively the same as that for the original d-wave superconductivity in cuprate superconductors.
Several experimental methods are considered for probing the Floquet topological superconductivity demonstrated in this paper.The first is the measurement of the optical conductivity which directly observes the topological gap.Several studies have already investigated the optical conductivity of cuprate superconductors [40,41].The irradiation of CPL opens the gap at the Dirac nodes in cuprate superconductors.Therefore, during CPL irradiation, optical absorption at frequencies below the topological gap is expected to be suppressed.The second is to detect the Higgs mode associated with the id xy -wave component of the superconducting pairing amplitudes.The d x 2 −y 2 -wave superconductivity transforms into d x 2 −y 2 + id xy -wave superconductivity under CPL irradiation.The Higgs mode of d-wave superconductivity in cuprate superconductors was investigated using terahertz pump and optical probe techniques [42,43].Thus, observing the Higgs mode of the additional id xy -wave component would provide evidence of Floquet topological superconductivity.Another approach involves measuring the Hall response.CPL irradiation breaks time-reversal symmetry and induces a finite Berry curvature.The realization of a QAHI in graphene through CPL irradiation was recently observed using a laser-triggered photoconductive switch [17].Similar methods may capture the Hall response in the present system.

Schrieffer-Wolff transformation
In this section, we derive Eq. ( 6).For clarity, we rewrite each term in ĤHub as where Td (t) is the kinetic energy operator that changes the double occupancy by d: [ D, Td (t)] = d Td (t) with D being the double occupancy operator.The explicit form of Td (t) is given by T−1 (t) where ↑ ≡↓, ↓ ≡↑.
In order to obtain the transformed Hamiltonian ĤSW (t) = e i Ŝ (t) ĤHub (t)e −i Ŝ (t) − e i Ŝ (t) (i∂ t e −i Ŝ (t) ) with no charge excitations (i.e.[ D, ĤSW (t)] = 0), let us determine Ŝ (t) order by order in terms of the hopping amplitude.We can write down the first-order term of the transformed Hamiltonian as where the superscript (n) for Ŝ (t) and ĤSW (t) denotes the order of the hopping amplitude, and we have chosen Ŝ (0) (t) = 0.In order to satisfy [ D, Ĥ(1) SW (t)] = 0, S (1) (t) = S (1) +1 (t) + S (1) −1 (t) have to be chosen such that with which the transformed Hamiltonian is given as Ĥ(1) SW (t) = − T0 (t).The solution of the above operator equation for Ŝ (1) ±1 (t) is given as the same form as T±1 (t) but with s i j,±1 (t) instead of t i j e −iA(t)•R i j , where for the time-periodic hopping amplitude.In a similar manner, the second-order term can be obtained as We obtain Eq. ( 6) by inserting Eqs. ( 25), (26), and (29) to the above expression.In the present method, we perform the Schrieffer-Wolff transformation up to the second-order of hopping.This generates three-site term which induces topological superconductivity.Moreover, if terms of fourth order in Schrieffer-Wolff transformation are considered, scalar spin chirality terms scalar spin chirality terms develop, which also induce topological superconductivity [26].We note that the expression for Ŝ (1) ±1 (t) should be slightly modified when the field amplitude A 0 becomes timedependent [44,45], while the correction is negligible when the change of the amplitude is slow enough.The generalized expression for Eq. ( 29) is given as By expanding t (m) i j in the Taylor series around t ′ = t, we obtain The correction term is smaller by a (dimensionless) factor of ∼ | Ȧ0 |a/U, and is indeed negligible for the present calculation.

Gutzwiller approximation
In this section, we describe the detailed derivation of Eqs. ( 12) and ( 13).First, we evaluate the denominator of Eq. ( 11) by site diagonal expectation values, which gives [35,36] where δ = 1 − 1 N i ni↑ + ni↓ is the hole doping rate, f = (1 − δ)/2 is the probability to have a singly-occupied site in the Gutzwiller wave function PG |Ψ BCS ⟩, and f = 1 − f .Here we have used n! ∼ (n/e) n on the second line.We evaluate the expected value of each term in Eq. ( 6) in a similar manner.For the hopping term, we obtain while for the Heisenberg interaction we have The spin part of the three-site term, PG (ĉ As for the density operator, there is arbitrariness in whether to renormalize this as an operator or to replace it as the density at each site.In the present case, in which δ has a finite value, the variational Monte Carlo calculation shows that the densitydensity term has a few effects in contrast to the half-filling case, which allows us to replace PG ni PG with 1 − δ [46].We thus discard the second order fluctuations ( around the expected value of the density in the Hamiltonian here, like PG ni PG ≃ ni δ + 2 f 2 with ni = ni↑ + ni↓ [26].Then the remaining terms in the Hamiltonian can be evaluated as with i k.In this way, we obtain the effective Hamiltonian Eq (12).Evaluation of the first term of Eq. ( 5) incurs additional calculation of the derivative of the variational parameters.Using the relation we can rewrite the time derivative in an operator form as with which the first term of Eq. ( 5) is written as The expectation value PG ĉ † i↑ ĉ j↑ / PG can be evaluated within the Gutzwiller approximation as Note that both doubly-occupied and singly-occupied configurations at site j are allowed as the initial configuration due to the absence of the projection operator.The other expectation value can also be evaluated in the same way as PG ni↑ / PG ≃ ni↑ , with which Eq. ( 13) follows.

Time-periodic Hamiltonian and Floquet Hamiltonian
In this section, we describe how to obtain the time-periodic BdG Hamiltonian and the associated Floquet Hamiltonian (22).The superconducting pairing amplitude of each bond τ can be computed from the gap equation, This implies that the pairing amplitude is transformed as ∆ τ (t) → ∆ ′ τ (t) = ∆ τ (t)e 2iϵt under the unitary transformation ⃗ ψ k → ⃗ ψ ′ k ≡ e iϵtτ z ⃗ ψ k , which keeps the relative phase of the pairing amplitudes intact.When the pairing amplitudes satisfy ∆ τ (t + T ) = ∆ τ (t)e iα , the transformed amplitude becomes time-periodic when −ϵT = α/2.Namely, under this condition the transformed BdG Hamiltonian H ′ (k, t) ≡ e iϵtτ z H(k, t)e −iϵtτ z − ϵτ z becomes time-periodic, and we can define the corresponding Floquet Hamiltonian, where the time evolution of its eigenstates are expressed as the time-periodic wave function multiplied by a plane wave factor.Using the relation for the transformed time evolution operator, U ′ (t, t ′ ) = e iϵtσ z U(t, t ′ )e −iϵt ′ σ z , we can compute the Floquet Hamiltonian as U ′ (t 0 + T, t 0 ) : = e iϵ(t 0 +T )σ z U(t 0 + T, t 0 )e −iϵt 0 σ z .(47) By setting t 0 = 0, we arrive at Eq. ( 22) with ⃗ ψ ′ k (t 0 ) = ⃗ ψ k (t 0 ).Occupation of Floquet band Figures 7(c,d) show the occupation probabilities of the lower Floquet band after the time evolution with the CPL driving.For comparison, we replot the Berry curvature in Figs.5(h,i) as Figs.7(a,b).In the high-frequency case (ω = 5.8t 0 , E = 5.4t 0 /a), the lower Floquet band is almost fully occupied except for the regions where the point nodes were originally located.The lack of occupancy in these regions is attributed to the choice of our initial state in which the energy gap closes.These unoccupied regions should be filled when the scattering processes neglected in the present mean-field treatment are taken into account.In the low-frequency case (ω = 5.8t 0 , E = 5.4t 0 /a), other unoccupied regions arise due to the foldings of copies of the Floquet band.When the external electric field is turned on slowly as in the present case, the occupation of Floquet bands is almost determined by the original distribution function in the absence of the external field.In those cases, the boundaries of the occupied and unoccupied areas of the Floquet bands appear at the position of band crossings arising from band folding in the Floquet Brillouin zone as seen in Fig. 7.In general, the occupation of Floquet bands strongly depends on the time profile of the electric field and presence or absence of scatterings.

FIG. 1 .
FIG. 1. Schematics of Floquet topological superconductivity.Illuminating circularly polarized light to cuprate superconductors gives rise to topological superconductivity of d x 2 −y 2 + id xy -wave pairing from many body effects.

FIG. 4 .
FIG. 4. Driving field dependence of the id xy -wave component of the pairing amplitude.(a) The time average of the id xy -wave component of the pairing amplitude ∆ idxy in the steady state obtained from the time-dependent Gutzwiller simulation.We plot ∆ idxy as a function of the electric field E and the frequency ω in the high-frequency region.(b) The id xy -wave component of the pairing amplitude obtained from the high-frequency expansion ∆ HFE idxy .We plot ∆ HFE idxy as a function E and ω with the formula in Ref. [26].∆ idxy and ∆ HFE idxy show qualitatively similar behavior in the high-frequency region.(c) Comparison of ω dependence of ∆ idxy and ∆ HFE idxy with fixing Ea/ω = 1.In the low-frequency region, the difference between ∆ idxy and ∆ HFE idxy becomes significant, while they show a qualitative agreement in the high-frequency region.∆ idxy shows a sign change around ω = t 0 that is not captured by HFE.(d) (E, ω) dependence of ∆ idxy in the low-frequency region.The sign change appears around E = 1.05t 0 /a, indicating a topological phase transition by increasing the electric field of the CPL.
), we plot the pairing amplitudes with the phase rotation, Re[∆ τ e −iθ ] and Im[∆ τ e −iθ ].Figures 3(f,g) are the blowups of shaded areas of Figs.3(d,e), which corresponds to the approximate steady state under the driving.Here, the d x 2 −y 2 -wave component is given by (Re[∆ [1,0] e −iθ ]−Re[∆ [0,1] e −iθ ])/2 and remains finite in the presence of the driving field.Now let us look at the id xy -wave component of the pairing amplitude.We define the id xy -wave component of superconducting pairing amplitudes as b)  show that the id xy -wave components obtained in the present approach are generally consistent with the HFE in the high-frequency region of ω > 5t 0 .

FIG. 5 .
FIG. 5. Topological properties of Floquet energy bands for a d-wave superconductor without the CPL driving (a,d,g), with the CPL driving of high frequency (b,e,h), and with the CPL driving of low frequency (c,f,i).(a-c) Energy dispersion of Floquet bands with blowups around the point nodes.Without driving point nodes exist in the nodal line direction (a).Those nodes are gapped out when the CPL driving is applied (b,c).(d-f) The phase of gap functions F(k) in the Floquet Hamiltonian.The black curves indicate the k points where the diagonal (normal) component of the Floquet Hamiltonian is zero.Without driving (d), the gap function is real and has nodal lines.The black curve coincides with the Fermi surface in the normal state.In the high-frequency driving case (e), the phase winds twice in the counterclockwise direction along the black curve.In the low frequency driving case (f), the phase winds twice in the clockwise direction along the black curve.Additional black curves arise from the folding of the Floquet bands.(g-i) Berry curvatures of the Floquet bands with and without the CPL driving.The driving frequency ω and the field amplitude E are chosen as ω = 5.8t 0 , E = 5.4t 0 /a for (b,e,h) and ω = 1.03t 0 , E = 0.96t 0 /a for (c,f,i).

FIG. 6 .
FIG. 6. Phase diagram of Floquet topological superconductors with the CPL driving.The color plot of Chern number as a function of (E, ω) (a) in the high-frequency region and (b) in the low-frequency region.The overall feature of the phase diagram is consistent with the sign of the id xy -wave component of the pairing amplitude ∆ idxy .Jumps of the Chern number appear due to the folding of the Floquet bands in the low-frequency region.