Generation and decay of Higgs mode in a strongly interacting Fermi gas

We investigate the life cycle of the large amplitude Higgs mode in strongly interacting superfluid Fermi gas. Through numerical simulations with time-dependent density functional theory and the technique of the interaction quench, we verify the previous theoretical predictions on the mode’s frequency. Next, we demonstrate that the mode is dynamically unstable against external perturbation and qualitatively examine the emerging state after the mode decays. The post-decay state is characterized by spatial fluctuations of the order parameter and density at scales comparable to the superfluid coherence length scale. We identify similarities with FFLO states, which become more prominent at higher dimensionalities and nonzero spin imbalances.


I. INTRODUCTION
Spontaneous breaking of continuous symmetries inevitably leads to the appearance of associated collective modes of a physical system.In the case of U(1) symmetry these are either gapless Goldstone mode, associated with oscillation of the phase of the order parameter, or the Higgs mode, which is associated with oscillations of the amplitude (for this reason, it is also called the amplitude mode).The Mexican-hat-like potential is typically drawn to visualize these modes: the Goldstone mode corresponds to motion around the hat in the minimum energy valley.In contrast, the Higgs mode is related to oscillations around the minimum in the perpendicular direction and consequently has much higher energy.Many studies related to the Higgs mode have been conducted in condensed matter (see eg. [1] for review).In recent years, advances in cooling techniques also allowed investigating this phenomenon in highly controllable environments of ultracold atoms [2,3], as well as developments in THz spectroscopy for superconductors [4][5][6][7][8][9].Simultaneously, many theoretical considerations have been presented concerning the Higgs mode in ultracold atomic systems: properties of small and large amplitude oscillations [10][11][12][13][14], impact of the trapping potential on the mode characteristics [15][16][17], angle-resolved photoemission spectroscopy (ARPES) [18][19][20] and various other ways to induce collective modes [21][22][23][24].Moreover, many analytical results have been derived in the weak-coupling limit [13,[25][26][27][28].The most well-known property is that the frequency of small-amplitude oscillations is related to the equilibrium value of the pairing gap by Ω = 2∆.Through this, one can deduce the pairing gap by measuring the Higgs mode frequency [3,29].
The most widely discussed method of generation of the amplitude mode is through the interaction quench [13,[30][31][32][33][34][35][36].This idea exploits the fact that the equilibrium value of the pairing field (order parameter) depends on the interaction strength, which is typically characterized by the dimensionless value ak F , where a is the scattering length and k F = (3π 2 n) 1/3 is the Fermi wave-vector for a gas of density n.Thus, by preparing the system in the ground state for a selected initial interaction strength a (initial) k F and then changing the interaction regime rapidly to the new value a (final) k F , one can induce oscillations of the order parameter within the range specified by [∆(a (initial) k F ), ∆(a (final) k F )]. Theoretical studies of such scenarios have been performed by use of the Bogoliubov de-Gennes (mean-field) method, which is justified for weakly interacting systems with attractive interaction: |ak F | 1 and a < 0. On the other hand, Fermi superfluids produced in laboratories are typically strongly interacting |ak F | 1.A separate aspect concerns the stability of the amplitude (Higgs) mode.In the case of weak couplings and small amplitude oscillations, by utilizing linear response theory, Dzero, Yuzbashyan, and Altshuler demonstrated that the mode is unstable [37].The mode decays into a new state with spatially nonuniform order parameter.Numerical simulations performed for the attractive Fermi-Hubbard model [30,38] indicated that the mode is indeed dynamically unstable, which was consistent with the theoretical predictions.
The purpose of the presented article is twofold.First, we examine the stability of the Higgs mode in strongly interacting Fermi gases -the systems that are of particular experimental interest.By means of time-dependent density functional theory, we simulate the whole process: from the generation of the Higgs mode via the interaction quench, through the decay dynamics until the equilibration of the final (post-decay) state.Second, we investigate the properties of the post-decay state.We show that it is characterized by spatially inhomogeneous order parameter, bearing similarities to Larkin-Ovchinnikov (LO) [39] or Fulde-Ferrel (FF) [40] type oscillations.Upon introducing another degree of freedom to the problem in the form of spin imbalance, the post-decay state indeed acquires properties as predicted for spin-imbalanced systems [41][42][43][44], although the final state corresponds to an excited state.

II. METHOD AND FRAMEWORK
We employ the Density Functional Theory (DFT) formalism to study the properties of the amplitude mode.Many variants of DFT methods exist; here we utilize the one known as Superfluid Local Density Approximation (SLDA), which has been specifically designed to simulate fermionic superfluid systems where interparticle interactions have short range [45,46].It is a microscopic theory where the system is described in terms of fermionic quasi-particles, defined through Bogoliubov amplitudes: ϕ nσ (r, t) = [u nσ (r, t), v n−σ (r, t)] T , where σ = ± indicates spin projection.Precisely, v nσ and u nσ stand for amplitude probabilities for n-th state to be occupied by a particle and a hole with spin projection σ.Formally, the equations of motion have the same structure as Bogoliubov-de Gennes (BdG), however they are obtained as a result of energy minimization expressed as a functional of several densities characterizing fermionic superfluid: The energy density E depends on the following densities (we skip position and time dependence for brevity): • normal densities • kinetic densities • current densities • anomalous density The Fermi-Dirac distribution function f β (E) = 1/[1 + exp (βE)] accounts for temperature T (β −1 = k B T ), although in this work we mainly focus on the T → 0 limit.The densities are constructed from the quasiparticle amplitudes up to a certain energy cut-off E c .This procedure accounts for regularization of divergences appearing due to short-range interactions [47].The minimization yields static equations for the quasi-particle amplitudes (hereafter we use units where m = = 1) Note that there are two sets of equations for σ = ±.The time-dependent variant is obtained by replacing E nσ → i ∂ ∂t .Chemical potentials, denoted as µ σ , are used to control particle numbers N σ = n σ (r)dr associated with a given spin component.The time-dependent formulation conserves the total number of particles of each species N σ (t) = const.The single particle Hamiltonian h σ and the pairing potential ∆ are defined through functional derivatives of the energy: The pairing potential is a complex field ∆(r, t) that serves as the order parameter.In addition we have added an external pairing potential δ∆ ext that will be used to generate a perturbation when studying the stability of the amplitude mode.Due to symmetry of Eqs. ( 6) one needs effectively a solution for either σ = + or σ = −.The other one can be obtained through symmetry transformation (see eg. [48]).The energy functional defines the system.In almost all cases devoted to the studies of the Higgs mode in Fermi systems, the mean-field Bogoliubov-de Gennes equations were used.They are valid for weakly interacting Fermi systems, and upon specific choice of the energy density functional the formalism presented here becomes identical to the BdG.The single particle Hamiltonian and the pairing potential are then of the form: The advantage of a description based on DFT is that it allows to incorporate so-called "beyond mean-field" effects, while keeping the numerical complexity at a level comparable to the mean-field methods.Here we use the functional known as SLDAE [49], which has the generic form where n = n + +n − and e F = is the associated Fermi energy.The dimensionless coupling functions A, B and C are constructed in such a way to ensure the correct reproduction of thermodynamics quantities for a uniform system (equation of state, chemical potential, pairing gap, effective mass) over the entire interaction regime from weak (BCS regime) to strong (unitary Fermi gas regime).Their explicit forms are given in paper [49].The SLDAE functional has been designed to provide accurate results for spin-symmetric systems (N + = N − ), however it can be extended to spin-imbalanced systems as well.In this paper we will show a few selected results for spin-imbalanced systems (N + = N − ).Since we are interested in qualitative aspects for such cases, we use the functional without further adjustments (in the same spirit as it is usually done with BdG approach).Note that further optimizations of the functional towards spin-imbalanced systems can be performed in a similar way as in Refs.[46,50], but from the perspective of this work it is not needed.
A quantity of special importance in the context of this work is the order parameter ∆(r, t).It is a complex function, and both attributes (magnitude and phase) carry physical information.The amplitude/Higgs mode corresponds to uniform oscillations of the absolute value across the whole system.Thus, we consider a uniform system (no external trapping potential) in which we induce the mode through the interaction quench method (see next section for details).At the beginning of each simulation the density distributions (2)-( 5) are uniform and may depend on time only.In order to study the stability, we solve the equations of motion on spatial grid of size N x × N y × N z with lattice spacing dx = dy = dz = 1 (definition of the length unit).The lattice spacing sets a natural cut-off energy scale 2 .Most of the calculations will be presented for the constrained case, where we impose translation symmetries along y and z directions.It implies that quasiparticle wave-functions have a generic structure that can be written as where k y and k z take discrete values of multiples of 2π/N y and 2π/N z , respectively.Under this assumption during the evolution (breaking of translation symmetry), the densities and the order parameter can acquire dependence along the x direction only.To test the robustness of conclusions obtained from quasi-1D simulations, we have also performed additional simulations for quasi-2D case: as well as fully unconstrained simulations in 3D.We specifically focus on the large amplitude mode induced in the strongly interacting regime (|ak F | → ∞), as the stability for this case has not been discussed in the literature to date.For the computation we use W-SLDA Toolkit [51][52][53].All technical parameters, needed to restore the calculations presented below are included into reproducibility packs [54].They contain also detailed information about the numerical setup and computation process.

III. GENERATION OF THE HIGGS MODE VIA THE INTERACTION QUENCH
Rapid quenching of the interaction strength a (initial) k F → a (final) k F is the standard method of inducing the Higgs mode [10,15,17,25,29,55].After the generation of the Higgs mode, it is interesting to test its stability with respect to random perturbation.The perturbation is realized by adding the weak external potential δ∆ ext (r, t).Consequently, in our simulations the emerging time evolution can be divided into 4 parts.
The homogeneous gas characterized by s-wave scattering length a 0 and (normal) density n 0 = k 3 F /(3π 2 ) is in its ground state.We consider spin-symmetric simulations (N ↑ = N ↓ ), and only in some cases we start simulations from states corresponding to spin-imbalanced (N ↑ = N ↓ ) systems.

Quenching of interaction:
The initial state is driven adiabatically towards the new interaction strength a 0 k F → a 1 k F in the time period t 0 < t ≤ t 1 , and then quasi-instantaneously driven back to the initial value a 1 k F → a 0 k F in the time interval t 1 < t ≤ t q , where (t 1 − t 0 ) (t q − t 1 ).Explicitly, the time-dependent s-wave scattering length is parametrized as follows: such that a(t) and the first derivative ∂ t a(t) are continuous, leading to t 1 − t 0 = π/ω 1 and t q − t 1 = π/ω 2 .Using this parametrization, the limit ω 1 e F corresponds to an adiabatic process, i.e. a small energy transfer to the system.On the other hand ω 2 e F generates the rapid quench.
Higgs oscillations: (t q < t ≤ t d ) The interaction quench induces the Higgs oscillations with wave vector k = 0 (uniform oscillations) and frequency Ω.Thus, observables do not exhibit position dependence, for example for the pairing field we have |∆(x, t)| → |∆(t)|.In order to investigate the stability of this mode, we add a weak perturbation to the system δ∆ ext (r, t) once the oscillations have developed.We have tested the stability of the Higgs mode with respect to two types of perturbation: local in momentum space and spatially localized.The explicit forms of perturbations read as: where e F (weak perturbation) and σe F 1 (applied at much shorter times scale as the expected time scale of the Higgs oscillation).In numerical realization the Dirac function δ(x) is approximated by narrow Gaussian function.

Higgs mode decay: (t d < t)
After some time t d we find that the Higgs mode decays and an inhomogeneous phase emerges.Namely, the observables depend now on position and time, for example for the quasi-1D calculations: n(x, t) = n 0 and ∆(x, t) = ∆(t).
In Fig. 1 we present an example time evolution of the absolute value of the order parameter for selected point in space |∆(x = 0, t)|.The described stages of the numerical experiment can be easily identified.
As it will be shown in this paper, the Higgs mode is dynamically unstable, i.e small perturbations amplify and destroy the mode.We note that even in the case of δ∆ ext (x, t) = 0 (no external perturbation) the mode decays spontaneously after some time t max .In such case, the mode destabilizes due to numerical noise, since we integrate the equations of motion numerically with some precision.When measuring the lifetime of the Higgs mode, defined as the time from the weak perturbation to the decay, we limit our considerations to time intervals t < t max to rule out complications related to imperfections of the numerical integration.We recognize that in the case of studies of the dynamical stability of physical phenomena by means of computational physics, code quality is an important issue.We note that the accuracy of time integrator implemented within W-SLDA, which is multistep Adams-Bashforth-Moulton of 5th order, has been tested in work [56].Namely, the dependence of the results with respect to the perturbations of the initial states, for cases where no dynamical instability is expected, was tested.Moreover, the W-SLDA Toolkit has been applied to a variety of problems.Results were documented in papers [42,44,51,[57][58][59][60].We have confirmed stability of the results with respect to size of the integration time step dt (all presented here results were obtained with dt = 0.0025e −1 F ).We have also checked the stability of the results with respect to the lattice size.For the quasi-1D calculations (see sections V and VI), we have checked that conclusions are stable if we vary lattice size from 256 to 1024.Thus we regard the code as well tested, and the observed instability in the simulations is expected to be due to intrinsic properties of the studied problem, not due to artifacts generated by the numerical implementation.

IV. PROPERTIES OF THE HIGGS MODE OSCILLATIONS
We start by examining the properties of the amplitude mode.The aim of this section is to demonstrate that the interaction quench method induces the Higgs mode, having the well known properties.The mode is typically regarded as as oscillating system in the effective potential in the form of a Mexican hat: The solutions for a classical particle moving in such potential are well known [1,61].In the case of small amplitude oscillations around the minimum of the potential ∆ 0 = ±µ/ √ we have: where ω is the oscillator's frequency, φ is its phase, and A(t) is the amplitude, which in general may depend on time [13,26].Moreover, for the small amplitude Higgs mode, the oscillation frequency is related to the equilibrium value of the pairing gap ω = 2∆ 0 .The oscillation of arbitrary amplitude are expressed by delta amplitude Jacobi elliptic function dn: where ∆(t 0 ) = ∆ 0 = 0.4425e F and the elliptic modulus parameter k is related to maximum and minimum values of the paring field during the oscillations: In Fig. 2 we compare the numerically extracted time dependence of the pairing field ∆(t) after the interaction quench with analytic predictions.The time-dependent variant of Eq. ( 6) has been solved numerically on a lattice of size 256 × 32 × 32 within quasi-1D geometry, see Eq. (11).In both cases, small and large amplitude oscillations, we observe a close match with analytic predictions.In particular, the well known property of the small amplitude Higgs mode, ω = 2∆ 0 , has been reproduced with very good accuracy.The large amplitude mode oscillates with lower frequency Ω ≈ ∆ 0 and the parameter k = 0.9935 provided by the fit also agrees well with the expected value k = 0.9952 obtained by applying Eq. 17.These findings are fully consistent with previous works, such as [10,13,26].

V. LIFETIME OF THE HIGGS MODE
The Higgs mode is attributed to oscillations of the amplitude of the order parameter, while the density of the system remains unchanged.Indeed, when the mode is induced in a uniform system, n(r, t) stays constant.The appearance of inhomogenities in the density is associated with the decay of the Higgs mode.To quantify the magnitude of spatial variations of the density we compute: FIG. 3. Time evolution of the normalized standard deviation of density σ[n](t) and amplitude of the order parameter in selected point in space ∆(0, t) after the applied perturbation in form ∆ (1) . In these simulations the large amplitude mode is induced by rapidly quenching the interaction from a1kF = −0.1 to a0kF = −10.The perturbation function parameters are /eF = 10 −3 and σeF = 0.5.Each panel corresponds to various wave-vectors k, indicated on the respective labels.Calculations were executed with quasi-1D geometry on a lattice of size 512 × 32 × 32.
The decay time t d is defined as the time for which the normalized standard deviation exceeds a threshold value σ[n](t > t d ) > γ = 0.01.The dynamics of real system is always attributed by stochastic fluctuations (for example, due to thermal effects) that may amplify in time in case of unstable modes.In the numerical scenario, we introduce the weak external perturbation δ∆ ext and check if induced fluctuations amplify or decay.In the case of the mode studied here, we find that it is unstable.One can also induce fluctuations by adding weak external potential that couples to the density.The final qualitative result does not depend on the perturbation method.We define the lifetime of the Higgs mode as the time from the external perturbation t p to the decay t d .In Fig. 3 we have shown induced decay by the perturbation δ∆ (1) ext (x, t).We clearly observe a sensitivity of the decay time with respect to the spatial modulation wavelength λ ∼ 1/k.The perturbation with k = 0 (top panel) corresponds to preserving translational symmetry and therefore it does not have any impact on the dynamics of the system.The mode starts decaying due to the amplification of numerical errors only, after a time t max e F 150. On the other hand, perturbations corresponding to k > 0 induce a decay whose onset time t d depends on the wavelength of the spatial modulation, see Fig. 4(a).Clearly, there exist optimal value of k (resonance) that trigger the instability almost immediately.The location of the optimal value depends on the interaction strength quench.For tested cases it is located around k/k F ≈ 0.2 − 0.4.In Fig. 4(b) we show the Fourier spectra of the order parameter after applying the perturbation, but before the decay.It is seen that the uniform Higgs mode (k = 0) first decays into two modes with +k d and −k d .The diagram representing schematically the decay of small amplitude mode is shown in inset of Fig. 4(b).
Instead of using a perturbation with well defined wave length 2π/k, one can instead use a spatially localized perturbation such as δ∆ (2) ext ∼ δ(x) = dk 2π e ikx .In numerical realization, we model the δ(x) function by a narrow Gaussian.Similarly, we observe that the amplitude mode decays and the state after the decay is nonuniform.In Fig. 5 we compare spatio-temporal evolution of the order parameter after applying different types of perturbation.The decay dynamics depend on the perturbation type; however there are gross properties that remain unchanged.Namely, in the Fourier spectra of the order parameter |∆(k, t)| we observe dominant wave-vectors ±k d into which the The results obtained for the δ-kick need extra clarification (top right panel of Fig. 5).In general, one can expect that the fluctuation induced by a localized perturbation should propagate with speed not exceeding the speed of sound.Indeed in such case, we should see the perturbation only in the causality cone.It is visible only for times (t − t p )e F 75.For later times one may conclude from the graph that the causality is violated.Such a conclusion is incorrect.The observed effect is the numerical artifact.Namely, we model δ(x) function by narrow gaussian, as the numerical scheme assumes that all derivatives are continuous and smooth.This implies that our perturbation is not localized in practice but affects the system in the whole spatial domain.
Dzero, Yuzbashyan & Altshuler [37] demonstrated that the small amplitude Higgs mode is unstable.Their result has been obtained within mean-field (BdG) treatment, which is valid for weakly interacting superfluid Fermi gases.In such case, the presence of dynamical instability was demonstrated analytically within the framework of the linear response theory.In the present study, we have released these constraints by considering large amplitude and strongly interacting regime and we have demonstrated the instability of the mode.Namely, we have shown that similarly to the small amplitude regime, the mode decays predominantly into a state where the spatial modulation of the pairing field is given by ∼ 1/k d , i.e. ∆(x, t) = ∆(t) + A(t)[e −ik d x + e ik d x ], where A(t) is the amplitude of the perturbation which grows rapidly in time.Consequently, this implies that after the decay, the order parameter exhibits fluctuations of the type ∆(x) ∼ cos(k d x).Indeed, in Fig. 5 (right panel) one can clearly see lines where |∆(x, t)| = 0 (called nodal lines) for later times.These are lines where the order parameter changes sign.The spatial fluctuations of ∆(x) ∼ cos(kx) are typically regarded as a fingerprint of Larkin-Ovchinnikov (LO) state [39], if they emerge in spin-imbalanced system in a ground state.None of these requirements is satisfied in the cases discussed so far.However, it is interesting to note that the state emerging from the decay process shares similarities with the exotic type of superfluidity.This issue will be discussed in more detail in next section.

VI. PROPERTIES OF THE POST-DECAY STATE
We now consider the properties of the state after the decay.To better visualize the typical characteristics of the post-decay state, we switch to a quasi-2D geometry, ansatz (12).The calculations were executed on a 64 × 64 × 16 lattice.Qualitatively, we observe the same phenomena as reported above for quasi-1D cases: the mode decays to a state that is attributed by the non-homogeneous distribution of the order parameter.The structure of ∆(x, y) consists of regions where the phase changes by π, which are separated from each other by nodal lines defined as |∆(x, y)| = 0.The development of inhomogeneities of the pairing field is accompanied with spatial modulations of density distribution which are induced simultaneously.One needs to emphasize that it is a transient state, which evolves on much longer time scale (see Fig. 6(a,b) with comparison at different times from the perturbation).For   later times (larger time scale than reachable in numerical computation), we expect the system to thermalize again to the state with uniform density distribution and the order parameter.
Interestingly, the transient state that emerges from the decay shares similarities to states as predicted for spinimbalanced systems [41][42][43].The main difference is that if the imbalance is present (N ↑ = N ↓ ), the majority particles tend to accumulate in the nodal regions if such are developed in the system.It can be understood from the point of view of the system's energetics: not all particles can form Cooper pairs (superfluid component), and the energetically lowest price is paid if we store them close to the nodal lines where the pairing correlations vanish.The local spin polarization, defined as p(r) = [n ↑ (r) − n ↓ (r)]/[n ↑ (r) + n ↓ (r)], acts as a stabilizer and converts the states consisting of nodal lines into meta-stable structures.Such structures, stabilized by the spin polarization, were recently the subject of studies where they were called ferrons [42] or ring solitons [64].The size of these grains depend on the global population imbalance of the system P = For small imbalances P 10% their sizes can by of the order ∼ 10ξ, and decrease to ∼ ξ as the imbalance increases, and finally converting to well known LOFF state [41].Suppose we induce the Higgs mode in the spin-imbalanced system and let it decay, precisely in the same manner as we did for the spin-symmetric case.Then, we indeed find that the post-decay state is stabilized by the spin-polarization, which at time scales larger then the time scale associated with the Higgs mode decay evolves towards state consisting of many spin-polarized impurities, see Fig. 6(c,d).Note that just after the decay, the typical modulation wave-length for the order parameter (and also for the density) is close to the the wave-length of the fastest decaying mode.Precisely, the average size of the structures we observed in Figs 6 (a) and (c) translate into wave-vectors k/k F = 2π/λk F ≈ 0.2 for both P = 0% and P = 10% at the onset of the oscillation.Clearly, these values to the minimum location in Fig. 4. Once the disordered structure is set by the decayed Higgs mode properties, it is driven to the new configuration by effects related to dynamics of the spin polarization.They operate on longer time scales, and in large time scales (te F ∼ 10 3 ) the polarized system evolves towards to the known from past studies configurations [41][42][43].
Finally, we demonstrate that the observations hold for fully 3D cases, as shown in Fig. 6(e).These results were obtained for a lattice of size 64 × 64 × 64.Similarly to the quasi-2D case, we introduce the spin-imbalance to the  [63] system in order to stabilize the final state.The post-decay state, at large times, consists of many bubble-like structures.Inside each bubble, the phase of the order parameter is shifted by π.The overall properties of the post-decay state remain fully consistent with results obtained for simplified geometries (quasi-1D and quasi-2D).The insensitivity of the results at the quantitative level with respect to the computation dimensionality demonstrates the robustness of derived conclusions.

VII. CONCLUSIONS
We have studied the life cycle of an ultracold atomic gas after the interaction quench from weak to strong coupling.The quench applied to uniform system induces the large amplitude Higgs/amplitude mode.This mode turns out to be dynamically unstable and it decays to a state with spontaneously broken spatial symmetry.It is of a different type than discussed so far in literature [13,26,65], where the decay of amplitude of Higgs oscillations was inspected while the uniformity of the system was maintained.
The relation between frequency of the small amplitude Higgs mode and strength of the pairing correlations ( ω = 2∆ 0 ) makes it a very valuable tool from the perspective of experimental measurements.For this reason, recent works (such as [55]) focus on the stabilization of this mode.Here, and also in [37], we demonstrate that the nonuniform state emerging from the decay of the mode can provide access to new class of states with order parameter that is periodically modulated in space.Such states are frequently associated with exotic types of superfluidity, generically related to Fulde-Ferrell-Larkin-Ovchinnikov phase -a state that is routinely discusses in the context of spin-imbalanced systems.Indeed, when introducing a spin-imbalance to the configuration, the states emerging from the decay consist of many spin-polarized droplets [41,42].Surprisingly, experiments devoted to Higgs modes in strongly interacting Fermi gases may contribute to deeper understanding of exotic states emerging in superfluids.
Tunable ultracold atomic gases confined in box traps can be potentially used to verify the findings of this work [66].Such setups have already demonstrated their high capabilities for detecting the system's nonuniformities.For example, there were successfully used to visualize density modulations due to the propagation of sound waves up to wavelengths of the order 0.1 k F [67] -the resolution that should be sufficient to detect predicted here inhomogeneities in the postdecayed state.

FIG. 1 .
FIG.1.The numerical experiment investigated in this work.The Higgs mode is induced in the uniform system by quenching the interaction strength in the time interval t0(= 0) < t ≤ tq.The mode persists for some time (tq < t ≤ t d ), after which it decays into a nonuniform state (t > t d ).The red line presents the absolute value of the pairing field for a selected point at the center of the lattice, |∆(x = 0, t)|, normalized to the initial value.The presented case corresponds to a large amplitude Higgs mode.

FIG. 4 .
FIG. 4. Panel (a): Decay time as a function of the modulation wave vector k/kF for different stregths of the quench; values on the label indicate the values of a1kF and a0kF = −10 is fixed.Panel (b): Fourier spectra of the order parameter just before the decay induced by the external modulation k/kF = 0.4 for the case with a1kF = −0.1.Inset: The decay diagram for the small amplitude Higgs mode.

FIG. 5 .
FIG.5.Evolution of the order parameter in coordinate space (top row) and its Fourier spectra (bottom row) for perturbation δ∆

( 2 )
ext (x, t) (right column).The ∆(k = 0, t) component was removed from the Fourier spectra for better visibility.The amplitude of the perturbation was /eF = 10 −3 and evolution is shown only after the perturbation is applied.Maps generated by Matplotlib library[62].

FIG. 6 .
FIG. 6.(a) and (b): 2D configurations for P = N ↑ −N ↓ N ↑ +N ↓ = 0%, respectively after a short and long time from the perturbation is applied.(c) and (d): 2D configurations for P = 10%.Panels on the bottom row (rainbow palette) show the absolute value of the order parameter, normalized to eF ; panels on the top row (heat palette) show the phase of the order parameter.e) 3D configuration for P = 5%.Upper colorbar indicates local spin polarization p(r) on the right panel; colorbar on the right shows the absolute value of the order parameter |∆(r)|, normalized to eF .Inner inset shows the contour surface for which the order parameter is zero.The dimensionless size of the box is LkF = 64 in each direction.Maps and 3D view created by VisIt software[63]