Asymptotic correlation functions and FFLO signature for the one-dimensional attractive Hubbard model

We study the long-distance asymptotic behavior of various correlation functions for the one-dimensional (1D) attractive Hubbard model in a partially polarized phase through the Bethe ansatz and conformal field theory approaches. We particularly find the oscillating behavior of these correlation functions with spatial power-law decay, of which the pair (spin) correlation function oscillates with a frequency $\Delta k_F$ ($2\Delta k_F$). Here $\Delta k_F=\pi(n_\uparrow-n_\downarrow)$ is the mismatch in the Fermi surfaces of spin-up and spin-down particles. Consequently, the pair correlation function in momentum space has peaks at the mismatch $k=\Delta k_F$, which has been observed in recent numerical work on this model. These singular peaks in momentum space together with the spatial oscillation suggest an analog of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state in the 1D Hubbard model. The parameter $\beta$ representing the lattice effect becomes prominent in critical exponents which determine the power-law decay of all correlation functions. We point out that the backscattering of unpaired fermions and bound pairs within their own Fermi points gives a microscopic origin of the FFLO pairing in 1D.


Introduction
With addition of an on-site interaction term to the tight binding hamiltonian, the Hubbard model successfully provides a paradigm for condensed matter physics [1,2,3,4]. In contrast to its simple form, this model exhibits diverse features of many-body systems, such as a Mott phase, high T c superconductivity, quantum phase transition, Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase, spin-charge separation etc. Recently this model has received further attention in the experimental developments of trapping ultracold atoms on optical lattices, which provides a prominent opportunity to verify theoretical predictions of this kind [5,6,7,8,9,10]. The one-dimensional (1D) Hubbard model can be exactly solved through the application of the Bethe ansatz [11,12], a thorough study of which may also help to understand key aspects of many-body physics in higher dimensions.
The Hubbard model with attractive interaction is considered as a promising candidate to explain high T c superconductivity. To this end, understanding the pairing mechanism in 1D is of significant importance. According to the Bardeen-Cooper-Schrieffer (BCS) theory, a Cooper pair is formed by electrons with opposite spins and momenta and total momentum zero. This balance between Fermi energies breaks down in the presence of strong magnetic field, so that a novel superconductive state -the FFLO state -appears [13,14]. The FFLO pair carries non-zero centre-of-mass momentum, originating from the strong magnetic field. The superconducting order parameter and density of spins in the FFLO state exhibit a periodic oscillation in the spatial coordinate. The experimental observation of the FFLO state in various materials has been sought for decades. Within this scenario, more evidence has been found in heavy-fermion systems [15,16].
The pair mechanism of the 1D attractive Hubbard model has been investigated [17,18,19]. Cooper pairs in 1D exist in the region B < B c1 below the critical magnetic field B c1 , where the average distance between pairs is much larger than the average pair size. This means that the single-particle Green's function decays exponentially and the single-pair correlation function decays as a power of distance. Once exceeding the critical magnetic field (B > B c1 ), the field begins to break up the pairs, and both the above correlation functions decay as a power of distance.
The FFLO state in the 1D attractive Hubbard model has been extensively studied by various methods, such as the numerical approaches of the density matrix renormalization group (DMRG) [20,21,22,23,24] and quantum Monte Carlo (QMC) [25,26]. It has been found that the pair correlation complies with a power-law decay, i.e., n pair ∝ cos (k F F LO |x|) /|x| α , and its corresponding momentum distribution has peaks in the position of k F F LO = π (n ↑ − n ↓ ) [20]. As far as we know, the confirmation of the FFLO state in 1D mainly relies on numerics, although the application of conformal field theory (CFT) has provided a definite answer for the FFLO signature of the 1D attractive Fermi gas [27]. In this paper we focus on asymptotic correlation functions of the 1D attractive Hubbard model in the partially polarized phase. In particular, we investigate the long-distance asymptotics of the single-particle Green's function, the charge density correlation function, the spin correlation function and the pair correlation function. We present a theoretical confirmation of the existence of the FFLO state in the 1D attractive Hubbard model. In contrast to the continuum case of the 1D attractive SU(2) Fermi gas, here the lattice effects characterized by a parameter β become prominent in the critical exponents of all correlation functions.
The paper is organized as follows. In Sec. 2, we make a brief introduction to the exact solution of the 1D attractive Hubbard model and derive the finite-size correction to the ground state. In Sec. 3, we discuss three types of elementary excitations and the dressed charge integral equations. The dressed charge matrix is derived in the low density limit, which affords access to the conformal dimensions. The long-distance asymptotics of the various correlation functions are studied in Sec. 4, along with the discussion of the FFLO signature in the partially polarized phase. Sec. 5 is reserved for the conclusion.

Exact solution of the Hubbard model and ground state
The 1D fermionic Hubbard model is described by the hamiltonian whereĉ † j,a andĉ j,a , withn j,a =ĉ † j,aĉ j,a , are the creation and annihilation operators of fermions with spin a (a = ↑ or a = ↓) at site j on a 1D lattice of length L. The chemical potential and magnetic field are denoted by µ and B, respectively. Meanwhile u represents the on-site interaction between particles (u > 0 for repulsion and u < 0 for attraction). The hamiltonian (1) is exactly solvable with periodic or open boundary conditions [4,11,12,28,29,30,31,32].
By means of the (nested) Bethe ansatz, the diagonalization of hamiltonian (1) leads to a set of nonlinear algebraic equations, known as the Lieb-Wu equations and written as [11,12] Here Λ is the spin rapidity and k is the quasimomenta of fermions. N and M are the total number of fermions and the number of fermions with down spins. The energy E and momentum P of this model are The magnetization per site m = N −2M 2L . The distributions of the quasimomenta {k i } with i = 1, 2, . . . , N and spin rapidities {Λ β } with β = 1, 2, . . . , M in the thermodynamic limit comply with the string hypothesis. Accordingly, all roots of the Lieb-Wu equations are divided into three categories. These are single real k, k-Λ strings and Λ-Λ strings [33,34,35]. Namely single real k's, the α-th k-Λ string of length m, with for which there are 2m k's accompanied by m spin-rapidities in Λ space, with where j = 1, 2, . . . , m and Λ ′ α m is the real center of the k − Λ string. There is also the β-th Λ-Λ string of length m where j = 1, 2, . . . , m, and Λ m β is the real center of the Λ string. The Λ strings represent the spin wave bound states in the spin sector.
The situation for the ground state is much simplified, where only single real k's and k-Λ strings of length one (m = 1) are permitted. The Λ-Λ strings are suppressed due to the ferromagnetic ordering. This means that the quasimomenta of fermions in bound pairs and the excess unpaired fermions can be respectively written as k b η,± = arcsin(k b γ ± i|u|) and k u β , with η = 1, 2, . . . , M . For simplicity, we have denoted k as k u β and Λ as k b γ . Substituting the above simplified string hypothesis into Eqs. (2) and (3) and then taking logarithms yields the discrete Bethe ansatz equations in the form 2π where N u (N b ) is the number of unpaired fermions (bound pairs) satisfying N u + 2N b = N . The quantum numbers I u j and I b α take the values The ground state energy and momentum are explicitly expressed as We define counting functions, y u L (k u j ) = 2πI u j /L and y b L (k b β ) = 2πI b β /L, which are monotonic increasing functions and satisfy the equations Here ρ γ L (k γ ) (γ = u, b) is the root density of the corresponding quasimomentum and a n (x) = 1 π n|u| (nu) 2 +x 2 . To obtain finite-size corrections in terms of the above root densities when L ≫ 1, we utilize the Euler-MacLaurin formula to obtain the integral equations (up to high orders) where a ′ n (x) is the derivative of a n (x) and Q γ ± (γ = u, b) denote the Fermi points.
It is necessary to introduce the thermodynamic Bethe ansatz (TBA) equations for the derivation of finite-size corrections to the ground state and low-lying excitations. The TBA equations describe the full thermodynamics of the model in the whole temperature regime [33,34,35]. Building on these equations, the equation of states can be derived in terms of densities of single fermions and bound states [35]. In the ground state, the phase diagram consists of five states of the bound pairs and excess fermions in the B-µ plane, see Fig. 1. The T = 0 phase diagram is presented in the left panel of Fig. 1, whereas the right panel presents the low temperature phase diagram obtained via the compressibility Wilson ratio, which is defined by κ Cv /T [35] in terms of the Boltzmann constant k B , the compressibility κ and the specific heat C v . This dimensionless ratio measures the competition between quantum and thermal fluctuations, which exhibits an enhancement near the quantum critical points and thus serves as a powerful tool in determining the low temperature phase diagram. Hereafter we only concentrate on the partially polarized phase IV, which, at zero temperature, presents a novel FFLO-like state on a 1D lattice. At low temperature, the low energy physics of this phase reveals a subtle two-component Tomonaga-Luttinger liquid (TLL), composed of bound pairs and of unpaired fermions.
The ground state TBA equations for the 1D attractive Hubbard model are written as The conformal invariance of a 1D many-body system at T = 0 provides a universality class of criticality in terms of the central charge C of the underlying Virasoro algebra. Indeed, the dimensionless central charge classifies the finite-size scaling form of energies in low-lying excitations. In particular, the C = 1 universality class gives rise to a systematic calculation of the critical exponents which govern the power-law decay of correlation functions in long-distance [39,40,41]. With the help of the root densities presented in Eqs. (16) and (17) and the TBA equations (18) and (19), we derive the finite-size correction to the ground state energy in the form where C = 1 is the central charge for both branches of excitations. In the above equations v u,b are the sound velocities of the unpaired fermions and bound pairs. The sound velocities are defined by Here we have denoted ±Q γ as the Fermi points of corresponding single fermions and bound pairs and the momenta p γ (k γ ) = lim L→∞ y γ L (k γ ).
, and B c2 = 2 + 2|u|, see analytical phase boundaries in [35]. The right panel is a contour plot of the compressibility Wilson ratio R κ W with interaction u = −1 and temperature T = 10 −3 . Both diagrams present five quantum phases. These are I: vacuum; II: fully polarized phase, less than half-filling; III: half-filling phase and fully polarized; IV: partially polarized phase, less than half-filling; V: fully paired phase and less than half-filling.

Low-lying excitations and dressed charge matrix
In order to obtain the conformal dimensions and the critical exponents, the finitesize corrections to the low-lying excitations need to be calculated. The 1D attractive Hubbard model has two branches of excitations in terms of the unpaired fermions and bound pairs. The method used in this work follows the scheme established in Refs. [4,36,37,38]. In general, the low-lying excitations can be realized by combination of three types of elementary excitations, all of which involve the distortion of the Fermi points due to the transformation or variation of the Fermi points. Such changes can be characterized by the changes in the quantum numbers given in Eq. (11).
Excitations of Type I move particles inside the Fermi sea to location j outside the Fermi sea, known as a particle-hole excitation. The lowest particle-hole excitation is described by the change of quantum numbers I γ j close to I γ Here the superscript u stands for the unpaired fermions and b for the bound pairs, respectively. The changes of the total momentum and energy are given by Figure 2: A sketch of distributions for the quantum number I γ j (γ = u, b) in the ground state and excited states. There are three types of elementary excitations: (a) the symmetric distribution for the ground state; (b) a particle-hole excitation near the right Fermi point, the Type I excitation; (c) one more particle is added to the right Fermi point, the Type II excitation; (d) two particles near the left Fermi point are moved to the right Fermi point, the Type III excitation. One may find that all of the possible vacancies for quantum numbers in Type II and Type III excitations are occupied by particles, which are different from the Type I excitation.
where N γ + > 0 (N γ − > 0) stands for the change of the quantum number for corresponding particle near the right (left) Fermi point.
Type II excitations originate from the change of particle numbers of unpaired fermions and bound pairs in the Fermi seas. It is easy to see that the particle number is given by N γ = I γ + − I γ − , and a Type II excitation is characterized by quantum number ∆N γ = N γ e − N γ g where subscripts e and g respectively represent the excited state and ground state.
Type III excitations are caused by the backscattering process, where particles from one Fermi point move to the other one. They are characterized by quantum number In light of the 'continuous' Fermi sea in the Type II and Type III excitations, quantum number ∆N α and ∆D α characterize the changes of the sizes of Fermi seas and the displacements of the center of Fermi seas, respectively.
The calculation of finite-size corrections for Type I excitations is straightforward. Calculations for the Type II and Type III excitations can also be carried out in a systematic way. These calculations are given in the Appendix. Here we summarize the results for the three types of elementary excitations. For convenience in the following calculation of the conformal dimensions, the excitations can be cast into unified finite-size scaling forms of the energy and total momentum, which read For both results there are additional higher order corrections. In these equations we have introduced the notation The superscript t in the above equations stands for matrix transpose, and the quantum numbers ∆D β (β = u, b) obey the relations The dressed charges at the Fermi points Q γ + (γ = u, b) are obtained from the elements of the dressed charge matrix, which satisfy the integral equations We note that the form of the dressed charges is quite different from those for the 1D repulsive Hubbard model [4].
In the ground state, phase V is gapped in the spin sector due to the existence of the bound pairs. However, the system becomes gapless if the magnetic field is greater than the lower critical field, at which bound pairs break. Consequently, phase IV consists of both bound pairs and unpaired fermions. The conformal invariant symmetry enables one to obtain the finite-size scaling forms Eqs. (20) and (23). In what follows we will calculate the conformal dimensions which determine the critical exponents of twopoint correlation functions between primary fields Ô † (x, t)Ô(x ′ , t ′ ) .
We focus on the asymptotics of correlation functions in phase V. For this phase, one expects a power-law decay of the correlation function at T = 0. Meanwhile, at T > 0, the correlation functions should decay exponentially. The conformal scaling dimensions can be read off as where N α ± (α = u, b) characterizes the descendent field from the primary field. It follows that the long-distance asymptotics of the two point correlation functions are given by The dressed charge equations can be simplified in the low density regime, i.e., with small integration boundaries Q γ ≪ 1 (γ = u, b). Here we replace Q γ ± by ±Q γ in the dressed charge matrix, whose elements are calculated for the ground state. Obviously, the dressed charge equations can be separated into two sets of coupled integral equations, composed of Eqs. (28) and (29), and of Eqs. (30) and (31), respectively. By analysing the order of Q γ in the dressed equations Eqs. (28) -(31) we can further obtain asymptotic forms of these equations.
To begin, we substitute Eq. (29) into Eq. (28) to give We further substitute this equation into Eq. (29) to readily obtain Similarly, we have Thus we obtain the dressed charge matrix to leading order, namelŷ By virtue of the TBA equations (18) and (19) with the condition ε γ (±Q γ ) = 0 (γ = u, b), and using standard thermodynamic relations, one can express the cut-off quasi-momenta in terms of particle densities [35] where the density n γ = N γ /L (γ = u, b) must satisfy both conditions n γ /β ≪ 1 and n γ ≪ 1. In the above equations, the lattice parameter is defined by Furthermore, in the strong coupling regime, we have β ≈ 2/|u|. Without losing generality, the approximation used here requires that n γ (γ = u, b) is less than the order of 1/|u|. Meanwhile the condition n γ /β ≪ 1 is required. For the weak coupling regime, numerical calculation enables the confirmation of the asymptotic behaviour of the correlation functions. In Fig. 3, we show the numerical solution of the dressed charge equations (28)(29)(30)(31). We then substitute Eqs. (40) and (41) into Eq. (39). Using the leading order of n γ (γ = u, b) gives the dressed charge in the form With the help of this dressed charge matrix, the conformal dimensions given in Eqs. (32) and (33) in the low density regime can be approximated as These results provide a direct calculation of the asymptotics of the correlation functions.

Asymptotic behavior of correlation functions at zero temperature
We study four types of correlation functions. These are the single-particle Green's function G ↑ (x, t) = ĉ † x,↑ (t)ĉ 0,↑ (0) , the charge density correlation function G nn (x, t) = n † x (t)n † 0 (0) , the spin correlation function G z (x, t) = ŝ z x (t)ŝ z 0 (0) , and the pair cor- . Each of these correlation functions is accessible through choosing suitable quantum numbers ∆N γ (γ = u, b) with respect to the low-lying excitations.
The single-particle Green's function decays exponentially if the magnetic field B < B c1 , for which the external field does not provide enough energy to break up bound pairs. However, if B c1 < B < B c2 then the excess unpaired fermions appear in this gapless phase. In this regime every correlation function satisfies power-law decay [39,40,41,42,42,43]. The single-particle Green's function is determined by the quantum numbers (∆N u , ∆N b ) = (1, 0), which results in (∆D u , ∆D b ) ∈ (Z + 1/2, Z + 1/2). Using Eq. (34) the leading terms of the single-particle Green's function where the critical exponents are given by The leading order term is associated with the quantum numbers (∆ u , ∆ b ) = ±(1/2, −1/2), with the next term coming from (∆ u , ∆ b ) = ±(1/2, 1/2). The coefficients A ↑,1 and A ↑,2 cannot be derived from the CFT approach, yet this does not impede our understanding of the long-distance asymptotic behavior of the correlation functions. Here we have introduced particle densities n ↑ = (N u + N b )/L and n ↓ = N b /L for the particles with up and down spins. We now turn to the charge density correlation function and the spin correlation function, both of which are characterized by quantum numbers (∆N u , ∆N Z). The leading terms are expressed as where the critical exponents are given by Here the constant terms n 2 and m 2 z originate from quantum numbers (∆ u , ∆ b ) = (0, 0), while the second, third, and fourth terms come from (∆ u , ∆ b ) = ±(1, 0), ±(0, 1) and ±(−1, 1), respectively.
Last but not least we discuss the pair correlation function G p (x, t), which is described by the quantum numbers (∆N u , ∆N b ) = (0, 1), allowing (∆D u , ∆D b ) ∈ (Z + 1/2, Z). We find that the leading terms for the pair correlation function are with critical exponents Here the first and the second terms are associated with (∆ u , ∆ b ) = ±( 1 2 , 0) and ±( 1 2 , 1), respectively. The asymptotic behaviour of the correlation functions reveals an important manybody correlation nature and lattice effect, which is apparent in the lattice parameterdependent critical exponents. Moreover, the leading order terms of the pair correlation function and spin correlation function reveal spatial oscillating behavior in their longdistance asymptotics. The pair correlation function oscillates with a wave number ∆k = π(n ↑ −n ↓ ). So does the spin correlation function with 2∆k. Our analytic results provide a confirmation of the previous numerical observations of this oscillatory nature [20,21,22,23,24]. The oscillation stems from the backscattering processes in the two Fermi seas, where the imbalance between the densities of spin-up and spin-down particles results in the mismatch of their Fermi surfaces. It is interesting to see that the spatial oscillation in the 1D attractive Hubbard model is the feature of the Larkin-Ovchinnikov phase predicted in [14]. We find that the oscillation terms in the spin and pair correlation functions arise from the Type III elementary excitations (backscattering process). Our theoretical result for the wave number shows good agreement with the numerics [20], where the numerical wave number is almost ∆k = π(n ↑ − n ↓ ). This comparison implies that the coefficient A p,1 is much larger than A p,2 . Additionally, with respect to the coefficients in the spin correlation function, A z,1 is much larger than A z,2 and also A z,3 . Finally, we point out that the parameter β defined in Eq. (42) represents the lattice effect, and thus distinguishes the Hubbard model from its continuum limit -the 1D attractive SU(2) Fermi gas [35].
Applying Fourier transformation to the above correlation functions allows the derivation of their counterparts in momentum space [38]. For the equal-time correlation function, where ∆ ± = ∆ u ± + ∆ b ± . The Fourier transformation is thus Here the conformal spin and the exponent are given by s = ∆ + − ∆ − and ν = 2(∆ + + ∆ − ) − 1. Consequently, the Fourier transforms of the equal-time correlation functions near the singularities k ≈ k 0 are expressed as where 2s z = 2s nn ≈ 0, ν z = ν nn ≈ 1, We plot the pair correlation function in momentum space in Fig. 4. Notably, the correlation functions in momentum space shown in Eqs. (55) to (58) are valid only in the vicinity of the wave numbers k 0 , i.e., k ≈ k 0 . Fig. 4 shows that G p (k) has a singularity in non-zero momentum in a partially polarized phase. This qualitatively agrees with the numerical result given in [20]. One should notice that this plot is correct only if k ≈ π(n ↑ − n ↓ ), where the extrapolation is used for the purpose of better visualization.

Conclusion
In this paper we have investigated four types of correlation functions for the 1D attractive Hubbard model at zero temperature. The finite-size corrections to the ground state and the low-lying excitations are derived explicitly. Based on these corrections to the momentum and energy, we have applied CFT in the study of long-distance asymptotic behavior of the correlation functions. The critical exponents have been obtained explicitly in this way. In contrast to the Fermi gas, these asymptotics of correlation functions essentially depend on the parameter β which represents the lattice effect. We have found that the spin and pair correlation functions have spatial oscillations with frequencies 2π(n ↑ − n ↓ ) and π(n ↑ − n ↓ ), respectively. From the perspective of CFT, this type of oscillation is induced by the backscattering process of unpaired fermions and bound pairs among the Fermi points. This gives a microscopic origin of the FFLO pair correlation in the 1D system. Using the Fourier transform, we have also derived the correlation functions in momentum space. The pair correlation is singular at the mismatch point k = π(n ↑ − n ↓ ), which confirms the frequency found by numerical methods [20]. Meanwhile the correlation functions at zero temperature display powerlaw decay in the partially polarized phase IV. This suggests an analog of long-range order in a 1D many-body system and thus demonstrates the existence of a superconducting state. We further point out that the dressed charge matrix in this CFT approach can be numerically resolved for arbitrary interaction strength u < 0 (see fig. 3). It follows that one can calculate the conformal dimensions and critical exponents for arbitrary interaction strength. Our results provide benchmark physics of 1D strongly correlated fermions on a lattice which may be testable in current ultracold atomic experiments [5,6,7,8,9,10].
Then noticing that the counting function connects the momentum k α and quantum number I α , we similarly introduce a Taylor expansion in the last equation, with result where y α = lim L→∞ y α L (k α ) = p α (k α ) is the counting function in the thermodynamic limit, and we have used the definition for the sound velocity.
In light of the zero momentum in the ground state, we can write down the change in total momentum, where N α + = I α +,e −I α −,g and N α − = I α −,g −I α +,e arise from the change in distribution of quantum numbers close to the right and left Fermi points, respectively. The subscripts e and g represent the excited state and ground state.
On the other hand, the derivations for the changes in energy and total momentum of Type II and III elementary excitations are rather complicated. We give an outline of these calculations here.
Prior to discussion of Type II and III excitations, we write down the root density in the thermodynamic limit, where for simplicity we have denoted π −π dk a 1 (k b − sin k), and introduce the matrix expression for the integral kernels. We then use a simpler expression for Eq. (A.5), where ⊗ stands for an integral on k β over interval [Q β − , Q β + ]. The low-lying excitation corresponds to small variations in the Fermi points from the ground state case. This allows expansion of the energy per site with respect to the changes of the Fermi points to leading order, is the energy per site, and ±Q α denotes the ground state Fermi points. The subscript G implies Q u ± = ±Q u and Q b ± = ±Q b . In this expansion, the first order derivative ∂eE vanishes due to the fact that the Fermi points minimize the energy in the ground state.
With the help of the TBA equations and integral equations of root densities, the second order derivative is given by Substituting this result and the sound velocity with additional higher order correction terms. In order to rewrite Q γ + −Q γ and Q γ − +Q γ in terms of ∆N γ and ∆D γ , we introduce the notation ν γ = stands for the number of unpaired fermions (bound pairs). D γ is the position of the center of the Fermi sea for unpaired fermions (bound pairs). Their explicit expressions are given by with θ(x) = 2 arctan(x).
The total differential of ν γ (γ = u, b) with respect to Q β ± (β = u, b) in the vicinity of the ground state is We furthermore denote Then Eq. (A.15) can be rewritten as a vector equation, Here we need to calculate According to Eq. (A.8), it is straightforward to derive which together with the dressed charge equations yields We also introduce the notation In the next stage, we similarly consider the total differential of δ γ with respect to Q β ± in the vicinity of the ground state, Similarly, we introduce the notation in terms of which Eq. (A.25) can be rewritten as the vector equation We further calculate Now we introduce a set of new integral equations similar to the dressed charge equations, of the form where η = u, b and the integral kernel matrixK T is written aŝ With the help of the above integral equations for σ αβ (k β ) and Eq. (A.19), we obtain We now further defineσ whose inverse is given by Finally, we express Eq. (A.11) in terms of dQ γ ± (γ = u, b), which gives, subject to higher order correction terms, In light of the definitions ofŜ,ρ, ∆ N , ∆ D, with ν γ = N γ /L and δ γ = D γ /L (γ = u, b), this last equation gives our final result, for the leading finite-size correction term.
The change in total momentum caused by Type II and III excitations is given by After the above lengthy calculations, the changes in the energy and total momentum of the three types of elementary excitations have been derived. The change in energy shown in Eqs. (22) and (A.49) is summarized as Eq. (23). The change in total momentum shown in Eqs. (21) and (A.50) is summarized as Eq. (24).