Light-induced switching between singlet and triplet superconducting states

While the search for topological triplet-pairing superconductivity has remained a challenge, recent developments in optically stabilizing metastable superconducting states suggest a new route to realizing this elusive phase. Here, we devise a testable theory of competing superconducting orders that permits ultrafast switching to an opposite-parity superconducting phase in centrosymmetric crystals with strong spin-orbit coupling. Using both microscopic and phenomenological models, we show that dynamical inversion symmetry breaking with a tailored light pulse can induce odd-parity (spin triplet) order parameter oscillations in a conventional even-parity (spin singlet) superconductor, which when driven strongly can send the system to a competing minimum in its free energy landscape. Our results provide new guiding principles for engineering unconventional electronic phases using light, suggesting a fundamentally non-equilibrium route toward realizing topological superconductivity.


Introduction
Topological superconductors are elusive unconventional superconducting phases [1][2][3] that can host topologicallyprotected Majorana boundary modes and non-Abelian vortex excitations [4,5], which are of fundamental as well as tremendous practical interest as a route towards fault-tolerant quantum computing [6].Spin-triplet superconductors with finite angular momentum Cooper pairs [7,8] have long been regarded as particularly promising candidates, with degeneracies between nodal order parameters expected to favor a chiral topological superconducting state [9].However, spin-triplet pairing remains rare in nature and signatures of chiral topological order remain inconclusive, despite several candidate compounds such as Sr 2 RuO 4 [10][11][12] or UTe 2 [13,14] being placed under exceptional experimental scrutiny.
At the same time, a series of pioneering pump-probe experiments have established irradiation with light as an alternative and fundamentally non-equilibrium tool for interrogating and manipulating superconducting phases on ultrafast time scales, ranging from time-resolved probes of Higgs [15][16][17][18] and Leggett [19][20][21][22] mode oscillations in conventional and multi-gap superconductors to the light-induced enhancement or induction of long-lived superconducting signatures in the fullerides [23][24][25].With the underlying mechanisms still under substantial debate, these observations coincide with broader experimental [26][27][28][29] and theoretical efforts [30][31][32][33][34] in exploring thermal and non-thermal pathways to suppress or control competing ordered phases with light.
These results immediately raise the tantalizing question of whether elusive topological spin-triplet superconducting states can instead emerge as metastable phases upon irradiating a conventional superconductor with light.Consider an inversion-symmetric material with a conventional s-wave superconducting phase that preempts a closely competing topo-ence, the material can refuse to relax back to its equilibrium phase, instead becoming trapped in a metastable spin-triplet superconducting phase, with thermalization to the equilibrium spin-singlet phase suppressed due to restored inversion symmetry after the pulse.
In this work, we illustrate this mechanism as both a new route towards engineering spin-triplet and topological superconducting phases out of equilibrium as well as a generic probe of competing odd-parity instabilities in conventional superconductors.The protocol is summarized in Figure 1.We first identify ultrafast and two-color pulses as two complementary routes to dynamically break inversion symmetry, and demonstrate using a minimal model of quasiparticle dynamics that they can conspire with spin-orbit coupling to transiently induce odd-parity order parameter oscillations in a system with even-parity order.We then derive an effective time-dependent Ginzburg-Landau theory dictated via symmetry and present an explicit switching protocol for driving the system to settle into a metastable odd-parity superconducting state.Remarkably, we find that the coupling between equilibrium conventional s-wave order and a competing spintriplet order parameter necessarily scales linearly with the field strength, in contrast to ordinary Higgs mode excitations which scale quadratically with the field.We illustrate the proposed mechanism for light-induced switching to a triplet superconductor by example of two different lattice models, one of which features a metastable chiral topological superconducting state, and discuss implications for real systems such as dilutely-doped 1T' WTe 2 .Our results reveal new guiding principles for engineering metastable unconventional superconducting states using light.

Dynamical inversion symmetry breaking with light
While conventional centrosymmetric superconductors typically comprise local spin-singlet Cooper pairs with even parity and net spin S = 0, spin-triplet superconductivity entails odd-parity pairs with net spin S = 1 and three allowed values m s = 0, ±1 for the spin component.Spin and parity are linked by necessity; since the combination of inversion and spin-exchange has the net effect of exchanging two fermions, even-parity pairs must be singlets while odd-parity pairs must be triplets.Coupling singlet and triplet orders therefore immediately requires the breaking of both inversion and SU (2) spin rotation symmetries.The latter is readily broken in superconductors with strong spin-orbit coupling such as heavyfermion compounds.Most notably, monolayer 1T ′ WTe 2 has been observed to host both a quantum spin hall insulating phase [40][41][42][43] due to strong spin-orbit coupling and a proximal superconducting phase for dilute doping [42] with a rich array of predicted competing superconducting states with different pairing symmetries [44], rendering it a prime candidate to search for an out-of-equilibrium topological superconducting phase.
In contrast, inversion-breaking turns out to be more subtle, requiring careful consideration of the optical driving scheme.
A simple monochromatic light wave is not sufficient, since it preserves a more general dynamical inversion symmetry defined by a combination of parity and time translation by half the wave period T , This dynamical symmetry can be strongly broken in two ways: (a) via the envelope of an ultrafast pulse, or (b) via a two-color pulse E(t) = E 1 cos(ω 1 t) + E 2 cos(ω 2 t) if the two constituent frequencies are not odd harmonics ω 1 = pω 0 , ω 2 = qω 0 of a common frequency ω 0 [45][46][47].Dynamical symmetries have been extensively studied as a means of high-harmonic generation in atomic and molecular systems [48][49][50] and are more recently being used for optically controlling solids [51][52][53][54].
We first illustrate the ramifications of breaking these two symmetries for a minimal mean-field model of a centrosymmetric honeycomb lattice superconductor with effective attractive local (U ) and nearest-neighbor (U ′ ) interactions.Importantly, the inclusion of spin-orbit coupling λ via spin-dependent next-nearest-neighbor hopping with phases ν ij = ±1 for left or right turns [55] reduces spin rotation symmetry to U (1), permitting a coupling between singlet and m s = 0 triplet pairs ∼ (ĉ −k↑ ĉk↓ ∓ ĉ−k↓ ĉk↑ )/ √ 2. A standard BCS mean-field decoupling of the interaction in the Cooper channel introduces the superconducting gap function ∆ αβ (k) = i f i αβ (k)∆ i with sublattice indices α, β which can be decomposed into pairing channels classified in terms of the irreducible representations of the crystal point group with form factors Suppose now that a conventional s-wave superconducting phase in equilibrium is irradiated with a weak but wide pump pulse, which couples to electrons via the Peierls substitution with a vector potential A(t) polarized in the x direction.The pulse is parameterized via the dimensionless field strength |A| = ea 0 E 0 /ℏω where E 0 , ω and a 0 denote the electric field amplitude, frequency, and the lattice constant, respectively.For numerical expediency, in Fig 2 we use effective interactions U = U ′ = −2 and β = 10, in units of hopping.Strikingly, the onset of p x -wave order parameter oscillations for weak light pulses is linear in the field strength, shown in Fig. 2a-b for a two-color pulse, and in stark contrast to ∼ A 2 scaling for ordinary amplitude mode oscillations.Furthermore, the amplitude of p x order scales linearly with λ, completely vanishing in the SU(2) symmetric limit λ = 0. To illustrate and linear-in-λ (spin-orbit parameter) coupling between an equilibrium order parameter ∆s and a px-wave order parameter ∆p x .A two-color x-polarized pulse that breaks dynamical inversion symmetry is applied with various small magnitudes |A|, and example time-plots of the resulting fluctuations in |∆p x (t)| are plotted in b.The maxima of these fluctuations are plotted versus |A| in a, for five different magnitudes of the spin-orbit coupling strength λ in units of hopping.c Time dynamics resulting from a two-color pulse that breaks dynamical inversion symmetry (with a 2 : 1 frequency ratio).d Time dynamics resulting from a two-color pulse that preserves dynamical inversion symmetry (with a 3 : 1 frequency ratio).Both c and d start from purely even superconducting order (indicated by blue curves), and the onset of odd orders (indicated by red and pink curves) is suppressed when dynamical inversion symmetry is unbroken (see legend in bottom right).e Short-time deviations in the s-wave and px-wave order parameters for weak driving with x-polarized light as a function of frequency Ω, revealing a BS mode for the px-wave order parameter that is sub-gap (indicated by the unshaded region).
the role of inversion symmetry breaking, Fig. 2c and d depict the order parameter response to a two-color pulse with 2 : 1 and 3 : 1 frequency ratio, respectively.A 3 : 1 frequency ratio preserves a dynamical inversion symmetry [Eq.(1)] with time translation t → t + π/ω; p-wave order oscillations are consequently dramatically suppressed at short times.Conversely, broken dynamical inversion symmetry efficiently excites oddparity order already for short times [Fig.2c].Note that, due to the nonlinear nature of the quasiparticle equations of motion, these heuristics only apply for time scales of only a few cycles of the pump pulse.Central to efficient switching to triplet order, Fig. 2e reveals a BS mode resonance for subdominant pwave order parameter oscillations that crucially lies below the ordinary Higgs mode and pair breaking excitations.Spectroscopic observation of this mode would immediately provide an experimentally accessible handle to probe the existence of a subdominant pairing channel, and would importantly suggest that stronger ultrafast excitation of this mode can potentially nudge the system well beyond linear order parameter oscillations and into a metastable competing phase with triplet pairing.
Optical switching to a metastable state Insight into whether strongly driving a triplet BS mode can allow for light-induced switching to a metastable odd-parity superconductor can be readily gleaned from an effective timedependent Ginzburg-Landau (TDGL) description, which en-codes the coupling of multiple order parameters to light and importantly accounts for relaxation.In this picture, a suitably tailored pulse liberates the superconducting order parameter from its global free energy minimum and brings it close enough to a proximal local minimum that it relaxes into a metastable opposite-parity phase.A minimal Lagrangian that describes this process reads and includes a kinetic contribution with damping coefficients ij and inertial coefficients Γ (2) ij .The equilibrium free energy F is taken to generalize the usual Ginzburg-Landau action to N order parameters that crucially include subdominant orders not stabilized in equilibrium, and formally reads Here, the bar denotes complex conjugation, and summation over repeated indices is implied.Coefficients A ij , B ijmn , and D µν ij are tensorial generalizations of the usual Ginzburg-Landau coefficients for quadratic, quartic, and gradient contributions.We use subscript Latin indices (i, j, m, n = 1, . . ., N ) to index different order parameters, and superscript Greek indices (µ, ν) to index spatial directions.The theory accurately represents a multi-dimensional free energy landscape in the vicinity of the critical temperature T c for the equilib-rium even-parity instability.Finally, we couple the superconductor to light by introducing minimal coupling to a gauge field in velocity gauge, ∇ µ −→ ∇ µ + i 2e ℏ A µ (t), where −e < 0 is the electron charge.We discard subsequent gradient terms by considering only spatially homogeneous irradiation and order parameters (∇∆ i = 0).From this, one can derive Euler-Lagrange equations of motion ∂L ∂∆i − ∂ t ∂L ∂(∂t∆i) = 0 that describe light-induced dynamics of competing orders.
The structure of the TDGL action is dictated solely by parity and angular momentum conservation, and importantly permits a coupling between even-and odd-parity order parameters already to linear order in A µ (t).If the lattice has C n rotational symmetry with n ≥ 2, order parameters for s-, p-, d-wave instabilities can be enumerated by their angular momentum eigenvalues e i2πl/n , with l = 0, l = ±1, l = ±2 (modulo n) respectively.An appealing selection rule permits light-induced coupling C µ ij between superconducting orders i, j with ∆l = ±1 at linear order in the field.Notably, this dictates that the Lagrangian couples an equilibrium s-wave order parameter to odd-parity p-wave order already at linear order in A µ (t), in agreement with the quasiparticle dynamics of Fig. 2(a).Conversely, second order in A µ (t) contributions couple same-parity order parameters with ∆l = 0 or ∆l = ±2, capturing the excitation of the conventional amplitude mode.This observation has intriguing consequences for Higgs mode experiments (see Summary and outlook).
To address the existence of a metastable triplet superconductor as a local minimum of the free energy landscape, we explicitly compute the multi-component Ginzburg-Landau coefficients for microscopic Hamiltonians, thus connecting the top-down phenomenological approach to the bottom-up microscopic approach of the previous section.Starting from a generic multiband Hamiltonian with effective attractive pairing interactions, the generalized Ginzburg-Landau coefficients can be computed diagrammatically as with correlation functions depicted in Fig. 3, and Matsubara Green's functions for particles and holes given by In the above equations, Tr ≡ 1 βL d k,ωn tr, indicating a sum over momenta, orbital indices, and Matsubara frequencies.The gradient terms C µ ij and D µν ij are calculated from expanding Π (2) ij (q) in powers of q.We first demonstrate key qualitative features by example of the conceptually simplest single-band model on a rectangular lattice [Fig.4a] that captures the essential physics for switching to metastable odd-parity superconductor.We choose a FIG. 3. Effective Action via Diagrams.Diagrammatic representation of Eqs. ( 7) and ( 8), which summarize how to calculate the generalized Ginzburg-Landau coefficients as an effective field theory starting from a microscopic Hamiltonian.
hopping anisotropy δt = 0.2, chemical potential µ = −0.1, and interaction parameters v s = −1.75,v px = v py = −2.5, in units of hopping.As spin-orbit interactions in centrosymmetric crystals require at least a two-orbital unit cell, we break SU (2) via a small Zeeman magnetic field ∆ Zeeman = 0.1 [see Methods].The free energy parameters are computed from the microscopic model, and the phenomenological inertial/damping coefficients are set to γ = 1.0, η = 0.1.The temperature T = 0.1 lies below the critical temperature for the s-wave and p x -wave channel, with s-wave pairing stabilized in equilibrium.Crucially, the coupling between these channels is first-order in A(t).For simplicity, we assume that pairing interactions in the d-wave and extended s-wave channels vanish.Furthermore, an x-polarized pulse couples solely to p x order, yielding an appealingly simple minimal action where we abbreviate A ii ≡ a i , B iiii ≡ b i , and B iijj ≡ b ij , and assume equal inertial (γ) and damping (η) coefficients without loss of generality.The first two lines describe decoupled TDGL actions for s and p ≡ p x order parameters.The third line describes their quartic coupling that dictates the emergence of a metastable minimum, and the last line couples s and p orders to linear order in the light field A µ .In the rectangular model, this coupling keeps the s-and p x -wave order parameters completely real, allowing the dynamics to be completely captured by a visualizable two-dimensional free energy contour plot (see Figure 4b).
Intriguingly, as a function of pulse width and fluence, one finds a contiguous parametric family of Gaussian pulses that result in the final state settling into the pure ∆ px stationary point of the free energy, as in Figure 4e.The presence of a clear threshold for the width of the Gaussian pulse is evidence that dynamical inversion symmetry breaking, which occurs more strongly with a tighter pulse width, is crucial for efficient coupling between opposite-parity orders.We find switching to be most reliable for strong driving near the Bardasis-Schrieffer mode frequency [Fig.2d], which reads and is determined by deviations from the global free energy minimum in the direction of the target instability [see Methods].Light polarized in the x-direction drives the order parameter strongly in the direction of a closely competing p x minimum in its free energy landscape, and for Gaussian pulses centered at this frequency with widths on the order of a few resonant periods (see Figure 4e), the system is efficiently switched to the target metastable state.Figures 4d and 4e together reveal a trade-off whereby the pulse needs to be short enough to allow for strong dynamical inversion symmetry breaking while being wide enough (i.e.spectrally sufficiently narrow) to avoid excitation of quasiparticles across the gap (excitations that are not captured by TDGL).The latter constraint leads us to detune the central frequency of the Gaussian pulse to be slightly below Ω BS , decreasing the spectral over-lap with the quasiparticle continuum while still maintaining a large overlap with the BS mode.[See the Supplementary Material for further details.] We now turn to a realistic model of a conventional centrosymmetric superconductor with strong spin-orbit coupling on the honeycomb lattice, with a putative closely-competing triplet pairing instability [Fig.5a].This models the lowenergy physics near the Dirac points of several systems of recent interest, including kagome superconductors [56][57][58], honeycomb materials with large spin-orbit coupling and recently-reported superconductivity [59][60][61], and many moiré heterostructures of transition metal dichalcogenides [62][63][64][65][66]. Rotation symmetry and energetics dictate that the chiral triplet states ∆ p±ip ≡ 1 √ 2 ∆ px ± i∆ py are stable local free energy minima, while the nodal ∆ px and ∆ py states individually are not.Figure 5c depicts the free energy landscape with an equilibrium s-wave minimum and proximal p ± ip local minima.This is therefore a promising example of a truly metastable non-equilibrium chiral superconducting state in a system with conventional s-wave order in equilibrium.
The key parameter allowing singlet-triplet switching at first order in the gauge field A(t) is the coefficient C µ ij in Eq. ( 5).For the honeycomb model, up to first order in the spin-orbit coupling strength λ, this matrix element becomes a function of the (matrix-valued) non-Abelian Berry connection A µ kmn ≡ i u mk |∂ kµ u nk .Since the energy eigenvalues of the Hamiltonian are at order λ 2 and higher, the dominant contribution to C µ for small λ is therefore purely geometric in nature, controlled by the Berry connection.Notably, this effect is distinct from quantum-geometric corrections to the superfluid weight [67,68], and can be equivalently expressed via the Fermi surface Berry connection at low temperatures (see Methods).
To approximate the energetics in conventional but strongly spin-orbit coupled superconducting systems, [69,70] we choose parameters so that the nearest neighbor hopping amplitude is 1 eV, the spin-orbit strength is 0.1 eV, the system is hole-doped with chemical potential µ = −0.4eV and the critical temperatures of the s and p ± ip states are ∼ 10 K for effective renormalized interactions v s = −1.25 eV, v p = −2.3eV.In this scenario, the superconducting gaps are between 1 − 10 meV.The inertial coefficients Γ (2) ij in Eq. ( 4) are set so that the frequency of small amplitude oscillations of ∆ i matches the expected Higgs mode frequency, (eq) i |/ℏ.This sets the time scale of interest to be on the order of picoseconds, with relevant frequencies between 1 − 10 THz.The time scale for relaxation (parametrized by Γ (1) ij in Eq. ( 4)) is longer than the Higgs mode period by the dimensionless factor β|∆ (eq) i | [31], which in our case is about 10. Lastly, assuming a lattice constant on the order of Angstroms, our results call for peak field strengths on the order of ∼ 100 kV/cm for ultrafast pulses.
Figure 5e reveals a broad region of pulse widths (∼ 0.1 ps) and fluences (∼ 1.0 mJ/cm 2 ) that lead to successful switching from s to a chiral triplet p + ip state via an ultrafast circularly polarized pulse.The resulting order parameter trajectory in a subspace of the free energy landscape is shown in Fig. 5c and illustrates key features of the switching protocol, whereby the dynamical inversion symmetry breaking kicks the order parameter in the vicinity of a chiral instability.To allow the system to settle into one of the chiral p ± ip states within TDGL theory, we set a slight difference in their relaxation rates (η ± = (1 ± δ) η with δ = 0.5).Though introduced as a phenomenological parameter, one can expect this to arise microscopically from pump-induced time-reversal symmetry breaking of the environmental degrees of freedom that provide dissipation, an effect which is not captured by TDGL theory.The chiral states are stable against fluctuations in the amplitudes and relative phases of all competing order parameters that are supported by nearest-neighbor pairing interactions, providing a promising proof-of-concept for engineering a metastable topological superconducting state.

Summary and outlook
We present light-induced dynamical inversion symmetry breaking as a generic route to explore competing triplet pairing instabilities in conventional superconductors with strong spin-orbit coupling.Identifying a sub-gap odd-parity BS mode that couples linearly to light, we find that driving this mode with a tailored two-color or ultrafast light pulse can switch the system to a metastable unconventional superconducting state.We illustrate this mechanism for minimal models of superconductors with closely competing pairing insta-bilities.
Immediate future directions include searching for experimental signatures of a sub-gap opposite-parity Bardasis-Schrieffer mode, in addition to applying this methodology to realistic tight-binding models of candidate materials, including 1T ′ WTe 2 and moire heterostructures of transition metal dichalcogenides.Following experiments that use pump-probe microscopy to detect the conventional Higgs mode via a linear scaling between pump intensity and the amplitude of resulting gap fluctuations [16][17][18], a simple optical protocol may be to search for gap fluctuations scaling with the square root of the pump intensity, indicating an amplitude mode coupling linearly to light.Additionally, an intriguing open problem is under what general conditions does a subleading stationary point in the free energy become a local minimum.We have found that the only examples of metastable states are Kramers partners with chiral pairing, i.e. p x ± ip y order parameters in lattices where the p x and p y orders are degenerate.It would be fruitful to further understand the necessary conditions for a system to support a metastable superconducting instability.Another interesting direction would be to model heating effects from irradiation with an ultrafast pulse, accounting for microscopic mechanisms for heat dissipation not included in our methods.Lastly, throughout this work we use effective renormalized interaction strengths v i as phenomenological parameters for emergent pairing at low energy scales.Detailed modeling of these parameters in candidate materials, via renormalization group treatments and incorporating RPA effects [39], is an important direction for future work.

Gap function equations of motion via quasiparticle dynamics
Quasiparticle equations of motion for a superconductor coupled to light are derived for a generic BCS-like Hamiltonian where ĉ † kασ creates an electron in the state with momentum k, orbital index α ∈ {1, . . ., N }, and spin σ ∈ {↑, ↓}.The dependence of the one-body Hamiltonian h σ αβ (k) on spin encodes spin-orbit coupling that reduces the SU (2) spin rotation symmetry to spin-z conservation.L denotes the linear system size.We decompose the four-fermion interaction vertex V αβα ′ β ′ (k, k ′ ) in terms of a set of orthonormal basis functions A standard mean-field decoupling of the interaction in the Cooper channel introduces the superconducting gap function where the components ∆ i are classified in terms of the irreducible representations of the crystal point group.The meanfield decomposition of Eq. ( 13) can be written in Bogoliubovde Gennes (BdG) form as where The equal-time Nambu Green's function obeys equations of motion (setting ℏ = 1) where n F-D (X) ≡ (1 + exp(βX)) −1 is the Fermi-Dirac distribution function.Combined with the instantaneous selfconsistency equation for order parameters, this yields equations of motion for the gap function in terms of the Heisenberg dynamics of the Bogoliubov quasiparticles.

Symmetry analysis of the Lagrangian
Suppose we demand that the Lagrangian be invariant under some symmetry transformation on the order parameters ∆ i and the gauge field A µ (t), for tensors D and R. Let R be orthogonal ((R ⊤ ) µν R νρ = δ µρ ) and choose a basis for is a product of two (not necessarily distinct) eigenvalues of R. For example, if one considers C n rotational symmetry on a lattice with n ≥ 2, the constraints given by Eq. ( 22) yield the l selection rules discussed in Results.

Effective action from the path integral
We start with the partition function written in terms of an imaginary-time path integral over the Grassmann fields ψ, ψ, with a free-electron action given by, and an interacting action given by, As before, α, β denote orbital indices, σ denotes a spin index, and v i denotes the interaction decomposed into pairing channels.A Hubbard-Stratonovich transformation [71] decouples S int in the Cooper channel in terms of auxiliary fields where the Nambu spinor Ψ and Gor'kov Green's operator G −1 are defined in terms of their momentum components and Matsubara frequency dependence as follows, Note that here, k can be interpreted as the internal momentum of a Cooper pair, while q corresponds to the external momentum.One can then perform the path integral over the Grassmann fields Ψ, Ψ, which gives, (29) Using ln det M = tr ln M, this results in an effective action Expanding this to fourth order in ∆ yields, which, comparing to Eq. ( 5), allows one to compute the generalized Ginzburg-Landau coefficients for completing orders as described in Results.

Rectangular lattice model
Consider a rectangular lattice with hopping amplitudes t x > t y > 0, on-site interaction U , and nearest neighbor interaction U ′ .As discussed in the main text, breaking SU (2) symmetry in a single-band model requires including a small the Zeeman splitting due to a z-aligned magnetic field.The single-particle dispersion reads When the Zeeman splitting energy ∆ Zeeman is nonzero, there are separate spin-up and spin-down Fermi surfaces, making it difficult to define singlet/even order parameters and triplet/odd order parameters in this model.For this reason, we will consider the applied magnetic field to be zero initially, only to be adiabatically turned on before the optical pulse and adiabatically turned off after the optical pulse.We note that a small ∆ Zeeman has a negligible effect on the Ginzburg-Landau free energy, but must necessarily be present to allow coupling between singlet and triplet orders.Conversely, multiband models discussed below allow for the inclusion of spin-orbit coupling, obviating a magnetic field.
The form factors f i (k) for this model up to nearestneighbor interactions are tabulated in Table I.There are three order parameters in the A irrep (s-wave) and two in the B irrep (p-wave).In equilibrium, the system settles into a combination of the three s-wave orders dictated by the relative strengths of the on-site and nearest neighbor interactions.For simplicity, we re-express these order parameters in the eigenbasis of Π (2) ij (0) [Eq.( 7)], defining ∆ s in the main text to be the order parameter with the lowest eigenvalue.We then define |∆ stotal | in Figure 4c to be ∆ 2 slocal + ∆ 2 sext,x + ∆ 2 sext,y .

Honeycomb lattice model with spin-orbit coupling
We now consider a honeycomb lattice with Kane-Mele spin-orbit coupling [55] as well as effective on-site and nearest-neighbor attractive interactions, with nearest and next-nearest-neighbor hopping Here, λ parameterizes the strength of spin-orbit coupling, µ denotes the chemical potential, t denotes the hopping parameter, {d i } denotes the three nearest-neighbor lattice vectors, and {a i } denotes the three next-nearest-neighbor lattice vectors (a . For the interaction term, we consider an on-site attraction U and nearest-neighbor attraction U ′ , decomposing V αβα ′ β ′ (k, k ′ ) as in Eq. ( 14) using the basis functions tabulated in Table II.Setting U = U ′ , the ∆ s state is found to have the highest critical temperature, followed by two degenerate p x and p y states.These have a d-wave k-space structure, but are odd under sublattice exchange and have angular momentum l = ±1.Chiral superpositions are abbreviated as 1 √ 2 ∆ px ± i∆ py ≡ ∆ p±ip .

BS mode frequency and free energy barrier between two order parameters
The Lagrangian [Eq.( 4)] results in Euler-Lagrange equations of motion where N µ ij ≡ i 2e ℏ C µ ij and M µν ij ≡ (i 2e ℏ ) 2 D µν ij are defined to absorb factors of i 2e ℏ from minimally substituting ∇ µ −→ ∇ µ + i 2e ℏ A µ (t) into Eq.( 5) and assuming spatiallyhomogeneous order parameters and irradiation.
Focusing on the simplest case of two competing order parameters that are only coupled at linear order in A µ (t) (e.g.s-wave and p x -wave), the free energy for two order parameters can be written concisely as, where a i ≡ A ii , b i ≡ B iiii , and b 12 equals any fourthorder coefficient coupling two ∆ 1 's and two ∆ 2 's (these are all equivalent in this case).This is in agreement with Ref. [72].Assuming the system is initially in the equilibrium state, = 0 (assume overall phase equals zero without loss of generality), the equation of motion for ∆ 2 up to linear order in ∆ 2 and A µ (t) is given by, where 22 (we assume the inertial and damping coefficients are diagonal), and M 2 is given by, One can formally integrate this equation of motion as, which yields a resonant response when ω 2 ≡ Ω 2 = M 2 /γ 2 .We identify this as the BS mode frequency Ω BS .One can also derive from Eq. (36) the free energy barrier between the two states by finding the saddle point surrounded by the ∆ 1 minimum, the ∆ 2 minimum, and the ∆ 1,2 = 0 maximum.The free energy difference between this saddle point and the global minimum is given by, where ϕ ≡ arg (∆ 2 /∆ 1 ) is the relative phase between the two order parameters.For successful switching, the time-integrated pulse must supply enough kinetic energy ∼ 1 2 γ 2 |∂ t ∆ 2 | 2 to overcome this free energy barrier.Since the matrix element N µ ij ∝ C µ ij depends linearly on the spin-orbit coupling strength λ (for small λ), we expect the requisite fluence for switching to scale as ∼ (t/λ) 2 (where t is the nearestneighbor hopping amplitude).

C µ ij in terms of the Fermi surface Berry connection
Assuming the chemical potential exists in a band with Bloch functions |u k ⟩, and dispersion ξ k (measured with respect to the chemical potential) with a large gap ∆ gap to all other bands that we formally take to infinity, In the β → ∞ limit, the integrand diverges at the Fermi surface, meaning the dominant contribution to C µ ij at low temper-atures and low spin-orbit strengths comes from the geometry of the Bloch states at the Fermi surface.

Data availability
Data sets are available from the corresponding author on request.
Irrep.Pairing In Figures 4 and 5 of the main text, we report parameter regions as a function of inverse pulse width and fluence that lead to successful switching in our time-dependent Ginzburg-Landau (TDGL) simulations.In Figure 5e, we comment that the window of pulse widths was chosen such that fraction of fluence with frequencies in the quasiparticle continuum (defined as ω > 2∆ (eq) , where 2∆ (eq) is the equilibrium superconducting gap) never exceeds 10%.Here we provide plots of this fluence fraction as a function of the pulse parameters.We compute this via the ratio The results are plotted in Figure 1.Here we use a path integral to derive the Ginzburg-Landau effective action presented in the main text from an arbitrary multiband Bloch Hamiltonian.We start with the partition function for the microscopic model written in terms of a path integral over the electronic Grassmann fields ψ, ψ with a free-electron action given by, and an interaction given by, We write the momentum indices in the "center of mass" frame of the Cooper pair and decompose the interaction function where i indexes some subset of the irreducible representations (irreps) of the crystal symmetry group, and f i αβ (k) is the momentum-space form factor associated to irrep i.We now perform a Hubbard-Stratonovich transformation to bring this action to a form that is quadratic in ψ, ψ.We decouple in the Cooper channel in terms of order parameters ∆i,q , ∆ i,−q , e −S int [ ψ,ψ] Here, q corresponds physically to the external momentum of a Cooper pair, in contrast to k, which corresponds to the internal momentum.The action can be compactly written using Nambu fields in which case the partition function takes the form, with the inverse Gor'kov Green's function Ĝ−1 .In Matsubara frequency representation, Ψ(τ ) = 1 √ β ωn Ψ ωn e −iωnτ and therefore ∂ τ → −iω n , β 0 dτ → β ωn .(We assume that ∆ has no imaginary time dependence.)This allows writing, where the matrix elements Ĝ−1 ωn k,k ′ of the inverse Green's function operator are given by Ĝ−1 ωn Here, we use a condensed notation in which ĥkσ is the matrix with elements h σ αβ (k), and fk,i is the matrix with elements f i αβ (k).Evaluating the path integral over the Grassmann fields Ψ, Ψ, one obtains with an effective action given by S eff [ ∆, ∆] = −tr ln Ĝ−1 − i q ∆i,q We expand S eff as a functional Taylor series around the saddle point ∆ = 0, S eff [ ∆, ∆] = S eff | ∆=0 + i i ′ q q ′ ∆i,q δ 2 S eff δ ∆i,q δ∆ i ′ ,q ′ ∆=0 ∆ i ′ ,q ′ + 1 4 ijmn ∆i,0 ∆m,0 δ 4 S eff δ ∆i,0 δ∆ j,0 δ ∆m,0 δ∆ n,0 ∆=0 ∆ j,0 ∆ n,0 The first-order term vanishes by virtue of ∆ = 0 being a saddle point.The zeroth-order term is simply the free-electron action, which we can neglect since it just contributes an overall constant.The second-order and fourth-order terms read ijmn (0) δ q0 ∆m,0 ∆ n,0 ∆ j,−q , (13) in terms of tensors Π (2) ij (q) and Π ijmn (0), described in the main text.

Second-order terms
The second-order in ∆ terms (including the gradient terms), are computed in terms of Π (2) ii ′ (q).Comparing Eqs. ( 12) and ( 13), we find, Π ii ′ (q) = − 1 βL d δ 2 δ ∆i,q δ∆ i ′ ,q ′ tr ln Ĝ−1 where we use the shorthand Ĝ0 ≡ Ĝ ∆=0 . For the honeycomb lattice considered in the main text, the matrix elements of the Green's function Ĝ0 are given by Ĝ0 ωn where we define Ĝ↑ k and Ĝ↑ k (non-script G's) for later convenience (and we will sometimes suppress the iω n dependence to declutter notation).With these definitions in order, one can carefully compute Π (2) ii ′ (q) making use of the matrix elements of these operators, Here, q ′ must equal −q, as imposed by the Kronecker deltas.

FIG. 2 .
FIG.2.Dynamical inversion symmetry breaking and odd-parity Bardasis-Schrieffer (BS) modes.a-b Demonstration of a linear-in-|A| and linear-in-λ (spin-orbit parameter) coupling between an equilibrium order parameter ∆s and a px-wave order parameter ∆p x .A two-color x-polarized pulse that breaks dynamical inversion symmetry is applied with various small magnitudes |A|, and example time-plots of the resulting fluctuations in |∆p x (t)| are plotted in b.The maxima of these fluctuations are plotted versus |A| in a, for five different magnitudes of the spin-orbit coupling strength λ in units of hopping.c Time dynamics resulting from a two-color pulse that breaks dynamical inversion symmetry (with a 2 : 1 frequency ratio).d Time dynamics resulting from a two-color pulse that preserves dynamical inversion symmetry (with a 3 : 1 frequency ratio).Both c and d start from purely even superconducting order (indicated by blue curves), and the onset of odd orders (indicated by red and pink curves) is suppressed when dynamical inversion symmetry is unbroken (see legend in bottom right).e Short-time deviations in the s-wave and px-wave order parameters for weak driving with x-polarized light as a function of frequency Ω, revealing a BS mode for the px-wave order parameter that is sub-gap (indicated by the unshaded region).

FIG. 4 .
FIG. 4. Engineering metastable triplet superconductors with light.a Schematic of a minimal rectangular tight-binding model of competing singlet and triplet pairing instabilities.b Free energy landscape plot of Re(∆p x ) vs Re(∆s), with the time-dependent trajectory in c plotted in green.c Time plots of Ax(t) and |∆i(t)| for stotal [Methods], px, and py order parameters subjected to an example Gaussian pulse.d Schematic plot showing the spectral profile of the pulse in c in relation to the Bardasis-Schrieffer mode frequency ΩBS and the edge of the quasiparticle continuum (QPC).e A parametric plot of fluence versus inverse pulse width showing the final state induced by a family of Gaussian pulses, with frequency Ω = 0.75 ΩBS centered slightly below the BS mode frequency.A star indicates the pulse parameters used in b-d.

FIG. 5 .
FIG. 5. Switching between trivial s and topological triplet p ± ip pairing in honeycomb superconductors.a Schematic of a minimal honeycomb-lattice conventional s-wave superconductor with a competing chiral triplet pairing instability.b Order parameter trajectory in the free energy landscape as a function of the magnitudes of the s and p+ip order parameters, demonstrating the stable switching of the equilibrium order parameter to a metastable chiral state.A small difference in the relaxation rates for the two chiral states is assumed.c Dynamics of A x (t), A y (t), and |∆i(t)| for ∆s and ∆p±ip.d Schematic plot of the spectral profile of the pulse in c, showing strong overlap with the BS mode frequency ΩBS and minimal overlap with the quasiparticle continuum (QPC).e Parametric plot in terms of inverse pulse width (in THz) and fluence (in mJ/cm 2 ) showing which Gaussian pulses centered at Ω = 0.75 ΩBS lead to successful switching to the metastable p + ip state.The window is chosen so that no more than 10% of the fluence overlaps with the QPC.A star indicates the pulse parameters used in b-d.

Figure 1 :
Figure 1: The amount of pulse fluence contributing to quasiparticle excitations as a function of the pulse parameters used in Figure 4 (Left, Rectangular model) and Figure 5 (Right, Honeycomb model).The upper plots calculate the fluence.The lower plots calculate the log of the ratio between the quasiparticle continuum fluence and the full fluence of the pulse.