1 Introduction

The discovery of a Higgs boson at the LHC [1, 2] completed the search for the building blocks of the standard model of particle physics. While the mass of the Higgs boson is in principle a free parameter of the standard model, it has long been known that it is not necessarily an arbitrary parameter. For instance, assuming that the description of fundamental physics in terms of standard-model degrees of freedom is valid at a high energy scale \(\Lambda \) and that the theory is sufficiently weakly coupled, the range of possible Higgs boson masses is restricted by a finite range, the so-called infrared (IR) window [39]. The fact that this IR window shrinks for increasing \(\Lambda \) can be traced back to the fixed-point structure of the RG flow [6, 10] which connects the mass of the Higgs boson to that of the heaviest quark.Footnote 1

The edges of the IR window [7, 1239], i.e., the upper and lower admissible values, for the Higgs mass are actually not sharply defined, but depend on a number of additional assumptions. This is most obvious for the upper “triviality bound”, as the Higgs sector becomes strongly coupled at high scales for large values of the Higgs mass. Perturbative estimates of this bound, e.g., depend on an ad hoc choice of coupling value up to which perturbation theory is trusted. Non-perturbative methods have shown that this upper bound relaxes considerably if one allows the system to start microscopically with a strong Higgs self-coupling [40, 41].

As is less well appreciated, a similar fuzziness also holds for the lower edge of the IR window on which we concentrate in the present work. This indeterminacy arises from the (implicit) assumptions imposed on the precise form of the microscopic theory at the high scale \(\Lambda \). For instance, by considering only renormalizable operators in conventional perturbative estimates of the lower bound, the couplings of all higher-order operators are implicitly chosen to vanish at the high scale \(\Lambda \). Strictly speaking, this corresponds to fixing infinitely many further parameters, in addition to the parameters of the standard model.Footnote 2

By contrast, if one defines the standard model more agnostically in terms of its symmetries, field content and measured IR parameters, the microscopic interactions quantified by the bare action at the high scale \(\Lambda \) remains largely unspecified. It can contain infinitely many operators and corresponding coupling constants which would be computable if the underlying theory was known. It is reasonable to assume that these couplings if measured in units of the scale \(\Lambda \) are of order \(\mathcal {O}(1)\). Still, an even more strongly coupled UV regime is not excluded by observation. The fact that the perturbative description of electroweak collider data works so well indicates that colliders so far probe a regime where Nature is close to the Gaußian fixed point of the renormalization group (RG). Here, power-counting arguments hold and higher-dimensional operators are irrelevant, exerting a negligible influence on IR observables. In turn, the IR observables dominantly constrain only the marginal and relevant (renormalizable) operators of the bare action and put hardly any bound on the irrelevant couplings.

In a series of works, it has recently been shown that these unconstrained higher-dimensional operators in fact can relax the lower edge of the IR window, i.e., can lower the lower (stability) bound on the Higgs mass without introducing metastability [4042]. Comparatively simple modifications of the bare action, e.g., in terms of a dimension-six operator at the Planck scale can lower the lower mass bound by \(\sim \)1 GeV, while preserving absolute stability on all scales [42]. The mechanism behind this relaxation of the lower edge has two aspects: first, negative couplings of renormalizable operators which seem to introduce an instability from a perturbative viewpoint can still be associated with a fully stable potential in the presence of the higher-dimensional operators. Second, the higher-dimensional operators take influence on the running of the renormalizable part over a range of scales thus modifying the approach to the perturbative region while leaving the IR region itself intact as addressed by perturbation theory.

These results can already be observed in a simple mean-field study (large-N limit) as well as in extended mean-field approximation (1/N corrections), but unfold more comprehensively in a controlled non-perturbative study using the functional RG. They have been confirmed by lattice simulations in the range of scales accessible to current lattice sizes [4346]. Recent functional RG studies also including higher-order fermionic operators have shown the same features [4749]; for further studies of higher-dimensional operators in this context, see, e.g., [5055].

While controlled quantitative results have so far been obtained for a small class of operators represented by simple low-order polynomials of the field, a possible metastable regime with competing vacua has not been explored so far.

The present work is devoted to a first step in this direction, namely to study the full functional renormalization of the Higgs effective potential as a function of both field amplitude and RG scale. This is facilitated by the development of pseudo-spectral methods for functional flows [56, 57]. If applied to the fully stable regimes studied before, our results confirm the lowering of the Higgs-mass bounds to a high accuracy. In addition, the full potential solver allows us to address the fate of the RG flow in potential metastable regimes and the renormalization of the tunnel barrier.

Our present study is performed within a simple Higgs–Yukawa model with a discrete chiral symmetry, which has proved useful for addressing the qualitative properties of the IR window. Our main results are read off from the properties of the Higgs effective potential as a function of field amplitude \(\phi \) and RG scale k. A first illustration for such a fully stable flow is shown in Fig. 1, where the dimensionless potential u is depicted as a function of the dimensionless field amplitude \(\rho \sim \phi ^2\) for various values of k ranging from the UV toward the IR (blue to black) where a vacuum expectation value has developed.

Fig. 1
figure 1

Example of the RG flow (from blue to black) of the dimensionless effective potential u in dependence of the dimensionless field amplitude \(\rho \sim \phi ^2\). Here, the UV potential has been chosen as linear in \(\rho \); see Eq. (19) below for conventions. The flow results from a \(\beta \) functional of the potential u, defined for each value of the field amplitude \(\beta [u](\rho )\); see Eq. (21)

The paper is organized as follows: in Sect. 2, we introduce our toy model and briefly summarize some recent results which are of relevance for our present study. Section 3 is devoted to an extensive mean-field study, for the first time also including the metastable regime. Here, we also clarify the fate of a recently discovered pseudo-stable phase [42] and discuss the connection with the effective single-scale potential obtained from conventional perturbative constructions. In Sect. 4 we investigate the full RG flow of the entire scalar potential by means of the functional RG. We compare the local behavior around the electroweak minimum as well as the global behavior with the results obtained by polynomial expansion of the potential and the mean-field estimates. While we still span the bare potential in terms of a few polynomial operators in this work, the techniques used here will facilitate future studies of general classes of bare actions.

2 Higgs–Yukawa model with discrete symmetry

2.1 The model

All mechanisms relevant for the present work can already be studied in a simple Higgs–Yukawa model corresponding to the reduction of the standard model to the top quark \(\psi \) with the largest Yukawa coupling, and a real scalar Higgs degree of freedom \(\phi \). The classical Euclidean action of this model is given by

(1)

The model features a discrete chiral \(\mathbb {Z}_2\) symmetry, \( \psi \rightarrow e^{i\frac{\pi }{2}\gamma _5} \psi \), \(\bar{\psi } \rightarrow \bar{\psi } e^{i\frac{\pi }{2}\gamma _5} \), \( \phi \rightarrow -\phi , \) mimicking the global part of the electroweak symmetry group, and protecting the fermion against acquiring a mass term. No massless Goldstone bosons appear after spontaneous symmetry breaking as the symmetry is discrete. Hence, the particle spectrum is gapped in the broken phase as in the standard model. This toy model was intensively discussed in the context of stability of the effective potential in the literature, e.g., [40, 5861].

In order to make semi-quantitative contact with the standard model, we impose Coleman–Weinberg renormalization conditions [62] on the effective potential obtained after integrating out all fluctuations down to the IR,

$$\begin{aligned} U_{\text {eff}}'(\phi _0)=0, \quad m_{\text {H}}^2= & {} {U_\text {eff}}''(\phi _0), \quad m_{\text {t}}^2=\phi _0^2 h^2, \end{aligned}$$
(2)

where \(\phi _0\) denotes the (renormalized) field value at the minimum of the potential,Footnote 3 and all couplings are also considered to be renormalized at a suitable renormalization point \(\mu \), e.g., \(\mu =\phi _0\). In the present work, we choose for the observable parameters: \(m_{\text {t}}=173\) GeV for the top mass, and \(\phi _0\equiv v = 246\) GeV. The Higgs mass \(m_{\text {H}}\) then is treated as a function of the cutoff and a functional of the bare action, \(m_{\text {H}}=m_{\text {H}}[S_\Lambda ;\Lambda ]\).

Despite this apparent physical fixing, the simplified model, of course, deviates quantitatively from the standard model in essential aspects: for instance, whereas the center of the IR window for the Higgs mass is near \(\sim \)150 GeV for a Planck scale cutoff in the standard model [37], it is near \(\sim \)215 GeV for the present simple model at high energy scales [47] mainly due to the absence of the gauge sectors.

2.2 Perturbative effective single-scale potential

In order to make contact with the conventional perturbative treatment, we briefly sketch the standard line of argument to obtain an estimate of the effective potential. For simplicity, we consider only the one-loop level. Perturbatively, only the renormalizable operators of the bare potential are considered,

$$\begin{aligned} U_{\Lambda } = \frac{m_{\Lambda }^2}{2}\phi ^2 + \frac{\lambda _{2,\Lambda }}{8}\phi ^4 , \end{aligned}$$
(3)

featuring the bare mass parameter \(m_{\Lambda }^2\) and bare \(\phi ^4\) coupling \(\lambda _{2,\Lambda }\). The estimate for the effective potential is based on the \(\beta \) function for the renormalized running coupling \(\lambda _2\),

$$\begin{aligned} k \frac{\mathrm{d}\lambda _2}{\mathrm{d}k} \equiv \partial _t\lambda _2 = \frac{1}{16\pi ^2}(9\lambda _2^2 + 8h^2\lambda _2 - 16h^4), \end{aligned}$$
(4)

depending as well on the renormalized running Yukawa coupling h. For the present line of argument, it suffices to ignore the running of h (it will be fully included in our detailed studies later). The discussion can even be simplified further by noting that the \(\lambda _2\)-terms in Eq. (4) are small compared to the \(h^4\) term for small Higgs masses and large top masses. In this limit which corresponds to ignoring scalar fluctuations, the integration of the \(\beta \) function yields

$$\begin{aligned} \partial _t\lambda _2&= -\frac{h^4}{\pi ^2} \quad \Rightarrow \quad \lambda _2(k) = \lambda _{2,\mu } - \frac{h^4}{2 \pi ^2} \ln {\frac{k^2}{\mu ^2}}, \end{aligned}$$
(5)

with \(\mu \) denoting the renormalization point for \(\lambda _2\).

The conventional perturbative estimate of the effective potential is then inspired by the Coleman–Weinberg form of the effective potential [62]. One assumes that the effective potential is well approximated by identifying the dependence of the integrated scalar self-coupling on the RG scale k with the scalar field itself, \(\lambda _2(k=\phi )\). We emphasize that the identification \(k=\phi \) mixes momentum scale information k with the field amplitude. In general, the full effective action in field theory would provide separate information as regards the two scales which need not be the same. By this identification, we obtain a single-scale potential which in our simple approximation reads

$$\begin{aligned} U_{\text {eff}}^{\text {S}}(\phi )&= \frac{1}{2}m_{\mu }^2\phi ^2 + \frac{\lambda _2(k=\phi )}{8}\phi ^4 \nonumber \\&=\frac{1}{2} m_{\mu }^2\phi ^2 + \frac{\lambda _{2,\mu }}{8}\phi ^4 - \frac{h^4\phi ^4}{16\pi ^2}\ln {\frac{\phi ^2}{\mu ^2}}. \end{aligned}$$
(6)

Imposing the renormalization conditions (2) together with the choice \(\mu =\phi _0=v\), we can write the single-scale potential as

$$\begin{aligned} U_{\text {eff}}^{\text {S}}(\phi ) =&- \frac{1}{4}\left[ m_{\text {H}}^2 + \frac{m_{\text {t}}^4}{2\pi ^2 v^2} \right] \phi ^2 + \frac{1}{8}\left[ \frac{m_{\text {H}}^2}{v^2} + \frac{3m_{\text {t}}^4}{4\pi ^2 v^4} \right] \phi ^4 \nonumber \\&- \frac{m_{\text {t}}^4 \, \phi ^4}{16\pi ^2 v^4} \ln {\frac{\phi ^2}{v^2}}. \end{aligned}$$
(7)

Note also that the bare potential (3) remains completely unspecified in this derivation. The implicit use of only renormalizable operators together with the limit \(\Lambda \rightarrow \infty \) permitted by perturbative renormalizability seems to suggest that the details of the bare potential are irrelevant.

Clearly, this single-scale potential develops an instability for large Yukawa couplings, i.e., large \(m_{\text {t}}\). For the present choice of parameters, the instability occurs at a scale of \(\sim \)10\(^7\) GeV in our toy model; see Fig. 2. This instability is related to the running of \(\lambda _2(k)\), which turns negative at sufficiently large k; cf. Eq. (5).

Fig. 2
figure 2

Conventional effective single-scale potential \(U_{\text {eff}}^{\text {S}}\) as a function of the field amplitude \(\phi \). While the potential looks stable around the electroweak minimum, it develops an instability at large field values within our toy model. This instability seems to be driven by top fluctuations which turn the scalar self-coupling negative at large scales; cf. Eq. (5)

In the full standard model, the corresponding instability scale is of order \(\sim \)10\(^{10}\) GeV. Though current state-of-the-art calculations [34, 37, 39] determine the single-scale potential to NNLO precision, including two-loop threshold corrections, and self-consistent resummations [63, 64], the present rather cartoon-like presentation in a toy model still captures the essence of the origin of the instability occurring in the perturbative estimate of the single-scale potential.

A qualitative difference arises in the standard model from the electroweak gauge fluctuations, which render the \(\phi ^4\) coupling positive again at even higher scales such that the single-scale potential becomes bounded from below and a second minimum arises beyond the Planck scale which turns out to be the global one. Therefore, the absolute instability of the single-scale potential is a particularity of our model. Below, this will actually be useful to make one of our main points more transparent.

3 Mean-field effective potential and stability

In the following, we use mean-field methods to study the effective potential. We stick to the same simplifications as before, ignoring bosonic fluctuations and the running of the Yukawa coupling, but keep track of all scales involved, the momentum scale of fluctuations k, the field amplitude \(\phi \) and the UV cutoff scale \(\Lambda \). Parts of this discussion follows [40, 41], where also more technical details can be found. Here, we focus on the new aspects arising for un-/metastable scenarios.

3.1 Mean-field potential

With these prerequisites, the mean-field potential is directly related to the fermion determinant. More precisely, working with an explicit UV cutoff \(\Lambda \) and an IR regulator scale k, the mean-field potential reads

(8)

where \(\Omega \) denotes the spacetime volume, and irrelevant field independent constants are ignored. If we introduced N fermion flavors, the mean-field potential would become exact in the limit \(N\rightarrow \infty \). The notation \(\det _{\Lambda ,k}\) indicates that the determinant is regularized and includes momentum modes p in the range \(k^2\le p^2\le \Lambda ^2\). The result is regularization dependent. As long as we do not send \(\Lambda \rightarrow \infty \), this dependence is physical and can be viewed as a model for the details of the embedding into a more fundamental underlying UV complete theory. For a close contact with later sections, we use a piece-wise linear regulator familiar from functional RG studies [65, 66]; we emphasize that all conclusions remain the same also for a sharp momentum cutoff, propertime or zeta-function regularization, see [41]. We obtain,

$$\begin{aligned} U_k^{\text {MF}} = U_{\Lambda } - \frac{h^2 (\Lambda ^2-k^2)\phi ^2}{16\pi ^2} + \frac{h^4 \phi ^4}{16 \pi ^2} \ln {\frac{\Lambda ^2+h^2\phi ^2}{k^2+h^2\phi ^2}}, \end{aligned}$$
(9)

which makes all scale dependencies explicit. By varying the RG scale k, we can observe how the mean-field effective potential as a full function of the field amplitude \(\phi \),

$$\begin{aligned} U_{\text {eff}}^{\text {MF}}(\phi ) = U_{k=0}^{\text {MF}}(\phi )= & {} \frac{1}{2}\left( m_\Lambda ^2- \frac{h^2 \Lambda ^2}{8\pi ^2} \right) \phi ^2 + \frac{\lambda _{2,\Lambda }}{8}\phi ^4 \nonumber \\&+\, \frac{h^4 \phi ^4}{16 \pi ^2} \ln {\left( 1+\frac{\Lambda ^2}{h^2\phi ^2}\right) }, \end{aligned}$$
(10)

is built up from fermionic fluctuations renormalizing the bare potential \(U_{\Lambda }\) while running from \(k=\Lambda \) to \(k\rightarrow 0\).

Apart from the induced mass term \(\sim h^2\Lambda ^2\phi ^2\), the whole interaction part of the determinant [2nd line of Eq. (10)] is positive. The bare mass term \(m_{\Lambda }^2\) can now be fixed by the renormalization condition \(U_{\text {eff}}^{\text {MF}}{}'(\phi _0=v)=0\), fixing the Fermi scale,

$$\begin{aligned} m_{\Lambda }^2&= \frac{h^2 \Lambda ^2}{8\pi ^2} - \frac{h^4v^2}{8\pi ^2}\left[ 2\ln {\left( 1+\frac{\Lambda ^2}{h^2v^2}\right) } - \frac{\Lambda ^2}{\Lambda ^2 + h^2v^2} \right] \nonumber \\&\quad - \frac{1}{2}\lambda _{2,\Lambda }v^2. \end{aligned}$$
(11)

Inserting Eq. (11) into (10) yields a globally stable effective potential for any value of the UV cutoff \(\Lambda \) and any admissible non-negative value of the bare \(\phi ^4\) coupling \(\lambda _{2,\Lambda }\ge 0\); cf. solid black line in Fig. 3. It is important to stress that a bare potential of quartic type, Eq. (3), with negative \(\lambda _{2,\Lambda }\) would be inconsistent right from the beginning, as the functional integral over the scalar field would be ill-defined.

Fig. 3
figure 3

Comparison between the effective potential where the cutoff is kept finite (black solid line, \(\Lambda =\Lambda _{\text {cr}}=1.22 \times 10^7\) GeV) and the single-scale potential (red dashed line). Both approaches describe the same low-energy physics around the Fermi scale as they should, while at high energies a seeming instability appears in \(U_{\text {eff}}^{\text {S}}\)

For completeness of the presentation, we recall that the mass of the scalar particle now becomes a function of the cutoff and \(\lambda _{2,\Lambda }\), cf. [40, 41],

$$\begin{aligned} m_{\text {H}}^2&= U_{\text {eff}}^{\text {MF}}{}''(v) \nonumber \\&= \frac{m_{\text {t}}^4}{4\pi ^2 v^2} \left[ 2 \ln \left( 1+ \frac{\Lambda ^2}{m_{\text {t}}^2} \right) - \frac{3\Lambda ^4+2m_{\text {t}}^2\Lambda ^2}{(\Lambda ^2+m_{\text {t}}^2)^2}\right] \nonumber \\&\quad +\, v^2 \lambda _{2,\Lambda }. \end{aligned}$$
(12)

This demonstrates that a lower bound for the Higgs mass is obtained by the physical restriction that the bare potential of \(\phi ^4\)-type at a given UV cutoff \(\Lambda \) must be bounded from below, i.e., \(\lambda _{2,\Lambda }\ge 0\). Thus, the lower bound (lower edge of the IR window) is given by \(\lambda _{2,\Lambda } = 0\) for this class of bare potentials. This way of determining the lower bound has been suggested in [58, 59], and has been used in full non-perturbative lattice simulations [6770]. Generically, one observes a strong quantitative agreement with mean-field theory for this lower bound. In this fashion, strong constraints on the existence of a heavy fourth generation arise [7174].

For the purpose of the present work, we reverse the line of argument: for a given Higgs mass of, say \(m_{\text {H}}= 125\) GeV, this implies that a maximal scale of UV extension \(\Lambda \) is obtained. Choosing the minimal admissible value \(\lambda _{2,\Lambda } = 0\) a cutoff of \(\Lambda _{\text {cr}} = 1.22 \times 10^7\) GeV is obtained by writing \(\Lambda =\Lambda (m_{\text {H}}^2,\lambda _{2,\Lambda })\). For larger values of the UV cutoff, \(\Lambda >\Lambda _{\text {cr}}\) no physical (mean-field) RG trajectory can be found that connects an admissible bounded bare potential to an IR Higgs mass of 125 GeV. As long as \(\Lambda \le \Lambda _{\text {cr}}\), the bare potential as well as the effective potential do not exhibit an instability. Figure 3 shows a comparison between the mean-field potential (solid/black line) where the cutoff is kept finite and the single-scale potential (red/dashed line) where the cutoff has implicitly been sent to infinity. The single-scale potential approximation starts to break down for field amplitudes, where \(h\phi /\Lambda \gtrsim \mathcal {O}(1)\), i.e., where terms which are dropped in the implicit \(\Lambda \rightarrow \infty \) limit are actually sizable.

Fig. 4
figure 4

Approach of the mean-field potential to the single-scale potential if we blindly allow for cutoffs \(\Lambda \) larger than the critical value of \(\Lambda _{\text {cr}} = 1.22 \times 10^7\) GeV (black solid) with IR physics kept fixed; \(\Lambda = 2 \times 10^7\) GeV (blue dotted line), and \(\Lambda = 5 \times 10^7\) GeV (orange dotted line). For \(\Lambda > 10^8\) GeV, there is no visible difference between the mean-field potential and the single-scale potential (red dashed line) in this plot

It is, of course, possible to reduce the multi-scale mean-field potential to the single-scale potential. First, we blindly enforce all renormalization conditions. In particular the first condition in Eq. (2) for large cutoffs implies

$$\begin{aligned} \lambda _{2,\Lambda } = \frac{m_{\text {H}}^2}{v^2} - \frac{h^4}{2\pi ^2} \left[ \ln \frac{\Lambda ^2}{m_{\text {t}}^2} -\frac{3}{2} \right] + \mathcal {O}\left( \! \frac{1}{\Lambda ^2} \! \right) . \end{aligned}$$
(13)

Inserting Eqs. (13) and (11) into the mean-field effective potential finally leads to a potential with the requested minimum at v and Higgs mass of \(m_{\text {H}}\) by construction. The cutoff remains still a free parameter. In the naive large cutoff limit, we obtain

$$\begin{aligned} U^{\text {MF}}_{\text {``}\Lambda \rightarrow \infty \text {''}} =&- \frac{1}{4}\left[ m_{\text {H}}^2 + \frac{m_{\text {t}}^4}{2\pi ^2v^2} \right] \phi ^2 + \left[ \frac{m_{\text {H}}^2}{v^2} + \frac{3m_{\text {t}}^4}{4\pi ^2 v^4} \right] \frac{1}{8} \phi ^4 \nonumber \\&- \frac{m_{\text {t}}^4 \, \phi ^4}{16\pi ^2 v^4} \ln { \left( \frac{\Lambda ^2}{\Lambda ^2+h^2\phi ^2} \frac{\phi ^2}{v^2} \right) } + \mathcal {O}\left( \! \frac{1}{\Lambda ^2} \! \right) . \end{aligned}$$

For a cutoff larger than the critical value \(\Lambda _{\text {cr}}\), the potential develops an instability and rapidly approaches the single-scale potential, see Fig. 4. For \(\Lambda >10^8\)GeV, the difference between the mean-field effective potential with a finite cutoff and the single-scale potential with implicit limit \(\Lambda \rightarrow \infty \) becomes very small. Taking the naive limit \(\Lambda \rightarrow \infty \), the single-scale potential (7) is obtained as expected.

We emphasize that the consistency condition that the bare potential should be bounded from below for a well-defined generating functional is no longer fulfilled for all \(\Lambda >\Lambda _{\text {cr}}\). This can be directly read off from expression (13): \(\lambda _{2,\Lambda }\) has to be chosen negative for \(\Lambda >\Lambda _{\text {cr}}\), and thus already the bare potential is unstable. At this point, we conclude that the apparent instability of the single-scale potential appears due to an inconsistent UV boundary condition for the theory. As long as the consistency condition \(\lambda _{2,\Lambda }\ge 0\) is fulfilled, no instability can be found within the class of quartic bare potentials.

3.2 Generalized bare potentials

As already emphasized in [4042], these observations do not imply that in- or metastabilities are completely excluded. Whether or not an in-/metastability occurs is not a matter of the fermionic fluctuations but has to be seeded by the microscopic underlying theory. A specific example from string phenomenology is given in [75].

From the perspective of the standard model as an effective field theory, the embedding into a UV complete theory is parametrized by the bare action at the cutoff \(\Lambda \). Of course, the bare action is expected to host all operators compatible with the symmetry with couplings of order \(\mathcal {O}(1)\) in units of the cutoff \(\Lambda \).

In the following, we consider the simplest extension of the bare potential by including a higher-dimensional \(\phi ^6\) operator as an example,

$$\begin{aligned} U_{\Lambda } = \frac{m_\Lambda ^2}{2} \phi ^2 + \frac{\lambda _{2,\Lambda }}{8} \phi ^4 + \frac{\lambda _{3,\Lambda }}{48\Lambda ^2} \phi ^6. \end{aligned}$$
(14)

Within the same mean-field approximation as used before, we can straightforwardly compute the mass of the Higgs boson in our model as a function of \(\Lambda \) and the parameters \(\lambda _{2,\Lambda }\) and \(\lambda _{3,\Lambda }\), cf. Eq. (12),

$$\begin{aligned} m_{\text {H}}^2= & {} \frac{m_{\text {t}}^4}{4\pi ^2 v^2} \left[ 2 \ln \left( 1+ \frac{\Lambda ^2}{m_{\text {t}}^2} \right) - \frac{3\Lambda ^4+2m_{\text {t}}^2\Lambda ^2}{(\Lambda ^2+m_{\text {t}}^2)^2}\right] \nonumber \\&+ v^2 \lambda _{2,\Lambda } + \frac{v^4}{2\Lambda ^2}\lambda _{3,\Lambda }. \end{aligned}$$
(15)

It is obvious that the previous lower bound of Eq. (12) can be relaxed by a negative value for \(\lambda _{2,\Lambda }\), while a positive \(\lambda _{3,\Lambda }\) can stabilize the bare potential. For small negative \(\lambda _{2,\Lambda }\) and sufficiently large \(\lambda _{3,\Lambda }\) the effective potential as well as the potential at intermediate scales k are globally stable and have a unique minimum. In this regime, it is easily possible to obtain Higgs masses below the perturbative lower bound, i.e., decrease the edge of the IR window.

For even smaller \(\lambda _{2,\Lambda }\), i.e., larger absolute values of a negative \(\lambda _{2,\Lambda }\), the effective potential \(U^\text {MF}_k\) starts to develop a second minimum toward lower RG scales k and becomes metastable, while the bare potential \(U_\Lambda \) is still stable. For even smaller values of \(\lambda _{2,\Lambda }\), also the bare potential can become metastable.

For an illustration, let us assume a fixed cutoff \(\Lambda =10^7\) GeV. Within the class of quartic bare potentials (3), the lowest Higgs mass according to Eq. (12) is given by \(m_{\text {H}}=123.8\) GeV for \(\lambda _{2,\Lambda }=0\). Stabilizing the more general class of bare potentials (14) with a fixed value of \(\lambda _{3,\Lambda }=3\), we can choose negative values of \(\lambda _{2,\Lambda }\), yielding also smaller values of the Higgs mass; see Fig. 5. The resulting mean-field potentials are stable with a unique (electroweak) minimum on all scales (blue solid line) until we reach a value for the bare quartic coupling of \(\lambda _{2,\Lambda }=-0.065\). For even smaller values of \(\lambda _{2,\Lambda }\), a second minimum arises in the course of the mean-field flow, while the bare potential still has a unique minimum. This second minimum is a local minimum only for a small range of \(\lambda _{2,\Lambda }\) values, \( -0.0671< \lambda _{2,\Lambda }< -0.065\). For \(\lambda _{2,\Lambda }<-0.0671\), the second minimum becomes the global one (blue dashed line), which renders the electroweak minimum in the effective potential metastable. Within this regime of metastability, the Higgs mass can be made arbitrarily small by a suitable choice of parameters even without any metastability in the bare potential.

Fig. 5
figure 5

Higgs masses for the class of generalized bare potentials for \(\Lambda =10^7\) GeV. The bare potential is stabilized by \(\lambda _{3,\Lambda }=3\). The horizontal black solid line marks the lower Higgs-mass consistency bound within quartic bare potentials. The blue solid line indicates values for \(\lambda _{2,\Lambda }\) where the IR potential is stable while for the blue dashed line a metastability occurs

Fig. 6
figure 6

Mean-field potential (black solid) as the difference between the bare potential (blue dashed) and the absolute value of the fermion determinant (red dotted). The Fermi minimum at \(\phi =246\) GeV is hardly visible on the scale of the plot. Left panel the quartic bare potential always exceeds the contributions from the fermion loop for field values above the Fermi scale (\(\Lambda =10^7\) GeV, \(\lambda _{2,\Lambda }=0\), \(\lambda _{3,\Lambda }=0\)). Right panel a metastability seeded by the bare potential develops in the course of the RG flow (\(\Lambda =10^7\) GeV, \(\lambda _{2,\Lambda }=-0.15\), \(\lambda _{3,\Lambda }=3\))

Fig. 7
figure 7

Full mean-field potential (black solid line) and the potential approximated by a Taylor expansion (red dashed line) around the origin up to \(\phi ^8\) for different values of the RG scale k. The left panel shows the bare potential for \(\Lambda =10^9\) GeV (\(\lambda _{3,\Lambda }=3\), and \(\lambda _{2,\Lambda }\) is chosen such that \(m_{\text {H}}=125\) GeV), where the Taylor approximation fits the full potential as it should. The middle plot shows the scalar potential slightly below the UV cutoff, \(k_1=2.5\times 10^8\) GeV \(<\Lambda \), where the second minimum is built up. Toward the IR, \(k_2=5\times 10^7\) GeV, the second minimum settles while it disappears within the Taylor expansion (right plot)

It is important to emphasize that the metastability observed here in this model is a consequence of the shape of the bare potential encoded in both renormalizable and non-renormalizable operators. In the present model, this metastability remains invisible in the perturbatively estimated single-scale potential which would predict complete instability. We conclude that metastability properties of the model can only be reliably calculated if the bare potential at a UV scale is known. The single-scale potential is not sufficient as a matter of principle.

In the present model, this conclusion becomes obvious as the single-scale potential does not even exhibit a metastable region. This is different from the standard model, where the single-scale potential itself predicts metastability for light Higgs masses, as the single-scale potential is stabilized by electroweak fluctuations again at high field amplitudes. Still, the same conclusion about the reliability of the metastability estimate of the single-scale potential holds as for the simple model.

The fact that the metastability in the effective potential is seeded in the bare potential is illustrated in Fig. 6. Here, the effective mean-field potential (black solid line) is shown as the difference between the bare potential (blue dashed line) and the absolute value of the fermion determinant (red dotted line). The left panel depicts the case with stable bare as well as effective potential (initial parameters: \(\Lambda =10^7\) GeV, \(\lambda _{2,\Lambda }=0\), \(\lambda _{3,\Lambda }=0\)). By contrast, the right panel shows the case where a second minimum arises in the effective potential (initial parameters: \(\Lambda =10^7\) GeV, \(\lambda _{2,\Lambda }=-0.15\), \(\lambda _{3,\Lambda }=3\)). One clearly sees how the modified structure of the generalized bare potential with a negative \(\lambda _{2,\Lambda }\) is responsible for the second minimum at large scales besides the electroweak one (the latter at \(\phi =246\) GeV is hardly visible on the scales of the plot). We emphasize again that there is no possibility for the mean-field potential to develop a second minimum for the case of quartic bare potentials because the bare potential always exceeds the fermion determinant.

With the full mean-field potential at hand, we can also clarify the nature of the pseudo-stable phase observed in a polynomial expansion of the effective potential in [42]. In this approximation, RG flows were observed that start at \(k=\Lambda \) with a globally stable bare potential, then run through a metastable regime with two minima and finally end up in the IR \(k=0\) with one stable minimum at the Fermi scale. In the same spirit, we now expand the mean-field effective average potential (9) around the minimum at the origin and follow its flow in comparison with the flow of the full mean-field effective potential. This is depicted in Fig. 7. Indeed, the potential approximated by a polynomial expansion shows the same pseudo-stable behavior as observed in [42]. A second minimum appears but disappears again after a short RG time. The polynomial expansion thus looks stable again in the IR. This is in contrast to the full mean-field potential where the second minimum survives the RG flow toward the IR. We conclude that the pseudo-stable phase is an artifact of the finite convergence radius of the polynomial expansion. The global effective mean-field potential exhibits a metastability also in this phase.

4 Non-perturbative flow of the scalar potential

4.1 Functional renormalization group

While the mean-field approximation is highly convenient for first analytically controllable estimates, we have to go beyond for quantitative accuracy. The functional RG is an ideal tool for this task, as it resums large classes of higher-order diagrams, automatically accounts for threshold corrections and provides information as regards the global RG flow of the effective potential, i.e., all relevant scales can be studied independently. We use the Wetterich equation [76] in order to compute the RG flow of the effective action \(\Gamma _k\),

$$\begin{aligned} \partial _t\Gamma _k = \frac{1}{2} \text {STr}\left[ \partial _tR_k\Big (\Gamma _k^{(2)}+R_k\Big )^{-1} \right] , \quad t = \ln \frac{k}{\Lambda }. \end{aligned}$$
(16)

The solution to this equation interpolates between the bare action \(S_\Lambda =\Gamma _{k=\Lambda }\) and the full effective action \(\Gamma =\Gamma _{k\rightarrow 0}\), which accounts for all quantum fluctuations. The regulator function \(R_k\) implements the shell-by-shell integration at the momentum scale k. \(\Gamma _k^{(2)}\) denotes the Hessian of the effective average action with respect to the fields in the theory, \((\hat{\phi },\psi ,\bar{\psi })\). For detailed reviews, see [7781].

We solve the Wetterich equation within a systematic derivative expansion of the action at next-to-leading order,

(17)

The Wetterich equation (16) provides flow equations for the potential, the Yukawa coupling, and the wave function renormalizations for the fields \(Z_{\phi ,k}\) and \(Z_{\psi ,k}\). The latter can be summarized by the anomalous dimensions of the fields,

$$\begin{aligned} \eta _{\phi } = -\frac{\partial _tZ_{\phi ,k}}{Z_{\phi ,k}}, \quad \eta _{\psi } = -\frac{\partial _tZ_{\psi ,k}}{Z_{\psi ,k}}. \end{aligned}$$
(18)

In Eq. (17), we have introduced the unrenormalized scalar field \(\hat{\phi }\), which is related to the renormalized field by \(\phi =Z_{\phi ,k}^{1/2} \hat{\phi }\). At mean-field level, the distinction is irrelevant as \(Z_{\phi ,k}|_{\text {MF}}=1\). In terms of dimensionless renormalized quantities,

$$\begin{aligned} \rho= & {} Z_{\phi ,k} k^{2-d}\, \frac{1}{2}\hat{\phi }^2 = \frac{k^{2-d}}{2} \phi ^2, \quad u = k^{-d}U_k, \end{aligned}$$
(19)
$$\begin{aligned} h^2= & {} Z_{\phi ,k}^{-1}Z_{\psi ,k}^{-2} k^{d-4} \bar{h}_k^2, \end{aligned}$$
(20)

the non-perturbative flow equations in agreement with [40, 41, 82] read in d spacetime dimensions

$$\begin{aligned} \partial _t\, u&\equiv \beta [u]= -d\, u + (d-2+\eta _{\phi }) \rho u'\nonumber \\&\quad + 2v_{d} \, \Big [ l_{0}^{d} \left( u' + 2\rho u''; \eta _{\phi } \right) - d_{\gamma } \, l_{0}^{(F)d}\left( 2\rho h^{2}; \eta _{\psi } \right) \Big ], \end{aligned}$$
(21)
$$\begin{aligned} \partial _t\, h^2&= \left[ \eta _{\phi } + 2\eta _{\psi } + d - 4 \right] h^{2}\nonumber \\&\quad + 8v_{d} h^{4} \Big [ l_{1,1}^{(FB)d}\left( 2 h^2 \kappa ,u' + 2\kappa u'';\eta _{\psi },\eta _{\phi }\right) \nonumber \\&\quad - \left( 6 \kappa u'' + 4 \kappa ^{2} u''' \right) \, l_{1,2}^{(FB)d}\left( 2 h^2 \kappa ,u' + 2\kappa u'';\eta _{\psi },\eta _{\phi } \right) \nonumber \\&\quad - 4 h^{2} \kappa \, l_{2,1}^{(FB)d}\left( 2 h^2 \kappa ,u' + 2\kappa u'';\eta _{\psi },\eta _{\phi }\right) \Big ]_{\rho =\kappa }, \end{aligned}$$
(22)
$$\begin{aligned} \eta _{\phi }&= \frac{8v_{d}}{d} \bigg [ \kappa \left[ 3 \, u'' + 2 \kappa u''' \right] ^{2} m_{4,0}^{d}\left( u' + 2 \kappa u'';\eta _{\phi }\right) \nonumber \\&\quad + d_{\gamma } \, h^{2} \Big [m_{4}^{(F)d}\left( 2h^2 \kappa ;\eta _{\psi }\right) -2 h^2 \kappa \, m_{2}^{(F)d}\left( 2 h^2 \kappa ;\eta _{\psi }\right) \Big ] \bigg ]_{\rho =\kappa }, \end{aligned}$$
(23)
$$\begin{aligned} \eta _{\psi }&= \frac{8 v_{d}}{d} h^2 m_{1,2}^{(FB)d}\left( 2 h^2 \kappa ,u' + 2\kappa u'';\eta _{\psi },\eta _{\phi }\right) \Big |_{\rho =\kappa }, \end{aligned}$$
(24)

where primes denote derivatives with respect to \(\rho \). Here, we extract the flow of the Yukawa coupling and the anomalous dimensions at the (running) minimum of the potential \(\kappa = \rho _{\text {min}} \sim \phi _0^2\); see [83] for an extended flow in the Yukawa sector. The threshold functions l and m encode the decoupling of massive modes. Evaluated for a piece-wise linear regulator function [65, 66], these are listed, for instance, in [40]. Of course, the perturbative \(\beta \) functions can be recovered from these non-perturbative flow equations by neglecting resummation and threshold effects. The flow equation also contains mean-field theory as a simple limit: ignoring the anomalous dimensions as well as the flow of the Yukawa coupling, and dropping the bosonic fluctuations \(\sim l_0^d\) in Eq. (21), the integration of the potential \(\partial _tu\) precisely yields the mean-field potential (9). Hence, all known limits are unified in the functional RG framework, which we can now use to go beyond the perturbative/mean-field studies.

Previous work on Higgs boson mass bounds has solved the potential flow by means of a polynomial expansion [40, 41, 48, 49] about the flowing minimum. More concretely,

$$\begin{aligned} \begin{array}{ll} u &{}= \displaystyle \sum _{n=1}^{N_{\text {p}}} \frac{\lambda _n}{n!} \rho ^n \,\, \text {(SYM)} \quad \text {or} \\ u &{}= \displaystyle \sum _{n=2}^{N_{\text {p}}} \frac{\lambda _n}{n!} (\rho -\kappa )^n \,\, \text {(SSB)} \end{array} \end{aligned}$$
(25)

approximate the potential by a polynomial of degree \(N_{\text {p}}\) in the symmetric (SYM) or symmetry-broken (SSB) regime. For the present problem, the quality of this expansion has been confirmed for potentials [40] with a single minimum at any given scale. As the example of the seeming pseudo-stable phase above has shown, a proper description of metastability doubtlessly requires a PDE solver for the RG flow of the full potential (21) as a function of both k and \(\rho \).

Solvers for such types of Yukawa systems have been successfully developed and applied in various functional RG studies from low-energy QCD models [84, 85], critical phenomena [8688] to ultracold-gas systems [89]. The particular difficulty in the present case arises from the necessity to run the RG over many orders of magnitude in the presence of a relevant operator \(\phi ^2\) of canonical dimension 2. This requires substantial precision.

Another challenge is the approach to convexity which is expected to hold for the full effective potential [90, 91], but is spoiled by both the perturbative single-scale potential and the mean-field approximation. Whereas the approach to convexity is less interesting for the case of a single minimum, it may become essential for metastable scenarios as the tunneling lifetime depends on the shape of the manifestly non-convex tunnel barrier.

4.2 Pseudo-spectral flows

In general, the derivative expansion of the functional RG results in a system of coupled partial differential equations (PDEs). In the present case, we have to deal with one such PDE (21) coupled to an ODE (22) and two algebraic equations (23) and (24).

In order to obtain global information with high accuracy, we advocate pseudo-spectral methods [92] which span the solutions in terms of global basis functions of high polynomial (or rational) degree. Under mild analyticity assumptions, the convergence to the exact result is exponential [92]. In particular in comparison with finite difference methods, pseudo-spectral methods are memory minimizing, as a certain accuracy requires only a comparatively small number of grid points.

In the following, we briefly sketch our methods; for more details, see [57]. It is worth mentioning that pseudo-spectral methods have already successfully been applied to various problems in physics in general [93]; for first applications in the context of the functional RG, see [94, 95]. The present code development is based on a highly accurate pseudo-spectral solver for computing global fixed-point solutions within the functional RG [56].

In principle, pseudo-spectral methods include any kind of basis function system. In the present work, we use Chebyshev polynomials of the first kind. Boundary effects can be controlled since the polynomials are defined on a finite interval. The behavior of higher-order coefficients provides an error estimate of the absolute remainder of the solution. As a main advantage also in practice, function values, derivatives, and integrals of the objective function are easily accessible from the Chebyshev coefficients by means of recursive algorithms yielding a high precision calculation.

We apply pseudo-spectral methods in both the field and the RG time direction. For this purpose, we subdivide the time interval into patches, apply a Newton–Raphson iteration to each patch, and solve for the coefficients. For minimizing the number of coefficients and increasing the resolution in physically interesting regions, we use multiple domains in field direction. All these patches are connected by additionally demanding matching conditions for the objective function and a sufficient number of its derivatives. The Newton–Raphson iteration provides an error estimate for the solution of the equations.

As a result, this method is stable over many orders of magnitude in k. This enables us to choose high UV cutoffs. This choice is solely restricted by the number of digits needed for tuning the IR quantities. As the flow of the present problem includes a scalar mass parameter with canonical scaling of dimension 2, we need to tune approximately twice as many digits as the number of orders of magnitude between the UV scale and the Fermi scale. All full potential computations have been done with long double. Thus, we restricted ourselves to a maximal UV cutoff of \(10^{10}\)GeV for the full potential calculation. In principle, higher values for \(\Lambda \) are straightforwardly accessible by using a higher accuracy for the floating-point arithmetics.

Fig. 8
figure 8

Higgs boson mass as a function of the UV cutoff for various bare potentials. The filled circles are obtained by solving the full PDE system. These match perfectly with the Higgs masses computed within a polynomial expansion (25) of the scalar potential for the class of \(\phi ^4\)-type bare potentials (black, conventional “lower bound” \(\lambda _{2,\Lambda }=0\)) as well as generalized bare potentials for the case where the effective potential is stable for \(\Lambda \gtrsim 10^{4}\) GeV (red, \(\lambda _{2,\Lambda }=-0.1\), \(\lambda _{3,\Lambda }=3\)) or develops a metastability (orange, \(\lambda _{2,\Lambda }=-0.15\), \(\lambda _{3,\Lambda }=3\))

4.3 Higgs-mass bounds

As a first benchmark, we perform a comparison to local polynomial solutions of the flow equation. For this purpose, we compute Higgs masses for different initial values for the flow equations over a large range of cutoff values. In Fig. 8, we depict the resulting IR Higgs mass as a function of the UV cutoff \(\Lambda \): for the restricted class of \(\phi ^4\) bare potentials, the black data shows the resulting lowest possible Higgs mass, i.e., the conventional lower bound for \(\lambda _{2,\Lambda }=0\). Examples within the class of generalized bare potentials that lead to a relaxation of the lower bound are shown in red (\(\lambda _{2,\Lambda }=-0.1\), \(\lambda _{3,\Lambda }=3\)) and orange (\(\lambda _{2,\Lambda }=-0.15\), \(\lambda _{3,\Lambda }=3\)). The solid lines mark the Higgs masses computed within the polynomial truncation, and filled circles indicate the Higgs masses resulting from the pseudo-spectral full potential computation of this work. The black and red line agree with the results of [40], and illustrate the relaxation of the conventional lower bound by higher-dimensional operators. The orange data corresponds to a potential that develops a metastability (i.e., seems pseudo-stable in the polynomial expansion). For all cases, the pseudo-spectral data lies perfectly on top of the polynomial results. The full numerical PDE solution thus provides a strong confirmation that the polynomial expansion is suitable for extracting local information such as the Higgs mass (\(\sim \) curvature of the potential at \(\phi =v\)).

Fig. 9
figure 9

Comparison between the effective potential of the full calculation (black, \(k_\text {IR}=113.9\) GeV), the mean-field result (green dashed) and polynomial truncated potentials (dotted) up to fourth order (orange), sixth order (red) and eighth order (purple) for the stable case. The cutoff is chosen to be \(\Lambda = 10^9\) GeV and bare couplings are tuned such that the Higgs mass is 159.4 GeV in all cases for reasons of comparison. The vertical gray-dashed line indicates the position of the Fermi minimum at \(\phi =246\) GeV

4.4 Global flows

Let us start with a closer look at the global behavior of the flow for the class of the \(\phi ^4\) bare potentials. Obviously, the polynomial truncation lacks in describing the asymptotic behavior of the potential which can be seen in Fig. 9. This is not surprising since the flow equations suggest the asymptotic behavior to be that of the UV potential \(\sim \phi ^4\) because fluctuations for large field amplitudes are suppressed. By contrast, the asymptotic behavior of the polynomial expansion by construction is fixed to the highest power of the field which is taken into account in the truncation, \(\sim \phi ^{2N_{\text {p}}}\). These higher-order couplings are generated during the RG flow, even if the bare potential is of \(\phi ^4\) type. Therefore, considering only terms up to \(\phi ^4\), (accidentally) displays the asymptotic behavior best.

Naively, the polynomial truncation up to sixth order in \(\phi \) seems to suggest an instability; however, the inflection point is beyond the radius of convergence of the polynomial expansion around the Fermi scale. This radius of convergence is approximately of the order of the curvature around the electroweak minimum [40]. For large field values the polynomial expansion behaves like an asymptotic series with alternating signs between the coefficients. Incidentally, an alternating series is also obtained from the polynomial expansion of the mean-field effective potential. As long as the \(\phi ^4\) class of bare potentials is considered, no hint for an in-/metastability can be found within the radius of convergence of the polynomial expansion. This is confirmed by the fully stable potential obtained from the global pseudo-spectral flow.

We observe that the mean-field potential agrees quite well with the results for the full potential, for small as well as for larger field values, see green dashed curve in Fig. 9. The fluctuations of the bosons appear to play a minor role in this parameter regime near the lower edge of the IR window, where the fermionic contributions dominate. Moreover, neglecting the anomalous dimensions and the flow of the Yukawa coupling does not have a significant effect. Thus, the simple mean-field effective potential is an effective tool to get first insights into the global behavior of the scalar potential, justifying the seemingly severe approximations of Sect. 3.

The qualitative picture remains the same for full flows also in the class of generalized bare potentials. For those cases where no second minimum emerges during the polynomial truncated flow, we observe that also the full flow does not develop an outer minimum for field values larger than the Fermi scale. Our new results hence confirm the existence of a class of stable bare potentials giving rise to Higgs masses below the conventional lower bound as, for instance, depicted by the red curve in Fig. 8 for \(\Lambda \gtrsim 10^{4}\) GeV. Therefore, the mechanism of lowering the lower bound for completely stable potentials remains active beyond the polynomial expansion and the mean-field analysis. Thus higher-order operators can diminish the lower Higgs-mass bound of the standard model.

Let us take a closer look at the inner workings of the equations: for large field amplitudes in the asymptotic regime of the potential, all fluctuations are suppressed as the threshold functions approaching zero. In this regime, the flow is dominated by dimensional scaling, i.e., the first two terms in Eq. (21). This holds for the \(\phi ^4\) as well as for the generalized class of bare potentials irrespectively of the stability properties.

Deviations from the mean-field limit require relevant bosonic fluctuations. In the small coupling (i.e., small Higgs-mass) regime, this can indeed occur in the full flow due to threshold effects of the following type: if a second minimum emerges seeded by a suitable bare potential, the curvature near the top of the barrier between the minima is negative, such that the bosonic threshold function \(l_0^d(u'+2\rho u''; \eta _\phi )\) is enhanced. This type of bosonic enhancement is only visible in a full potential flow. For a meaningful comparison between mean-field and full flow, we tune the bare potential such that the IR physics including the Higgs mass is kept fixed. For illustrative purposes, we choose a UV cutoff of only 5 TeV, a Higgs mass of 25GeV and \(\lambda _{3,\Lambda } = 1\). The resulting potentials end up in the metastable regime as depicted in Fig. 10. The mean-field potential (dashed line) differs from the full solutions (solid lines) integrated down to an IR value of \(k_{\text {IR}}\sim 100\) GeV. The full solutions correspond to defining the anomalous dimensions at the local Fermi minimum (black) or at the second global minimum (red). These curves differ as fields and couplings are renormalized differently within the two schemes, with the black curve corresponding to a renormalization choice better resolving the Fermi minimum and the red curve better suited for the second minimum. In other words, the axes for the different solid lines have a different meaning due to the different field rescaling during the flow. If the anomalous dimensions were artificially set to zero, mean-field and full potential results would still match rather well.

Fig. 10
figure 10

Comparison between the effective potential for a metastable case obtained by the full pseudo-spectral calculation where the anomalous dimensions are computed at the electroweak vev (black/lower solid line, \(k_\text {IR}=105.5\) GeV), at the outer minimum (red/upper solid line, \(k_\text {IR}=93.3\) GeV) and the mean-field result (green dashed line). Setting the anomalous dimensions by hand to zero would yield a rather good agreement with the mean-field potential. The Higgs mass is tuned to \(m_{\text {H}}=25\) GeV where the UV cutoff is 5 TeV and \(\lambda _{3,\Lambda } = 1\)

4.5 Convexity of the effective potential

By definition of the effective action as a Legendre transform of the Schwinger functional, we expect the full effective action and particularly the effective potential to be convex functions of the field. This convexity property cannot be seen in the perturbative construction of the single-scale potential nor in the mean-field approximation. The Wetterich equation does have this convexity property in the limit \(k\rightarrow 0\) for the bosonic potential, see e.g. [77, 91]. However, at finite k, the regulator term \(\sim R_k\) sources a non-convex contribution which vanishes in the limit \(k\rightarrow 0\).

For potentials with a single minimum, it is well known that convexity of the running potential sets in rather late in RG time, i.e., convexity is driven by the very deep IR modes which are often no longer relevant for the IR observables. For instance in the examples above, we have stopped the flow at scales \(k_{\text {IR}}\sim 10\ldots 100\) GeV, where the IR Fermi scale observables have already settled to their physical values. Still, for these values of k, the approach to convexity has not fully set in yet. Whereas this demonstrates that convexity is not important for the static observables, it is an interesting question as to whether the approach to convexity can be important for estimates of the tunneling rate between two different minima. The relevance of this question becomes obvious from the fact that any tunneling barrier in a convex potential is (naively) exactly zero by construction.

In order to understand the onset of convexity in the present model, let us start with the simpler case with only one minimum at the Fermi scale. Here, convexity only affects the inner region of the potential with \(\phi <v\). Technically, convexity of the effective potential is generated by singularities in the bosonic propagators entering the threshold functions. In the present case, the bosonic threshold function \(l_0^d\) corresponding to the regularized propagator is proportional to

$$\begin{aligned} l_0^d \sim \frac{1}{1+u'+2\rho u''}, \end{aligned}$$
(26)

exhibiting a singularity at \(u'+2\rho u'' \rightarrow -1\), or \(u'\rightarrow -1\) for small \(\rho \) or small \(u''\). The flow avoids this singularity by renormalizing the negative curvature of the potential in the inner region \(0\ge U_k''(\phi )\sim k^2( u'+ 2\rho u'') \gtrsim - k^2\rightarrow 0\) with \(k\rightarrow 0\). This establishes convexity for \(k\rightarrow 0\).

In comparison to purely bosonic models, fermionic fluctuations delay convexity since they enter the flow equation with an opposite sign, cf. the last term in Eq. (21). Thus, bosonic fluctuations have to exceed the fermionic fluctuations first. As convexity also introduces a non-analyticity at \(\phi =v\), its onset becomes numerically first visible in higher derivatives. Therefore, we consider the first derivative of the potential \(u'\) in the following. The balancing between bosons and fermions also implies that the onset of convexity becomes more pronounced if the boson coupling \(\lambda _2\) is enhanced relative to the fermion coupling h. In terms of physical parameters, this implies that convexity should become more prominent for larger Higgs-to-top mass ratios. In Fig. 11, we plot \(u'\) for three different ratios \(m_{\text {H}}/m_{\text {t}}\). The flow has been stopped at a scale \(k_{\text {IR}}\) such that the potentials have the same distance from the singularity of the bosonic propagator \(1-|u'(0)|=0.01\). The faster approach to convexity then is directly visible in terms of the position of non-analyticity \(\rho _{\text {kink}}\) which we observe to move toward larger field amplitudes if \(m_{\text {H}}/m_{\text {t}}\sim \lambda _2/h\) increases. By contrast, if \(m_{\text {H}}/m_{\text {t}}\) is small, the characteristic flat region of \(u'\) is hardly visible at this particular scale \(k_{\text {IR}}\) and would increase only toward even smaller scales.

It should be emphasized that convexity is a notoriously difficult problem for pseudo-spectral methods, since it introduces a non-analyticity which violates the assumptions on the function space for which exponential convergence can be proved. Hence, the numerics will unavoidably face a singularity problem in the very deep IR. For further adapted methods, see [96, 97].

Fig. 11
figure 11

The approach to convexity of the potential is faster if the relation between bosonic and fermionic coupling \(m_{\text {H}}/m_{\text {t}}\) increases (from green/top to blue/bottom \(m_{\text {H}}/m_{\text {t}}= 0.23, 0.66, 1.12\)). Here, the first derivative of the potential as a function of the dimensionful field invariant is depicted. All potentials exhibit a minimum (i.e. \(u'=0\)) at \(\phi =246\) GeV (\(\phi ^2/2{=}30258\,\text {GeV}^2\)). The approach to convexity becomes manifest by a characteristic flattening of the inner region and \(u'\) approaching \(u'\rightarrow -1\), see text. We have chosen \(\Lambda = 10^3\) GeV and \(m_{\text {H}}=39.7\) GeV for the green/top curve stopping the flow at \(k_\text {IR}=33.4\) GeV, and \(\Lambda = 10^6\) GeV and \(m_{\text {H}}=113.6\) GeV for the orange/middle curve stopping at \(k_\text {IR}=80\) GeV. The blue curve is added for illustration; here \(\Lambda = 6.5\times 10^4\) GeV, \(\text {vev}=246\) GeV, \(m_{\text {t}}=426.3\) GeV, \(m_{\text {H}}=476.6\) GeV, and \(k_\text {IR}=264.5\) GeV, such that the curve is not tuned to the physical top mass in terms of the renormalization conditions (2)

Fig. 12
figure 12

Onset of convexity at a scale \(k\sim 100\) GeV for the case of two minima. Whereas the small-field region is still dominated by fermionic fluctuations, the bosons control the flow for larger field amplitudes, especially in the region of the minimum of the field-dependent mass term \(u'+2\rho u''\). Left panel mass term and the first derivative of the potential \(u'\) over a wide field range. The vertical dashed lines mark the location of both minima and the maximum of the tunnel barrier of the potential in between. The minimum of the mass term approaches the singularity \(u'+2\rho u''\rightarrow -1\), triggering the approach to convexity. This is shown in the right panel in detail for decreasing scale \(k\in \{112.5,110,109.5\}\) GeV from top (green) to bottom (blue). For this example flow, we have used \(m_{\text {H}}=24.1\) GeV, \(\lambda _{3,\Lambda } = 1\), and \(\Lambda =5\) TeV

Let us now turn to the more interesting case of two minima which is numerically more challenging since the field-dependent “mass term” becomes negative, \(u'+2\rho u''<0\), not only for small fields but also for a second region at larger fields. As an illustrative example, we choose a similar flow as above. We plot this mass term in the region of both minima of the tunnel barrier; cf. Fig. 12 (left panel). For comparison we show \(u'(\rho )\) as well. The dashed vertical lines indicate the position of both minima and the maximum in between. Convexity becomes first visible for larger fields at the minimum of the mass term which tends to \(u'+2\rho u'' \rightarrow -1\) after k has dropped below the scale of the top quark; cf. Fig. 12 (right panel). For the current example, this minimum of the mass term is located in between the maximum and the outer minimum of the potential, but this relative position may vary depending on the scale and the precise choice of parameters. As the maximum of the potential is situated within the region where \(u'+2\rho u''<0\) for larger fields, the flat region eventually extends beyond the maximum, significantly affecting the tunnel barrier for small scales \(k\lesssim 100\,\text {GeV}\). For small fields, the fermions still control the flow at \(k\sim 100\) GeV. However, for decreasing scale k the bosonic fluctuations win out over the fermionic ones and convexity sets in as well, similar to Fig. 11. We emphasize that the approach to convexity appears to set in at different scales for large field amplitudes than for small ones.

In the present example, convexity affects the tunnel barrier at scales k which are more than an order of magnitude smaller than the field amplitude of the barrier and the outer minimum. A calculation of the tunnel rate which is dominated by the latter scales hence is expected to be only weakly influenced by the approach to convexity. As a general rule, we conclude that the standard recipes for calculations of the tunnel rate [98, 99] remain unaffected as long as the fermion fluctuations dominate the renormalization of the potential. Whether or not this is the case at the relevant scales of interest will in general depend on the details of the scale-dependent potential and thus also on the details of the bare potential. As soon as the boson fluctuations become important, the approach to convexity has also be accounted for in estimates of the tunneling rate.

In the functional RG context, a proposal for this has been worked out for scalar models in [100]. A formalism that can also systematically deal with further radiative effects in the resulting inhomogeneous instanton background on top of a radiatively generated potential has recently been developed with the help of a self-consistent functional scheme based on the 2PI effective action [101103].

We would like to emphasize the necessity of a simultaneous consistent treatment of the renormalization flow of the potential together with the fluctuation contributions in a tunnel-rate calculation—even if the bare potential was known exactly. Of course, unknown higher-dimensional operators then further add to the indeterminacy of the vacuum decay rate [42, 104107]. For instance, the influence of gravity-induced higher-dimensional operators has been studied in [108111].

4.6 Quantum phase diagram of the Higgs–Yukawa model

Can the outer global minimum be used to define the electroweak vacuum? If the occurrence of metastability is rather generic in presence of higher-dimensional operators, could it be possible to fix physical parameters with respect to the global minimum as the Fermi scale? In order to address these questions, we now reconsider the model from a more general viewpoint.

So far, we have fixed the model with the help of the renormalization conditions (2) applied to the first or innermost minimum. Instead, let us now start from a fixed UV cutoff \(\Lambda \) with some bare potential bounded from below and read off the IR phases from the effective potential at some IR scale k where all modes have decoupled (apart from the approach to convexity). We are most interested in this quantum phase diagram as a function of the (super-)renormalizable operators \(\sim m_\Lambda ^2 \phi ^2\) and \(\sim \lambda _{2,\Lambda }\phi ^4\), as the electroweak precision data tells us that the standard model is sufficiently close to the Gaußian fixed point, where perturbation theory based on these operators works very well. In other words, higher-dimensional operators do not take a momentarily measurable influence on collider data.

In the language of critical phenomena, the standard model appears to be close to a second-order quantum phase transition that effectively allows one to push the UV cutoff to large values (compared to collider scales). The natural candidate in the standard model is the electroweak (quantum) phase transition represented by the order–disorder phase transition of discrete chiral symmetry in our simple model. It is, in fact, straightforward to verify by means of perturbation theory, mean-field theory or the functional RG that this phase transition is of second order for \(\phi ^4\) type bare potentials in the stable regime. The “control parameter” for the quantum phase transition is the bare mass term \(m_\Lambda ^2\).

Fig. 13
figure 13

Quantum phase portraits of the IR effective potentials with possible metastabilities seeded from the bare action. As an example, the large amplitude \(\rho \sim \phi ^2\) region is stabilized by a \(\phi ^6\) operator in the UV. The phase portraits are sketched for various initial values for a negative bare \(\lambda _{2,\Lambda }\) and suitable mass parameters \(m_\Lambda ^2\) near the critical regions. In the dotted region the effective potential is metastable, while it is stable in the white and gray-shaded region. In the gray-shaded region only the outer minimum exists

In the following, we perform this investigation for the class of generalized bare potentials. For this purpose, we fix \(\lambda _{3,\Lambda }=1\) as a representative of a higher-dimensional operator that induces absolute stability. We expect the following results to hold also for other polynomial operators that ensure absolute stability for large field amplitudes. For technical simplicity, we keep the Yukawa coupling \(h^2\sim \mathcal {O}(1)\) fixed and also neglect the anomalous dimensions, as both do not induce qualitative differences. Still, we keep the full bosonic-fluctuation contribution to the flow of the potential.

Choosing \(\lambda _{2,\Lambda }\) negative but with a small absolute value, the potential will still show only one minimum and the phase transition as a function of \(m_{\Lambda }^2\) still is of second order as for the \(\phi ^4\)-class; cf. left-hand side of Fig. 13. Increasing the absolute value of a negative \(\lambda _{2,\Lambda }\) a bit, and starting with a large value of \(m_\Lambda ^2\), the sufficiently negative \(\lambda _{2,\Lambda }\) may seed a local higher minimum at large field amplitudes. Nevertheless, the system is in the symmetric phase with the global minimum at \(\phi =0\) (upper left part of Fig. 13). (On the left-hand side of this figure, we do not further distinguish between the existence or non-existence of a further local outer minima; potentials with a local outer minimum shown here only represent possible examples.)

Decreasing \(m_\Lambda ^2\), we indeed observe a second-order phase transition to a broken phase driven by fermion fluctuations where the order parameter \(\phi _0=v\) is switched on continuously; cf. the white region in Fig. 13. A local higher minimum at larger field amplitudes may arise by decreasing \(m_\Lambda ^2\) or persists if it already existed. It is this second-order phase transition which can serve to define a “continuum limit” essentially establishing cutoff independence.

Decreasing \(m_\Lambda ^2\) further first leads to a lowering of the outer minimum such that the inner minimum becomes metastable (dotted region). The phase transition between the two cases is of first order (dashed lines). For even smaller mass parameters \(m_\Lambda ^2\), the inner minimum vanishes discontinuously while the outer remains (gray-shaded region in Fig. 13). We also classify this discontinuous change of the system as a first-order transition, even though it would not correspond to a thermal phase transition: on both sides of the lower dashed line the system is dominated by the global minimum in the thermodynamic limit.

This analysis demonstrates that only the transition between the symmetric and the symmetry-broken phase, which is driven by fermion fluctuations, is of second order. Therefore, only this transition can be used to separate the cutoff scale from the IR physics in this model. This forces us to associate the Fermi scale with the innermost minimum driven by fermion fluctuations. The outer one seeded by the bare action cannot be used for a definition of the Fermi scale as it is highly unlikely to match with the perturbative description of electroweak precision data. In our flows, this mismatch becomes visible in the practical difficulty to push the cutoff beyond \(\sim 1\) TeV while satisfying all renormalization conditions with v corresponding to the outer minimum.

For negative \(\lambda _{2,\Lambda }\) with an even larger absolute value (right-hand side of Fig. 13), the phase portraits are similar in the sense that only the transition driven by the fermionic fluctuations in the inner region of the potential is a second-order transition. The difference is that this transition occurs only after the outer minimum seeded by the bare action has become the global one. As a consequence, both the symmetric phase for larger \(m_\Lambda ^2\) and the fermion driven broken phase (inner minimum) are metastable (dotted region) on both sides of the transition. In this regime of bare potentials, the separation of IR physics from the UV scale hence goes along with a metastability.

We can also compare the phase portraits for fixed \(m_\Lambda ^2\) and use \(\lambda _{2,\Lambda }\) as control parameter. For instance, the transition marked by the gray thick arrow is a first-order broken-to-broken transition. This is likely to correspond to an equivalent transition first observed in lattice simulations of a similar chiral model [44].

We emphasize that the phase portraits determined here correspond to quantum phase transitions with control parameters corresponding to parameters of the bare action. This is a priori unrelated to the nature of finite temperature phase transitions in the same model, even though a relation might be established dynamically because of a thermal decoupling of the fermions. For recent lattice studies, see [46, 112].

5 Conclusions

We have investigated the renormalization flow of the effective scalar potential keeping the full dependence on all relevant scales: the field amplitude \(\phi \), an RG scale k and the UV cutoff scale \(\Lambda \). This is necessary to overcome the limitations of conventional approximation schemes, relying on identifications such as \(k=\phi \), or implicit perturbative limits \(\Lambda \rightarrow \infty \). The advantage of keeping the full scale dependence becomes obvious for the analysis of metastabilities.

Using a simple Higgs–Top–Yukawa model as an example, our analysis demonstrates that metastabilities are not primarily induced by fermion fluctuations, but have to be seeded by suitable structures in the bare potential. In particular, metastabilities cannot occur if a well-defined bare action is restricted to contain only renormalizable operators. Upon the inclusion of suitable higher-dimensional operators, metastabilities generically occur for small Higgs masses and large cutoffs—at least within the class of simple polynomial bare potentials studied in this work.

Whereas the latter conclusions are in part reminiscent of results obtained within perturbative estimates based on a single-scale effective potential (\(k=\phi \)), it is worthwhile to note some decisive differences: our approach allows for addressing arbitrary bare potentials, defining the model purely in terms of its symmetry, field content and a minimal set of IR parameters. No assumption as regards the manifest absence of an infinity of irrelevant operators is made. While the occurrence of higher-dimensional operators is conventionally interpreted as particle physics beyond the standard model (e.g., induced by integrating out heavy particle multiplets), our approach allows us to also associate such operators to any origin that can be parametrized in terms of any effective action, e.g., a coarse-grained continuum action in a spacetime arising from discrete building blocks.

In this work, we have studied global RG flows of the potential using pseudo-spectral methods. This facilitates to study the full RG evolution of competing minima, to analyze the quantum phase diagram of the model, and to quantify the approach to convexity. The latter is not accessible in perturbation theory or mean-field theory. The global flows also serve to confirm earlier results from local flows evaluated at the running Fermi minimum to a high accuracy, such as, for instance, the relaxation of the perturbative lower bound on the Higgs mass. On the other hand, the global flow also reveals the limitations of the local flow in metastable regimes as competing minima turn out to be beyond the radius of convergence of local flows. Further, the global flows also demonstrate the usefulness of the mean-field approximation in the small-Higgs-mass regime.

Finally, we emphasize that a full determination of consistency bounds for the IR observables of the standard model as a function of the cutoff \(\Lambda \) as the scale of maximum UV extension has not yet been completed. For this purpose, the mapping of a wide range of bare actions to the IR observables would have to be computed with the RG, technically corresponding to an extremization problem in an infinite-dimensional space. The capability of handling global flows and extending the current studies to nonpolynomial interactions will be a necessary prerequisite for this.