ABSTRACT 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.

continuous symmetries inevitably leads to the appearance of associated collective modes of a physical system. In the case of a 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 e.g.1 for review). In recent years, advances in cooling techniques also allowed investigating this phenomenon in highly controllable environments of ultracold atoms2,

3, as well as developments in THz spectroscopy for superconductors4,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 oscillations10,11,12,13,14, impact of the trapping potential on the mode characteristics15,16,17, angle-resolved photoemission

spectroscopy (ARPES)18,19,20 and various other ways to induce collective modes21,22,23,24. Moreover, many analytical results have been derived in the weak-coupling limit13, 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 \(\hbar \Omega =2\Delta\). Through this, one can

deduce the pairing gap by measuring the Higgs mode frequency3, 29. The most widely discussed method of generation of the amplitude mode, in ultracold atomic systems, is through the

interaction quench13, 30,31,32,33,34,35,36. This idea relies on 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_{\textrm{F}}\), where _a_ is the scattering length and \(k_{\textrm{F}}=(3\pi ^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^{(\text {initial})}k_{\textrm{F}}\) and then changing the interaction regime

rapidly to the new value \(a^{(\text {final})}k_{\textrm{F}}\), one can induce oscillations of the order parameter within the range specified by \([\Delta (a^{(\text

{initial})}k_{\textrm{F}}),\Delta (a^{(\text {final})}k_{\textrm{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_{\textrm{F}}|\ll 1\) and \(a<0\). On the other hand, Fermi superfluids produced in laboratories are

typically strongly interacting \(|ak_{\textrm{F}}|\gg 1\). A separate aspect concerns the stability of the amplitude (Higgs) mode. In the case of weak couplings and small amplitude

oscillations, Dzero, Yuzbashyan, and Altshuler demonstrated by utilizing linear response theory that the mode is unstable37. The mode decays into a new state with spatially nonuniform order

parameter. Numerical simulations performed for the attractive Fermi–Hubbard model30, 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. Secondly, 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 systems41,42,43,44, although the final state corresponds to an excited state. 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 range45, 46. It is a

microscopic theory where the system is described in terms of fermionic quasi-particles, defined through Bogoliubov amplitudes: \(\varphi _{n \sigma }(\varvec{r},t)=[u_{n \sigma

}(\varvec{r},t), v_{n -\sigma }(\varvec{r},t)]^T\), where \(\sigma =\pm\) indicates spin projection. Precisely, \(v_{n \sigma }\) and \(u_{n \sigma }\) stand for amplitude probabilities for

_n_-th state to be occupied by a particle and a hole with spin projection \(\sigma\). 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: $$\begin{aligned} E = \int {\mathcal {E}}(n_\sigma ,\tau

_\sigma ,\varvec{j}_\sigma ,\nu _\sigma )d\varvec{r}. \end{aligned}$$ (1) The energy density \({\mathcal {E}}\) depends on the following densities (we skip position and time dependence for

brevity): * Normal densities $$\begin{aligned} n_{\sigma } = \sum _{|E_{n \sigma }|<E_c}|{v_{n \sigma }}|^2 f_{\beta }(-E_{n \sigma }), \end{aligned}$$ (2) * Kinetic densities

$$\begin{aligned} \tau _{\sigma } = \sum _{|E_{n \sigma }|<E_c}|{\nabla v_{n \sigma }}|^2 f_{\beta }(-E_{n \sigma }), \end{aligned}$$ (3) * Current densities $$\begin{aligned}

\varvec{j}_{\sigma } = \sum _{|E_{n \sigma }|<E_c}\text {Im}[v_{n \sigma }\nabla v_{n \sigma }^{*}] \,f_{\beta }(-E_{n \sigma }), \end{aligned}$$ (4) * Anomalous density $$\begin{aligned}

\nu _{+-} = \nu = \sum _{|E_{n \sigma }|<E_c} u_{n -}v_{n +}^{*} \frac{ f_{\beta }(-E_{n +} )-f_{\beta }(E_{n -})}{2}. \end{aligned}$$ (5) The Fermi–Dirac distribution function

\(f_{\beta }(E) = 1/[1 + \exp \left( \beta E \right) ]\) accounts for temperature _T_ (\(\beta ^{-1}=k_B T\)), although in this work we mainly focus on the \(T\rightarrow 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

interactions47. The minimization yields static equations for the quasi-particle amplitudes (hereafter we use units where \(m=\hbar =1\)) $$\begin{aligned} \begin{pmatrix}h_{\sigma }-\mu

_{\sigma } &{} \sigma (\Delta + \delta \Delta _{\text {ext}}) \\ \sigma (\Delta ^*+ \delta \Delta ^*_{\text {ext}}) &{} -h^*_{-\sigma }+\mu _{-\sigma }\end{pmatrix}

\begin{pmatrix}u_{n \sigma } \\ v_{n -\sigma }\end{pmatrix}= E_{n \sigma }\begin{pmatrix}u_{n \sigma } \\ v_{n -\sigma }\end{pmatrix}. \end{aligned}$$ (6) Note that there are two sets of

equations for \(\sigma =\pm\). The time-dependent variant is obtained by replacing \(E_{n\sigma }\rightarrow i\frac{\partial }{\partial t}\). Chemical potentials, denoted as \(\mu _\sigma\),

are used to control particle numbers \(N_\sigma = \int n_\sigma \,(\varvec{r}) d\varvec{r}\) associated with a given spin component. The time-dependent formulation conserves the total

number of particles of each species \(N_\sigma (t)=\text {const}\). The single particle Hamiltonian \(h_\sigma\) and the pairing potential \(\Delta\) are defined through functional

derivatives of the energy: $$\begin{aligned} h_{\sigma } = -\varvec{\nabla }\frac{\delta E}{\delta \tau _\sigma }\varvec{\nabla } + \frac{\delta E}{\delta n_\sigma }-\dfrac{i}{2}\left\{

\frac{\delta E}{\delta \varvec{j}_\sigma },\varvec{\nabla } \right\} ,\quad \Delta =-\frac{\delta E}{\delta \nu ^*}. \end{aligned}$$ (7) The pairing potential is a complex field \(\Delta

(\varvec{r},t)\) that serves as the order parameter. In addition we have added an external pairing potential \(\delta \Delta _{\text {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 \(\sigma =+\) or \(\sigma =-\). The other one can be obtained

through symmetry transformation (see e.g.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 BdG

equations were used. They are valid for weakly interacting Fermi systems, and upon specific choice of the energy density functional $$\begin{aligned} {\mathcal {E}} = \sum _{\sigma =\pm

}\frac{\tau _\sigma }{2}+4\pi a |\nu |^2 \end{aligned}$$ (8) the formalism presented here becomes identical to the BdG. The single particle Hamiltonian and the pairing potential are then of

the form: $$\begin{aligned} h_{\sigma } = -\frac{1}{2}\varvec{\nabla }^2,\quad \Delta =-4\pi a \nu . \end{aligned}$$ (9) 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 SLDAE49, which

has the generic form $$\begin{aligned} {\mathcal {E}} = \sum _{\sigma =\pm } \frac{A(ak_{\textrm{F}})}{2}\left( \tau _\sigma -\frac{\varvec{j}^2_\sigma }{n_\sigma }\right) +

\frac{3}{5}B(ak_{\textrm{F}})\,n\,e_{\textrm{F}}+ \frac{C(ak_{\textrm{F}}) }{n^{1/3}} |{\nu }|^2 + \sum _{\sigma =\pm } \frac{\varvec{j}^2_\sigma }{2n_\sigma }, \end{aligned}$$ (10) where

\(n=n_{+} + n_{-}\) and \(e_{\textrm{F}}=\frac{k_{\textrm{F}}^2}{2}=\frac{1}{2}[3\pi ^2n]^{2/3}\) 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 paper49. 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_{+}\ne 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 \(\Delta (\varvec{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 \times N_y \times N_z\) with lattice spacing

\(dx=dy=dz=1\) (definition of the length unit). The lattice spacing sets a natural cut-off energy scale \(E_c\approx \frac{k_c^2}{2}=\frac{\pi ^2}{2}\). Most of the calculations will be

presented for the constrained case, where we impose translational symmetries along _y_ and _z_ directions. It implies that quasiparticle wave-functions have a generic structure that can be

written as $$\begin{aligned} \psi _n(\varvec{r},t) \rightarrow \psi _n(x,t)\frac{1}{\sqrt{N_y}}e^{ik_y y}\frac{1}{\sqrt{N_z}}e^{ik_z z}, \end{aligned}$$ (11) where \(k_y\) and \(k_z\) take

discrete values of multiples of \(2\pi /N_y\) and \(2\pi /N_z\), respectively. Under this assumption during the evolution (breaking of translational 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: $$\begin{aligned} \psi _n(\varvec{r},t) \rightarrow \psi _n(x,y,t)\frac{1}{\sqrt{N_z}}e^{ik_z z}, \end{aligned}$$ (12) as well as fully unconstrained simulations in 3D. We

specifically focus on the large amplitude mode induced in the strongly interacting regime (\(|ak_{\textrm{F}}|\rightarrow \infty\)), as the stability for this case has not been discussed in

the literature to date. For the computation we use W-SLDA Toolkit51,52,53. All technical parameters, needed to restore the calculations presented below are included into reproducibility

packs (See supplementary files added to this submission). They also contain detailed information about the numerical setup and computation process. GENERATION OF THE HIGGS MODE VIA THE

INTERACTION QUENCH Rapid quenching of the interaction strength \(a^{(\text {initial})}k_{\textrm{F}}\rightarrow a^{(\text {final})}k_{\textrm{F}}\) is the standard method of inducing the

Higgs mode10, 15, 17, 25, 29, 54. 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 \(\delta \Delta _{\text {ext}}(\varvec{r},t)\). Consequently, in our simulations the emerging time evolution can be divided into 4 parts. * _Initial

configuration:_ (\(-\infty < t \le t_0\)) The homogeneous gas characterized by _s_-wave scattering length \(a_0\) and (normal) density \(n_0 = k_{\textrm{F}}^3 / (3\pi ^2)\) is in its

ground state. We consider spin-symmetric simulations (\(N_{\uparrow }=N_{\downarrow }\)), and only in some cases we start simulations from states corresponding to spin-imbalanced

(\(N_{\uparrow }\ne N_{\downarrow }\)) systems. * _Quenching of interaction:_ (\(t_0< t \le t_q\)) The initial state is driven adiabatically towards the new interaction strength

\(a_0k_{\textrm{F}}\rightarrow a_1k_{\textrm{F}}\) in the time period \(t_0<t\le t_1\), and then quasi-instantaneously driven back to the initial value \(a_1k_{\textrm{F}}\rightarrow

a_0k_{\textrm{F}}\) in the time interval \(t_1<t\le t_q\), where \((t_1-t_0) \gg (t_q-t_1)\). Explicitly, the time-dependent _s_-wave scattering length is parametrized as follows:

$$\begin{aligned} a(t) = {\left\{ \begin{array}{ll} a_0 + \dfrac{a_1-a_0}{2}[1 - \cos \omega _1(t-t_0)], &{} t_0< t \le t_1 \\ a_1 + \dfrac{a_0-a_1}{2}[1 - \cos \omega _2(t-t_1)],

&{} t_1 < t \le t_q \end{array}\right. } \end{aligned}$$ (13) such that _a_(_t_) and the first derivative \(\partial _t a(t)\) are continuous, leading to \(t_1 - t_0 = {\pi }/{\omega

_1}\) and \(t_q - t_1 = {\pi }/{\omega _2}\). Using this parametrization, the limit \(\omega _1 \ll e_{\textrm{F}}\) corresponds to an adiabatic process, i.e. a small energy transfer to the

system. On the other hand \(\omega _2 \gtrsim e_{\textrm{F}}\) generates the rapid quench. * _Higgs oscillations:_ (\(t_q < t \le t_d\)) The interaction quench induces the Higgs

oscillations with wave vector \(k=0\) (uniform oscillations) and frequency \(\Omega\). Thus, observables do not exhibit position dependence: for example, for the pairing field we have

\(\vert \Delta (x,t) \vert \rightarrow \vert \Delta (t)\vert\). In order to investigate the stability of this mode, we add a weak perturbation to the system \(\delta \Delta _{\text

{ext}}(\varvec{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: * 1. \(\delta \Delta _\text {ext}^{(1)}(x,t) = \epsilon \exp (ikx) \exp (-(t-t_p)^2/2\sigma ^2)\) (induced decay by spatial

modulation) * 2. \(\delta \Delta _\text {ext}^{(2)}(x,t) = \epsilon \delta (x) \exp (-(t-t_p)^2/2\sigma ^2)\) (induced decay by Dirac-delta perturbation) where \(\epsilon \ll

e_{\textrm{F}}\) (weak perturbation) and \(\sigma e_{\textrm{F}}\lesssim 1\) (applied at much shorter times scale as the expected time scale of the Higgs oscillation). In numerical

realization the Dirac function \(\delta (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) \ne n_0\) and \(\Delta (x,t) \ne \Delta

(t)\). In Fig. 1 we present an example time evolution of the absolute value of the order parameter for a selected point in space \(|\Delta (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 \(\delta \Delta _\text {ext}(x,t) = 0\) (no external perturbation) the mode decays spontaneously after some time \(t_{\text {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_{\text {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 work55. 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 papers42,

44, 51, 56,57,58,59. 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_{\textrm{F}}^{-1}\)). We have also checked the stability of the results with respect to the lattice size. For the quasi-1D calculations (see "Lifetime of the Higgs mode and

Properties of the post-decay state" Sections), 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. 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: $$\begin{aligned} V(\Delta ) =

-\frac{1}{2}\mu ^2|{\Delta }|^2 + \frac{1}{4}\epsilon |{\Delta }|^4. \end{aligned}$$ (14) The solutions for a classical particle moving in such potential are well known1, 60. In the case of

small amplitude oscillations around the minimum of the potential \(\Delta _0=\pm \mu /\sqrt{\epsilon }\) we have: $$\begin{aligned} \Delta (t) = \Delta _0 + A(t)\sin (\omega t+\phi ),

\end{aligned}$$ (15) where \(\omega\) is the oscillator’s frequency, \(\phi\) is its phase, and _A_(_t_) is the amplitude, which in general may depend on time13, 26. Moreover, for the small

amplitude Higgs mode, the oscillation frequency is related to the equilibrium value of the pairing gap \(\hbar \omega =2\Delta _0\). The oscillation of arbitrary amplitude are expressed by

delta amplitude Jacobi elliptic function \(\text {dn}\): $$\begin{aligned} |\Delta (t)| = |\Delta (t_0)| \text {dn}(\Omega t+\phi ,k), \end{aligned}$$ (16) where \(\Delta (t_0)=\Delta

_0=0.44e_F\) (for strongly interacting Fermi gas) and the elliptic modulus parameter _k_ is related to maximum and minimum values of the paring field during the oscillations:

$$\begin{aligned} k \approx 1-\left( \frac{\min [\Delta (t)]}{\max [\Delta (t)]}\right) ^2 \end{aligned}$$ (17) In Fig. 2 we compare the numerically extracted time dependence of the pairing

field \(\Delta (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 \times 32 \times

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, \(\omega =2\Delta _0\), has been reproduced with very good accuracy. The large amplitude mode oscillates with lower frequency \(\Omega \approx \Delta _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 as10, 13, 26. 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(\varvec{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: $$\begin{aligned} \sigma [n](t) = \frac{\sqrt{\langle n^2 \rangle (t) - \langle n \rangle ^2(t)}}{\langle n \rangle (t

= 0)}, \quad \text {with} \quad \langle n^k \rangle (t) = \frac{1}{V}\int n^k(\varvec{r},t)\,\textrm{d}\varvec{r}. \end{aligned}$$ (18) The decay time \(t_d\) is defined as the time for

which the normalized standard deviation exceeds a threshold value \(\sigma [n](t>t_d)>\gamma =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 \(\delta \Delta _\text {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 \(\delta \Delta _\text {ext}^{(1)}(x,t)\). We clearly observe a sensitivity of the decay time with respect to

the spatial modulation wavelength \(\lambda \sim 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_{\text {max}}e_{\textrm{F}}\gtrsim 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. 4a. Clearly, there exist an optimal value of

_k_ (resonance) that triggers 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_{\textrm{F}}\approx\) 0.2 to 0.4. In Fig. 4b 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. 4b. Instead of

using a perturbation with well defined wave length \(2\pi /k\), one can instead use a spatially localized perturbation such as \(\delta \Delta ^{(2)}_\text {ext}\sim \delta (x) = \int

\frac{dk}{2\pi }e^{ikx}\). In numerical realization, we model the \(\delta (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 \(|\Delta (k,t)|\) we observe dominant wave-vectors \(\pm k_d\) into which

the mode decays. These are typically in the range of \(0.1\lesssim k_d/k_{\textrm{F}}\lesssim 0.5\). It is comparable to the wave-vector associated with the coherence length scale \(k_{\xi }

= 1/\xi = \Delta /k_F\). Namely, for results presented in Figs. 3, 4 and 5 the wave-vector associated with the coherence length reads \(k_{\xi }/k_{\textrm{F}}\approx 0.22\), which

translates into \(0.5\lesssim k_d/k_{\xi }\lesssim 2.3\). The results obtained for the \(\delta\)-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\lesssim 75\). For later times one may conclude from the graph that the causality is violated. Such a conclusion is incorrect. The observed

effect is a numerical artifact. Namely, we model \(\delta (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 et al.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 \(\sim 1/k_{d}\), i.e. \(\Delta (x,t)=\Delta (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 \(\Delta (x) \sim \cos (k_d x)\). Indeed, in Fig. 5 (right panel) one can

clearly see lines where \(|\Delta (x,t)|=0\) (called nodal lines) for later times. These are lines where the order parameter changes sign. The spatial fluctuations of \(\Delta (x) \sim \cos

(k x)\) are typically regarded as a fingerprint of Larkin-Ovchinnikov (LO) state39, 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. 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\times 64 \times 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 \(\Delta (x,y)\) consists of

regions where the phase changes by \(\pi\), which are separated from each other by nodal lines defined as \(|\Delta (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. 6a,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

spin-imbalanced systems41,42,43. The main difference is that if the imbalance is present (\(N_\uparrow \ne N_\downarrow\)), 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(\varvec{r})=[n_\uparrow

(\varvec{r})-n_\downarrow (\varvec{r})]/[n_\uparrow (\varvec{r})+n_\downarrow (\varvec{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 ferrons42 or ring solitons63. The size of these grains depend on

the global population imbalance of the system \(P=[N_\uparrow -N_\downarrow ]/[N_\uparrow +N_\downarrow ]\). For small imbalances \(P\lesssim 10\%\) their sizes can by of the order \(\sim

10\xi\), and decrease to \(\sim \xi\) as the imbalance increases, and finally converting to well known Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state41. 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. 6c,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 Fig. 6a,c translate into wave-vectors \(k/k_{\textrm{F}}=2\pi /\lambda k_{\textrm{F}}\approx 0.2\) for both \(P=0\%\) and \(P=10\%\) at the

onset of the oscillation. Clearly, these values correspond 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_{\textrm{F}}\sim 10^3\)) the polarized system

evolves towards to the known from past studies configurations41,42,43. Finally, we demonstrate that the observations hold for fully 3D cases, as shown in Fig. 6e. These results were obtained

for a lattice of size \(64 \times 64 \times 64\). Similarly to the quasi-2D case, we introduce the spin-imbalance to the 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 \(\pi\). 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. 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 literature13, 26, 64, 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 (\(\hbar \omega = 2\Delta _0\)) makes it a very valuable

tool from the perspective of experimental measurements. For this reason, recent works (such as54) focus on the stabilization of this mode. Here, and also in37, 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 FFLO phase—a state that is routinely discussed 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 droplets41, 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 work65. 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_{\textrm{F}}\)66—the resolution that should be sufficient to detect predicted here

inhomogeneities in the post-decayed state. DATA AVAILIBILITY The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.

REFERENCES * Pekker, D. & Varma, C. Amplitude/Higgs modes in condensed matter physics. _Annu. Rev. Condens. Matter Phys._ 6, 269–297 (2015). Article  CAS  ADS  Google Scholar  * Endres,

M. _et al._ The ‘Higgs’ amplitude mode at the two-dimensional superfluid/mott insulator transition. _Nature_ 487, 454–458 (2012). Article  CAS  PubMed  ADS  Google Scholar  * Behrle, A. _et

al._ Higgs mode in a strongly interacting fermionic superfluid. _Nat. Phys._ 14, 781–785 (2018). Article  CAS  Google Scholar  * Matsunaga, R. & Shimano, R. Nonequilibrium BCS state

dynamics induced by intense terahertz pulses in a superconducting NbN film. _Phys. Rev. Lett._ 109(18), 187002 (2012). Article  PubMed  ADS  Google Scholar  * Matsunaga, R. _et al._ Higgs

amplitude mode in the BCS superconductors \({{\rm Nb}}_{1-x}{{\rm Ti}}_{x} {{\rm N}}\) induced by terahertz pulse excitation. _Phys. Rev. Lett._ 111(5), 057002 (2013). Article  PubMed  ADS 

Google Scholar  * Matsunaga, R. _et al._ Light-induced collective pseudospin precession resonating with Higgs mode in a superconductor. _Science_ 345(6201), 1145–1149 (2014). Article 

MathSciNet  CAS  PubMed  MATH  ADS  Google Scholar  * Katsumi, K. _et al._ Higgs mode in the \(d\)-wave superconductor \({{\rm bi}}_{2}{{\rm sr}}_{2}{{\rm cacu}}_{2}{{\rm o}}_{8+x}\) driven

by an intense terahertz pulse. _Phys. Rev. Lett._ 120, 117001 (2018). Article  CAS  PubMed  ADS  Google Scholar  * Chu, H. _et al._ Phase-resolved Higgs response in superconducting cuprates.

_Nat. Commun._ 11(1), 1793 (2020). Article  CAS  PubMed  PubMed Central  ADS  Google Scholar  * Shimano, R. & Tsuji, N. Higgs mode in superconductors. _Annu. Rev. Condens. Matter Phys._

11, 103–124 (2020). Article  CAS  ADS  Google Scholar  * Bulgac, A. & Yoon, S. Large amplitude dynamics of the pairing correlations in a unitary Fermi gas. _Phys. Rev. Lett._ 102,

085302 (2009). Article  PubMed  ADS  Google Scholar  * Liu, B., Zhai, H. & Zhang, S. Evolution of the Higgs mode in a fermion superfluid with tunable interactions. _Phys. Rev. A_ 93,

033641 (2016). Article  ADS  Google Scholar  * Barankov, R. A. & Levitov, L. S. Dynamical selection in developing fermionic pairing. _Phys. Rev. A_ 73, 033614 (2006). Article  ADS 

Google Scholar  * Yuzbashyan, E. A., Dzero, M., Gurarie, V. & Foster, M. S. Quantum quench phase diagrams of an \(s\)-wave BCS-BEC condensate. _Phys. Rev. A_ 91, 033628 (2015). Article 

ADS  Google Scholar  * Ojeda Collado, H. P., Usaj, G., Lorenzana, J. & Balseiro, C. A. Nonlinear dynamics of driven superconductors with dissipation. _Phys. Rev. B_ 101, 054502 (2020).

Article  CAS  ADS  Google Scholar  * Tokimoto, J., Tsuchiya, S. & Nikuni, T. Higgs mode in a trapped superfluid Fermi gas. _J. Low Temp. Phys._ 187, 765 (2017). Article  CAS  ADS  Google

Scholar  * Bruun, G. M. Long-lived Higgs mode in a two-dimensional confined Fermi system. _Phys. Rev. A_ 90, 023621 (2014). Article  ADS  Google Scholar  * Tokimoto, J., Tsuchiya, S. &

Nikuni, T. Excitation of Higgs mode in superfluid Fermi gas in BCS-BEC crossover. _J. Phys. Soc. Jpn._ 88, 023601 (2019). Article  ADS  Google Scholar  * Kemper, A. F., Sentef, M. A.,

Moritz, B., Freericks, J. K. & Devereaux, T. P. Direct observation of Higgs mode oscillations in the pump-probe photoemission spectra of electron-phonon mediated superconductors. _Phys.

Rev. B_ 92, 224517 (2015). Article  ADS  Google Scholar  * Schwarz, L. _et al._ Classification and characterization of nonequilibrium Higgs modes in unconventional superconductors. _Nat.

Commun._ 11(1), 287 (2020). Article  CAS  PubMed  PubMed Central  ADS  Google Scholar  * Schwarz, L., Fauseweh, B. & Manske, D. Momentum-resolved analysis of condensate dynamic and Higgs

oscillations in quenched superconductors with time-resolved ARPES. _Phys. Rev. B_ 101, 224510 (2020). Article  CAS  ADS  Google Scholar  * Krull, H., Bittner, N., Uhrig, G., Manske, D.

& Schnyder, A. Coupling of Higgs and Leggett modes in non-equilibrium superconductors. _Nat. Commun._ 7(1), 11921 (2016). Article  CAS  PubMed  PubMed Central  ADS  Google Scholar  *

Nosarzewski, B., Moritz, B., Freericks, J. K., Kemper, A. F. & Devereaux, T. P. Amplitude mode oscillations in pump-probe photoemission spectra from a \(d\)-wave superconductor. _Phys.

Rev. B_ 96, 184518 (2017). Article  ADS  Google Scholar  * Murotani, Y., Tsuji, N. & Aoki, H. Theory of light-induced resonances with collective Higgs and Leggett modes in multiband

superconductors. _Phys. Rev. B_ 95, 104503 (2017). Article  ADS  Google Scholar  * Gao, H., Schlawin, F. & Jaksch, D. Higgs mode stabilization by photoinduced long-range interactions in

a superconductor. _Phys. Rev. B_ 104, L140503 (2021). Article  CAS  ADS  Google Scholar  * Barankov, R. A., Levitov, L. S. & Spivak, B. Z. Collective Rabi oscillations and solitons in a

time-dependent BCS pairing problem. _Phys. Rev. Lett._ 93, 160401 (2004). Article  CAS  PubMed  ADS  Google Scholar  * Barankov, R. A. & Levitov, L. S. Synchronization in the BCS pairing

dynamics as a critical phenomenon. _Phys. Rev. Lett._ 96(23), 230403 (2006). Article  CAS  PubMed  ADS  Google Scholar  * Yuzbashyan, E. A., Tsyplyatyev, O. & Altshuler, B. L.

Relaxation and persistent oscillations of the order parameter in fermionic condensates. _Phys. Rev. Lett._ 96, 097005 (2006). Article  PubMed  ADS  Google Scholar  * Zhou, T.-G. & Zhang,

P. Attenuating dynamics of strongly interacting fermionic superfluids in SYK solvable models. arXiv:2303.02422, (2023). * Scott, R. G., Dalfovo, F., Pitaevskii, L. P. & Stringari, S.

Rapid ramps across the BEC-BCS crossover: A route to measuring the superfluid gap. _Phys. Rev. A_ 86, 053604 (2012). Article  ADS  Google Scholar  * Huang, B., Yang, X., Xu, N., Zhou, J.

& Gong, M. Dynamical instability with respect to finite-momentum pairing in quenched BCS superconducting phases. _Phys. Rev. B_ 99, 014517 (2019). Article  CAS  ADS  Google Scholar  *

Kombe, J., Bernier, J.-S., Köhl, M. & Kollath, C. Finite-duration interaction quench in dilute attractively interacting Fermi gases: Emergence of preformed pairs. _Phys. Rev. A_ 100,

013604 (2019). Article  CAS  ADS  Google Scholar  * Harrison, T. _et al._ Decay and revival of a transient trapped Fermi condensate. _Phys. Rev. Res._ 3, 023205 (2021). Article  CAS  Google

Scholar  * Szymańska, M. H., Simons, B. D. & Burnett, K. Dynamics of the BCS-BEC crossover in a degenerate Fermi gas. _Phys. Rev. Lett._ 94, 170402 (2005). Article  PubMed  ADS  Google

Scholar  * Yi, W. & Duan, L.-M. Dynamic response of an ultracold Fermi gas near the Feshbach resonance. _Phys. Rev. A_ 73, 013609 (2006). Article  ADS  Google Scholar  * Foster, M. S.,

Yuzbashyan, E. A. & Altshuler, B. L. Quantum quench in one dimension: Coherent inhomogeneity amplification and “supersolitons’’. _Phys. Rev. Lett._ 105, 135701 (2010). Article  PubMed 

ADS  Google Scholar  * Ojeda Collado, H. P., Usaj, G., Lorenzana, J. & Balseiro, C. A. Fate of dynamical phases of a BCS superconductor beyond the dissipationless regime. _Phys. Rev. B_

99, 174509 (2019). Article  CAS  ADS  Google Scholar  * Dzero, M., Yuzbashyan, E. A. & Altshuler, B. L. Cooper pair turbulence in atomic Fermi gases. _EPL_ 85, 20004 (2009). Article  ADS

  Google Scholar  * Chern, G.-W. & Barros, K. Nonequilibrium dynamics of superconductivity in the attractive Hubbard model. _Phys. Rev. B_ 99, 035162 (2019). Article  CAS  ADS  Google

Scholar  * Larkin, A. I. & Ovchinnikov, Y. N. Nonuniform state of superconductors. _Zh. Eksperim. i Teor. Fiz._ 47, 9 (1964). Google Scholar  * Fulde, P. & Ferrell, R. A.

Superconductivity in a strong spin-exchange field. _Phys. Rev._ 135, A550–A563 (1964). Article  ADS  Google Scholar  * Tüzemen, B., Zawiślak, T., Wlazłowski, G. & Magierski, P.

Disordered structures in ultracold spin-imbalanced Fermi gas. _New J. Phys._ 25, 033013 (2023). Article  ADS  Google Scholar  * Magierski, P., Tüzemen, B. & Wlazłowski, G. Spin-polarized

droplets in the unitary Fermi gas. _Phys. Rev. A_ 100, 033613 (2019). Article  CAS  ADS  Google Scholar  * Tüzemen, B., Kukliński, P., Magierski, P. & Wlazłowski, G. Properties of

spin-polarized impurities: Ferrons, in the unitary Fermi gas. _Acta Phys. Pol. B_ 51(3), 595 (2020). Article  MathSciNet  ADS  Google Scholar  * Magierski, P., Tüzemen, B. & Wlazłowski,

G. Dynamics of spin-polarized impurity in ultracold Fermi gas. _Phys. Rev. A_ 104, 033304 (2021). Article  CAS  ADS  Google Scholar  * Bulgac, A. Local-density-functional theory for

superfluid fermionic systems: The unitary gas. _Phys. Rev. A_ 76, 040502 (2007). Article  ADS  Google Scholar  * Bulgac, A., Forbes, M. M. & Magierski, P. The unitary Fermi gas: From

Monte Carlo to density functionals. In _The BCS-BEC Crossover and the Unitary Fermi Gas_ (W. Zwerger, ed.), Lecture Notes in Physics, pp. 305–373, Springer (2012). * Bulgac, A. & Yu, Y.

Renormalization of the Hartree–Fock–Bogoliubov equations in the case of a zero range pairing interaction. _Phys. Rev. Lett._ 88, 042504 (2002). Article  PubMed  ADS  Google Scholar  *

Magierski, P., Wlazłowski, G., Makowski, A. & Kobuszewski, K. Spin-polarized vortices with reversed circulation. _Phys. Rev. A_ 106, 033322 (2022). Article  MathSciNet  CAS  ADS  Google

Scholar  * Boulet, A., Wlazłowski, G. & Magierski, P. Local energy density functional for superfluid Fermi gases from effective field theory. _Phys. Rev. A_ 106, 013306 (2022). Article 

CAS  ADS  Google Scholar  * Bulgac, A. & Forbes, M. M. Unitary Fermi supersolid: The Larkin-Ovchinnikov phase. _Phys. Rev. Lett._ 101, 215301 (2008). Article  PubMed  ADS  Google Scholar

  * Wlazłowski, G., Sekizawa, K., Marchwiany, M. & Magierski, P. Suppressed solitonic cascade in spin-imbalanced superfluid Fermi gas. _Phys. Rev. Lett._ 120, 253002 (2018). Article 

PubMed  ADS  Google Scholar  * Bulgac, A., Forbes, M. M., Kelley, M. M., Roche, K. J. & Wlazłowski, G. Quantized superfluid vortex rings in the unitary Fermi gas. _Phys. Rev. Lett._ 112,

025301 (2014). Article  PubMed  ADS  Google Scholar  * W-SLDA Toolkit. https://wslda.fizyka.pw.edu.pl/. * Lyu, G., Xi, K.-T., Yoon, S., Chen, Q. & Watanabe, G. Exciting long-lived Higgs

mode in superfluid Fermi gases with particle removal. _Phys. Rev._ 107, 023321 (2022). Article  Google Scholar  * Bulgac, A., Abdurrahman, I. & Wlazłowski, G. Sensitivity of

time-dependent density functional theory to initial conditions. _Phys. Rev. C_ 105, 044601 (2022). Article  CAS  ADS  Google Scholar  * Wlazłowski, G., Xhani, K., Tylutki, M., Proukakis, N.

P. & Magierski, P. Dissipation mechanisms in fermionic Josephson junction. _Phys. Rev. Lett._ 130, 023003 (2023). Article  PubMed  ADS  Google Scholar  * Barresi, A., Boulet, A.,

Magierski, P. & Wlazłowski, G. Dissipative dynamics of quantum vortices in fermionic superfluid. _Phys. Rev. Lett._ 130, 043001 (2023). Article  CAS  PubMed  ADS  Google Scholar  *

Hossain, K. _et al._ Rotating quantum turbulence in the unitary Fermi gas. _Phys. Rev. A_ 105, 013304 (2022). Article  CAS  ADS  Google Scholar  * Tylutki, M. & Wlazłowski, G. Universal

aspects of vortex reconnections across the BCS-BEC crossover. _Phys. Rev. A_ 103, L051302 (2021). Article  CAS  Google Scholar  * Arefeva, I. Y., Volovich, I. V. & Piskovskiy, E. V.

Rolling in the Higgs model and elliptic functions. _Theor. Math. Phys._ 172(1), 1001–1016 (2012). Article  MathSciNet  MATH  Google Scholar  * Hunter, J. D. Matplotlib: A 2d graphics

environment. _Comput. Sci. Eng._ 9(3), 90–95 (2007). Article  Google Scholar  * Childs, H., Brugger, E., Whitlock, B., Meredith, J., Ahern, S., Pugmire, D., Biagas, K., Miller, M., Harrison,

C., Weber, G. H., Krishnan, H., Fogal, T., Sanderson, A., Garth, C., Bethel, E. W., Camp, D., Rübel, O., Durant, M., Favre, J. M. & Navrátil, P. Visit: An end-user tool for visualizing

and analyzing very large data. In _High performance visualization–enabling extreme-scale scientific insight_, pp. 357–372, (2012). * Barkman, M., Samoilenka, A., Winyard, T. & Babaev, E.

Ring solitons and soliton sacks in imbalanced fermionic systems. _Phys. Rev. Res._ 2, 043282 (2020). Article  CAS  Google Scholar  * Volkov, A. F. & Kogan, S. M. Collisionless

relaxation of the energy gap in superconductors. _Sov. Phys. JETP_ 38, 1018 (1974). ADS  Google Scholar  * Navon, N., Smith, R. P. & Hadzibabic, Z. Quantum gases in optical boxes. _Nat.

Phys._ 17, 1334–1341 (2021). Article  CAS  Google Scholar  * Patel, P. B. _et al._ Universal sound diffusion in a strongly interacting Fermi gas. _Science_ 370, 1222–1226 (2020). Article 

MathSciNet  CAS  PubMed  MATH  ADS  Google Scholar  Download references FUNDING This work was supported by the Polish National Science Center (NCN) under Contracts No.

UMO-2017/26/E/ST3/00428 (GW, ABa, ABo) and UMO-2021/43/B/ST2/01191 (PM). The calculations were executed on Piz Daint supercomputer based in Switzerland at Swiss National Supercomputing

Centre (CSCS), PRACE allocation No. 2021240031. We also acknowledge the Global Scientific Information and Computing Center, Tokyo Institute of Technology for resources at TSUBAME3.0 (project

ID: hp210079). AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Faculty of Physics, Warsaw University of Technology, Ulica Koszykowa 75, 00-662, Warsaw, Poland Andrea Barresi, Antoine Boulet, 

Gabriel Wlazłowski & Piotr Magierski * ISMANS CESI, 44 Avenue Frédéric Auguste Bartholdi, 72000, Le Mans, France Antoine Boulet * Department of Physics, University of Washington,

Seattle, WA, 98195-1560, USA Gabriel Wlazłowski & Piotr Magierski Authors * Andrea Barresi View author publications You can also search for this author inPubMed Google Scholar * Antoine

Boulet View author publications You can also search for this author inPubMed Google Scholar * Gabriel Wlazłowski View author publications You can also search for this author inPubMed Google

Scholar * Piotr Magierski View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS A.Ba., A.Bo. and G.W. contributed to the calculations and data

analysis. All authors contributed to research planning, interpretation of the results and manuscript writing. CORRESPONDING AUTHOR Correspondence to Gabriel Wlazłowski. ETHICS DECLARATIONS

COMPETING INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION PUBLISHER'S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published

maps and institutional affiliations. SUPPLEMENTARY INFORMATION SUPPLEMENTARY INFORMATION 1. SUPPLEMENTARY INFORMATION 2. SUPPLEMENTARY INFORMATION 3. SUPPLEMENTARY INFORMATION 4.

SUPPLEMENTARY INFORMATION 5. RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing,

adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons

licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a

credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted

use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT

THIS ARTICLE CITE THIS ARTICLE Barresi, A., Boulet, A., Wlazłowski, G. _et al._ Generation and decay of Higgs mode in a strongly interacting Fermi gas. _Sci Rep_ 13, 11285 (2023).

https://doi.org/10.1038/s41598-023-38176-9 Download citation * Received: 03 April 2023 * Accepted: 04 July 2023 * Published: 12 July 2023 * DOI: https://doi.org/10.1038/s41598-023-38176-9

SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to

clipboard Provided by the Springer Nature SharedIt content-sharing initiative