Collective modes in pumped unconventional superconductors with competing ground states

TL;DR

This paper explores how ultrafast pump-probe experiments can reveal collective excitations like Bardasis-Schrieffer and Higgs modes in unconventional superconductors with competing ground states, highlighting control via pulse polarization.

Contribution

It demonstrates the presence and controllability of collective modes in pumped superconductors with multiple pairing channels using a theoretical model and spectral analysis.

Findings

Identification of Bardasis-Schrieffer modes in spectral gap dynamics

Control of collective mode signatures through pulse polarization

Strong signatures of collective modes in quasiparticle distribution dynamics

Abstract

Motivated by the recent development of terahertz pump-probe experiments, we investigate the short-time dynamics in superconductors with multiple attractive pairing channels. Studying a single-band square lattice model with spin-spin interaction as an example, we find the signatures of collective excitations of the pairing symmetries (known as Bardasis-Schrieffer modes) as well as the order parameter amplitude (Higgs mode) in the short-time dynamics of the spectral gap and quasiparticle distribution after an excitation by a pump pulse. We show that the polarization and intensity of the pulse can be used to control the symmetry of the non-equilibrium state as well as frequencies and relative intensities of the contributions of different collective modes. We find particularly strong signatures of the Bardasis-Schrieffer mode in the dynamics of the quasiparticle distribution function. Our work shows the potential of modern ultrafast experiments to address the collective excitations in unconventional superconductors and highlights the importance of sub-dominant interactions for the non-equilibrium dynamics in these systems.

Figures

Equations

$$H=-t\sum _{⟨ij⟩,\sigma }{c}_{i\sigma }^{†}{c}_{j\sigma }+J\sum _{⟨ij⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$J_{\mathbf{k},{\mathbf{k}}^{\prime }}\equiv$$

$$+V_{p_x}{\gamma }_{\mathbf{k},p_x}{\gamma }_{{\mathbf{k}}^{\prime },p_x}+V_{p_y}{\gamma }_{\mathbf{k},p_y}{\gamma }_{{\mathbf{k}}^{\prime },p_y},$$

$$H_{\text{MF}}\simeq \sum _{\mathbf{k}\sigma }{\xi }_{\mathbf{k}}{c}_{\mathbf{k}\sigma }^{†}{c}_{\mathbf{k}\sigma }+\sum _{\mathbf{k},l}\left[{\Delta }_l{\gamma }_{\mathbf{k},l}{c}_{\mathbf{k}\uparrow }^{†}{c}_{-\mathbf{k}\downarrow }^{†}+\text{h.c.}\right]$$

$${\mathbf{s}}_{\mathbf{k}}=\frac{1}{2}\left(c_{\mathbf{k}\uparrow }^{†},c_{-\mathbf{k}\downarrow }\right)\mathbf{\sigma }\left(\begin{array}{c}{c}_{\mathbf{k}\uparrow }\\ c_{-\mathbf{k}\downarrow }^{†}\end{array}\right),$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k}}⟩={\mathbf{B}}_{\mathbf{k}}\times ⟨{\mathbf{s}}_{\mathbf{k}}⟩.$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k},s}⟩$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k},d_{x^2-y^2}}⟩$$

$$H=\sum _{\mathbf{k}}{\mathbf{B}}_{\mathbf{k}}\cdot {\mathbf{s}}_{\mathbf{k}}+({\xi }_{\mathbf{k}+\frac{e}{c}\mathbf{A}}-{\xi }_{\mathbf{k}-\frac{e}{c}\mathbf{A}})(c_{\mathbf{k}\uparrow }^{†}{c}_{\mathbf{k}\uparrow }+c_{-\mathbf{k}\downarrow }{c}_{-\mathbf{k}\downarrow }^{†}),$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k}}⟩=\frac{i}{\hslash }⟨[H,{\mathbf{s}}_{\mathbf{k}}]⟩={\mathbf{B}}_{\mathbf{k}}\times ⟨{\mathbf{s}}_{\mathbf{k}}⟩.$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k},s}⟩$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k},d_{x^2-y^2}}⟩$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k},d_{xy}}⟩$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k},g_{xy(x^2-y^2)}}⟩$$

$${\xi }_{\mathbf{k},\mathbf{A},s}$$

$${\xi }_{\mathbf{k},\mathbf{A},d_{x^2-y^2}}$$

$${\xi }_{\mathbf{k},\mathbf{A},d_{xy}}$$

$${\xi }_{\mathbf{k},\mathbf{A},g_{xy(x^2-y^2)}}$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k},s}⟩$$

$$\frac{d}{dt}⟨{\mathbf{s}}_{\mathbf{k},d_{x^2-y^2}}⟩$$

$$-i\omega \delta {\mathbf{s}}_{\mathbf{k}}={\mathbf{B}}_{\mathbf{k}}^{\text{eq}}\times \delta {\mathbf{s}}_{\mathbf{k}}+\delta {\mathbf{B}}_{\mathbf{k}}\times {\mathbf{s}}_{\mathbf{k}}^{\text{eq}},$$

$$\delta {\mathbf{s}}_{\mathbf{k}}=$$

$$\left(\begin{array}{cc}{\chi }_{xx}& {\chi }_{xy}\\ {\chi }_{yx}& {\chi }_{yy}\end{array}\right)\left(\begin{array}{c}\delta {\Delta }_l^{\prime }\\ \delta {\Delta }_l^{\prime \prime }\end{array}\right)$$

$${\chi }_{xx}$$

$${\chi }_{xy}$$

$${\chi }_{yy}$$

Full text

Ultrafast pump-probe techniques became recently a powerful tool to probe the temporal evolution of symmetry broken states and relaxation in conventional and unconventional superconductors.Pashkin et al. (2010); Beck et al. (2011); Fausti et al. (2011); Matsunaga and Shimano (2012); Dal Conte et al. (2012); Beck et al. (2013); Mansart et al. (2013); Matsunaga et al. (2014a); Hu et al. (2014); Casandruc et al. (2015); Matsunaga et al. (2017); Katsumi et al. (2018) An intense pulse couples non-linearly to the Cooper pairs of the superconductor and, as was argued theoretically, should lead to a coherent excitation of the Higgs amplitude mode, i.e. $|\Delta(t)|$ performs a damped oscillation with frequency $\omega_{H}=2|\Delta(\infty)|$.Volkov and Kogan (1974); Amin et al. (2004); Barankov et al. (2004); Yuzbashyan et al. (2005, 2006); Barankov and Levitov (2006); Papenkort et al. (2007); Krull et al. (2014); Yuzbashyan et al. (2015); Tsuji and Aoki (2015); Murotani et al. (2017); Chou et al. (2017) Nonlinear terahertz spectroscopy has enabled the observation of the Higgs mode in conventional superconductors in the form of a free or forced oscillation and the resulting third-harmonic generationMatsunaga and Shimano (2012); Matsunaga et al. (2014b, 2017). Interestingly, this technique has been also recently applied to the unconventional superconductors such as high-Tc cuprates with the $d$-wave symmetry of the superconducting gap Katsumi et al. (2018); Chu et al. (2019) where some additional oscillations have been reported.Chu et al. (2019)

In contrast to conventional superconductors, where the pairing is driven by the attractive electron-phonon interaction, the pairing interaction in unconventional superconductors is most likely of repulsive nature. To overcome the net repulsion among the quasiparticles, the superconducting gap has to change its sign across different parts of the Fermi surface, which typically yields the superconducting gap of a lower symmetry than an isotropic $s$-wave. For example, it is generally known that the antiferromagnetic spin fluctuations peaked near wave vector ${\bf Q}_{AF}=(\pi,\pi)$ within a single-band model on a square lattice give rise to a $d_{x^{2}-y^{2}}$-wave symmetry of the superconducting gap, yet states having other symmetries, such as strongly anisotropic sign-changing (extended) $s$-wave symmetry and the $d_{xy}$-symmetry, are closely competing. As a result, the temporal dynamics of single-band unconventional superconductors might be significantly richer than that of the conventional ones.Barankov et al. (2004); Papenkort et al. (2007); Yuzbashyan et al. (2005); Tsuji and Aoki (2015)

In this manuscript we analyze the short-time dynamics in a single-band unconventional superconductor with multiple competing pairing symmetries. In particular, we consider a single-band model of fermions on a square lattice interacting via spin-spin interaction. The interaction can be decoupled into various pairing channels with different symmetry. Varying the band filling we find two competing even-parity superconducting states forming a typical phase diagram of an unconventional superconductor where different ground states can be accessed by doping. Studying the system driven out of equilibrium by a laser pulse, we show how the collective signatures of symmetries different from a given ground state symmetry, known as Bardasis-Schrieffer modesBardasis and Schrieffer (1961); Maiti and Hirschfeld (2015); Maiti et al. (2016); Böhm et al. (2018) in the context of $s$-wave ground state in the equilibrium, evolve as a function of doping and the polarization direction. Depending on the polarization direction of the incoming light the tetragonal symmetry is broken, which necessarily leads to a mixing of $s$- and $d$-wave symmetries in non-equilibrium. Furthermore, we show that the particle distribution acquires an additional $d_{x^{2}-y^{2}}$-character and oscillates dominantly with the Bardasis-Schrieffer mode frequency, which may be observed in time-resolved ARPES experiments. Our study highlights the important role of sub-dominant pairing states in the short-time dynamics of unconventional superconductors and we identify the signatures of the resulting collective modes that can be observed in future pump-probe experiments.

Our starting point is a model of fermions interacting via a spin-spin interaction

[TABLE]

where $c_{i,\sigma}^{(\dagger)}$ are fermionic annihilation (creation) operators on site $i$ and spin $\sigma$, $t$ is the hopping integral between nearest neighbors and $S_{i}^{\alpha}=\frac{1}{2}\sum_{s,s^{\prime}}c_{is}^{\dagger}\sigma_{ss^{\prime}}^{\alpha}c_{is^{\prime}}$ are the spin-1/2 operators. This model was considered before in the context of cuprates, pnictides and heavy fermion systems.Davis and Lee (2013) In the momentum space the tight-binding energy dispersion is given by ${\xi_{\mathbf{k}}=-2t(\cos(k_{x})+\cos(k_{y}))-\mu}$. The spin-spin interaction can also be transformed into momentum space and decoupled into a number of superconducting channels

[TABLE]

where ${V_{s}=V_{d_{x^{2}-y^{2}}}=-3J/2}$ is the even-parity spin-singlet interaction with ${\gamma_{\mathbf{k},s}=(\cos(k_{x})+\cos(k_{y}))/2}$ and ${\gamma_{\mathbf{k},d_{x^{2}-y^{2}}}=(\cos(k_{x})-\cos(k_{y}))/2}$ form factors, respectively. The remnant components are repulsive odd-parity spin-triplet coupling constants ${V_{p_{x}}=V_{p_{y}}=J/2}$ with the $p_{x}$- and $p_{y}$-wave form factors ${\gamma_{\mathbf{k},p_{x(y)}}=\sin(k_{x(y)})}$. As we show below, the light pulse we consider does not couple the odd- and even- parity channels. As we focus on the regime, where the ground state has even parity, the odd and even parity solutions cannot mix both in- and out of equilibrium. In the mean-field approximation the Hamiltonian reduces to

[TABLE]

where $l=s,d_{x^{2}-y^{2}}$ and the ${\Delta_{l}=-V_{l}\sum_{\mathbf{k}}\gamma_{\mathbf{k},l}\langle c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow}\rangle}$ are the $s$- and $d_{x^{2}-y^{2}}$-wave component of the total superconducting order parameter ${\Delta_{\mathbf{k}}=\Delta_{s}\gamma_{\mathbf{k},s}+\Delta_{d_{x^{2}-y^{2}}}\gamma_{\mathbf{k},d_{x^{2}-y^{2}}}}$. In what follows, we work at $T=0$. Minimizing the energy, we obtain the equilibrium values of the superconducting order parameters $\Delta_{s},\;\Delta_{d_{x^{2}-y^{2}}}$ as a function of $n$, shown in Fig. 1. We concentrate on the region of the phase diagram where $s$-wave and $d$-wave symmetries are neighbors, but avoid the regime of coexistence of both order parameters, i.e. a possible $s+id$ state, as this state (and the similar case of $d+id$ state) and its non-equilibrium dynamics were discussed previously, Maiti and Hirschfeld (2015); Müller et al. (2018); Kirmani and Dzero (2019). Instead, we focus on pairing fluctuations of the sub-dominant symmetry close to the transition points.

To simplify further calculations, we introduce the Anderson pseudospin notationAnderson (1958)

[TABLE]

where $\bm{\sigma}=\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)^{T}$ are the Pauli-matrices. Using this vector, one can recast the Hamiltonian $H=\sum_{\mathbf{k}}\mathbf{B}_{\mathbf{k}}\cdot\mathbf{s}_{\mathbf{k}}$, which has the form of a set of (pseudo-)spins $\mathbf{s}_{\mathbf{k}}$ coupled to a (pseudo-)magnetic field ${\mathbf{B}_{\mathbf{k}}=\left(2\Delta^{\prime}_{\mathbf{k}},2\Delta^{\prime\prime}_{\mathbf{k}},2\xi_{\mathbf{k}}\right)^{T}}$, with the notation $\Delta_{\mathbf{k}}=\Delta_{\mathbf{k}}^{\prime}-i\Delta_{\mathbf{k}}^{\prime\prime}$. At zero temperature the thermal expectation values of the pseudospin components are given by ${\left\langle s_{\mathbf{k}}^{x}\right\rangle=-\frac{\Delta^{\prime}_{\mathbf{k}}}{2E_{\mathbf{k}}}}$, $\left\langle s_{\mathbf{k}}^{y}\right\rangle=-\frac{\Delta^{\prime\prime}_{\mathbf{k}}}{2E_{\mathbf{k}}}$, and ${\left\langle s_{\mathbf{k}}^{z}\right\rangle=-\frac{\xi_{\mathbf{k}}}{2E_{\mathbf{k}}}}$ with the quasiparticle energy dispersion $E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^{2}+|\Delta_{\mathbf{k}}|^{2}}$.

To investigate the collective modes in our model, we study the equations of motion for the pseudospin expectation values ${\mathbf{s}_{\mathbf{k}}}$ that have the form of Bloch equations:

[TABLE]

This equation together with the self-consistency equation for the superconducting gaps $\Delta_{s,d}(t)=-\sum_{\mathbf{k}}V_{s,d}\gamma_{\mathbf{k},s,d}(\left\langle s_{\mathbf{k}}^{x}\right\rangle(t)-i\left\langle s_{\mathbf{k}}^{y}\right\rangle(t))$ yields a closed set of coupled differential equations, which defines the temporal evolution of all relevant quantities. To drive the system out of equilibrium we model the electric field of a laser pulse by including a time-dependent vector potential $\mathbf{A}(t)$ via the Peierls substitution. This results in ${\mathbf{B}_{\mathbf{k}}=\left(2\Delta^{\prime}_{\mathbf{k}},2\Delta^{\prime\prime}_{\mathbf{k}},\xi_{\mathbf{k}+\frac{e}{c}\mathbf{A}}+\xi_{\mathbf{k}-\frac{e}{c}\mathbf{A}}\right)^{T}}$ (see Supplementary Material for details).

To consider the possibility of exciting order parameter symmetries, different from the ground state one it is instructive to split the equations of motion in Eq. (5) into different symmetry channels. As the pulse temporarily breaks $C_{4}$-rotational symmetry, different symmetry representations may mix. In particular, we can decompose the pseudomagnetic field into all even parity irreducible representation for the tetragonal $D_{4h}$ symmetry as follows: ${\mathbf{B}_{\mathbf{k}}=\mathbf{B}_{\mathbf{k},s}+\mathbf{B}_{\mathbf{k},d_{x^{2}-y^{2}}}+\mathbf{B}_{\mathbf{k},d_{xy}}+\mathbf{B}_{\mathbf{k},g_{xy(x^{2}-y^{2})}}}$, where we define $\mathbf{B}_{\mathbf{k},l}=\left(2\Delta_{l}^{\prime}\gamma_{\mathbf{k},l},2\Delta_{l}^{\prime\prime}\gamma_{\mathbf{k},l},\xi_{\mathbf{k},\mathbf{A},l}\right)^{T}$, where we have introduced the symmetrized notations for $B_{\mathbf{k}}^{z}=\xi_{\mathbf{k}+\frac{e}{c}\mathbf{A}}+\xi_{\mathbf{k}-\frac{e}{c}\mathbf{A}}\equiv\xi_{\mathbf{k},\mathbf{A},s}+\xi_{\mathbf{k},\mathbf{A},d_{x^{2}-y^{2}}}+\xi_{\mathbf{k},\mathbf{A},d_{xy}}+\xi_{\mathbf{k},\mathbf{A},g_{xy(x^{2}-y^{2})}}$. There appears no odd-parity component of the pseudomagnetic field $\mathbf{B_{\mathbf{k}}}$, as can be seen immediately from the definition above.

Let us now discuss the symmetry mixing in non-equilibrium. Denoting the ground state symmetry as $l_{0}$, one finds that $\frac{d}{dt}\langle s_{\mathbf{k},l^{\prime}}^{x}\rangle=0$ and $\frac{d}{dt}\langle\mathbf{s}_{\mathbf{k},l^{\prime}}^{y}\rangle=\mathbf{B}_{\mathbf{k},l}\times\langle\mathbf{s}_{\mathbf{k},l_{0}}\rangle$, where $l^{\prime}=l\otimes l_{0}$ is the symmetry of a product of functions with symmetries $l$ and $l_{0}$. Therefore, a finite $l$-symmetry component of $\mathbf{B}_{\mathbf{k}}$ field is needed to induce pairing correlations of symmetry $l^{\prime}$. Moreover, the self-consistency equations imply that inducing an out-of equilibrium order parameter $\Delta_{l^{\prime}}$ additionally requires a non-zero $V_{l^{\prime}}$, making the consideration of the sub-dominant pairing interaction crucial.

In our case, the pseudospin expectations values $\langle\mathbf{s}_{\mathbf{k}}\rangle$ at $t=0$ are the equilibrium values, with $s$- or $d_{x^{2}-y^{2}}$-wave symmetry. Consequently, one requires $B_{\mathbf{k},d_{x^{2}-y^{2}}}^{z}=\xi_{\mathbf{k},\mathbf{A},d_{x^{2}-y^{2}}}$ to be finite at non-zero t to induce a finite $d_{x^{2}-y^{2}}$-wave order parameter in the $s$-wave ground state and vice versa. In particular, a finite $\xi_{\mathbf{k},d_{x^{2}-y^{2}}}\sim\cos(A_{x})-\cos(A_{y})$ (see Supplementary Material for details) can be induced by a vector potential $\mathbf{A}$ with polarization at $\phi\approx 0$, where $\phi$ is the polar angle in momentum space. Importantly, the effect manifestly depends on the polarization of $\mathbf{A}$, suggesting a possibility to control the induced order parameter symmetry in non-equilibrium. As $\xi_{\mathbf{k},\mathbf{A},d_{xy}},\xi_{\mathbf{k},\mathbf{A},g_{xy(x^{2}-y^{2})}}$ and the odd-parity components are identically zero, only $s$- and $d_{x^{2}-y^{2}}$-components remain in Eq. (5) leading to

[TABLE]

We note that this result does depend on the absence of spatial variations of ${\bf A}$. However, for THz light used in the experiments,Matsunaga et al. (2013) the wavelength is of the order $\sim$100 $\mu$m which is much larger than, e.g., superconducting coherence length that rarely exceeds $\sim$100 nm (as estimated by the upper critical field $H_{c2}$Yin et al. (2009); Shan et al. (2011); Grissonnanche et al. (2014)), justifying the assumption.

We discuss now the short-time dynamics described by the equations (Collective modes in pumped unconventional superconductors with competing ground states). We integrate the equations numerically using Runge-Kutta method with a momentum grid of $513\times 513$ points for fixed $V_{s}=V_{d}=-0.4t$ varying $n$ and the polarization of ${\bf A}(t)$. Let us consider first the situation of an extended $s$-wave ground state symmetry for $\Delta_{\mathbf{k}}$, i.e. $\Delta_{s}\neq 0$ and $\Delta_{d_{x^{2}-y^{2}}}=0$ at $n=0.411$. The transition into a $s+id_{x^{2}-y^{2}}$-state occurs at the filling $n\approx 0.415$. We start the analysis by driving the model out of equilibrium with a vector potential in $k_{x}$-direction. The vector potential $\mathbf{A}$ is for simplicity chosen as $\mathbf{A}=\mathbf{A}_{0}\theta(t-\tau)\theta(2\tau-t)$ simulating two pulses at $t=\tau$ and $t=2\tau$. Here we choose $\frac{e}{c}|\mathbf{A}_{0}|=0.01$ and $\tau\Delta=0.05$. Typically a Gaussian pulse shape is used in similar approachesPapenkort et al. (2007), with a half-duration $t_{1/2}\ll\Delta^{-1}$ to work in the non-adiabatic regime. Our choice of the pulse is equivalent to this approach in the non-adiabatic limit, but we found it to be advantageous for computations. The results of the numerical calculations are shown in Fig. 2. As one can see, the vector potential (pump pulse) induces a non-zero $d_{x^{2}-y^{2}}$-wave order parameter $\Delta_{d_{x^{2}-y^{2}}}$ in the extended $s$-wave ground state, which shows undamped oscillations at a frequency $\omega_{\text{BS}}$. Note that $|\Delta_{d_{x^{2}-y^{2}}}|$ shown in Fig. 2 oscillates at twice the frequency of $\Delta_{d_{x^{2}-y^{2}}}$ due to the absence of sign changes in the former. According to above considerations a finite $\Delta_{d_{x^{2}-y^{2}}}$ can only be induced if $\xi_{\mathbf{k},\mathbf{A},d_{x^{2}-y^{2}}}$ is finite for some period of time. Thus by applying the pulse along the diagonal $k_{x}=k_{y}$, for which $\xi_{\mathbf{k},d_{x^{2}-y^{2}}}$ remains zero, there is no induced oscillations of the $d_{x^{2}-y^{2}}$-wave order parameter. In both cases the $s$-wave order parameter oscillates with its Higgs mode $\omega_{H}=2\Delta_{s,\text{max}}$, where $\Delta_{s,\text{max}}=\max\{\Delta_{s}\gamma_{\mathbf{k}}|\mathbf{k}\in FS\}$ is the maximum gap size on the Fermi surface. The frequency $\omega_{\text{BS}}$ is smaller than $\omega_{H}$ due to residual attractive interaction and therefore undamped, while $\omega_{H}$ shows the weak dampingBarankov et al. (2004) observed for isotropic s-wave states. We attribute this to the gap being nearly isotropic along the Fermi surface at the doping we consider. In Fig. 3 we show the fluence dependence of $\omega_{BS}$ and $\omega_{H}$. It can be seen that in the case where $\omega_{BS}$ and $\omega_{H}$ are well separated (Fig. 3 b) ), $\omega_{BS}$ is not directly affected by the light pulse, unlike $\omega_{H}$. However, once the Higgs mode frequency becomes close to $\omega_{BS}$, both start being suppressed by the increasing fluence and eventually become indistinguishable. Note that instead of using a pulse one could also by hand set a finite $\Delta_{d_{x^{2}-y^{2}}}\neq 0$ as initial value, as it was shown above that the pulse also generates a finite $d_{x^{2}-y^{2}}$-wave component. We find that the quench scenario yields results equivalent to Fig. 2

To verify whether the oscillations at $\omega_{BS}$ are due to the existence of the Bardasis-Schrieffer mode in the sub-dominant pairing channel we analyze Eq. (5) in the linear regime (i.e. weak driving) and compare the frequencies of the resulting modes to the ones observed in the non-equilibrium (see Fig. 2). In particular, we linearize the Eq. (5) around the equilibrium state at $T=0$ and perform a subsequent Fourier transform ${-i\omega\delta\mathbf{s}_{\mathbf{k}}=\mathbf{B}^{\text{eq}}_{\mathbf{k}}\times\delta\mathbf{s}_{\mathbf{k}}+\delta\mathbf{B}_{\mathbf{k}}\times\mathbf{s}^{\text{eq}}_{\mathbf{k}}}$. Numerical solution of the homogeneous part of the resulting equations (see Supplementary Material for details) gives three types of solutions: (i) $\omega=0$ mode for the ground state pairing symmetry independent of doping; (ii) $\omega=2\Delta_{l,\text{max}}$ modes for the both types of the ground states $l=s$ and $l=d_{x^{2}-y^{2}}$ (Higgs mode), and (iii) Bardasis-Schrieffer mode, which frequency depends on the proximity to the phase boundary between extended $s$-wave and $d_{x^{2}-y^{2}}$-wave symmetries in the equilibrium. We show the frequency of the Bardasis-Schrieffer mode in the $s$-wave ground state in Fig. 4 together with the results of the numerical solution of the full non-linear equation, Eq. (5) as a function of doping. For values of $n$ far away from the transition this mode merges into the Higgs mode, while close to the phase transition it softens as is expected, signalling the transition to the new ground state. The results of linearized calculation of $\omega_{\text{BS}}$ are in good agreement with the results for the full non-linear dynamics.

We also considered the excitation of the Bardasis-Schrieffer mode in the $d_{x^{2}-y^{2}}$-wave ground state, which should also exist once the $s$-wave ground state is close enough. As the $d_{x^{2}-y^{2}}$-wave gap is nodal, all excited modes are strongly damped. However, we still see a clear damped mode at $\omega<2\Delta_{d,max}=\omega_{H}$. Due to the finite size effects and the limited accuracy of the numerical integration, the error bar is quite broad. Therefore, we only show the frequency of this mode close to the $s+id$ transition in the inset of Fig. (4). To obtain these frequencies we can again equivalently use both, a pulse or a quench, for example, if one sets $\Delta_{s}=0.1\Delta_{d_{x^{2}-y^{2}}}$. Note, for a better numerical accuracy all frequencies in the inset of Fig. (4) are obtained via quench.

The presence of the signatures of the Bardasis-Schrieffer mode in the order parameter dynamics raises the possibility of an enhanced third-harmonic generation when the pump frequency roughly matches $\omega_{\text{BS}}$. Unlike the Higgs mode, where the issues of quasiparticle contribution is still under debate,Cea et al. (2016, 2018); Tsuji and Aoki (2015); Matsunaga et al. (2017); Murotani and Shimano (2019); Kumar and Kemper (2019) the signature of a Bardasis-Schrieffer mode would be definitive as $\omega_{\text{BS}}<2\Delta_{0}$, below the edge of the quasiparticle continuum, at least for fully gapped $s$-wave superconductors. At the same time, the amplitude oscillation of the order parameter (Higgs oscillation) was also predicted to show up in the tr-ARPES experiments.Nosarzewski et al. (2017) The quasiparticle distribution function that can be measured in these experiments can be addressed in our model. In particular, the quasiparticle distribution function is determined by the $s_{\mathbf{k}}^{z}$ component of the pseudospin $s_{\mathbf{k}}^{z}=\frac{1}{2}\left(\langle c_{\mathbf{k}\uparrow}^{\dagger}c_{\mathbf{k}\uparrow}\rangle+\langle c_{\mathbf{k}\downarrow}^{\dagger}c_{\mathbf{k}\downarrow}\rangle-1\right)$. In equilibrium this quantity is equal to $-\xi_{\mathbf{k}}/(2E_{\mathbf{k}})$ and thus is $C_{4}$-symmetric. Due to the perturbation via an electric field, the tetragonal symmetry is temporarily broken down to $C_{2}$-symmetry and $s_{\mathbf{k}}^{z}$ develops a finite $d_{x^{2}-y^{2}}$-wave component $s_{\mathbf{k}}^{z,d}$, i.e. the quasiparticle distribution along $x$- and $y$-axis becomes asymmetric. Therefore, one expects that the non-equilibrium particle distribution acquires the information on $\omega_{\text{BS}}$ and $\omega_{H}$. To investigate this in detail we define the quantity $I(t)=\int_{-\pi}^{\pi}dk_{x}(2s_{\mathbf{k}=(k_{x},0)}^{z,d}(t))$, which describes the integrated $d_{x^{2}-y^{2}}$-wave component of the quasiparticle distribution along the $k_{y}=0$ cut and is equivalent to the number of particles with $k_{y}=0$. In Fig. 5 we show $I(t)$ for the same parameters as in Fig. 2 for $\mathbf{A}$ oriented along $x$-axis. One can readily see that it oscillates mostly with $\omega=\omega_{\text{BS}}$, where we find the amplitude of the oscillations increases with an increasing pump strength. The integration along the $k_{x}=0$ cut leads to similar results but shifted by a phase $\varphi=\pi$.

To conclude we analyze the short-time dynamics in single band unconventional superconductor with multiple competing pairing symmetries. Driving the system out of equilibrium with a light pulse (modeled as a time-dependent vector potential), we show how the collective signatures of symmetries different from a given ground state symmetry, known as Bardasis-Schrieffer modesBardasis and Schrieffer (1961); Maiti and Hirschfeld (2015); Maiti et al. (2016); Böhm et al. (2018) in the context of $s$-wave ground state in the equilibrium, evolve as a function of doping and the polarization direction. Depending on the polarization direction of the incoming light the tetragonal symmetry is broken, which necessarily induces a coupling of $s$- and $d$-wave symmetries in a non-equilibrium. Furthermore, we show that the particle distribution acquires an additional $d_{x^{2}-y^{2}}$-character, due to coupling to the vector potential and that this quantity shows a dominant signature of the Bardasis-Schrieffer mode frequency, which might be observed in time-resolved ARPES experiments. Therefore we conclude that taking the sub-dominant pairing channels and the corresponding interactions into account is important, when discussing polarization dependent excitation of unconventional superconductors.

Acknowledgements: M.A. M. and I.M. E. were supported by the joint DFG-ANR Project (ER 463/8-1). P.A.V. acknowledges the support by the Rutgers University Center for Materials Theory Postdoctoral fellowship. I. P. is supported by the ANR grant (ANR-15-CE30-0025).

I Symmetrized Equations of Motion

Here we give detailed information about the derivation of the equations of motion and its symmetrization. After performing a minimal substitution to include the vector potential the Hamiltonian has the following form

[TABLE]

where ${\mathbf{B}_{\mathbf{k}}=\left(2\Delta^{\prime}_{\mathbf{k}},2\Delta^{\prime\prime}_{\mathbf{k}},\xi_{\mathbf{k}+\frac{e}{c}\mathbf{A}}+\xi_{\mathbf{k}-\frac{e}{c}\mathbf{A}}\right)^{T}}$. The last term transforms as $\frac{1}{2}\left(c_{\mathbf{k}\uparrow}^{\dagger},c_{-\mathbf{k}\downarrow}\right)\sigma_{0}\begin{pmatrix}c_{\mathbf{k}\uparrow}\\ c_{-\mathbf{k}\downarrow}^{\dagger}\end{pmatrix}$, where $\sigma_{0}$ is the unit matrix and thus commutes with all other pseudospins. Therefore the structure of the equations of motions remains unaffected

[TABLE]

Now one can rewrite $\mathbf{s}_{\mathbf{k}}$ and $\mathbf{B}_{\mathbf{k}}$ in Eq. (S2) as a sum of all even parity irreducible representations of the tetragonal $D_{4h}$ symmetry group as discussed in the main text. Using the product rules for the irreducible representations of the $D_{4h}$ symmetry group one can split this equations into four orthogonal parts

[TABLE]

For the $z$-component of the pseudomagnetic field $B_{\mathbf{k}}^{z}=\xi_{\mathbf{k}+\mathbf{A}}+\xi_{\mathbf{k}-\mathbf{A}}\equiv\xi_{\mathbf{k},\mathbf{A}}^{s}+\xi_{\mathbf{k},\mathbf{A}}^{d_{x^{2}-y^{2}}}+\xi_{\mathbf{k},\mathbf{A}}^{d_{xy}}+\xi_{\mathbf{k},\mathbf{A}}^{g_{xy(x^{2}-y^{2})}}$ one can integrate out each symmetry. For the chosen dispersion they have the form

[TABLE]

In our model the occurring order parameters only have either $s$- or $d_{x^{2}-y^{2}}$-wave symmetry, which simplifies equations (I) into

[TABLE]

as those two representations form a subgroup.

II Short-time dynamics in the d-wave ground state

To verify that the existence of a sub-dominant $s$-wave channel in a $d_{x^{2}-y^{2}}$-wave superconductor also gives rise to an additional collective mode in the short-time dynamics we solve Eq. (6) numerically for the $d_{x^{2}-y^{2}}$-wave ground state. Due to the nodal gap structure more quasiparticles are excited, which leads stronger damping. However, this also means that the total error bars in the main text are due to both finite accuracy of the Fourier transform and numerical errors due to finite grid size effects.

We consider a $d_{x^{2}-y^{2}}$-wave ground state at $n\approx 0.437$ and drive the system out of equilibrium by quenching the $s$-wave order parameter into $\Delta_{s}=0.1\Delta_{d}$. The result in Fig. S1 shows an induced Higgs mode at $\omega_{H}=2|\Delta_{d,\text{max}}(\infty)|$ and also a mode $\omega_{\text{BS}}<\omega_{H}$.

The error bars in the main text are given by the frequency, which is closest to the half width of the peak. Due to the computational error and the stronger damping, the peak of the Fourier transform is broader than in Fig. 2. Therefore in addition to the discrete resolution of the Fourier transform the error bars are increased due to computational accuracy.

III Linearization of Equations of Motion

Here, we perform the linearization of the equations of motion around the equilibrium state at $T=0$. We start with linearizing Eq. (5) and performing a subsequent Fourier transform to obtain

[TABLE]

with the equilibrium pseudomagnetic field $\mathbf{B}_{\mathbf{k}}^{\text{eq}}=(2\Delta_{l_{0}}^{\prime}\gamma_{\mathbf{k},l_{0}},2\Delta_{l_{0}}^{\prime\prime}\gamma_{\mathbf{k},l_{0}},2\xi_{\mathbf{k}})^{T}$ and the equilibrium pseudospin expectation value $\mathbf{s}^{\text{eq}}_{\mathbf{k}}=-\mathbf{B}_{\mathbf{k}}^{\text{eq}}/(2E_{\mathbf{k},l_{0}})$. We use $l_{0}$ to denote the groundstate symmetry. After some rearrangement of this expression one obtains a matrix equation

[TABLE]

where $\Delta_{\mathbf{k}}^{\text{0}}=\Delta_{l_{0}}\gamma_{\mathbf{k},{l_{0}}}$, which corresponds to the equilibrium order parameter in either $s$-wave or $d_{x^{2}-y^{2}}$-wave ground state depending on band filling. Eq. (8) can be then plugged into the definition of the superconducting gap amplitude $\Delta_{l}=-\sum_{\mathbf{k}}\gamma_{\mathbf{k},l}V_{l}(\delta s_{\mathbf{k}}^{x}-i\delta s_{\mathbf{k}}^{y})$. Consequently, one obtains a linear equation for the real and imaginary part of the gap $\Delta_{l}$, where the vector $\bm{\lambda}_{l}(\omega)$ includes necessary information of the perturbating vector potential in $\delta B^{z}_{\mathbf{k}}(\omega)$.

[TABLE]

The entries of the response function $\hat{\chi}(\omega)$ are given by

[TABLE]

for each $l$ respectively. The entries of $\bm{\lambda}_{l}(\omega)$ are given by ${\lambda_{l}^{\prime}=V_{l}\sum_{\mathbf{k}}\frac{2\xi_{\mathbf{k}}\gamma_{\mathbf{k},l}\Delta_{\mathbf{k}}^{0}\delta B_{\mathbf{k}}^{z}}{2E_{\mathbf{k},l_{0}}(\omega^{2}-4E^{2}_{\mathbf{k},l_{0}})}}$ and ${\lambda_{l}^{\prime\prime}=V_{l}\sum_{\mathbf{k}}\frac{2i\omega\gamma_{\mathbf{k},l}\Delta_{\mathbf{k}}^{0}\delta B_{\mathbf{k}}^{z}}{2E_{\mathbf{k},l_{0}}(\omega^{2}-4E^{2}_{\mathbf{k},l_{0}})}}$. Collective modes of the system are excited for values of $\omega$ for which the determinant of $\hat{\chi}$ vanishes, i.e., $\det(\hat{\chi})=0$ independent of the strength of perturbation.