# Exact renormalization group equation for lattice Ginzburg-Landau models
adapted to the solution in the local potential approximation

V. I. Tokar1,2 1Université de Strasbourg, CNRS, IPCMS, UMR 7504, F-67000
Strasbourg, France  2G. V. Kurdyumov Institute for Metal Physics of the N.A.S.
of Ukraine, 36 Acad. Vernadsky Boulevard, UA-03142 Kyiv, Ukraine

###### Abstract

The Wilson Green’s function approach and, alternatively, Feynman’s diffusion
equation and the Hori representation have been used to derive an exact
functional RG equation (EFRGE) that in the course of the RG flow interpolates
between the interaction part of the lattice Ginzburg-Landau Hamiltonian and
the logarithm of the generating functional of the S-matrix. Because the
S-matrix vertices are the amputated correlation functions of the fluctuating
field, it has been suggested that in the critical region the amputation of the
long-range tails makes the S-matrix functional more localized and thus more
amenable to the local potential approximation (LPA) than the renormalized free
energy functional used in Wilson’s EFRGE.

By means of a functional Legendre transform the S-matrix EFRGE has been
converted into an EFRGE for the effective action (EA). It has been found that
the field-dependent part of EA predicted by the equation is the same as
calculated within the known EA EFRGE approaches but in addition it is
accurately accounts for the field-independent terms. These are indispensable
in calculation of such important quantities as the the specific heat, the
latent heat, etc.

With the use of the derived EFRGE a closed expression for the renormalization
counterterm has been obtained which when subtracted from the divergent
solution of the Wetterich equation would lead to a finite exact expression for
the EA thus making two approaches formally equivalent.

The S-matrix equation has been found to be simply connected with a generalized
functional Burgers’ equation which establishes a direct correspondence between
the first order phase transitions and the shock wave solutions of the RG
equation.

The transparent semi-group structure of the S-matrix RG equation makes
possible the use of different RG techniques at different stages of the RG flow
in order to improve the LPA solution.

Keywords: exact renormalization group equations, local potential
approximation, $n$-vector spin-lattice models, Burgers’ equation, shock wave
solutions, first order phase transitions

##  1 Introduction

Exact functional renormalization group equations (EFRGEs) were introduced by
Wilson [1] as a prospective method of dealing with problems that cannot be
solved by other techniques. An important problem of this kind is the solution
of models with strong coupling between the fluctuating fields for which
rigorous solution methods similar to the perturbative techniques of the weak
coupling case are unavailable. The EFRGE derived in [1], however, was too
complicated to be solved beyond the perturbation theory. Therefore, using the
flexibility of the RG approach, simpler EFRGEs were derived in [2, 3, 4, 5, 6,
7, 8, 9] and used in the solution of various field-theoretic problems (see
review papers [10, 11, 12, 13] and references to earlier literature therein).

To the purposes of the present study the most important will be two kinds of
EFRGEs : those dealing with the renormalization of the interaction functional
suggested in [5, 4] and the equations introduced in [7, 8, 14, 9] which
renormalize the effective action (EA). Our interest to these two approaches
stems from the fact that they proved to be successful in application to the
lattice models of the Ginzburg-Landau type [15, 16, 17, 18, 19]. Moreover,
when applied to the same problem and solved in the local potential
approximation (LPA) [20, 10, 11, 21, 13, 22, 15, 16, 6, 5, 17] EFRGEs of both
kinds produce very similar results provided that the same cutoff in the
momentum space in the form of the step function [5, 15] is being used in
implementing the LPA. This suggests that both kinds of equations may be
equivalent in some sense and so will give similar results in any
approximation.

Another possible explanation is that LPA is very accurate and all EFRGEs would
give similar solutions if solved within this approximation. This, however,
does not seem to be plausible because when LPA is used with RG equations that
do not belong to the above two categories, the results may disagree
significantly as, e.g., between the solutions of the Blume-Capel model in [23]
and in [19]. Besides, a rigorous assessment of the LPA accuracy meets with
serious difficulties. The derivative expansion (DE) that is usually invoked as
a way of justifying and correcting LPA (see extensive bibliography on the
subject in [24] and [13]) is not a systematic approach [24]. For the purposes
of the present study the most important conclusion of [24] was that DE can be
valid only at moderate interaction strengths. This was confirmed in [16] by
comparison of LPA solution using Wilson-type cutoff with the Monte Carlo (MC)
simulations of [25] where the agreement did worsen with the growing
interaction strength. But the spin-lattice models we will be interested in
formally correspond to the case of infinitely strong coupling and though in
the course of the RG flow it considerably weakens, in the transient region [1]
the large values may persist long enough to make the resort to DE ungrounded.
In view of this, the remarkable accuracy of the values of the phase transition
temperatures obtained within LPA in [6, 15, 16, 23, 17, 18, 19] in strongly
coupled lattice models requires explanation.

The aim of the present paper is to clarify the above issues and to develop
techniques facilitating the use of nonperturbative RG methods in the solution
of strongly coupled lattice models. Of special interest to us will be the
spin-lattice case because of the prominent role that the Ising model (IM), the
classical Heisenberg model, the Blume-Capel model and many others play in
statistical physics.

Specifically, in this paper we will derive EFRGE for the renormalization of
the interaction functional [5, 4] using a general cutoff or, equivalently, the
regulator function [13] and will show that a fully renormalized interaction
functional coincides with the generating functional of the logarithm of the
S-matrix with minus sign. Therefore, RG equations of this kind will be called
the S-matrix equations to distinguish them from the equations of [7, 9, 11,
13] which will be called the EA equations. We will show that EFRGEs of the two
kinds are connected by a functional Legendre transform (LT) and thus formally
are equivalent, but with a caveat.

In Wilsonian RG approach it is habitual to neglect the field-independent
multiplicative factors that contribute only to the normalization of the field
probability distribution [1] or, equivalently, to the field independent terms
(f.i.t.) in the renormalized Hamiltonian. According to [9], this practice was
also accepted in the derivation of EA EFRGEs in [7, 9, 11, 13] which means
that only the field-dependent part of the EA can be assuredly correct when
calculated with the use of the EA EFRGE schemes of the above references. In
many cases this is completely sufficient which explains the successful use of
the RG equation of [7, 9, 11, 13] in quantitative description of the equation
of state of magnetic systems [26] and to the calculation of the magnetization
and the critical temperatures in lattice systems [15, 16]. The reason was that
only the field derivatives of EA were needed in the calculations and so the
f.i.t. disappeared from all expressions.

However, in other physically important cases the knowledge of the f.i.t. is
indispensable. For example, the Onsager solution of the 2D IM is a f.i.t.; the
specific heat and the latent heat of the first order phase transitions (FOPT)
(see, e.g., [19]) both depend on f.i.t..

An advantage of the EA EFRGE derivation by the LT of the S-matrix equation is
that no normalization factors have been omitted which has made possible to
correct the RG approaches of [7, 9, 11, 13] in such a way that they treated
the f.i.t. exactly. This may have important consequences for the theories
based on the Wetterich equation in which f.i.t. at the lower limit of
integration is divergent [7, 9, 11]. To regularize the divergence various ways
of choosing the ultraviolet cutoff have been suggested in [27, 28, 29, 30, 31]
which, however, relied on heuristic ad hoc assumptions which are difficult to
rigorously justify. To overcome this difficulty, by analogy with the
renormalization in the relativistic quantum field theories a divergent
counterterm will be inserted into the initial condition of the Wetterich
equation which will render f.i.t. in the EA to be finite and unambiguously
defined.

Besides the latent heat, discontinuities in the order parameter are important
characteristics of FOPT. The analytic structure of the S-matrix EFRGE can be
helpful in their description because when differentiated w.r.t. the field
variable the equation may be qualitatively described as the generalized
functional Burgers’ equation (BE) [32]. In the LPA it can be used to describe
FOPT as the shock wave solutions of the $n$-vector Ginzburg-Landau model for
all natural $n$ in contrast to the EA equation where the shock waves were
found only in the $n\to\infty$ limit [33]. The shock-wave picture puts on a
firmer footing the mechanism of the suppression of the van der Waals loops
plaguing the mean field (MF) theories discovered in the LPA RG approach in [6]
.

In the present study of major interest will be the RG description of non-
universal quantities, such as the phase transition temperatures, because of
their truly non-perturbative nature. Universal quantities, such as the
critical exponents, though not accurately predicted by the LPA do not depend
on the interaction strength and in the three-dimensional systems that we will
consider in this paper can be efficiently calculated by perturbative
techniques within simplified Hamiltonians characterized by a few parameters
[1, 34]. In contrast, non-universal quantities may depend on an unlimited
number of parameters and are particularly difficult to calculate beyond the
perturbation theory.

Obviously, the universal quantities are also needed for a comprehensive
description of critical phenomena in strongly interacting systems so it would
be reasonable to accommodate the available rigorous results [1, 34] to the LPA
solutions. To this end in the present paper a multi-step renormalization
technique will be developed that will make possible the use of different
renormalization methods at different stages of the RG flow. For example, it
will be possible to use LPA within the transient region at an early stage of
renormalization when the interactions are strong [1] and to switch to the
perturbative treatment in the critical region where the conventional
perturbation theory or the DE expansion become efficient [24, 35, 36].

The presentation of material will be as follows. After introducing the
necessary notation in the next section, in section 3 we will derive the EFRGE
for the interaction functional that underlies the LPA equations of [5, 18, 17,
19]. It will be shown that the RG flow interpolates between the initial local
potential and the generating functional of the logarithm of the S-matrix. It
will be argued that in conjunction with a self-consistency (SC) condition the
S-matrix functional is more local than the FE functional of the Wilson
approach [1] so the S-matrix EFRGE should be more accurately represented in
the LPA than the Wilson one [1]. In B several forms of a LT will be introduced
that connect the S-matrix and the EA RG equations both the exact and in the
LPA form including the case of the multi-step renormalization. Because EFRGEs
are the evolution equations, their solutions depend on the initial conditions.
In sections 5.1, 5.2 and F it will be shown that for S-matrix EFRGE analytic
initial conditions can be easily established even in the spin-lattice case
despite that it formally corresponds to an infinitely strong coupling. The
renormalization of the divergence in the Wetterich equation will be discussed
in B.1.

In section 3.4 we will show that the Hori representation makes the (semi)group
structure of the S-matrix EFRGE fully transparent which simplifies the
stepwise renormalization needed in the preliminary exact renormalization in
the spin models and also allows one to switch to a perturbative technique in
the critical region.

Further, in section 4 qualitative arguments will be given in support of the
standpoint that the use of the step function cutoff for the elimination of
high momenta [5, 15, 16, 18, 17, 19] should be more accurate in the lattice
case than in the continuum models. The FOPT will be discussed in section 6 and
in D.

In final section 7 a brief summary will be given.

##  2 Definitions and notation

To simplify notation and to facilitate comparison with [1] all derivations in
the present paper will be done for the Landau-Ginzburg model of the Ising
universality class, that is, with the scalar field variable. Generalization to
the $n$-vector models can be achieved along the same lines as in [5, 18] and
in simple cases is straightforward (see section 5.2 below).

Thus, we will be interested in the calculation of the partition functional

$Z(\vec{h})=\int d\vec{s}\,\rme^{-H_{0}(\vec{s})+\vec{h}^{\dagger}\vec{s}}$
(1)

by means of the $N$-dimensional integral over the fluctuating scalar field
which will be usually denoted as a column vector $\vec{s}=[s_{i}]$ (similarly,
$\vec{h}=[h_{i}]$) defined at sites $i$ of a periodic lattice of size $N$;
$\int d\vec{s}\equiv\prod_{i}\int_{-\infty}^{\infty}ds_{i}$. The factor
$1/k_{B}T$ will be assumed to be included in the parameters of the
dimensionless Hamiltonian of the Landau-Ginzburg model

$H_{0}(\vec{s})=\frac{1}{2}\sum_{ij}\left(\epsilon_{ij}+r\delta_{ij}\right)s_{i}s_{j}+U_{0}(\vec{s})=\frac{1}{2}\vec{s}\,^{\dagger}(\hat{\epsilon}+r\hat{I})\vec{s}+U_{0}(\vec{s})$
(2)

where $\hat{\epsilon}=[\epsilon_{ij}]$ is the matrix of the pair interactions
which in the thermodynamic limit $N\to\infty$ (which will be always implicitly
assumed) is translationally invariant; the strength of the pair interactions
will be characterized by a dimensionless constant $K$ with
$\hat{\epsilon}\propto K$; $\hat{I}$ is the unit matrix, $U_{0}$ the
interaction potential and $\vec{h}$ is the source and/or the external field.

The dagger symbol used in (2) denotes Hermitian conjugation which for real
fields means the transposition but for the Fourier transformed fields the
complex conjugation should be additionally applied.

For convenience we will use The symmetric Fourier transform, e.g., in

$s_{\bf k}=N^{-1/2}\sum_{j}\rme^{-\rmi{\bf k}\cdot
j}s_{j}=\sum_{j}\left(\hat{{\cal F}}\right)_{{\bf k},j}s_{j};$ (3)

the use of the unitary matrix $\hat{{\cal F}}$ ($\hat{{\cal
F}}^{-1}=\hat{{\cal F}}^{\dagger}$) will allow us to formally consider the
real-space vector and its Fourier transform as the same abstract vector in two
reference frames connected by a unitary rotation. This makes possible the use
of the same symbol for vectors both with the lattice and with the momentum
components.

The bare interaction $U_{0}$ in (2) can in principle be any functional of
$\vec{s}$, but in the present paper we will assume the interaction functional
to be the sum of site-local potentials

$U_{0}(\vec{s})=\sum_{i}u(s_{i},0)$ (4)

where the second argument of $u$ is $t$,—a scalar variable that will
parametrise the RG flow. It can be chosen arbitrarily so we define it to be
varying from $t=0$ corresponding to the system in its initial or “bare” state
to $t=t^{\mathrm{R}}$ in the fully renormalized state. The bare potential in
(4) will be assumed to be analytic in the field variable but as we will see
later the analyticity of the renormalized potential at $t=t^{\mathrm{R}}$ can
be broken by phase transitions.

Here and below by superscript “R” we will mark all fully renormalized
quantities that appear also in only partially renormalized form at arbitrary
$t$.

Further, because the separation of Hamiltonian in the quadratic and the
interaction parts is not unique, the ambiguity has been used to define the
first term in (2) in such a way that, first, the Fourier transform of
$\hat{\epsilon}$ had the small-momentum asymptotic

$\epsilon({\bf k})_{k\to 0}\simeq ck^{2}$ (5)

where $k=|{\bf k}|$ and $c$ is a constant proportional to $K$. This can be
achieved by subtracting an appropriate site-diagonal term in the quadratic
part of (2) and adding it to $U_{0}$. Second, an arbitrary diagonal term
$r\hat{I}$ has been added to the quadratic part and subtracted from $U_{0}$.
Obviously, neither $H_{0}$ nor the exact partition function (1) depend on $r$
but in an approximate solution the independence may be broken and the
arbitrary parameter can be used to improve the accuracy. This will be done
below with the use of a SC condition (section 3.2).

The quantities that we are going to calculate in this paper will be the
Helmholtz FE

$F^{\mathrm{R}}(\vec{h})=-\ln Z(\vec{h})$ (6)

and its two derivatives: the magnetisation

${m}_{i}\equiv\langle s_{i}\rangle=-\frac{\partial
F^{\mathrm{R}}({\vec{h}})}{\partial{h_{i}}},$ (7)

and the pair correlation function

$G_{ij}^{\mathrm{R}}=-\frac{\partial
F^{\mathrm{R}}({\vec{h}})}{\partial{h_{i}}\partial{h_{j}}}=\langle{s}_{i}{s}_{j}\rangle-{m}_{i}{m}_{j}.$
(8)

Throughout the paper the arrows will denote the $N$-dimensional lattice
vectors; Fourier momenta and $n$-vectors will be boldface and for simplicity
the same letters but without arrows or subscripts, such as $h$ and $m$, will
denote the scalar values of the homogeneous external field and of the on-site
magnetisation, respectively. In particular, according to (3) the homogeneous
external field will have vector components

$h_{\bf k}^{h}=\sqrt{N}\delta_{\bf k,0}h.$ (9)

##  3 Exact RG equations

In the derivation of S-matrix EFRGE we will use the Wilson Green’s function
approach (see [1], ch. 11) which will be slightly modified in order to
facilitate comparison with the EA approach as presented in review articles
[11, 13]. To begin with let us consider a$N$-dimensional Gaussian kernel

${\cal
G}(\vec{s},\vec{s}\,^{\prime},t)={\det{[\hat{{R}}(t)}/{2\pi}]}^{\frac{1}{2}}\exp\left[-\frac{1}{2}\left(\vec{s}-\vec{s}\,^{\prime}\right)^{\scriptscriptstyle
T}\hat{{R}}(t)\left(\vec{s}-\vec{s}\,^{\prime}\right)\right]$ (10)

where $\hat{R}(t)$ is the regulator matrix [11, 13] which is assumed to be a
non-negative definite, symmetric, translationally invariant and that its site-
diagonal matrix elements behave at $t\to 0$ as

$\hat{{R}}(t\to 0)\to\infty$ (11)

in order to satisfy the conventional initial condition

${\cal
G}(\vec{s},\vec{s}\,^{\prime},t=0)=\delta(\vec{s}-\vec{s}\,^{\prime})\equiv\prod_{i}\delta(s_{i}-s_{i}^{\prime}).$
(12)

In the Wilson approach the regulator $\hat{R}$ in (10) should be chosen in
such a way that the convolution of the Green function with the Boltzmann
factor $\rme^{-H_{0}+\vec{h}^{\dagger}\vec{s}}$ at $t^{\mathrm{R}}$ satisfied

$Z(\vec{h})=C\int d\vec{s}\,^{\prime}{\cal
G}(\vec{s},\vec{s}\,^{\prime},t^{\mathrm{R}})\rme^{-H_{0}(\vec{s}\,^{\prime})+\vec{h}^{\dagger}\vec{s}\,^{\prime}}$
(13)

with some proportionality coefficient $C$ which may depend on $\vec{s}$,
$\vec{h}$, $t^{\mathrm{R}}$, etc., with the only requirement that it could be
calculated exactly.

Replacing $t^{\mathrm{R}}$ in (13) by $t\in(0,t^{\mathrm{R}})$ one obtains an
expression that interpolates between the Boltzmann factor and the fully
renormalized partition function. Because the Green function satisfies a
diffusion-type differential equation, by choosing a suitable form of the
regulator one can derive a Wilsonian EFRGE in differential form, as shown in
[1] with a particular choice of the regulator. In this paper we will derive
along similar lines a different, S-matrix EFRGE, adopted for the use in
conjunction with the LPA.

But before proceeding it is pertinent to point out that similar to the Wilson
equation [1], the S-matrix EFRGE deals with the calculation of the Helmholtz
FE. On the other hand, the EA in application to statistical-mechanical
problems aims at calculating the Gibbs FE which necessitated the use of a
modified LT in the derivation of EA EFRGE [13]. Of course, the two fully
renormalized thermodynamic potentials are connected by the conventional LT
[37].

###  3.1 S-matrix EFRGE

The advantage of S-matrix EFRGE derived below in comparison with the Wilson-
type RG equations [1, 2, 3] and with EA EFRGEs of [7, 8, 14, 9] is that from a
qualitative standpoint they are not suitable for the LPA in lattice models
even with local interaction potentials. In both cases EFRGEs deal with
renormalization of the total Hamiltonian which contains both the local
potential and the non-local quadratic part. Because EFRGEs are essentially
nonlinear, in the course of the RG flow the two terms are intermixed which, in
particular, spreads the interaction potential over many lattice sites thus
making the LPA along the whole RG trajectory dubious from a qualitative point
of view.

In such cases LPA can be justified only formally within the DE [24, 13].
However, for DE to converge the effective expansion parameter should be at
least smaller than unity. This is possible in the continuum models with a
single cutoff parameter $\Lambda$. The lowest order correction in the
dimensionless DE parameter in the isotropic systems will be of order
$(k/\Lambda)^{2}\leq 1$. In the lattice case, however, the cutoff is
determined by the lattice constant $a_{\mathrm{l}at}$ but because in the
corners of the Brillouin zone (BZ) of, e.g., sc lattice the squared momentum
can be as large as $(ka_{\mathrm{l}at})^{2}\approx 3\pi^{2}\gg 1$ the use of
DE becomes meaningless at the early stages of the RG flow. (Henceforth we
chose $a_{\mathrm{l}at}$ to be our length unit.)

Therefore, to justify LPA in lattice models at least qualitatively it seems
reasonable to derive an EFRGE for the interaction part of the Hamiltonian and
develop approximations that would keep the renormalized potential during the
RG flow as local as possible. To this end let us consider functional

$S(\vec{s},t)=\rme^{-U(\vec{s},t)}$ (14)

subject to the initial condition $U(\vec{s},t=0)=U_{0}(\vec{s})$. In the
Wilson Green’s function approach the RG evolution of $S$ is described by the
equation

$S(\vec{s},t)=\int{\cal
G}(\vec{s},\vec{s}\,^{\prime},t)\rme^{-U_{0}(\vec{s}\,^{\prime})}\,d\vec{s}\,^{\prime}.$
(15)

From explicit form of ${\cal G}$ (10) it is easy to surmise that the partition
function can be obtained from the fully renormalized $S^{\mathrm{R}}$ or,
equivalently, from $U^{\mathrm{R}}$ if instead of (70) we impose on the
regulator the following condition

$\hat{{R}}^{\mathrm{R}}=\hat{R}(t^{\mathrm{R}})=\hat{\epsilon}+r\hat{I}.$ (16)

With this regulator one gets at $t=t^{\mathrm{R}}$ after expanding the
quadratic form in (10)

$\rme^{-U^{\mathrm{R}}(\vec{s})}=\sqrt{\det{\left(\hat{{R}}^{\mathrm{R}}/{2\pi}\right)}}\exp\left(-\frac{1}{2}\vec{s}\hat{{R}}^{\mathrm{R}}\vec{s}\right)Z(\hat{{R}}^{\mathrm{R}}\vec{s}).$
(17)

Next choosing $\vec{s}$ to satisfy

$\hat{{R}}^{\mathrm{R}}\vec{s}=\vec{h}$ (18)

we obtain the partition functional

$Z(\vec{h})=\sqrt{\det{\left({2\pi}\hat{{G}}^{\mathrm{P}}\right)}}\exp\left(\frac{1}{2}\vec{h}\hat{{G}}^{\mathrm{P}}\vec{h}\right)\rme^{-U^{\mathrm{R}}(\hat{{G}}^{\mathrm{P}}\vec{h})}$
(19)

where

$\hat{{G}}^{\mathrm{P}}=\left(\hat{{R}}^{\mathrm{R}}\right)^{-1}=\left(\hat{\epsilon}+r\hat{I}\right)^{-1}$
(20)

is the bare propagator. The superscript “P” has been used instead of “R” to
distinguish $\hat{{G}}^{\mathrm{P}}$ from the fully renormalized pair
correlation function $G^{\mathrm{R}}_{ij}$ defined in (8).

The FE functional (6) corresponding to (19) reads

$F^{\mathrm{R}}(\vec{h})=U^{\mathrm{R}}(\hat{G}^{\mathrm{P}}\vec{h})-\frac{1}{2}\vec{h}^{\dagger}\hat{{G}}^{\mathrm{P}}\vec{h}-\frac{1}{2}\Tr\ln{\left(2\pi\hat{{G}}^{\mathrm{P}}\right)}.$
(21)

###  3.2 Self-consistency condition

Standard perturbative analysis (see, e.g., [1, 38]) shows that renormalized
interaction vertices in $F^{\mathrm{R}}(\vec{h})$ (21) are the connected
correlation functions coupled to the source field $\vec{h}$ with the external
legs in the diagrammatic representation corresponding to the exact pair
correlation functions $\hat{G}^{\mathrm{R}}$. Here by interaction vertices we
mean the terms in $F^{\mathrm{R}}$ of the third and higher orders in the
source field which obviously are all contained in the renormalized potential
$U^{\mathrm{R}}$. But according to [1] (see, e.g., equation (7.45)), factors
$\hat{G}^{\mathrm{R}}$ describe the leading long-distance behaviour at the
critical point when the nonlocality is the most pronounced. Thus, if instead
of $\hat{G}^{\mathrm{P}}$ in the argument of $U^{\mathrm{R}}$ in (21) we had
$\hat{G}^{\mathrm{R}}$ then functional $U^{\mathrm{R}}(\vec{s})$ would be more
local than $F^{\mathrm{R}}$ and so better representable by the LPA.

Obviously that with only one parameter $r$ at hand it is in general not
possible to make $\hat{G}^{\mathrm{P}}$ to be equal to $\hat{G}^{\mathrm{R}}$
because in the exact pair correlation function the mass operator will not be
local in nontrivial models. However, in ferromagnetic systems the long-
distance behaviour is governed by the smallest momenta in the Fourier
transformed pair correlation function so for our purposes it would be
sufficient to satisfy a simpler requirement

$G^{\mathrm{R}}({\bf k\to 0})\approx G^{\mathrm{P}}({\bf k\to 0})=1/r$ (22)

where according to (20)

$G^{\mathrm{P}}({\bf k})=\frac{1}{\epsilon({\bf k})+r}.$ (23)

It should be pointed out that here we have neglected the modification of the
long-distance behaviour due to the critical exponent $\eta$. It can be
expected, however, that it is important only when $r\to 0$. Besides, the
exponent is equal to zero in the LPA so anyway we could not account for it in
this approximation. A possibility to correct this within our RG approach will
be discussed in section 5.3.

Now substituting (21) into (8) one finds that condition (22) in terms of the
fully renormalized potential reads

$\left.\frac{\partial^{2}U^{\mathrm{R}}(\vec{s})}{\partial{s_{\bf
k}}\partial{s_{\bf-k}}}\right|_{{\bf k\to 0},s_{\bf k=0}}=0$ (24)

where we assume that the solution has been obtained in a translationally-
invariant system in homogeneous external field $h$ so that
$\vec{s}=\hat{G}^{\mathrm{P}}\vec{h}$ in (21) in this case should be
calculated with the use of (9) and (23) as

$s_{\bf k=0}=\left.G^{\mathrm{P}}({\bf k})h^{h}_{\bf k}\right|_{\bf
k=0}=\sqrt{N}h/r$ (25)

The SC condition (24) will be used everywhere below in the solutions of the
LPA RG equations.

###  3.3 S-matrix EFRGE in differential form

Differentiating (15) w.r.t. $t$, using in the derivative of ${\cal G}$ (10)
the commutativity of translationally-invariant matrices and Jacobi’s formula
for invertible matrices one arrives at the exact linear RG equation

$\partial_{t}S(\vec{s},t)=\frac{1}{2}\sum_{ij}\partial_{t}{{G}}_{ij}\frac{\partial^{2}S(\vec{s},t)}{\partial{s_{i}}\partial{s_{j}}}$
(26)

where $\partial_{t}\equiv\partial/\partial t$ and

$\hat{G}(t)=[G_{ij}(t)]=\hat{{R}}^{-1}(t).$ (27)

Substituting $S$ from (14) in (26) the non-linear EFRGE for the interaction
potential in the lattice coordinates is obtained as

$\partial_{t}{U}=\frac{1}{2}\sum_{ij}{\partial_{t}{G}}_{ij}^{(0)}\left(\frac{\partial^{2}U}{\partial{s_{i}}\partial{s_{j}}}-\frac{\partial
U}{\partial{s_{i}}}\frac{\partial U}{\partial{s_{j}}}\right).$ (28)

But in implementation of the Wilsonian renormalization we will need it in the
momentum representation [5, 4]

$\partial_{t}U(\vec{s},t)+\frac{1}{2}\sum_{\bf k}\partial_{t}G({\bf
k},t)\frac{\partial U}{\partial s_{\bf-k}}\frac{\partial U}{\partial s_{\bf
k}}=\frac{1}{2}\sum_{\bf k}\partial_{t}G({\bf
k},t)\frac{\partial^{2}U}{\partial s_{\bf-k}\partial s_{\bf k}}.$ (29)

Here the terms have been rearranged in such a way that after differentiation
w.r.t. $s_{\bf q}$ (29) acquired the structure of a $N$-dimensional
generalized BE [32, 39] for $\vec{\cal U}_{s}=\left[\partial U/\partial s_{\bf
q}\right]$. It will greatly simplify under LPA and in this form will be useful
in the description of FOPTs in section 6.

###  3.4 Feynman’s diffusion equation as an EFRGE

In our notation Feynman’s diffusion equation [40] is obtained from (26) by
assuming that $G_{ij}(t)=tG^{\mathrm{P}}_{ij}$, $0\leq t\leq
t^{\mathrm{R}}=1$:

$\partial_{t}S(\vec{s},t)=\frac{1}{2}\sum_{ij}{G}^{\mathrm{P}}_{ij}\frac{\partial^{2}S(\vec{s},t)}{\partial{s_{i}}\partial{s_{j}}}.$
(30)

Thus, (30) is an EFRGE with a particular choice of the cutoff matrix
$\partial_{t}\hat{G}(t)$.

The Hori representation [40, 38, 5] can be obtained by formally integrating
(30) or, more generally, (26) as

$S(\vec{s},t)=\exp\left(\frac{1}{2}\sum_{ij}{G}_{ij}(t)\frac{\partial^{2}}{\partial{s_{i}}\partial{s_{j}}}\right)S_{0}(\vec{s}).$
(31)

It can also be obtained independently of Feynman’s equation [38] so the latter
as well as the S-matrix EFRGE (29) can be straightforwardly derived from (31)
by simple differentiation w.r.t. $t$ [5, 17].

An important property of the Hori representation is that its semi-group RG
structure becomes completely transparent if $\hat{G}(t)$ in the exponential is
written as $\int_{0}^{t}\partial_{t^{\prime}}\hat{G}(t^{\prime})dt^{\prime}.$
Now from the properties of the integral it is seen that the renormalization
can be performed in several finite steps, for example, first from $t=0$ to
$t=t_{0}$ and then from $t_{0}$ to $t$. At the second step $S(\vec{s},t_{0})$
should be taken as the initial condition and $\hat{G}(t)$ in the exponential
should be replaced by

$\hat{\Delta}(t,t_{0})=\hat{G}(t)-\hat{G}(t_{0}).$ (32)

This possibility will be used in sections 4 and 5.3 below. The multi-step
renormalization can be also obtained in the Green’s function approach (see A),
though not as straightforwardly as in the Hori representation.

###  3.5 Legendre transform of S-matrix EFRGE

In B explicit LTs of the S-matrix equation have been performed both exactly
and in the LPA. It has been shown that the EA equation (69) from [7, 9, 11,
13] can be obtained from (29) by neglecting f.i.t. in the exact equation. We
use somewhat different notation in these equations because it is not clear
whether our LT in B and the modified LT in [7, 9, 11, 13] are the same at
intermediate stages of the RG flow. But the equations are mathematically
equivalent so with the same initial conditions will lead to the same
solutions.

The EA EFRGE (68) takes into account all contributions and so can be used in
statistical-mechanical calculations. It is to be noted that because with
condition (16) $G({\bf k},t)$ in (68) is not divergent neither at $t=0$ nor at
$t^{\mathrm{R}}$ which makes this EFRGE easier to deal with computationally
than equations (69) and (81). It is this EA EFRGE in LPA that was used in all
calculations in [17, 18, 19]. An equation similar to (68) was obtained in [9]
but in the derivation the normalization factors and the vacuum contributions
that correspond to f.i.t. in FE case have been dropped.

##  4 LPA

The step-function cutoff suggested in [5] for continuum models proved to be
very accurate in application to LPA solutions of lattice models [6, 15, 16,
23, 17, 18, 19]. This success can be qualitatively understood within the
Kadanoff picture of critical phenomena [41, 1] as follows.

Formally, LPA consists in assuming that functional $U(\vec{s},t)$ preserves
the local structure of $U_{0}$ (4) throughout the whole evolution from $t=0$
to $t=t^{\mathrm{R}}$:

$\displaystyle U(\vec{s},t)\stackrel{{\scriptstyle
LPA}}{{\approx}}N\sum_{l,\\{{\bf k}_{i}\\}}N^{-\frac{l}{2}}u_{l}(t)s_{{\bf
k}_{1}}s_{{\bf k}_{2}}\dots s_{{\bf k}_{l}}\delta^{L}_{{{\bf k}_{1}}+{{\bf
k}_{2}}+\dots+{{\bf k}_{l}}},$ (33)

i.e., the Fourier-transformed coefficients of the expansion do not depend on
the momenta apart from the lattice Kronecker symbols $\delta^{L}$. The latter
differs from the conventional Kronecker delta in that the momenta may sum not
only to zero but to any of the reciprocal lattice vectors, as can be seen from
the formal definition (103).

The LPA ansatz (33) allows one to establish a one-to-one correspondence
between functional $U$ and (in the Ising universality class) a function of a
real variable $x$

$u(x,t)=\sum_{l=0}^{\infty}u_{l}(t)x^{l}$ (34)

which can be interpreted as the on-site potential. The general $n$-vector case
is treated similarly with using $n$-vector ${\bf
x}=(x_{1},x_{2},\cdots,x_{n})$ instead of $x$ [5].

The evolution equation for $u$ should be obtained by substituting (33) into
(29). However, in general case the equation will be incompatible with the LPA
ansatz so the task is to satisfy it approximately. The difficulty causes the
momentum dependence of the cutoff function $\partial_{t}G({\bf k},t)$. Indeed,
if it was equal to unity

$\partial_{t}G=1,$ (35)

then equation (29) could be satisfied exactly because in the real space the
inverse Fourier transform of unity is $\delta_{ij}$ which would make the RG
equation separable, that is, representable as the sum over the sites of
identical local potentials. Feynman’s diffusion equation (30) in this case
would admit solution in quadratures which in terms of $u$ would read

$\rme^{-u({x},t)}=\frac{1}{\sqrt{2\pi
t}}\int_{-\infty}^{\infty}d{x}_{0}\;\rme^{-({x-x}_{0})^{2}/{2t}}\rme^{-u(x_{0},0)}.$
(36)

The layer-cake renormalization scheme [5, 17, 18] has been devised with the
aim of emulating (35)–(36) with a maximum possible accuracy in the case of
momentum-dependent cutoff function. To this end the latter was chosen as

$\partial_{t}G({\bf k},t)=\theta[G^{\mathrm{P}}({\bf k})-t]$ (37)

where $\theta$ is the step function that is equal to unity when the Fourier
momentum ${\bf k}$ is inside the region $\Omega(t)$ defined by the condition

$G^{\mathrm{P}}({\bf k}\in{\Omega(t)})\geq t$ (38)

(see figure 1) and according to (37) $\partial_{t}G=0$ outside of $\Omega$.
For the purposes of the present paper it will be sufficient to keep in the
fully renormalized quantities only the zero-momentum component $s_{\bf k=0}$.
Therefore, in the course of the RG flow the components outside $\Omega$ can be
set equal to zero [1]: $s_{\bf k}|_{{\bf k}\in\bar{\Omega}}=0$ where
$\bar{\Omega}=\mbox{BZ}\backslash\Omega$ is the part of BZ external to
$\Omega$. Next one may rescale the momenta in (33) so that $\Omega$ regained
the same volume as the original BZ and the rescaled system could be considered
as describing the same lattice model approximately renormalized by the
Kadanoff real-space procedure [41].

Figure 1:  Illustration of the layer-cake renormalization scheme. The stacked
narrow rectangles are cross-sections of $dG=\theta[G^{\mathrm{P}}({\bf
k})-t]dt$ along the diagonal of the Brillouin zone on which
$G^{\mathrm{P}}({\bf k})$ reaches its minimum value $t_{0}$ defined in (47);
$\bar{t}=t-t_{0}$; $\Lambda$ is an approximate cutoff momentum in the region
where $G^{\mathrm{P}}({\bf k})$ is almost isotropic (for further explanations
see the text).

Further details on the rescaling can be found in [5] but in this paper the RG
LPA equation of [5] in the scaling form will be used only for the calculation
of critical exponents. For the calculation of non-universal quantities,
however, the rescaling has been omitted because it introduces a large but
trivial Lyapunov exponent equal to the spatial dimension of the system which
makes the RG equation numerically unstable. Instead, we simply substitute (33)
in (29) and observe that the second derivative on the r.h.s. makes all terms
in (33) independent of ${\bf k}$ so they all acquire the common factor [17]

$p(t)=\frac{1}{N}\sum_{\bf k}\theta\left[G({\bf
k})-t\right]=\int_{0}^{t^{-1}-r}dE\rho(E)$ (39)

where $\rho(E)$ is the density of states corresponding to the dispersion
$\epsilon({\bf k})$ [15, 18, 17].

In the second term on the l.h.s. of (29), however, the summation over ${\bf
k}$ is lifted by one of the two Kronecker deltas so the dependence on the
momentum remains in the cutoff function (37). In the LPA we neglect it by
assuming that in the critical region the most important are small momenta for
which (35) holds. In E it has been argued that this approximation should be
better fulfilled in the lattice case than in the continuum Landau-Ginzburg
model [5]. Thus, the LPA RG equation for the local potential (34) in the Ising
universality class reads

$u_{t}+\frac{1}{2}u_{x}^{2}=\frac{1}{2}p(t)u_{xx}.$ (40)

To make notation more compact, here we have made use of the fact that in the
Ising universality class in the LPA there is no need for the vector component
subscripts so henceforth they will be used to denote partial derivatives. In
this notation the free energy (21) per site in LPA reads

$f(h)=u^{\mathrm{R}}(h/r)-h^{2}/2r-\frac{1}{2N}\Tr\ln{\left(2\pi\hat{{G}}^{\mathrm{P}}\right)}$
(41)

where use has been made of (19) and the argument of $u^{\mathrm{R}}$ has been
found by substituting (25) in (33) at the end of renormalization when only
${\bf k=0}$ component contributes to $U^{\mathrm{R}}$ which means that at
$t^{\mathrm{R}}$ the argument of $u^{\mathrm{R}}$ is

$x=h/r.$ (42)

The equation of state is obtained as

$m(h)=-\frac{df}{dh}=\frac{h}{r}-\frac{1}{r}u^{\mathrm{R}}_{x}|_{x=h/r}$ (43)

and the SC condition (24) in the LPA takes the form

$\left.u^{\mathrm{R}}_{xx}\right|_{x=h/r}=0.$ (44)

The general $n$-vector case can be treated similarly [5] and in fully
$O(n)$-symmetric case the LPA equation simplifies to [18]

$u_{t}+\frac{1}{2}(\nabla u)^{2}=\frac{1}{2}p(t)\nabla^{2}u$ (45)

where the differential operators act in the space of $n$-vectors ${\bf x}$.

##  5 Illustrative calculations

###  5.1 Critical temperatures in $\varphi^{4}$ model

The first example we consider is the $\varphi^{4}$ model on sc lattice [25,
16]. In notation of section 2 the interaction potential reads

$U_{0}(\vec{s})=\sum_{i}\left[(1-3K-r/2)s_{i}^{2}+\lambda(s_{i}^{2}-1)^{2}\right]$
(46)

(cf. equation (2) in [25]). The LPA solution has been obtained as follows.
First $U_{0}$ in (46) has been mapped onto $u(x,0)$ by replacing $s_{i}$ with
$x$ in the summand. Then equation (40) has been solved numerically for $h=0$
(the critical value of the field) and iterated to satisfy the SC condition
(44) (for details of the procedure see [17]) to determine $r$ corresponding to
chosen parameters $K$ and $\lambda$. The critical value $K_{c}$ corresponds to
$r=0$. Calculations within the method of lines have been performed for 11
values of $\lambda$ in the range 0.1—2.5 that were used in MC simulations in
[25]. The calculated $K_{c}$ agreed with the MC values with accuracy better
than 0.25% similar to the LPA calculations in [16] based on EA EFRGE [7, 9,
11, 13] and the step function cutoff. Moreover, even the systematic change
with $\lambda$ of the discrepancy sign has been similar in both calculations,
presumably, as a consequence of equivalence of the exact equations and of the
same cutoff used in the LPA in [16] and in the present paper.

###  5.2 LPA for spin-lattice models

The $O(n)$-symmetric spin-lattice models in principle can be solved in the
same way as the $\varphi^{4}$ model, only using equations (45) instead of
(40). To this end it would be sufficient to replace in (46) $s_{i}$ by
$n$-dimensional vectors ${\bf s}_{i}$ and to let $\lambda\to\infty$ to
suppress fluctuations of the vector length. However, numerical solution of
(45) with infinite initial condition would be technically very difficult. But
in the partition function (1) the $\lambda\to\infty$ limit amounts to
appearance in the integrand of the factor $\prod_{i}\delta({\bf s}_{i}^{2}-1)$
[42] which is site-local and despite being singular it can be exactly
renormalized with the use of (36) which has been generalized to the $O(n)$
case in F. The exact renormalization has been done in the interval (see figure
1)

$0\leq t\leq t_{0}=\min_{\bf k}G^{\mathrm{P}}({\bf k})=\frac{1}{r+\max_{\bf
k}\epsilon({\bf k})}$ (47)

where (35) is satisfied exactly because the integral in (39) saturates to
unity above the upper limit of the dispersion. Explicit expressions for
$u^{(n)}({\bf x},t_{0})$ from (109) can be used as the initial conditions for
LPA RG equation (45). In the symmetric phase the equation depends only on the
radial coordinate and its numerical integration for $n>1$ has been as
straightforward as in $n=1$ case [18]. The critical temperatures in table 1
have been found as the points where the SC values of $r$ interpolated to zero.
The accuracy similar to that in table 1 was also obtained in RG calculations
for IM [23] and additionally for $XY$\- and the Heisenberg models on the cubic
lattice in [15]. In the later case, however, it was found to be necessary to
resort to a non-functional technique and to introduce heuristic modifications
in the partially renormalized EA which makes the good accuracy obtained less
convincing than in the present approach.

Table 1: Dimensionless inverse critical temperatures of the $n$-vector spin models on cubic lattices calculated in the LPA. The errors have been estimated by comparison with the MC simulations data [43] for $n=1$ (IM) and with the high temperature expansion for $n=2$, 3 and 4 [42]. $n$ | Lattice | $K_{c}$ | Error
---|---|---|---
1 | fcc | 0.1023 | 0.2%
1 | bcc | 0.1579 | 0.3%
1 | sc | 0.2235 | 0.8%
2 | bcc | 0.3225 | 0.6%
2 | sc | 0.4597 | 1.2%
3 | bcc | 0.4905 | 0.8%
3 | sc | 0.7025 | 1.4%
4 | bcc | 0.6608 | 0.8%
4 | sc | 0.9488 | 1.4%

As is seen, the values of critical temperatures have been calculated with the
accuracy better than 2% which should be sufficient for most practical purposes
because in realistic lattice models microscopic Hamiltonians are rarely known
with better accuracy.

###  5.3 Critical exponents

The critical exponents have been calculated with the use of the scaling form
of (45) [5]. Because in the LPA there is only one independent exponent, for
comparison purposes exponent $\nu$ has been used. The values found were
0.65(0.63) for $n=1$, 0.71(0.67) for $n=2$, 0.76(0.71) for $n=3$ and
0.80(0.75) for $n=4$; in parentheses are shown the rounded values from precise
calculations taken from [12]. The LPA values were closest to those calculated
in [11] in the lowest order of the derivative expansion, though being
systematically larger on $\sim 0.01$. This apparently is a consequence of a
similar non-perturbative RG approach used by the authors while the difference
could be attributed to the fact that $\eta$ in the calculations in the above
reference was not equal to zero but obtained from RG equations.

Though not very large, the errors in the LPA values of critical exponents can
be seen in the discrepancy between experimental and theoretical curves in the
disordered phase in [17] and can be even larger in other experiments, e.g., on
the specific heat where exponent $\alpha$ is small and the LPA error is about
100% in the Ising universality class.

Going beyond LPA within lattice nonperturbative RG, however, is not an easy
task mainly because of the large values of the momenta within BZ which
precludes the use of DE; seemingly more appropriate expansion in circular
harmonics, on the other hand, worsens the accuracy of critical temperatures
which is undesirable [44, 15]. From a practical standpoint it seems reasonable
to exploit the fact that in the critical region the interactions in the course
of the RG flow diminish to the point where the perturbative treatment becomes
justifiable [1]. This makes possible to invoke the multi-step renormalization
technique to switch from LPA to a perturbative treatment in a consistent
manner at an appropriately chosen late-stage point $t_{1}\lesssim
t^{\mathrm{R}}$ (see figure 1). Using the Hori representation (31) the usual
diagrammatic techniques can be used in the calculations with the only
modification that the conventional propagator that behaves at small $k$ as
$G^{\mathrm{P}}(k)\simeq 1/(ck^{2}+r)$ (see (5)) should be replaced according
to (32) by

$\Delta(k,t^{\mathrm{R}},t_{1})\simeq\frac{1}{ck^{2}+r}-\frac{1}{c\Lambda^{2}+r}$
(48)

where $\Lambda$ is the momentum cutoff corresponding to $t_{1}$ (see
discussion after equation (104)). As is seen, when $r,k\to 0$ the conventional
first term in (48) dominates so the leading terms in the perturbative
expansion will coincide with the standard theory. The diagrams containing the
second term will differ from the leading contribution in that the propagator
will be replaced by a constant and so the diagram should be easier to
calculate than the main contribution. It may be hoped that when $t_{1}$ is
sufficiently close to $t^{\mathrm{R}}$ the corrections to the critical
temperature will be small and good precision of $K_{c}$ can be preserved. The
critical exponents in this approach will be the same as those calculated
within the chosen perturbative technique [1, 34, 35, 36].

##  6 First order phase transitions

For simplicity we restrict our discussion of FOPTs to the ferromagnetic IM
where transitions occur below the critical temperature when the external
magnetic field $h$ changes between infinitesimally small values $0^{-}$ and
$0^{+}$. During the change the spontaneous magnetization jumps from $-m_{0}$
to $m_{0}>0$ or vice versa so the susceptibility $dm/dh$ is infinite inside
the coexistence region but in the immediate vicinity of $h=0$ it should be
finite. This behaviour means that the free energy per site $f(h)$ contains a
singular contribution of the form $f_{sing}\simeq m_{0}|h|$ [45, 46] with the
respective term in the local potential

$u_{sing}^{\mathrm{R}}=-rm_{0}|x|.$ (49)

Though mathematically simple, this term is not straightforward to obtain as a
solution of RG equations because the local potential (34) is assumed to be an
analytic function of the field variable.

The origin of the singularity can be elucidated by differentiating LPA
equation (40) w.r.t. $x$ to arrive at the generalized BE [32]

$\mu_{t}+\mu\mu_{x}=(p/2)\mu_{xx},$ (50)

where

$\mu(x,t)=u_{x}(x,t).$ (51)

According to (49), $\mu^{\mathrm{R}}(x)$ should be discontinuous at $x=0$. As
is known [47], the discontinuous shock waves appear in the inviscid BE which
in our case would mean vanishing $p$. According to (39), $p=0$ at the end of
renormalization when $t^{\mathrm{R}}=1/r$. It is important to note that in
models studied in [17, 19] and in the present paper $p(t)$ additionally
satisfied the condition $(dp/dt)_{t^{\mathrm{R}}}=0$ due to $\rho(E=0)=0$ in
(39). It seems that the latter condition should also be satisfied for the
presence of FOPTs because the calculations with a model which did not satisfy
it the discontinuous solutions could not be found.

To describe FOPTs with the use of the equation of state (43) it would be
sufficient to solve (50) and (51) using the method of lines [48]. However,
because $\mu^{\mathrm{R}}(x)$ is discontinuous, its first and the second
derivatives entering (50) will be difficult to deal with numerically.
Therefore, in the actual calculations in [17, 19] the difficulties were
alleviated by using the integrated form of (50), that is, the original LPA
equation (40). Furthermore, additional regularization has been found to be
useful that can be achieved with the use of the partial LT (86)–(87). By
solving the transformed equation (88) instead of (40) it can be seen from (86)
that $x=h/r$ at $t^{\mathrm{R}}$ so the jump in $u_{x}$ at $x=0$ is
transformed in a finite interval of $y$ values. Thus, the only singularities
in $v^{\mathrm{R}}(y)$ will be two kinks at points
$y^{\pm}_{0}=-\bar{t}^{\mathrm{R}}u^{\mathrm{R}}_{x}|_{x=0^{\pm}}$ while
inside the interval $\mu^{\mathrm{R}}(y)=v^{\mathrm{R}}_{y}(y)$ will be a
linear function of $y$, as can be seen from (90) for $h=0$. So the strongest
singularities in equation (88) will be two discontinuities in the second
derivative at the kink points corresponding to the jump of the susceptibility
at the boundary separating the ordered phase and the coexistence region.

The mechanism of appearance of the linear in $y$ segment in $\mu(y,t)$ can be
qualitatively understood from the equation obtained by differentiating (88)
w.r.t. $y$:

$\mu_{t}=\frac{p(t)\mu_{yy}}{2(1+\bar{t}\mu_{y})^{2}}.$ (52)

As is seen, when the term in parentheses vanishes the r.h.s. becomes singular.
Because $\bar{t}\geq 0$ this may happen only for negative values of $\mu_{y}$
which appear when $v(y,t)$ has a negative curvature in some region of $y$
values. The latter appears in the initial potentials (109) $u$ or $v$ which
are equal due to (87) at sufficiently low temperatures. Further, the
coefficient before $\mu_{yy}$ in (52) is always positive and so can be
considered as a space and time dependent diffusivity. At vanishing denominator
the singularity becomes non-integrable so the maximum negative slope tolerated
by equation (52) is $-1/\bar{t}^{\mathrm{R}}$; steeper slopes will be smeared
out by the diffusion because the diffusivity will diverge in such cases. In
numerical calculations in [18, 17, 19] and in the exactly solvable model below
it was found that the fully renormalized solution exhibited the universal
behaviour $\mu^{\mathrm{R}}(y)=-y/\bar{t}^{\mathrm{R}}+(f.i.t.)$ corresponding
to the maximum allowable steepness. Substitution in (90) shows that the change
of $\mu^{\mathrm{R}}(x)$ within the coexistence region is confined to a single
point $x=0$, as expected.

The above qualitative reasoning can be illustrated with an exactly solvable
example provided by the infinite-range IM (IRIM) also known as the Husimi-
Temperley [49] and the Curie-Weiss model [50] which is often considered to be
exactly solvable in the MF approximation [50]. A straightforward application
of the MF theory, however, leads to unphysical van der Waals loops in the FOPT
region which suggests that the MF solution below $T_{c}$ is flawed. In [49,
51] it was shown that the loops are replaced by the shock wave solutions of
the BE that can be derived for that particular model without resort to RG.

In D it has been shown that exactly the same BE as in [49, 51] is obtained for
this model as the S-matrix RG equation in the LPA. In the IRIM case $p(t)$ in
(40) becomes equal to $1/N$ so that the generalized BE (50) turns into the
conventional BE with constant viscosity $1/2N$. To our purposes the most
important is the fact that the jumps in $\mu^{\mathrm{R}}(x)$ as the external
field crosses zero are between the MF spontaneous magnetisation values $\pm
m_{0}$. Thus, though $\mu=u_{x}$ in (50) and $m(x)$ in [49] do not coincide,
from (43) it can be seen that at $h=0$ the functions in both cases are the
same. Now taking into account that our LPA RG solution in D coincides with the
MF solution [50] outside the coexistence region and unifying it with the
rigorous treatment in [49] of discontinuities in the solution of BE which is
equally valid for our $h=0$ case we conclude that the LPA RG equation solves
exactly both the IRIM and the problem of the van der Waals loops. Finally, in
order to check whether the numerical solution by the method of lines used in
[18, 17, 19] gives in the IRIM case the same solution as the analytic
approach, (88) has been solved with $p=1/2N$ for $N=1000$ and the exact
solution has been reproduced with the accuracy $\Or(10^{-3})$.

In connection with the numerical solutions it is pertinent to note that though
(88) in the FOPT region has been found to be easier to deal with than with
equation (40) [17, 19], it may still be worthwhile to adopt the techniques of
[48] to the case of (40). The reason is that the phase transitions in
$u^{\mathrm{R}}(x)$ take place at a single value of argument $x$ while in (88)
$x$ is mapped on an interval of $y$ values of length $\Or(1)$. But
$v^{\mathrm{R}}(y)$ within the interval behaves linearly and in the symmetric
case can be fully characterized by a single parameter, the slope. However, for
good numerical accuracy the discretization step has to be chosen to be
$\Or(10^{-3})$ or smaller [22, 17]. Thus, in the method of lines $\Or(10^{3})$
points carry essentially the same physical information. Because in practice
the maximally possible number of lines is restricted, the use of (40) may
prove to be more suitable for achieving better accuracy.

##  7 Conclusion

In the extensive review of the nonperturbative RG approach [13] the authors
motivated the need for the modern implementation of Wilson’s RG suggested in
[7, 8, 14, 9] by the complexity of previously derived EFRGEs [1, 2, 3] which
impeded the development of reliable approximate computational schemes. The
results of the present paper can be considered as further advance in this
direction in the case of the classical lattice models.

The new advancements include such important achievement as the derivation of
the S-matrix EFRGE and of its Legendre transformed form which account for all
contributions to the free energy and EA functionals, not only for the field-
dependent terms as is conventional in the Wilsonian RG [1, 9]. This has made
possible the RG calculation of all thermodynamic quantities in lattice
systems. Moreover, by LT of the S-matrix EFRGE it has been possible to
identify the cause of the divergence inherent in the conventional formulation
of the Wetterich equation [7, 9, 11, 13]. An explicit expression for a
divergent conterterm has been derived which having being subtracted from the
fully renormalized solution leads to the exact finite expression for the Gibbs
FE. This allows one to obtain with the use of the Wetterich equation the same
results as with the S-matrix EFRGE.

Further, within the S-matrix RG approach a thermodynamic description of the
FOPT of general kind may be achieved. By the latter are meant the FOPT not
only due to the symmetry breaking as in the IM below $T_{c}$ but also those
that are defined by the crossing of the FE curves, such as in the Blume-Capel
model below the tricritical point [19]. The description is further facilitated
by the structure of the S-matrix EFRGE, especially in the LPA, which can be
straightforwardly connected with the generalized BE through its
differentiation w.r.t. the field variable. This has made possible to
generalize on the nonperturbative RG the observation made in connection with
the MF solution of IRIM in [49, 51] that FOPT correspond to the shock-wave
solutions of RG equations. In the case of EA EFRGEs this could be
approximately justified only in the large-$n$ limit [33] which excludes the
majority of physically important small-$n$ models, such as the IM, the Blume-
Capel model, etc.. The connection of S-matrix EFRGE with the generalized
Burgers’ equation has made possible to apply in the solution of RG equations
the concepts and techniques developed in the field of nonlinear differential
equations [47, 32, 39].

Some useful improvements in the S-matrix approach to the lattice models in
comparison with the approach of [15, 16] can be mentioned. For example, the
problem of establishing the initial condition does not arise in the S-matrix
case; the need to calculate improper integrals is absent because the momenta
much larger than the inverse lattice spacing do not appear in the formalism.
All necessary LTs have been found in explicit analytic form and the need for
their numerical calculation do not arise. These improvements are important in
the solution of the evolution equations of RG type because the positive
Lyapunov exponents inherent in RG in the critical region may enhance small
errors in the initial condition resulting in degraded accuracy.

The absence of large unphysical momenta in the S-matrix EFRGE is also
important from the point of view of qualitative soundness of the approximate
solutions obtained. In the absence of rigorous methods of solution of strongly
coupled field-theoretic models, qualitative arguments play a significant role.
Therefore, the use in the S-matrix approach of the SC condition and of the
amputated correlation functions in the renormalized interaction potential
makes LPA a natural approximation from a qualitative standpoint. Furthermore,
the use of the cutoff in the form of the step function [5, 15] in the spin-
lattice case leads to the exact renormalization at the earliest stage of the
RG flow when the Fourier momenta are not yet eliminated. In the present paper
it has been argued that at an early stage of the momenta elimination the
gradual loss of LPA accuracy caused by the widening gaps in the periodic zone
scheme is accompanied by the weakening of interactions. So that at the late
stage of the RG flow they may become sufficiently weak for LPA to be
justifiable by DE [24]. This qualitative reasoning agrees with the results of
[16] where the use of the step function cutoff in the strong coupling case
resulted in considerably more accurate values of critical temperatures than
with the Wilson-type cutoff. The values of the critical temperatures
calculated in [16] within the EA LPA in $\varphi^{4}$ model are very similar
to those of the present paper based on the S-matrix equation. Apparently this
is a consequence of the equivalence of the EA and the S-matrix equations also
in the LPA. The high accuracy of the obtained results further supports the
qualitative correctness of the picture of the RG flow described by the
S-matrix EFRGE.

Further, the possibility to easily implement a multi-step renormalization
within the S-matrix formalism makes possible to use in the critical region the
accurate perturbative results obtained, for example, with the use of DE [35,
36]. Such an approach may be useful in correcting the deviations of the LPA
predictions from the experimental data caused by errors in the LPA critical
exponents [17].

To sum up, the main advantage of the RG approach based on the S-matrix EFRGE
is that it provides conceptually simpler and computationally easier techniques
of solution of spin-lattice models than the alternative methods [15, 23] while
preserving the same accuracy in the calculation of the critical temperatures
as was achieved in [15, 16]. Besides, the possibility to calculate correct
values of the f.i.t. in FE within the S-matrix formalism allows for its use in
the description of FOPT of any kind [19].

This research was supported by the National Academy of Sciences of Ukraine
under contract No. 22/20-H. I am particularly indebted to Hugues Dreyssé for
his support and encouragement. I expresses my gratitude to Université de
Strasbourg and IPCMS for their hospitality.

##  Appendix A Stepwise renormalization in the Green’s function formalism

Renormalization in several finite steps can be performed recursively by adding
steps one at a time so only a two-step case need be considered. This can be
done with the use of the identity satisfied by the Gaussian kernel

$\displaystyle\int_{-\infty}^{\infty}d\vec{s}_{0}\;{\det\left(\frac{\hat{X}}{2\pi}\right)}^{\frac{1}{2}}\rme^{-\frac{1}{2}(\vec{s}-\vec{s}_{0})^{\dagger}\hat{X}(\vec{s}-\vec{s}_{0})}\;{\det\left(\frac{\hat{Y}}{2\pi}\right)}^{\frac{1}{2}}\rme^{-\frac{1}{2}(\vec{s}_{0}-\vec{s}\,^{\prime})^{\dagger}\hat{Y}(\vec{s}_{0}-\vec{s}\,^{\prime})}$
$\displaystyle={\det\left(\frac{\hat{X}\hat{Y}}{2\pi(\hat{X}+\hat{Y})}\right)}^{\frac{1}{2}}\rme^{-\frac{1}{2}(\vec{s}-\vec{s}\,^{\prime})^{\dagger}\frac{\hat{X}\hat{Y}}{(\hat{X}+\hat{Y})}(\vec{s}-\vec{s}\,^{\prime})}$
(53)

where all matrices are considered to be translationally-invariant, hence,
commuting with each other.

Now by assuming that the second kernel on the first line is ${\cal
G}(\vec{s}_{0},\vec{s}\,^{\prime},t_{0})$ from (10) with
$\hat{Y}=\hat{R}(t_{0})$ and that on the second line we want to obtain (10),
we conclude that

$\frac{\hat{X}\hat{R}(t_{0})}{\hat{X}+\hat{R}(t_{0})}=\hat{R}(t).$ (54)

Solving this one finds that the regulator matrix for the second step is

$\hat{X}(t,t_{0})=\left[\hat{R}^{-1}(t)-\hat{R}^{-1}(t_{0})\right]^{-1}=\left[\hat{G}(t)-\hat{G}(t_{0})\right]^{-1}$
(55)

As is seen, $\hat{X}$ satisfies the initial condition for Green’s functions
$\hat{X}(t\to t_{0},t_{0})\to\infty$; besides, because $\hat{R}(0)=\infty$,
$\hat{X}(t,0)=\hat{R}(t)$ so (55) is valid for all $t$ and $t_{0}\leq t$.

##  Appendix B Legendre transforms

By analogy with the LTs of [52, 18] we introduce a new fluctuating field
$\vec{\phi}$ and a new interaction potential $W$ as

$\displaystyle\phi_{\bf k}({\vec{s}},t)$ $\displaystyle=$ $\displaystyle
s_{\bf k}-G({\bf k},t)\frac{\partial U({\vec{s}},t)}{\partial{s_{\bf-k}}}$
(56) $\displaystyle W(\vec{\phi},t)$ $\displaystyle=$ $\displaystyle
U({\vec{s}},t)-\frac{1}{2}\sum_{\bf q}\frac{\partial
U({\vec{s}},t)}{\partial{s_{\bf-q}}}G({\bf q},t)\frac{\partial
U({\vec{s}},t)}{\partial{s_{\bf q}}}$ (57)

where we assume $G({\bf k},t)$ to be a non-negative even function of ${\bf k}$
satisfying $G({\bf k},0)=0$.

To transform the exact RG equation (29) to the new variables we first take the
derivatives of (56) and (57) w.r.t. $s_{{\bf k}^{\prime}}$

$\displaystyle\frac{\partial\phi_{\bf k}}{\partial s_{{\bf
k}^{\prime}}}=\delta_{{\bf k}{\bf k}^{\prime}}-G({\bf
k},t)\frac{\partial^{2}U}{\partial{s_{\bf-k}}\partial s_{{\bf k}^{\prime}}}$
(58) $\displaystyle\sum_{\bf k}\frac{\partial W}{\partial{\phi_{\bf
k}}}\frac{\partial\phi_{\bf k}}{\partial s_{{\bf k}^{\prime}}}=\frac{\partial
U}{\partial{s_{{\bf k}^{\prime}}}}-\sum_{\bf q}\frac{\partial
U({\vec{s}},t)}{\partial{s_{\bf q}}}G({\bf
q},t)\frac{\partial^{2}U}{\partial{s_{\bf-q}}\partial s_{{\bf k}^{\prime}}}$
(59)

and observe that by substituting (58) in (59) and changing the summation
subscript from ${\bf k}$ to ${\bf q}$ one arrives at the system of linear
equations for the row vector $\left({\partial W}/{\partial{\phi_{\bf
q}}}-{\partial U}/{\partial{s_{\bf q}}}\right)^{\dagger}$

$\sum_{\bf q}\left(\frac{\partial W}{\partial{\phi_{\bf q}}}-\frac{\partial
U}{\partial{s_{\bf q}}}\right)^{\dagger}\left(\delta_{{\bf q}{\bf
k}^{\prime}}-G({\bf q},t)\frac{\partial^{2}U}{\partial{s_{\bf-q}}\partial
s_{{\bf k}^{\prime}}}\right)=0.$ (60)

Because $G({\bf q},t)$ can be arbitrary, the matrix of the system in general
is not degenerate which means that the solution should be trivial and the
following useful equality should hold

$\frac{\partial W}{\partial{\phi_{\bf q}}}=\frac{\partial U}{\partial{s_{\bf
q}}}.$ (61)

In particular, with its use transformations (56) and (57) are easily reversed.
For example, from (56) one gets

$s_{\bf k}(\vec{\phi},t)=\phi_{\bf k}+G({\bf k},t)\frac{\partial
W}{\partial{\phi_{\bf-k}}}(\vec{\phi},t)$ (62)

The next step in transforming the RG equation is to differentiate (56) and
(57) w.r.t. $t$:

$\displaystyle\partial_{t}\phi_{{\bf k}}=-\partial_{t}G({\bf
k},t)\frac{\partial U}{\partial{s_{\bf-k}}}-G({\bf
k},t)\frac{\partial^{2}U}{\partial{s_{\bf-k}}\partial t}$ (63)
$\displaystyle\partial_{t}W+\sum_{\bf k}\frac{\partial W}{\partial{\phi_{\bf
k}}}\partial_{t}\phi_{{\bf k}}=\partial_{t}U-\sum_{\bf q}\frac{\partial
U}{\partial{s_{\bf q}}}G({\bf
q},t)\frac{\partial^{2}U}{\partial{s_{\bf-q}}\partial t}$
$\displaystyle-\frac{1}{2}\sum_{\bf q}\frac{\partial
U}{\partial{s_{\bf-q}}}\partial_{t}G({\bf q},t)\frac{\partial
U}{\partial{s_{\bf q}}}$ (64)

By substituting (63) in (B) and using (61) one can express the terms on the
l.h.s. in (29) in terms of $W(\vec{\phi},t)$ as

$\partial_{t}U+\frac{1}{2}\sum_{\bf q}\frac{\partial
U}{\partial{s_{\bf-q}}}\partial_{t}G({\bf q},t)\frac{\partial
U}{\partial{s_{\bf q}}}=\partial_{t}W.$ (65)

The r.h.s. in (29) can be transformed by differentiating (61) and (62) w.r.t.
$\phi_{{\bf-k}^{\prime}}$ as

$\displaystyle\sum_{\bf q}\frac{\partial^{2}U}{\partial s_{{\bf
k}}\partial{s_{\bf-q}}}\left[\delta_{{\bf q},{\bf k}^{\prime}}+G({\bf
q},t)\frac{\partial^{2}W}{\partial{\phi_{\bf
q}\partial\phi_{{\bf-k}^{\prime}}}}\right]=\frac{\partial^{2}W}{\partial{\phi_{\bf
k}\partial\phi_{{\bf-k}^{\prime}}}}$ (66)

Now denoting the matrix in the brackets by $\hat{B}=[B_{\bf q,k^{\prime}}]$
(66) can be solved as

$\frac{\partial^{2}U}{\partial s_{{\bf k}}\partial{s_{\bf-k}}}=\sum_{\bf
k^{\prime}}\frac{\partial^{2}W}{\partial{\phi_{\bf
k}\partial\phi_{{\bf-k}^{\prime}}}}\left(\hat{B}^{-1}\right)_{\bf
k^{\prime},k}.$ (67)

This expression can be substituted in (29) to give the exact RG equation in
terms of the Legendre-transformed quantities $\vec{\phi}$ and $W$

$\partial_{t}W=\sum_{\bf k^{\prime},k}\partial_{t}G({\bf
k},t)\frac{\partial^{2}W}{\partial{\phi_{\bf
k}\partial\phi_{{\bf-k}^{\prime}}}}\left(\hat{B}^{-1}\right)_{\bf
k,k^{\prime}}$ (68)

which structure is very similar (but not identical) to that of equation (3.31)
in [9].

###  B.1 The Wetterich equation

The EA EFRGE derived in [7, 11] (the Wetterich equation) reads

$\partial_{t}\tilde{\Gamma}(\vec{\varphi},t)=\frac{1}{2}\Tr\left\\{\partial_{t}\hat{{R}}^{(0)}\left(\hat{\tilde{\Gamma}}_{\varphi\varphi}+\hat{{R}}^{(0)}\right)^{-1}\right\\}$
(69)

where $\hat{{R}}^{(0)}$ in addition to the requirement to all regulators (11)
also satisfies the condition

$\hat{R}^{(0)}(t=t^{\mathrm{R}})=0$ (70)

(the superscript (0) is used to distinguish such regulators). The RG flow in
(69) is initiated by

$\tilde{\Gamma}(\vec{\varphi},t=0)=H_{0}(\vec{\varphi}).$ (71)

To establish its connection with (68) we first derive an expression for the
Gibbs FE that is equal to the fully renormalized EA. To this end using Fourier
transformed (7) and (21) one finds

$m_{\bf k}=G^{\mathrm{P}}({\bf k})h_{\bf k}-G^{\mathrm{P}}({\bf
k})\left.\frac{\partial
U^{\mathrm{R}}({\vec{s}})}{\partial{s_{\bf-k}}}\right|_{s_{\bf
k}=G^{\mathrm{P}}({\bf k})h_{\bf k}}.$ (72)

Solving this w.r.t. the derivative

$\left.\frac{\partial
U^{\mathrm{R}}({\vec{s}})}{\partial{s_{\bf-k}}}\right|_{s_{\bf
k}=G^{\mathrm{P}}({\bf k})h_{\bf k}}=h_{\bf k}-\left[G^{\mathrm{P}}({\bf
k})\right]^{-1}m_{\bf k},$ (73)

substituting it in (57) at $t^{\mathrm{R}}$ where $\phi_{\bf k}=m_{\bf k}$ and
using (20) one gets in the vector notation

$W^{\mathrm{R}}(\vec{m})+\frac{1}{2}\vec{m}\,^{\dagger}(\hat{\epsilon}+r\hat{I})\vec{m}=U^{\mathrm{R}}(\hat{G}^{\mathrm{P}}\vec{h})-\frac{1}{2}\vec{h}^{\dagger}\hat{{G}}^{\mathrm{P}}\vec{h}+\vec{m}\,^{\dagger}\vec{h}.$
(74)

Comparing this with (21) and adding the last term in (21) to both sides

$W^{\mathrm{R}}(\vec{m})+\frac{1}{2}\vec{m}\,^{\dagger}(\hat{\epsilon}+r\hat{I})\vec{m}-\frac{1}{2}\Tr\ln{\left(2\pi\hat{G}^{\mathrm{P}}\right)}=F^{\mathrm{R}}(\vec{h})+\vec{m}\,^{\dagger}\vec{h}$
(75)

one sees that the r.h.s. is the $\vec{h}-\vec{m}$ Legendre transform of the
Helmholtz FE, so, by definition, it is equal to the Gibbs FE which we will
denote $\Gamma^{\mathrm{R}}(\vec{m})$. The expression is naturally generalized
at arbitrary $t$ as

$\Gamma(\vec{\phi},t)=W(\vec{\phi},t)+\frac{1}{2}\vec{\phi}\,^{\dagger}(\hat{\epsilon}+r\hat{I})\vec{\phi}-\frac{1}{2}\Tr\ln{\left(2\pi\hat{G}^{\mathrm{P}}\right)}.$
(76)

Further, because at $t=0$ the LT (56)–(57) is trivial, the initial condition
for (76) is obtained by simply changing the variable from $\vec{s}$ to
$\vec{\phi}$ so taking into account (2) one gets

$\Gamma(\vec{\phi},t=0)=H_{0}(\vec{\phi})-\frac{1}{2}\Tr\ln{\left(2\pi\hat{G}^{\mathrm{P}}\right)}$
(77)

which differs from the condition (71) of EA approach by the second term.
Because the latter is constant in $\vec{\phi}$ and $t$, it remains unaltered
in the course of the RG flow. Therefore, the constant term instead of the
initial condition (77) can be added to the solution at the end of the flow
which below will be used in the renormalization of
$\tilde{\Gamma}^{\mathrm{R}}$ in the Wetterich equation.

The EFRGE for $\Gamma$ can be obtained as follows. First, in (68) we represent
matrix $\hat{B}$ as the product $\hat{G}\hat{A}$ with $\hat{G}=[G({\bf
k},t)\delta_{\bf k,k^{\prime}}]$ and

$\hat{A}=\left[W_{\phi_{\bf k}\phi_{\bf-k^{\prime}}}+{R}({\bf k},t)\delta_{\bf
kk^{\prime}}\right]\equiv\hat{W}_{\phi\phi}+\hat{{R}}.$ (78)

We remind that according to (27) $\hat{R}=\hat{G}^{-1}$, hence, in (68)
$\partial_{t}\hat{G}=-\hat{R}^{-1}\partial_{t}\hat{R}\hat{R}^{-1}$. Next, in
the matrix identity

$\hat{W}_{\phi\phi}\hat{B}^{-1}=(\hat{A}-\hat{R})\hat{A}^{-1}\hat{R}=\hat{R}-\hat{R}\hat{A}^{-1}\hat{R}.$
(79)

the regulator function can be represented in the form

$\hat{R}=\hat{\epsilon}+r\hat{I}+\hat{{R}}^{(0)}$ (80)

compatible with both (16) and (70). Now substituting (79) in (68), using
$t$-independence of the last two terms in (76) and the invariance of the trace
under the cyclic permutations we arrive at an EFRGE

$\partial_{t}\Gamma(\vec{\phi},t)=\frac{1}{2}\Tr\left\\{\partial_{t}\hat{R}^{(0)}\left(\hat{\Gamma}_{\phi\phi}+\hat{R}^{(0)}\right)^{-1}\right\\}-\frac{1}{2}\partial_{t}\ln\det\left(\hat{\epsilon}+r\hat{I}+\hat{{R}}^{(0)}\right)$
(81)

where $\hat{\Gamma}_{\phi\phi}=\left[\Gamma_{\phi_{\bf
k}\phi_{\bf-k^{\prime}}}\right]$.

As can be seen, (81) differs from the Wetterich equation (69) only by the last
term which, however, plays an significant role by regularizing the equation at
$t\to 0$ when $\hat{{R}}^{(0)}\to\infty$. Because in this limit $\Gamma$ in
(77) is finite and so is negligible in comparison with the regulator, the two
terms on the r.h.s. in (81) cancel to zero at $t=0$ so unlike (69) the
equation is integrable near $t=0$ as well as everywhere along the RG
trajectory. This ensures the use of (81) in determining the exact Gibbs FE.

Because the cancellation of singularities at $t=0$, is very important for the
consistency of (81), it cannot be rearranged in the WE form through a simple
transfer of the last term on the r.h.s. on the l.h.s.. To overcome this
difficulty we first regularize the equation in the conventional way by
starting the evolution from small positive value of $t=t_{0}$ (which is
equivalent to starting from large finite value of $\Lambda=\Lambda_{0}$ or $k$
in the momentum parametrization [7, 9, 11, 13]). Now it is easily seen, that
(81) can be cast in the form of WE (69) for a modified EA

$\tilde{\Gamma}(\vec{\phi},t,t_{0})={\Gamma}(\vec{\phi},t)+\frac{1}{2}\Tr\ln{\left(2\pi\hat{G}^{\mathrm{P}}\right)}+\frac{1}{2}\ln\det\left.\left(\hat{\epsilon}+r\hat{I}+\hat{{R}}^{(0)}(t^{\prime})\right)\right|_{t^{\prime}=t_{0}}^{t^{\prime}=t}$
(82)

where the second term on the r.h.s. has been added in order that instead of
(77) $\tilde{\Gamma}$ satisfied at $t_{0}\to 0$ the initial condition (71).

At the end of the evolution at $t=t^{\mathrm{R}}$ after some rearrangement the
fully renormalized Gibbs FE can be found with the use of the WE solution (82)
as

$\Gamma^{\mathrm{R}}(\vec{m})=\lim_{t_{0}\to
0}\left\\{\tilde{\Gamma}(\vec{m},t^{\mathrm{R}},t_{0})-\frac{1}{2}\ln\det{[2\pi\hat{G}(t_{0})]}\right\\}$
(83)

where the second term in the braces is the counterterm with $\hat{G}$ equal to
the inverse of the regulator (80).

It is to be pointed out that when differentiated w.r.t. $\phi$ the equations
(69) and (81) and their respective initial conditions (71) and (77) becomes
identical which proves that both the equations and their solutions differ only
by f.i.t.. It is this differentiated form of WE that has been used by
Wetterich et al. in concrete RG calculations [53, 26].

Thus, the Wetterich equation (69) solved with the conventional initial
condition (71) will give the exact solution of the EA EFRGE both for the
field-dependent part and for f.i.t. provided the divergence in the EA will be
removed according to (83).

###  B.2 LT in a two-step renormalization

As could be noted, the only property of $G({\bf k},t)$ that was used in the
derivation of RG equation (68) was that $G$ had the same time derivative as in
(29). But from (32) one can see that $\partial_{t}G=\partial_{t}\Delta$. So
for the second-step renormalization one can derive RG equation similar to (68)
by using the LT (56) and (57) with $G$ replaced by $\Delta$

$\displaystyle y_{\bf k}({\vec{s}},t)$ $\displaystyle=$ $\displaystyle s_{\bf
k}-\Delta({\bf k},t,t_{0})\frac{\partial U({\vec{s}},t)}{\partial{s_{\bf-k}}}$
(84) $\displaystyle V(\vec{y},t)$ $\displaystyle=$ $\displaystyle
U({\vec{s}},t)-\frac{1}{2}\sum_{\bf q}\frac{\partial
U({\vec{s}},t)}{\partial{s_{\bf-q}}}\Delta({\bf q},t,t_{0})\frac{\partial
U({\vec{s}},t)}{\partial{s_{\bf q}}}.$ (85)

However, unless $t_{0}=0$, functional $V$ and variables $\vec{y}$ do not
acquire immediate physical meaning at the end of renormalization (for
$t_{0}=0$ they will coincide with $W$ and $\vec{\phi}$). To see this we note
that the equation for $V$ similar to (68) has as the initial condition
$V(\vec{y},t_{0})=U(\vec{y},t_{0})$, as follows from (84)–(85) with
$\Delta(t_{0},t_{0})=0$ (see (32)). But according to (29) $U(\vec{y},t_{0})$
depends on the arbitrary values of $G({\bf k},t)$ in the interval $0\leq
t<t_{0}$ while in the RG equation for $V$ and in the LT (84)–(85) only $G({\bf
k},t)$ for $t_{0}\leq t<t^{\mathrm{R}}$ contribute. Because of this
arbitrariness, both $V^{\mathrm{R}}$ and $\vec{y}$ cannot be given a physical
meaning at $t^{\mathrm{R}}$. However, similar to (61) the equality $\partial
V/\partial{y_{\bf q}}={\partial U}/{\partial{s_{\bf q}}}$ holds so the
transform (84)–(85) can be reversed and the physical quantities
$U^{\mathrm{R}}$ and $\vec{s}$ in (18) can be expressed in parametric form in
terms of $V^{\mathrm{R}}$ and $\vec{y}$. If needed, the Gibbs FE can be found
with the use of (75).

##  Appendix C Legendre transforms in LPA

The formalism just described considerably simplifies in the LPA. First we note
that in the layer-cake renormalization scheme of section 4 visualized in
figure 1 the fluctuating field “sees” only the flat part of $G({\bf k},t)=t$
(38) because beyond the cutoff the field has been set to zero. Thus,
$\Delta^{LPA}(t,t_{0})=t-t_{0}=\bar{t}$ (see figure 1). In terms of the local
potential $v(y,t)$ corresponding to $V$ which now depends on a single scalar
variable $y$ the LT (84)–(85) simplifies to

$\displaystyle y$ $\displaystyle=$ $\displaystyle x-\bar{t}u_{x}$ (86)
$\displaystyle v$ $\displaystyle=$ $\displaystyle u-\bar{t}u_{x}^{2}/2$ (87)

where for brevity the arguments of $y(x,t)$, $v(y,t)$ and $u(x,t)$ have been
omitted and the subscript notation for partial derivative has been used.

The exact equation for $V$ is obtained from (68) via substitutions $W\to V$,
$\vec{\phi}$ to $y$ and $G\to\Delta$. In the LPA it simplifies to

$v_{t}=\frac{p(t)v_{yy}}{2(1+\bar{t}v_{yy})}.$ (88)

It is possible also to derive it directly from (40) by using (86)–(87) and
repeating the steps from B [18].

Similar to (61) in the exact case, the equality

$v_{y}(y,t)=u_{x}(x,t)$ (89)

holds also in LPA so at the end of renormalization this can be used to find
the equation of state in parametric form. Thus, from (86) one gets at
${t}^{\mathrm{R}}$ where according to (42)

$h/r=y+\bar{t}^{\mathrm{R}}v^{\mathrm{R}}_{y}.$ (90)

The expression for the magnetisation can be expressed through $y$ by replacing
in (43) $h/r$ by the r.h.s. of (90) which gives

$m(y)=y-t_{0}v^{\mathrm{R}}_{y}.$ (91)

As is seen, $y=m$ only when $t_{0}=0$ in which case the equation of state is
obtained by replacing $y$ in (90) by $m$. When $t_{0}\not=0$ (90) and (91)
define the equation of state parametrically.

##  Appendix D LPA solution of IRIM

In the IRIM all spin pairs interact with the same dimensionless strength $K/N$
so the sites can be numbered in an arbitrary order which means that the model
is structureless. Because the use of the conventional lattice Fourier
transform is not warranted in this case, all calculations will be carried out
in the space of the site numbers.

As is easy to see, the pair interaction matrix in IRIM can be cast in the form

$\hat{\epsilon}=K(\hat{I}-\hat{E})$ (92)

where the idempotent matrix $\hat{E}$, $\hat{E}^{2}=\hat{E}$, has matrix
elements ${E}_{ij}=N^{-1}$ and, as a consequence, matrix $\hat{I}-\hat{E}$ is
also idempotent. The spectrum of idempotent matrices consists of only two
points: 0 and 1, so the spectrum of $\hat{\epsilon}$ in (92) is $(0,K)$ so
similar to the lattice case the lowest eigenvalue of $\hat{\epsilon}$ is zero.

In the absence of the momentum representation the density of states can be
found according to the formula
$\rho(E)=-\pi^{-1}\Im\Tr(E+\mathrm{i}\varepsilon-\hat{\epsilon})^{-1}$. Now
using the idempotence of $\hat{\epsilon}$ it is easy to find that

$\left(\frac{1}{z-\hat{\epsilon}}\right)_{ii}=\frac{1}{N}\frac{1}{z}-\left(1-\frac{1}{N}\right)\frac{1}{z-K}.$
(93)

With $z=E+\mathrm{i}\varepsilon$ one gets

$\rho(E)=N^{-1}\delta(E)+(1-N^{-1})\delta(E-K)$ (94)

so that according to (39)

$p(t)=N^{-1}+(1-N^{-1})\theta({t^{-1}-r}-K).$ (95)

As is seen, when

$0\leq t\leq t_{0}=1/(r+K)$ (96)

$p(t)=1$ and in this range the RG equation can be solved by (36) with
$u(x,t_{0})$ for IM (which includes IRIM as a special case) calculated
explicitly in (110). By reminding that the LT (86)–(87) at $t_{0}$ is trivial,
on the basis of (110) we obtain the explicit expression

$v(y,t_{0})=\frac{y^{2}}{2t_{0}}-\ln\left(\cosh\frac{y}{t_{0}}\right)+\frac{1}{2}\ln(2\pi
t_{0}).$ (97)

For $t>t_{0}$ $p(t)=1/N$ in (95) so the RG equation (50) for IRIM becomes the
conventional BE with the constant viscosity $1/2N$ and in the thermodynamic
limit $N\to\infty$ it reduces to the inviscid BE with the discontinuous shock
wave solutions corresponding to FOPTs [49, 51]. The transformed equation (88)
acquires a simple form in this limit

$v_{t}|_{N\to\infty}=0.$ (98)

In this form, however, it produces the MF solution exhibiting the undesirable
van der Waals loops in the coexistence region. To see this we note that (98)
implies that the initial condition does not change with $t$ so at
$t^{\mathrm{R}}=1/r$ the renormalized $v^{R}(y)=v(y,t_{0})$. Substituting it
in (90) and (91) on gets after some rearrangement the MF equation of state in
parametric form

$\displaystyle h$ $\displaystyle=$ $\displaystyle y/t_{0}-K\tanh(y/t_{0})$
(99) $\displaystyle m$ $\displaystyle=$ $\displaystyle\tanh(y/t_{0})$ (100)

from which immediately follows the exact IRIM MF equation [50]

$m=\tanh(Km+h).$ (101)

The MF free energy of IRIM [50]

$f=Km^{2}/2-\ln\left[2\cosh(Km+h)\right].$ (102)

is obtained from (41) with the use of inverted transform (87), (97), (100) and
noticing that the last term in (97) is cancelled by the last term in (41) as
is easily calculated with the use of definitions of $\hat{G}^{\mathrm{P}}$
(20) and (92). As is seen, the arbitrary parameter $r$ completely disappears
from the solution given by equations (101) and (102). The MF solution,
however, is not fully consistent because FE (102) is not convex and as a
consequence the van der Waals loops will appear in the equation of state. As
explained in the main text, these deficiencies can be avoided if the
thermodynamic limit in (98) is taken more carefully.

##  Appendix E Particularity of the lattice LPA

The LPA ansatz (33) in general case does not satisfy the S-matrix EFRGE (29)
because of the second term on the r.h.s.. However, the continuum and the
lattice cases differ from each other due to the difference between the
Kronecker symbols in (33) which on the lattice takes the form

$\delta^{L}_{{{\bf k}_{1}}+{{\bf k}_{2}}+\dots+{{\bf k}_{l}}}=\delta_{{{\bf
k}_{1}}+{{\bf k}_{2}}+\dots+{{\bf k}_{l}},0}+\sum_{{\bf K}\not=0}\delta_{{{\bf
k}_{1}}+{{\bf k}_{2}}+\dots+{{\bf k}_{l},{\bf K}}}$ (103)

where on the r.h.s. the deltas are the conventional Kronecker symbols and in
the summation over the reciprocal lattice vectors ${\bf K}$ we singled out
term ${\bf K=0}$ which is the only one that is present in the continuum theory
[1, 5].

Substituting (33) in the nonlinear term in (29) and using the first Kronecker
delta to lift the summation over ${\bf k}$ and the definition (103) one
arrives at the sum of quadratic in $u_{l}$ contributions of the form

$\displaystyle u_{l}(t)u_{l^{\prime}}(t)\sum_{\bf k}\theta\left[G({\bf
k},t)-t\right]\delta^{L}_{\sum_{j=1}^{l-1}\mathbf{k}_{j}-\mathbf{k}}\delta^{L}_{\mathbf{k}+\sum_{j^{\prime}=1}^{l^{\prime}-1}\mathbf{k^{\prime}}_{j^{\prime}}}$
(104) $\displaystyle=$
$\displaystyle\delta^{L}_{\sum_{j=1}^{l-1}\mathbf{k}_{j}+\sum_{j^{\prime}=1}^{l^{\prime}-1}\mathbf{k^{\prime}}_{j^{\prime}}}u_{l}(t)u_{l^{\prime}}(t)\left.\sum_{\bf
K}\theta\left[\Lambda({\bf k}/|{\bf k}|,t)-|{\bf k}|\right]\right|_{{\bf
k}=\sum_{j=1}^{l-1}\mathbf{k}_{j}-\mathbf{K}}$

where $\Lambda$ in the argument of the step function on the second line should
be determined from equation (38) defining the momentum cutoff. In the lattice
case it depends on the crystal anisotropy, i.e., on the momentum direction
${\bf k}|{\bf k}|$. In the isotropic case $\Lambda$ would coincide with the
conventional momentum cutoff [1, 5].

If the step functions in (104) were equal to unity (the LPA) than the terms
(104) would sum up to the second term on the l.h.s. in (40). But this would be
valid only if the arguments of the step functions were always positive. In the
isotropic continuum theory when only ${\bf K=0}$ contributes to (104) this
would mean that the sum of ${\bf k}_{j}$ is smaller in absolute value than
$\Lambda$. But this can be guaranteed only for $l\leq 2$ because, by our
definition of the renormalization procedure, all individual field momenta
reside within $\Omega(t)$ defined in (38). However, for $l>2$ the absolute
value of the sum may exceed $\Lambda$ in which case the step functions will
dependent on ${\bf k}_{j}$. As a result, many contribution in the momenta
space will be lost from contributions $l,l^{\prime}>2$ and the errors
introduced by the LPA will be enhanced in the strong coupling case when
$u_{l>2}$ are large.

The problem alleviates in the lattice models when all ${\bf K}$ contribute in
(104). This can be visualized by considering the step functions in the
periodic zone scheme where each ${\bf K}$ becomes the centre of a region
$\Omega_{\bf K}$ which is the same $\Omega$ as in (38) only shifted on vector
${\bf K}$. Now the step function in (104) is equal to unity within all
$\Omega_{\bf K}$ and zero at the outside. In the course of the renormalization
this will create a muffin-tin structure where at the early stages of RG flow
there will much more volume in momentum space where the LPA conjecture is
fulfilled. For example, in the region $0\leq t\leq t_{0}$ in figure 1
$\Omega=$ BZ so in the periodic zone scheme the whole momentum space becomes
covered and the step function is unity everywhere. This makes LPA exact which
has allowed us to perform the exact renormalization in (36) and (105). As $t$
grows above $t_{0}$ all $\Omega_{\bf K}$ start to shrink but close to $t_{0}$
only narrow gaps between $\Omega_{\bf K}$ will appear so the violation of the
LPA assumption $\theta=1$ in (104) will be relatively small and can be roughly
assessed by the value of $\kappa(t)={\rm vol}[\bar{\Omega}(t)]/{\rm
vol}[{\Omega}(t)]$. We note that the estimate is not based on the smallness of
the momentum, as in DE, or on the values of $u_{l}$. In fact, it relies on the
closeness of $\Omega\simeq$ to BZ while $|{\bf k}|$ may be as large as
$\sim\sqrt{3}\pi$ on sc lattice. As $t$ advances farther toward
$t^{\mathrm{R}}$ $\Omega$ shrinks and $\kappa(t)$ grows to large values which
signals the breakdown of the LPA assumption (35) in a large part of BZ.
However, because under Wilsonian renormalization $u_{l}$ with large $l$ fast
attenuate towards small values [1], it is expected that the terms quadratic in
$u_{l}$ violating the LPA assumption will become small. Besides, closer to
$t^{\mathrm{R}}$ $\Omega(t)$ tends to be almost spherical, the lattice nature
of the system smears out and $k/\Lambda$ becomes an acceptable DE parameter to
justify LPA on the late stages of renormalization. It is to be stressed that
the above reasoning relied on the shape of the cutoff function similar to the
step function and is not applicable to the Wilson-type momenta elimination in
large-$k$ shells [1]. This seems to be confirmed by the RG calculations in
[16] where the use of the Wilson-type cutoff led to poor agreement with the MC
simulations at large interaction strengths while the use of the step function
resulted in a perfect agreement with the MC data.

##  Appendix F Initial condition in the spin models

In $O(n)$-symmetric case the $n$-dimensional integral in the exact initial
renormalization (36) takes the form

$\rme^{-u({\bf x},t_{0})}=\frac{1}{(2\pi t_{0})^{n/2}}\int d{\bf
x}_{0}\,\exp{\left(-\frac{({\bf x-x}_{0})^{2}}{2t_{0}}\right)}\rme^{-u({\bf
x}_{0},0)}.$ (105)

For spin-lattice models we assume that for all $n$ the spin length is equal to
unity

$\rme^{-u({\bf x}_{0},0)}=\delta({\bf x}_{0}^{2}-1).$ (106)

The integration in (105) is conveniently carried out in hyperspherical
coordinates [54] in which the integration over $|{\bf x}_{0}|$ is trivial due
to the delta-function in (106) and with the choice of the direction of ${\bf
x}$ along the first axis: ${\bf x}=(x\cos\vartheta,0,0,\dots,0)$, $x=|{\bf
x}|$ (105) reduces to

$\displaystyle\rme^{-u^{(n)}({\bf
x},t_{0})}=C_{n}\rme^{-\frac{t_{0}a^{2}}{2}}\int_{0}^{\pi}\rme^{\frac{x}{t_{0}}\cos\vartheta}\sin^{n-2}\vartheta\,d\vartheta$
$\displaystyle=C_{n}\rme^{-\frac{t_{0}a^{2}}{2}}\int_{-1}^{1}\rme^{az}(1-z^{2})^{\frac{n-3}{2}}\;dz\equiv
C_{n}\rme^{-\frac{t_{0}a^{2}}{2}}b(n,a)$ (107)

where in $C_{n}$ are gathered all ${\bf x}$-independent factors from (105) and
from the spherical volume [54] which can be easily recovered if needed;
$a=x/t_{0}$, the integral on the second line is obtained by the change of
variables $z=\cos\vartheta$. Functions $b(n,a)$ have been introduced in order
to make use of a recursion relation they satisfy. The latter is obtained by
integrating by parts twice in (107) with respect to $d(\rme^{az})$ to get

$b\left(n,a\right)=\frac{n-3}{a^{2}}\left[(n-5)b\left(n-4,a\right)-(n-4)b\left(n-2,a\right)\right].$
(108)

As is seen, functions $b(n,a)$ can be calculated recursively for all $n\geq 6$
if they are known for $n=2-5$. The latter are given by the terms in brackets
in the expressions calculated directly from (107):

$\displaystyle\rme^{-u^{(2)}({\bf x},t_{0})}$ $\displaystyle=$ $\displaystyle
C_{2}\rme^{-\frac{t_{0}a^{2}}{2}}\left[\pi I_{0}\left(a\right)\right]$
$\displaystyle\rme^{-u^{(3)}({\bf x},t_{0})}$ $\displaystyle=$ $\displaystyle
C_{3}\rme^{-\frac{t_{0}a^{2}}{2}}\left[\frac{2}{a}\sinh a\right]$
$\displaystyle\rme^{-u^{(4)}({\bf x},t_{0})}$ $\displaystyle=$ $\displaystyle
C_{4}\rme^{-\frac{t_{0}a^{2}}{2}}\left[\frac{\pi}{a}I_{1}\left(a\right)\right]$
$\displaystyle\rme^{-u^{(5)}({\bf x},t_{0})}$ $\displaystyle=$ $\displaystyle
C_{5}\rme^{-\frac{t_{0}a^{2}}{2}}\left[\frac{4}{a^{3}}\left(a\cosh a-\sinh
a\right)\right]$ (109)

where $I_{0,1}$ are the modified Bessel functions of the first kind.

Case $n=1$ is not covered by (107) but using (106) multiplied by 2 to accord
with the conventional definition of IM, the partially renormalized potential
can be straightforwardly calculated as

$u^{(1)}({x},t_{0})=\frac{t_{0}a^{2}}{2}-\ln(2\cosh a)+\frac{1}{2}\ln(2\pi
t_{0})$ (110)

where we took into account field-independent terms needed in D.

## References

  * [1] Wilson K G and Kogut J 1974 Phys. Rep. 12 75–199
  * [2] Wegner F J and Houghton A 1973 Phys. Rev. A 8 401–412
  * [3] Nicoll J F, Chang T S and Stanley H E 1976 Phys. Rev. A 13 1251–1264
  * [4] Polchinski J 1984 Nucl. Phys. B 231 269–295
  * [5] Tokar V I 1984 Phys. Lett. A 104 135–139
  * [6] Parola A, Pini D and Reatto L 1993 Phys. Rev. E 48 3321–3332
  * [7] Wetterich C 1993 Phys. Lett. B 301 90–94
  * [8] Bonini M, D’Attanasio M and Marchesini G 1993 Nucl. Phys. B 409 441–464
  * [9] Morris T R 1994 Int. J. Mod. Phys. A 09 2411–2449 ; arXiv:hep-ph/9308265
  * [10] Bagnuls C and Bervillier C 2001 Phys. Rep. 348 91–157
  * [11] Berges J, Tetradis N and Wetterich C 2002 Phys. Rep. 363 223 – 386
  * [12] Pelissetto A and Vicari E 2002 Phys. Rep. 368 549–727
  * [13] Dupuis N, Canet L, Eichhorn A, Metzner W, Pawlowski J M, Tissier M and Wschebor N 2021 Phys. Rep. 910 1–114
  * [14] Ellwanger U 1993 Z. Phys. C 58 619–627
  * [15] Machado T and Dupuis N 2010 Phys. Rev. E 82 041128
  * [16] Caillol J M 2012 Nucl. Phys. B 865 291–307
  * [17] Tokar V I 2021 J. Stat. Mech.: Theory Exp. 2021 013215
  * [18] Tokar V I 2019 Effective medium approach in the renormalization group theory of phase transitions (Preprint 1910.05123)
  * [19] Tokar V I 2021 Renormalization group approach to unified description of continuous and the first order phase transitions: application to the Blume-Capel model (Preprint 2103.12647)
  * [20] Bervillier C 2013 Nucl. Phys. B 876 587
  * [21] Bervillier C, Boisseau B and Giacomini H 2008 Nucl. Phys. B 789 525–551
  * [22] Caillol J M 2012 Nucl. Phys. B 855 854–884
  * [23] Parola A and Reatto L 2012 Mol. Phys. 110 2859–2882
  * [24] Morris T R and Tighe J F 1999 J. High Energy Phys. 1999 007–007
  * [25] Hasenbusch M 1999 J. Phys. A: Math. Gen. 32 4851–4865
  * [26] Berges J, Tetradis N and Wetterich C 1996 Phys. Rev. Lett. 77 873–876 arXiv:hep-th/9507159
  * [27] Berges J, Tetradis N and Wetterich C 1997 Phys. Lett. B 393 387–394 arXiv:hep-ph/9610354
  * [28] Rançon A and Dupuis N 2011 Phys. Rev. B 83 172501
  * [29] Rançon A and Dupuis N 2011 Phys. Rev. B 84(17) 174513
  * [30] Rançon A and Dupuis N 2012 Phys. Rev. A 85(6) 063607
  * [31] Rançon A, Kodio O, Dupuis N and Lecheminant P 2013 Phys. Rev. E 88(1) 012113
  * [32] Doyle J and Englefield M J 1990 IMA J. Appl. Math. 44 145–153
  * [33] Grossi E and Wink N 2019 Resolving phase transitions with discontinuous galerkin methods (Preprint 1903.09503)
  * [34] Pelissetto A and Vicari E 2002 Phys. Rep. 368 549–727
  * [35] De Polsi G, Balog I, Tissier M and Wschebor N 2020 Phys. Rev. E 101(4) 042113
  * [36] Shalaby A M 2021 Eur. Phys. J. C 81 87
  * [37] Zia R K P, Redish E F and McKay S R 2009 Am. J. Phys. 77 614–622
  * [38] Vasiliev A N 1998 Functional Methods in Quantum Field Theory and Statistical Physics (Amsterdam: Gordon and Breach)
  * [39] Gao Q and Zou M 2017 Appl. Math. Model. 45 255–270
  * [40] Hori S 1952 Prog. Theor. Phys. 7 578–584
  * [41] Kadanoff L P 1966 Physics 2 263–272
  * [42] Butera P and Comi M 2000 Phys. Rev. B 62 14837–14843
  * [43] Lundow P, Markström K and Rosengren A 2009 Phil. Mag. 89 2009–2042
  * [44] Dupuis N and Sengupta K 2008 Eur. Phys. J. B 66 271–278
  * [45] Nienhuis B and Nauenberg M 1975 Phys. Rev. Lett. 35 477–479
  * [46] Fisher M E and Berker A N 1982 Phys. Rev. B 26 2507–2513
  * [47] Whitham G 2011 Linear and Nonlinear Waves Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts (Wiley)
  * [48] Mousa M M 2015 Z. Naturforsch. A 70 47–58
  * [49] Brankov J G and Zagrebnov V A 1983 J. Phys. A Math. Gen. 16 2217–2224
  * [50] Salinas S 2001 Introduction to Statistical Physics Graduate Texts in Contemporary Physics (Springer New York)
  * [51] Choquard P and Wagner J 2004 Journal of Statistical Physics 116 843–853
  * [52] Morris T 2005 J. High Energy Phys. 0507 027
  * [53] Adams J, Tetradis N, Berges J, Freire F, Wetterich C and Bornholdt S 1995 Mod. Phys. Lett. A 10 2367–2379 arXiv:hep-th/9507093
  * [54] Wikimedia Foundation 2021 Spherical volume element https://en.wikipedia.org/wiki/N-sphere (accessed 2/10/2021)