Collective Excitations and Nonequilibrium Phase Transition in Dissipative Fermionic Superfluids

We predict a new mechanism to induce collective excitations and a nonequilibrium phase transition of fermionic superfluids via a sudden switch-on of two-body loss, for which we extend the BCS theory to fully incorporate a change in particle number. We find that a sudden switch-on of dissipation induces an amplitude oscillation of the superfluid order parameter accompanied by a chirped phase rotation as a consequence of particle loss. We demonstrate that when dissipation is introduced to one of the two superfluids coupled via a Josephson junction, it gives rise to a nonequilibrium dynamical phase transition characterized by the vanishing dc Josephson current. The dissipation-induced collective modes and nonequilibrium phase transition can be realized with ultracold fermionic atoms subject to inelastic collisions.

In this Letter, we theoretically investigate collective excitations and a nonequilibrium phase transition of fermionic superfluids driven by a sudden switch-on of two-particle loss due to inelastic collisions between atoms. By formulating a dissipative BCS theory that fully incorporates a change in particle number, we find that dissipation fundamentally alters the superfluid order parameter and induces collective oscillations in its amplitude and phase. In particular, we elucidate that a coupling between the order parameter and dissipation leads to a chirped phase rotation, in sharp contrast to the case of an interaction quench in closed systems [see Fig. 1(a)].
To experimentally observe the collective phenomena induced by dissipation, we propose introducing a particle loss in one of two coupled superfluids to induce a relativephase oscillation analogous to the Leggett mode [2, [33][34][35][36][37][38] [see Fig. 1(b)]. The phase mode causes an oscillation of a Josephson current around a nonvanishing dc component. Remarkably, when dissipation becomes strong, the coupled system undergoes a nonequilibrium phase transition characterized by the vanishing dc Josephson current, which can be regarded as a generalization of a dynamical phase transition [11,12,72,73] to dissipative quantum systems. Our findings can experimentally be tested with ultracold atoms through introduction of dissipation via a photoassociation process [50,59].
Dissipative BCS theory.-We consider ultracold fermionic atoms described by the three-dimensional attractive Hubbard model (a) (b)  1. (a) Schematic illustration of the amplitude and phase modes in a Mexican-hat free-energy potential as a function of the complex order parameter ∆, when either the interaction UR or the dissipation γ is suddenly switched on. A sudden quench of the interaction UR and that of the dissipation γ kick ∆ in a direction parallel and perpendicular to the radial direction, respectively. Note that a finite change of γ excites both the phase and amplitude modes. (b) Two superfluids coupled via a Josephson junction, where one superfluid (system 2) is subject to two-body loss.
σ fermion with momentum k (at site i). When the system is subject to inelastic collisions, scattered atoms are lost to a surrounding environment, resulting in dissipative dynamics as observed experimentally [50,55,56,58].
Here, we study the time evolution of the density matrix ρ which is described by the Lindblad equation [39,40] where L i = c i↓ c i↑ is a Lindblad operator that describes two-body loss with loss rate γ > 0. We note that the kinetic energy of lost atoms is large because of large internal energy of atoms before inelastic collisions. Under such situations, atoms after inelastic collisions are quickly lost into the surrounding environment and the Born-Markov approximation is justified [74][75][76]. We first study how the standard BCS theory is generalized in open dissipative systems by formulating a timedependent mean-field theory in terms of a closed-timecontour path integral [77,78]. We start with a generating functional defined as with an action where the subscripts + and − denote forward and backward paths, H α = kσ kckσα c kσα − U R ic i↑αci↓α c i↓α c i↑α , L iα = c i↓α c i↑α , andL iα = c i↑αci↓α (α = +, −). Note that the action has U(1) symmetry under c iσα → e iθ c iσα though the particle number is not conserved [79,80]. By introducing auxiliary fields via the Hubbard-Stratonovich transformation, we rewrite the action in a quadratic form of fermionic Grassmann fields as [49,81] Here ∆ is the superfluid order parameter which can be determined from the requirement that the action be extremal as [81] where U = U R + iγ/2 is an effective complex coupling constant including a contribution from the atom loss [49], and N 0 is the number of sites. Importantly, the order parameter includes the loss rate γ, which leads to dissipation-induced collective modes as discussed below. The action (5) describes the mean-field time-evolution equation of the density matrix as where Ψ k = (c k↑ , c † −k↓ ) t is the Nambu spinor. In the Supplemental Material [81], we show that Eq. (7) can be derived from two different methods, i.e. the mean-field theory for the Lindblad equation and the time-dependent Bogoliubov-de Gennes analysis. While Eq. (7) appears to describe unitary evolution, it is consistent with the original Lindblad equation (2) as a consequence of the time-dependent BCS ansatz [81].
We use Anderson's pseudospin representation [1, 7-12, 14, 32] defined by σ k = 1 2 Ψ † k · τ · Ψ k and H eff = 2 k b k · σ k , where τ = (τ x , τ y , τ z ) is the vector of the Pauli matrices. The pseudospins satisfy the commutation relations [σ j k , σ k k ] = i jkl σ l k . For simplicity of notation, we omit the bracket and regard σ k as the expectation value of the pseudospin operator. By using the commutation relation of the pseudospins, Eq. (7) is mapped to the Bloch equation: Equation (9) shows that the superfluid dynamics is characterized by precession of a pseudospin in an effective magnetic field b k . Here, the order parameter is determined self-consistently from the pseudospin expectation value as It is noteworthy that the norm of the pseudospin is conserved by the Bloch equation (9). The time evolution of the particle number due to particle loss is obtained from Eq. (7) as which reflects the dynamics of the order parameter.
As an initial state, we prepare a BCS ground state with γ = 0, whose pseudospin representation is given by The single-particle energy k is measured from the Fermi energy of the initial state. The bandwidth W is defined by the energy difference between the upper and lower edges of the energy spectrum with a constant density of states. We then switch on the atom loss γ at t = 0. The results shown in Fig. 2 are obtained by the second-order Runge-Kutta method. In the long-time limit, the amplitude of the superfluid order parameter ∆ is suppressed due to dissipation, indicating a decay of superfluidity [see Fig. 2(a)]. We note that the order parameter decays in the long-time limit due to a decrease of the particle number [see Fig. 2(b)], and such behavior has no counterpart in the quench in isolated systems [10,11]. Remarkably, after the dissipation γ is introduced, the U(1) phase of the order parameter rotates and shows chirping, i.e., its angular velocity increases with time [see Fig. 2(a), (b)] as a consequence of the dynamical shift of the Fermi level [81]. This property is unique to the dissipative superfluid and distinct from the usual dynamics in isolated systems where the U(1) phase stays constant [10-12]. The phase rotation is understood from an initial-state free energy as a function of ∆ [see Fig. 1(a)]. When dissipation is introduced, the sudden quench of the imaginary part of U in Eq. (11) pushes the order parameter towards the direction perpendicular to the radial direction irrespective of the initial choice of the gauge. Another way to understand the phase rotation is to introduce an effective chemical potential as ∆(t) = exp(−2i t 0 µ eff (t)dt)Ω(t) (Ω ∈ R). By performing a global gauge transformation from ∆ to Ω, the Bloch equation is written in the Larmor frame on which the energy dispersion is given by ξ k (t) = k − µ eff (t). This gauge transformation indicates that the phase rotation corresponds to a decrease of the effective chemical potential, which is consistent with the behaviors ofθ and N in Fig. 2(b). This result can naturally be understood since the phase and the particle number are conjugate variables.
We also find amplitude oscillations in |∆| as shown in Fig. 2(a). The amplitude oscillations are more pronounced when the interaction and the dissipation are simultaneously quenched [81]. The mechanism behind the oscillations is that the quench of the imaginary part of U changes the absolute value of ∆ [see Fig. 1(a)]. The frequency of the amplitude oscillation is close to 2∆ 0 at an early stage, and increases as time evolves. This behavior is distinct from that of an isolated system, where the amplitude mode is characterized by the constant frequency. Such behavior can be observed from the measurement of the time-dependent particle number via Eq. (12).
Collective excitations: Leggett mode.-To observe the chirped phase rotation of the superfluid order parameter that is a unique feature of dissipative superfluids, we propose that the phase rotation induced by dissipation can be detected when two superfluids are connected via a Josephson junction, which has been realized in ultracold atoms [82][83][84][85][86]. As the phase difference in the two superfluid order parameters is gauge-invariant, it leads to an observable Josephson current. We introduce dissipation to one of the two superfluids as schematically illustrated in Fig. 1(b) and assume that they are coupled via a tunneling Hamiltonian [2, 38] where V > 0 is the amplitude of Cooper-pair tunneling between system 1 without dissipation and system 2 with two-particle loss. By performing a meanfield analysis, we can write the system Hamiltonian as In the pseudospin respresentation, the Hamiltonian is written as H i = 2 k b ik · σ ik with an effective magnetic field b ik = (Re∆ i , −Im∆ i , ik ), which yields the Bloch equation dσ ik /dt = 2b ik × σ ik . The self-consistent conditions for the order parameters read where N 0 is the number of sites of each system. Here, the relations  satisfied. Then, the Josephson current between the two superfluids is given by the rate of change in the particle number of system 1: where δ = tan −1 (−γ/2U R ) is the phase shift due to the sudden switch-on of the atom loss. We numerically solve the coupled Bloch equations for σ ik . We assume that dissipation γ and tunneling V are turned on at t = 0 for the BCS ground state. The numerical results for weak dissipation are shown in Fig. 3(a1)-(d1). In Figs. 3(a1) and (b1), the dynamics of two superfluids almost synchronize with each other because the time scale of particle loss is comparable with the inverse tunneling rate. In the pseudospin picture, the dynamics of particle numbers shown in Fig. 3(c1) can be interpreted as the nutation of pseudospins. Importantly, we see that, although the particle number of the system decreases in time, the corresponding amplitude of the order parameter stays almost constant. This implies that the condensate fraction against the total particle number becomes larger than that of the initial state. As inferred from Fig. 3(d1), the Josephson current oscillates around its dc component. Such behavior is reminiscent of Shapiro steps in a Josephson junction under irradiation of a microwave [87]; however, in the present case, the Josephson current oscillate spontaneously without any external field. Moreover, from Fig. 3(d1), the frequency of the oscillation of the phase difference between the two systems is close to that of the relative-phase mode known as the Leggett mode [2, 38] whose dispersion relation is given by ω L = 2 (λ 12 + λ 21 )|∆ 1 ||∆ 2 |/detλ, where λ 11 = λ 22 = U R /W , λ 12 = λ 21 = V /W and detλ = λ 11 λ 22 − λ 12 λ 21 . We note that ω L includes the effect of loss through the order parameters. The Leggett mode with frequency ω L has been discussed in the context of a collective mode in a multiband superconductor irradiated by light [38]. The agreement between the frequencies of the relative-phase modes in very different situations can be understood as follows. When dissipation is weak, the time evolution of an order parameter is given by with an effective chemical potentials µ ieff . Then, by performing a global gauge transformation from c ikσ to c ikσ exp(i t 0 i µ ieff dt/2), we can linearize the Bloch equation with respect to the relative phase difference between ∆ i 's by following Ref. [38].
Nonequilibrium phase transition.-In the presence of strong dissipation, the order parameter of system 2 oscillates faster than that of system 1 [see Fig. 3(a2), (b2)] and the phase difference monotonically increases in time [see Fig. 3(d2)]. This is because the dissipation rate larger than the tunneling rate makes system 1 fail to follow the decay of system 2, resulting in the dynamics similar to that of a single superfluid shown in Fig. 2. In particular, the chirped phase rotation of the superfluid order parameter of system 2 can be detected from the Josephson current [ Fig. 3(d2)]. As the superfluidity of system 2 is suppressed, the Josephson current also decays, and the particle number of system 1 settles to a constant after some transient time [see Fig. 3(c2)]. The latter behavior is attributed to the continuous quantum Zeno effect [49,74,75,[88][89][90], which states that strong dissipation prevents tunneling and inhibits loss in system 1. In fact, an effective decay rate of system 1 is given by γ eff ≡ |V eff | 2 /γ with an effective tunneling rate Phase difference between the two systems (blue) and particle numbers of system 1 (red) and system 2 (yellow) after a sufficiently long time (t f = 97.9/∆0). The parameters used are UR = 3.06∆0, V = 0.02∆0, and W = 5.11∆0.
The two dynamically distinct regimes of superfluid behaviors suggest the existence of dynamical phases of matter [72,73] in dissipative superfluids.
The qualitative change in the superfluid behaviors with respect to the dissipation strength highlights a dynamical phase transition characterized by the vanishing dc Josephson current [ Fig. 4(a)], where the dc component of the Josephson oscillation is defined by 14)] after a sufficiently long time evolution with t f = 97.9/∆ 0 . We emphasize that the dynamical phase transition in dissipative superfluids is essentially distinct from the phase transition between ground states in a non-Hermitian BCS superfluid [49]. The former is caused by a change in particle number in the long-time dynamics, whereas the latter is caused by an exceptional point of a non-Hermitian BCS Hamiltonian, which is relevant to the short-time dynamics during which the number of particles does not change [91]. From Fig. 4(b), we see that the phase difference θ 2 − θ 1 starts to increase monotonically at the critical point and that the difference in particle number (N 2 − N 1 )/N 0 becomes much larger. The behavior of the phase difference is reminiscent of the localization-diffusion transition of a quantum-mechanical particle moving in a washboard potential in the presence of frictional force [92][93][94]. However, the origin of the transition shown in Fig. 4 is essentially different from frictional force, since it cannot change the particle number. In fact, as shown in the Supplemental Material [81], the dynamical phase transition in Fig. 4 is triggered by the competition between the Josephson coupling and particle loss. Moreover, as the steady state is a vacuum due to the particle loss, the dynamical phase transition is observed only in the transient dynamics, and thus distinct from steady-state transitions.
Conclusions.-We have investigated the loss-quench dynamics of fermionic superfluids, and have demonstrated that the dynamics exhibits amplitude and phase modes with chirped oscillations, the latter of which is a salient feature of a dissipative superfluid. To observe the chirped phase rotation, we have proposed a Josephson junction comprised of dissipative and nondissipative superfluids. We have shown that the relative-phase Leggett mode can be detected from the Josephson current for weak dissipation. Remarkably, when dissipation becomes strong, the superfluids exhibit the unique nonequilibrium phase transition triggered by particle loss. Our prediction can be tested with ultracold atomic systems of 6 Li [83,84], for example, by introducing dissipation using photoassociation processes [50,59]. It is of interest to explore how the dimensionality or confinement by a trap potential affects the dynamics and associated collective modes [15][16][17][18].
We are grateful to Yuto Ashida, Philipp Werner, Shuntaro Sumita, and Yoshiro Takahashi for fruitful discussions. This work was supported by KAKENHI

Supplemental Material for "Collective Excitations and Nonequilibrium Phase Transition in Dissipative Fermionic
Superfluids"

Detailed calculations of the time-dependent dissipative BCS theory
We explain the details of the Hubbard-Stratonovich transformation that is used for the derivation of the timedependent dissipative BCS theory in the path-integral formalism. The action (4) in the main text is rewritten as We perform the Hubbard-Stratonovich transformation for each term in the second line of Eq. (S1) with auxiliary fields ∆ α (α = +, −, ±) as which yield whereψ kα = c k↑α , c −k↓α t and ψ kα = c k↑α ,c −k↓α t (α = +, −). Then, from the saddle-point condition ∂S/∂∆ α = ∂S/∂∆ α = 0 (α = +, −, ±), we obtain where N 0 is the number of lattice sites and · · · is the expectation value for fixed ∆ α and∆ α . Then, we can reduce the number of the auxiliary fields by using c −k↓α c k↑α = tr(c −k↓ c k↑ ρ) (α = +, −) and tr(A † ρ) = [tr(Aρ)] * [78], giving Finally, the action (S5) is rewritten as where the superfluid order parameter of the system is given by Operator formalism of the BCS theory for a dissipative superfluid Here, we explain the operator formalism of the BCS theory for a dissipative superfluid that leads to Eq. (7) in the main text. First, we note that an operator |ψ + ψ − | acting on the Hilbert space of the system can be mapped to a tensor-product state |ψ + ⊗ |ψ − in the doubled Hilbert space H + ⊗ H − [95,96]. Using this mapping, we can rewrite the Liouvillian [see Eq.
(2) in the main text for definition] as where L iα = c i↓α c i↑α , and c iσα (c kσα ) with α = +, − is the fermion annihilation operator in the real (momentum) space that acts on the Hilbert space H α . The fermion operator with subscript + (−) corresponds to the fermion field in the forward (backward) path in the path-integral formalism. The BCS Hamiltonian equivalent to Eq. (1) is given by and H α is defined as By applying the mean-field approximation to L, we obtain the mean-field Liouvillian as where Ψ k = c k↑ , c † −k↓ t is the Nambu spinor. As we can see from Eq. (S13) and Eq. (S15), the Liouvillian is invariant under the U(1) gauge transformations c kσ+ → e iθ c kσ+ and c kσ− → e iθ c kσ− . Moreover, under the exchange of forward and backward operators, P c kσ+ P −1 = c kσ− and P c † kσ+ P −1 = c † kσ− with P 2 = 1, the Liouvillian has the following symmetry (S17) By imposing the same symmetry on the mean-field Liouvillian as P (iL MF ) † P −1 = −iL MF , we obtain the relations for the order parameters as which coincide with those obtained in the path-integral formalism [see Eqs. (S9) and (S10)]. Finally, by rewriting the superfluid order parameter ∆ + as we obtain the equations for the density matrix asρ which are the same as Eqs. (7) and (8) in the main text.
Generalization of the Bogoliubov-de Gennes analysis with a time-dependent BCS state We explain that Eq. (7) (Eq. (9) in the psedouspin representation) in the main text is equivalent to the Bogoliubovde Gennes equation with a time-dependent BCS state [3,9], which describes the unitary evolution of the density matrix.
We introduce the time-dependent BCS state of the effective Hamiltonian H eff = k Ψ † k k ∆ ∆ * − k Ψ k as follows: where |0 is the vacuum of fermions. Here, the superfluid order parameter ∆ is rewritten as Suppose that the density matrix is given by Then, the time-evolution equation [Eq. (7) in the main text] is equivalent to the Bogoliubov-de Gennes equation with the time-dependent BCS state By defining f k (t) and g k (t) as Eq. (S28) is rewritten as These equations take the same forms as those for closed systems [3,9]. Finally, by defining the psedouspins as we obtain the same Bloch equation as discussed in the main text: We note that the dynamics described by Eq. (S27) conserves the purity tr[ρ 2 ] as In general, the purity should decrease during the time evolution described by the quantum master equation [Eq.
(2) in the main text]. This fact is consistent with the time-dependent BCS ansatz (S23) as follows. Since |Ψ BCS (t) can be expanded in terms of N -particle states the density matrix is written as Then, for a gauge-invariant observable O, its expectation value is given by where is a mixed state of different particle numbers. Therefore, concerning gauge-invariant observables, the time-dependent BCS state (S26) is indistinguishable from the mixed state (S42) with tr[ρ 2 ] < 1. Since any physically observable quantity should be gauge invariant, the time-dependent BCS ansatz (S23) can describe the time evolution of the density matrix consistently with the quantum master equation (2).

Dynamics of pseudospins after switch-on of dissipation
The physical origin of the chirping of the U(1) phase can be understood from the pseudospin picture. As shown in Fig. S1, when the sign of σ z k changes from positive to negative, the magnitudes of σ x k and σ y k increase due to the norm conservation of pseudospins. This indicates that the Cooper-pair amplitude at specific momenta rapidly changes when atoms at those momenta are lost from the system. Since Cooper pairs are formed near the Fermi surface, a loss of Cooper pairs leads to a downward shift of the Fermi level.
To see the effect of the dynamics of pseudospins on the collective phase mode, we calculate the angular velocity of the order parameter. From Eq. (11) in the main text, the real and imaginary parts of the order parameter are written as By differentiating Eqs. (S43) and (S44) with respect to time, we obtain Then, by substituting the Bloch equation (S35) into Eq. (S45), we arrive at in which the last term increases due to the shift of the Fermi level, leading to chirping of the U(1) phase. Here, we note that the first term on the right-hand side of Eq. (S46) also increases as the particle number decreases; however, it is much smaller than the second term.
FIG. S1. Dynamics of pseudospins [σ x k (light green), σ y k (blue), σ z k (orange), |σ k | (violet)] after the atom loss with γ = 2.81∆0 is switched on for the initial state with UR = 12.2∆0 and the Fermi energy F = 0 for = −18.7∆0 (left), −9.35∆0 (center), and 0 (right), where is the single-particle energy in the band −23.4∆0 ≤ ≤ 23.4∆0. Figure S2 shows the dynamics after the dissipation γ is introduced at t = 0 and the interaction strength is simultaneously changed from U R = 8.4∆ 0 to U R = 16.8∆ 0 . In Figs. S2(a) and (b), we see an amplitude oscillation larger than that in the loss quench dynamics shown in Figs. 2(a) and (b) in the main text. This behavior is due to the fact that a change in the real part of U causes a large initial shift of the amplitude of the order parameter [see Fig. 1(a) in the main text]. The U(1) phase rotates with an increasing angular velocity due to dissipation as in Figs. S2(a) and (b). The amplitude oscillation of the order parameter can be detected through monitoring of the particle number and from Eq. (12) in the main text. As shown in Figs. S2(c) and (d), the amplitude of the particle-number oscillation is a few percent of the initial particle number, which can be detected with current experimental techniques [83]. Since the oscillation in the particle number cannot appear in isolated systems, the dissipation-induced dynamics can be used as a unique signature for the amplitude mode of the superfluid order parameters.

A simplified model for understanding the nonequilibrium phase transition
Here, we explain how the nonequilibrium phase transition associated with the vanishing dc Josephson current can be understood by considering a simplified model of the Josephson junction with particle loss. We consider two fermionic superfuids coupled via a Josephson junction, where two-body loss is introduced to one of them, as shown in Fig. 1(b) in the main text. The Josephson current flowing between the two systems is given by where, from Eq. (14) in the main text, ∆θ = θ 2 − θ 1 , I 0 = 4V |∆ 1 ||∆ 2 |/U R |U |, and we have neglected the phase shift δ because δ ∆θ. Taking into account the particle loss [see Eq. (12) in the main text], we obtain the rate of change in the particle number of system 1 and that of system 2 as We assume that the time evolution of the phase difference ∆θ is given by the effective chemical-potential difference ∆µ eff between the two systems as where W is a bandwidth and we assume a constant density of states for simplicity (we can also understand this equation from the phenomenological time-dependent Ginzburg-Landau theory, which is explained in the last part of this section). We obtain the equation of motion for ∆θ by using Eqs. (S48), (S49), and (S50) as This system is regarded as a Josephson junction with shunt resistance R = +∞, capacitance C = 1/2W and an external force F = γ|∆ 2 | 2 /|U | 2 , which is described as That is, the time evolution of ∆θ is equivalent to that of a particle moving in a washboard potential The condition for the extremum of V wash is given by dV wash /d∆θ = 0, giving The solution to this equation does not exist for γ|∆ 2 | 2 /I 0 |U | 2 > 1 and the time evolution of ∆θ becomes unstable. If we assume |∆ 1 | |∆ 2 | when the time evolution is sufficiently slow, we obtain the critical strength of the atom loss as which is of the same order of magnitude as that in the main text (γ c 3V ). Thus, the system exhibits a dynamical phase transition from the state in which ∆θ oscillates around an extremum of V wash for γ < γ c to the state in which ∆θ slips down the washboard potential for γ > γ c . Thus, the dynamical phase transition caused by the particle loss is the one between a trapped state and a running state. We note that the particle loss γ acts as an external force F rather than friction R in Eq. (S52). The loss-induced dynamical phase transition occurs spontaneously without any external fields, and has an essentially different origin from the localization-delocalization transition induced by friction R [92][93][94]. These features are also obtained from a phenomenological introduction of two-body loss [74,75], under which the rate of change in the particle number of system 1 and that of system 2 are given by where κ 2 is the two-body loss rate. By using Eqs. (S50), (S56), and (S57), we obtain the equation of motion for ∆θ as This system is regarded as a Josephson junction with resistance R = +∞, capacitance C = 1/2W and an external force F = κ2 2 N2 N0 2 .
As the system has the washboard potential the condition for the extremal of V wash is given by If we assume that the variation of the parameters is sufficiently slow and approximate them as constant, the critical strength of the loss rate where the solution of ∆θ becomes unstable is given by We can numerically solve Eq. (S56), (S57), and (S58), and the results are shown in Fig. S3. We see that the results shown in Fig  0 dtµ ieff (t))|∆ i (t)|), which reflects the conjugate nature of the particle number and the phase. From a more phenomenological point of view, Eq. (S47) can be understood using the time-dependent Ginzburg-Landau (GL) theory. Here we note that, strictly speaking, the GL theory cannot be applied to zero-temperature superfluids considered in our study, since the Taylor expansion of the free energy requires that the system should be close to the transition temperature. In nonequilibrium situations, the time-dependent GL theory cannot be applied to gapped superfluids and nonadiabatic regimes, since the time scale of the order-parameter dynamics becomes shorter than the lifetime of quasiparticles and quasiparticle contributions cannot be neglected [8,32]. Moreover, for the application of the time-dependent GL theory, it is required that the deviation from equilibrium is sufficiently small, which cannot be satisfied in our dissipative superfluids which are driven far from equilibrium. However, we assume below that the GL theory could be phenomenologically used to investigate the dynamics of the dissipative superfluids. We start from the phenomenological GL equation of the superfluid condensate, which is expanded in a series of the order parameter where Γ is a positive constant, ϕ is the scalar potential, e is the charge, and we have assumed that the order parameter is spatially uniform. Here, we have introduced complex-valued coefficients a and b, which might describe the effect of the atom loss. However, we remark on the validity of this treatment below. In our study, the scalar potential does not exist. Instead, we introduce an effective chemical potential µ eff that is determined from the total particle number of the system as follows: −Γ ∂∆ ∂t + 2iµ eff ∆ = a∆ + b|∆| 2 ∆.
By rewriting the equation in terms of the phase θ and the amplitude |∆| of the order parameter (the latter is proportional to the square root of the superfluid density), we obtain the equation of motion for the phase as ∂θ ∂t = −2µ eff − Im(a) + Im(b)|∆| 2 Γ .
For two superfluids connected via a Josephson junction, we similarly obtain where the subscript 1 and 2 label the two superfluids, and we have introduced the effect of loss to superfluid 2. If the coefficients a and b are real, we arrive at Eq. (S50): We note, however, that Eqs. (S48) and (S49) cannot be obtained from the GL theory, since N 1 and N 2 are the total particle numbers of the two superfluids, which are not directly related to the amplitude of the order parameter for dissipative superfluids. In the case of complex coefficients a and b, one might at first sight think that the effect of loss can be described by regarding the second term in Eq. (S65) as an effective change in the chemical potential, which reads ∆µ W (N − N 0 )/N 0 ∝ |∆| 2 − |∆ 0 | 2 . However, this result of the phenomenological treatment of the time-dependent GL theory has a serious problem. To see this, let us differentiate this equation with respect to time as dN/dt ∝ |∆|d|∆|/dt. As the amplitude of the order parameter oscillates as shown in Fig. 2(a) in the main text, the equation dN/dt ∝ |∆|d|∆|/dt indicates that the total particle number of the system can increase since the oscillating part d|∆|/dt can take a positive value. Such an increase of the particle number is unphysical since our system has no particle gain. As the correct equation obtained from the microscopic theory is dN/dt ∝ −|∆| 2 [Eq. (12) in the main text], which indicates that the total particle number of the system monotonically decreases, Eq. (S65) does not correctly describe the dynamics associated with the loss in the total particle number. Thus, it is highly nontrivial to consistently describe the dynamic evolution of the order parameter, collective excitations, and the dynamical phase transition caused by the particle loss in dissipative fermionic superfluids on the basis of the phenomenological GL theory.
If we wish to derive the time-dependent GL equation from the formalism based on the closed-time-contour path integral, we should start from a generating functional like Eq. (3) in the main text, and introduce the order parameter by performing the Hubbard-Stratonovich transformation. Then, by integrating out the fermionic degrees of freedom, we arrive at the action with respect to the order parameter. By expanding the exact action S(∆, ∆ * ) around its extremal value S(∆ 0 , ∆ * 0 ), where ∆ 0 is usually the equilibrium value that satisfies the BCS gap equation, the timedependent GL equation is obtained from a saddle-point equation ∂S(∆, ∆ * )/∂∆ * (r, t) = 0 [97]. However, we have to pay careful attention to the condition that is assumed in the standard derivation: the deviation of the order parameter ∆ around its extremal value ∆ 0 should be small enough for the Taylor expansion to be justified. We also note that the standard derivation assumes that the total particle number of the system does not change. In contrast, as the particle number significantly changes in our system, the phase of the order parameter rotates and largely deviates from the stationary value as shown in Fig. 2 in the main text. Thus, a time-dependent GL theory that fully incorporates the change in particle number needs highly nontrivial consideration, which deserves further study.