Algebraic approach to Bose Einstein Condensation in relativistic Quantum Field Theory. Spontaneous symmetry breaking and the Goldstone Theorem

We construct states describing Bose Einstein condensates at finite temperature for a relativistic massive complex scalar field with $|\varphi|^4$-interaction. We start with the linearized theory over a classical condensate and construct interacting fields by perturbation theory. Using the concept of thermal masses, equilibrium states at finite temperature can be constructed by the methods developed in arXiv:1306.6519 and arXiv:1502.02705. Here, the principle of perturbative agreement plays a crucial role. The apparent conflict with Goldstone's Theorem is resolved by the fact that the linearized theory breaks the $U(1)$ symmetry, hence the theorem applies only to the full series but not to the truncations at finite order which therefore can be free of infrared divergences.


Introduction
In this paper we shall analyze the perturbative construction of a Bose-Einstein condensate for a relativistic charged scalar field theory at finite temperature.
The first experimental realization of Bose-Einstein condensation (BEC) in dilute vapours of alkali atoms has been obtained few years ago [And95,Bra95,Dav95]. These works have pushed a lot the theoretical and experimental investigation of this phenomenon.
Usually, Bose-Einstein condensation is discussed in the realm of non relativistic quantum theories. (See e.g. [DGPS99] and references therein.) In that case, a Bose Einstein condensate is the effect where weakly interacting Bosons populate the ground state and thus dominate over the higher excitations. Below a critical value of the temperature, in the limit where the number of bosons tends to infinity keeping the density finite, almost every Boson will be located in the ground state. In the large N limit, the wave function of the ground state satisfies the Gross-Pitaevskii equation [Gr61,Pi61]. A rigorous analysis of this fact can be found for instance in [LSY01b,LSY01a].
Instead of considering a quantum system in a box and the limit where the box tends to cover the full space, we deal directly with quantum fields over the full spacetime. We shall work in a relativistic setting because the construction of the algebra of interacting observables in this case is free from infrared divergences 1 thanks to the causal property of the theory, see e.g. [BB94,BR14] 2 . Furthermore, there are physical systems where the relativistic nature of fundamental physics is manifest together with the phenomenon of condensation. Here we have in mind the possible existence of Boson stars at cosmological level [BMZ16], the condensation of relativistic particles like photons [So13] recently detected in a photon gas and condensation phenomena in high energy physics as e.g. the quark gluon plasma [BN77,Sa11].
The traditional construction of non-zero temperature equilibrium states (KMS states) for perturbatively defined interacting quantum field theories suffers even in the massive case from spurious infrared divergences at higher loop order (see [Al90,St95]). It was recently shown how these divergences can be circumvented [Li13,FL14]. This construction amounts to an adaptation of formulas derived in rigorous statistical mechanics by Araki [Ar73] to the framework of perturbative QFT.
To apply these methods to BEC at finite temperature one has to deal with the problem that correlations in the unperturbed state decay slowly. In the case of a massless scalar field this problem could be avoided by taking into account that the interaction produces at finite temperature a thermal mass. If this term is included in the free theory, the correlations of the unperturbed state decay sufficiently fast [DHP17]. In the case of BEC this is in conflict with the existence of a Goldstone mode induced by the spontaneous breakdown of the U (1) symmetry of the model.
We solve this problem in the following way: We linearize the theory around a classical solution. The linearized theory breaks the U (1) symmetry and shows nonvanishing thermal masses, hence the perturbative construction works as in the massless model treated in [DHP17]. The U (1) symmetry is recovered for the full theory which then has a massless mode in agreement with the Goldstone theorem.
The paper is organised as follows. In the next section we shall briefly recall the framework of pAQFT, the key steps in the construction of the KMS state performed in [FL14] and few facts about the principle of perturbative agreement [HW05] by which we can move the thermal mass into the free theory [DHP17]. The third section contains the perturbative analysis of the massive complex scalar field with a |ϕ| 4 interaction, expanded around a solution with a nonvanishing condensate. We shall then discuss the construction of the interacting state at finite temperature 1 The infrared divergences at zero temperature in the nonrelativistic case have been discussed by Benfatto [Be94] and by Di Castro group [PCDS04]. 2 The construction of an algebra of interacting nonrelativistic bosons was a longstanding open problem which was recently solved by Buchholz [Bu18] using the concept of the so-called resolvent algebra [BG08] over the condensate and we analyze its adiabatic limit. The forth section contains the discussion of the formation of the condensate in connection to the spontaneous symmetry breaking. We shall actually see that the symmetry is effectively broken in the background theory while it is recovered in the exact theory where it is thus spontaneously broken.

Equilibrium states for interacting quantum field theory
In this section we briefly review the formalism of perturbative algebraic quantum field theory. This framework combines renormalized perturbation theory with concepts of algebraic quantum field theory [HK64,Ha92]. The basic step is an assignment of the algebra generated by the observables of the theory to each region of the spacetime. Physical informations, like locality, are then stored in the relations between algebras labelled by different regions of the spacetime. In the case of an interacting theory treated with perturbation theory, the elements of the algebras associated to each region are given as formal power series in the coupling constant with values in suitable * −algebras. States are constructed in a second step as linear functionals on the algebra of observables which via the GNS construction then provide representations of the elements of the algebra by operators on a state space. Renormalization in this framework is automatically independent of the state; moreover, infrared problems do not occur in the construction of the algebra. They may become visible in the construction of states where they indicate physical properties of the system.

Perturbative construction of the interacting quantum field theory and the adiabatic limit
We shall here briefly recall the basic elements of the perturbative construction of the φ 4 scalar field theory propagating on a four dimensional Minkowski spacetime (M, η) where the metric η has the signature (−, +, +, +). The Lagrangian is where m is a positive mass and λ the coupling to the non linear perturbation. In a first step we construct the algebra corresponding to the free theory (λ = 0). We label the elements O of the algebra by functionals on the classical configuration space which in our case is the space of smooth functions φ on Minkowski space, where f n is a compactly supported distribution on M n which is symmetric under permutations of the arguments. The spacetime support supp O of O is the smallest closed subset G of Minkowski space such that suppf n ⊂ G n for all n. O is called regular if all f n are smooth, and it is local if all f n are of the form with a partial differential operator P with smooth compactly supported coefficients.
The product of regular observables is defined in terms of the commutator function ∆, which is the retarded minus advanced fundamental solution of ( − m 2 )φ = f . In the sense of generating functionals it is given by where φ(f ) = φ(x)f (x)dx for a compactly supported test function f . Due to singularities of ∆, the product cannot be extended to nonlinear local functionals as e.g. f (x)φ(x) n dx with n > 1 and a test function f . This well known problem is, as usual, circumvented by replacing these functionals by so-called normal ordered functionals. In our framework this means the following. Consider the Klein-Gordon equation with spacetime dependent mass m(x) on a globally hyperbolic spacetime. Let H be a symmetric bisolution of the form where U, V and W are smooth symmetric functions of 2 spacetime points x and y. U and V depend only on the geometry and on the mass near to the geodesic connecting the arguments, and σ (x, y) = σ(x, y) + i (t(x) − t(y)). σ is the square of the geodesic distance between x and y, equipped with the appropriate sign for spacelike and timelike separation, respectively, t is a time function and ξ a lenghtscale. Such a bisolution is called a Hadamard function [KW91]. According to Radzikowski [Ra96], a Hadamard function H can be characterized as a symmetric bisolution with the property that the wave front set of H + i 2 ∆ satisfies a positivity condition (microlocal spectrum condition [BFK96]). Examples of Hadamard functions are the symmetric parts of the 2 point functions of vacuum and KMS states.
Given a Hadamard function H, normal ordering is a linear map defined by The -product then is interwined with the so-called Wick product H , Due to the smaller wave front set of H + i 2 ∆ compared to ∆, the Wick product can be extended to a larger class of functionals O. The product among these objects is well defined as soon is the wave front set [Ho03] of the distribution f n ∈ E (M n ) and where V ± n denotes the closure of the forward/backward light cones in T * M n , These functionals are called microcausal, and we denote their set by F µc . It contains in particular the local functionals and their pointwise (classical) products, the multilocal functionals. We can now extend the algebra of observables by normal ordered microcausal functionals and define their product by Note that the enlarged algebra A does not depend on the choice of the Hadamard function. Only the labeling by the functional depends on H. Note furthermore that the normal ordered functionals are, in general, no longer functionals on the configuration space due to the singularities of the Hadamard function. The standard normal ordering in Fock space with respect to annihilation and creation operators is obtained if one chooses the symmetrized vacuum 2 point function ∆ 1 as the Hadamard function. It has the nice feature that the vacuum expectation value ω 0 (: O : ∆ 1 ) of a normal ordered functional : O : ∆ 1 coincides with the evaluation of the functional O at φ = 0, A corresponding formula holds for any quasifree (Gaussian) state ω, if one uses its symmetrized 2-point function as the Hadamard function, in particular for KMS states.
This choice of normal ordering, however, is problematic when one wants to identify observables in different states or under the process of renormalization. Then a preferred Hadamard function with W = 0 is better behaved, as first discussed by Hollands and Wald [HW01]. In the case of a generic curved background that preferred function is in general not well defined but in the case of Minkowski space it was explicitly constructed in appendix A of [BDF09]. It is unique up to the choice of a length scale.
The field equation is not yet implemented into the algebra A. The algebra contains instead an ideal generated by normal ordered functionals which vanish on solutions. The so-called on shell algebra is obtained by taking the quotient with respect to this ideal. This quotient is faithfully represented on Fock space and coincides with the standard algebra of the free field.
The off shell algebra A, however, is better behaved under the time ordered product which is used for the incorporation of interaction.
Interacting fields can be constructed by means of causal perturbation theory, a method of renormalization elaborated by Epstein and Glaser [EG73] on the basis of ideas of Stueckelberg [St51,St49] and Bogoliubov [BS76]. It was further developed by Scharf and collaborators (see e.g. [Sc89]). It is also the basis for a treatment of interactions on curved spacetimes [BF00,HW01] where other versions of renormalization did not work. Its main idea is the construction of time ordered products of interaction Lagrangians. In the work of Epstein and Glaser these products are operator valued distributions on Fock space. One uses the fact that the time ordered product for noncoinciding points agrees with the operator product in the appropriate order. Renormalization then consists in extending these distributions to coinciding points. By induction with respect to the number of factors one can show that this extension is always possible and unique up to the addition of a further interaction Lagrangian in each order. This corresponds precisely to the freedom in the choice of renormalization conditions known from other versions of reormalization.
In pAQFT one uses a version of causal perturbation theory which is independent of the choice of a state space. There the time ordering operator is a linear map T from the algebra of multilocal functionals (with respect to the pointwise (classical) product) 3 to the algebra of A. On regular functionals it is determined by with the Dirac propagator ∆ D (the mean of retarded and advanced propagator) and the Feynman propagator ∆ H F = ∆ D − iH associated to the Hadamard function H. It satisfies the causal factorization condition if the support of F does not intersect the past of the support of G. To extend T to multilocal functionals one fixes a Hadamard function H and characterizes the extensions by the initial conditions T 1 = 1 and T F =: F : H for local functionals F vanishing at φ = 0. The causal factorization condition then fixes the map T on n-local functionals by its values on k-local functionals for k < n up to a local functional. In order to reduce the ambiguity in this extension we choose a Hadamard function with W = 0 as discussed before, so only the dependence on the scale ξ remains for local functionals.
Given an interaction Lagrangian L I where L I (ϕ) is a translation invariant section of the jet bundle constructed over the field configuration ϕ, we consider the local functional with the test function g.
where the inverse is understood in terms of the product, and the interacting fields are given in terms of the Bogoliubov map (also called Møller operator) which extracts the contributions of S V (µF ) linear in µ 3 The algebra of multilocal functionals with respect to the pointwise product is isomorphic to the symmetric tensor algebra of local functionals vanishing at ϕ = 0. This fact has been proved in [FK13].
Interacting observables can now be represented as elements of the algebra generated by the relative S matrices S V (F ) where the square brackets denote the set of linear combinations of products of the elements inside the brackets. We observe that the association of spacetime regions O to algebras A I (O) forms a net of subalgebras of A[[λ]] in the sense of the Haag-Kastler axioms of algebraic quantum field theory. In the following we omit the symbol for the product within A[[λ]] and replace it by juxtaposition.
The last step in constructing the interacting theory is the removal of the cutoff g from the interaction lagrangian V in (3). This can be done taking the adiabatic limit g → 1. At algebraic level, as soon as the interacting observables are supported on a compact region O, this limit can be taken over larger and larger regions where the cutoffs are equal to 1, further details can be found in [BF00]. Here we can make this construction more explicit making use of the time slice axiom and the causal properties of the S−matrix. Actually, both the S−matrix and the relative S−matrix satisfy the following causal factorisation property valid for A, B, C ∈ F loc where A C means that A is later than C in the sense that supp(A) ∩ J − (supp(C)) = ∅, where J ± (O) denotes the causal future/past of O. This causal factorisation property implies that If g, g coincide on J + (O) ∩ J − (O) the correponding local functionals V, V differ by The limit g → 1 can now be taken at the algebraic level.
Consider now a Cauchy surface Σ = t −1 (0) where t is the time coordinate of a standard Minkowski coordinate system which is fixed once and forever. An neighborhood of the Cauchy surface Σ is Interacting fields satisfy the time slice axiom, see [CF09], namely, for every A ∈ A I (O) it exists a C ∈ A I (Σ ∩ J(O)) such that where W ∈ F µc vanishes on solutions φ. Here J(O) = J + (O)∪J − (O). Hence the on shell algebras B I (O), obtained by taking quotients with respect to the ideal generated by these elements, are subsets of B I (Σ ). Thus to construct a state for the interacting algebra B I (M ) it suffices to construct it for B I (Σ ). We therefore choose a cutoff function g of the form where now the time cutoff is realized by χ(t) which is a smooth function which is equal to 1 for t > − and 0 for t ≤ −2 . Furthermore, h is a space cutoff which is compactly supported on Σ.
To obtain a state in the abiabatic limit, it is sufficient to consider the limit where h tends to 1 keeping fixed the time cutoff χ.

Interacting KMS states and the adiabatic limit
Equilibrium states are characterized by the Kubo-Martin-Schwinger (KMS) condition, see [HHW67]. This condition yields canonical Gibbs states when they are well defined, but remains meaningful also for infinitely extended systems where the Gibbs formula can no longer be used [HHW67,Ha92]. We recall here the definition. A state ω over a C*-algebra B satisfies the KMS condition with respect to the one parameter group of * −automorphisms τ t at inverse temperature β if for every A, B ∈ B ω(Aτ t B) is an analytic function for Im(t) ∈ (0, β), continuous at the boundary and if ω(Aτ iβ B) = ω(BA).
A state which satisfies the KMS condition at inverse temperature β is called β-KMS state. It is automatically invariant under τ t and satisfies similar relations for n-point functions.
If one uses the concept of KMS states for general * -algebras one has to enrich the definition by some of these properties. See e.g. Definition 1 in [FL14] for an extended discussion.
Let τ t denote the one parameter group of * −automorphisms of B induced by the action of Minkowski space time translations Consider now the following two-point function The quasifree state ω β constructed out of this two-point function for the free theory is a KMS state at inverse temperature β with respect to time translations. This state is easily described by using normal ordering with respect to its symmetrization H β as The interacting time evolution whereas the free evolution is τ t (S V (F )) = S Vt (F t ). To construct the KMS for the interacting theory we have to relate the free and interacting time evolution. The causal factorisation property (4) implies that from positive real numbers to unitary elements of B I defines a cocycle which intertwines the free and interacting time evolution. The cocycle relation and its infinitesimal generators are where, in the case of V as in (3) with g as in (7) it turns out that Hence H I and H I play the role of the interacting Hamiltonian and the interacting Hamiltonian density. We stress that due to the smearing in time, H I (x) is a well defined formal power series with coefficients contained within the algebra of the free theory. For any β-KMS state ω β of the free theory, like the quasifree state (9) constructed with the two-point function (8) , and any spatial cutoff described by the test function h on Σ we obtain a β−KMS state of the theory with interaction H I (h) with respect to the evolution τ V t observing that t → ω β (AU (t)) for every A ∈ B I (Σ ) can be analytically continued to Imt ∈ [0, β]. Hence, defines a β−KMS state with respect to τ V t , as proved in [FL14]. Furthermore, the expectation values in the state ω β,V can be computed by the following formula Here ω β T denotes the truncated functional associated to ω β . As shown in [FL14], the limit h → 1 can now be taken provided the truncated n−point functions decay sufficiently fast for large spatial separations. Furthermore, the obtained state does not depend on χ anymore. In this way one obtains the correlation functions for an interacting field in thermal equilibrium in the case of a massive theory.

Principle of perturbative agreement -Massless case
If the linearized theory is massless, the limit h → 1 in (11) cannot be taken because the decay of the n−point function is too slow. However, in the case of a φ 4 theory, it is possible to use a similar construction [DHP17]. The idea is to modify the splitting by adding an artificial mass M to the background theory and subtracting it in the interacting Lagrangian, namely We see that the interaction Hamiltonian density in the KMS-state contains a mass term with a positive coefficient as long as M 2 < M 2 β with the thermal mass Under this condition the interaction Lagrangian remains convex and possesses a single stationary point at φ = 0. We recall that the thermal mass is is a renormalization constant and it depends on the length scales ξ in (2). If ξM is equal to 1 then M 2 β vanishes in the limit β → ∞ and M 2 β = T 2 /12 + O(M 2 ). We finally observe that the theories constructed with the two different splittings are equivalent thanks to the principle of perturbative agreement, which has been shown to hold in [HW05], see also [DHP17]. Further details about the validity of this principle are collected in Appendix A. We finally recall that, if the interaction Lagrangian is quadratic in the field Q = δm 2 ϕ 2 dx 4 and if it corresponds to a perturbation of the mass m of the free theory to √ m 2 + δm 2 , the equilibrium state constructed as in (10) is the KMS state at inverse temperature β with perturbed mass. This last observation has been proved in Theorem 3 of [Dr19] and it shows that perturbative agreement is compatible with the construction of equilibrium states discussed in [FL14].

Massive complex scalar field
We discuss the equilibrium states at finite temperature β when a chemical potential µ is present in the theory. We start discussing the free theory (λ = 0), afterwards we analyze the same kind of states when a condensate for the interacting theory is chosen.

Condensate in the free theory
Let us start discussing the condensate for a free massive complex scalar quantum field theory propagating in a Minkowski spacetime. Let us denote by ϕ the associated field configuration. Its equilibrium states with inverse temperature β > 0 and chemical potential µ, |µ| < m are quasifree ("Gaussian") states with the 2-point functions and ω β,µ (ϕ(x)ϕ(y))) = ω β,µ (ϕ * (x)ϕ * (y))) = 0. They satisfy the KMS condition with respect to the time evolution where e 0 denotes the unit vector in time direction. For chemical potentials µ = ±m there exist many equilibrium states at fixed temperature. The pure phases can be obtained from ω β,±m by an automorphism which is generated by where c is a harmonic function, ∆c = 0, with the spatial Laplacian ∆. Notice that the states ω ± β,c .
= ω β,±m • γ ± c are KMS states w.r.t. τ t,±m given in (16). The theory possesses an internal U (1) symmetry which is spontaneously broken in ω ± β,c . Hence, an interesting observable for these states is the current density where this is seen as an element of A. We observe that, in view of the symmetry of the Hadamard coefficients U , V i of V = V n σ n in (2), and W i of W = W n σ n = (ω 0 − H − i∆/2), ω 0 being the two-point function of the vacuum state, we have that :φ * ϕ − ϕ * φ : H =:φ * ϕ − ϕ * φ : H+W . Hence, J 0 (f ) can be seen as the current normal ordered with respect to vacuum state. At zero temperature, all the states with chemical potential |µ| < m coincide with the vacuum. Hence, the vacuum expectation value of the smeared current J 0 (f ), or better of the charge density J 0 (x), vanishes in that limit. Thus any non vanishing charge density in the limit T = β −1 → 0 of vanishing temperature requires a condensate. We find in particular that Here the first term ρ cr (β) is the critical charge density at a given temperature T = β −1 , and the second term is the condensate. Actually we might have a non vanishing condensate only if |ω ± β,c (J 0 )| − ρ cr (β) = 2m|c| 2 > 0. We furthermore observe that ρ cr (β) > 0, it is monotonically decreasing in β, it diverges for β = 0 and tends to 0 for large β. Hence, at a fixed charge density ω ± β,c (J 0 ), there is a critical temperature T cr > 0 obtained as ρ cr (1/T cr ) = ω ± β,c (J 0 ) such that only for T < T cr a condensate can be formed.

Massive complex scalar field with ϕ 4 interaction over the condensate
In this section we start discussing the perturbative construction of the ϕ 4 interacting theory over the condensate in the Minkowski spacetime. The Lagrangian of the theory we are considering is thus where ϕ is a complex scalar field. Following a similar procedure presented in section III of [ABS08], we expand L around a real classical solution φ which represents the condensate. Hence where µ is again the chemical potential, x 0 is a fixed Minkowski time and ψ is a complex scalar field which describe the perturbations. Its real and imaginary parts are denoted by ψ 1 and ψ 2 and thus ψ = ψ 1 + iψ 2 .
The Lagrangian density can now be written as a sum of contributions homogenous in the number of fields ψ as follows The term L 1 vanishes because φ is chosen to be a stationary point for the classical action L 0 dµ. In the following, we shall choose a non vanishing φ to describe the condensate, we discuss the quantization of the linearized theory (L 2 ) and finally we use perturbation theory over the linearized theory to take into account L 3 + L 4 .

The condensate in the vacuum theory
We look for the case of a translation invariant background φ. Then, the kinetic term in L 0 has no effect and φ is a stationary point for and only one real, positive and translational invariant background solution φ is thus available for µ 2 > m 2 . We notice that for fixed µ 2 > m 2 the background value of the field φ is of order 1/ √ λ. In this case, we observe that L 2 does not depend on λ, L 3 is of order √ λ while L 4 is of order λ. In the next we shall construct the interacting field theory with perturbation methods considering L I = (L 3 + L 4 ) the interaction Lagrangian. Hence, the solution we shall obtain will be a formal power series in √ λ.

Linearized theory
The first step to construct the quantization of ϕ is the analysis of the linearized equations of motion for the fluctuations (ψ 1 , ψ 2 ) around φ. They have the form Notice that if (21) holds, M 2 1 = 2(µ 2 − m 2 ) and M 2 2 = 0. Hence, we assume M 1 > M 2 ≥ 0. Let us introduce We observe that the equations (22) forψ = (ψ 1 , ψ 2 ) can be written in a compact form Dψ = 0 where D is given in terms of the standard Pauli matrices Notice that DD = (− + M 2 ) 2 − (δM 2 ) 2 + 4µ 2 ∂ 2 0 I. The retarded and advanced propagators of the theory can be obtained as where (DD) R and (DD) A are the retarded and advanced fundamental solutions of DD. Let us thus study hence a massless mode is present in this system as expected by the Goldstone theorem. However, if the linearized theory is not in a ground state, it could happen that the normal ordered interaction Lagrangian with respect to the state, as in (14), contains quadratic terms that could contribute to the masses of the fluctuations. If we use the formula valid for pairwise different x 1 , . . . , x n , a couple of times, we get We can construct ∆ A just changing i → −i , while the Feynman propagator ∆ F is obtained substituting (p 0 + i ) 2 with p 2 0 + i and multiplying by i. Finally, the commutator function is With ∆ at disposal the quantum product can be given as in section 2.1, in this way we obtain the * −algebra of field observables A 0 . The analog of the Hadamard singularity H (2) for this theory can be given. The form of some of the corresponding Hadamard coefficients are discussed in Appendix B. The extended * −algebra of field observables A containing Wick polynomials normal ordered with respect to H is obtained as in 2.1.

KMS states for the linearized theory
In view of the decomposition of the field ϕ given in (20) the action of τ t on ψ as time translation is equivalent to the action of τ t,µ on ϕ as given in (16). Hence, having the causal propagator of the linearized theory at disposal, we can construct the two-point function of the quasifree β−KMS state with respect to time translation τ t of the ψ fields as we have that the two-point function of the quasifree β−KMS state ω β,ψ is in position space Recalling the form of ω ± in (25) we notice that if M 1 > M 2 > 0 this means that no infrared divergences are present in ω β,ψ if M 2 > 0. The two-point function of the ground state of the ψ i theory (keeping the condensate φ = 0) can be obtained taking the limit β → ∞ of (27). Hence, to study expectation values in the state ω β,ψ of observables normal ordered with respect to the vacuum ω ∞,ψ we consider W = ω β,ψ − ω ∞,ψ and we obtain We observe that in the coinciding point limit, the off diagonal expectation values are vanishing and, introducing 2δω 2 .

Thermal masses
Having analyzed the equilibrium state of the free theory on field observables A, the next step in the construction of an equilibrium state for the interacting theory will be an application of the analysis given in [FL14] and summarized in section 2.2, namely to use (10) starting with a quasifree state whose two-point function is given in (27). However, we expect that the limit h → 1 cannot be directly taken because, as discussed above, if (21) holds, the mass M 2 given in (23) vanishes, hence, for vanishing spatial momentum, ω 2 − is also vanishing. This implies that, the various propagators of the linearized theory diverge for p → 0. Hence, in agreement with Goldstone theorem a massless mode is present in this case. This implies a slow decay in the connected n−point functions constructed with ω β,µ given in (27).
In order to cure this problem we use a different splitting of the Lagrangian into the free and interacting part. Actually, we add a virtual mass m 2 v to the linearized fields and we remove them in the interaction Lagrangian. More precisely, the Lagrangian of the free theory is now while the modified interaction Lagrangian is The elements of the interacting algebra are now given in terms of two parameter λ and m v . More precisely, keeping µ fixed, as in (21), they are formal power series in √ λ with coefficients depending on m 2 v , which can be understood as a partial resummation of the original perturbative expansion. The advantage of this new expansion is in the fact that the coefficients remain finite in the adiabatic limit, when they are evaluated in the state representing the condensate at finite temperature. We furthermore observe that the principle of perturbative agreement discussed below implies that the final theory does not depend on this extra parameter m v .
Let H ψ,β be the symmetrized two-point function of the β-KMS given in (27), we observe that if m v is chosen to be sufficiently small, the interaction Lagrangian normal ordered with respect to H ψ,β is again convex. To see this in detail, let T be a time ordering operator such that T F =: F : H where H is the distinguished Hadamard function constructed in (50). We have up to a choice of renormalization freedom (the lenghtscale ξ in (50) chosen in such a way that : |ψ| 4 : H =: |ψ| 4 : H ∞,ψ where H ∞,ψ is the symmetrized two-point function of the vacuum obtained taking the limit β → ∞ in (27)) where C is a constant which can de discarded and the two thermal massese m β,i have been computed above in (29) and (30) hence T ( 1 4 |ψ| 4 − m 2 v 2 |ψ| 2 ) remains convex, provided m 2 v < λ(3m 2 β,1 + m 2 β,2 ) and m 2 v < λ(3m 2 β,2 + m 2 β,1 ). For this reason it is expected that the stability properties of the theory are not altered adding the virtual masses m v in the free theory.

Condensate and perturbative agreement
We need to check that the Wick monomials in the interaction Lagrangian originally constructed over the linearized theory L 2 , are not corrected because of the new splitting. In other words, we prove that the principle of perturbative agreement holds also when a condensate is present. Let us recall the form of the equation of motion forψ = (ψ 1 , ψ 2 ) given in (24) and consider the preferred Hadamard function H M 2 ,δM 2 ,µ , with a lengthscale ξ, associated with this operator constructed in Appendix B.
We prove now that the time ordering operator T M,δM,µ (F ) =: F : H M 2 ,δM 2 ,µ satisfies the principle of perturbative agreement. To this end consider the 2 × 2 matrix Ψ = {ψ i ψ j } i,j∈{1,2} , following the discussion presented in Appendix A we want to prove that ∆Ψ = γT 0,0,µ Ψ − T M,δM,µ Ψ vanishes, where γ is the map which intertwines T 0,0,µ (ψ i (x)ψ j (y)) to T M,δM,µ (ψ i (x)ψ j (y)). Formally, indicating with the supscript c the quantities referred to the condensate (M, δM, µ) and with the supscript 0 those referred to the vacuum (0, 0, m) we have to compute where H c/0 F are the time ordered/Feynman propagator associated to the Hadamard functions H c/0 . By power counting we notice that all the contributions larger than order two in (D c − D) are removed from ∆Ψ by renormalization. In order to check if there is a finite reminder after this renormalization we analyze the form of the Hadamard singularity H M 2 ,δM 2 ,µ given in Appendix B, we remove the contributions of order lower than the third in x i from H M 2 +x 1 ,δM 2 +x 2 ,µ+x 3 before computing the coinciding point limit. Let us recall the form of some Hadamard coefficient given in Appendix B. From equations (52) and (53) we have U = cos(µx 0 )I − iσ 2 sin(µx 0 ) and hence, in H M 2 +x 1 ,δM 2 +x 2 ,µ+x 3 , U/σ does not depend on x 1 and x 2 and the contributions in x 3 larger than the second order vanish in the coinciding point limits. Similarly, the contributions larger than second order vanish also in V log( σ ξ 2 ). As for the mass perturbations we thus have that the Wick monomial Ψ n computed with respect to the Hadamard parametrix do not change under the action of the map which intertwines the time ordering constructed with two different sets of parameters M, δM, µ.

Construction of the condensate, cluster estimates
To construct the state at finite temperature over the condensate, we follow the construction given in [FL14] and summarized in section 3.2. In particular, at fixed spatial cutoff h, the equilibrium state at inverse temperature β can be constructed as in (11). Using the spatial translation invariance of the interacting hamiltonian and denoting by τ t,x the * -automorphisms realizing a spacetime translation of step (t, x) we have where ω β T denotes the truncated n−point function of the state ω β,ψ . Hence, in order to discuss the limit h → 1 we need to control the decay for large spatial directions of the truncated n point functions. We have actually the following theorem = ω T (A 0 ; τ iu 1 ,x 1 (A 1 ); . . . ; τ iun,x 1 (A n )).

It exists a constant C such that
for r > 4cR, uniformly in u for 0 < u 1 < · · · < u n < β with β − u n ≥ β n+1 . Proof. Thanks to the decay property for large spatial separations of the locally smeared twopoint functions given in Proposition C.2, the proof of this theorem can be done in a similar way as the proof of Theorem 3 of [FL14]. We recall here the main steps of that proof and we adapt them to the case studied here.
The truncated n−point functions can be written as a sum over all possible connected graphs joining n points. We shall denote the set of connected graphs, without tadpoles, with n + 1 vertices V = {0, . . . , n} as G c n+1 . Furthermore, for any G ∈ G n+1 , E(G) denotes the set of edges of G. For any l ∈ E(G), s(l) and r(l) denote the source and the range of l. A graph G is considered to be in G n+1 only if, for every l, s(l) < r(l). Finally l ij (G) is the number of lines connecting i, j ∈ V . With these definitions where sym(G) = i<j l ij (G)! is a numerical factor and Furthermore, with the integral kernel K(x − y) = ω β,ψ (ψ(x)ψ(y)), given in terms of the thermal two-point function of the background theory (27). Furthermore, ψ j = ψ j 1 + iψ j 2 is the field configuration in the j−th factor of the tensor product and the functional derivative δ δψ j acts on the j−th factor of the tensor product. We have that where U = (u 0 , . . . , u n ), X = (x 0 , . . . , x n ) with u 0 = 0 and x 0 = 0, while P = (p 1 , . . . , p |E(G)| ) and where λ + and λ − are the positive and negative frequency part hence, separating the positive and negative contributions in F G we get where the sum is taken over all possible partitions of E(G) in up to two sets {E + (G), E − (G)} ∈ P 2 (E(G)). We proceed now splitting again these contributions over the two possible frequencies ω ± . Hence, denoting by The functionΨ is an entirely analytic function which grows at most polynomially in every direction, we might thus integrate over all possible p 0 to get Since, by hypothesis, and r(l) > s(l), we have that e −(β−(u r(l) −u s(l) ))ω σ ≤ e − n n+1 βω σ .
is rapidly decreasing, in every direction, because, F G is a microcausal functional and Φ(P) is the restriction on a particular subdomain ofΨ(−P, P ) which is an entire analytic function which grows at most polynomially. Hence, the negative frequencies are exponentially suppressed and if directions containing only positive frequencies are considered, they are also rapidly decreasing by Proposition C.1. The integral over P can now be taken and we may apply Proposition C.2 to estimate the decay of the result of that integral, actually, we obtain where the constant c does not depend on u i . In the last inequality we used the fact that G is a connected graph and thus every x i can be reached from the origin (x 0 = 0), hence thus concluding the proof.
To apply the results of Theorem 3.1 we observe that, if R is sufficiently large H I (0) ∈ A I (B R ), furthermore, the form of the integration domain of the u variables as given in (37) is such that Using the KMS condition we might restrict attention to the case where, β − u n ≥ β n+1 . In fact, if this is not the case, it must exist an m for which u m − u m−1 ≥ β n+1 . Actually, for A i ∈ A I (O) by the KMS condition we have that hence we might now consider F (v 1 , y 1 ; . . . ; v n , y n ) . = ω β T (K; τ iv 1 ,y 1 (K); . . . ; τ iv n−m ,y n−m (K); τ iv n−m+1 ,y n−m+1 (A 0 ); τ iv n−m+2 ,y n−m+2 (K); . . . τ ivn,yv (K)) at the place of F . In fact the previous equality obtained with the KMS condition together with translation invariance of the state implies that The arguments of the function F have the desired property, actually β − v n = u m − u m−1 ≥ β/(n + 1). We might thus use F at the place of F . Since the integration over the u variables is over a compact set and since the points where u i = u j for some i = j forms a zero measure set. Hence, the Theorem 3.1 implies that the integral over x i can be taken for all i to conclude the proof.

Spontaneous symmetry breaking and the Goldstone Theorem
The model we are considering possesses an internal U (1) symmetry. Actually, the Lagrangian is invariant under transformations where θ ∈ [0, 2π]. However, in the state which describes the condensate, this symmetry is spontaneously broken. From the Goldstone theorem we expect that a massless (gapless) mode is present in the model. This observation is in contrast with the analysis discussed in the previous section. Actually there, all the fields in the linearized theory have considered to be massive. Notice that, once the background φ is fixed in the decomposition ϕ = e −iµx 0 (φ + ψ), the Lagrangian for the linearized theory, is not invariant under U (1) transformations hence, the fact that both linearized fields ψ i are massive, is not in contrast with the Goldstone theorem. Furthermore, if the Goldstone theorem holds for the full theory, this would imply that at least one gapless mode should exists if the full perturbation series is considered.
In the case of thermal theories, the proof of Goldstone theorem is not completely straightforward as for theories at zero temperature because the original proof makes use of Lorentz invariance [GSW62] (see also the work of Jona-Lasinio using effective action methods in [Jo64]). The equilibrium states are however not Lorentz invariant because of the presence of a preferred time direction in the KMS condition. Furthermore, even if a gapless mode exists the particle content of the gapless mode is not immediately evident, as discussed by Bros and Buchholz in [BB98].
The presence of Goldstone modes at finite temperature, has been discussed in [Ko87] using effective action methods. Based on the analysis of Swieca [Sw67], a proof of the Goldstone theorem without using Lorentz invariance has been given by Morchio and Strocchi in [MS87], see also the book [St08] for the application of similar ideas for the analysis of the case of finite temperature. Furthermore, the analysis of the slow decay of large spatially separated correlation functions in the presence of spontaneous symmetry breaking is discussed in [JW11]. However, when a non trivial background is present as for the case of Bose Einstein condensation, we don't expect that the presence of a gapless mode is directly related to the clustering properties of the correlation functions for large spatial separation. As an example, consider the two-point function of the state ω ± β,c discussed in section 3.1 in the limit of vanishing temperature namely β → ∞. The obtained state is the composition of the massive vacuum ω 0 with the map γ ± c given in (17). Even if one of the modes in the two-point function of ω ± β,c is gapless, the clustering properties of ω ± β,c are equivalent to the one of the vacuum because γ ± c does not change the localization of the observables.
The mentioned proofs cannot be directly applied for perturbatively constructed theories, for this reason in the next section we shall give a proof of the validity of Goldstone theorem which can hold in our setting. We shall actually follow Swieca's proof without making use of Lorentz invariance.

Proof of the Goldstone theorem
Here we would like to give a proof of the validity of Goldstone theorem at finite temperature in the presence of a condensate. For this purpose, we observe that the U (1) invariance of the Lagrangian density L for ϕ = ϕ 1 + iϕ 2 is such that where R is a rotation of an angle θ. Its infinitesimal version is where t is the antisymmetric metric Associated to the symmetry which is spontaneoulsy broken in the state ω there is a current J which is conserved. This current is defined as By Noether theorem the action possesses the desired U (1)−symmetry if and only if the current J is conserved, namely if ∇ µ J µ = 0. Following [KRS66] we can now introduce a regularized charge operator associated to the current density J 0 introduced above. Let f ∈ C ∞ 0 (R) be a time cutoff with suppf ∈ (− , ), f ≥ 0 and f 1 = 1. Furthermore, g ∈ C ∞ 0 (R 3 ) is a space cutoff, g(x) = 1 for x < 1. The regularized charge operator associated to J can be seen as the large R limit of The charge operator can be used to implement the infinitesimal U (1) transformation of the field 4 lim R→∞ [Q R , ϕ m (t, y)] = t mn ϕ n (t, y).
Hence, in a state where the symmetry is spontaneously broken, namely when ω(ϕ n (0)) = φ n = 0 for some n, lim Notice that, in view of the support properties of f , g, and of the conservation of the current J, for R sufficiently large we have that ω([Q R , ϕ n (0)]) is constant. Hence, the limit R → ∞ can be safely taken and the final result does not depend on the particular form of g.
We are now ready to state the Goldstone theorem in the following form Theorem 4.1. (Goldstone) Consider a complex scalar quantum field ϕ whose Lagrangian density L possesses an U (1)−symmetry generated by the current J given in (40). Consider the distribution In perturbation theory equation (42) can be proved starting from the master ward identity, using the current conservation and the causal properties of the commutator. For further details in the case of QED we refer to [DF98].
If the symmetry generated by J is spontaneously broken in the state ω, namely ω(ϕ n (0)) = φ n = 0 for some n, then in the spectrum of G n there is a zero frequency (gapless) contribution at vanishing momentum, namelyĜ n , the Fourier transform of G n , is such that in the sense of distribution.
Proof. The support properties of the distribution F i n (x) , ϕ n (0)]) and invariance under spatial rotations imply that it exists a distribution F n such that where f and g are chosen as in the definition of regularized charge Q R given in (41). Causality implies that for R sufficiently large, the left hand side does not depend on R. Hence dp 0f (p 0 )p i F n (p 0 , |p|) must be bounded near p = 0 and as a consequence of this fact it holds that Hence, current conservation furnishes a condition for G n , namely Another constraint on the form of G n is given by (43), namely Both conditions imply that (44) holds in the sense of distributions.
As discussed in [BB98], no direct particle interpretation can be inferred from (44). Actually, the singularity inĜ n can proved to exists only at p = 0 and not on the whole null cone. For a particle interpretation of this fact we refer to the paper [BB98].

Analysis of the validity of Goldstone theorem in perturbation theory
We have seen that Goldstone theorem in the form stated in Theorem 4.1 holds for a quantum scalar field theory if the corresponding Lagrangian density L is invariant under the U (1) transformations given in (39), and if this symmetry is spontanously broken in a state ω. Hence, in a generic quantum field theory, we can apply Goldstone theorem if the current J given in (40) is conserved and if ω(ϕ n (0)) = φ n = 0 for some n.
We now check that at linear order some of the desired hypotheses are not satisfied. In particular, we see that in the linearized theory the internal U (1) symmetry is explicitly broken. Actually, the current J defined as in (40), whose time component has the explicit form (31), is not conserved. To check this fact notice that where we have used the equation of motion given at linear order (22). If now µ 2 > m 2 we have that φ = 0 because λφ 2 = (µ 2 − m 2 ). Furthermore, M 2 1 = (m 2 − µ 2 ) + 3λφ 2 , and M 2 2 = (m 2 − µ 2 ) + λφ 2 and thus even if M 2 = 0, we have that M 2 1 − M 2 2 = 0 and hence ∂ µ J µ = 0.
We pass now to analyse the interacting case. We observe that, if the full interacting equation of motion is used in evaluating ∂ µ J µ , namely taking into account L 3 and L 4 , we have that equation (45) needs to be changed as Now both M 2 = 0 and M 2 1 − M 2 2 − 2λφ 2 = 0 and hence the symmetry is not explicitly broken in the full classical theory. Notice that this analysis does not depend on the splitting between free and interacting Lagrangian, hence the eventual thermal mass contributions do not alter this analysis.
We now check that conservation of the current J holds in the case of interacting quantum field theory treated with perturbation methods. The very same analysis holds up to the quotient with respect to the free equation of motion. Actually, the Schwinger-Dyson equation implies that where S is the action constructed with the full Lagrangian density L and S 0 is the action of the linearized theory, constructed with the lagrangian density L 2 . If expectation values on a quantum state are considered we have We now analyze in the same spirit the conservation of the current density J. To this end, we observe that where the last equality holds because of the properties of the time ordered product and S = S. We stress that the divergences present in ϕ(x)· T δS/δϕ(x) and proportional to lim y→x ∆ F (x, y) δ 2 S δϕ(y)δϕ (x) are cancelled by the subtraction of ϕ · T δS/δϕ. To proceed we need the Master Ward Identity which can be shown to hold for this theory without anomalies. Master Ward identities have been discussed in [DF03,DF04,Ho07]. For our purposes we start with equation (70) of [FK13].
Rewriting it in our context, we obtain the desired equation 5 Notice that, the divergences on the right hand side of the previous equality due to the pointwise multiplication with the quantum product vanishes because δS 0 /δϕ is in the ideal of the linear equation of motion. Hence, in any quantum state We conclude that the current J constructed with the interaction quantum scalar field is conserved in the sense of perturbation theory up to the ideal describing the equation of motion. Finally we observe that even if this last observation has been given in terms of the field ϕ, since relations (46) and (47) are algebraic relations, the very same analysis holds for the fluctuations ψ 1 , ψ 2 . We also observe that on the state we are considering ω(ϕ 1 ) = 0 in the sense of perturbation theory. Actually φ = 1 √ λ µ 2 − m 2 and thus L 3 is of order √ λ while L 4 is of order λ (m 2 v needs to be chosen smaller than λc). This implies that ω(R V ψ 1 ) is at least of order O(1) in λ and hence it cannot totally cancel/compensate φ.
We thus have that the hypotheses of Theorem 4.1 holds in the sense of perturbation theory. Furthermore, we also notice that all the identities in the proof holds in the sense of perturbation theory. Hence we conclude that the thesis of that theorem holds in the sense of perturbation theory. 5 A direct proof of (47) can be obtained studying where the limit is taken in a direction where y is always in the future of x. The first equality is a consequence of the causal factorisation property and the last equality is the Schwinger-Dyson equation

Acknowledgments
NP thanks the ITP of the University of Leipzig for the kind hospitality during the preparation of part of this work and DAAD for supporting that visit with the program "Research Stays for Academics 2017".

A The principle of perturbative agreement
Thanks to the principle of perturbative agreement [HW05,DHP17] the theories obtained with the two different splittings given in (12) and (13) of the Lagrangian density are equivalent. In particular, the principle of perturbative agreement holds because it exists a map γ which intertwines the time ordered products constructed in theories with different masses m and m which will be later set to 0. We shall here denote by T m the time ordering operator of a massive theory − m 2 . Hence, where here F mloc is the commutative algebra of multilocal functionals (with respect to the pointwise (classical) product). The non trivial nature of γ can be seen from the fact that the time ordered products T m and T m are different. At the same time, the perturbative agreement requires that the common part of L I and L I , or more generally all Wick powers, are left invariant up to ordinary renormalization freedom by the map γ. We shall here see that this is the case provided the time ordering operator T m are constructed with the preferred symmetrized Hadamard singularity H m of − m 2 as in (2) with W = 0 for some fixed lengthscale ξ, namely T m (F ) =: F : Hm . For completeness we discuss this invariance under γ in the simple case of the Wick square We want to compare T m Φ 2 (f ) with T m Φ 2 (f ). To this end we observe that and we study Suppose that, on regular functionals, T m and T m are constructed starting with the Feynman propagators ∆ F,m (x, y) and ∆ F,m (x, y), hence we have Furthermore, the principle of perturbative agreement (see e.g. [DHP17] Theorem 3.2 where γ is denoted by β) implies that γϕ = ϕ and that γT m (ϕ(x)ϕ(y)) = T m (ϕ(x)ϕ(y)) and thus It remains to compare ∆ F,m with ∆ F,m . Considering the form of the Hadamard singularities, we observe that where V m is the standard Hadamard V coefficient, σ is the regularized squared geodesic distance, µ is a length scale and W is a smooth reminder. In the coinciding point limit (V m − V m ) is proportional to δm 2 = m 2 − m 2 , hence the limits present in the previous equation are divergent. However, the kind of divergences present in that coinciding point limit are not different than the standard divergences present in R Q (ϕ 2 ) with if the time ordered product of T m (Qϕ 2 ) is constructed without making use of a correct renormalization prescription. We may thus avoid (renormalize) them in a similar way as the divergences present in the naive construction of T m are resolved. To this end we notice that, if we consider the Feynman operator as operator on functions, we have that We observe, by power counting, that the divergent contributions in the sum are the contributions n = 0 and n = 1. The contributions n = 0 is removed by the subtraction of ∆ F,m while the divergences present in the contribution n = 1 are similar to the divergences present in T m (QQ) and thus they can be treated with renormalization theory. Actually, at order 1 in δm 2 the remaining contribution ∆Φ 2 (1) in the difference γT m Φ 2 − T m Φ 2 is thus We regularize the product of Feynman propagators in the following way where a is a parameter which takes into account the known renormalization freedom. If we fix it to 0 we have that where in the last step we used the fact thatĝ tends to δ(p) in the limit g → 1. Higher order contributions can be computed directly and they give a finite result.
Instead of performing these explicit computations, fur our purposes, we need to understand the origin of these contributions with the following observation. Actually, the essential steps in the previous computation are the following: we have computed the expansion in powers of δm 2 of ∆ F,m (x − y), we have removed the contributions of order 0 and 1 and we have eventually taken the coinciding point limits. In formulae We prove now that if one starts with the time ordered propagators constructed with the preferred Hadamard functions H instead of the Minkowski vacuum ∆Φ 2 vanishes, we have actually that and expanding V in powers of σ we have Hence, we conclude that The same argument essentially holds also on any curved background, so in the case of mass renormalization we see that there is a choice of renormalization freedom (a = 0 in (48)) such that where now the Wick powers T m Φ 2 =: Φ 2 : Hm are constructed (regularized) with respect to the distinguished Hadamard singularity H m and not with respect to the symmetric part the twopoint function of a state. Finally, we observe that, the same results holds also when Wick monomial of higher order are considered, namely γT m Φ n = T m Φ n .

B First Hadamard Coefficients for D
We compute in this section the first Hadamard coefficients for the operator D given in (24). We recall that in this case the Hadamard singularity has the following structure similar to (2) with vanishing W , where now U, V, V i are 2 × 2 matrices and where σ is again one half of the regularized geodesic distance. In the case of Minkowski σ(x, y) = 1 2 (x − y) µ (x − y) µ . The requirement that DH is smooth implies the following transport equations The first two equations give considering integrals along the geodesic γ joining x, y and indicating by r its affine parameter (γ(0) = x, γ(1) = y), that equation gives Integrating along γ we get The first transport equation can be solved once the initial condition U (x, x) = I is fixed (as we obtained in the previously computed vacuum state). To find its solution, due to translation invariance, it is enough to study a solution which satisfy the desired initial condition is hence, since U −1 = cos(µx 0 )I + iσ 2 sin(µx 0 ) We can now expand V as n≥0 V n σ n . The equation DV = 0 and the knowledge of V 0 , permits to compute V n recursively. We are in particular interested in [V 1 ](x) = V (x, x) because this coefficient is proportional to the trace anomaly of the stress tensor of the linearized theory [Wa78, Mo03,HW05] and in particular enters in the expressions We observe that furthermore we can expand V 0 = U X for some X hence we have Summarizing this analysis We observe that [V 1 ] is diagonal and constant, hence, we see that this anomaly is not visible in the conservation of the charge J µ given in (45).

C Some technical propositions
We report here a technical proposition, similar to Proposition 9 in [FL14], and adapted to the case studied here.
Proof. Since A i ∈ A I , is a microcausal functional we have that and the same holds with W − at the place of W + , where p l = 0 and hence, since for all l p l ∈ V + , for every l such that s(l) = 0 p l = 0. Furthermore, If we have already proved that p l = 0 for every s(l) = j < i, we can prove it also for s(l) = i. Actually, in that case we already know that for every l such that r(l) = i, s(l) < i and thus p l = 0. Hence, the direction P we are analysing can be of non rapid decrease only if s(l)=i p l = 0 which implies that for every l such that s(l) = i, p l = 0. We have thus proved that all these directions are of rapid decrease because the direction P can be of non rapid decrease only if P = {0} and the zero section does not intersect W F (Ψ).
In the following proposition we prove the exponential decay of the two-point function of the KMS states for the linearized theory studied in section 3.3. This proposition is used in the proof of the clustering estimate necessary for the adiabatic limit.
whereD is as in (26), see also (24), and ω ± as in (25). Assume M − > 0, it holds that Hence, we analyze separately the following functions I is then a linear combination of I i with constant coefficients. We prove with some details the decay of I 2 σ for σ ∈ {+, −}, analogous results holds also for the other components. Without loosing generality, let us assume that x = rn with n = (1, 0, 0). We shall discuss the decay for large r. Notice that We evaluate the integral in p 1 with the help of complex analysis considering p 1 = z = x+iy with x, y ∈ R a complex variable. We shall furthermore, consider a particular contour of integration in the upper half plane. The integral over one branch of the contour we shall chose corresponds to I 2 σ , furthermore the contour will be extended to infinity and chosen in such a way to avoid the poles in 1/(ω 2 + − ω 2 − )2ω σ and the brunch cuts present in ω σ . Actually, we observe that since f is a compactly supported distribution, and its support is contained in the disc centered in 0 and of radius R, by the Paley Wiener theorem, its Fourier transform is an entire analytic function. Furthermore, it grows at most polinomially in every real direction and exponentially in complex directions. Hence, it exist two constants C > 0, C > 0 and a N > 0 such that hence, if it is composed with ω σ and if it is seen as a function of p 1 , it is analytic everywhere up the branch cuts which are present in the principal squares of ω ± = w 2 + 2µ 2 ± 4µ 4 + 4µ 2 w 2 + (δM 2 ) 2 .
This implies that the integrand of I 2 σ , seen as a function of p 1 is analytic every where up to the poles of 1/(ω 2 + − ω 2 − )2ω σ and the branch cuts mentioned above. To describe the contour of integration we analyze the location of the branch cuts and the poles.
The p 1 integration can evaluated by choosing a contour in the upper half plane. The contour can the taken in such a way that both the the poles and the branch cuts lie outside it. Furthermore, the non trivial part of the contour is chosen to be at imaginary part larger than |p ⊥ | 2 + M 2 − δM 2 . The contour ξ is formed by the following paths ξ 1 , ξ 2 , ξ 3 , ξ 4 , ξ 5 where ξ 1 is the real line, ξ 2 is part of the semicircle centered in the origin of the complex plane and with radius B tending to infinity. Furthermore the points with |Rez| < 2µ are removed from the circle.