Light-induced topological superconductivity via Floquet interaction engineering

Hossein Dehghani ,1,2 Mohammad Hafezi,1,2 and Pouyan Ghaemi3,4 1Joint Quantum Institute, College Park, Maryland 20742, USA 2The Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, Maryland 20742, USA 3Physics Department, City College of the City University of New York, New York, New York 10031, USA 4Physics Program, Graduate Center of the City University of New York, New York, New York 10031, USA

In this paper, we show that the extension of these ideas could lead to creation and manipulation of topological superconducting phases. In particular, we show that optical pumping of electrons in such two-dimensional (2D) semiconductors can generate topologically nontrivial chiral SC [27,28] which hosts topologically protected chiral edge states in the prethermal regime of our driven system. The idea is illustrated in Fig. 1(a), where we apply a circularly polar- ized laser field, in the presence of an external bath to create the population imbalance, required for the development of a nonequilibrium superconducting phase.
The key underlying mechanism for the development of unconventional SC in our system is the following. By varying the pump's frequency, we excite photocarriers of select momentum classes. Due to the optical valley polarization, this leads to an asymmetric occupation distribution around the resonance surfaces in the two valleys, as in Fig. 1(b). This nonequilibrium occupation creates an effective pairing population inversion around one of the valleys [ Fig. 1(c)], which leads to an interband pairing of electrons for a repulsive density-density interaction; that is, the pairing population inversion effectively changes the interaction sign. We should note that although superconductivity due to population inversion has not yet been observed, recent advancements in the development of transient light-induced quantum phases of matter make the realization of such forms of superconductivity much more conceivable.
To study this pairing, the bare density-density interaction should be projected into the band basis, composed of Bloch wave functions. Due to the nonzero Berry curvature of Bloch wave functions around each valley, the effective interaction has a chiral nature and can be decomposed into different angular momentum channels, where each channel has a different dependence on the momentum of electrons. Combined with the fact that the momentum distribution of the excited photocarriers is controlled by the pump's frequency, our setup allows for engineering the dominant form of electron-electron (e-e) interaction. Consequently, we find frequency regimes where a chiral p-wave pairing becomes more favorable than an s-wave pairing. Therefore our results indicate that the periodic drive could be a powerful tool to not only engineer a band but also control the form and strength of the interaction [29][30][31].

II. MODEL
The system considered in this paper consists of a 2D semiconductor with honeycomb lattice structure, such as a single layer of hexagonal boron nitrate (h-BN) [32,33] or transition metal dichalcogenides (TMDs) [34,35]. The electronic band structure consists of two degenerate valleys, and the broken inversion symmetry leads to a gap at two Dirac points K and −K at the corners of the Brillouin zone (BZ), which are labeled by η = (±), respectively. The semiconductor is driven by a laser beam, whose frequency is slightly larger than the semiconductor gap. The Hamiltonian describing the system is H s = H K + H e-e , where H K is a driven kinetic term and H e−e is an e-e interaction. The driven Hamiltonian for the semiconductor has the form where c η † a,k is the electron creation operator of sublattice type a in the vicinity of valley ηK, τ i with i = {x, y, z} is the Pauli matrix acting on the sublattice space in the unit cell, and μ k is the chemical potential, which we assume to be momentum independent. The low-energy Hamiltonian of the two valleys is given by d η k = (vk x , ηvk y , m − κk 2 ), where k denotes the deviation from the K or −K points in the BZ; m, v, and κ correspond to the band gap, the Fermi velocity, and the band curvature, respectively; and k 2 = k 2 x + k 2 y . The optical driving of the system with a circularly polarized laser field is described by minimal coupling (k → k + eA), where the laser field's vector potential is A(t ) = A 0 (cos ωt, sin ωt, 0), with A 0 and ω labeling the amplitude and frequency of the pump, respectively, and we seth = 1. Up to linear order in A 0 the associated Rabi vector of the optical pump is For simplicity, we ignore the effect of the physical spin, which only affects our results when the spin-orbit coupling is comparable to the semiconductor gap [3].
In the following, we denote the corresponding eigenenergies and eigenstates of the undriven Hamiltonian by η α,k and |u η α,k , where the valence and conduction bands are labeled by α = {v, c}, respectively.
For the e-e interaction, we consider a repulsive densitydensity potential U (r − r ), with the corresponding Hamiltonian where ψ † a (r) represents the electronic creation operator with the sublattice index a. To study the possibility of Cooper pairing between electrons, we suppose that the dominant form of the interaction is a screened Coulomb interaction [36]. Therefore, in passing to the momentum space, such interactions are treated as a constant coupling. Denoting the Fourier transform of Coulomb potential by U kk , this implies that U kk = g/N, where g is the interaction strength and N stands for the number of particles in the unit cell. We also restrict our interactions to intravalley scatterings such that in U kk , k, and k belong to the same valley.
To create an effective pairing population inversion, we need a thermal bath. Our bath can have a fermionic or bosonic nature; however, here we only consider a bosonic bath composed of photons or phonons, which is experimentally more feasible, and leave the study of the fermionic bath to Appendix B 2. Specifically, here we assume that our bath can induce relaxation processes between the valence and conduction bands via absorption and emission of photons.

III. MASTER EQUATION
To examine the out-of-equilibrium nature of SC in the presence of a thermal reservoir at temperature T , we use the master equation approach. Assuming that the system-bath coupling is sufficiently weak and the bath has a short autocorrelation time, we employ the Born-Markov approximation to trace out photons from the equations of motion (EOMs). We also consider large pump frequencies compared with Rabi frequencies, which allows us to use the rotating wave approximation (RWA). As a result, we obtain an effective static master equation for the density matrix of the system ρ s [37], 023039-2 where the dissipator operators are L v,k = L † c,k = c † c,k c v,k . Associated with these dissipators we suppose momentumindependent decay rates α corresponding to an effective population n B such that v = n B and c = (1 + n B ), see the Appendices.

IV. ROTATING FRAME TRANSFORMATION
By applying the RWA to the time-dependent term (t ) in Eq. (1), this term in the rotating frame is replaced by the static vectors¯ (+) k and¯ (−) k around the K and −K Dirac points, whose magnitudes are given by |¯ The modification of the e-e interactions becomes transparent through a mode decomposition of the field operator ψ a (r) = α={v,c};η,k 1 √ S u η,a α,k e i(ηK+k).r c η α,k in Eq. (2), where S is the quantization area. The resulting projected Hamiltonian has a contribution in the Cooper channel for the interband pairing as [38] H e-e = k,k ,η=±Ū where the projected density-density interaction We see that the Bloch wave functions which encompass the topological characteristics of the system control the form of e-e interactions. The crucial effect of the Berry curvature of the band structure on the e-e interactions is embedded in the Bloch wave function overlaps inŪ kk which after inserting for the eigenstates can be decoupled in three channels according to their angular mo- , and d k = |d k |, see the Appendices. Correspondingly, for momenta close to the corners of BZ, for small k and k this givesŪ kk g[1 + F 4 (2k.k − 2ẑ.k × k − k 2 − k 2 )]/N, where F = v 2 /m 2 is the Berry curvature at the Dirac point. Recently, such topologically induced modifications of the e-e interaction have been associated with modification of the excitons' spectrum [39,40].

V. MEAN-FIELD ANALYSIS
To study the possibility of the Cooper pair condensation, we use a mean-field (MF) approximation for the e-e interaction and express it as Notice that we introduce two pairing order parameters (±) k , depending on whether the valence electrons around the K valley and the conduction electrons around the −K valley are bound to each other or vice versa. Employing the MF expression above, we can use Eq. (3) to write a closed set of equations for the occupation numbers n η α,k = tr(ρ s c η † α,k c η α,k ), the polarization p η k = tr(ρ s c η † c,k c η v,k ), and the anomalous pair- . This approach leads to legitimate results at the onset of the SC phase transition, where the distinction between the Bogoliubov quasiparticles and electrons is negligible. The EOMs for n η α,k and p η k in the absence of pairing are familiar and known as the optical Bloch equations in the literature [41,42], and we leave their derivation to Appendix B. Here, we only present the EOM of the anomalous pairing, which is less familiar, where we define the total decay rate as t = v + c and the total kinetic energy as t,k = v,k + c,k . We note that on the right-hand side of this equation, the two occupation probabilities in the parentheses belong to two different valleys, which is a manifestation of the Cooper pairing. In the steady state, where the right-hand side of Eq. (B34d) vanishes, s η k satisfies s η k = −i η k n η sc,k /(i t,k + t /2), where we define the interband pairing population as n η sc,k = 1 − n η v,k − n −η c,−k . Since n η v,k and n η c,−k can be independently populated due to the optical valley selection rules, under nonequilibrium conditions the pairing population can acquire a finite value.
Using the MF definition of the anomalous pairing s k , the self-consistency equation at the onset of the phase transition becomes η k − k Ū kk s η k . However, since in dissipative systems, s k may have a finite imaginary part, this equation should be modified by a more accurate Keldysh approach [3,25]. For weak decay rates, we get where the real part in brackets indicates the dissipative suppression of the pairing and ensures that η k remains real. Correspondingly, in the steady state this equation gives where we remark that n η sc,k effectively determines the interaction sign [3,25]. The asymmetry of the Rabi frequencies at the two valleys results in the corresponding steady-state occupation probabilities differing significantly as depicted in Fig. 1(b). For the polarization we have chosen, this leads to a positive (negative) value for n (+) sc,k (n (−) sc,k ) as is displayed in Fig. 1(c). Thus, with a repulsive interaction, we can have a SC instability by developing a nonvanishing value for the order parameter (−) k , while (+) k remains vanishing. Thus, in what follows we drop the valley superscript in k for the nonvanishing order parameter.
We also highlight that the form ofŪ kk crucially determines the form of the resulting gap. In fact, the self-consistency equation above can be solved using a simple ansatz for the gap function of the form where l = {0, 1, 2} should be ascribed to the angular momentum of the s-, p-, and d-wave pairing modes and f (l ) 's play the role of the SC form factors. Using this ansatz and inserting t,k = 2μ, the linearized gap equations for the three different types of pairing decouple, and the critical coupling strength for each channel reads where Since we only consider direct optical transitions, ζ η k is essentially the pairing population in each momentum class, in the weak-drive limit.
As illustrated in Fig. 1(c), the main contribution to n sc,k comes from around the resonance energy ring, denoted by k r , where the detuning frequency vanishes, d,k r = 0. Hence we consider a given UV cutoff E around this surface and show that the resulting phenomena are independent of the exact value of this parameter.
The frequency of the pump determines which momentum classes are resonantly excited. Equation (7) indicates that only the states with negative n sc,k can form pairing. Moreover, since the projected Coulomb interaction has momentum dependence, the critical value of the coupling strength depends on the pump frequency. This behavior is depicted in Fig. 2(a), where we notice that the preferred form of pairing transforms from s wave to p wave as the drive frequency increases. This transition is associated with the momentum dependence of the SC form factors. Since n sc,k is peaked around the resonant surface, we only need to consider the momentum dependence of the form factors around this surface. Consequently, for a given frequency the radius of the resonance ring k r (ω) is obtained, which can be inserted to determine the form factors f (l ) [k r (ω)]. These functions are displayed in the inset of Fig. 2(a). Here, we observe a similar behavior to that in the main plot of Fig. 2(a), where by increasing the frequency the initially dominant s-wave form factor becomes subdominant with respect to the p-wave form factor.
Other than the pump's frequency, critical coupling depends on the amplitude of the pump, too. This is displayed in Fig. 2(b), where the horizontal axis has been chosen to be the dimensionless parameter evA 0 / , which appears in the gap equation through ζ k . We notice that the critical coupling strength of all the three SC modes always decreases as the pump amplitude increases. Specifically, in the lowintensity limit (evA 0 ), the critical coupling is inversely proportional to the intensity ∝ ( evA 0 ) 2 , which could be associated with the fact that the peak value of the excited population (|n (±) sc,k | ζ (+) k ) increases linearly with the intensity. At higher intensities, where evA 0 , this behavior changes, since the width of n (±) sc,k , i.e., the number of momentum classes participating in pairing, keeps increasing, and therefore we do not observe a saturating behavior. Furthermore, in the inset of Fig. 2(b), we depict g crit as a function of the energy cutoff, which shows that once E becomes comparable to the band gap, the cutoff dependence of g −1 crit becomes insignificant. Finally, note that g crit can be controlled directly by the chemical potential and its minimum value is reached when μ = t /2. This feature provides a high tunability in choosing the other parameters of our system, see the Appendices.

VI. SIGNATURES OF TOPOLOGICAL SC
Given the diversity of proposals to realize topological phases in driven-dissipative systems, it is crucial to present the experimental signature of the light-induced topological phase in our proposal. Here, we show that the interband p-wave pairing hosts edge supercurrents which are experimentally detectable. To analyze the edge modes, the SC behavior of the system can be described by an effective MF Hamiltonian which includes a kinetic energy and a SC pairing between the valence (conduction) electrons at K (−K). Due to the nonequilibrium nature of superconductivity in our system, only electrons around the resonance ring participate in SC. Hence the steady-state Hamiltonian corresponds to electrons with momentum close to the resonance surface with momenta k ≈ k r + δk, where δk k r denotes the momentum deviation from the resonance surface. Thus the effective MF Hamiltonian should be projected into these states. This structure is comparable to conventional superconductors, where electrons at the Fermi energy control the properties of the superconductor and edge modes [43,44]. To examine the existence of localized states at a hard boundary parallel to the y direction, we set k y = 0 and replace the momentum deviation from the resonant point, δk x , with −i∂ x . The resulting Hamiltonian becomes where the SC order has a p-wave structure ( k r = 0 k r ). The effective Hamiltonian commutes with σ z τ z , and consequently, the eigenstates are in two independent sectors corresponding to σ z τ z = ±1. We put the boundary at x = 0, and the sample is in the x < 0 region. The two chiral states resulting from momenta close to each resonant momentum have energies E ± = ± eA 0 v 2 and are separated with a gap of roughly 0 from the continuum states close to the resonant momenta. Their corresponding wave functions are where λ = m 2 Spinors |s σ z (|s τ z ) correspond to eigenspinors of σ z (τ z ) Pauli matrices with eigenvalue s. The ± sign corresponds to the two sectors with σ z τ z = ±1. It is noteworthy that these edge states have opposite chirality in the two sectors as shown in Fig. 1(a) and can solely emerge when the pairing has a p-wave structure. Furthermore, these states carry supercurrents with no potential drop, which distinguishes them from normal edge currents and was experimentally detected before in graphene [14]. To verify the development of the edge states outlined above, we can use the method used before to distinguish the edge and bulk supercurrents in the steady-state superconducting phase [45]. To this end, we consider a TMD sample where part of the sample (including part of the edge) is illuminated by an appropriate laser beam. We then measure the supercurrent carried through the edge which is illuminated by the laser field. It is also shown that in a magnetic field, fluxoid quantization generates a periodic modulation of the edge condensate which is observable as a fast-mode oscillation of the critical current versus the magnetic field. It should be noted that the fast-mode frequency is distinct from the conventional Fraunhofer oscillation displayed by the bulk supercurrent and the frequency of such oscillations should increase with the superconducting area [45]. Interestingly, in our setup the superconducting area can be readily increased by an increase in the laser beamwidth, which makes such measurements convenient for our setup.
The edge states' supercurrent is contrasted with the bulk states' current through its dependence on the size of the superconducting region and an external magnetic field. Such measurements have been realized through experiments on samples with different sizes [45], which in our setup is simply controlled by the width of the pump beam.

VII. EXPERIMENTAL FEASIBILITY
Finally, we provide an estimate for the pump's amplitude based on typical energy scales in 2D two-valley semiconductors. Specifically, to verify the feasibility of the realization p-wave SC in our model, we need to estimate the required critical coupling constant. Promising 2D semiconductors to realize the phenomena outlined here are h-BN or TMDs with a band gap of the order of 5 eV [46,47]. In these materials the screened Coulomb interaction g typically is comparable to the band gap of these materials [36] or can be enhanced via Coulomb engineering [48]. From Fig. 2(b), we deduce that to obtain log 10 (g crit /m) 0, one requires an electric field such that the ratio evA 0 / would be of the order of 10 2 . Since typically the inverse decay rate is of the order of picoseconds [49], this implies that the required Rabi frequency for our proposal should be 10 THz. For a typical semiconductor, this Rabi frequency corresponds to an electric field of 5 × 10 6 V/m. Recently, strong fields with an electric field of 4 × 10 7 V/m have been used in generating a light-induced Hall effect in graphene [14]; therefore we believe that our proposal is within experimental reach.

VIII. OUTLOOK
The Floquet engineering of the interaction described here is a versatile effect which can be generalized to other lattice symmetries where the band structure has a nontrivial topological structure. Correspondingly, interesting directions to explore in the future are the study of similar effects in multilayer twisted semiconductors and the generation of other interacting topological phases. While here it suffices to study the prethermal regime of the system, which at the high frequencies considered in our proposal can last for many periods of the drive [50,51], for longer time scales scattering processes with acoustic phonons, which mediate momentum and tend to thermalize the system, become important. It would be fascinating to explore under what circumstances the superconducting state proposed here survives in the presence of such processes [52,53]. Finally, for large decay rates a mean-field approach is not satisfactory, and the effects of the fluctuations should be investigated [54][55][56].
where we have dropped the valley index in d k and d z,k , which have the same form in the two valleys and are defined as d η ±,k = d η k,x ± id η k,y . To apply the rotating wave approximation to a nondiagonal Hamiltonian, we need to first transform the Hamiltonian into the energy basis where it is diagonal and then apply a time-dependent rotation to the two energy levels so that the time dependence of the two transformed eigenstates becomes approximately the same. The first transformation is done through a similarity transformation by the unitary matrix where we have used the eigenstates of the undriven Hamiltonian, and the second transformation is realized by the diagonal time-dependent transformation diag(e iωt/2 , e −iωt/2 ). The combination of these two transformations is R η k (t ) = (|u η v,k e iωt/2 |u η c,k e −iωt/2 ). Therefore, denoting the electronic spinors in the laboratory frame and the rotating frame as T We apply this transformation to all of the terms in the Hamiltonian system. While the undriven Hamiltonian is trivially diagonalized, the pump Hamiltonian should be obtained after averaging over time. To evaluate the temporal average of the drive term, it would be convenient to decompose the time-dependent terms as where . Correspondingly, the expression that must be averaged over time is The final result of this calculation becomes Similarly, we need to transform the electron-electron interaction term by rotating the electronic creation and annihilation operators and averaging over time. The final result of this calculation, in addition to the right-hand side of Eq. (8) in the main text, has other contributions which include the overlap of the valence and conduction wave functions at close momenta which makes such terms negligible.
Here, we mention that in order to integrate the gap equation, we consider an energy cutoff with respect to the resonance surface. The resonance ring in the BZ is defined by ω = 2d k=k r , demanding that where k r is the radius of the resonance ring. We can also use the above equation to define the integral bounds of the radial momentum through the UV energy cutoff E as follows: Furthermore, we note that we have used this relation in the main text to define the form factors as a function of the frequency f (l ) (ω). In particular, to determine the frequency where a transition from the s-wave SC to p-wave SC occurs, we should satisfy where From here, we can observe that in order to satisfy this condition, it is desirable to have a positive band curvature κ so that d z,k and correspondingly f (0) k would decrease with the momentum. Therefore, since f (1) k increases monotonically with the momentum, this condition can be satisfied with smaller values of the momentum deviation from the valley center.
Regarding the numerical parameters chosen in the main text, we should mention that in Fig. 2 the magnitude of v is chosen such that it results in a large value for k r 2/a, which may not be accessible in lattice models. A more realistic parameter regime is obtained by increasing v, which reduces the required values of k r and still captures the same competition between the p-wave and s-wave superconducting states but now with a larger value for the critical coupling log 10 (g crit /m) 1 which has been numerically verified in our simulations. This issue can also be resolved, and it does not affect the feasibility results of our paper since the magnitude of the critical coupling strength can be decreased by choosing a smaller value for the chemical potential closer to the optimum value of the chemical potential, i.e., μ = t /2 as mentioned in the main text.
Finally, having a finite deviation from the Dirac points has an additional dynamical effect which helps our proposal indirectly. We note that in our proposal, to develop a superconducting state, it is essential to have a large Coulomb interaction which can occur as a result of screening the Coulomb interaction. However, this screening can also lead to large exciton binding energy in TMDs. Therefore it might appear that exciton formation could compete with the formation of SC. This issue can be avoided by a strong optical pumping of the system well above the band gap. More precisely, for the formation of excitons, the excited electrons should relax to the bottom of the conduction band. The effective energy associated with inverse relaxation rate through optical phonons is 2 meV [57], which is two orders of magnitude smaller than the Rabi frequency being more than 100 meV. Correspondingly, in our derivation, we ignore the exciton formation processes.

Bosonic bath
In this section, we obtain the equations of motion (EOMs) in the presence of a bosonic bath. In general, the bath used in our formalism and its coupling to our system can be described by the following Hamiltonian: Using such a bath leads to a master equation where the action of the Lindbladian superoperator L with a quantum jump operator L on the density matrix ρ is defined as While in general the decay rates depend on the coupling constant λ k and the density of states of the bath, to simplify our formalism, we consider constant decay rates, v = n B and c = (1 + n B ), where n B denotes the effective Bose-Einstein population of the bath, which we assume is momentum independent.
Here, we are interested in the EOMs for the occupation probability n η α,k , the polarization σ η k , and the anomalous pairing s η k = tr(ρ s c −η c,−k c η v,k ). For convenience of labeling we define the notation n η αα,k ≡ n η α,k and n η cv,k ≡ σ η k . To derive the EOMs for an arbitrary operator O = tr(ρ sÔ ) in the Schrodinger equation, we use ∂ t O = tr(Ô∂ t ρ). To write down the EOMs for these quantities, we also need the following We can study the contributions of the Hamiltonian and Lindbladian in the time evolution separately. For the kinetic Hamiltonian part with H K = c † α,k h αβ,k c η β,k we can use the following identities: Similarly, for the interband pairing we get We can also compute the commutators of the anomalous pairing and the electron-electron interaction, which is The corresponding commutator becomes For the Lindbladian contributions we employ Eq. (B5) to obtain After combining the Hamiltonian and Lindbladian contributions, we get Notice that should we want to take exciton formation into account, in Eq. (B11c), which is the EOM for the polarization, we need to consider the Hartree-Fock contribution of the electron-electron interaction in the particle-hole channel. This adds a term −i η k Ū kk σ η k on the right-hand side of this equation. We can show that in our system, exciton formation and the Cooper instability do not compete with each other. Therefore, even in the presence of a finite density of excitons, we can still have a phase transition into a superconducting state. Thus, in our derivation, we drop such terms to simplify our analysis. From Eqs. (B11c) and (B11d) we can obtain the anomalous pairing and the polarization These relations can be inserted in the first two EOMs in the steady state, where we have defined The equations at the two valleys should be solved together. This gives where the effective Rabi frequency and pairing amplitude are given by respectively. The resulting equations can be rewritten in a matrix form, ⎛ Here, we are mainly interested in studying the onset of the SC phase transition, which implies that we can ignore the pairing amplitude in the above equations so that the matrix on the left becomes block diagonal. In the limit where the effective Rabi frequency around the −K valley, i.e.,¯ (−) k , is negligible, after using the conservation of the particle densities in the valence and conduction bands of the two valleys separately (n (η) v,k + n (η) c,k = 1), these probability populations become Let us further assume that n B = 0, which results in γ v = 0 and γ c = 1. In this limit, it is evident that we can have an effective SC population inversion around one of the valleys because We should hint that in the weak-drive limit, the right-hand side reduces to +ζ (+) k and −ζ (+) k . More generally, after defining the interband pairing population, n (η) sc,k ≡ 1 − n (η) v,k − n (−η) c,−k , we have where we have used the fact that at every momentum we have n (+) v,k + n (+) c,k = 1. These equations lead to the linearized gap equation in the main text.
Finally, in deriving the final form of the gap equation from the mean-field solution, we note that in dissipative superfluid or superconducting systems, in general it is possible that the condensate attains a time dependence [55]. Let us decompose the total Hamiltonian into its system and system-bath components H = H s + H s-b , where the former has kinetic and interaction contributions as H s = H k + H int . After integrating out the reservoir degrees of freedom the effective Hamiltonian that is obtained for the system's degrees of freedom is quadratic and non-Hermitian. Therefore the total quadratic Hamiltonian that is obtained from summing this contribution and H k would be non-Hermitian, too. Consequently, if we try to use the corresponding non-Hermitian free energy to find saddle-point solutions, it would lead to inconsistency as the contributions of the quadratic terms are non-Hermitian while those of the interaction terms are Hermitian. This problem is solved by using the Keldysh action, which has a forward and backward temporal contour such that by differentiating between the retarded, advanced, and Keldysh Green's functions, the Hermiticity of the action is always built into it. The detailed Keldysh calculation is explained in Ref. [3]. In the limit where the dissipation rate is small, this result agrees with the modified mean-field approximation.

Fermionic bath
Here, we show that we can obtain similar results with a fermionic bath at a fixed temperature T , where α = {v, c}. We consider a system-bath coupling which allows the exchange of particles between the system and the reservoir, Starting with the system-bath coupling term, we assume a thermal Fermi-Dirac distribution for the bath degrees of freedom (DOFs) at temperature T , so that these DOFs can be traced out. After applying the RWA and eliminating the oscillating terms, we arrive at the following master equation for the density matrix of the driven semiconductor: where n F α,k = n F ( α,k ) is the Fermi-Dirac distribution. The decay rates α (k) = 2π a |t a | 2 ν( α,k )u α * a,k u α a,k , where ν( ) represents the density of states of the bath's electrons at energy .
For the pairing amplitude between the valence α = v and conduction band α = c, this yields For the Lindbladian part we get Ô , c η β,k ρ , (B32) where we have used the creation and annihilation operators in the rotating frame. Consequently, we can assume that oscillating terms in the rotating frame can be ignored. This way, we can time average over the Lindbladian, which results in considering only the diagonal terms in the above with α = β.
Without loss of generality, in the rest of this section, we assume momentum-independent dissipation rates, and we label its diagonal components as α . The terms obtained from expanding the right-hand side are similar to the terms obtained in the bosonic case. The final result of this expansion reads where as before we have ignored the terms which are relevant in exciton formation. Next, we derive the steady-state solution by assuming constant densities and pairing amplitudes in the rotating frame. We start by obtaining the equations for the anomalous pairing, where we have defined t,k = c,k + v,k and t = c + v . In the next step, we consider the EOM for the polarization σ k , where we have defined d,k = c,k − v,k − ω. Inserting these two equations in the occupation probabilities, we get 023039-9 where the effective Rabi frequency and the effective pairing amplitudes are (B38) As in the bosonic case, we need to solve four equations simultaneously, where t = v + c , γ α = α / . We can rewrite these equations in a matrix form, ⎛ In general, one needs to invert the matrix on the left to find the solutions for the occupation probabilities. As the first step, we consider the linearized gap equation where we only consider the solutions of the above equation in the zeroth order of k . Furthermore, we consider the zero-temperature limit where n F,v/c k = 0, 1, where (−) −k = 0 for k around the K Dirac cone. This yields We can further simplify these relations in the limit that the Rabi frequency around the K point is negligible, The above relations can be employed for the interband pairing, which can be used to derive the gap equation. To perform this task, we need to write the self-consistency definition of meanfield order parameter. The result of this calculation yields where we have used the definition n η sc,k = 1 − n η v,k − n −η c,k . As in the bosonic bath case, we can see that this equation can only be satisfied around one of the valleys, which for our choice of the laser's polarization will be the K valley. Therefore we can drop the valley index and rewrite this equation as Since the only difference between this gap equation and the gap equation in the bosonic bath case is in the effective value of n η sc,k , we can use the same ansatz for the pairing amplitude as before, Using this ansatz, we can evaluate the critical value of the coupling constant g numerically. After employing the same integration method, we obtain a similar behavior for g crit as a function of the frequency of the pump, and we observe that a transition from a s-wave SC pairing to a p-wave pairing is possible. This shows that the phenomenon we observe is due to the specific form of the electron-electron interaction that we engineer and independent of the type of bath that we use in our model.

APPENDIX C: CHIRAL EDGE STATES
Here, we investigate the possibility of having protected edge states which can carry supercurrent for our system [58]. In particular, we study the possibility of having localized modes in the presence of a hard-wall boundary which is parallel toŷ which requires the wave function to vanish at x = 0.
We assume that we are in the regime where the SC pairing order parameter has acquired a significant value and is no longer negligible as we previously imagined in the course of obtaining the dominant form of the SC order parameter. Besides, we recall that the main effect of the dissipation in our system is to allow the formation of out-of-equilibrium steady states where pairing is possible for electrons around the resonant ring corresponding to momentum k ≈ k r + δk, where |δk| < |k r |. Similar to our calculations in the main text, this goal is achievable with a relatively small value of the system-bath coupling. Therefore, for our purpose, which is to study the topological properties of the system after reaching this state, we can ignore the system-bath coupling and only consider the Hermitian terms of the effective Hamiltonian in our system.
In the rotating frame the kinetic energies of the valence and conduction bands around the two valleys are where on the right we set α = ±1, which corresponds to the conduction and valence band energies, respectively. Correspondingly, the second-quantized form of the kinetic term reads H K = α,k α,k c η † α,k c η α,k . Up to quadratic order, d k ≈ m + k 2 2m (v 2 − 2mκ ). For states with momentum close to resonant momentum k r = 1 2 ω 2 −4m 2 v 2 −2mκ , we get d k r +δk ≈ ω 2 + k r ·δk m (v 2 − 2mκ ). Us-ing this approximation, the kinetic energies read η α={v,c},k = μ + αṽ The light-induced modification of the band structure is implemented through the Rabi vectors in the two valleys η k as found in Eqs. (A3) and (A4). To simplify our study, we consider k y = 0. The resonance surface is reduced to the resonance points at k r = ±k rx . Without loss of generality we also set the band curvature to zero (κ = 0). The magnitude of the Rabi vectors then reads as where 0 = eA 0 v. We note that after thermalization, the pairing obtains a finite value in the basis above where the kinetic energies are diagonalized. Thus, as the next step, we add the Bogoliubovde Gennes (BdG) pairing Hamiltonian of the system to the Hamiltonian above.
Next, we transform the pairing Hamiltonian of the system in the momentum space in the mean-field limit where (−) k ≡ k is finite and (+) k is vanishing. The pairing Hamiltonian is where the superconducting order parameter k = 0 (k x ± ik y ) corresponds to a chiral p wave. As outlined above, we initially consider k y = 0.
(b)−1 sector. In the other sector we project the Hamiltonian to the subspace where τ x σ x = −1. This subspace is spanned by the following states: |3 = |σ z = 1, τ z = −1 and |4 = |σ z = −1, τ z = 1 , where as before we use ξ i 's to denote the Pauli matrices in this basis. Consequently, the Hamiltonian required for the Volovik approach becomes which should be compared with Eq. (C13). As before we apply a phase shift to the wave function as |ψ → |ψ = e −iγ kr x |ψ . The required phase is given by Similar to the +1 sector, we can find a localized state around x = 0, by superposing the localized states in the vicinity of k r and −k r . The final resulting state is |ψ = sin(γ k r x)e λ kr x |ξ y = 1 = sin(γ k r x)e λ kr x (|σ z = 1, τ z = −1 + i|σ z = −1, τ z = 1 ), (C20) whose energy is − 0 /2. The analysis above demonstrates that in the absence of the matrix τ x σ 0 in the Hamiltonian H eff , we have two chiral localized states on the edge with energies ± 0 . To show that the states are chiral with opposite chirality, we need to consider small momentum k y parallel to the x = 0 edge. Such small momentum corresponds to the addition of a term of the form 0 k y τ x σ y . One can readily see that the edge states |ψ and |ψ are eigenstates of k y τ x σ y with eigenvalue ±1, correspondingly. As a result, they disperse linearly with parallel momentum with opposite velocities and form oppositely chiral supercurrents. Now the presence of pairing elements associated with τ x σ 0 in general can mix the corresponding eigenstates in the two sectors. This can be easily shown by evaluating the matrix elements of τ x σ 0 between the localized states and continuum states in the two sectors. For large system sizes due to the exponential decay of the localized states at the boundaries, the hybridization of the localized states from different sectors is negligible. However, we should still consider the possibility of the hybridization of the localized states of one sector with the bulk states of the other sectors. To answer this question, let us first find the functional form of the bulk states. Here, we only consider the +1 sector, and the bulk states of the other sector can be obtained in a similar manner. The eigenvalues are conveniently obtained by diagonalizing the corresponding Hamiltonian, which (for the two sectors) gives The associated eigenstates are where N = [ δk x + √ v 2 m 2 k r 0 ] 2 + 2 k is the normalization factor. As discussed above, τ x σ 0 mixes the two sectors (both for localized and continuum states). Since γ and γ for nonzero μ and k r are independent, the matrix elements between states from the two sectors (localized or continuum) vanishes. As a result, the projection of the τ x σ 0 term into the states close to the resonance vanishes, and these terms do not contribute to the effective Hamiltonian close to resonant momenta.
However, as mentioned previously, these matrix elements are negligible because in order to match the traveling component of the corresponding localized states and bulk states with different momenta, we need to incorporate states with momenta largely different from the resonant momentum. Consequently, since such states do not contribute to superconductivity, their corresponding matrix elements can be ignored.