The $O(4)$ $\phi^4$ model as an effective light meson theory: A lattice-continuum comparison

TL;DR

This paper compares lattice Monte Carlo simulations and continuum methods for the 4D $O(4)$ $^4$ model, assessing its viability as an effective light meson theory and exploring triviality bounds and cutoff effects.

Contribution

It provides a detailed lattice-continuum comparison of the $O(4)$ $^4$ model, highlighting limitations in cutoff independence for thermal quantities.

Findings

Triviality bounds restrict lattice spacings for the model.

Cutoff independence is practically unachievable for thermal quantities.

Matching critical lines shows consistency between lattice and continuum approaches.

Abstract

We investigate the possibility of using the 4 dimensional $O(4)$ symmetric $\phi^4$ model as an effective theory for the sigma-pion system. We carry out lattice Monte Carlo simulations to establish the triviality bound in the case of explicitly broken symmetry and to compare it with results from continuum functional methods. In case of a physical parametrization we find that triviality restricts the possible lattice spacings to a narrow range, therefore cutoff independence in the effective theory sense is practically impossible for thermal quantities. We match the critical line in the space of bare couplings in the different approaches and compare vacuum physical quantities along the line of constant physics (LCP).

Tables

Table 1. Table 1: Field and curvature mass values in units of the cutoff at the minimum of the potential of the LCP shown in Figs.  2 and 4 . The points are denoted by P i subscript 𝑃 𝑖 P_{i} with i ∈ 1 , … , 10 𝑖 1 … 10 i\in 1,\dots,10 in increasing order from left to right of the LCP. For the first six (fixed h ℎ h ) points the values comes from the direct numerical solution of ( 30 ), while for the last four (fixed λ 𝜆 \lambda ) points the values come from the solution obtained using the expansion ( 32 ) with an extrapolation to N g = ∞ subscript 𝑁 𝑔 N_{g}=\infty .

Point	$m_0^2{a}^2$	$g_0$	$H_0{a}^3$	$a$[fm]	$\overline{\varphi }/\Lambda$	${\hat{M}}_T/\Lambda$	${\hat{M}}_L/\Lambda$	$\Lambda \left[\mathrm{MeV}\right]=\frac{6.923}{a\left[\mathrm{fm}\right]}$	$\overline{\varphi }$[MeV]	${\hat{M}}_T$[MeV]	${\hat{M}}_L$[MeV]
P1	-2.087$E+1$	1.648$E+2$	2.599$E-1$	1.022$E+0$	6.471$E-2$	1.100$E-1$	2.712$E-1$	1335.804	86.438	146.983	362.300
P2	-2.356$E+1$	1.909$E+2$	1.944$E-1$	9.375$E-1$	5.934$E-2$	9.938$E-2$	2.465$E-1$	1457.094	86.459	144.808	359.174
P3	-2.690$E+1$	2.254$E+2$	1.298$E-1$	8.164$E-1$	5.182$E-2$	8.691$E-2$	2.113$E-1$	1673.228	86.701	145.425	353.486
P4	-2.922$E+1$	2.484$E+2$	1.038$E-1$	7.632$E-1$	4.849$E-2$	8.033$E-2$	1.953$E-1$	1789.874	86.784	143.775	349.481
P5	-3.635$E+1$	3.167$E+2$	8.389$E-2$	7.130$E-1$	4.503$E-2$	7.493$E-2$	1.820$E-1$	1916.000	86.280	143.566	348.727
P6	-3.990$E+1$	3.528$E+2$	6.458$E-2$	6.475$E-1$	4.125$E-2$	6.869$E-2$	1.640$E-1$	2109.794	87.026	144.931	345.914
P7	-5.259$E+1$	4.746$E+2$	5.463$E-2$	6.305$E-1$	3.870$E-2$	6.517$E-2$	1.533$E-1$	2166.640	84	141	332
P8	-6.699$E+1$	6.122$E+2$	5.152$E-2$	5.952$E-1$	3.753$E-2$	6.417$E-2$	1.497$E-1$	2295.396	86	147	344
P9	-1.024$E+2$	9.534$E+2$	4.431$E-2$	5.742$E-1$	3.392$E-2$	6.257$E-2$	1.446$E-1$	2379.292	85	149	344
P10	-1.474$E+2$	1.387$E+3$	4.203$E-2$	5.637$E-1$	3.532$E-2$	6.095$E-2$	1.558$E-1$	2423.379	86	148	377

Equations

$$S_{\mathrm{C}}=$$

$$\int d^{4}x\left\{\frac{1}{2}{\partial }_{\mu }{\vec{\varphi }}_0{\partial }_{\mu }{\vec{\varphi }}_0+\frac{m_0^2}{2}{\vec{\varphi }}_0{\vec{\varphi }}_0+\frac{g_0}{24}{\left({\vec{\varphi }}_0{\vec{\varphi }}_0\right)}^2-{\vec{H}}_0{\vec{\varphi }}_0\right\},$$

$$S_{\mathrm{L}}=$$

$$a{\vec{\varphi }}_0$$

$$g_0$$

$$a^3{\vec{H}}_0$$

$$a^2{m}_0^2$$

$$\frac{m_{\sigma }a}{{\bar{\varphi }}_{R}a}$$

$$\frac{m_{\pi }a}{{\bar{\varphi }}_{R}a}$$

$$\vec{s}(t)=\frac{1}{N_S^3}\sum _{\vec{x}}\vec{\phi }(t,\vec{x}),$$

$$C_{ij}(t)=\frac{1}{N_T}\sum _{\tau }{s}_i(t)s_j(t-\tau )$$

$$f(t)=A+B\left(\exp (-mt)+\exp (-m(N_T-t))\right),$$

$$C_{\pi }(t)=\frac{1}{3}\sum _{i=2}^4{C}_{ii}(t).$$

$$⟨C_{\sigma ,\mathrm{c}}(t)⟩=⟨C_{\sigma }(t)⟩-⟨M_1{⟩}^2,$$

$$\vec{M}=\frac{1}{N_T{N}_S^3}\sum _x\vec{\phi }(x).$$

$$⟨C_{\sigma ,\mathrm{c}}(t)⟩=⟨C_{\sigma }(t)-M_1^2⟩.$$

$$⟨M_1⟩=\frac{a{\bar{\varphi }}_0}{\sqrt{2\kappa }},$$

$${\bar{\varphi }}_R=\frac{{\bar{\varphi }}_0}{\sqrt{Z}}.$$

$$G_{R,\pi }^{-1}(p=0)\stackrel{!}{=}{m}_{\pi }^2.$$

$$G_{R,\pi }^{-1}(p=0)=\frac{H_R}{{\bar{\varphi }}_R}=Z\frac{H_0}{{\bar{\varphi }}_0},$$

$$Z=2\kappa m_{\pi }^2⟨M_1⟩h^{-1}.$$

$$a{\bar{\varphi }}_R=a_{\mathrm{phys}}\times 93\mathrm{MeV}=a_{\mathrm{phys}}\times (93/197.327){\mathrm{fm}}^{-1}.$$

$$g_0^{-1}(a)=g_1^{-1}-({\beta }_1+{\beta }_2{g}_1)\log (a_1/a),$$

$${\beta }_1=\frac{1}{3}\frac{N+8}{8{\pi }^2},{\beta }_2=-\frac{1}{3}\frac{3N+14}{(16{\pi }^2{)}^2}$$

$$g_R=\frac{3m_R^2}{v_R^2},$$

$$\displaystyle\log\big{(}m_{R}a_{\rm min}\big{)}=\frac{1}{\beta_{1}g_{R}}+\frac{\beta_{2}}{\beta_{1}^{2}}\log(\beta_{1}g_{R})-1.9(1)\,,$$

$$-\frac{m_{0,c}^2}{{\Lambda }^2}=\frac{N+2}{6}\left[\frac{g_0}{(4\pi {)}^2}-\frac{2}{3}\frac{g_0^2}{(4\pi {)}^4}+\mathcal{O}(g_0^4)\right]$$

$$\displaystyle\bar{M}^{2}_{\phi=0}(K)=m_{0}^{2}+\frac{N+2}{6}\left[g_{0}{\cal T}[\bar{G}]-\frac{g^{2}_{0}}{3}{\cal S}[\bar{G}](K)\right].\ \$$

$${\tilde{M}}^2$$

$$\frac{m_{0,c}^2}{{\Lambda }^2}$$

$${\partial }_k{\Gamma }_k[\varphi ]=\frac{1}{2}\mathrm{Tr}\left[{\partial }_k{R}_k{\left(\frac{{\delta }^2{\Gamma }_k[\varphi ]}{\delta {\varphi }_i\delta {\varphi }_j}+R_k{\delta }_{ij}\right)}^{-1}\right],$$

$$\displaystyle\Gamma_{k}[\phi]=\int d^{d}x\left(\frac{1}{2}\big{(}\partial_{x}\phi_{i}\big{)}^{2}+U_{k}(\rho)\right),$$

$$\displaystyle\partial_{k}U_{k}(\rho)=\frac{k^{5}}{32\pi^{2}}\left(\frac{N-1}{k^{2}+\hat{M}^{2}_{\rm T}(k)}+\frac{1}{k^{2}+\hat{M}^{2}_{\rm L}(k)}\right),\$$

Full text

I Introduction

The $\phi^{4}$ scalar model with an internal $O(4)$ symmetry in 4 space-time dimensions has been long used as a model for spontaneous chiral symmetry breaking Lee72 . The direction of the symmetry breaking is associated with the sigma meson, while the pions are the Goldstone bosons emerging as a result of the spontaneous symmetry breaking. It is also widely known, that as a field theory it is trivial, it has no finite ultraviolet limit with non-zero coupling strength Luscher:1988uq . Although this property is still discussed (see e.g. Shrock:2017zuk ) we accept it as a fact and investigate what is the bound set by triviality to the quantitative applicability of the model. Based on a calculation carried out by Lüscher and Weisz (LW) in the same model applied to the Higgs particle Luscher:1988gc one can estimate the lowest lattice spacing that can be reached in a parametrization adjusted to light mesons. This turns out to be $a_{\rm min}^{\rm LW}=0.40(4)$ fm, which corresponds to a maximal cutoff in momentum representation to a few times $500$ MeV. This foreshadows that a scaling region of physical quantities as a function of $a$ on the lattice is unlikely to be found without getting too close to the triviality bound, and therefore cutoff independence, even in the effective theory sense is not feasible.

The above estimate was derived in a specific renormalization scheme for the case without explicit symmetry breaking. It is interesting to see to what extent it changes when, compared to Luscher:1988gc , a different renormalization scheme is employed in the case when the pions are massive. At the same time, experience shows Marko:2013lxa that the use of continuum functional methods is less restricted in the shadow of triviality and can retain some predictivity. To study this in more detail we use two continuum methods: the functional renormalization group (FRG) Wetterich:1992yh ; Berges:2000ew in the local potential approximation (LPA) and the two-loop and ${\cal O}(g_{0}^{2})$ truncations of the two-particle irreducible approach (2PI) Cornwall:1974vz ; Berges:2005hc ; Marko:2013lxa . Treating the model as a cutoff theory, we solve it using the same bare couplings as in the lattice version along the line of constant physics (LCP). Then, to compare the values of physical quantities, we need the relation between the lattice spacing $a$ and the cutoff $\Lambda$. This is determined by matching the critical line of the model at zero temperature with the one determined by Lüscher and Weisz in Luscher:1988uq using the hopping parameter expansion.

The paper is structured as follows. In Sec. II we introduce notations for the model and summarize the details of the lattice simulations. In Sec. III we define the LCP and describe how the triviality bound is obtained. We also discuss the immediate consequences of the value of the minimal lattice spacing. In Sec. IV we compare the lattice results with those obtained in the continuum approximations, and finally in Sec. V we summarize our findings.

II Generalities

We discuss the $O(N)$ symmetric, Euclidean $\phi^{4}$ model specifically for $N=4$, described in terms of bare quantities denoted by the subscript [math] by the continuum action (omitting the obvious $x\equiv(t,\vec{x})$ dependencies)

[TABLE]

where $\vec{\phi}_{0}$ is the $N=4$ component field, $m_{0}$ is the mass, and $g_{0}$ is the quartic coupling. In the explicit symmetry breaking term the external field $\vec{H}_{0}$ is chosen to point in the direction of the first component of the scalar field with a length of $H_{0}$, independent of $x$. Discretization on a periodic, four-dimensional cubic lattice consisting of $N_{T}\times N_{S}^{3}$ sites (using a forward derivative), and rewriting in terms of the hopping parameter $\kappa$ leads to the well known lattice action Luscher:1988uq :

[TABLE]

where $a$ is the lattice spacing and $\hat{\mu}$ is the usual four-dimensional unit vector. The connections between the continuum and the lattice parameters are

[TABLE]

We use Monte Carlo integration with importance sampling to evaluate path integrals. Configuration generation is done by using a poor man’s heat bath algorithm, in which each site is updated using ten metropolis steps before its neighbors are updated in order to make the new field value at the chosen site practically independent of its initial value. Between two heat bath sweeps we also include two overrelaxation sweeps in order to sample a much larger part of the phase space using the same number of configurations.

III Line of constant physics

III.1 Observables defining the LCP

The explicitly broken $O(4)$ symmetric $\phi^{4}$ model has three parameters: the hopping parameter $\kappa$, the quartic coupling $\lambda$ and the external field $h$. In order to define a continuum limit111Even though a true continuum limit is not possible due to the triviality of the model, we follow the standard procedure which would allow to define it if it existed. we give two physical prescriptions, which restrict our parameter space to the LCP, along which the lattice spacing $a$ tends to zero in physical units, at least in principle. The two prescriptions are

[TABLE]

where $m_{\sigma,\pi}$ are the respective pole masses and $\bar{\phi}_{R}$ is the expectation value (denoted by the bar) of the $\sigma$ component of the renormalized field, which takes the role of the pion decay constant in the linear sigma model (LSM). We choose a lower sigma mass (300 MeV) than what is generally agreed upon ($\approx 450$ MeV) Pelaez:2015qba . Our choice is limited on the one hand by the fact that higher sigma masses are barely reachable in approximate continuum solutions of the LSM Patkos:2002xb ; Andersen:2004ae and on the other hand by the fact that we want to retain the kinematic possibility of the $\sigma\to 2\pi$ decay.

To obtain the pole masses we measure time slice correlators. Let us define a time slice as

[TABLE]

and then

[TABLE]

is the time slice correlator matrix for one configuration. The ensemble average of $C_{ij}(t)$ is the time slice correlator. By our choice of $\vec{h}$ the $\sigma$ direction is $i=1$; therefore $C_{\sigma}(t)\equiv C_{11}(t)$ is dominated by $m_{\sigma}$, while $C_{ii}(t)\,,i\neq 1,$ are all dominated by $m_{\pi}$. We do a least squares fit using the function222Contributions of excited states, if they exist, are invisible within our precision: the remnant obtained by subtracting the fitted form from the data is consistent with zero.

[TABLE]

with parameters $A\,,B$ and $m$ to $C_{\sigma}(t)$ as well as to the average of the three pion directions333Averaging over the three pion directions lowers the statistical error of $m_{\pi}$ compared to $m_{\sigma}$.

[TABLE]

The average and error of the fit parameters and in particular the masses are obtained by a jackknife analysis. The fit is carried out on each jackknife sample, leaving out the $t=0$ point of the correlator from the data in order to lower the distortions caused by possible higher excitations.

In the case of the sigma mass, one must take care of the disconnected part of the correlator. The connected part of the correlator is

[TABLE]

where $M_{1}$ is the first component of the average field over one configuration

[TABLE]

To subtract the correlated errors from the connected sigma correlator, instead of (12) we use another prescription (the two definitions differ only in a constant),

[TABLE]

The definition (12) has a bad signal to noise ratio due to correlated errors which are canceled in (14) leading to a better signal. We show the reduction of error achieved by using the definition in (14) in Fig. 1.

The measurement of $\bar{\phi}_{R}$ goes as follows. The ensemble average of the first component (the sigma direction) of $\vec{M}$ is

[TABLE]

where the [math] index on the right-hand side denotes that $\phi_{0}$ is a bare fields; that is, wave function renormalization is still needed. Then the renormalized vacuum expectation value is

[TABLE]

We obtain $Z$ by prescribing the value of the zero-momentum inverse pion propagator to be the pion pole mass:

[TABLE]

Through a Ward identity ZinnJustin:2002ru the inverse two-point function can be rewritten as

[TABLE]

which, in terms of lattice quantities and combined with the renormalization prescription, leads to

[TABLE]

The value of $Z$ is slowly changing between 0.74 and 0.8 along the LCP in the measured range of $a$.

III.2 Determining the LCP

To obtain a point of the LCP curve, we fix the value of one of the parameters (usually $h$, but in the region where triviality strongly influences the LCP we fix $\lambda$) and measure the ratios of observables appearing in (7a) and (7b) on an appropriate grid in the plane of the remaining two parameters, $\kappa-\lambda$ and $\kappa-h$ planes respectively. The physical values of ratios define contour lines and the LCP is obtained as the intersection of two contour lines, each belonging to one surface. Points of the contour lines are obtained by linear interpolation between grid points and are fitted with parabolas. The intersection of the two parabolas is one point of the LCP corresponding to the $h$ or $\lambda$ where the grid was defined. The error is estimated by a bootstrap resampling using $10^{4}$ samples. A detailed description of this procedure is relegated to Appendix A. The original ratios of observables were obtained using $16\times 16^{3}$ lattices with $10^{5}$ field configurations. While $16\times m_{\pi}a>6$ even for the smallest lattice spacing, we tested volume independence by checking at the last LCP point that the observable ratios of (7a) and (7b) only change within errors when $20\times 20^{3}$ and $24\times 24^{3}$ lattices are used.

With the method outlined above, we obtain the points of the LCP shown in Fig. 2. The conversion of $a$ to physical units is done using

[TABLE]

In Fig. 2 we see that the LCP follows the critical line in the $\lambda-\kappa$ plane. In a theory with a proper continuum limit the LCP should run into the critical line at least at infinite coupling. While in our case $h\neq 0$ for all LCP points, $h$ tends to zero along the line as it should to approach criticality. Triviality appears here by seeing that even at $\lambda\to\infty$ the LCP does not converge to the critical line, meaning that $a$ remains finite. This means that the bare $\phi^{4}$ coupling $g_{0}$ must have a pole as a function of $a$, at the minimal value of the lattice spacing. The results for $g_{0}$ are shown in Fig. 3 and are in compliance with the generally accepted view on the triviality of the $\phi^{4}$ model. Fitting the data shown in Fig. 3 with the second order perturbative $\beta$ function

[TABLE]

where

[TABLE]

are the standard $\beta$-function coefficients Smit:2002ug , we estimate $a_{\rm min}=0.52(2)$ fm. This leads to an estimate for the minimal value of the lattice $\sigma$ mass, $am_{\sigma}=0.79(3)$.

The result for the minimal lattice spacing can be compared to the one which can be given based on Luscher:1988gc . In the renormalization scheme of Lüscher and Weisz

[TABLE]

where $g_{R}$ is the renormalized quartic coupling, $m_{R}$ is the renormalized mass, which we identify with the sigma mass for the sake of the estimate, and $v_{R}\equiv\bar{\phi}_{R}$ takes the value of the pion decay constant as in our case, although the $Z$ factor, which we do not need here, is defined differently. In Luscher:1988gc the renormalization trajectories are described, and taking the $\lambda\to\infty$ limit in them yields a relation between $m_{R}a_{\rm min}$ and $g_{R}$:

[TABLE]

where the number $1.9(1)$ is the result of a numerical calculation at a high order of the hopping parameter expansion. Plugging $m_{R}=300$ MeV and $v_{R}=93$ MeV into (23) and (24) yields $a_{\rm min}^{\rm LW}=0.4$ fm already mentioned in Sec. I. We see that our result $a_{\rm min}=0.52(3)$ is even more restrictive.

An important implication of the largeness of $a_{\rm min}$ is that on the lattice the maximal temperature that can be simulated is $T=(N_{t}\cdot a_{\rm min})^{-1}|_{N_{t}=1}\approx 420$ MeV. Furthermore if one is interested in a “continuum limit” in the effective theory sense the feasible temperature range is definitely below $50$ MeV. This limits the comparison of the continuum methods practically to vacuum quantities.

IV Comparison with 2PI and FRG results

Since according to the previous section, a comparison between lattice and continuum physical quantities is not feasible at a finite temperature, we remain at $T=0$ and using continuum functional methods we determine the masses along the LCP shown in Fig. 2. In order to compare to a lattice result determined at a fixed lattice spacing, we need to treat the continuum version of the model as a cutoff theory. Hence we need the relation between the continuum cutoff $\Lambda$ and the lattice spacing $a$, i.e. we need $c=\Lambda a$. This relation was studied in Brahm:1994zv ; Besjes where the conversion factor $c\approx 4.9$ was calculated analytically for the four-dimensional (4D) hypercubic lattice and obtained also by fitting the perturbative continuum result (using (4) and (6))

[TABLE]

to the critical line $m^{2}_{0,c}(g_{0})$ obtained by Lüscher and Weisz in Luscher:1988uq . The above equation comes from the condition of the vanishing curvature mass at the vanishing field value at second order in the perturbation theory. We can see in Fig. 4 that at $\mathcal{O}(g_{0}^{2})$ it reproduces the LW critical line only at small values of the coupling. One expects that this behavior changes if one uses a more sophisticated approximation.

In Fig. 4 we also show the LCP points corresponding to the “fixed h” data points of Fig. 2. While doing the transformation of the bootstrap data from the $\lambda-\kappa$ plane to the $g_{0}-m_{0}^{2}$ plane, we have also performed a principal component analysis of the correlation matrix. As indicated in Fig. 4, the bootstrap sample in the $g_{0}-m_{0}^{2}$ plane is much more elongated in one direction than it was in the $\lambda-\kappa$ plane. Based on the standard deviations of $g_{0}$ and $m_{0}^{2}$ over the bootstrap sample, shown for the rightmost LCP point, one could be incorrectly led to think that the LCP is compatible with the critical line. This is not the case, of course, as the bootstrap sample is well separated from the critical line, as was the case already in the $\lambda-\kappa$ plane.

IV.1 The critical line $m^{2}_{0,c}(g_{0})$ in the 2PI framework

In the 2PI framework, the curvature mass at the vanishing field value is given at the $\mathcal{O}(g_{0}^{2})$ level of truncation of the effective action by $\hat{M}^{2}_{\phi=0}=\bar{M}^{2}_{\phi=0}(K=0)$, where the gap mass satisfies the self-consistent equation Marko:2015gpa

[TABLE]

The tadpole and setting-sun integrals involve the propagator $\bar{G}(K)=1/(K^{2}+\bar{M}^{2}_{\phi=0}(K)).$ The critical line is determined from the condition of the vanishing curvature mass: $\hat{M}^{2}_{\phi=0}=\bar{M}^{2}_{\phi=0}(K=0)=0$.

The nontrivial momentum dependence makes (26) rather hard to solve, however, the solution $\bar{M}^{2}_{\phi=0}(K=0)$ can be approximated by using a localized propagator with momentum independent mass gap $\bar{M}^{2}=m_{0}^{2}+\frac{N+2}{6}g_{0}{\cal T}[\bar{G}_{\rm loc}]$, where $\bar{G}_{\rm loc}(K)=1/(K^{2}+\bar{M}^{2}).$ This approximation corresponds in fact to the two-loop 2PI truncation. In this approximation, the tadpole can be explicitly computed with a 4D cutoff $\Lambda$ and the condition of the vanishing curvature mass can be written as

[TABLE]

where $\tilde{\cal S}(\tilde{M}^{2})$ is a perturbative setting-sun integral at vanishing external momentum and we used the tilde for a quantity scaled by appropriate powers of $\Lambda$.

For a given $g_{0}$, one then solves (27a) for $\tilde{M}^{2}$ and using this solution one has $m_{0,c}^{2}$ from (27b). The critical line obtained in this way is shown in Fig. 4. It still deviates from the LW curve, but remains closer to it in a wider range of the coupling than the $\mathcal{O}(g_{0}^{2})$ perturbative curve. We mention that (25) can be obtained by first using (27a) in the first term on the right-hand side of (27b) and then taking the $\tilde{M}^{2}\to 0$ limit, in which $\tilde{\cal S}(\tilde{M}^{2})\to 2/(4\pi)^{4}$.

Now let us discuss the determination of the critical line by solving (26) without further approximation. $m^{2}_{0,c}(g_{0})$ could be obtained in principle by approaching it from the symmetric phase: fixing $g_{0}$, the equation is solved for increasing values of $|m_{0}^{2}|$ and $m^{2}_{0,c}(g_{0})$ is obtained by extrapolating the determined values of $\bar{M}^{2}_{\phi=0}(K=0)$ to zero. As detailed in Appendix C, Eq. (26) is solved by treating the setting-sun ${\cal S}(K)$ as a double convolution: a convolution of the propagator with a bubble integral, where the latter is itself a convolution of two propagators. It turned out that the solution to (26) is lost for a value of $m^{2}_{0}$ where $\bar{M}^{2}_{\phi=0}(K=0)$ is nonzero (see Fig. 8). This loss of solution, which seems to be a feature of the $\mathcal{O}(g_{0}^{2})$ 2PI gap equation, and was investigated in detail in Marko:2015gpa , prevents us from the direct determination of the critical line at this order of the 2PI truncation scheme and furthermore from a comparison along the LCP.

As Fig. 4 shows, the simpler two-loop approximation indeed has a critical line determined by Eqs. (27a) and (27b). Nevertheless, a loss of solution can also happen in this approximation in the broken phase (that is at $\phi\neq 0$) depending on the parameters Marko:2015gpa . We found that the usual iterative procedure to solve the broken phase two-loop equations (which were written down and solved as detailed in Marko:2013lxa with little modifications to accommodate for the use of nonrenormalized equations and approximating the $T\to 0$ limit numerically) breaks down close to the critical line in comparison to where the points of the LCP are; therefore in the LCP points no solution exists and no comparison can be made. We checked that this loss of solution persists in the even simpler localized two-loop approximation which we detailed in Marko:2015gpa . We conclude that in the considered approximations the 2PI formalism cannot be compared to the lattice LCP results.

IV.2 Determination of observables using the FRG method

Another functional method from which one can calculate curvature masses along the LCP is the functional renormalization group method. The flow equation describing the evolution of the scale-dependent average action $\Gamma_{k}$ from the ultraviolet (UV) scale $k=\Lambda$, where the microscopic theory is defined through the bare action, down to the deep infrared (IR), where the usual quantum effective action is obtained in the $k\to 0$ limit, is Wetterich:1992yh

[TABLE]

where $R_{k}$ is a regulator function, that is, in momentum space it suppresses the IR modes, while $\partial_{k}R_{k}$ regulates the integral in the UV. In the LPA the Ansatz

[TABLE]

is used, where $\rho=\vec{\phi}^{2}/2$ is $O(N)$ invariant and it is customary to choose the LPA-optimized regulator Litim:2001up $R_{k}(q)=(k^{2}-q^{2})\Theta(k^{2}-q^{2})$ ($q$ is the Euclidean four-momentum). Then, using $\displaystyle\frac{\partial^{2}U_{k}(\rho)}{\partial\phi_{i}\partial\phi_{j}}=U_{k}^{\prime}(\rho)(P^{\rm L}_{ij}+P^{\rm T}_{ij})+2\rho U_{k}^{\prime\prime}(\rho)P^{\rm L}_{ij}$ with $P^{\rm L/T}$ being the longitudinal/transverse projectors, the integral can be performed and, at zero temperature and $d=4$, one obtains

[TABLE]

where $\hat{M}^{2}_{\rm T}(k)=U_{k}^{\prime}(\rho)$ and $\hat{M}^{2}_{\rm L}(k)=M^{2}_{\rm T}(k)+2\rho U_{k}^{\prime\prime}(\rho).$ This equation is solved numerically by integrating it down to $k=0$ (in practice to some $k_{\rm end}>0$, due to the flattening of the potential) starting at scale $k=\Lambda,$ where the initial condition for the potential is given in terms of the couplings $m_{0}^{2}$ and $g_{0}$ as $U_{k=\Lambda}(\rho)=m_{0}^{2}\rho+g_{0}\rho^{2}/6.$

In the so-called grid method $U_{k}(\rho)$ is discretized using $N_{\rho}$ grid points so that (30) transforms into a system of $N_{\rho}$ coupled ordinary differential equations. We solve this system using the Runge-Kutta-Fehlberg algorithm with adaptive step-size control provided by the GNU Scientific Library (GSL)GSL . We work in units of the cutoff, denoting with tilde a quantity scaled with the cutoff, and we choose $N_{\rho}=5000$ equidistant values of $\tilde{\rho}=\rho/\Lambda$ in the range between 0 to $\tilde{\rho}_{\rm max}=0.026$. The flow was stopped at $\tilde{k}_{\rm end}=1.28\times 10^{-2}$ where all the monitored quantities became practically constants. At each point of the grid the first and second order derivatives of the potential are calculated with ${\cal O}(\Delta\tilde{\rho}^{4})$ finite difference formulas. The minimum of the potential is obtained with spline interpolation, while the curvature masses at the minimum are obtained fitting a sixth order polynomial to the potential in an appropriate $\tilde{\rho}$ interval which has the minimum as its left end point.

The transverse and longitudinal curvature masses are obtained at $k_{\rm end}$ as $\hat{M}_{\rm T}=\sqrt{U_{k}^{\prime}(\bar{\rho})}$ and $\hat{M}_{\rm L}=(\hat{M}^{2}_{\rm T}+2\bar{\rho}U_{k}^{\prime\prime}(\bar{\rho}))^{1/2}$, where $\bar{\rho}=\bar{\phi}^{2}/2$ is the minimum of the potential. In the LPA they can be regarded as approximations to the pole masses due to the simplicity of the Euclidean propagator. In order to compare $\hat{M}_{\rm T/L}$ along the LCP with the values $M_{\sigma}=300$ MeV and $M_{\pi}=138$ MeV which in the lattice simulation are constant along the line, we need to know what is the relation between the cutoff scale $\Lambda$ and the lattice spacing $a$. This relation is obtained by matching the critical curve $\tilde{m}_{0,c}^{2}(g_{0})$ determined in the FRG case to the one obtained by Lüscher and Weisz in Luscher:1988uq using the hopping parameter expansion to 14th order. We determine $\tilde{m}_{0,c}^{2}(g_{0})$ working at fixed $g_{0}$ and using dichotomy on $\tilde{m}_{0}^{2}$, as shown in Fig. 5, where the quantity that distinguishes between the broken and symmetric phases is $U^{\prime}_{k}(0)/k^{2}$.

Once we match $\tilde{m}_{0,c}^{2}$ and $m_{0,c}^{2}a^{2}$ at some value of $g_{0}$, finding the relation

[TABLE]

the entire critical curve determined using FRG agrees with the one obtained by Lüscher and Weisz, as shown in Fig. 4. The very good agreement of the two critical curves is in line with the findings of Ref. Caillol:2012dt , where it was reported that in the one component $\phi^{4}$ model the critical line, obtained in the LPA with the Litim regulator and lattice discretization, compares well with the one determined with Monte Carlo simulations.

Having obtained the relation between $a$ and $\Lambda$, we can now solve the flow equation (30) and determine the curvature masses for the fixed $h$ data points of Fig. 2. The results shown in the first four rows of Table 1 in units of the cutoff can be used in two ways. In the first case, shown in the last four columns of Table 1, one can determine for each point of the LCP the value of the cutoff from the lattice spacing using (31). Then $\bar{\phi}$ is smaller than $f_{\pi}=93$ MeV by $\sim 8\%,$ $\hat{M}_{\rm T}$ is $5\%-8\%$ larger than $M_{\pi}=138$ MeV, while $\hat{M}_{\rm L}$ is $15\%-20\%$ larger than $M_{\sigma}=300$ MeV. In the second case one can require $\bar{\phi}$ to be $f_{\pi}$. In this case, due to the larger value of the cutoff, one finds that $\hat{M}_{\rm L}$ is $22\%-30\%$ larger than the sigma values used to determine the LCP, while $\hat{M}_{\rm T}$ is larger by around $10\%$ than the pion mass. The deviation from the lattice results decreases for smaller $a$.

IV.2.1 Solution of a modified flow equation for $m_{0}^{2}<-\Lambda^{2}<0$

With the chosen quartic potential at the initial value of the scale, $k=\Lambda,$ the flow equation (30) cannot be solved for $m_{0}^{2}<-\Lambda^{2}<0$ due to a singularity in the equation. One could either change the initial condition by including higher order, perturbatively nonrenormalizable terms in the potential or, as we do it here following Barnafoldi:2016tkd , try to circumvent the problem by modifying the flow equation expanding in power series to some order $N_{g}$ the fractions appearing in the right-hand side of (30),

[TABLE]

where $\xi=(M_{0}^{2}-\hat{M}^{2}_{\rm L/T}(k))/(k^{2}+M_{0}^{2})$ with $M_{0}$ some large parameter, i.e. $M_{0}>\Lambda$, which for numerical reasons has to be chosen appropriately.444We use the same values of $M_{0}$ and $N_{g}$ for both fractions in (30).

First, keeping the numerical framework used so far, that is changing only the right-hand side of (30) according to (32), we tested the method in a case where a direct solution to (30) exists and then we applied it for the fixed $\lambda$ data points of the LCP shown in Fig. 2 (points P7-P10 in Table 1). In the latter case the solution is regarded as an approximation to the solution of the original Wetterich equation (28), assumed to exist for an appropriate form of the effective action at scale $\Lambda$.

In the case of point P1, it turns out that in order to reproduce the available direct solution of (30) with the expansion method, one has to go to rather high orders in the expansion. Also, for the method to work, the first and second derivatives of the potential at $\tilde{\rho}_{\rm max}$ had to be kept fixed as a function of $k$, however, the chosen values were practically arbitrary. We fixed the derivatives to their values calculated at $k=\Lambda$.

At a given order of the expansion the deviation from the direct result increases with $M_{0}.$ Among the studied quantities, $\hat{M}_{\rm L}$, presented in Fig. 6, shows the slowest convergence rate with $N_{g}$ at a fixed value of $M_{0}$. For $\tilde{M}_{0}^{2}=2$ and $N_{g}=50$ the deviation from the direct result is around $10\%$. To estimate the result of the curvature masses and the minimum of the potential we fitted with $f(x)=a+b/(x-c)^{d}$ the data obtained at various $N_{g}$ with the expansion method. For P1 one can practically recover the direct results from a dataset obtained with up to $N_{g}\simeq 100$ terms in the expansion, but as $|m_{0}^{2}|/\Lambda^{2}$ increases we need larger $M_{0}^{2}$ and larger values of $N_{g}$ to maintain the quality of the fit. Eventually, numerical errors prevent us from going above a certain value of $N_{g}$. All these features are illustrated in Fig. 6 and the results obtained with the expansion values are given in the last four rows of Table 1. Based on the variation of the extrapolated results on the fitting $N_{g}$-interval, one can estimate the error of $\hat{M}_{\rm L}$ to be $1\%-2\%$ for P7 and P8 and $5\%-10\%$ for P9 and P10. For the other two quantities the error of the extrapolation to $N_{g}=\infty$ is smaller.

V Conclusion

We studied the four component Euclidean $\phi^{4}$ model in four dimensions. In the presence of an explicit symmetry breaking term, we determined with Monte Carlo simulations the line of constant physics (LCP) in the bare parameter space of the model based on ratios involving the pion and sigma masses and the expectation value of the field. In this process we brought further evidence in support of the triviality of the model in a renormalization scheme which is different from the one usually used by the lattice community (see Siefert:2014ela for a recent study).

Using the bare couplings of the LCP, we solved the model with two continuum functional methods (the 2PI formalism and the FRG method) in an attempt to compare the vacuum masses and expectation value obtained with these continuum methods to the corresponding input values of the lattice study of the model. The manifestation of triviality prevented us from a meaningful comparison of finite temperature quantities. It turned out that the comparison at $T=0$ can be done only with the FRG, since the 2PI is hindered by the loss of solution to the propagator equation. The needed relation between the lattice spacing and the cutoff, used in the lattice and continuum versions of the model, respectively, was obtained by matching the critical line of the parameter space determined originally by Lüscher and Weisz using hopping parameter expansion.

Acknowledgements.

The authors would like to thank Sz. Borsányi, S. D. Katz, G.  Fejős, A. Jakovác and Cs. Török for helpful and interesting discussions. In the 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. acknowledges support by the János Bólyai Research Scholarship of the Hungarian Academy of Sciences.

Appendix A Finding the LCP point and its error

We detail here the procedure used to find the points of the LCP. The seemingly overly complex procedure summarized in Sec. III.2 is needed because the observable ratios (7a) and (7b) change very similarly over the three-dimensional parameter space. Therefore, the search for the LCP becomes the mathematical problem of finding the intersection of two noisy surfaces which are almost parallel.

Our generic approach to find a point of the LCP is to fix one parameter and scan the remaining two for the pair of parameters where the observable ratios have both the prescribed physical values. For points of the LCP that are far from the triviality bound, that is $a_{\rm latt}$ is relatively large, we fix $h$ and scan the $\kappa$-$\lambda$ plane. As the lattice spacing becomes smaller, we change to fixing $\lambda$ and scan the $\kappa$-$h$ plane.

We discuss in what follows the case when we fix $h$ and search for the intersection point of two contour lines in the $\kappa$-$\lambda$ plane. The two curves are marked out by (7a) and (7b), which are almost parallel locally, as mentioned earlier. We scan the plane by doing simulations in a set of $(\kappa,\lambda)$ pairs forming a square grid. In each point of the grid we independently measured the observable ratios ($R_{1}$ and $R_{2}$) with known errors ($\Delta_{1}$ and $\Delta_{2}$). To estimate the location of the LCP point in the plane we generate $10^{4}$ bootstrap configurations on the grid, where the value of both observable ratios at each grid point is randomly taken from independent normal distributions with means $R_{1,2}$ and standard deviations $\Delta_{1,2}$. Each bootstrap configuration can then be thought of as two surfaces $R_{1,2}(\kappa,\lambda)$ over the parameter plane sampled on the grid points.

For each bootstrap configuration we linearly interpolate the observable surfaces along the edges of the square grid. We locate the points which belong to the physical contour lines of the respective surfaces (physical values of $R_{1}$ and $R_{2}$). These two sets of points are then fitted by parabolas. The intersection of the two parabolas is one bootstrap realization of the LCP point. The bootstrap realizations outline the probability distribution of the LCP point in the $\kappa$-$\lambda$ plane, and hence its error can be estimated.

The procedure is illustrated in Fig. 7 for $h=0.1$ (P3 of Table 1), which shows the bootstrap realization of the LCP point and its estimated average value and error. The bands correspond to the averages and $1\sigma$ deviation of the parabolas approximating the physical contours. The distribution shows correlation between the $\kappa$ and $\lambda$ coordinates of the LCP point and it is elongated along one of the contour lines. This elongation is further enhanced when transformed to $m_{0}^{2}$ and $g_{0}$, where the distribution becomes practically one dimensional. As a result we show error bands in the direction of the major axis of the distribution for each LCP point in Fig. 4.

By appropriately choosing the grid and including more points in it, as well as increasing the precision of the lattice simulations, one could in principle further confine the location of the LCP points.

Appendix B The $T=0$ setting-sun integral as a Hankel transform

Using the Fourier transform, the convolution of two momentum-dependent functions can be written as

[TABLE]

where we used the shorthands $\int_{k}=\int\frac{d^{4}k}{(2\pi)^{4}}$ and $\int_{x}=\int d^{4}x$. Working with spherical coordinates in 4D, the angular integral can be performed analytically by exploiting the rotation invariance. Choosing $q$ to point in the fourth direction, such that $q\cdot x=QX\cos\theta_{1}$, where $Q=|q|$ and $X=|x|,$ and introducing $\tau=\sin\theta_{1},$ the nontrivial part of the angular integration gives a Bessel function of the first kind,

[TABLE]

such that (33) becomes

[TABLE]

where we used the Hankel transform of order 1 of the function $A(X)=F_{1}(X)F_{2}(X)/X$, with $F_{i}(X)=Xf_{i}(X)$.

A similar calculation shows that $F_{i}(X)=Xf_{i}(X)$ appearing in (35) can be given as the inverse Hankel transform (denoted in what follows by a tilde) of order 1 of $\tilde{F}_{i}(P)=P\tilde{f}_{i}(P)$ ($i=1,2$):

[TABLE]

So, in terms of Hankel transforms, the convolution (33) can be written as

[TABLE]

Then, writing the momentum dependent setting-sun as a convolution of a propagator and a bubble integral $B[G](q-k)=\int_{p}G(p)G(q-k-p)$, both the bubble integral and the setting-sun integral $S[G](q)=\int_{k}G(k)B(q-k)$ can be written in terms of Hankel transforms. The discrete version of the Hankel transform (DHT) and its inverse are implemented in the GSL package GSL . In the discretized case it is understood that all momenta and combinations of momenta are cut by the cutoff $\Lambda.$ Working in units of the cutoff and using $3\times 2^{10}$ sampling points the smallest momentum on the grid is $|\tilde{k}|_{\rm min}=3.96\times 10^{-4}.$

Appendix C On the solution of (26)

As mentioned in Sec. IV.1, the setting-sun ${\cal S}[G](K)$ can be regarded as a double convolution. The convolution integral can be computed using the Fourier transform, which for rotational invariant functions leads at $T=0$ to the use of the Hankel transform, as detailed in Appendix B. As discussed in Sec. V.E.1 of Marko:2012wc , calculating a convolution using discrete Fourier transform is not accurate if the function does not decrease fast enough in the UV. We expect this behavior in the case of the Hankel transform as well. Although to determine the critical line we are interested in an IR quantity, namely $\bar{M}^{2}_{\phi=0}(K=0),$ since a momentum integral is involved in its calculation, the discretization error in the UV will influence this quantity. To check the method that uses the Hankel transform, we also computed on a nonuniform momentum grid 555We use $N_{K}=128-256$ values of momenta: $K_{i}=K_{\rm min}+(\Lambda-K_{\rm min})\,(i/(N_{K}-1))^{2.5}$, $i=0,\dots,N_{K}-1$ with $\tilde{K}_{\rm min}=3.96\times 10^{-4}$. the convolution as a double integral using Eq. (A1) of Fejos:2011zq . The inner integral in that expression is calculated numerically after a Tanh-Sinh transformation (see Appendix A of Reinosa:2011cs for details) using the splined $\bar{M}^{2}(K)$ in the propagator. The bubble integral is calculated in this way on a grid, then splined and used for the calculation of ${\cal S}(K).$

The iterative solution of (26) obtained for $g_{0}=150$ and $m_{0}^{2}/\Lambda^{2}=-0.5275$ using underrelaxation method Berges:2004hn with parameter $\alpha=0.1$ is shown in Fig. 8. The upper part of the figure shows what happens if the solution obtained with DHT is used as an initial propagator in the solver that computes the convolutions using adaptive integration routines on a grid with 256 momentum values. We see that $\bar{M}^{2}_{\phi=0}(K)$ obtained in the first iteration deviates by $5\%-8\%$ from the solution obtained with DHT as a result of the fact that, as anticipated, the setting-sun calculated with DHT is not accurate. As the iteration progresses, $\bar{M}^{2}_{\phi=0}(K)$ departs even more from the used initial function and hence the converged solution is substantially different than the one obtained with DHT.

In the lower part of Fig. 8 we show $\bar{M}^{2}_{\phi=0}(K_{\rm min})$ as function of $|m_{0}^{2}|$ at four values of $g_{0}$. The difference between the curves obtained with the two ways of treating the convolution increases with the value of the coupling. This is due to the fact that the numerical error made in computing the convolution with DHT is magnified when the setting-sun is multiplied with a larger coupling. More importantly, the shape of the curves is compatible with the fact that the solution of (26) is lost at some value of $m_{0}^{2}$ where $\bar{M}^{2}_{\phi=0}(K_{\rm min})$ is still finite. As a result the critical line cannot be determined.