Pad\'e approximants and analytic continuation of Euclidean Phi-derivable approximations

TL;DR

This paper explores the effectiveness of Padé approximants for analytically continuing Euclidean propagator data to extract spectral functions and physical information, tested on benchmarks and applied to the O(4) model.

Contribution

It demonstrates the use of Padé approximants for analytic continuation in Phi-derivable approximations and compares zero-momentum and pole masses in a nonperturbative context.

Findings

Padé approximants successfully access spectral functions from Euclidean data.

The method reveals differences between zero-momentum and pole masses.

Application to the O(4) model shows potential for nonperturbative analysis.

Abstract

We investigate the Pad\'e approximation method for the analytic continuation of numerical data and its ability to access, from the Euclidean propagator, both the spectral function and part of the physical information hidden in the second Riemann sheet. We test this method using various benchmarks at zero temperature: a simple perturbative approximation as well as the two-loop Phi-derivable approximation. The analytic continuation method is then applied to Euclidean data previously obtained in the O(4) symmetric model (within a given renormalization scheme) to assess the difference between zero-momentum and pole masses, which is in general a difficult question to answer within nonperturbative approaches such as the Phi-derivable expansion scheme.

Figures

Equations

$${\int }_Q\equiv {\int }_{-\infty }^{\infty }\frac{d{q}_0}{2\pi }\int \frac{d{q}^3}{(2\pi {)}^3}.$$

$${\bar{M}}^2(Q)=m_0^2+\frac{{\lambda }_2}{2}{\varphi }^2+\frac{{\lambda }_0}{2}\mathcal{T}[\bar{G}]+\frac{{\lambda }^2}{2}{\varphi }^2\mathcal{I}[\bar{G}](Q),$$

$$\mathcal{T}[G]$$

$$\mathcal{I}[G_1,G_2](K)$$

$$0=\bar{\varphi }\left(m_2^2+\frac{{\lambda }_4}{6}{\bar{\varphi }}^2+\frac{{\lambda }_2}{2}\mathcal{T}[\bar{G}]+\frac{{\lambda }^2}{6}\mathcal{S}[\bar{G}]\right),$$

$$\mathcal{S}[G_1,G_2,G_3]$$

$${\bar{M}}^2(Q)$$

$$0$$

$${\mathcal{T}}_{\mathrm{F}}[G]$$

$${\mathcal{I}}_{\mathrm{F}}[G](Q)$$

$${\mathcal{S}}_{\mathrm{F}}[G]$$

$$\delta G(Q)$$

$$G_{\mathrm{r}}(Q)$$

$$G(Q)={\int }_0^{\infty }\frac{ds}{2\pi }\frac{i\rho (s)}{Q^2-s+i\epsilon },$$

$${\int }_0^{\infty }\frac{ds}{2\pi }\rho (s)=1.$$

$$\Im {\mathcal{I}}_{\mathrm{F}}[G](Q)={\int }_0^{\infty }\frac{d{s}_1}{2\pi }{\int }_0^{\infty }\frac{d{s}_2}{2\pi }\rho (s_1)\rho (s_2)\Im {\mathcal{I}}_0[G_1,G_2](Q),$$

$$\Im{\cal I}_{0}[G_{1},G_{2}](Q)=-\frac{\sqrt{\lambda(Q^{2},s_{1},s_{2})}}{16\pi Q^{2}}\Theta\big{(}Q^{2}-(\sqrt{s_{1}}+\sqrt{s_{2}})^{2}\big{)},$$

$$\displaystyle\Re{\cal I}_{\rm F}[G]\left(\sqrt{s}\right)={\cal I}_{\rm F}[G]+\frac{s}{\pi}{\cal P}\int_{0}^{\infty}ds^{\prime}\frac{\Im{\cal I}[G]\big{(}\sqrt{s^{\prime}}\big{)}}{s^{\prime}(s^{\prime}-s)}\,,$$

$${\mathcal{I}}_{\mathrm{F}}[G]$$

$${\mathcal{I}}_{0,\mathrm{F}}[G_1,G_2]=\frac{(s_1-s_2{)}^{-1}}{16{\pi }^2}\left(s_1\ln \frac{s_1}{e{M}_0^2}-s_2\ln \frac{s_2}{e{M}_0^2}\right).$$

$$\rho (s)=2\pi Z\delta (s-{\bar{M}}_{\mathrm{p}}^2)+\sigma (s),$$

$${\bar{M}}_{\mathrm{p}}^2={\bar{M}}^2({\bar{M}}_{\mathrm{p}}^2).$$

$$\Re {\mathcal{I}}_{\mathrm{F}}[G]\left(\sqrt{s}\right)={\mathcal{I}}_{\mathrm{F}}[G]+\frac{s}{\pi }\mathcal{P}{\int }_{s_{\mathrm{th}}}^{\infty }d{s}^{\prime }\frac{\Im \mathcal{I}[G]\left(\sqrt{s^{\prime }}\right)}{s^{\prime }(s^{\prime }-s)},$$

$$G_{\mathrm{E}}(Q_{\mathrm{E}})={\int }_0^{\infty }\frac{ds}{2\pi }\frac{\rho (s)}{Q_{\mathrm{E}}^2+s}.$$

$$m^2+\frac{{\lambda }_{\mathrm{c}}}{2}{\mathcal{T}}_{\mathrm{F}}[{\bar{G}}_{\varphi =0}]+\frac{{\lambda }_{\mathrm{c}}^2}{6}{\mathcal{S}}_{\mathrm{F}}[{\bar{G}}_{\varphi =0}]=0.$$

$$C_N(z)=\frac{a_1}{1+}\frac{a_2(z-z_i)}{1+}\cdots \frac{a_N(z-z_{N-1})}{1},$$

$$g_1(z_i)$$

$$g_p(z)$$

$${\mathcal{B}}_{0,\mathrm{F}}[G_M](Q_{\mathrm{E}})=\frac{1}{16{\pi }^2}\left(2-\log \frac{M^2}{M_0^2}+\sqrt{1+\frac{4M^2}{Q_{\mathrm{E}}^2}}\log \frac{\sqrt{1+\frac{4M^2}{Q_{\mathrm{E}}^2}}-1}{\sqrt{1+\frac{4M^2}{Q_{\mathrm{E}}^2}}+1}\right).$$

$$\displaystyle{\cal Q}=\int_{\sqrt{s}_{\rm min}}^{\sqrt{s}_{\rm max}}dq\big{(}F(q)-{\cal F}(q)\big{)}^{2},$$

$${\mathcal{G}}_{\mathrm{E}}^{-1}(x,y)={\frac{{\delta }^2\Gamma [\varphi ,G]}{\delta \varphi (x)\delta \varphi (y)}|}_{G={\bar{G}}_{\varphi }}-\frac{1}{2}{\int }_{z_1,z_2,z_3,z_4}\Lambda (x;z_1,z_2){\bar{G}}_{\mathrm{E}}(z_1,z_3){\bar{G}}_{\mathrm{E}}(z_4,z_2)V(z_3,z_4;y),$$

$${\frac{{\delta }^2\Gamma [\varphi ,G]}{\delta \varphi (x)\delta \varphi (y)}|}_{G={\bar{G}}_{\varphi }}={\mathcal{G}}_{\mathrm{E},0}^{-1}(x-y)+\frac{1}{2}\left[{\lambda }_4{\varphi }^2+{\lambda }_2{\bar{G}}_{\mathrm{E}}(0)\right]\delta (x-y)-\frac{{\lambda }^2}{6}{\bar{G}}_{\mathrm{E}}^3(x-y),$$

$${\hat{M}}_{\mathrm{E}}^2(Q_{\mathrm{E}})=\delta Z_2{Q}_{\mathrm{E}}^2+m_2^2+\frac{{\lambda }_4}{2}{\varphi }^2+\frac{{\lambda }_2}{2}{\mathcal{T}}_{\mathrm{E}}[{\bar{G}}_{\mathrm{E}}]-\frac{{\lambda }^2}{2}{\varphi }^2\mathcal{B}[{\bar{G}}_{\mathrm{E}}](Q_{\mathrm{E}})-\frac{{\lambda }^2}{6}{\mathcal{S}}_{\mathrm{E}}[{\bar{G}}_{\mathrm{E}}](Q_{\mathrm{E}}).$$

Full text

I Introduction

In recent years, much progress was achieved in accessing the properties of interacting quantum fields in equilibrium. Part of the success stems from the development of continuum nonperturbative tools Braun:2007bx ; Braun:2014ata ; Contant:2017gtz ; Reinhardt:2017pyr that allow to implement some of the resummations of perturbative diagrams that are required in this context. Most of these approaches are formulated in Euclidean space (or imaginary-time formalism) because, on general grounds, the equations are easier to solve than the corresponding ones in Minkowski space (or real-time formalism). For instance, to date, there is no complete Minkowski space study of the temperature driven symmetry restoration in the four dimensional $O(N)$ model which could compete in precision with its Euclidean counterpart. Yet, certain quantities of interest, such as transport coefficients, pole masses and decay rates, require the use of the real-time formulation, and, therefore, a constant effort is put into extending the previous methods (not only formally but also at a computational level) to Minkowski space. In particular, the self-consistent scalar propagator equation was solved in Minkowski space at zero and finite temperatures, both in the symmetric and broken phases (see e.g. Sauli:2001mb ; VanHees:2001pf ; Cooper:2004rs ; Sauli:2004bx ; Arrizabalaga:2007zz ; Jakovac:2006dj ; Jakovac:2006gi ). These spectral function based numerical solutions were obtained in various nonperturbative approximations. More recently, the spectral function was accessed using an analytic continuation of the Functional Renormalization Group (FRG) equations Kamikado:2013sia ; Strodthoff:2016pxx .

Another possible route towards real-time quantities is not to employ the above approaches directly in Minkowski space, but rather, to develop methods that allow for the analytic continuation of Euclidean data. A standard method used to extract the spectral function from Euclidean data emploies Bayesian inference to invert the integral relation between Euclidean correlation functions and spectral functions. In order to overcome some difficulties of the Maximal Entropy Method (MEM) Jarrell:1996rrw , a novel Bayesian method was developed recently in Ref. Burnier:2013nla and applied for example to extract spectral properties red from bottomonium correlators in Ref. Kim:2014iga and solutions of the quark Dyson-Schwinger equation obtained in Landau gauge Fischer:2017kbq . In the context of the $\Phi$-derivable expansion scheme (also referred to as the two-particle irreducible (2PI) formalism), knowledge of the external propagator, which fully displays the (linear) symmetries of the theory, requires the resolution of a Bethe-Salpeter equation. Since the latter equation is more easily solved in Euclidean space (see for instance Carrington:2013jta ), analytic continuation techniques are of great help in this case.

In this paper, we would like to investigate the analytic continuation method of Ref. VS_Pade which constructs multipoint Padé approximants in the form of finite continued fractions. We test this method using various benchmarks. First, we consider the two-loop $\Phi$-derivable approximation for a one-component scalar field at zero temperature for which Euclidean data are already available Fejos:2011zq and Minkowski data are relatively easy to generate using dispersion techniques. Using this benchmark, we test to what extent the pole mass and more generally the spectral function determined directly in Minkowski space can be reproduced from the Euclidean solution. Our test shows also that we have some control on the choice of the appropriate analytic continuation (among the infinitely many possible choices). As a second test, we use simple one-loop perturbative formulas and investigate the ability of the continuation method to extract (physical) information hidden in the second Riemann sheet.

We then apply the Padé analytic continuation method on the previously obtained Euclidean space solution of Refs. Marko:2013lxa ; Marko:2015gpa to assess how much pole masses differ from zero-momentum masses. A similar question was recently addressed in Ref. Helmboldt:2014iya in the context of the quark-meson model using the FRG approach. In our case, such a comparison is motivated by the fact that, in our previous Euclidean studies of the $O(4)$ model, the parametrization was done at $T=0$ based, for simplicity, on the pion and sigma zero-momentum masses of the external propagator which are directly accessible from the effective potential in the imaginary-time formalism. Since a more correct parametrization should involve the corresponding pole masses, we would like to quantify how much the parametrization was deformed by the use of the zero-momentum masses. We mention also that, recently, the Padé analytic continuation method was applied in Ref. Pilaftsis:2013xna to compute the pole mass of the internal propagator of the symmetry improved 2PI (SI2PI) formalism. In Ref. Marko:2016wtw we argued that, at least at two-loop order, the SI2PI formalism contains an inherent untamed infrared sensitivity in the broken phase, according to which the Euclidean propagator is not defined above some value of the volume of the system.111This infrared sensitivity can be avoided for some specific (sharp) UV regulators. However, the specificity of these UV regulators makes this removal of the IR sensitivity artificial in a sense, see Ref. Marko:2016wtw for more details. Here we use the Padé analytic continuation method to investigate to which extent this IR sensitivity affects the pole mass of the corresponding Minkowski propagator.

The organization of the paper is as follows. In Sec. II, using dispersion techniques, we generate Minkowski data within the two-loop $\Phi$-derivable approximation for a $\varphi^{4}$ model at zero temperature, with or without spontaneous symmetry breaking. In Sec. III, we present the multipoint Padé continuation method and perform various tests using both the Minkowski space solution of Sec. II and simple perturbative formulas. The continuation method is then used in Sec. IV to assess the difference between pole and zero-momentum masses in the context of the $O(4)$ model, as well as the IR sensitivity of the pole mass in the two-loop SI2PI approximation. The conclusions of our study are presented in Sec. V and some technicalities can be found in the Appendices.

II Minkowskian two-loop $\Phi$-derivable approximation in the $\varphi^{4}$ model

In this section, we consider a relatively simple situation where the $\Phi$-derivable equations can be solved directly and accurately in Minkowski space. The obtained results shall then be used as a benchmark for the Padé analytic continuation method to be presented in the next section. In what follows $Q=(q_{0},\vec{q})$ denotes a four-momentum in Minkowski space and we introduce the notation

[TABLE]

II.1 Generalities

The 2PI formalism gives typically access to a self-consistent equation for the two-point function of a given model, from which various observables can be determined, in principle. In contrast to the infinite tower of Dyson-Schwinger equations, the equation for the two-point function in the 2PI formalism is not coupled to higher $n$-point functions. The prize to pay is, however, that the equation contains infinitely many terms coming from the loop expansion of the 2PI effective action and then it needs to be truncated for any practical purpose. In the case of the $\varphi^{4}$ model, the two-loop truncation of the 2PI effective action leads to the following gap equation for the momentum dependent mass function $\bar{M}(Q)$ of the self-consistent propagator $\bar{G}(Q)\equiv i/(Q^{2}-\bar{M}^{2}(Q)+i\epsilon)$ in the presence of an arbitrary field expectation value $\phi$:

[TABLE]

where the tadpole and bubble integrals are defined respectively as

[TABLE]

with the shorthand notation ${\cal I}[G](K)\equiv{\cal I}[G,G](K)$. The physical value of the field $\phi$, denoted $\bar{\phi}$, is given in terms of the field equation. In the two-loop order truncation, it reads

[TABLE]

where the setting-sun integral at vanishing external four-momentum is defined as

[TABLE]

with a similar shorthand notation as above, ${\cal S}[G]\equiv{\cal S}[G,G,G]$. Let us mention that the field equation may have multiple solutions. The solution $\bar{\phi}$ we are interested in is the one that minimizes the effective potential. The latter can also be computed within the 2PI formalism, but we shall not recall its expression here.

As discussed in Ref. Berges:2005hc , a total of five counterterms are needed in order to absorb the divergences present in Eqs. (1) and (4). These counterterms can be expressed in terms of two renormalized quantities, a mass parameter $m^{2}$ and a coupling $\lambda,$ as well as a renormalization scale. Expanding the propagator around an auxiliary propagator $G_{0}(Q)=i/(Q^{2}-M_{0}^{2}+i\varepsilon),$ in which $M_{0}$ plays the role of the renormalization scale, the counterterms were explicitly determined in Ref. Patkos:2008ik and employed numerically in Ref. Fejos:2011zq . Similar counterterms, although determined at some fixed temperature in the symmetric phase, have been used in Arrizabalaga:2007zz ; Marko:2012wc ; Marko:2013lxa ; Marko:2015gpa . In the present study, we use the counterterms of Ref. Patkos:2008ik and, in Appendix A, we show the formal relation between the renormalization method of Ref. Patkos:2008ik and that of Ref. Berges:2005hc , which is formulated in terms of various two-point and four-point functions. The explicitly finite gap and field equations obtained after renormalization are

[TABLE]

The finite parts of the integrals, denoted by the subscript ‘F’, read

[TABLE]

where ${\cal I}[G_{0}]\equiv{\cal I}[G_{0}](Q=0)$ and, as a result of expanding $G$ around $G_{0}$, one has

[TABLE]

with $M^{2}_{\rm l}\equiv m^{2}+\frac{\lambda^{2}}{2}\left(\phi^{2}+{\cal T}_{\rm F}[G]\right)$. Each of the three terms in the r.h.s. of Eq. (7c) is UV finite.

II.2 Use of dispersion relations

To put the system (6) in a convenient form in view of its numerical resolution, we now make use of dispersion relations. We first focus on the momentum dependent bubble integral ${\cal I}_{\rm F}[G](K)$. A standard method for dealing with such an integral in Minkowski space is to use the spectral representation for the propagator and a subtracted dispersion relation (see e.g. Ref. Bjorken_Drell_RQF ) relating the real and imaginary parts of the integral. The spectral representation for the propagator is

[TABLE]

with the spectral function defined as the (real) function $\rho(Q^{2})=2\Im[iG(Q)].$ In the two-loop approximation considered here, the leading contribution to the propagator at large $Q$ is entirely given by the tree-level term, $G(Q)\sim i/Q^{2}$. It then follows from Eq. (9) that the spectral function obeys the following sum rule222This relation still holds beyond the present approximation for the spectral function associated with the bare propagator, in the presence of an ultraviolet regulator. It needs to be modified by a renormalization factor for the spectral function associated to the renormalized propagator.

[TABLE]

We note also that $\rho(s)\sim 1/s^{2}$ at large $s$ in the approximation at hand (if $\bar{\phi}\neq 0$).

Using Eq. (9) in a particular integral containing the self-consistent propagator is helpful because the integral over the momentum becomes perturbative and can be carried out using standard techniques, yielding an integral kernel in the rewritten expression of the original integral. In particular, we can use this rewriting in order to evaluate the imaginary part of the bubble integral ${\cal I}_{\rm F}[G](Q)$. Using Eq. (9) in the integral for ${\cal I}[G](Q)$, carrying out the momentum integration, and taking the imaginary part, while noting that $\Im{\cal I}_{\rm F}[G](Q)=\Im{\cal I}[G](Q)$, one obtains

[TABLE]

where $G_{j}(Q)=i/(Q^{2}-s_{j}+i\varepsilon),$ $j=1,2$. The UV finite kernel $\Im{\cal I}_{0}[G_{1},G_{2}](q)$ is the imaginary part of the perturbative bubble with mass squares $s_{1}$ and $s_{2}.$ It is given in terms of the Källén function $\lambda(x,y,z)=(x-y-z)^{2}-4yz$ as

[TABLE]

with $\Theta(x)$ the Heaviside step function.

We could obtain the real part of ${\cal I}_{\rm F}[G](Q)$ in the same fashion. It is however more efficient to make use of the dispersion relation that connects it to its imaginary part. Since the real part of the finite bubble integral grows logarithmically in the approximation at hand, we need a once-subtracted dispersion relation (see Appendix B for more details):

[TABLE]

where we have again used that $\Im{\cal I}_{\rm F}[G](Q)=\Im{\cal I}[G](Q)$ as well as $\Re{\cal I}_{\rm F}[G]={\cal I}_{\rm F}[G]$, with our notational convention ${\cal I}[G]={\cal I}[G](Q\equiv 0)$. Each term in the previous relation is finite. The subtracted piece ${\cal I}_{\rm F}[G]$ can be conveniently computed using Eq. (9) for $G$. We find

[TABLE]

where ${\cal I}_{0,\rm F}[G_{1},G_{2}]={\cal I}_{0}[G_{1},G_{2}]-{\cal I}[G_{0}]$ and we have used the sum rule (10) to bring ${\cal I}[G_{0}]$ under the integral. The difference ${\cal I}_{0,\rm F}[G_{1},G_{2}]$ of pertubative bubble integrals is UV finite. It can then be computed in any regularization scheme, with the result

[TABLE]

We now turn our attention to the principal value integral appearing in Eq. (13). Since this integral depends on which interval the imaginary part of the bubble has a support, we have to say a few words on the form of the spectral function. In the present section we restrict ourselves to spectral functions which have a singular part, corresponding to a real pole of the propagator. In this case one can write

[TABLE]

where the continuum part $\sigma(s)$, starting at the two-particle threshold, is given in Eq. (44) and the pole mass $\bar{M}_{\rm p}$ is defined as $\bar{G}^{-1}(Q=\bar{M}_{\rm p})=0$, or, equivalently, in terms of the gap mass $\bar{M}^{2}(Q)$, as

[TABLE]

Using that the imaginary part of the bubble integral, along with $\sigma,$ is non-zero only for $s>s_{\rm th}=4\bar{M}^{2}_{\rm p}$, one can write the real part of the finite bubble integral in the form

[TABLE]

with ${\cal I}_{\rm F}[G]$ given in (14).

Finally we discuss the evaluation of the finite tadpole and setting-sun integrals ${\cal T}_{\rm F}[G]$ and ${\cal S}_{\rm F}[G]$ defined in (7a) and (7c). The simplest way to compute them is by using the Euclidean propagator obtained through analytic continuation from (9) as

[TABLE]

For instance, the unsubtracted tadpole ${\cal T}[G]$ that enters ${\cal T}_{\rm F}[G]$ can be written, owing to a Wick rotation, as ${\cal T}[G]=\int_{Q_{\rm E}}G_{\rm E}(Q_{\rm E})\equiv{\cal T}_{\rm E}[G_{\rm E}]$. Similarly, one finds $\smash{{\cal S}[G]=-\int_{Q_{\rm E}}\int_{K_{\rm E}}G_{\rm E}(Q)G_{\rm E}(K)G_{\rm E}(Q+K)\equiv-{\cal S}_{\rm E}[G_{\rm E}]}$. With the Euclidean version of (8), ${\cal T}_{\rm F}[G]$ and ${\cal S}_{\rm F}[G]$ are then evaluated with an appropriate numerical cutoff.

II.3 Results

The iterative numerical algorithm for the resolution of the system (6) is described in detail in Appendix C. Here we concentrate on describing the results obtained in the broken symmetry phase of the model. We explore in Sec. II.3.2 the region of the parameter space ($m^{2},\lambda$) where such type of solution exists. All dimensionfull quantities are measured in units of $M_{0}$ throughout this paper, that is in the numerical code $M_{0}=1.$

II.3.1 Minkowskian solution

A typical broken phase solution of the gap and field equations (6) can be seen in Fig. 1: the left panel shows the real and imaginary part of the self-energy, while the right panel shows the pole and the continuum part of the corresponding spectral function. We compare these quantities to perturbative ones obtained using a tree-level spectral function $\rho(Q)=2\pi\delta(Q^{2}-M^{2})$ to evaluate the integrals, with $M^{2}\equiv\bar{M}^{2}(Q=0)$. We see that even for relatively large value of the coupling the 2PI self-energy is not much different from a perturbative one. This is due to the fact that in the two-loop approximation the momentum dependence of the self-energy is logarithmic and that the mass $M^{2}\equiv\bar{M}^{2}(Q=0))$ is then very close to $\bar{M}_{\rm p}^{2}.$ Note however that resummation of a perturbative series is needed in order to know the value of $M^{2}$, hence the comparison tells that the momentum dependence of the self-energy is similar to a perturbative one.333Another possibility is to parameterize the system directly in terms of $M^{2}$, with however the important subtlety that two physically distinct systems (one in the broken phase and one in the symmetric phase) can lead to the same parameters $(M^{2},\lambda)$, see the discussion in Appendix C.

II.3.2 Parametrization

In previous studies Marko:2012wc ; Reinosa:2011ut , we analyzed the parameter space of the model at $\smash{T=0}$ and how it was divided into regions corresponding to systems that displayed a symmetric phase and systems that displayed a broken phase. In those studies, the renormalization was carried out in the symmetric phase at a large enough fixed temperature $T_{\star}$ that played the role of the renormalization scale. It is interesting to see how this discussion appears in the present scheme where the renormalization is performed directly at $T=0$, where the parameters $m^{2}$ and $\lambda$ do not have the same meaning as in previous studies, and where the role of the renormalization scale is played by $M_{0}$ rather than $T_{\star}$.

Apart from the occurrence of the above mentioned two regions of the parameter space, the detailed investigation of Ref. Marko:2015gpa shows that, within a given 2PI approximation, one could have physical and unphysical branches of solutions to the gap equation (as $\phi$ is varied), which could merge for some values of the field. In this case, there is a region of $\phi$ where the gap equation admits no solution. This could result in the absence of solution to the coupled system of gap and field equations if the would-be $\bar{\phi}$ were to be engulfed by the region of $\phi$ over which the gap equation has no solution. We have shown in Ref. Marko:2015gpa that a localized approximation to the momentum dependent gap mass is useful to investigate the loss of solution. We now show that an unphysical solution to the full 2PI equations can indeed be found. However, for the sake of simplicity, we use the localized approximation introduced in Ref. Marko:2015gpa to find the region of parameters where a loss of solution happens.

Our analysis reveals the existence of three curves which, as shown in the left panel of Fig. 2, delimit different regions of the $m^{2},\,\lambda$ parameter space. The physically most relevant one, denoted $\lambda_{\rm c}(m^{2})$, separates the region where only a symmetric phase exists from the region where there is potentially a broken phase solution. This curve is a critical line in parameter space, in the sense that, for these parameters, the physical solution is at $\bar{\phi}=0$ with no $\bar{\phi}\neq 0$ solution and the curvature of the potential at $\phi=0$ is vanishing. The condition of vanishing curvature thus defines the $\lambda_{\rm c}(m^{2})$ curve through the relation

[TABLE]

The second line, denoted $\lambda_{-}(m^{2})$, delimits the region where a loss of solution occurs from the region where the coupled gap and field equations do admit a solution. We expect such a region based on our previous investigations in Ref. Marko:2015gpa . As an example, the left panel of Fig. 2 shows a value of the coupling ($\lambda=34$) at which the solution is lost for $m^{2}<-0.12,$ while for $m^{2}>-0.12$, $\bar{M}_{\bar{\phi}\neq 0}^{2}(m^{2})$ is multivalued. Multivaluedness is usually the signal of the appearance of unphysical branches, which can collide with the physical branch and lead to a loss of solution. In our case, the physical branch is always the upper one. This branch continuously connects to the $\phi=0$ solution (in case of parameters where the solution exists for any $\phi\geq 0$), where the unphysical branch disappears and only the physical solution remains.444Actually at all $\phi$, including the vanishing field limit, there exists yet another unphysical branch, where the gap mass is exponentially large. We disregard this solution of the gap equation when discussing physical and unphysical branches. It is important to note here that, in order to find the unphysical branch of solution, it was very helpful to reparametrize our equations in terms of the local part of the self-energy $\bar{M}^{2}(Q=0)$, rather than $m^{2}$, as explained in Appendix C around (49). For further details on unphysical branches of solutions and loss of solution we direct the reader to Ref. Marko:2015gpa . In principle, the $\lambda_{-}(m^{2})$ curve can be obtained in the full 2PI case by searching, at each value of $m^{2}$, for the value of the coupling where a loss of solution to the coupled equations occurs, as shown in the right panel of Fig. 2. Since this is a tedious procedure, as for large $\lambda$ the convergence of the iterations is slow and strongly initial condition dependent, we resort to the much simpler localized approximation used in Ref. Bordag:2000tb ; Reinosa:2011cs ; Marko:2015gpa , which qualitatively reproduces the solution of the full 2PI equations. In this approximation the full self-energy is replaced by its zero momentum value resulting in a tree-level type propagator with a self-consistent mass. Finally we define $\bar{\lambda}_{\rm c}(m^{2})$ as the curve where the gap mass at $\phi=0$ vanishes. For $\lambda<\bar{\lambda}_{\rm c}$ the gap equation admits no solution at $\phi=0$ which also makes the potential unreachable there. Using the gap equation (6a) at $\phi=0$ and that {\cal T}_{\rm F}[\bar{G}]\big{|}_{\bar{M}^{2}=0}=M_{0}^{2}/(16\pi^{2}) one can easily see that $\bar{\lambda}_{\rm c}(m^{2})=-32\pi^{2}m^{2}/M_{0}^{2}\,.$

III Analytic continuation of Euclidean solutions using Padé approximants

In this section we would like to use Padé approximants in order to analytically continue functions known only at a finite number of points in the complex plane. We employ the multipoint Padé approximant (see Refs. VS_Pade and Baker ) calculated for a function known at $N$ complex points $z_{i}$, $f(z_{i})=u_{i},$ $i=1\dots,N,$ using a finite continued fraction. Using the notation $\frac{1}{1+}x\equiv\frac{1}{1+x}$, the finite continued fraction is written in the form

[TABLE]

and the task is to determine its $N$ coefficients $a_{i}$ from the conditions $C_{N}(z_{i})=u_{i}.$ An elegant and efficient way to achieve this is by recursion: the coefficient are obtained as $a_{i}=g_{i}(z_{i})$, by defining555In the second relation it is understood that $z$ is part of the set of discrete points $z_{i}.$

[TABLE]

Working out explicitly the condition $a_{i}=g_{i}(z_{i})$ for a few values of $i$ one sees that basically one needs to construct a triangular matrix $t_{i,j}$ using the recursion $t_{i,j}=(t_{i-1,i-1}/t_{i-1,j}-1)/(z_{j}-z_{i-1}),$ for $j=2,\dots,N$ and $i=2,\dots,j$, starting from its first row $t_{1,j}=u_{j},$ $j=1,\dots,N$.

The finite continued fraction can be written as a rational function: $C_{N}(z)=A(z)/B(z)$, which we do not give here because we do not use that form. We only mention that the polynomials $A(z)$ and $B(z)$ are both of order $(N-1)/2$ for $N$ odd and of order $(N-2)/2$ and $N/2$, respectively, for $N$ even. This means that with increasing $N$ two Padé sequences of the type $P_{k}^{k}(z)$ (for $N=2k+1,k\geqslant 0$) and $P_{k+1}^{k}(z)$ (for $N=2(k+1),k\geqslant 0$) are generated. When all the continued fraction coefficients are nonnegative it is known (see Ref. Bender ) that $P_{k}^{k}(z)$ decreases, while $P_{k+1}^{k}(z)$ increases monotonically with $k$. The first sequence has a lower bound, while the second sequence has an upper bound which for $k\to\infty$ is a Stieltjes function $F(z)$, that is a function of the form $F(z)=\int_{0}^{\infty}dt\rho(t)/(1+zt),$ with $\rho(t)\geqslant 0$ in the domain of integration.

We shall use the above method to obtain the propagator $\bar{G}(Q)$ for $Q^{2}>0$, from the knowledge of the same propagator at a finite number of negative $Q^{2}=-Q_{\rm E}^{2}$ values in the Euclidean domain. In practice we fix $\vec{q}$ (to zero in our case), define the original Padé approximant as a function of the Matsubara frequency $\omega_{n}$ and, to carry out the analytic continuation, we evaluate it at $\omega_{n}=-i\omega+\varepsilon$, where the $\varepsilon\to 0$ limit can be safely taken. We should of course keep in mind that the analytic continuation of a finite number of data is not unique and leads to various solutions that differ by their asymptotic behaviors as $|z|\to\infty$.666Knowing some data $d_{1},\dots,d_{n}$ for a finite set of points $z_{1},\dots,z_{n}$, and given an analytic continuation $f(z)$ of these data, which is such that $f(z_{i})=d_{i}$, we can construct another analytic continuation of the same data as, for instance, $g(z)=f(z)+\prod_{j=1}^{n}(\exp\{i2\pi z/z_{j}\}-1)$. The ability of the method to reproduce the expected propagator needs then to be tested and we will do so by using the explicit Minkowskian solution obtained in the previous section. Another subtle point is that, by definition, the Padé approximant has no branch cut and thus there is only one Riemann sheet to be considered. It is then a question how the method can allow us to access physical information usually hidden in the second Riemann sheet of the propagator. A related issued it the fact that the Padé approximant does not obey the Schwarz reflection property, whereas the analytic propagator (39) obeys this property (${\cal G}(s^{*})^{*}=-{\cal G}(s)$ with our conventions) in the first Riemann sheet, as it is easily checked.

III.1 Test I: ability to access the spectral function

To test the quality of the Padé analytic continuation we use two test functions: the Euclidean 2PI self-energy extracted from Eq. (19) using the spectral function obtained as the solution of (6) and the zero temperature finite perturbative Euclidean bubble with square mass $M^{2}=0.5$:

[TABLE]

In both cases we compare the real and imaginary parts of the analytically continued quantity to the one computed in the Minkowski space, which in the latter case is $\bar{M}^{2}(Q)$. For the comparison in the 2PI case we use the parameters $m^{2}=-0.2611$ and $\lambda=5$, for which one obtains $M^{2}\equiv\bar{M}^{2}(Q=0)=0.5$ and $\bar{\phi}=0.5566.$ These agree with the values obtained with the numerical code used in Ref. Fejos:2011zq to solve the Euclidean version of the model in the present renormalization scheme.

In Fig. 3, we show the Padé analytically continued real and imaginary parts of the Euclidean 2PI self-energy and the corresponding spectral function compared to the quantities computed directly in Minkowski space. We see that while using $N=10$ points in (21) leads to a noticeable error almost everywhere, except for the real part of the self-energy at small $\sqrt{s}$, the multipoint Padé approximation using $N=50$ points has difficulties reproducing the Minkowski result only in a narrow neighborhood of the threshold. However, even with $N=10$ points, the analytical continuation gives a good approximation for the pole mass, $\bar{M}_{\rm p}=0.70387$, compared to the value $\bar{M}_{\rm p}=0.70392$ obtained from a direct calculation. Based on this observation, we can trust our results presented in Sec. IV concerning the comparison of zero-momentum and pole masses.

In order to quantify the quality of the analytic continuation over a significant region of momentum, we define the following Minkowski space integral

[TABLE]

where $F$ is the original function and ${\cal F}$ its approximation obtained using analytic continuation, and with $\sqrt{s}_{\rm min/max}$ chosen appropriately. In Fig. 4 we show this quantity ${\cal Q}$ as a function of the sampled region with fixed sampling frequency (bottom $x$-axis), and as a function of the number of sample points on fixed region (top $x$-axis). We see convergence both by increasing the sampled region and by increasing the number of points. While a perfect continuation would mean ${\cal Q}=0$ the curves do not tend to zero, which we attribute to the error arising from finite number representation and the inherent differences in the analytic properties of the approximated functions and the Padé approximants. Also note that the combined results have some implications on the continuation of finite temperature data. The results corresponding to the bottom axis tell us that, at a certain temperature, there is a limiting Matsubara-frequency, over which including new frequencies does not improve the quality of the continuation. At the same time, the results associated to the top axis indicate that with fewer frequencies in the relevant range, i.e. by increasing the temperature, the quality of the continuation deteriorates.

III.2 Test II: ability to access physical information in the second Riemann sheet

We show now that the Padé analytic continuation is capable of finding physically relevant complex poles in the case of an $O(4)$ symmetric model. Since we have 2PI results in Minkowski space only in the one-component case, we use as a benchmark the $T=0$ Minkowskian results of Chiku:1998kd and Hidaka:2003mm obtained within perturbation theory at one loop. We calculate the integrals of the sigma self-energy both in Euclidean and Minkowski spaces with the renormalization prescription of Chiku:1998kd and then carry out the analytic continuation (cf. Ch. 6.3 of Brown ) of the Euclidean bubble integrals using Padé approximants. We construct one Padé approximant for each bubble integral, as it turns out that with this procedure we get better agreement with the Minkowskian results than in the case where a single Padé approximant is associated to the entire Euclidean self-energy.

The physically relevant complex pole of the sigma propagator is on the 2nd Riemann sheet (see e.g. Sec. 1.3 of Badalyan ), which is accessed by crossing the real axis in between the thresholds of the pion and the sigma bubbles. In the present case this range is given in terms of the tree-level pion and sigma masses as $\sqrt{s}\in[2m_{\pi,0},2m_{\sigma,0}]$.

In search for the complex pole of the $\sigma$ propagator, we continue the Padé approximant (21) to complex values of its argument. With the parameters given in the first line of Table I in Ref. Chiku:1998kd , we find a complex pole at $\sqrt{s}=569.25-i\cdot 119.02$ MeV, using the exact formulas of the bubbles analytically continued to the 2nd Riemann sheet, and at $\sqrt{s}=569.32-i\cdot 118.83$ MeV, using the multipoint Padé approximation with $N=200$ points in the $[0,3]$ GeV interval of Euclidean momentum. This finding is in line with the low temperature location of the pole II in the 2nd Riemann sheet shown in Fig. 6 of Ref. Hidaka:2003mm .777However, we could not reproduce the results shown in Fig. 2 of Ref. Hidaka:2002xv , where the scale of the imaginary axis seems to be off. As in Fig. 3 of Ref. Hidaka:2002xv , the real part of this pole is close to the maximum of the spectral function, which is accurately reproduced with a Padé approximant using $N=200$ points, in line with the findings of the previous subsection. However, a word of caution is in order, since as we show in Fig. 5, far from the real axis the contours $\Re G_{\sigma}^{-1}(\sqrt{s})=0$ and $\Im G_{\sigma}^{-1}(\sqrt{s})=0$ are strongly distorted due to spurious poles of the Padé approximants fitted to the bubble integrals. In the shown range there are eight fake poles, four to the left and four to the right of the physical one. While in our benchmark case we are lucky, as the position of the physical pole is well captured with the Padé analytic continuation, one should always check whether the result obtained in the complex plane is influenced by non-physical poles of the Padé approximant.

The fact that the Padé analytic continuation gives so accurately the coordinates of a complex pole on the 2nd Riemann sheet located far away from the real axis is a surprise to us in view of the fact that, as already mentioned, the Padé approximant does not know about the existence of different Riemann sheets. When using Padé approximants one can in principle mimic to some extent the procedure of the analytic continuation based on a known functional form. In order to do this, we continue the Euclidean bubbles to the real axis and then we construct new Padé approximants for the real and imaginary part of the pion bubble for $\sqrt{s}>2m_{\pi,0}$ and a new Padé approximant for the the sigma self-energy in the range $\sqrt{s}\in[2m_{\pi,0},2m_{\sigma,0}]$ where it is purely real. One can then continue these Padé approximants to complex values. Actually, Fig. 5 was obtained in this way, and in the present case no significant change in the coordinates of the complex pole and of the zero contours was observed for $\Re\sqrt{s}<2m_{\sigma,0}$.

IV Comparison of pole and zero-momentum masses in the O(N) model

In our earlier works Marko:2013lxa ; Marko:2015gpa , we have studied the parametrization of the $O(4)$-model in regard to its application to light meson phenomenology. The parametrization was done using zero-momentum masses for simplicity. However, a correct parametrization requires, instead, the use of pole masses. Here we use the Padé analytic continuation method in order to assess how different these masses can be. We should mention that this comparison only makes sense in a given renormalization scheme because the zero-momentum masses are scheme dependent whereas the pole masses are not.

Another subtlety comes from the fact that, in a given 2PI approximation, one can define two types of propagators, the so-called internal and external propagators. The internal propagator is the one considered in Sec. II for the case of the two-loop approximation. However, one can define an external propagator, as we recall in the next subsection. Since the parametrization of Ref. Marko:2013lxa ; Marko:2015gpa was done using the zero-momentum mass of the external propagator, in what follows we compare this mass with the pole mass of the same propagator. In doing so, we shall use a well motivated approximation for the external propagator which leads to the same expression in the two-loop and ${\cal O}(\lambda^{2})$ truncations of the 2PI effective action, the two approximations used in Refs. Marko:2013lxa ; Marko:2015gpa . For completeness, we shall also do this comparison in the case of the internal propagators which are different in the two types of truncations.

IV.1 External vs internal propagators

In the 2PI formalism, the internal propagator $\bar{G}_{\phi}$ in the presence of a given field expectation value is obtained by solving the equation $0=\delta\Gamma[\phi,G]/\delta G_{G=\bar{G}_{\phi}}$, where $\Gamma[\phi,G]$ is the so-called 2PI effective action. Knowing $\bar{G}_{\phi}$, one can reconstruct the usual 1PI effective action as $\Gamma[\phi]=\Gamma[\phi,\bar{G}_{\phi}]$. From the latter, it is also possible to compute the propagator as ${\cal G}^{-1}=\delta^{2}\Gamma/\delta\phi^{2}$. Although these two definitions for the two-point function are equivalent in the absence of approximations, they differ in practice whenever a truncation of $\Gamma[\phi,G]$ is considered and the second definition is referred to as the external propagator. We refer to Berges:2005hc for details on the relation between the two types of propagators and their respective renormalization. Here we recall the expression for the external propagator in the two-loop and ${\cal O}(\lambda^{2})$ truncations, in the $N=1$ case. In Euclidean configuration space, one obtains

[TABLE]

with

[TABLE]

where we have restricted to homogeneous fields $\phi(x)=\phi$ and where $\Lambda(x;z_{1},z_{2})$ and $V(z_{1},z_{2};x)$ are proportional to $\phi$ and have tree-level contributions equal to $\lambda_{2}\phi\,\delta(x-z_{1})\delta(x-z_{2})$. In this work, we shall neglect all contributions to $\Lambda$ and $V$ beyond tree-level, which is the minimum needed if we want the two-loop external propagator to contain at least the same diagrams than the two-loop internal one. Higher order contributions, since they are proportional to $\phi^{2}$, lead only to logarithmic modifications of the momentum dependence of the self-energy. We thus expect those contributions not to modify by much the values of the pole and zero-momentum masses. This is an important simplification because, in general, obtaining the function $V$ requires to solve a Bethe-Salpeter equation, which we avoid here without affecting our conclusions in a dramatic way. Putting these pieces together and Fourier transforming, we arrive at ${\cal G}^{-1}_{\rm E}(Q_{\rm E})\equiv Q_{\rm E}^{2}+\hat{M}_{\rm E}^{2}(Q_{\rm E})$, with

[TABLE]

We mention that, in the ${\cal O}(\lambda^{2})$ truncation, the internal propagator has exactly the same expression as the approximate external propagator introduced above, since the gap mass $\bar{M}_{\rm E}^{2}(Q)$ is given by the right hand side of (27).

In fact, we can do a little bit better concerning the external propagator, as the zero-momentum value $\hat{M}_{\rm E}^{2}\equiv\hat{M}_{\rm E}^{2}(Q_{\rm E}=0)$ does not need to be approximated, but can be computed instead from curvature of the effective potential associated to the effective action. Proceeding this way, one completely takes into account the contribution of the Bethe-Salpeter equation to $\hat{M}_{\rm E}^{2}$ and away from $Q_{\rm E}=0$ we can approximate $\hat{M}_{\rm E}^{2}(Q_{\rm E})$ by

[TABLE]

The counterterm $\delta Z_{2}$ is needed to absorb divergences in the difference of setting-sun integrals and is determined numerically as explained in Appendix A.2 of Ref. Marko:2015gpa , using instead of $T_{\star}$ a small value of the temperature, $T=0.05T_{\star}.$ If we denote by $\hat{M}^{2}(Q)$ the analytic continuation of $\hat{M}^{2}_{\rm E}(Q_{\rm E})$, the pole mass of the external propagator corresponds to $\hat{M}^{2}_{\rm p}=\hat{M}^{2}(\hat{M}^{2}_{\rm p})$.

The generalization of Eq. (28) to the $O(4)$ case is straightforward. To get the approximated $\hat{M}^{2}_{\rm L,E}(Q)$ and $\hat{M}^{2}_{\rm T,E}(Q)$ in the corresponding longitudinal and transverse external propagators, one can use the expressions given in Eq. (15) of Ref. Marko:2015gpa , subtract from the right hand side the zero momentum expression and add the corresponding curvature mass determined numerically using Eq. (16) of that reference. We emphasize that the approximated expressions are the same in the two-loop and ${\cal O}(\lambda^{2})$ truncations, but the internal propagators entering the integrals are different in the two cases.

IV.2 On the quality of the curvature mass based parametrization of the O(4)

model

We now turn to the analytic continuation of the $N=4$ Euclidean results obtained in Marko:2013lxa ; Marko:2015gpa at a small value of the temperature, $T=0.05T_{\star}$. Applying the method based on Padé approximants we find similar results in the two investigated truncations for any parameters where the curvature masses (zero-momentum masses of the external propagator) are physical. The main result is that although the difference between the pole masses of the longitudinal external and internal propagators can be of 40%, the difference between the pole and zero momentum mass of a given propagator (internal or external) is only within 1%. For the transverse mode the differences are typically smaller, which is in line with the results of Strodthoff:2016pxx ; Helmboldt:2014iya obtained using the FRG approach. We already knew from previous studies that the difference between $\hat{M}_{\rm L}$ and $\bar{M}_{\rm L}$ could be more than 30% for some parameters (see Fig. 3 of Marko:2015gpa ) and now we see that both $\bar{M}_{\rm L/T}-\bar{M}_{\rm L/T,p}$ and $\hat{M}_{\rm L/T}-\hat{M}_{\rm L/T,p}$ are very small in both truncations, for parameters where the curvature masses are physical. For such parameters the pole of the analytically continued sigma propagator becomes complex. Based on the study presented in Sec. III.B, we interpret this pole as being the physical pole on the 2nd Riemann sheet. Unfortunately, the imaginary part of the pole is almost zero: the ratio between real and imaginary part is $\sim 10^{-3}$. We illustrate this for a typical parameter set of the two-loop truncation in Fig. 6, where the plots are obtained by continuing the Euclidean external longitudinal propagator calculated with the approximation described above.

The smallness of the imaginary part of the complex $\sigma$ pole is in line with the findings of Ref. Patkos:2002xb , where it turned out that at leading order in a $1/N$ expansion its ratio to the real part of the pole only starts to grow for larger values of the coupling and becomes physically acceptable for $\lambda\in(300,400).$ We mention that taking into account the scaling by $N=4$ used in Refs. Marko:2013lxa ; Marko:2015gpa the coupling of Ref. Chiku:1998kd (also used in Sec. III.B) corresponds to $\lambda=292.$ Even for such a big value of the coupling constant $M_{\rm L}(Q=0)$ is only 1.7% larger than $M_{\rm L,p}$ given in Sec. III.B. Such a large value of the coupling, which is by a factor of 5 larger than the largest coupling used for parametrization in Refs. Marko:2013lxa ; Marko:2015gpa , cannot be reached in our 2PI investigations because the closeness of the Landau pole to the physical scales makes the solution of the propagator equation, if accessible at all iteratively, highly cutoff sensitive.

IV.3 Infrared sensitivity of the pole mass in the two-loop SI2PI approximation

In Ref. Marko:2016wtw , we investigated the symmetry improved 2PI (SI2PI) formalism of Ref. Pilaftsis:2013xna at two-loop order, using various types of UV regulators. In particular, we pointed out that, for generic smooth regulators, the solution for the internal propagator of the SI2PI framework possess an untamed infrared sensitivity in the broken phase that leads to a loss of solution for large enough volumes. This may not be a problem though in cases where the volume at which the solution disappears is many order of magnitudes higher than the physical volume of the system under study, and that the quantity under scrutiny presents a plateau behavior for a large range of volumes below the volume at which the solution is lost. In Ref. Marko:2016wtw , we tested this scenario on the zero-momentum mass of the Higgs propagator in the model of Ref. Pilaftsis:2013xna and we observed sensible changes with the volume as we also illustrate in Fig. 7. However, as pointed out to us by D. Teresi, in the same model, the Higgs pole mass, should be (almost) insensitive to the volume. This insensitivity was reported in Ref. Pilaftsis:2017enx and it is expected based on the fact that the Higgs pole mass is not related to the zero-momentum behavior of the propagator, which is the only region where the propagator is infrared sensitive in our case. However, we note that it is not clear which type of UV regulator was used in Refs. Pilaftsis:2013xna ; Pilaftsis:2017enx , so one could a priori think that the insensitivity could also originate from the use of a sharp regulator that tames artificially the IR sensitivity, as explained in Ref. Marko:2016wtw . We can now test these scenarios using the Padé continuation technique with various types of UV regulators. The result for the Higgs pole mass is compared to the zero-momentum Higgs mass in Fig. 7. Imposing the constraint of the SI2PI formalism as $\bar{M}_{\rm T}(|K_{\rm E}|=\kappa)=0$, we observe that, even with a smooth UV regulator, for which the internal Higgs propagator ceases to exist beyond some volume, the corresponding pole shows a plateau behavior over a large range of volumes.

V Conclusions

We searched for broken phase Minkowski space solutions in the $(m^{2},\lambda)$ parameter space of the one-component $\varphi^{4}$ model at zero temperature using the 2PI effective action truncated at two-loop level. Since the Euclidean solution can be obtained from the spectral function, we could use the Minkowskian solution as a benchmark for the analytic continuation of the Euclidean propagator with the method based one multipoint Padé approximants. The various tests we performed show that at zero temperature the Padé analytic continuation is an efficient tool to obtain not only the spectral function, but also the complex pole of the propagator on the second Riemann sheet. However, it is expected that as the temperature grows the applicability of the method is limited by the decreasing number of available Matsubara modes in a given frequency interval, as we have argued in Sec. III.1.

We applied the Padé analytic continuation on previously determined practically zero temperature Euclidean propagators of the $O(4)$ model and studied the relation between the pole and the gap mass at vanishing external momentum of both the internal and external 2PI propagators. We have searched for real poles of the transverse propagator and for complex poles of the longitudinal propagator on the 2nd Riemann sheet. For parameters where the longitudinal curvature mass (gap mass at vanishing external momentum of the external propagator) is physical, the real part of the corresponding pole is within 1% of the curvature mass, both at two-loop and at ${\cal O}(\lambda^{2})$ level truncations of the 2PI effective action. This shows that the error of the parametrization of Marko:2013lxa and Marko:2015gpa , where the curvature masses where used instead of the pole masses, is smaller than 1%. However, with the parameters determined in these references the imaginary part of the pole of the longitudinal external propagator turns out to be very small, although it increases with the value of the coupling constant. Unfortunately, the presence of the Landau pole prevents us from accessing regions of the parameter space with larger values of the coupling, as the solution of the model, which is increasingly difficult to find with an iterative method, becomes highly cutoff sensitive.

Acknowledgements.

We would like to thank D. Teresi for clarifying discussions about his work that motivated the analysis presented in Sec. IV.3. In case of G. M.  this work is part of Project no. 121064 for which support was provided from the National Research, Development and Innovation Fund of Hungary, financed under the PD_16 funding scheme. Zs. Sz. would like to thank École Polytechnique for hospitality and support during the last stages of this work.

Appendix A Renormalization

As mentioned in the Introduction above Eq. (6), we now discuss how to obtain the counterterms with a renormalization at $T=0,$ in a way that resembles the renormalization prescriptions at $T=T_{\star}$ which were used in some of our previous works. Since there is no physical renormalization condition giving exactly the counterterms determined in Patkos:2008ik , we define $m_{0}^{2}$ and $m_{2}^{2}$ by absorbing those quadratic and logarithmic divergences which arise when expanding the propagator $\bar{G}$ around $G_{0}(Q)=i/(Q^{2}-M_{0}^{2}+i\varepsilon)$ and do not depend on the environment ($\phi$ or ${\cal T}_{\rm F}[\bar{G}]$). One has

[TABLE]

where the environmental free divergences of the tadpole and setting-sun integrals are given by

[TABLE]

with $T_{\rm d}^{(I)}=-i\int_{Q}G_{0}^{2}(Q){\cal I}_{\rm F}[G_{0}](Q).$

In order to renormalize the couplings we proceed as usual at $\phi=0$ (see Ref. Berges:2005hc for details) and introduce the 2PI kernels and the related four-point functions, but formally replace in them and in the equations they satisfy $\bar{G}_{\phi=0}$ by $G_{0}.$ One has for instance

[TABLE]

where at the present order of truncation

[TABLE]

$\lambda_{0}$ is obtained from the condition $\bar{V}_{0}=\lambda$ in the form $\lambda_{0}^{-1}=\lambda^{-1}+{\cal I}[G_{0}]/2$, while $\lambda_{2}$ from the condition $V_{0}(0,K)=\lambda,$ which gives after a few lines of algebra $\lambda_{2}=\lambda_{0}-\lambda^{2}{\cal I}[G_{0}]-\lambda^{2}\lambda_{0}T_{\rm d}^{(I)}/2.$ The counterterm $\delta\lambda_{4}$ is obtained from the condition $\hat{V}_{0}=\lambda,$ using the equation for

[TABLE]

Appendix B Dispersion relations and other useful relations

Let us start recalling here Eqs. (9) and (10) can be understood from the sole assumption of the existence of an analytic propagator ${\cal G}(Q^{2})$ for complex values of $Q^{2}$ away from the positive real axis (in the first Riemann sheet), that decreases as $1/Q^{2}$ at large $Q^{2}$ and which gives back the Minkowski propagator when approaching the positive real axis: $G(Q)={\cal G}(Q^{2}+i\epsilon)$ with $Q^{2}>0$. Indeed, the propagator ${\cal G}(Q^{2})$ being analytic away from the positive real axis, one can apply Cauchy formula to a small circle not intersecting this axis and centered around a certain $Q^{2}$. Deforming this contour into a large circle $C_{R}$ whose radius $R$ tends to infinity and a contour $B_{R}$ with two parallel branches on each side and close to the positive real axis, one obtains

[TABLE]

Owing to the assumed behavior of ${\cal G}(s)$ at large $s$, the first term goes to [math] as $R\to\infty$. The second term can be expressed in terms of the spectral function and we finally obtain

[TABLE]

which provides a spectral representation for the analytic propagator that leads to Eq. (9) when approaching the positive real axis. The sum rule is obtained in a similar fashion by applying the argument to the function $(s-Q^{2}){\cal G}(s)$ instead of ${\cal G}(s)$. In this case the integral over the contour $C_{R}$ does not vanish and gives the $1$ in the right-hand-side of Eq. (10).

The previous argument can be repeated for any function ${\cal F}(Q^{2})$ analytic in $Q^{2}$ away from the positive real axis. In the case where ${\cal F}(s)$ does not go to zero fast enough, one should apply the argument to ${\cal F}(s)-\sum_{p=0}^{n}{\cal F}^{(p)}(s_{0})(s-s_{0})^{p}/p!$ with $n$ high enough, such that the integral over $C_{R}$ vanishes as $R\to\infty$. For instance, if ${\cal F}(s)$ grows logarithmically, we apply the argument to ${\cal F}(s)-{\cal F}(s_{0})$ with $s_{0}=Q_{0}^{2}$ away from the positive real axis. We obtain, for any $Q^{2}$ and $Q_{0}^{2}$ away from the positive real axis,

[TABLE]

In particular, if one is interested in $F(Q)\equiv{\cal F}(Q^{2}+i\epsilon)$, one has for $Q^{2}$ and $Q_{0}^{2}$ real,

[TABLE]

The previous formula applies in particular to the finite bubble integral considered in the present work. Indeed, using the spectral representation and its sum rule,

[TABLE]

from which it follows that this is analytic in the variable $Q^{2}$ away from the positive real axis. Moreover since the finite bubble grows like $\ln Q^{2}$, we need to apply the once-subtracted formula.

Let us now derive some other useful results Sauli:2004bx . We start from the identity $\bar{G}(Q)\bar{G}^{-1}(Q)=1$ and derive two expressions for the finite wave-function renormalization constant $Z,$ the expression for the continuum part of the spectral function and the sum rule it satisfies. For $\bar{G}(Q)$ we use its spectral representation (9), while for the inverse propagator we use the expression $\bar{G}^{-1}(Q)=-i(Q^{2}-\bar{M}^{2}(Q)+i\varepsilon).$ The gap mass is a complex valued function and $\Re\bar{M}^{2}(Q)$ and $\Im\bar{M}^{2}(Q)$ can be read off from (6a), remembering that only the bubble integral has imaginary part. Using the form of the spectral function given in (16) and the relation

[TABLE]

one obtains from $GG^{-1}=1$ two equations, one for the real part and one for the imaginary part. With some algebraic manipulations one derives from them the usual expression for the continuum part of the spectral function

[TABLE]

and the equation ($s_{\rm th}=4\bar{M}^{2}_{\rm p}$)

[TABLE]

Taking the limit $Q\to\infty$ in (45) one obtains

[TABLE]

which is the sum rule for a spectral function of the form (16), while taking the limit $Q\to\bar{M}_{\rm p},$ one obtains

[TABLE]

which is the residue of the pole $Q=\bar{M}_{\rm p}$ of the propagator.

Appendix C Numerical algorithm and its implementation

In this Appendix we describe the steps of the iterative process applied to solve sequentially the explicitly finite coupled field and gap equations (6) and present the numerical implementation of these steps.

Being interested in broken symmetry phase solutions, we start from an initial propagator $\bar{G}^{(0)}$ and express the nontrivial $(\bar{\phi}^{2})^{(1)}$ from the field equation. Then, using both $\bar{\phi}^{(1)}$ and $\bar{G}^{(0)},$ we evaluate $(\bar{M}^{2})^{(1)}$ from (6a), which in turn gives us $\bar{G}^{(1)}(Q)$. These steps are repeated until the relative change from one iteration to the next of both $\bar{\phi}$ and the pole mass $\bar{M}^{2}_{\rm p}$ obtained from (17) is less than our stopping parameter, chosen to be $10^{-7}.$ If $(\bar{\phi}^{2})^{(i+1)}<0$ we leave $\bar{\phi}^{(i+1)}=\bar{\phi}^{(i)}$. When for 5 consecutive iterative steps $\bar{\phi}^{2}$ is negative, we stop the iteration and say that no broken phase solution was found for the given parameters and/or initial values.

The explicit steps for one full iteration (from the $i$-th to $i+1$-th) are the following. The propagator $G^{(i)}$ is used to denote a given quantity in the $i$th iteration, even though it is clear from Sec. II.2 that the spectral function of the form (16) is the central object of the entire iterative procedure.

Using (7) and (8), evaluate ${\cal S}_{\rm F}[\bar{G}_{\rm E}^{(i)}]$ and ${\cal T}_{\rm F}[\bar{G}_{\rm E}^{(i)}]$ and update $\bar{\phi}^{(i+1)}$ from (6b). 2. 2.

Evaluate $\Im{\cal I}[\bar{G}^{(i)}](Q)$ using (11) and then $\Re{\cal I}_{\rm F}[\bar{G}^{(i)}](Q)$ using the dispersion relation (18). 3. 3.

Determine $(\bar{M}^{2}_{\rm p})^{(i+1)}$ solving the pole equation (17) and update $\sigma^{(i+1)}(s)$ using its expression given in (44), then update $Z$ using the sum rule (46).888Note that by changing $\bar{M}^{2}_{\rm p}$ in each iteration $s_{\rm th}$ also changes. Due to the structure of the bubble integral $(s_{\rm th})^{(i+1)}=4(\bar{M}^{2}_{\rm p})^{(i)}$.

We initialize $\bar{\phi}$ with some crude non-zero estimate of the solution (e.g. the classical value $\sqrt{-m^{2}/\lambda}$) and choose the initial spectral function to contain only a pole part with unit residue, with the pole at $\bar{M}_{\rm H}$, defined as the solution of the gap equation in the Hartree approximation

[TABLE]

with $\bar{G}_{\rm H}(Q)=i/(Q^{2}-\bar{M}_{\rm H}^{2}+i\varepsilon)$.

We note that one can introduce a reparametrization of the equations (6) in terms of the mass parameter

[TABLE]

that is the zero momentum value of the self-energy. While this introduces some intricacy connected to the fact that $M^{2}$ depends on the a priori unknown solution $\bar{\phi},$ it has several advantages. First, by fixing the local part of the self-energy, the number of iterations required for convergence is greatly reduced. Also, the tadpole diagram (2), one of the main sources of numerical error, only has to be evaluated to obtain $m^{2}$ corresponding to the fixed value of $M^{2}$. Second, this allows the introduction of an underrelaxation parameter $\alpha$ as

[TABLE]

This is needed to achieve convergence in parameter regions where otherwise our iterative method fails. This combined with the good initial guess provided by $M^{2}$ allowed us to find the unphysical solution shown in Fig. 2.

C.1 Numerical implementation

In the iterative method described above we have to store $\bar{\phi},\,\bar{M}^{2}_{\rm p},$ and discretized version of $\sigma(s)$. For the latter we use a grid of $N_{s}$ points (typically we use $N_{s}=500$) defined actually in the transformed variable $t=\frac{s}{1+s}$. Starting from $t(s=s_{\rm th})$ and ending in $1$ one has

[TABLE]

with $\kappa=\frac{s_{\rm th}}{s_{\rm th}+1}$ and $\Delta=\frac{1-\kappa}{(N_{s}-1)^{3}}$. The property $t(s=\infty)=1$ allows us to carry out the $s$-integrals of Sec. II.2 without introducing an explicit numerical cutoff.

We follow below the order of the iterative steps described above and give some technical details concerning the evaluation of the integrals.

The explicitly finite expression of the two local Euclidean integrals, the tadpole ${\cal T}_{\rm F}[\bar{G}_{\rm E}]$ and ${\cal S}_{\rm F}[\bar{G}_{\rm E}]$ is given in Appendix D, in (54) and (55)–(56) respectively, in terms of the Euclidean propagator (19). The integral of (19) is evaluated on a rather complicated looking momentum grid defined as

[TABLE]

with $\Delta(i)=\Lambda/(N_{s}-1)^{\alpha(i)},$ where $\Lambda$ is an appropriate numerical cutoff for which we typically use values in the range $(50-200).$999In principle we could also use a transformed variable similar to (51) and integrate up to infinity since the integrals are finite. However, the integrands of (54) and (55)–(56) only decrease fast enough due to cancellations, which break down when numerical number representation errors are comparable to the results. The reason for using such a grid is to sample with sufficient accuracy the Euclidean propagator both in the IR and UV regions. For integration we need to know $G_{\rm E}$ in between the grid points, but in order to improve on the UV behavior, we subtract $1/(q_{i}^{2}+M_{0}^{2})$ from $G_{\rm E}(q_{i})$ and fit the difference with a Steffen-spline.101010This is used because Steffen’s method leads to an interpolation function which is monotonic between the given data points. This proves to be accurate enough in the UV to see the apparent convergence of the integrals with increasing $\Lambda$. The spline is evaluated in the integrand called by the cquad adaptive integration routine of GSL library GSL .

To compute $\Im{\cal I}[\bar{G}](Q),$ we use (16) in (11). This leads to $\delta\times\delta$, $\delta\times\sigma$, and $\sigma\times\sigma$ type terms, which start giving non-vanishing contribution for $Q>2\bar{M}_{\rm p},$ $Q>3\bar{M}_{\rm p},$ and $Q>4\bar{M}_{\rm p},$ respectively. Since $\Im{\cal I}[\bar{G}](Q)$ will be used in the dispersion relation (18) giving $\Re{\cal I}_{\rm F}[\bar{G}](Q)$, we spline it in order to be able to evaluate it anywhere. We use the grid introduced for $\sigma$ in (51) in the same transformed coordinates. However, a spline interpolation cannot reproduce a functional form with infinite derivatives, such as $\Im{\cal I}[\bar{G}](Q),$ which behaves as a square root function around $Q=2\bar{M}_{\rm p}.$ Therefore, in order to avoid this problem, we chose to fit with a cubic spline the function $\left(\Im{\cal I}[\bar{G}](Q)\right)^{2}.$

$\Re{\cal I}_{\rm F}[\bar{G}](Q)$ is evaluated from the dispersion relation (18). In order to circumvent the problem of numerically computing integrals with a principal value prescription when $s=Q^{2}>4\bar{M}^{2}_{\rm p},$ we rewrite them with the usual method (see e.g. Shima ). First we split the interval of integration by introducing the point $2s-4\bar{M}^{2}_{\rm p},$ then we add and subtract an appropriate term in the integrand of the resulting integral over the interval $[4\bar{M}^{2}_{\rm p},2s-4\bar{M}^{2}_{\rm p}],$ exploiting the fact that ${\cal P}\int_{a}^{2b-a}dx/(x-b)=0$. In this way we obtain

[TABLE]

Actually, the integrand of the first integral on the r.h.s. is continuous at $s^{\prime}=s,$ so that the principal value prescription can be omitted when integrating it numerically with the cquad routine.

We end this part by mentioning two numerical intricacies. The first is encountered when fitting $\sigma$ with a spline on the grid (51). $\sigma$ grows out at the threshold as a square root function, similarly to $\Im{\cal I}[\bar{G}]$. However, contrary to $\Im{\cal I}[\bar{G}]$, which enters integrals which decrease fast enough in the UV (see the suppression by powers of $s$ in (53)), the convergence of the integrals with $\sigma$ (see e.g. (11) and (46)) is assured by the UV behavior of $\sigma$, which has to be preserved by the fit. This time, fitting $\sigma^{2}$ would not be enough because when computing integrals using the transformed variable $t$ (see the definition in the paragraph before (51)) the Jacobian $J(t)=(1-t)^{-2}$ of the transformation is divergent as $t\to 1.$ This would lead in the integrals involving $\sigma$ to a blowing-up of the oscillations of the spline when $t\to 1.$ Since $\underset{t\to 1}{\lim}\,J(t)\sigma(t)$ is finite, the way out is to fit the function $(J(t)\sigma(t))^{2}$ with a spline on the grid introduced in (51). In this way the spline respects the requirements of having appropriate behaviors around the threshold and in the UV.

The second intricacy is related to the coefficient $Z$ of the pole part of the spectral function. This has to be evaluated from iteration to iteration either from the sum rule (46) or as the residue of the pole (47). While in principle the two ways of computing $Z$ are equivalent, numerical differences could occur, as we show in Fig. 8. There we present the value of $Z$ computed in the above mentioned two ways at the end of the iterative method. It turned out that for the iterations to converge, it is actually very important to ensure that the sum rule is satisfied iteration by iteration, since this guarantees the proper asymptotic behavior of the propagator, which was also used to renormalize the real part of the bubble (see Sec. II.2). Hence, when solving the equations, $Z$ is always evaluated from the sum rule (46). The difference seen in Fig. 8 has purely a numerical reason. Increasing the precision of the numerical integration substantially reduces the difference which grows with increasing $\lambda.$

Appendix D Explicit form of some integrals

In this Appendix we give the expression of the finite Euclidean tadpole (7a) and setting-sun (7c) integrals after analytically performing as many of their angular integrals as possible. The cutoff regularization used is such that an integral involving at least two propagators is invariant with respect to shifts of the loop momenta (see e.g. Appendix D of Reinosa:2011cs ). The parts containing $G_{\rm r}$, that is (7a) and the third term of (7c), reads ($q=|Q_{\rm E}|$)

[TABLE]

[TABLE]

where

[TABLE]

with $\delta G_{\rm E}(q)=G_{\rm E}(q)-G_{\rm 0,E}(q)$ and $G_{\rm 0,E}(q)=(q^{2}+M_{0}^{2})^{-1}.$ Observing that the other two terms of (7c) can be combined, one obtains

[TABLE]

with $\lambda(q^{2},k^{2},p^{2})$ defined after (11).