Spontaneous symmetry breaking and Nambu-Goldstone modes in open classical and quantum systems

We discuss spontaneous symmetry breaking of open classical and quantum systems. When a continuous symmetry is spontaneously broken in an open system, a gapless excitation mode appears corresponding to the Nambu-Goldstone mode. Unlike isolated systems, the gapless mode is not always a propagation mode, but it is a diffusion one. Using the Ward-Takahashi identity and the effective action formalism, we establish the Nambu-Goldstone theorem in open systems, and derive the low-energy coefficients that determine the dispersion relation of Nambu-Goldstone modes. Using these coefficients, we classify the Nambu-Goldstone modes into four types: type-A propagation, type-A diffusion, type-B propagation, and type-B diffusion modes.


Introduction
Spontaneous symmetry breaking is one of the most important notions in modern physics.When a continuous symmetry is spontaneously broken, a gapless mode appears called the Nambu-Goldstone (NG) mode [1][2][3], which governs the low-energy behavior of the system.For example, a phonon in a crystal, which is the NG mode associated with translational symmetry breaking, determines the behavior of the specific heat at low temperatures, which is nothing but the Debye T 3 law.
The nature of the NG modes such as the dispersion relations and the number of modes is determined by symmetry and its breaking pattern.The relation between them was first shown for relativistic systems [1][2][3], where the number of NG modes is equal to the number of broken symmetries or generators.It was extended to isolated systems without Lorentz symmetry [4][5][6][7][8], where the number of NG modes does not coincide with the number of broken symmetries [9][10][11].A typical example is a magnon in a ferromagnet, which is the NG mode associated with spontaneous breaking of the spin symmetry O(3) to O (2).The number of broken symmetries, dim(O(3)/O(2)), is equal to two, but only one magnon with quadratic dispersion appears.In general, when a global internal symmetry G is spontaneously broken into its subgroup H, the number of NG modes is expressed as [4][5][6][7][8] where N BS = dim(G/H) is the number of broken symmetries, and ρ βα := − [i Qβ , ĵα0 (x)] is the Watanabe-Brauner matrix with the Noether charge Q β and the charge density j α0 (x) of G [12].The NG modes are classified into two types: type-A and type-B modes.The type-B mode is characterized by nonvanishing ρ βα which implies that the two broken charge densities j β0 (x) and j α0 (x) are not canonically independent [13].(rank ρ)/2 counts the number of canonical pairs of broken generators, and every pair of broken charges constitutes one NG mode.Therefore, the number of type-B modes, N B , is equal to (rank ρ)/2.The rest of the degrees of freedom, N A = N BS − rank ρ, become the type-A modes, which have the same property as NG modes in Lorentz-invariant systems.The sum of the number of type-A and type-B modes leads to Eq. ( 1).Both type-A and type-B NG modes are propagation modes with typically linear and quadratic dispersions, respectively.Similar to isolated systems, spontaneous symmetry breaking occurs in open systems.A well-known example is a diffusion mode in the synchronization transition of coupled oscillators, which describe chemical and biological oscillation phenomena [14,15].In the synchronization transition, U (1) phase symmetry is spontaneously broken, and the diffusion mode appears as the NG mode, which has different dispersion from that in isolated systems.Another example is ultracold atoms in an optical cavity [16].In the system, a laser and its coupling to radiation fields give rise to spontaneous emission and dissipation, where the internal energy does not conserve.Bose-Einstein condensation and the symmetry breaking can occur even in such a case, and they have been observed [16].Furthermore, the synchronization transition and the diffusive NG mode of ultracold atoms in the driven-dissipative setup are discussed in Refs.[17][18][19][20][21][22].The diffusive NG modes are characteristic of open systems.
Compared with isolated systems, the relation between the NG modes, the broken symmetries, and the dispersion relations has not been established in open systems.A crucial difference between isolated and open systems is the lack of ordinary conserved quantities such as the energy, momentum, and particle number, because of the interaction with the environment.Thus, we cannot naively apply the argument for isolated systems to that for open systems.In our previous work, we studied properties of the NG modes in open systems based on toy models [22], in which we found two types of NG modes: diffusion and propagation modes, whose poles have ω = −iγ|k|2 and (±a − ib)|k| 2 , respectively.Here, γ, a, and b are constants.In the model study, the nonvanishing Watanabe-Brauner matrix for the open system leads to the propagation NG modes, where a similar relation to Eq. ( 1) is satisfied [22].However, it is found that the propagation modes split into two diffusion modes by changing a model parameter in the study of time-translation breaking [23], where Eq. ( 1) is not always satisfied.Therefore, we need a model-independent analysis to understand the nature of NG modes for open classical and quantum systems.This is the purpose of the present paper.
Open classical and quantum systems can be uniformly described by the path integral formulation, called the Martin-Siggia-Rose-Janssen-De Dominicis (MSRJD) formalism [24][25][26][27] for classical systems, and the Keldysh formalism for quantum systems [28].Both path integral formulations are written in two fields.In the language of classical theories, they are called classical and response fields.In the Keldysh formalism, they correspond to the a linear combination of fields on the forward and backward paths.These doubled fields play an important role in the symmetry of open systems.
The notion of symmetry in open systems is slightly different from that in isolated systems [21][22][23].In isolated systems, when there is a continuous symmetry, there exists a physical Noether charge, which we call Q α R in this paper."Physical" means that the Noether charge is an observable like the energy and momentum.In contrast, in open systems, the Noether theorem does not necessarily lead to the physical conserved charge.Instead, another conserved charge, which we refer to as Q α A , arises in the path integral formulation.The doubled charges, Q α R and Q α A , relate to the fact that the path integral is written in doubled fields.Although Q α A itself is not a physical conserved quantity, it plays the role of the symmetry generator.By using Q α A , we can define the spontaneous symmetry breaking for open systems and derive the Ward-Takahashi identities [29,30].As is the case in isolated systems, spontaneous symmetry breaking implies the existence of gapless excitation modes.Our previous analysis based on the Ward-Takahashi identities [22] is limited in the zero-momentum limit.In this paper, we generalize it to analysis with finite momentum, and derive the low-energy coefficients for the inverse of the Green functions in the NG mode channel.Using these coefficients, we classify the NG modes, and discuss the relation between these modes and the broken generators.
The paper is organized as follows: In Sec. 2, we show our main result, in which we classify the NG modes into four types and discuss their dispersion relations.We also discuss how the NG modes can be observed in experiments.In Sec. 3, we review the path integral formulation in open classical and quantum systems.In Sec. 4, we discuss the concept of symmetry in open systems and provide two field-theoretical technique that we employ to show our main result.Section 5 shows the detailed derivation of our main result.Section 6 is devoted to the summary and discussion.
Throughout this paper, we use the relativistic notation in (d + 1)-dimensional spacetime with the Minkowski metric η µν = diag(1, −1, −1, • • • , −1), although the system need not have the Lorentz symmetry.We employ the natural units, i.e., c = = 1, where c is the speed of light and is the Planck constant over 2π.

Main result
In this section we summarize our main results first, because the derivation is technical and complicated.We show the relationship between the inverse of the retarded Green function and low-energy coefficients that determine the dispersion relation of NG modes.Using the low-energy coefficients, we classify the NG modes into four types: type-A propagation, type-A diffusion, type-B propagation, and type-B diffusion modes.We also discuss how these modes can be observed in experiments.

Green functions and low-energy coefficients
We consider a (d + 1)-dimensional open system that has a continuous internal symmetry G (the precise definition of symmetry in open systems is shown in Sec.4.1).Suppose that the symmetry is spontaneously broken into its subgroup H in a steady state.We assume that the steady state is unique and stable against small perturbations.The steady state may be not only the global thermal equilibrium but also a nonequilibrium steady state.We also assume that any spacetime symmetries of the open system such as the time and spatial translational symmetry are not spontaneously broken, which implies that the frequency and momentum are good quantum numbers.We are interested in the behavior of the retarded Green functions that contain NG modes [G π (k)] αβ , where α and β run from 1 to the number of broken symmetries N BS .As in the case of isolated systems, N BS is equal to dim(G/H).Our main result (generalization of the Nambu-Goldstone theorem) is that the inverse of the retarded Green function can be expanded as with coefficients: Here, µ and ν are spacetime indices, k µ = (ω, k) are the frequency and the wave vector, and d d+1 x = dtdx.The subscript c denotes the connected part of correlation functions.Q π is the projection operator that removes the contribution of NG modes from the operator [the definition of Q π (= 1 − P π ) is given in Eq. ( 110)].The expectation values are evaluated in the path integral formulation, Eq. ( 26).In this formulation, there are two types elementary of fields denoted as χ a R and χ a A , which correspond to the classical and fluctuation fields, respectively.The operators in Eqs. ( 4) and ( 6) are related to the symmetry transformation of the action S[χ a R , χ a A ]. To see this, consider an infinitesimal local transformation of fields under G: , where α (x) is the infinitesimal parameter depending on spacetime, and δ α A χ a i (x) is the infinitesimal transformation of χ a i (x).If χ a i (x) is a linear representation of G, it can be represented as , where T α is the generator of G. Since G is the symmetry, the action is invariant under δ A if α is constant.For a spacetime-dependent α (x), the action transforms as where In isolated systems, these are also symmetry transformations; however, they are not in open systems [21,22].Under this transformation, the action transforms as Since δ R transformation is not the symmetry of the action, h α R (x) exists.

4/29
δ β R j αµ A (x) and S βα;νµ (x) in Eqs. ( 4) and ( 6) are given through the infinitesimal local transformation of The low-energy coefficients in the inverse of retarded Green functions in Eqs. ( 3)-( 6) are expressed as one-or two-point functions of these operators.The ordinary Nambu-Goldstone theorem shows C βα = 0 in (3), which is derived from the only symmetry breaking pattern, and thus it is independent of the details of the underlying theory.This claims that there is at least one zero mode when a continuum symmetry is spontaneously broken.To determine the dispersion relation and the number of NG modes, we need additional data.If we impose Lorentz invariance on the system, we find C βα;µ = 0 and C βα;νµ = −η νµ g βα , where g βα is an N BS × N BS matrix with det g = 0.In this case, det G −1 π = 0 has 2N BS solutions with ω = ±|k|.Each pair of solutions ω = ±|k| gives one mode.Therefore, the number of NG modes is equal to N BS , whose dispersion is linear.This is the Nambu-Goldstone theorem in relativistic systems [3,31].Equations ( 4) and (6) give the data for general cases, not just for isolated systems without Lorentz invariance but also for open ones.In this sense, our formulae in Eqs. ( 2)-( 6) provide the generalization of the Nambu-Goldstone theorem.In the next subsection, we discuss the classification of NG modes and their dispersion relations.

Classification of NG modes
It is interesting to clarify the relation between the broken symmetries, the NG modes, and their dispersion relations in open systems.In isolated systems without Lorentz invariance, those relations have been made clear in Refs.[4][5][6][7].In this section, we generalize the relation obtained in isolated systems to that in open systems.The formula in isolated systems can be reproduced as a special case of our result.
Our formulae in Eqs. ( 2)-( 6) are quite general, so we need additional assumptions to classify NG modes.First, we assume that det(−iωρ − ω 2 ḡ) = 0 at an arbitrary small but nonzero ω, where we define ρ βα := C βα;0 and ḡβα := −C βα;00 .This assumption means that the low-energy degrees of freedom are contained in, at least, up to quadratic order in ω.We note that this is implicitly assumed in the analysis of NG modes in isolated systems [4].Second, we assume that the action of the underlying theory satisfies the reality condition All the models discussed in Sec. 3 satisfy this condition.The reality condition implies that all the coefficients C βα;µ••• in Eq. ( 2) are real.This property leads to the relation that if which makes our analysis simpler.In this section, we concentrate on systems that satisfy this condition.One can generalize our analysis by relaxing the this assumption.From these three assumptions, we can conclude that there are two types of modes: diffusion and propagation, and the poles of the retarded Green functions have the form: Here, a, b, and γ are real and positive due to the stability of the system and they vanish at k = 0 due to symmetry breaking.We note that the above three assumptions do not exclude the 5/29 possibility of the pole ω = 0 with finite momentum k or equivalently the pole with a = 0, b = 0, or γ = 0.The pole can be excluded from the assumption that the translational symmetry is not broken discussed in Sec.2.1.If such a pole exists, there exists a time-independent local operator.By operating the operator to a given steady state, we can generate another steady state with a finite momentum, which breaks the translational symmetry.This contradicts the assumption of broken translational invariance 1 .
Let us now focus on the relation between the number of modes and the coefficients ρ βα and ḡβα .The number of zero modes can be evaluated from the number of solutions of det(−iωρ − ω2 ḡ) = ω NBS det(−iρ − ωḡ) = 0 with ω = 0, which is equal to N BS + (N BS − rank ρ) = 2N BS − rank ρ.As in the case of isolated systems [4], we characterize the type-B modes by linearly independent row vectors in ρ. ρ generally has both symmetry and antisymmetric parts.This is different from isolated systems, where ρ is an antisymmetry matrix.The symmetric part plays the role of dissipation.The number of type-B degrees of freedom is equal to rank ρ.This characterization is different from our previous work [22], in which we proposed that the type-B mode was characterized by the antisymmetry part of ρ, and conjectured that the existence of the antisymmetric part implied the existence of a propagation mode.However, this conjecture is not always true [23]; it depends on the details of the underlying theory, as will be seen below.Therefore, we change the definition of type-B modes 2 .
Since propagation poles are always paired with the positive and negative real parts, we count the pair of poles as two.On the other hand, the pole of a diffusion mode is counted as one.This observation gives the relation between rank ρ and the number of diffusion and propagation modes as As mentioned above, the number of diffusion modes depends on the details of the theory.
To see this, let us consider a simple model with and ḡβα = 0. Here, κ 1 and κ 2 are parameters.We choose these to be positive.Since det ρ = κ 1 κ 2 + 1 = 0, rank ρ = 2, i.e., there are two solutions in det G −1 π = 0, which are given as If 4 > (κ 1 − κ 2 ) 2 , one propagation mode appears.On the other hand, if 4 < (κ 1 − κ 2 ) 2 , two diffusion modes appear.In both cases, Eq. ( 12) is, of course, satisfied; however, the diffusion or propagation depends on the parameters.The remaining of gapless degrees of freedom describe type-A modes, whose number is (N BS − rank ρ).Since they have ω 2 terms in G −1 π , we count each degree of freedom as two, and we find that the following relation is satisfied: One might think the existence of the type-A diffusion mode is unnatural.However, we cannot exclude this possibility at this stage.For example, let us consider G , which correspond to the type-A modes because ρ = 0.If ζ < 1, there is one propagation mode.On the other hand, if ζ > 1, there are two diffusion modes.These are nothing but the type-A diffusion modes in our classification.The possibility of type-A diffusion modes is excluded by assuming an additional condition, as will be seen later.
Combining Eqs. ( 12) and ( 15), we find the general relation: For isolated systems, where h β R vanishes, the transformation generated by δ β R is upgraded to symmetry whose Noether current is j βµ R .In this case, the symmetric part of ρ βα vanishes, and ρ βα turns into the Watanabe-Brauner matrix: . In isolated systems, there are no diffusion modes, i.e., N A-diffusion = N B-diffusion = 0, so the relation reduces to the standard one: Furthermore, let us consider that the system is isotropic.In this case, C αβ;ij can be expressed as g αβ δ ij .We assume that det g = 0, which simplifies the dispersion relations: Type-B diffusion mode In this case, there is no type-A diffusion mode (see Appendix A for the proof).
In addition to gapless solutions, gapped solutions in det G −1 π (ω, k) = 0 may exist.As is in the case of the gapless modes in Eq. ( 11), the gapped solutions can be classified into damping and gapped modes, with at k = 0.The total number of solutions of det G −1 π (ω, 0) = 0 is equal to the degree of the polynomial det , the degree is equal to (N BS + rank ḡ).There are 2N BS − rank ρ solutions with ω = 0, so that the remaining (rank ḡ + N BS ) − (2N BS − rank ρ) = rank ḡ + rank ρ − N BS solutions correspond to gapped ones.Therefore, we find the relation for the gapped modes to be rank ḡ + rank ρ − N BS = 2N gapped + N damping . ( If γ G or |a G − ib G | is much smaller than the typical energy scale, these gapped modes become the low-energy degrees of freedom.In isolated systems, these kinds of gapped modes are known as gapped partners [5,7,[32][33][34][35][36].An example is a Ferrimagnet, in which ferroand antiferro-order parameters coexist.The Ferrimagnet has one magnon with quadratic dispersion, which is the type-B mode, and one gapped mode. 7/29

Spectral function and experimental detection
The dispersion relation of NG modes can be experimentally observed by inelastic scattering processes [37][38][39].For example, in atomic Fermi superfluids, the NG mode is observed in the spectra with focused Bragg scattering [39].The differential cross section is proportional to the correlation function, which behaves as 2Im G π (ω, k)/ω =: S(ω, k) at small ω [40]3 .S(ω, k) will be multiplied by an additional factor depending on the processes.From the results in the previous subsection, the retarded Green functions in broken phases are written as The corresponding spectra respectively.As a comparison, we show the functional form of S = ρ(ω, k)/ω for type-A and -B modes in isolated systems as respectively.Figures 1 and 2 illustrate the ω-|k| dependence of S(ω, k) in open and isolated T q n g 8 z T j + B J C 9 g O j y p m 9 < / l a t e x i t > Fig. 2 The left and right panels show the ω-|k| dependence of S(ω, k) for type-A and B modes in an isolated system, respectively.The both have the sharp pair peaks.The parameters are chosen as a A = 1, b A = 0.5, a B = 1, and b B = 0.5.isolated system has sharp pair peaks.Thus, S(ω, k) has significantly different behaviors depending on the open or the isolated system.

Open classical and quantum systems
Here, we briefly review the path integral approach to open classical and quantum systems.Readers who are familiar with this approach may jump to Sec. 4. As will be seen below, in both classical and quantum systems, the expectation value of an operator can be expressed as the path integral Here, the subscript ρ denotes the contribution from the initial density operator, whose definition is given later.S[χ a R , χ a A , C, C] is the action with two types of physical degrees of freedom, χ a R and χ a A , corresponding to classical and fluctuation fields, respectively.In classical systems, χ a A is often called the response field.In quantum systems, χ a R and χ a A are fundamental fields of the Keldysh basis in the Keldysh path formalism [28].C and C are ghost fields, which are responsible for the conservation of probability.The existence of the ghost term depends on the theory.In the following, we show three examples for classical and quantum systems that can be expressed as Eq. ( 26).

Stochastic system
The first example is a stochastic system.The path integral formalism in this system is the MSRJD formalism [24][25][26][27].Let us consider a stochastic system whose dynamics is described by the Langevin type equation: where ξ represents Gaussian white noise that satisfies Here, κ represents the strength of the noise.The probability distribution is defined as where ρ[φ R (t I )] is the initial probability distribution.Since φ(t) follows Eq. ( 27), P [t F ; φ R ] can be expressed as where det(∂ t + γ + 3λφ 2 R ) is the Jacobian associated with the transformation from Here, we introduced the following notation: Using the Fourier representation of the delta function, we can write the probability distribution as The average of the noise can be evaluated as which leads to The determinant can be expressed as the Fermionic path integral: We eventually obtain the path-integral formula: with the action, The expectation value of an operator O[φ R , φ A ] is given as In the second line, we include the integral of φ R (t F ) at the boundary into ρ Dφ R .Here, we only considered a single variable with simple interaction and noise.Generalization to multi-component fields is straightforward.For more detailed derivation, see Refs.[43][44][45].

e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " y d c U 7 c L J X B 2 d x I 1 T 5 b p I j U u 9 D h g = " >
A e e / i M D / i K n 3 / W 6 j g 1 + l 7 a N M s D L T P K / r P l 9 M e / q g b N A e e / i M D / i K n 3 / W 6 j g 1 + l 7 a N M s D L T P K / r P l 9 M e / q g b N A e e / i M D / i K n 3 / W 6 j g 1 + l 7 a N M s D L T P K / r P l 9 M e / q g b N A e e / i M D / i K n 3 / W 6 j g 1 + l 7 a N M s D L T P K / r P l 9 M e / q g b N x S L u 6 3 y w g q s w j r 1 Y x t i k I A k Z J 0 X P Y c r u B Y C w p a w J + w P U g W P q 1 m C H y E k v g B i l Z M o < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x S L u 6 3 y w g q s w j r 1 Y x t i k I A k Z J 0 X P Y c r u B Y C w p a w J + w P U g W P q 1 m C H y E k v g B i l Z M o < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x S L u 6 3 y w g q s w j r 1 Y x t i k I A k Z J 0 X P Y c r u B Y C w p a w J + w P U g W P q 1 m C H y E k v g B i l Z M o < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Fig. 3 Closed time contour in the Keldysh formalism.

Open quantum systems
The second example is an open quantum system.The open system can be formally constructed from an isolated system.To see this, let us consider a quantum system coupled with an environment.The path integral formula can be obtained by integrating the environment out [46].The action of the total system consists of three parts: where S sys [φ], S env [B], and S int [B, φ] are the actions of the system, the environment, and the interaction between them, respectively.φ and B are the degrees of freedom of the system and the environment.We assume that the initial density operator at t = t I is the direct product of those of the system and the environment: ρ = ρsys ⊗ ρenv .In the path integral formalism, the expectation value of an operator O[φ 1 , φ 2 ] is given on the path shown in Fig. 3 as where the indices 1 and 2 represent the label of the forward and backward paths in Fig. 3, and we introduced the matrix elements of the density operators ρ sys [φ we can express the expectation value as where we defined and the notation of the subscript ρ by To see the connection between S eff and the action in Eq. ( 37) in the MSRJD formalism, let us move on to the Keldysh basis, which is defined as φ R := (φ 1 + φ 2 )/2 and φ A := (φ 1 − φ 2 ).φ R and φ A are called classical and quantum fields, respectively.Expanding the action with 11/29 respect to φ A , and keeping φ A up to quadratic order, we can obtain the MSRJD action.For example, we consider φ 4 theory in (0 + 1) dimension, whose action is given as In the Keldysh basis, S sys [φ 1 ] − S sys [φ 2 ] is expressed as The functional form of Γ[φ 1 , φ 2 ] depends on the details of the environment.For a simple case, we assume where ν and D are parameters coming from couplings to the environment.Then, the action reduces to For a classical treatment, we drop φ 3 A term, and neglect the ∂ t φ A ∂ t φ R term for slow dynamics.Then, we arrive at where we rescaled φ A → −φ A /ν and introduced γ = m 2 /ν, κ = D/ν 2 and λ = u/ν.We can see that the action in Eq. ( 49) is the same form as the MSR one in Eq. ( 37) except for the Jacobian term.

Lindblad equation
The third example is a system described by a Lindblad equation [47] (see Ref. [21] for a review of the path integral formalism of the Lindblad equation).The Lindblad equation with a single Lindblad operator L is given as Here, γ is a coefficient and L is a function of fields.The term proportional to γ describes the dissipation and fluctuation effects.By construction, the trace of ρ is time independent, ∂ t tr ρ = 0.When the degrees of freedom is a bosonic Schrödinger field ψ(t, x) in (d + 1) dimensions with the commutation relation [ ψ(t, x), ψ † (t, x )] = δ (d) (x − x ), the path integral for the expectation value of an operator O[ψ i , ψ † i ] is given as with the action, where H i := ψ i | Ĥ|ψ i , and For example, if we choose the Hamiltonian and the Lindblad operator as the action is written in the form, where m is the mass.In the Keldysh basis, it becomes As we have seen in the three examples, all the systems are described by the path integral.

Field-theoretical technique
In this section, we provide the concept of symmetry in open systems and two field-theoretical techniques to derive low-energy coefficients.One is the Ward-Takahashi identity that gives the relation between different Green functions.The other is the generating functional and effective action methods, which are often employed to prove the Nambu-Goldstone theorem in quantum field theories (see, e.g., Ref. [31]).

Symmetry of open systems
In isolated systems, the existence of a continuous symmetry implies the existence of conserved current.On the other hand, in open systems the energy, momentum, and particle number are generally not conserved due to the interaction between the system and the environment.Thus, one may wonder what symmetry of open systems means.Even in such a case, symmetry may exist [21,22].To see the symmetry of open systems, let us consider the concrete example with the action in Eq. (54).If γ = 0, this system corresponds to an isolated system.In this case, the action is invariant under ψ 1 → e iθ1 ψ 1 and ψ 2 → e iθ2 ψ 2 , where θ 1 and θ 2 are constant.In this sense, there is the U (1) 1 × U (1) 2 symmetry, whose Noether charges are respectively.Here, we defined Q 2 such that [ Q2 , ψ2 ] = + ψ2 in the operator formalism, whose sign is opposite to [ Q1 , ψ1 ] = − ψ1 .This relative sign is caused by the canonical commutation relation, [ ψ1 (t, x), ψ † 1 (t, x )] = δ (d) (x − x ) and [ ψ2 (t, x), ψ † 2 (t, x )] = −δ (d) (x − x ) in the operator formalism.This can be understood as the sign in front of the time derivative term in the Keldysh action d d+1 (ψ † 1 i∂ t ψ 1 − ψ † 2 i∂ t ψ 2 ).This kind of doubled symmetry is used for the construction of an effective field theory of fluids in the Keldysh formalism [48][49][50][51][52][53][54][55][56][57][58][59].
In the Keldysh basis, these are expressed as In isolated systems, these charges generate an infinitesimal transformation as where δ i ψj = −i[ Qi , ψj ] in the operator formalism.On the other hand, when γ = 0, one of the symmetries is explicitly broken.The residual symmetry is U (1) A , which is given by setting the parameters as θ 1 = θ 2 .The Noether charge Q A is still conserved, but Q R is not, which means that the charge in the usual sense is not conserved.More generally, we can consider a system with continuous symmetry G.The symmetry of open systems means that the action is invariant under the infinitesimal transformation of G, where χ a A (x) and χ a R (x) are field degrees of freedom, and α is constant.We can formally define the transformation δ α R χ i (x) such that These are generalizations of Eqs. ( 59) and ( 60) [22].In open systems, invariance of the action under δ α R depends on the theory.If the particle number is conserved, but energy is not, then, the action is invariant under the U (1) R transformation χ A → e iθ χ R and χ R → e iθ χ A /4, but it is not invariant under the time translation χ A (t) → χ R (t + c) and χ R (t) → χ A (t + c)/4, where θ and c are constants.
In addition, all the previous examples in Sec. 3 satisfy the reality condition: In this paper, we focus on theories satisfying the reality condition.

Ward-Takahashi identity
Symmetry plays an important role not only in isolated systems but also in open ones.Here, we show the Noether current in open systems and the Ward-Takahashi identity [29,30].Let us consider a theory with field variables χ a A and χ a R .Suppose the action S is invariant under G and the infinitesimal transformation is given as χ a i → χ a i + α δ α A χ a i , with an infinitesimal 14/29 constant α .When α (x) depends on the spacetime coordinate, the action transforms as because δ A S vanishes if α (x) is constant.The operator j αµ A (x) is called the Noether current.To see the Ward-Takahashi identity in the path integral formalism, we consider the expectation value of an operator O[χ a i ], which is given as Since χ a i is the integral variable, the integral is invariant under the relabeling χ a i → χ a i : If we choose χ a i (x) = χ a i (x) + α (x)δ α A χ a i (x), and assuming that the path-integral measure is invariant under this transformation: Dχ a i = Dχ a i , we find that Here, the local transformation of O[χ a i ] is given as The leading-order term in α (x) gives Differentiating Eq. ( 67) with respect to α (x), we obtain This is the Ward-Takahashi identity, which is the conservation law in the path integral formalism.

Generating functional and effective action
In order to show the Nambu-Goldstone theorem, it is useful to introduce the generating functional.We start with a path integral representation for the generating functional in (d + 1) spacetime dimensions: where the χ a i are elementary degrees of freedom in the Keldysh basis.φ a i = (φ a R , φ a A ) is a set of real scalar fields, which may be the elementary or a composite field of χ a i 5 .In general, the 5 When φ a R and φ a A are composite, these fields are defined by using 1/2 basis such that φ R = (φ 1 + φ 2 )/2 and φ A = φ 1 − φ 2 , where φ 1 and φ 2 are polynomials of χ 1 and χ 2 , respectively.15/29 stationary state need not be the thermal equilibrium state, i.e., a nonequilibrium steady state is allowed in this formalism.For the generating functional to be well defined, we assume that the stationary state is stable against any small perturbations.We also assume a stationary state that is independent of the choice of the initial density operator, so that we omit the subscript ρ in the following.
Connected Green functions are generated by differentiating ln Z[J] with respect to J: Here, ... J denotes the expectation value in the presence of the source J i a (x), which is defined as The subscript c denotes the connected part of the correlation function.The Green function without the external field is given as the limit of J i a (x) → 0: Since the J A a φ a A term in Eq. ( 69) is proportional to φ a A , it can be absorbed into the action as the external force, while J R a φ a R cannot.The external force does not change the conservation of the probability, so that the partition function becomes trivial if J R = 0: Z[J R = 0, J A ] = 1.This identity leads to Therefore, any correlation functions contracted from only the φ A vanish.In particular, φ a A (x) cannot be the order parameter.
We introduce here the effective action, which is defined as where ϕ a i (x) := φ a i (x) J = −iδ ln Z[J]/δJ i a (x).We assume that ϕ a i (x) = φ a i (x) J is invertible so that the effective action is well defined.The effective action is a functional of ϕ a i (x) rather than J i a (x), since δΓ In the absence of the external field, J i a (x) = 0, Eq. (75) gives the stationary condition of the effective action.The second functional derivative of Γ is Here, we defined the two-point Green function as To see the symmetry of the effective action, we apply the Ward-Takahashi identity in Eq. (67).By choosing O = exp i d d+1 J i a (x)φ a i (x) in Eq. ( 67), we obtain or equivalently, where we used Eq. ( 75).These equations play an essential role in the analysis of Nambu-Goldstone modes.

Spontaneous symmetry breaking and the Nambu-Goldstone theorem
In this section, we nonperturbatively establish the Nambu-Goldstone theorem in open systems.We will derive the formulae for the inverse of the retarded Green function for NG modes shown in Eqs. ( 2)-( 6).We consider an open system with elementary fields χ a i (x).The system has spacetime translational invariance, and a continuous internal symmetry G.The fields transform under an infinitesimal transformation of G: , where the T α are generators of G.We assume that the order parameter belongs to a set of real fields {φ a i } that transforms under G as In other words, the φ a i belong to a linear representation of G. φ a i may be elementary or composite.The representation of φ a i may be different from the χ a i .For technical reasons, we assume that φ a i transforms under the local transformation: This assumption will make the derivation of our formulae simple.We assume that the continuous symmetry G is spontaneously broken to its subgroup H.The symmetry breaking is characterized by the nonvanishing order parameter, δ α A φa i := δ α A φ a i (x) .We also assume that the translational symmetry is not broken, i.e., the order parameter is independent of time and space, which enables us to work in momentum space.
Our procedure consists of three steps: First, we show the existence of gapless excitations.Second, we will find the relation between the inverse of the retarded Green functions for {φ a i } and their low-energy coefficients.In general, {φ a i } contains not only NG fields but also fields with gapped modes that we are not interested in.Therefore, in the third step, we derive the inverse of the retarded Green functions with only NG fields by projecting out the contribution from gapped modes.

Existence of gapless excitations
The existence of gapless excitations can be shown by using the standard technique developed in quantum field theory [31].In the generating functional method, the order parameter is 17/29 given as δ α A φa i := lim J→0 δ α A φ a i (x) J .Taking α (x) to be constant in Eq. ( 79), we obtain Differentiating Eq. ( 82) with respect to ϕ b j (x ), we get where we used δ α A ϕ a i (x) = i[T α ] a b ϕ b i (x) and Eq.(76).Taking the limit J j a → 0, we obtain where In momentum space, we obtain This identity represents the eigenvalue equation with the zero eigenvalue, whose eigenvectors are δ α A φa R , and it implies the existence of gapless excitations.One can check that the independent number of δ α A φa R equals dim(G/H) =: N BS .In an isolated system with Lorentz symmetry, N BS is equal to the number of gapless excitations [31].However, this is not the case in open systems.To obtain information on finite frequency and momentum, we need the further data shown in the following subsections.

Low-energy coefficients
In order to obtain finite momentum information, let us go back to Eq. ( 79) and consider its functional derivative with respect to ϕ b j (x ): Here, we employed Eq. ( 7) on the left-hand side and defined the local transformation, to simplify the notation.This local transformation is obeyed from Eq. (81).We can also introduce the local transformation δ R ϕ a i (x) defined as Here, we introduced the symbol ε j i for the transformation defined from δ α R χ a i (x) = ε j i δ α A χ a j (x) in Eq. ( 8).The bar on the ¯ α is introduced to distinguish the transformation for δ A , where α is used.Multiplying Eq. (86) by − d d+1 x δ R ϕ b j (x ) and taking the J i a (x) → 0 limit, we obtain where we interchange the left-and right-hand sides.The left-hand side represents the inverse of the retarded Green function in the NG mode channel.We would like to express the righthand side in the language of symmetry.Naively, one might think that the right-hand side in 18/29 Eq. ( 89) has the form − δ R δ A S .This is not the case.In general, a transformation of fields and the expectation value are not commutative due to fluctuations, i.e., δ R O J = δ R O J , where we define δ R O := d d+1 x ¯ α (x) δO δχ a i (x) δ R O J := d d+1 x δ R ϕ a i (x) Recall that χ a i (x) is the elementary field of theory and all local operators are polynomials of χ a i .The equality holds if O is a linear function of φ a R and φ a A .To see the explicit relation between δ R O J and δ R O J , let us consider the Ward-Takahashi identity for δ R O, which leads to We note that since the action is not invariant under this transformation, i δ R S O J does not vanish even when ¯ β (x) is constant.For the trivial operator O = 1, Eq. ( 92) reduces to i δ R S J + i d d+1 x J i a (x)δ R ϕ a i (x) = 0.
In order to obtain a more compact expression, we introduce the projection operator defined as Here, we assumed that O does not explicitly depend on ϕ.We also introduce Q φ := 1 − P φ , which satisfies P 2 φ = P φ , Q 2 φ = Q φ , and P φ Q φ = Q φ P φ = 0.By construction Q φ ∆φ a i (x) = 19/29 0, that is, the projection operator Q φ removes the linear component of ∆φ a i (x) from the operator.The first two terms on the right-hand side of Eq. ( 96) are simply expressed as Noting that [D −1 ] ji ba (x , x; ϕ) = δJ j b (x )/δϕ a i (x), and using the chain rule, we obtain the last term of Eq. ( 96 (99) Eventually, we arrive at the expression, If we choose O = δ A S and take the J → 0 limit, Eq. ( 94) becomes lim Substituting Eq. (101) into Eq.( 89), we obtain This expression gives the relation between the inverse of the retarded Green function and expectation values of operators.The indices a and b in [D −1 ] AR ba (x − x) include not only NG fields but also other fields with gapped modes.Therefore, our next step is to find the inverse of the retarded Green function consisting of only the NG fields.which are consistent with Eq. (137) 9 . 9The index of α and β represents the index of broken generators, which is different from the index of fields.α = 1, 2, and 3 correspond to χ 2 , χ 1 and ψ 2 , respectively.26/29

Summary and discussion
We have derived low-energy coefficients of the inverse of retarded Green functions for Nambu-Goldstone modes in open systems, Eqs. ( 2)-( 6), which provide a generalization of the Nambu-Goldstone theorem.As is the case in isolated systems, we have classified the NG modes into type-A and type-B modes by using the coefficient of the single time-derivative term, ρ βα .These modes are further classified into diffusion and propagation modes.The relation between broken symmetries and the numbers of these modes are summarized in Eqs. ( 12), (15), and (16).In this paper, we have employed continuum models.It is straightforward to generalize our results to systems with discrete translational symmetry, where the momentum is replaced by the Bloch momentum.
We have not taken into account hydrodynamic modes, which are expected to appear when Q α R is conserved.The coupling between NG and hydrodynamic modes may change low-energy behaviors.The limit k → 0 in Eqs. ( 123) and (124) to obtain Eqs. ( 2)-( 6) might not be well defined, although our formulae in Eqs. ( 118)-(120) still give the correct results.In such a case, it is better to treat the hydrodynamic modes as the dynamical degrees of freedom.
It is also fascinating to consider a spontaneous breaking of spacetime symmetries.In particular, time-translational symmetry breaking will be interesting in the context of "quantum time crystal" [62].In isolated systems, it is shown that the time-translational symmetry cannot be broken in the ground state [63].On the other hand, in open systems, time-translationalsymmetry breaking can occur.For example, in reaction-diffusion systems, time-translational symmetry breaking is known as synchronization phenomena, and the dynamics of phase, which corresponds to that of the NG mode, is described by a nonlinear diffusion equation [14].A similar situation can occur in open quantum systems [23,64].
In the clarification of NG modes, we assumed G −1 π (ω, k) = G −1 π (ω, −k), which leads to C βα;i vanishing.It will be interesting to relax this condition.If there is a steady-state current that is an order parameter, C βα;i will be nonvanishing.In this case, a different type of NG mode, which cannot be classified in this paper, will appear.
Another direction is to consider spontaneous symmetry breaking of higher-form symmetries [65].Photons can be understood as the NG modes associated with the spontaneous breaking of U (1) one-form symmetry [65].As in the ordinary symmetry breaking, there also exists the type-B photon, whose dispersion is quadratic [66][67][68].It will be interesting to see how such a mode changes in an open system.We leave these issues for future work.

Fig. 1
Fig. 1 The left panel shows the ω-|k| dependence of S(ω, k) for a type-B diffusion mode, which has a single peak.The right panel shows S(ω, k) for a type-B propagation mode, which has blunt pair peaks.The parameters are chosen as γ B = 1, a B = 1, and b B = 0.5. k e s N k z p y 5 5 8 6 Z u b p j S s 8 n e u p Q O r u 6 e 3 r 7 + m M D g 0 P D 8 c T I 6 K 5 e s N k z p y 5 5 8 6 Z u b p j S s 8 n e u p Q O r u 6 e 3 r 7 + m M D g 0 P D 8 c T I 6 K 5 e s N k z p y 5 5 8 6 Z u b p j S s 8 n e u p Q O r u 6 e 3 r 7 + m M D g 0 P D 8 c T I 6 K 5 e s N k z p y 5 5 8 6 Z u b p j S s 8 n e u p Q O r u 6 e 3 r 7 + m M D g 0 P D 8 c T I 6 K 5 Oφ b j (x) J = δ O J /δiJ j b (x) + φ b j (x) J O J , and from Eq. (93), we can expressδ R O J as δ R O J = −i δ R S O J − i d d+1 x J i a (x )¯ α (x )i[T α ] a b ε j i δ δiJ j b (x ) O J + φ b j (x ) J O J = −i δ R S O c;J − i d d+1 x J i a (x )¯ α (x )i[T α ] denotes the connected diagram defined as O 1 O 2 c;J := ∆O 1 ∆O 2 J with ∆O := O − O J .We would like to eliminate the explicit dependence of J j b (x) in Eq. (94).For this purpose, consider the functional derivative of Eq. (93) with respect to ϕ b j (x ), which becomes δ δϕ b j (x )i δ R S J + i d d+1 x [D −1 ] ji ba (x , x; ϕ)δ R ϕ a i (x) + iJ i a (x )¯ α (x )i[T α ] a b ε j i = 0. (95)Substituting this into Eq.(94), we can expressδ R O J as δ R O J = −i δ R S O c;J + d d+1 x d+1 x d d+1 x [D −1 ] ji ba (x , x; ϕ)δ R ϕ a i (x) d+1 x d d+1 x ∆φ b j (x )[D −1 ] ji ba (x , x; ϕ) φ a i (x)O c;J .