Engineering Higgs dynamics by spectral singularities

We generalize the dynamical phase diagram of a Bardeen-Cooper-Schrieffer condensate, considering attractive to repulsive, i.e., critical quenches (CQ) and a non-constant density of states (DOS). We show that different synchronized Higgs dynamical phases can be stabilized, associated with singularities in the density of states (DOS) and different quench protocols. In particular, the CQ can stabilize an overlooked high-frequency Higgs dynamical phase related to the upper edge of the fermionic band. For a compensated Dirac system we find a Dirac-Higgs mode associated with the cusp singularity at the Fermi level, and we show that synchronized phases become more pervasive across the phase diagram. The relevance of these remarkable phenomena and their realization in ensembles of fermionic cold atoms confined in optical lattices is also discussed.

We generalize the dynamical phase diagram of a Bardeen-Cooper-Schrieffer condensate, considering attractive to repulsive, i.e., critical quenches (CQ) and a non-constant density of states (DOS). We show that different synchronized Higgs dynamical phases can be stabilized, associated with singularities in the density of states (DOS) and different quench protocols. In particular, the CQ can stabilize an overlooked high-frequency Higgs dynamical phase related to the upper edge of the fermionic band. For a compensated Dirac system we find a Dirac-Higgs mode associated with the cusp singularity at the Fermi level, and we show that synchronized phases become more pervasive across the phase diagram. The relevance of these remarkable phenomena and their realization in ensembles of fermionic cold atoms confined in optical lattices is also discussed.
Introduction. Many-body systems are characterized by the occurrence of correlated phenomena, which have no counterpart in the few-body realm. In this perspective, the spontaneous symmetry breaking (SSB) of any Hamiltonian symmetry by the establishment of a finite order parameter represents one of the fundamental examples [1][2][3][4]. Consequences of SSB include the appearance of superfluid and superconducting phases in condensed matter systems, as well as the occurrence of a finite mass of the intermediate vector bosons, the carrier particles of the weak force in the Standard Model [5][6][7][8].
The excitations on top of the SSB ground state, in systems with continuous (gauge) symmetries, consists of Nambu-Goldstone (phase) modes and massive Higgs (amplitude) modes. In the Standard Model, the latter manifest themselves as the Higgs boson, whose experimental observation led to the 2013 Nobel Prize in physics [9,10].
While the Lorentz invariance stabilizes the Higgs mode in the high-energy context, its signatures at low-energies are less sharp and crucially depend on the strength of the interactions [26,[35][36][37][38][39]. A comprehensive picture of the Higgs mode features across the BEC-BCS crossover has recently been obtained in a fermionic cold atom ensemble [40], reinvigorating the interest of the cold atom's community in the signatures of SSB and Higgs mechanism also in the few-body limit [41,42].
Upon temporal variation of a control parameter, many-body dynamics may display several peculiar features, which are reminiscent of the behavior of thermodynamic functions at transition points [43]. In particular, theoretical investigations uncovered the appearance of dynamical phase transitions following interaction quenches in strongly correlated systems [44][45][46][47][48]. These dynamical transitions occur both as order parameter modulations and as singularities in the Löschmidt echo dynamics [49,50], which have been observed in quantum optics experiments [51,52]. These two phenomena have been shown to be deeply intertwined both between each other and with the existing equilibrium transitions [53,54], except for few examples [55][56][57].
Following the current perspective, we consider a critical quench (CQ) of the BCS model, i.e., an attractive to repulsive interaction quench where the sign of the coupling constant is flipped. Then, the system, which is prepared in the superconducting equilibrium state for attractive interaction, evolves according to a repulsive Hamiltonian, whose equilibrium ground-state would be a normal (non-superconducting) gas. Thanks to the flexibility of our approach, we can target a wide range of different models parametrized by diverse density of states (DOS).
The CQ protocol displays the distinctive features of the Higgs mode dynamics found in pioneering works [58,59], including synchronized oscillations of the order parameter. However, contrary to the common belief, there is no a unique synchronized Higgs phase (SHP) but different SHPs can be stabilized depending on the model and the protocol. The oscillation frequency can be determined by any spectral singularity. This includes those singularities not connected with the SSB but present in the bare DOS. Thus, a combination of protocol and the optical lattice in a cold atom system allows engineering generalized synchronized Higgs phases on-demand.
Model. In order to prove our picture, we consider a weak-coupling fermionic condensate with s-wave pairing described by the BCS model and subject to a sudden quench arXiv:2205.06826v1 [cond-mat.quant-gas] 13 May 2022 of the pairing interaction. The Hamiltonian can be written as where ξ k = ε k − µ measures the energy from the Fermi level µ and the pairing interaction λ(t) = λ i θ(−t)+λ f θ(t) with θ the Heaviside step function. Hereĉ † kσ (ĉ kσ ) is the usual creation (annihilation) operator for fermions with momentum k and spin σ.
Due to the all-to-all interaction assumed in the last term of Eq. (1), a mean-field approach becomes exact in the thermodynamic limit. Thus, we shall consider the BCS mean-field Hamiltonian which can be written, using the Anderson pseudospin formulation [60] aŝ Here, b k (t) = (2∆ (t) , 0, 2ξ k ) represents an effective magnetic field vector for the 1 2 Without loss of generality, we have assumed that the equilibrium BCS order parameter ∆ is real, which remain true over time due to electron-hole symmetry. The instantaneous BCS order parameter is given by where symbols without hat denote the expectation value of operators in the time-dependent BCS state. At equilibrium, the 1 2 -pseudospins align in the direction of their local fields b i k = (2∆ i , 0, 2ξ k ) in order to minimize the system's energy.
The system is prepared with an initial interaction λ i and a gap parameter satisfying the equilibrium gap equation, where ρ(ξ) is the DOS and the bandwidth satisfies, W ≫ ∆ i ensuring the system is in the weak-coupling regime. The interaction is then suddenly changed to the final value, λ f and the gap parameter is studied as a function of time. The pseudospins evolve obeying the Bloch-like equation of motion interacting with all the others pseudospins via the selfconsistent order parameter ∆(t). We take N = 10 4 , 10 5 pseudospins within an energy range of W = 40∆ i , 60∆ i and 80∆ i around µ. We consider both a constant DOS and the DOS of a Dirac semimetal like graphene (see white background insets of Fig. 1). It is convenient to parameterize the quench by the variable δ ≡ W λ f − W λi , whose value is closely related to the one already used to parametrize non-critical quenches [59,61].
The out-of-equilibrium dynamics shows collective effects and single-particle (pseudospin) excitations. Similar to the case of periodic driving [62][63][64][65], we find that the latter are k +∆ 2 is the quasiparticle dispersion defined in terms of an average order parameter∆ computed on a large time window after the quench. The resulting DOS has edge singularities at energies ±∆ and ±W 2 for quasiparticles (2∆ and W for Ω L fluctuations, with ∆ ≪ W ).
Dynamical Phase diagram. In Fig. 1 we present the dynamical phase diagram for three different bandwidths and a constant DOS (left) and for the graphene-like DOS (right). Panels (a) and (b) show key values of the order parameter characterizing the dynamics, while panels (c) and (d) show the generalized SHP oscillation frequencies. Full lines were obtained exploiting the integrability of the model through a Lax roots analysis (see Supplementary Information, SI for details) and were checked by numerical simulations (circles). The frequencies of Higgs modes are labeled according to a singularity of the quasiparticle DOS to which we show below  to be associated, namely, lower edge (l), upper edge (u), Van Hove (V) and Dirac point (D). In Fig. 1 the curves to the right of the vertical dashed line in panel (a) and (c) reproduce the results of Ref. [59] where three dynamical phases were found for non-critical quenches λ f λ i > 0 (see SI Fig. S1 for details of the dynamics at colored dots): For small increase and decrease of the attraction in the quench (small δ ), the superconducting order parameter shows damped oscillations with a Higgs frequency associated with the lower edge of the quasiparticle DOS, ω l = 2∆ ∞ and saturates to a constant value ∆ ∞ at long times which therefore coincides with∆ (yellow shading). For large quenches there are two possible outcomes. Decreasing the pairing constant beyond a critical point (δ > δ c+ = π 2), the system goes to a gapless regime (∆ ∞ = 0) with an overdamped dynamics (red region). Increasing the pairing above a critical point (δ < δ c− = −π 2), the system synchronizes and the order parameter oscillates between the values ∆ + and ∆ − with a fundamental Higgs frequency equal to twice the average order parameter [66], ω l = 2∆ (magenta area). We will refer to this well known phenomena [58,59,61] as the "lower-edge SHP".
Using δ as the control parameter has the advantage that for large enough bandwidth, the results for λ f λ i > 0 become independent of the bandwidth, so the curves for different W fall almost on top of each other. Thus, the dynamical phases obtained after non-critical quenches display a certain degree of universality, which is correctly captured in terms of the variable δ. The upper scale (and the position of the vertical dashed line) corresponds to the particular case, W ∆ i = 40. The region close to the vertical dashed line has a final interaction in the strong coupling regime which is out of our scope, so data is omitted.
The region to the left of the vertical dashed line shows the result of the CQ. A different synchronized regime is found where the order parameter oscillates with symmetric amplitudes ∆ + and ∆ − around zero instead of having a finite average∆ [see SI Fig. S1(c) for the detailed dynamics at the colored dots]. This zero order-parameter average (ZOPA) behavior is reminiscent of the time-crystal phases found with periodic driving [63]. In contrast to the purely attractive interaction quench, the amplitudes and the Higgs frequency are strongly dependent on the bandwidth. In this case, the SHP frequency converges to the upper edge of the fluctuation DOS (ω u → W ), when the final interaction is repulsive and small (large negative δ). For δ < −8 (large repulsive λ f ) the amplitude becomes vanishing small and the dynamical phase can not be distinguished from a gapless state. Therefore, for large negative and positive δ (λ f → ±0) the system converges to the same gapless phase.
The occurrence of SHP associated with upper and lower edges of the fluctuation/quasiparticle DOS suggests that singularities act as nucleation centers in frequency space to stabilize the synchronized phases during the non-linear dynamics. We expect this mechanism to be ubiquitous, thus leading to the appearance of Higgs mode signatures in any conventional superconducting system with singular DOS. In order to confirm this expectation, we have studied the dynamical phase diagram of a Dirac system with a graphene-like DOS (right column) where two additional singularities are present already in the bare DOS: one at the Fermi level (Dirac point) and the Van Hove singularities at ξ V k ≈ ±6.66∆ i (see DOS in the white background inset of Fig. 1).
Interestingly, the phase diagram of the graphene-like model turns out to be quite different from the flat-DOS case. The damped regime (yellow region) is completely absent and synchronization occurs even for an arbitrary small quench [see zoomed region inset in Fig. 1(b) and the detailed dynamics at the light blue point in Fig. 2(a)]. Also, differently from the constant DOS model, the synchronization phenomenon takes place both for an increase and a decrease of the pairing interaction. In the latter case, decreasing enough λ f , the pseudospins effectively decouple from each other and the gapless regime is recovered (red region). A representative dynamics of this phase is shown in Fig. 2(a) in red using the parameters indicated with red circle in Fig. 1 (b).
Between the δ = 0 and the gapless phase, a quite rich  Fig. S1(b)]. The synchronized Dirac-Higgs phase, in contrast, appears much more harmonic.
Panel (c) and (d) of Fig. 2 exemplify the dynamics for the CQ corresponding to the matching full dots in Fig. 1 (b), (d), where the upper-edge SHP is excited. The overall appearance is very similar to the case of a flat DOS [SI Fig. S1 (b), (d)] showing ZOPA behavior. However, an extra weak modulation appears, which is revealed as a new frequency ω V in the Fourier transform (d). This frequency matches twice the Van Hove singularity in the DOS 2ξ V k , thus as expected, a synchronized Van Hove-Higgs mode can be excited. Its amplitude decreases in time, indicating that the mode is damped, although with quite a long coherence time, as witnessed by the narrow peak in the FT.
To fully characterizes the dynamical phases and the emergence of synchronization in the system, it is instructive to analyze the ξ k -resolved FT of the pseudospin dynamics. As we are interested in the pairing dynamics, we show the FT of the x-component of pseudospins for the graphene-like model in Fig. 3 (SI Fig. S2 shows the same information for the flat DOS model).
Single-particle (pseudospin) excitations appear as dispersive features, while synchronized collective modes appear as vertical lines. For the non-ZOPA synchronized phase (light blue dots in Fig. 1) single-particle excitations appear in the dynamics as Larmor precessions with frequency Ω L (k) = 2 ξ 2 k +∆ 2 as shown in panel Fig. 3(a) with the large black dots. The same panel shows that the lower edge Higgs mode is not determined only by the quasiparticles participating in the edge singularity of Ω L (k) at the frequency 2∆. Indeed, the vertical feature emerging from 2∆ witnesses that all quasiparticles in this ξ k window are synchronized and participate in the collective mode. Thus, the edge singularity of the dispersion can be seen as a nucleation center in frequency space given a "rhythm" which is followed by the rest of the quasiparticles due to the interactions.
In addition, of the main Ω L (k) dispersion, Floquet side bands appear analogous to the bands observed under periodic drive [63]. Here, of course, an external periodic drive is not present and the bands are self-generated by the action of the lower-edge Higgs mode with frequency ω l yielding replicas at Ω L (k) + nω l with n = −1, 0, 1, 2, ... and weaker features at −Ω L (k) + nω l with n = 0, 1.
Panel (b) shows that the upper edge singularity of the DOS is enough to trigger the SHP (appearing as a vertical feature) with frequency ω u ≃ W provided that the appropriate protocol, i.e. the CQ, is used. Thus, the 1 √ ω − 2∆ divergence present at the lower edge of the DOS is not a prerequisite to stabilize Higgs modes. No feature associated with the Van Hove-Higgs mode ω V can be seen here, since the numerical analysis has been performed on a time window larger than its coherence time in order to have high-frequency resolution.
In panel (b) the ZOPA manifests as a gapless linear dispersion of single particle excitations, Ω L (k) = 2ξ k (black dots). The same linear dispersion appears for the other ZOPA modes: the Dirac-Higgs shown in panel (c) and the gapless mode shown in (d). For the former, the collective nature of the synchronization is also evident from the vertical feature at ω D .
Conclusions. We have shown that different synchronized Higgs phases can be excited in a BCS system by choosing an appropriate quench protocol. For a given system, the frequency of the mode is determined by singularities in the DOS with small corrections due to quasiparticle interactions. The previously known lower-edge Higgs mode appears at the same frequency of a singularity in the equilibrium particle-particle response. The upper edge SHP is reminiscent of antibound states appearing in the equilibrium pairing response of repulsive systems. However, at equilibrium the mode is not present in particle-hole symmetric situations [67] while here the mode is stabilized in an out-of-equilibrium setting. Thus, generalized Higgs modes do not appear to have always an equilibrium counterpart.
Our findings provide an innovative interpretation to the Higgs mode dynamics, which appear as synchronized quasiparticles oscillations nucleated by DOS singularities. This, in turns, implies that any spectral singularity can give rise to Higgs-mode like oscillations given a suitable quench.
The observation of the previously reported lower-edge SHP in real superconductors is hindered by its decay in other excitations and its weak coupling to light [26,35,36,61]. Furthermore, signatures of Higgs dynamics in pump and probe experiments [20,21,68] cannot be clearly distinguished from other Raman active modes [69] with similar frequencies, but different underlying mechanisms. The experimental detection of other Higgs modes presented here will probably encounter similar difficulties in solid state systems, as our picture is based on the BCS model, whose integrable nature does not account for thermalization.
A proper understanding of the generalized Higgs dynamics and its decay at strong coupling, as well as the possible relation with the equilibrium characteristic of the SSB phase, may be obtained by direct observation in Fermi superfluids of cold atoms. Indeed, recent improvements in the control and observation of ultracold atoms in optical lattices [70,71] allowed the study of both equilibrium and transport properties of Fermi systems on the lattice [72][73][74][75], paving the way to the realization of the generalized Higgs dynamics described here. In particular, an artificial graphene-like lattice with tunable interactions has been realized [76]. Another route is to use cold atoms in an optical cavity, which has recently been proposed as a BCS simulator [77].
Interestingly, the search for more than one Higgs boson is a subject of intense search also in high-energy scattering experiments [78]. This kind of probes, however, are more akin to equilibrium responses in condensed matter. Instead, the strongly out-of-equilibrium physics investigated here may find parallels in the electroweak transition of relevance for early universe cosmology and baryogenesis [79].

LAX ROOTS ANALYSIS
Because of the integrability of the BCS model, the different dynamical phases described in the main text can be obtained by analyzing the integrals of motion. For this purpose, it is useful to define the so-called Lax vector [1,2] defined as a function of an auxiliary complex parameter y, where z is a unit vector along the z direction, and S k is the pseudospin texture before the quench. The square of the Lax vector is a conserved quantity under time evolution with the BCS Hamiltonian. Therefore, the complex roots of such vector (in the following Lax roots) are also conserved.
Since the square of the Lax vector is non-negative, all roots are complex-conjugated pairs. Furthermore, due to the electron-hole symmetry of the problem, it is easy to show that if the initial texture is particle-hole symmetric then if y is a Lax root −y is also a root. We choose the initial superconducting order parameter, ∆, to be real. Also, this property is preserved at all times due to particle-hole symmetry. As previously discussed [2], Lax roots provide information on the frequency spectrum in the Fourier transform of ∆(t) after a quench. Starting from the equilibrium pseudospin texture, a dense distribution of Lax roots appears along the real axis which, in the thermodynamic limit, define the continuous part of the spectrum. For t → ∞ this contribution vanishes in the Fourier transform of ∆(t) and only isolated pairs of complex-conjugated Lax roots (bound states) contribute corresponding to persistent oscillations. It has been shown [2] that the number of discrete frequencies k in the Fourier transform of ∆(t) is equal to m − 1. Here, m is the number of isolated pairs of complex-conjugated Lax roots. Thus, ∆ (t) shows persistent oscillations at long times To complement the numerical results, we constructed the Lax vector and computed its roots using the equilibrium pseudospin texture for λ i as the initial condition. For non-critical quenches (λ f λ i > 0) the Lax analysis was done in Refs. [1,2] and the possible isolated pairs of Lax roots are always purely imaginary: i) In the synchronized regime, there are two pairs of isolated roots (m = 2) namely, y = ±iu 1 and y = ±iu 2 that give information on the amplitude of persistent oscillation as ∆ ± = (u 1 ±u 2 )∆ i . In this case, ∆(t) shows only one fundamental frequency (k = m − 1 = 1) corresponding to the lower edge Higgs mode 2∆. ii) In the damped phase there is only one isolated pair of Lax roots (m = 1) namely, y = ±i∆ ∞ indicating there are not persistent oscillations (k = m − 1 = 0) but damped. iii) In the overdamped regime, the latter pair of roots collapse to the origin (∆ ∞ = 0) so no isolated roots are present (m = 0).
We now extend these results to the critical quench (CQ, λ f λ i < 0) considered in the main text. In this case, the Lax roots are complex numbers with a finite real part, giving precise information not only on the amplitude but also the oscillation frequency of ∆(t). Both in the case of the constant and the graphene like-DOS, the roots can be written as y = ±(u r ± iu i ) with the real part providing information on the upper edge Higgs mode frequency, ω u = 2u r and the imaginary part yielding the oscillation amplitude ∆ ± = ±2u i .
For the graphene-like DOS and λ f λ i > 0, there is a regime where 2∆ ≠ 0 for which the isolated Lax roots are purely imaginary, mimicking the non-ZOPA synchronized phase of the constant DOS model. Still in the synchronized regime, where the Dirac-Higgs mode emerges, the isolated Lax roots acquire a finite real part defining the mode frequency as ω D = 2u r and the amplitudes of oscillations ∆ ± = ±2u i . Figure S1 shows representative dynamics ∆(t) and their Fourier transform (FT) for each regime shown in the Fig. 1 (a) of the main text for the flat DOS case. The color of each curve matches the colored dots in Fig. 1 (a) of the main text encoding the simulation parameters. Panels (a) and (b) correspond to non-critical quenches λ f λ i > 0 while (c) and (d) correspond to the CQ. In the top panels, we can see examples of the behavior in the three dynamical phases previously studied in Ref. [3] for λ f λ i > 0: persistent oscillations (blue), damped oscillations (orange) and overdamped dynamics (red). The former corresponds to the lower-edge Higgs mode with frequency ω l = 2∆ as can be seen in panel (b) where also a second harmonic peak at 4∆ appears. In the damped regime, a peak at frequency 2∆ ∞ is also resolved due to the finite time window we consider computing the FT. For the CQ (panels (c) and (d)), similar to the graphene-like DOS case discussed in the main text, ∆(t) shows persistent oscillations between ∆ ± with the upper-edge Higgs mode frequency ω u . Figure Fig. S2 (a) shows the pseudospin resolved FT of the dynamics in the synchronized phase with λ f λ i > 0. The large black dots correspond to the Larmor frequency The quench parameters δ used in each case, have been pointed out in Fig. 1 of main text with colored dots. The FT were computed considering a time window t∆i ∈ [5,50].

SUPPLEMENTARY FIGURES
Ω L = 2 ξ 2 k +∆ 2 . In addition, the spectrum shows Floquet satellites at Ω L ± ω l . The vertical features show that all pseudospins respond with the Higgs frequency 2∆ and high harmonics. Figure S2 (c) shows the same information in the damped regime. Each pseudospin oscillates with an effective Larmor frequency 2 ξ 2 k + ∆ 2 ∞ which introduces dephasing giving rise to damped oscillations in ∆(t) [c.f. Fig. S1 (a)]. The vertical line at ω = 2∆ ∞ is due to the finite time window in which the FT is performed. In the long time limit, its spectral weight vanishes, indicating the absence of persistent oscillations.
By decreasing enough the final coupling, the system enters into the gapless phase, where each pseudospin oscillates with its own bare Larmor frequency 2ξ k [panel (d)].
In the case of the CQ [panel (b)] the pseudospin resolved FT shows both the individual Larmor frequency, 2ξ k and a vertical feature in the upper edge corresponding to the upper-edge Higgs mode discussed in the main text for a constant DOS. The frequency is ω u ≈ W = 40∆ i and the intensity increase with ξ k indicating that high energy pseudospin participate with larger amplitude.