potential $U(t,\sigma)$ with respect to the background field $\sigma$. It is
also possible to already start by directly deriving the flow equation for
$u(t,\sigma)$ via a projection on the one-point function
$\bar{\Gamma}_{t}^{(1)}(\sigma\,)$,

$\displaystyle\partial_{t}u(t,\sigma)=\vphantom{\Bigg{(}\Bigg{)}}$ (101)
$\displaystyle=\,$
$\displaystyle\big{(}\partial_{t}\bar{\Gamma}_{t}^{(1)}[\vec{\varphi}\,]\big{)}_{\varphi_{1}=\ldots=\varphi_{N-1}=0,\,\varphi_{N}=\sigma}=\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle=\,$
$\displaystyle\big{(}\delta_{\varphi_{N}}\,\partial_{t}\bar{\Gamma}_{t}[\vec{\varphi}\,]\big{)}_{\ldots}=\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle=\,$
$\displaystyle\Big{(}\delta_{\varphi_{N}}\,\mathrm{Tr}\Big{[}\big{(}\tfrac{1}{2}\,\partial_{t}R_{t}\big{)}\big{(}\bar{\Gamma}_{t}^{(2)}[\vec{\varphi}\,]+R_{t}\big{)}^{-1}\Big{]}\Big{)}_{\ldots}=\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle=\,$
$\displaystyle-\begin{gathered}\includegraphics{diagrams/derivativepotential_pi.pdf}\end{gathered}-\begin{gathered}\includegraphics{diagrams/derivativepotential_sigma.pdf}\end{gathered}=$
(104) $\displaystyle=\,$
$\displaystyle-\big{[}\tfrac{1}{2}\,\partial_{t}r(t)\big{]}\,\frac{N-1}{\big{[}r(t)+\frac{1}{\sigma}\,u(t,\sigma)\big{]}^{2}}\,\partial_{\sigma}\big{[}\tfrac{1}{\sigma}\,u(t,\sigma)\big{]}-\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle-\big{[}\tfrac{1}{2}\,\partial_{t}r(t)\big{]}\,\frac{1}{\big{[}r(t)+\partial_{\sigma}u(t,\sigma)\big{]}^{2}}\,\partial_{\sigma}^{2}u(t,\sigma)=\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle=\,$
$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\sigma}\,\left(+\begin{gathered}\includegraphics{diagrams/potential_pi.pdf}\end{gathered}+\begin{gathered}\includegraphics{diagrams/potential_sigma.pdf}\end{gathered}\right)\,.$
(107)

This corresponds to an interchange in the order of operations (evaluating the
FRG equation on the background field configuration and taking derivatives with
respect to the background field versus taking functional derivatives of the
FRG equation and afterwards evaluating on the background field) and it is non-
trivial (especially for flow equations for more complex models in higher
dimensions and with truncation beyond LPA) that the resulting equations are
identical.

Before we turn to the fluid-dynamical interpretation of the conservation laws
(99) and (100), we comment on the question whether one of the two formulations
(99) and (100) is preferable or even others should be considered. The answer
to this question is not yet settled. From our present understanding, a
formulation of the conservation equation in terms of $\sigma$ is preferable,
for reasons of numerical implementability. This is discussed at length in the
context of the PDE boundary conditions for the RG flow equation in Sub.Sec.
IV.4 and also in part III of this series of publications [31]. Therefore, our
discussion in the next sections is based on Eq. (100), and hence we identify
$\sigma$ with the spatial coordinate $x$ and
$u(t,\sigma)\equiv\partial_{\sigma}U(t,\sigma)$ as the conserved quantity.

####  IV.1.2 Advection-diffusion equation, irreversibility of RG flows, and
entropy production

This section is dedicated to the fluid-dynamical interpretation of the RG flow
equation (100). To this end, we split the flux (current) on the right-hand
side of the conservation law (100) and rewrite the whole equation in terms of
an advection-diffusion equation in one spatial dimension $x=\sigma$ and one
temporal dimension $t$,

$\displaystyle\partial_{t}u(t,x)+\tfrac{\mathrm{d}}{\mathrm{d}x}\,F[t,x,u(t,x)]=\tfrac{\mathrm{d}}{\mathrm{d}x}\,Q[t,\partial_{x}u(t,x)]\,.$
(108)

The pionic contributions to the RG flow,

$\displaystyle F[t,x,u(t,x)]=\,$
$\displaystyle-\big{[}\tfrac{1}{2}\,\partial_{t}r(t)\big{]}\,\frac{N-1}{r(t)+\frac{1}{x}\,u(t,x)}\,,\vphantom{\Bigg{(}\Bigg{)}}$
(109)

are identified with a non-linear, position-dependent advection flux, while the
contribution of the radial $\sigma$ mode,

$\displaystyle Q[t,\partial_{x}u(t,x)]=\,$
$\displaystyle+\big{[}\tfrac{1}{2}\,\partial_{t}r(t)\big{]}\,\frac{1}{r(t)+\partial_{x}u(t,x)}\,,\vphantom{\Bigg{(}\Bigg{)}}$
(110)

corresponds to a non-linear diffusion flux. This decomposition in $F$ and $Q$
can be understood as follows:

##### Advection:

If we ignore the contribution of the $\sigma$ mode for a moment (which – after
rescaling – corresponds to the large-$N$ limit of the $O(N)$ model [26, 28,
27, 31]), we can rewrite the left-hand side of Eq. (108) as follows,

$\displaystyle\partial_{t}u(t,x)+\tfrac{\mathrm{d}}{\mathrm{d}x}\,F[t,x,u(t,x)]=\vphantom{\bigg{(}\bigg{)}}$
(111) $\displaystyle=\,$
$\displaystyle\partial_{t}u(t,x)+\partial_{u}F[t,x,u(t,x)]\,\partial_{x}u(t,x)+\vphantom{\bigg{(}\bigg{)}}$
$\displaystyle+\partial_{x}F[t,x,u(t,x)]=0\,.\vphantom{\bigg{(}\bigg{)}}$

This is a non-linear advection equation for $u(t,x)$ in its primitive form
including an internal source term, where $\partial_{u}F[t,x,u(t,x)]$ is
identified with the velocity of the characteristics (the local $u$-dependent
flow velocity of the quantity $u$) and $\partial_{x}F[t,x,u(t,x)]$ acts like
an $x$ and $u$ dependent internal source term. Hence $F[t,x,u(t,x)]$ is not
purely advective nevertheless we will continue to refer to it as advection
term.

The identification of Eq. (111) as an advection equation is easily understood
by comparison with the one-dimensional linear advection equation,

$\displaystyle\partial_{t}u(t,x)+v\,\partial_{x}u(t,x)=0\,.$ (112)

Here $v$ corresponds to the a (constant) fluid velocity. Another prominent
example is the inviscid Bateman-Burgers equation [154, 155, 78, 79, 150, 156],
where the velocity is itself proportional to $u$.

However, for our RG flow, the local flow velocity is highly non-linear in $t$,
$x$ and $u$ and explicitly reads

$\displaystyle\partial_{u}F[t,x,u(t,x)]=\,$
$\displaystyle\frac{(N-1)\,\tfrac{1}{2}\,\partial_{t}r(t)}{x\,\big{[}r(t)+\frac{1}{x}\,u(t,x)\big{]}^{2}}\,.\vphantom{\Bigg{(}\Bigg{)}}$
(113)

Considering for example the exponential regulator shape function (8), one
finds that the advection velocity $\partial_{u}F[t,x,u(t,x)]$ is always
negative (positive) for $x>0$ ($x<0$). In a fluid-dynamical picture, this
means that the conserved quantity $u(t,x)$ is always propagated from larger
values of $|x|$ towards the point $x=0$ by advection. Furthermore, the closer
the fluid $u(t,x)$ is to $x=0$, the faster the fluid moves, due to the factor
$\frac{1}{x}$. Since $u(t,x)$ is anti-symmetric in $x$ [because of the $O(N)$
symmetry of $U(t,\vec{\varphi}\,)$], this implies that “waves” of positive and
negative $u(t,x)$ collide with huge velocity at $x=0$ and annihilate. At large
$|x|$, the fluid velocity tends to zero.

We also observe that the advection velocity (113) is proportional to the
number of pions, $N-1$. Hence, in the large-$N$ limit the system is completely
advection driven, while for small $N$ the diffusive contributions (110) gain
in importance. In the case $N=1$, there is no advection at all and the
dynamics of the fluid $u(t,x)$ is purely diffusive. Both limiting cases $N=1$
and $N\rightarrow\infty$ are discussed at length in the subsequent parts II
and III in this series of publications [33, 31] respectively.

It is also well known [78, 79, 150, 156] that systems which involve non-linear
advection tend to exhibit non-analytical behavior in the form of shock and
rarefaction waves. We will return to this issue below.

##### Diffusion:

Next, we turn to the contribution of the radial $\sigma$ mode to the RG flow.
We find that t enters the conservation law (108) as a non-linear diffusion
flux (110), because it is overall of second order in spatial derivatives of
$u(t,x)$. The characteristic property of diffusive processes is that they
transport a quantity, in this case $u(t,x)$, from regions where its density or
concentration is high to regions where it is low [78, 79, 150, 156]. Diffusive
processes are therefore usually important in regions of high gradients and
smear out cusps, shocks etc., which might form via advection. Besides this,
diffusive processes are generically undirected, which is also the case for
(110) and therefore also propagate the quantity $u(t,x)$ in both directions,
depending on the local gradients of $u(t,x)$, which is especially relevant for
models in their symmetry-broken phase with rather weak advection (small $N$).
The effective transport velocities via diffusion are usually much slower than
those via advection, which is, due to the non-linearity, not necessarily true
for RG flow equations. A famous example of a purely diffusive process is heat
transport via the heat equation [77, 78, 79]. The diffusion flux (108) can
indeed be formulated as a non-linear time-dependent realization of the heat
equation. By performing the spatial derivative in the advection-diffusion
equation (108) for the purely diffusive ($N=1$) case, one finds

$\displaystyle\partial_{t}u(t,x)=\,$
$\displaystyle-\frac{\tfrac{1}{2}\,\partial_{t}r(t)}{[r(t)+\partial_{x}u(t,x)]^{2}}\,\partial_{x}^{2}u(t,x)\,,$
(114)

where

$\displaystyle-\frac{\tfrac{1}{2}\,\partial_{t}r(t)}{[r(t)+\partial_{x}u(t,x)]^{2}}\,,$
(115)

plays the role of a non-linear time dependent, strictly positive diffusion
coefficient. The positivity of the diffusion coefficients ensures that
$u(t,x)$ is only dispersed and never accumulates locally, i.e., that $u(t,x)$
tends to equilibrate towards a linear function in space. A positive diffusion
coefficient also ensures stability and uniqueness of (numerical) weak
solutions, see, e.g, [78, 79, 157, 150].

Directly comparing these findings with the linear heat equation with constant
diffusion coefficient, we can already qualitatively predict the behavior of
the diffusion transport for the RG flow of $u(t,x)$, as long as $N$ is small
and the system is diffusion-dominated. At a constant RG time $t$, we find that
the diffusion coefficient is much larger in regions where the gradient
$\partial_{x}u(t,x)$ is negative with a large absolute value, compared to
regions where it is positive, because in the first case the denominator of Eq.
(115) is smaller than in the second case. This plays a crucial role for
systems that involve symmetry breaking, where $\partial_{x}u(t,x)$ is negative
for at least some small $|x|$, while asymptotically for $|x|\rightarrow\infty$
the sign of $\partial_{x}u(t,x)$ is always positive. Hence, for diffusion-
dominated problems in RG flow equations (small number $N$ of pions), the
symmetry restoration is driven by the negative gradients $\partial_{x}u(t,x)$
at small $|x|$. Furthermore, we find that for $t\rightarrow\infty$, the
numerator of the diffusion coefficient (115) tends to zero such that the
diffusion stops, the system equilibrates and the dynamics freezes, even though
there are still gradients in $u(t,x)$. This would not happen for the linear
heat equation. The same is true for $t=0$, where the diffusion coefficient is
suppressed by $1/\Lambda$. However, the $t$ dependence of the diffusion
coefficient of $O(1)$ models strongly depends on the space-time dimension,
which will be discussed elsewhere.

##### Irreversibility and entropy production:

In a fluid-dynamical setting, it is very easy to understand the role of the
radial $\sigma$ mode: Due to its diffusive character, it is directly
responsible for the irreversibility of the RG flow and RG transformations in
general. Diffusion is a particular example of a dissipative process, which is
irreversible and increases the entropy of the system191919Interestingly, Ref.
[147] comes to the same conclusion arguing in reverse order: “Some of the
information on the ultraviolet behavior of the field theory is lost under
renormalization transformations with $t>0$, since in the field theory it is
not legitimate to examine correlations at scales smaller than the cutoff. We
would therefore expect that a motion of the space $Q$ [a change of the set of
all couplings] under the influence of the renormalization group would become
an ‘irreversible’ process, similar to the time evolution of dissipative
systems.” We remark that also Ref. [148] stated that a term of second order in
field space derivatives in related RG flow equations “[$\ldots$] corresponds
to a dissipation in the flow and is responsible for the semi-group property of
the RG.”. The dissipative and irreversible character can be seen as a
“thermodynamic” version of the irreversible Kadanoff block-spin
transformations [158, 109, 104]. Hence, the dissipation clearly singles out
the RG time $t$ as a temporal direction, because it introduces a
“thermodynamic arrow of time” and “thermodynamic time asymmetry” via entropy
production [159]. This also explains why

$\displaystyle
t=-\ln\big{(}\tfrac{k}{\Lambda}\big{)}\vphantom{\bigg{(}\bigg{)}}$ (116)

(including the minus sign!) is a natural choice for a temporal coordinate in
higher dimensions, see also Refs. [152, 160, 148, 149, 26].

Interestingly, the irreversibility and the dissipative character of the system
is lost if one does not include the full field dependence of the effective
potential in the flow equation, but instead uses a truncated system like the
Taylor expansion (97). Then, the system of coupled ODEs for the vertices can
theoretically be integrated in either direction in RG time, as long as it
consists of a finite number of couplings202020In momentum space this enables
an integration to higher energy scales, which corresponds to a reversion of
the coarse-graining in position space. More generally speaking, this implies
that it is possible to resolve the microphysics from the macrophysics. Both is
physically not possible and solely an artifact of the truncation.. The most
extreme examples are the RG flows of one single $t$ dependent coupling, e.g.,
the quartic coupling of $\phi^{4}$ theory or the QCD $\beta$ function [161,
162, 163, 164], see also the textbooks [40, 39]. Here the integration to both
higher and smaller RG scales is possible, which is the well-known result for
the universal one-loop $\beta$ function and is an artifact of the restriction
(truncation) to a finite number of couplings [109]. However, this
reversibility of RG transformations is not possible for the field-dependent
effective potential, which is obvious from the advection-diffusion equation
(108), where entropy increases and the information about the initial condition
in the UV cannot be recovered from the IR anymore.

This point of view was already shared, presented, and discussed by K. G.
Wilson: In Ref. [109] he pointed out the differences between his “coarse-
graining” version of the (F)RG, which is also applicable in highly non-
perturbative regimes, and the RG flow equations used by C. Callan, K.
Symanzik, M. Gell-Mann, F. Low, G. t’Hooft, S. Weinberg, H. Georgi, D.
Politzer et al. to calculate the running of a single (or small number of)
coupling constants, which solely describes a system correctly in a
perturbative regime.

The irreversibility of the RG flow and entropy production is also directly
related to the presence of discontinuities in the solution, which can arise
from the advective contributions to the flow. As shown in Refs. [118, 26, 27,
31] for the large-$N$ limit, a shock wave arises when the weak solution of the
PDE is multi-valued. The correct solution is usually constructed by means of
the Rankine-Hugoniot condition [165, 166, 78, 79, 150, 156]. This would lead
to ambiguities when one tries to invert the flow (integrating backwards in
time) in the presence of a shock. Hence, shock formation is an irreversible
process and produces entropy. In summary, these are further strong arguments
why the assumption of expandability of the effective average action in terms
of vertices as well as the truncation of the system should in general be
considered with care.

Therefore, it would be extremely interesting to explicitly construct an
entropy function for the flow equation, i.e., a quantity that is either non-
decreasing or non-increasing under the RG transformations during the RG flow
(depending on the sign convention), and that is a functional of the quantity
$u(t,x)$. The entropy for the flow equation will be a helpful instrument to
design a stable numerical scheme for generic truncations [78, 79, 150] and
will also highlight general properties of the RG flow. In this context we also
have to mention the recent publication [93] by J. Cotler and S. Rezchikov who
were able to interpret the Polchinski equation as an “optimal transport
gradient flow of a field-theoretic relative entropy” thus establishing a firm
and explicit connection between an information-theoretic entropy and (F)RG
flows.

Additionally, a numeric entropy (function) might provide a direct link to the
$\mathcal{C}$ theorem ($\mathcal{A}$ theorem) [147, 167, 168, 169, 170, 171,
172, 74], which states that in certain QFTs there exists some positive real
function $\mathcal{C}(\\{g_{i}\\},t)$, which depends on all coupling constants
of the QFT and which is monotonically increasing212121It can also be defined
as a monotonically decreasing function. This flip of sign corresponds to the
difference of the mathematicians’ and physicists’ definition of entropy. We
chose to the “thermodynamic convention” of increasing entropy for this and
subsequent publications.  during RG flows (transformations), while it stays
constant at (critical) fixed points,

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}\,\mathcal{C}(\\{g_{i}\\},t)\geq
0\,.$ (117)

Here, $\\{g_{i}\\}$ denotes the set of all (possibly infinitely many)
dimensionless coupling constants. In contrast to previous formulations [173,
174, 175, 176, 177, 178, 179], a non-local version, which is directly linked
to the numerical entropy function (similar to versions presented in Refs.
[160, 148, 74] for related field-dependent flow equations), would not rely on
expandability in the couplings or vertices and could naturally display the
dissipative character of RG transformations, which was already described by
Refs. [147, 148]. Fixed-point solutions of the RG flow would directly
correspond to steady-state or thermal-equilibrium solutions [78] in the fluid-
dynamical picture222222This actually brings up the interesting question
whether previous studies about global fixed-point solutions for field-
dependent flow equations, which seemed to deliver interesting results, e.g.,
Refs. [180, 148, 91], should be reanalyzed from the fluid-dynamical steady-
flow perspective, especially regarding their interpretation and the spatial
discretization methods [78]. . A caveat at this point is that a $\mathcal{C}$
function is based on the rescaled dimensionless RG flow equations. Hence, also
a numerical entropy should be formulated in this framework, if one seeks a
direct link to a $\mathcal{C}$ function. The dimensionless flow equations in
the LPA can be recast in terms of conservation laws, which might be a good
starting point.

An explicit discussion of (numerical) entropy for the zero-dimensional $O(1)$
model as well as possible links to $\mathcal{C}$ functions is discussed in
great detail in part II of this series of publications [33]. The situation for
the $O(N)$ model in the limit $N\rightarrow\infty$ is discussed in App. E of
part III in this series of publications [31]. The construction of an explicit
(numerical) entropy has proven to be elusive in the case of finite $N>1$ for
the $O(N)$ model [33, 31] due to the explicit position dependences in Eqs.
(100) and (99) and the related internal source terms, cf. Eq. (111).

##### Generalizations:

Before we turn to the numerical implementation of the RG flow equation for
$u(t,\sigma)$, we briefly comment on the generalization of the fluid-dynamical
picture to RG flow equations in higher-dimensional QFTs, systems with more
(field-dependent) couplings, and RG flow equations that involve fermions.

In higher-dimensional QFTs, the fluid-dynamical interpretation of the RG flow
of the effective potential survives, see for example Ref. [26, 27, 140, 32]. A
difference is that in higher dimensions the RG time enters as the negative
logarithm of the ratio of the RG scale $k$ and the UV reference scale
$\Lambda$, see Eq. (116), while in zero dimensions $t$ just parametrizes some
mass-like scale $r(t)$, see Eq. (8). Furthermore, the fluxes gain further $t$
dependent prefactors via the momentum integrals of the trace in the FRG
equation. This leads to a different time scaling but does not affect the
overall discussion. The inclusion of further field-independent but scale-
dependent couplings (such as a scale-dependent Yukawa coupling) adds ODEs to
the advection-diffusion equation for the effective potential, which does not
spoil its conservative fluid-dynamical character. It is currently investigated
by us and collaborators [141, 27] whether the inclusion as well as the
conservative formulation of further field-dependent couplings (such as a
field- and scale-dependent wave-function renormalization
$Z(t,\vec{\varphi}\,)$ in higher-dimensional models) is possible. However,
this analysis is beyond the scope of the present work and will be presented
elsewhere. In any case, simply adding fermions in the LPA does not destroy the
fluid-dynamical character of the RG flow equation at all: On the level of the
LPA for the RG flow equation of the effective potential, the contributions
from fermion loops can be interpreted as a source/sink term, which only
depends on $\sigma$, i.e., the spatial position $x$. For a detailed discussion
of such fermionic source/sink terms at zero and non-zero temperature and
especially quark chemical potential we refer the interested reader to Ref.
[32]. Another possible generalization concerns models with more than one
invariant of the underlying symmetry group of the model and respective
condensation directions in field space, see, e.g., Refs. [181, 182, 183, 184,
185, 186]. Here, the fluid-dynamical framework should still be applicable.
However, a suitable identification of a complete basis of field-space
directions with “spatial directions” of the fluid-dynamical problem and a
clear separation of the single contributions into advection, diffusion, and
source terms might be challenging and calls for future investigations –
especially when it comes to an actual numerical implementation. For first
attempts of generalizing our findings to a quark-meson-diquark model, we refer
to Ref. [187].

Summarizing we find that the fluid-dynamical interpretation of flow equations
has tremendous benefits, because it allows for a rather intuitive
understanding of the dynamics of the system. Furthermore, it allows for a
novel, physically intuitive interpretation of the RG flow and provides an
understanding of its irreversibility. Finally, it opens up the opportunity to
employ extremely powerful numerical tools from computational fluid dynamics,
which are discussed in the next two sections. Phenomena discussed within this
section are observed, quantified, and visualized in Sec. V, in recent and
parallel works [26, 27, 33], and in parallel and upcoming publications [153,
32, 34].

###  IV.2 Finite-volume method

In this section we discuss numerical solution schemes for advection-diffusion
equations232323Oftentimes, such equations are also referred to as “convection-
diffusion equations”. The semantically correct term is nevertheless
“advection-diffusion equation” because “convection” includes also diffusive
processes besides the transport by bulk motion (advection), see also Ref.
[78]. with source terms of the generic type

$\displaystyle\partial_{t}u(t,x)+\tfrac{\mathrm{d}}{\mathrm{d}x}\,F[t,x,u(t,x)]=\tfrac{\mathrm{d}}{\mathrm{d}x}\,Q[t,x,u(t,x),\partial_{x}u(t,x)]+S[t,x,u(t,x)]\,.$
(118)

Whether $S$ acts as a source or sink in the dynamics of $u(t,x)$ depends on
its explicit form. Nevertheless we will refer to $S$ as source term for
convenience for the scope of this paper. In the following, we occasionally
suppress the $t$ and $x$ dependences of $u$, $F$, $Q$, and $S$ for the sake of
simplicity. Equation (118) is a partial differential equation describing the
evolution of the conserved quantity $u\equiv u(t,x)$ in one spatial ($x$)
direction and one temporal ($t$) direction. Depending on the problem at hand
these two directions are not necessarily identical with physical spatial and
temporal dimensions of reality, but for the following discussion we denote
them as such. The function $F[u]\equiv F[t,x,u(t,x)]$ is a (non-linear)
advection flux, $Q[u,\partial_{x}u]\equiv Q[t,x,u(t,x),\partial_{x}u(t,x)]$ is
a (non-linear) diffusion(dissipation) flux and $S[u]\equiv S[t,x,u(t,x)]$ is a
source term. The concepts discussed in the following apply directly to systems
of $M$ conserved quantities $u\rightarrow\\{u_{1},\ldots,u_{M}\\}$ and can be
generalized beyond one spatial dimension to $d+1$ dimensional space-time
$(x,t)\rightarrow(\vec{x},t)=(x_{1},\ldots,x_{d},t)$. Equations or systems
similar or even identical to Eq. (118) are often referred to as conservation
laws and appear in many areas of the natural sciences, engineering, and
economics. They are extensively studied in the field of computational fluid
dynamics.

Consider the advection-diffusion equation (118) with specified $F$, $Q$, and
$S$ in a finite computational domain $\Omega=\mathcal{V}\times[t_{0},t_{N}]$,
where $\mathcal{V}\subset\mathbb{R}^{1}$ denotes the spatial volume, with an
initial condition $u(t_{0},x)$ and Dirichlet (Neumann) boundary conditions
specifying $(\partial_{x})u(t,x)|_{x\in\partial\mathcal{V}}$. The question is
how to evolve the initial condition in time from $t_{0}$ to $t_{N}>t_{0}$ to
acquire a solution $u(t_{N},x)$ respecting the specified boundary conditions.
For most problems of the type (118) an analytic solution is not known or is
even known to be nonexistent. Strategies for finding numerical (weak)
solutions are required. Numerical schemes in the broad class of so-called
finite-volume (FV) methods are very popular for the numerical solution of PDEs
describing the conservation or balance of quantities. Alternative high-
resolution shock-capturing (HRSC) schemes in modern computational fluid
dynamics are among others finite-difference schemes including flux limiters
and numerical viscosity or finite-element methods.

The concept that all numerical FV methods share is a discretization of the
computational domain into space-time control volumes
$\mathcal{V}_{j}\times[t^{n},t^{n+1}]$, where the set of spatial control
volumes $\mathcal{V}_{j}$ covers the spatial computational domain
$\mathcal{V}$. Integrating Eq. (118) over such a control volume centered at
$x$, using the divergence theorem (Gauss-Ostrogradsky theorem) on the fluxes
and introducing the sliding cell average

$\displaystyle\bar{u}(t,x)\equiv\frac{1}{|\mathcal{V}_{j}|}\int_{\mathcal{V}_{j}}\mathrm{d}\xi\,u(\xi,t)\,,$
(119)

where $\mathcal{V}_{j}=\\{\xi:|\xi-x|\leq\Delta x/2\\}$ we arrive at an
equivalent integral form of Eq. (118),

$\displaystyle\bar{u}(t^{n+1},x)=\vphantom{\bigg{(}\bigg{)}}$ (120)
$\displaystyle=\,$ $\displaystyle\bar{u}(t^{n},x)-\tfrac{1}{\Delta
x}\bigg{(}\int_{t^{n}}^{t^{n+1}}\mathrm{d}\tau\,F\big{[}\tau,x+\tfrac{\Delta
x}{2},u\big{(}\tau,x+\tfrac{\Delta
x}{2}\big{)}\big{]}-\int_{t^{n}}^{t^{n+1}}\mathrm{d}\tau\,F\big{[}\tau,x-\tfrac{\Delta
x}{2},u\big{(}\tau,x-\tfrac{\Delta
x}{2}\big{)}\big{]}\bigg{)}+\ldots\,.\vphantom{\bigg{(}\bigg{)}}$

Here, the ellipsis denotes the corresponding integrals for the diffusion and
the source term. The solution of Eq. (120) presents the central challenge for
an explicit FV scheme.

A central aspect of each practical FV scheme is an appropriate and informed
choice of the space-time control volumes which, depending on the scheme and
problem at hand, might change during the time evolution. Given a set of
control volumes and a corresponding set of cell averages
$\bar{u}(t^{n},x_{j})\equiv\bar{u}^{n}_{j}$ the time evolution to
${t^{n+1}\equiv t^{n}+\Delta t}$ requires the solution of the Riemann problems
[188, 78, 79, 157, 150] at each cell interface. Part of these problems are the
fluxes through the cell boundaries. The computation of those fluxes requires a
reconstruction of the values of $u$ on the cell interfaces located at
$x_{j+\frac{1}{2}}$, which we denote as $u_{j+\frac{1}{2}}^{n}$, from the
given set of cell averages $\bar{u}_{j}^{n}$. This is usually done by means of
a polynomial approximation respecting the given cell averages of the
neighboring cells. The order of the chosen approximation is one of the
parameters contributing to the overall spatial order (of the error) of the
scheme at hand.

Given the cell averages $\bar{u}_{j}^{n}$ and fluxes through the cell
interfaces at $t=t^{n}$ it remains to solve the Riemann problems at the cell
interfaces. The solution of the Riemann problem amounts to the exact
evaluation of the flux integrals on the right-hand side of Eq. (120).
Depending on the complexity of the underlying conservation equation an exact
solution of the Riemann problems at the cell boundaries might be either
impossible or unfeasible. Most explicit FV schemes, especially those for
general advection-diffusion equations, either use approximate Riemann solvers
(e.g., the Roe [189] or the HLLE [190, 191] solver) or do not require Riemann
solvers at all (e.g., the KT [192] scheme). For a pedagogic introduction into
the broad field of FV methods and HRSC schemes in general we refer the
interested reader to Refs. [156, 78, 79, 157, 150] and references therein.

In the following section we will introduce a particular FV scheme, which we
have chosen for the numerical solution of the RG flow equations because of its
flexibility, efficiency, and relative simplicity.

###  IV.3 Kurganov-Tadmor (KT) central scheme

In this section we will summarize the central scheme presented in Ref. [192]
by A. Kurganov and E. Tadmor, which we will refer to in the following as KT
scheme. The KT scheme can be implemented and applied as a black-box solver for
systems of the type of Eq. (118). Apart from the PDE with its initial and
boundary conditions the only additional information about the PDE required for
its solution using the KT scheme is $\big{|}\frac{\partial F}{\partial
u}\big{|}$ or the spectral radius of the Jacobian of $F[u]$ when considering
systems of conserved quantities, see Eq. (124) and the related discussion. The
scheme does not require a Riemann solver of any kind and as such does not rely
on a characteristic decomposition of the advection flux.

The KT scheme provides a direct method for evaluating the flux integrals on
the right-hand side of Eq. (120). The main focus lies on the treatment and
implementation of the flux integrals for the advection flux $F[u]$, which is
kind of a natural approach, because the advection flux determines the
characteristic velocities in an advection-diffusion equation. The diffusion
and source terms are treated separately and will be discussed at the end of
this subsection.

The KT scheme admits a meaningful ${t^{n+1}-t^{n}\equiv\Delta t\rightarrow 0}$
limit in the context of Eq. (120) and is thus an improvement on it predecessor
the Nessyahu-Tadmor (NT) scheme [193] with which it shares its piecewise-
linear MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws)
reconstruction [194]. We will focus on the KT scheme in its so called semi-
discrete from – in the limit $\Delta t\rightarrow 0$ – which involves only an
explicit spatial discretization. The KT scheme is formally second-order
accurate in the spatial direction and as such an improved version of the
first-order accurate Lax–Friedrichs (LxF) scheme [195, 196]. A semi-discrete
form reduces the PDEs (118) or equivalently (120) to a set of coupled ODEs,
which can be solved by a large class of general-purpose ODE solvers. This is
especially useful when working on stiff problems or PDE systems coupled to
additional ODEs. We will proceed with the introduction of quantities involved
in the semi-discrete form (135) of the KT scheme. The following quantities are
especially relevant for the numerical advection flux (131).

Consider a set of volume averages $\bar{u}_{j}^{n}$ at $t^{n}$ based on an
equidistant242424The generalization of the KT scheme to non-uniform grids is
on a conceptual level straightforward and especially useful for higher-
dimensional extensions and for adaptive or moving mesh variants, see, e.g.,
Ref. [197]. Its implementation is of course much more challenging and not
needed in this work. However, in the context of FRG flow equations this might
be relevant for models with multiple condensate directions, see, e.g., Refs.
[181, 182, 183, 184, 185, 186].  grid of volume cells
$\mathcal{V}_{j}\equiv[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$, with $\Delta
x=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$.

The time evolution of the averages $\bar{u}_{j}^{n}$ at $t^{n}$ to averages at
$\bar{u}_{j}^{n+1}$ at $t^{n+1}$ on the same volume grid is a three-step
process:

  1. 1.

The piecewise-linear MUSCL reconstruction is computed from the cell averages:

$\displaystyle\tilde{u}(t^{n},x)=\vphantom{\bigg{(}\sum_{j}\bigg{)}}$ (121)
$\displaystyle=\,$
$\displaystyle\sum_{j}\big{\\{}\bar{u}_{j}^{n}+(\partial_{x}u)_{j}^{n}\,(x-x_{j})\big{\\}}\openone_{[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]}\,,\vphantom{\bigg{(}\sum_{j}\bigg{)}}$

where the sum runs over all volume cells. The reconstruction step is needed to
gain access to the function values $\tilde{u}(t^{n},x)$. The reconstruction
uses approximations to the exact derivatives $(\partial_{x}u)_{j}^{n}$ by
employing a scalar total variation-diminishing (TVD) reconstruction [78, 79,
198],

$\displaystyle(\partial_{x}u)_{j}^{n}=\frac{\bar{u}_{j+1}^{n}-\bar{u}_{j}^{n}}{\Delta
x}\,\phi\bigg{(}\frac{\bar{u}_{j}^{n}-\bar{u}_{j-1}^{n}}{\bar{u}_{j+1}^{n}-\bar{u}_{j}^{n}}\bigg{)}\,,$
(122)

with a TVD limiter $\phi(r)$. An overview of TVD flux limiters can be found,
e.g., on the webpage [199], in Refs. [78, 79], or in Sec. 9.3.1 of Ref. [150].
Here, we follow Ref. [192] and use the so-called minmod limiter [200]252525We
also implemented and tested other flux limiters in the context of FRG flow
equations, which however did not influence our numerical results very much. An
optimization of the choice of flux limiters with regard to the overall runtime
could be part of future work.,

$\displaystyle\phi(r)=\,$ $\displaystyle\max[0,\min(1,r)]\,.$ (123)

The limiter $\phi$ is used in Eq. (122) to limit the slopes during the
reconstruction process. This is crucial to prevent spurious oscillations
around discontinuities, e.g., shocks, in high-resolution schemes like the KT
scheme. The KT scheme can also be used with higher-order reconstruction
schemes262626Examples for such improvements are the use of the third-order
central weighted essentiallynon-oscillatory (C-WENO) reconstruction [201, 202]
outlined in Ref. [203], the fifth-order WENO scheme (WENO5) [204, 205]
employed in Ref. [206], or the fifth-order monotonicity-preserving (MP5)
reconstruction [207] used in Ref. [208]. WENO schemes were first introduced in
Ref. [209] and use polynomial reconstructions based on smoothness indicators.
In regions where the solution is estimated to be smooth high-order polynomials
are used, while in regions with large gradients lower-order polynomials are
employed in order to ensure an essentially non-oscillatory [201] interpolation
and a TVD reconstruction. The MP5 scheme was first introduced in Ref. [207]
and uses a limiting procedure together with a fourth-order polynomial
reconstruction. The MP5 reconstruction does not require smoothness indicators
and performs favorable when compared to WENO in terms of speed and accuracy
according to Ref. [210]. to increase the spatial accuracy of the scheme, which
is not needed for this work.

When using a piecewise-constant or -linear reconstruction the cell averages
$\bar{u}_{j}^{n}$ coincide with the midpoint values $u_{j}^{n}$. While we
employ a piecewise-linear reconstruction, we still maintain the distinction
between averages and midpoint values for the sake of clarity.

  2. 2.

The time step from $t^{n}$ to $t^{n+1}$ is performed by computing the flux
integrals on the right-hand side of Eq. (120) using the reconstruction
$\tilde{u}(t^{n},x)$ from Eq. (121) and carefully chosen control volumes
discussed below. In the limit ${t^{n+1}-t^{n}\equiv\Delta t\rightarrow 0}$
only the expressions for $a_{j+\frac{1}{2}}^{n}$, $u_{j+\frac{1}{2}}^{n,-}$,
and $u_{j+\frac{1}{2}}^{n,-}$ from Eqs. (124), (125), and (126) respectively
are relevant for the semi-discrete KT scheme. The other quantities discussed
for this second step of the KT scheme are however necessary to understand the
underlying algorithm.

At each cell interface $x_{j+\frac{1}{2}}$ the respective local speed of
propagation $a_{j+\frac{1}{2}}^{n}$ is estimated using

$\displaystyle a_{j+\frac{1}{2}}^{n}\equiv\vphantom{\bigg{(}\bigg{)}}$ (124)
$\displaystyle\equiv\,$ $\displaystyle\max\bigg{\\{}\bigg{|}\frac{\partial
F}{\partial
u}\Big{[}u_{j+\frac{1}{2}}^{n,+}\Big{]}\bigg{|},\bigg{|}\frac{\partial
F}{\partial u}\Big{[}u_{j+\frac{1}{2}}^{n,-}\Big{]}\bigg{|}\bigg{\\}}\,,$

with the left and right intermediate values $u_{j+\frac{1}{2}}^{n,\mp}$ of
$\tilde{u}(t^{n},x)$ at the cell interface $x_{j+\frac{1}{2}}$:

$\displaystyle u_{j+\frac{1}{2}}^{n,-}=\,$
$\displaystyle\bar{u}_{j}^{n}+\tfrac{\Delta
x}{2}\,(\partial_{x}u)_{j}^{n}\,,\vphantom{\bigg{(}\bigg{)}}$ (125)
$\displaystyle u_{j+\frac{1}{2}}^{n,+}=\,$
$\displaystyle\bar{u}_{j+1}^{n}-\tfrac{\Delta
x}{2}\,(\partial_{x}u)_{j+1}^{n}\,.\vphantom{\bigg{(}\bigg{)}}$ (126)

For multi-valued conserved quantities $u=\\{u_{1},\ldots,u_{M}\\}$ the KT
scheme of Ref. [192] is limited272727The KT scheme of Ref. [192] can be
improved using refined estimates for the local speed of propagation. In Ref.
[211] the KT scheme is refined by employing left- and right-sided local speeds
$a_{j+\frac{1}{2}}^{n,\mp}$, which are extracted from the eigenvalue spectrum
of the Jacobian $\frac{\partial F}{\partial u}$. A further improvement in
terms of estimates of local speeds of propagation engineered for non-convex
hyperbolic (systems of) conservation laws is presented in Ref. [206] using
further information about the eigensystem of the Jacobian $\frac{\partial
F}{\partial u}$. When an explicit evaluation of the Jacobian is impossible or
unfeasible numerical approximations can be employed [192, 212, 213]. For the
numerical applications in this paper, the rather simple estimate of Eq. (124)
has proven to be sufficient. For specific computations in part III in this
series of publications [31] we additionally employed the improved scheme of
Ref. [211]. to hyperbolic advection fluxes signaled by a non-degenerate
eigenvalue spectrum $\lambda_{1}<\ldots<\lambda_{M}$ of the Jacobian
$\frac{\partial F}{\partial u}$ for all $x$, $t$, and $u$, since in this case
the local speed of propagation is computed using the spectral radius
$\rho(M)\equiv\max_{i}|\lambda_{i}(M)|$ in Eq. (124).

Using the estimated local speed of propagation, a space-time control volume
${[x_{j+\frac{1}{2},l}^{n},x_{j+\frac{1}{2},r}^{n}]\times[t^{n},t^{n}+\Delta
t]}$ around each cell interface $x_{j+\frac{1}{2}}$ is chosen. The spatial
extent corresponds to the domain which is causally affected by information
propagating with the local velocities away from the cell interface at
$x_{j+\frac{1}{2}}$. The flux integrals of Eq. (120) are performed on these
space-time control volumes separately using the midpoint rule to approximate
the flux integrals and leading to averages $\bar{\omega}_{j}^{n+1}$ and
$\bar{\omega}_{j+\frac{1}{2}}^{n+1}$ based on the new intermediate spatial
grid spanned by the points

$\displaystyle x_{j+\frac{1}{2},l}^{n}=\,$ $\displaystyle
x_{j+\frac{1}{2}}-a_{j+\frac{1}{2}}^{n}\Delta t\,,\vphantom{\bigg{(}\bigg{)}}$
(127) $\displaystyle x_{j+\frac{1}{2},r}^{n}=\,$ $\displaystyle
x_{j+\frac{1}{2}}+a_{j+\frac{1}{2}}^{n}\Delta t\,.\vphantom{\bigg{(}\bigg{)}}$
(128)

In the regions $[x_{j-\frac{1}{2},r}^{n},x_{j+\frac{1}{2},l}^{n}]$ the
solutions underlying the computed averages $\bar{\omega}_{j}^{n+1}$ are smooth
while the solutions underlying the computed averages
$\bar{\omega}_{j+\frac{1}{2}}^{n+1}$ are non-smooth based on the regions
$[x_{j+\frac{1}{2},l}^{n},x_{j+\frac{1}{2},r}^{n}]$. Details of this step can
be found in Ref. [192].

  3. 3.

A MUSCL-type piecewise-linear reconstruction based on
$\bar{\omega}_{j+\frac{1}{2}}^{n+1}$ and $\bar{\omega}_{j}^{n+1}$ is used to
project these averages back onto the original uniform grid spanned by
$x_{j+\frac{1}{2}}$. This results in a fully discrete second-order central
scheme, see Eq. (3.9) of Ref. [192], which gives an algebraic expression for
$\bar{u}_{j}^{n+1}$ in terms of the averages

$\displaystyle\\{\bar{u}_{j-2}^{n},\bar{u}_{j-1}^{n},\bar{u}_{j}^{n},\bar{u}_{j+1}^{n},u_{j+2}^{n}\\}$
(129)

and $\\{a_{j\pm\frac{1}{2}}^{n}\\}$. A pictographic representation of the
multi-step evolution procedure with the involved quantities and grids can be
found in Fig. 3.2 of Ref. [192]. The numerical viscosity of this second-order
scheme is $\mathcal{O}(\Delta x^{3})$ and does not depend on $\Delta t$, which
is an improvement when compared to the $\Delta t$ dependent numerical
viscosities $\mathcal{O}(\Delta x^{2}/\Delta t)$ and $\mathcal{O}(\Delta
x^{4}/\Delta t)$ of the LxF and NT schemes, respectively [192].

The $\Delta t$ independent numerical viscosity allows for a controlled limit
$\Delta t\rightarrow 0$, resulting in a reduction to a practical semi-discrete
scheme in the conservative form [192], which can be implemented
straightforwardly:

$\displaystyle\partial_{t}\bar{u}_{j}=-\frac{H_{j+\frac{1}{2}}-H_{j-\frac{1}{2}}}{\Delta
x}+\ldots\,,$ (130)

where the ellipsis denotes the diffusion and source fluxes. The numerical
advection fluxes $H_{j+\frac{1}{2}}$ are given by

$\displaystyle H_{j+\frac{1}{2}}\equiv\,$
$\displaystyle\frac{F\big{[}t,x_{j+\frac{1}{2}},u_{j+\frac{1}{2}}^{+}\big{]}+F\big{[}t,x_{j+\frac{1}{2}},u_{j+\frac{1}{2}}^{-}\big{]}}{2}-\vphantom{\Bigg{(}\Bigg{)}}$
(131) $\displaystyle-
a_{j+\frac{1}{2}}\,\frac{u_{j+\frac{1}{2}}^{+}-u_{j+\frac{1}{2}}^{-}}{2}\,.\vphantom{\Bigg{(}\Bigg{)}}$

This semi-discrete scheme is second-order accurate in $\Delta x$ and can be
used in conjunction with various ODE time-step algorithms (in this work, we
use Mathematica’s NDSolve [138]). The KT scheme for a position-independent
advection flux is conservative, meaning detailed balance at the cell
interfaces is maintained. It is also total variation diminishing/non-
increasing (TVD/TVNI) [198, 78, 79]282828This property is important for the
resolution of shocks and discontinuities. A TVD scheme does not produce
spurious oscillations around discontinuities. For a TVD scheme
${\sum_{j}\big{|}\bar{u}_{j+1}^{n+1}-\bar{u}_{j}^{n+1}\big{|}\leq\sum_{j}\big{|}\bar{u}_{j+1}^{n}-\bar{u}_{j}^{n}\big{|}}$
holds between all time steps, which is why the more concise expression is
total variation non-increasing, cf. Sec. 9.2.2 of Ref. [150]. The total
variation
${\mathrm{TV}[\\{\bar{u}_{j}\\}]=\sum_{j}\big{|}\bar{u}_{j+1}-\bar{u}_{j}\big{|}}$
is a discrete measure for the arc-length of $u(t,x)$. In the subsequent parts
II and III in this series of publications [33, 31] we show that the TV can
serve as a (numerical) entropy for the zero-dimensional $O(N)$ model with
$N=1$ and $N\rightarrow\infty$ and might provide a link to the formulation of
the $\mathcal{C}$ function in truncated flow equations [147]. when used with
appropriate flux limiters like the minmod limiter (123).

So far we only considered the advection term $\partial_{x}F[u]$ in the
discussion of the KT scheme. The explicit treatment of source and diffusion
fluxes in the semi-discrete KT scheme is as follows:

  1. 1.

If the source term is independent of $u$ and of type
$S(t,x)=\partial_{x}s(t,x)$, it can be integrated directly and contributes
$(s_{j+\frac{1}{2}}-s_{j-\frac{1}{2}})/\Delta x$ with $s_{j+\frac{1}{2}}\equiv
s(t,x_{j+\frac{1}{2}})$ to the right-hand side of Eq. (130). Otherwise, if
$S[u]\equiv S[t,x,u(t,x)]$ suitable approximations, like
$S[t,x_{j},\bar{u}_{j}(t)]$ may be used. For a detailed discussion source/sink
terms arising in the FRG treatment of higher-dimensional QFTs we refer the
interested reader to Ref. [32] and especially App. E of the aforementioned
publication.

  2. 2.

When considering a non-linear diffusion flux $Q[u,\partial_{x}u]$ Eq. (118) is
a potentially strongly degenerate parabolic equation (system) admitting
potentially non-smooth solutions. In the KT scheme the hyperbolic and
parabolic part of the PDE (118) are treated simultaneously based on the strict
splitting between $F$ and $Q$. Kurganov and Tadmor [192] presented a
discretization of the diffusion flux based on a kind of central-difference
approximation,

$\displaystyle P_{j+\frac{1}{2}}=\,$
$\displaystyle\tfrac{1}{2}\,Q\big{[}t,x_{j},\bar{u}_{j},\tfrac{\bar{u}_{j+1}-\bar{u}_{j}}{\Delta
x}\big{]}+\vphantom{\bigg{(}\bigg{)}}$ (132)
$\displaystyle+\tfrac{1}{2}\,Q\big{[}t,x_{j+1},\bar{u}_{j+1},\tfrac{\bar{u}_{j+1}-\bar{u}_{j}}{\Delta
x}\big{]}\,.\vphantom{\bigg{(}\bigg{)}}$

An alternative second-order discretization like the one put forward in App. B
of Ref. [214] can also be successfully employed: If the diffusion flux can be
written as

$\displaystyle\qquad Q[t,x,u,\partial_{x}u]$
$\displaystyle=\varphi[t,x,u]\cdot\tilde{Q}[t,x,\partial_{x}u]\,.\vphantom{\bigg{(}\bigg{)}}$
(133)

then

$\displaystyle\qquad P_{j+\frac{1}{2}}=\,$
$\displaystyle\varphi[t,x_{j+\frac{1}{2}},\tfrac{1}{2}\big{(}u_{j+\frac{1}{2}}^{+}+u_{j+\frac{1}{2}}^{-}\big{)}]\cdot\vphantom{\bigg{(}\bigg{)}}$
(134)
$\displaystyle\quad\cdot\tilde{Q}\big{[}t,x_{j+\frac{1}{2}},\tfrac{\bar{u}_{j+1}-\bar{u}_{j}}{\Delta
x}\big{]}\,,$

which can be applied for the RG flow equation (100) under consideration with
$\varphi[t,x,u]=1$ and $\tilde{Q}[t,x,\partial_{x}u]=Q[t,\partial_{x}u]$,
since $Q$ has no explicit position dependence when using the formulation in
$\sigma$, cf. Eq. (110). Both approaches (using Eq. (132) or (134) for the
numerical diffusion flux $P_{j+\frac{1}{2}}$) were successfully employed in
the context of this work for the diffusion flux (110) without visible
differences in quality. Still, all numerical results in Sec. V are obtained
using the formulation (132). For improved KT-type schemes employing higher-
order reconstructions (like, e.g., C-WENO/WENO5 or MP5) higher-order
discretizations for the diffusion flux are used like the fourth-order one put
forward in Eqs. (4.9) and (4.10) of Ref. [203].

The full semi-discrete KT scheme including advection, diffusion, and source
flux is given by

$\displaystyle\partial_{t}\bar{u}_{j}=\vphantom{\bigg{(}}$ (135)
$\displaystyle=\,$
$\displaystyle-\frac{H_{j+\frac{1}{2}}-H_{j-\frac{1}{2}}}{\Delta
x}+\frac{P_{j+\frac{1}{2}}-P_{j-\frac{1}{2}}}{\Delta
x}+\mathrm{source}\,.\vphantom{\bigg{(}}$

Specific spatial boundary conditions
$(\partial_{x})u(t,x)|_{x\in\partial\mathcal{V}}$ manifest themselves in the
KT scheme in the choice of volume averages for so-called ghost cells. For the
evaluation of the fluxes for the first two ($\bar{u}_{0}$ and $\bar{u}_{1}$)
and last two ($\bar{u}_{n-2}$ and $\bar{u}_{n-1}$) cell averages ghost cells
outside the computational domain are required when considering $n$ volume
cells, see, e.g., Refs. [78, 79] for a detailed discussion. The ghost cells
are $u_{-2}$, $u_{-1}$, $u_{n}$, and $u_{n+1}$ centered at $x_{-2}$, $x_{-1}$,
$x_{n}$, and $x_{n+1}$, respectively. Depending on the problem and
computational domain at hand a multitude of different boundary conditions can
be implemented by an appropriate choice or reconstruction of these averages
for the ghost cells. The specific boundary conditions for the RG flow
equations in this work are discussed at length in Sub.Sec. IV.4 below.

At this point we have to remark that the original KT numerical scheme
presented in Ref. [192] was constructed for position- and time-independent
advection and diffusion fluxes. Since we employ the KT scheme in its semi-
discrete form a resolution of potentially highly complicated and non-linear
dynamics in $t$ is possible and ultimately outsourced to the ODE solver. The
spatial discretization of the advection term (109) in the semi-discrete KT
scheme seems to be able to cope with the explicit position dependence in
$F[t,x,u]$ when evaluating the numerical fluxes at the appropriate cell
interfaces, cf. Eq. (131). In the scope of this paper we could not trace any
practical problems back to the explicit position- and time-dependence of the
advection and diffusion fluxes. The comparisons in Sec. V between results
obtained from a direct computation of correlation functions using the
generating functional (77) and the results computed using RG flow equations
via the KT scheme (with $t$ and $x$ dependent fluxes) can be seen as hard
tests for both – the FRG methodology as well as the (slightly modified) KT
scheme – depending on the respective perspective. In total, the precision of
our results for the non-trivial test cases gives us some confidence that our
approach is generically justified and the KT scheme is suitable for our
purpose. As mentioned in the paragraph concerning entropy and irreversibility
of Sub.Sub.Sec. IV.1.2 the explicit position dependences in both RG flow Eqs.
(99) and (100) prevents a direct usage of the total variation TV as entropy
functional for the zero-dimensional $O(N)$ model at finite $N>1$. The internal
source terms in the primitive forms originating in the explicit position
dependences of $F$ and $Q$ before performing the spatial derivatives lead to
an increase in arc length during RG time evolution, rendering the TV useless
as an entropy functional. For more details see the subsequent parts II and III
in this series of publications [33, 31], especially App. E of part III.

###  IV.4 Boundary conditions and computational domain in FRG flow equations

In the form of the conservation law (99) or (100), the RG flow equation (84)
is a non-linear PDE which has contributions of parabolic (diffusion terms) as
well as hyperbolic type (advection terms). In this subsection, we specify the
boundary conditions for Eq. (99) or (100) in field space (the effective
spatial $x$ direction).

The correct implementation of field space (spatial) boundary conditions is
very important for the correct numerical solution of PDEs [78, 79]. For (non-
linear) PDEs of elliptic type, the problem of solving them is also referred to
as boundary-value problem. For (non-linear) PDEs of hyperbolic and parabolic
type, the spatial boundary conditions are needed (in addition to the initial
condition) to make finding a (weak) solution a well-defined problem. The
latter case is also referred to as Cauchy or initial-boundary-value problem.
Thus, without explicitly specifying the boundary conditions, e.g., of Neumann-
or Dirichlet-type, as well as the initial conditions, the problem of finding a
unique (weak) solution is actually ill-posed and therefore impossible to solve
– a well-known mathematical fact with particular and severe implications in,
e.g., classical electrodynamics [215], fluid dynamics [216], soliton and
instanton solutions of classical field equations [217, 218], general
relativity [219, 220, 221], and other fields of research. This also holds true
for the FRG. However, explicit boundary conditions and especially their
numerical implementation are rarely discussed the FRG literature, with ,e.g.,
Refs. [129, 130, 131, 176] as notable exceptions.

For the derivative of the effective potential $u(t,\sigma)$, we find that the
spatial boundary conditions must be imposed at $\sigma=\pm\infty$, because the
field space domain of $u(t,\sigma)$ is given by $\mathbb{R}$. Thus, when
considering the flow equation on the non-compact domain $(-\infty,\infty)$ the
problem represents a pure initial-value/Cauchy problem [156, 78, 79] and,
given the asymptotics of the flow equation and the initial condition, explicit
boundary conditions at $x\rightarrow\pm\infty$ are not required. However,
spanning a non-compact computational interval from $-\infty$ to $+\infty$ is
practically impossible on a finite computational grid. A possible solution is
a compactification [137] of $\mathbb{R}$ to the interval $[-1,+1]$, via a
suitable mapping $\sigma\mapsto x(\sigma)$ usually supplemented with a mapping
$u\mapsto v(u)$ rendering $v$ finite on $[-1,+1]$. Another popular solution is
a truncation of the computation interval at a large value
$\sigma_{\text{max}}\sim x_{\mathrm{max}}$ with a suitable boundary condition
[130, 129, 180, 137]. We will return to this issue below.

In any case, one of the boundaries at spatial infinity can already be replaced
by a finite value by making use of the $O(N)$ symmetry of the potential
$U(t,\vec{\varphi}\,)$ and the flow equations, which implies a
$\mathbb{Z}_{2}$ antisymmetry of $u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$,

$\displaystyle U(t,\sigma)=U(t,-\sigma)$ $\displaystyle\Longleftrightarrow$
$\displaystyle u(t,\sigma)=-u(t,-\sigma)\,.$ (136)

This reduces the spatial domain to the half-open interval
$\sigma\in[0,+\infty)$, but now we need an additional artificial boundary
condition at $\sigma=0$, see, e.g., Ref. [130]. In previous studies, the use
of the $O(N)$ symmetry was usually implemented right from the beginning by
replacing the variable $\vec{\varphi}$ by the $O(N)$ invariant
$\varrho=\tfrac{1}{2}\,\vec{\varphi}^{\,2}$, whose domain is already by
definition $[0,\infty)$.292929In any case, independent of the implementation
of the boundary condition itself, one should make use of symmetries of the
flow equations in numerical implementations. First of all, this leads to a
reduction of the number of computational grid points in spatial direction,
while keeping the spatial resolution fixed, which significantly speeds up the
calculations (independently of the specific numerical method for spatial
discretization). An additional consequence is the reduction of numerical
errors: It is highly unlikely that the numerical errors are symmetric in $x$,
if a symmetric interval around $x=\sigma=0$ is used. This might lead to an
artificial breaking of the $\mathbb{Z}_{2}$ (anti-)symmetry by unbalanced
numerical errors. Although these errors might be tiny and almost negligible
they can be easily circumvented by exploiting the symmetries. Using the
symmetries of a problem is a standard procedure in practical computations and
of particular importance in, e.g., numerical fluid dynamics and numerical
(general) relativity, see Refs. [222, 223, 224, 225]. In this case one has to
define

$\displaystyle
u(t,\varrho)\equiv\partial_{\varrho}U(t,\rho)=\tfrac{1}{\sigma}\,\partial_{\sigma}U(t,\sigma)=\tfrac{1}{\sigma}\,u(t,\sigma)\,,$
(137)

to obtain a flow equation for $u(t,\varrho)$ in a manifestly conservative
form, see Eqs. (99) and (100).

Before returning to the remaining boundary condition at $+\infty$, we first
consider the newly introduced artificial boundary condition at ${x=\sigma=0}$
or, correspondingly, at $\varrho=0$.

####  IV.4.1 Boundary condition at $\sigma=0$

At first sight it might be appealing to formulate the whole problem – the
conservation equation and the boundary condition at $\sigma=0$ – in the
variable $\varrho$. However, we believe that a formulation in $\sigma$ is more
suitable and easier to implement in our numerical setup.303030We do not claim
that it is impossible to formulate well-defined discrete boundary conditions
in $\varrho$ at $\varrho=0$, as can be seen for example in Refs. [26, 27, 31]
for the specific case of the large-$N$ limit of the $O(N)$ model. However, we
were not able to provide a suitable discretization of the boundary condition
at $\varrho=0$ in the implementation of the FV method for flow equations that
include diffusion via the radial $\sigma$ mode.  A key feature of (non-linear)
hyperbolic/parabolic conservation equations is that their weak solutions may
exhibit non-analyticities in the form of shock and rarefaction waves etc.,
which manifest themselves in the solution in cusps or discontinuities in
spatial direction during the time evolution. These effects can develop during
the time evolution even if the initial condition is smooth/analytic see, e.g.,
Refs. [192, 226, 154, 155, 78, 79, 137]. As demonstrated in Refs. [118, 142,
139, 26, 27, 141, 32] this also holds for FRG flow equations, where non-
analyticities are inherent properties of the effective IR potential
$U(t_{\mathrm{IR}},\sigma)$. These statements are also true for the point
$\sigma=0$, where $U(t,\sigma)$ and $u(t,\sigma)$ do not need to be analytic,
see Sub.Sec. V.4. Hence, there might be a scenario where the potential
$U(t,\sigma)$, although it is symmetric in $\sigma$, has a cusp at $\sigma=0$,
which would correspond to a jump in a weak solution for
$u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$ at $\sigma=0$. If formulated in
$\varrho$, any scenario (analytic or non-analytic at $\sigma=0$) merely
corresponds to some arbitrary value for
$u(t,\varrho)=\partial_{\varrho}U(t,\varrho)$ at $\varrho=0$, which seems to
be of great advantage, because one does not have to deal with possible
discontinuities in the conserved quantity $u$. Furthermore, the problematic
factors of $\frac{1}{\sigma}$ in the pion propagator and the advection flux
(109), which are diverging at $\sigma=0$, can be avoided when formulating the
flow equations in $\varrho$.

Nevertheless, a problem with the variable $\varrho$ becomes apparent when
turning to the discretized form of $u$ within the FV scheme presented in
Sub.Secs. IV.2 and IV.3: FV methods (and also other discretization schemes)
usually require ghost cells at the boundaries of the computational domain,
since the in- and out-flows for the $i^{\text{th}}$ cell are calculated from
the cell averages $\bar{u}$ of its neighboring cells, cf. Eq. (129). However,
initially these values are not specified for the cells at the boundaries of
the computational domain. Thus, artificial ghost cells must be introduced and
the numerical values for $\bar{u}$ in these ghost cells have to be implemented
by hand or reconstructed from the cells within the computational domain in
accordance with the boundary conditions [78, 79]. In the second-order
formulation of the one-dimensional KT scheme one needs two ghost cells at each
of the two spatial boundaries, cf. Eq. (129).

However, implementing ghost cells for $u(t,\varrho)$ at $\varrho=0$ is
conceptually difficult, because these ghost cells must be centered at negative
values for $\varrho$ outside the computational domain $[0,\infty)$, which by
definition do not exist due to the positivity of
$\varrho=\tfrac{1}{2}\,\sigma^{2}$. A priori, it is therefore not clear how
numerical values $\bar{u}(t,\varrho_{i})$ should be assigned to ghost cells at
negative $\varrho_{i}$, because symmetry arguments cannot be applied anymore.

Furthermore, it is also not a feasible option to move the ghost cells to
positive values of $\varrho_{i}$, such that the point $\varrho=0$ is no longer
part of the computational domain. Namely, having ghost cells centered at small
but positive $\varrho_{i}$ implies that one has to extrapolate the numerical
values $\bar{u}(t,\varrho_{i})$ to these ghost cells and to the point
$\varrho=0$ from the other ordinary cells of the computational domain.
However, the functional behavior of $u(t,\varrho)$ is unknown for small
$\varrho$ and is actually exactly what we want to calculate in the first place
by solving the PDEs. Thus, any extrapolation at small $\varrho$ can only be
considered an educated guess. It is especially dangerous, because the physical
point will be part of the extrapolated ghost cells if it is located at
$\varrho=0$, which is the case for all models in their symmetric phase [32],
irrespective of the dimensionality of space-time. Consequently, extrapolation
errors at the physical point have the potential to spoil the numerical values
of all $n$-point correlation functions, which are calculated at the physical
point via derivatives of $u$ and contain the physics of the model. Even if the
physical point is at finite non-zero $\varrho$ far away from the ghost cells
and the boundary at $\varrho=0$, any kind of extrapolation at small $\varrho$
leads to numerical errors, because the diffusive contributions of the radial
$\sigma$ mode will propagate this information from smaller to larger $\varrho$
and hence to the physical point. Similar problems in formulating appropriate
boundary conditions at $\varrho=0$ also exist in other discretization schemes
like finite-difference or finite-element methods.

There is one exception to this conclusion: In the large-$N$ limit of the
$O(N)$ model the flow equation for $u$ reduces to a pure advection equation.
Studying the characteristic velocities, which are given by $\partial
F/\partial u$, respectively, see Eq. (113), we find that these cannot change
their sign, and information (or the conserved quantity $u$) is always
propagated via advection in the direction of smaller $\varrho$ or $|\sigma|$.
In this scenario, ghost cells can be positioned at negative $\varrho_{i}$ and
the corresponding cell averages $\bar{u}_{i}$ in the ghost cells can take any
numerical value since information from the ghost cells is never propagated
back into the computational domain and cannot cause any errors, cf. Ref. [26,
27, 31]. Shifting the ghost cells into regions of positive $\varrho$ is still
not suitable for the reasons already discussed above.

In order to avoid all these difficulties when formulating the problem in the
variable $\varrho$, we suggest a formulation in $\sigma$ and an implementation
of the boundary condition at $\sigma=0$. The key argument for using $\sigma$
instead of $\varrho$ is that positioning ghost cells at negative $\sigma$
poses no problem at all, since negative $\sigma$ exist in the first place.
Furthermore, it is clear how the cell averages $\bar{u}(t,\sigma_{i})$ in the
ghost cells have to be chosen: Using the antisymmetry (136), one merely has to
mirror the last physical cells of the computational domain at $\sigma=0$ to
the ghost cells (including a flip in sign). The only issue that requires
careful consideration is the choice of the position of the first physical cell
$x_{0}$ next to $\sigma=0$: The flux term of our PDE contains factors
$\frac{1}{\sigma}$ via the pion propagators, which diverge if evaluated at
$\sigma=0$. Therefore, we must avoid evaluating the fluxes $F[t,x,u(t,x)]$ at
$x=\sigma=0$. However, inspecting the KT scheme, we find that the fluxes as
well as the Jacobian $\frac{\partial F[u]}{\partial u}$ must only be evaluated
at the cell boundaries $x_{j\pm\frac{1}{2}}$, cf. Eqs. (124) and (131).
Consequently, the natural choice for the position of the cell center $x_{0}$
of the first physical cell in the computational domain is at $x=\sigma=0$,
such that the in- and out-fluxes of this cell are evaluated at
$x_{\pm\frac{1}{2}}$, which is not problematic. Incidentally, this
automatically cures the problem of the possibility of non-analyticities in
$u(t,\sigma)$ at $\sigma=0$: Even if $u(t,\sigma)$ is discontinuous at
$\sigma=0$ we do not run into problems, because all numerical calculations are
performed on the level of cell averages $\bar{u}(t,\sigma_{i})$. The cell
average of an antisymmetric function in a cell that is centered at $\sigma=0$
must always vanish identically, independent of all other properties of the
function, see also Ref. [130, 131].

In summary, we switch from the open computational interval $(-\infty,+\infty)$
to the half-open computational interval $[0,+\infty)$ by means of the
$\mathbb{Z}_{2}$ (anti-)symmetry using $\bar{u}_{-2}(t)=-\bar{u}_{2}(t)$,
$\bar{u}_{-1}(t)=-\bar{u}_{1}(t)$ for the cell averages in the ghost cells
left of $x_{0}=0$ and $\bar{u}_{0}(t)=0$ for the cell average in the cell at
$x_{0}$. This effectively corresponds to reflective boundary conditions
frequently imposed in numerical fluid dynamics [78, 79], cf. Fig. 3.

Figure 3: Second-order accurate FV implementation of the spatial boundary
condition for $u(t,x)$ or $\bar{u}_{i}(t)$, respectively, at $x=0$. We use the
fact that $u(t,x)$ is an odd function in $x$ by positioning the first
computational cell $x_{0}$ at $x=0$, such that the cell average is exactly
zero, $\bar{u}_{0}=0$, which is true for $u(t,x)$ which are analytic (blue-
dashed) as well as non-analytic (green-solid) at $x=0$. The ghost-cell
averages can thus be fixed by setting $\bar{u}_{-2}=-\bar{u}_{2}$ and
$\bar{u}_{-1}=-\bar{u}_{1}$. Corresponding cell averages $\bar{u}_{i}$ are
depicted as horizontal bars (magenta-dashed and yellow-solid). This boundary
condition can be generalized to lower- and higher-order accurate FV schemes as
well as finite-difference or finite-element schemes.

####  IV.4.2 Boundary condition at $\sigma\rightarrow\infty$

Now we return to the boundary condition at $\sigma\rightarrow+\infty$. Without
loss of generality we discuss the interval ${\sigma\in[0,+\infty)}$ since the
situation in $\sigma\in(-\infty,0]$ follows from $\mathbb{Z}_{2}$ antisymmetry
of $u(t,\sigma)$.

We have already argued that there are no real boundary conditions at spatial
infinity on a non-compact domain. The behavior of $u$ at
$\sigma\rightarrow\infty$ is rather given by the asymptotics of the FRG
equation, which makes the PDE an pure initial-value problem. The boundary
condition at spatial infinity is actually fixed implicitly: As long as the
initial potential $U(t=0,\sigma)$ is bounded from below and grows faster than
$\sigma^{2}$ for $\sigma\rightarrow\infty$ both pion and sigma propagator tend
to zero for sufficiently large $\sigma$, such that the right-hand side of the
PDE (101) vanishes during the entire FRG flow. In the fluid-dynamical picture
this corresponds to vanishing advection and diffusion fluxes (109) and (110)
at $\sigma\rightarrow\infty$, which is a zero-influx boundary condition for
$u(t,\sigma)$. The derivative of the effective potential $u(t,\sigma)$ is
therefore fixed to its initial value $u(t=0,\sigma)$ at
$\sigma\rightarrow\infty$.

The limiting case, when the asymptotic behavior of the initial UV potential is
quadratic,

$\displaystyle\lim\limits_{\sigma\rightarrow\infty}U(t=0,\sigma)\sim\sigma^{2}\,,$
$\displaystyle\lim\limits_{\sigma\rightarrow\infty}u(t=0,\sigma)\sim\sigma\,,$
(138)

is more delicate. In this case, the advection and diffusion fluxes (109) and
(110) do not vanish for $\sigma\rightarrow\infty$ for all RG times. However,
for small RG times $t\approx 0$, the fluxes are actually independent of
$\sigma$ at large $\sigma$ due to the constant asymptotic slope of the initial
condition $u(t=0,\sigma)$. This in turn implies that the in- and out-flux for
all volume cells at large $\sigma$ only depend on $t$ and must cancel exactly,
such that the net flux of these cells vanishes. Therefore, also in this
scenario $u(t,\sigma)$ is fixed to its initial condition at
$\sigma\rightarrow\infty$ not only for small $t$, but rather for all RG times
$t$. For late RG times $t\rightarrow\infty$, the advection and diffusion
fluxes (109) and (110) vanish anyhow, due to the derivatives of the regulator
shape functions in the numerators, i.e.,
$\partial_{t}r(t)=-\Lambda\,\mathrm{e}^{-t}$. In the language of fluid
dynamics, initial conditions with quadratic asymptotics can therefore be
interpreted as boundary conditions with time-dependent but spatially constant
in-flux, cf. Examples 7 and 9 in Ref. [192].

However, both cases cannot be implemented directly on a finite computational
domain and we basically have two options:

  1. 1.

We could try to map the interval $[0,\infty)$ to a compact interval $[0,1]$
via a suitable map $\sigma\mapsto x(\sigma)$. This also includes a suitable
mapping of $u\mapsto v(u)$ to keep the values for the conserved quantity
finite on $[0,1]$. This option has the advantage that the correct asymptotic
behavior $u(t,\sigma)$ can be implemented as boundary conditions for $v(t,x)$
at $x=1$. However, the same question then arises as before in the discussion
of an appropriate choice of ghost cells for negative values of $\varrho$: It
is highly non-trivial how the cell averages $\bar{v}_{i}$ should be fixed for
ghost cells which no longer belong to the physical values of $x$ within the
interval $[0,1]$. Additionally, the two mappings would introduce at least two
new numerical functional-mapping parameters. A suitable choice of these
parameters is not obvious. Still, these mappings would have to ensure dense
grids and high resolution around the physical point and low resolution at
large field values All this is extremely hard to achieve. Therefore, we
propose and favor another option.

  2. 2.

The second option, which is our preferred choice, is to split the physical
domain $[0,\infty)$ into a compact domain $[0,\sigma_{\mathrm{max}}]$ and a
non-compact domain $[\sigma_{\mathrm{max}},\infty)$. Here,
$\sigma_{\mathrm{max}}$ should be chosen much larger than the physical scales
of the problem and the position of the physical point, see, e.g., Refs. [26,
141, 32, 130, 180, 129]. We will provide explicit tests for an appropriate
choice of $\sigma_{\mathrm{max}}$ later on in Sec. V. For the compact domain
$[0,\sigma_{\mathrm{max}}]$, we keep a direct identification of the field
$\sigma$ and the computational spatial variable $x$, thus $x=\sigma$. For
higher-dimensional models this might be replaced by a linear map of $\sigma$
to a dimensionless spatial variable $x$ via appropriate rescaling with some
characteristic dimensionful quantity, e.g., the UV cutoff $\Lambda$ or a non-
vanishing condensate. In any case, this allows for a direct implementation of
the boundary condition at $\sigma=0$.

In the compact domain $[0,\sigma_{\mathrm{max}}]$, we have to ensure a high
spatial resolution via a sufficiently large number of cells, in order to
capture all aspects of the dynamics around the physical point. Explicit tests
to find an appropriate spatial resolution are also presented in Sec. V.

For the non-compact domain $[\sigma_{\mathrm{max}},\infty)$, instead of using
a discretization scheme like the FV method, we suggest an expansion or
approximation of $u(t,\sigma)$ via polynomials or complete sets of functions
with $t$ dependent expansion coefficients, which account for the asymptotic
behavior of the initial condition $u(t=0,\sigma)$ for large $\sigma$. As
discussed before, it is expected that for large $\sigma$ the deviations of
$u(t,\sigma)$ from the initial condition $u(t=0,\sigma)$ are small during the
FRG flow, such that a finite amount of expansion coefficients should be
satisfactory to capture this minimal dynamics.

At the point $\sigma_{\mathrm{max}}$, the ghost cells for the FV method in
$[0,\sigma_{\mathrm{max}}]$ can therefore be fixed via the values
$u(t,\sigma)$ from the asymptotic expansion in the non-compact interval
$[\sigma_{\mathrm{max}},\infty)$.

Interestingly, our numerical tests showed that, as long as
$\sigma_{\mathrm{max}}$ is chosen sufficiently large, the fluxes at
$\sigma_{\mathrm{max}}$ are already negligibly small. As a consequence, the
deviation of $u(t,\sigma)$ from the initial condition in the non-compact
interval $[\sigma_{\mathrm{max}},\infty)$ is extremely small and can be
ignored. In this case, the computational boundary conditions for the ghost
cells at $\sigma_{\mathrm{max}}$ can be fixed via an extrapolation using the
asymptotics of the initial condition. For extremely high spatial resolution,
hence rather small $\Delta x$, even a simple linear extrapolation might be
sufficient.

On the other hand, choosing $\sigma_{\mathrm{max}}$ rather large while keeping
a high spatial resolution in the compact computational domain
$[0,\sigma_{\mathrm{max}}]$ requires a large number of cells. However, this
slows down the computations drastically. For problems where this issue becomes
relevant, we suggest to further divide the compact domain
$[0,\sigma_{\mathrm{max}}]$ into several smaller subdomains. In each of these
subdomains one can implement the FV method with different spatial resolution
$\Delta x$ for each domain. This ensures high resolution at small $\sigma$
next to the physical point and also allows to truncate the spatial interval at
large $\sigma_{\mathrm{max}}$, while keeping a decent and manageable total
number of cells [26, 27]. An alternative approach would be switching from
equally sized volume cells on a uniform grid to a non-uniform (potentially
even moving/time dependent) grid, see, e.g., Ref. [197]. However, in the
context of FRG flow equations this might be relevant for models with multiple
condensate directions, see, e.g., Refs. [181, 182, 183, 184]. In our test
cases for the zero-dimensional $O(N)$ model the subdivision of the compact
interval $[0,\sigma_{\mathrm{max}}]$ or a formulation on non-uniform grids
was, however, not necessary.

We close the discussion of spatial boundary conditions with some final
remarks: We have discussed that the correct implementation of spatial boundary
conditions is of utmost importance. Otherwise there is no guarantee for the
solution of a PDE and especially our FRG flow equations to converge against
the correct result. Additionally, we remark that our specific choice and
implementation might not be the best option at hand for arbitrary (higher-
dimensional) models and arbitrary initial conditions within the FRG framework.
In the current context of the zero-dimensional $O(N)$ model initial conditions
without a proper large-$|\vec{\phi}\,|$ asymptotics, e.g.,
$[2-\sin(\vec{\phi}^{\,2})]\,\vec{\phi}^{\,2}$ or even worse
$[2-|\sin(\vec{\phi}^{\,2})|]\,\vec{\phi}^{\,2}$, and/or periodic potentials
could be a very interesting topics for further research.

##  V Zero-dimensional field theory as testing ground for FRG

After our general discussion of the theoretical basis for the solution of FRG
flow equations, we shall discuss concrete applications in the following
subsections. To this end, we study the RG flow of various zero-dimensional
field theories which differ by distinct initial conditions. Our choices for
the initial conditions range from smooth potentials to extreme choices
featuring non-analyticities. Note that such extreme choices are not only
relevant to demonstrate the numerical performance and stability of our
implementation but also for phenomenological reasons. In fact, in higher
dimensions we expect non-analytic behavior to build up, e.g., in the IR limit,
as associated with spontaneous symmetry breaking and the emergence of
convexity of the effective action.

###  V.1 Test case I: Non-analytic initial condition

Consider the following initial UV potential,

$\displaystyle
U(\vec{\varphi}\,)=\begin{cases}-\tfrac{1}{2}\,\vec{\varphi}^{\,2}\,,&\text{if}\quad\varphi\leq
2\,,\vphantom{\bigg{(}\bigg{)}}\\\ -2\,,&\text{if}\quad 2<\varphi\leq
3\,,\vphantom{\bigg{(}\bigg{)}}\\\
+\tfrac{1}{2}\,(\vec{\varphi}^{\,2}-13)\,,&\text{if}\quad
3<\varphi\,,\vphantom{\bigg{(}\bigg{)}}\end{cases}$ (139)

where $\varphi=|\vec{\varphi}\,|$, cf. Fig. 4.

Figure 4: The plot shows the UV potential $U(\sigma)$ (red-dashed) and its
first derivative $u(\sigma)=\partial_{\sigma}U(\sigma)$ (blue, solid) of test
case Eq. (139) evaluated on the constant background field configuration.

The test case is designed this way for the following reasons:

  1. 1.

All parameters of the potential $U(\vec{\varphi})$ are by default
dimensionless and chosen to be approximately of order one, such that scales
can easily be compared with each other.

  2. 2.

The UV potential $U(\sigma)$ has non-analytical points at $\sigma=2$ and
$\sigma=3$, which correspond to discontinuities in its derivative $u(\sigma)$.
In the fluid-dynamical context such discontinuities present a Riemann problem
and result in rarefaction waves. In QFT and thermodynamics such
discontinuities can be associated with phase transitions, see App. B. The non-
analytical behavior of this potential makes commonly used techniques like the
FRG Taylor expansion inapplicable. Furthermore, naive discretizations that
rely on smoothness are doomed to fail. One has to choose numerical schemes
that can handle this numerically challenging dynamics.

  3. 3.

The potential is initialized in the symmetry-broken phase, with infinitely
many degenerate minima at $\sigma\in(2,\,3]$. Furthermore, the UV potential is
neither convex nor smooth. However, in the IR the $O(N)$ symmetry has to be
restored and there must only be one minimum at $\sigma=0$, due to the Coleman-
Mermin-Wagner-Hohenberg theorem. Furthermore, for $t\rightarrow\infty$ the
potential has to be convex and smooth, see App. B. This non-trivial transition
and complicated dynamics from the UV to the IR is a numerically challenging
test.

  4. 4.

Furthermore, we choose a potential which is asymptotically quadratic in
$\sigma$. This is to ensure that the large-$\sigma$ boundary condition for
$u(t,\sigma)$ is fully under control and errors can be excluded, see Sub.Sec.
IV.4. This allows for a high-precision analysis of all other sources of
numerical errors.

The reference values for the exact IR 1PI vertex functions $\Gamma^{(2n)}$ of
the $O(N)$ model are calculated numerically from the UV potential (139) via
the integral (77) using Eqs. (73) –(75). They are listed in Table 1 for
selected $N$.

Table 1: The table lists the (up to numerical-integration errors) exact results for the $\Gamma^{(2n)}$ of the $O(N)$ model with the initial UV potential (139) for selected $N$. They are obtained by a high-precision one-dimensional numerical integration of the expectation values $\langle(\vec{\phi}^{\,2})^{n}\rangle$ using Mathematica’s numerical integration routine NIntegrate [138] with a PrecisionGoal and AccuracyGoal of $10$. Here, we present the first ten digits only. $N$ | $\Gamma^{(2)}$ | $\Gamma^{(4)}$ | $\Gamma^{(6)}$
---|---|---|---
$1$ | $0.1768130358$ | $0.0520549107$ | $\phantom{-}0.0865733100$
$3$ | $0.3973542395$ | $0.1408641479$ | $\phantom{-}0.2249964322$
$10$ | $0.8451440328$ | $0.1519326337$ | $-0.0691341063$

####  V.1.1 General discussion of the FRG flow –
advection and diffusion

We start our analysis with a general discussion of the FRG flow with initial
condition (139).

To this end, we fix $O(N=3)$ to include both diffusive and advective
contributions via the radial $\sigma$ mode and two pions. For $N=3$ the number
of pions is reasonably small and the (diffusive) effects of the $\sigma$ mode
remain visible. Furthermore, we choose $[0,x_{\mathrm{max}}=10]$ as the
spatial computational domain with $800$ volume cells and use the KT scheme
from Sub.Sec. IV.3 for spatial discretization. The initial cell averages
$\bar{u}_{i}(t=0)$ are computed by exact averaging313131Using the exact
averages for $\bar{u}_{i}(t=0)$ has proven advantageous in terms of achievable
numerical precision in the IR compared to taking the mid-point values of the
exact derivative of
${\bar{u}_{i}(t=0)=\partial_{\sigma}U(\sigma)|_{\sigma=\sigma_{i}}}$ when
considering non-analytic initial conditions like Eqs. (139) or (150).

$\displaystyle\bar{u}_{i}(t=0)=\tfrac{1}{\Delta\sigma}\,\big{[}U\big{(}\sigma_{i+\frac{1}{2}}\big{)}-U\big{(}\sigma_{i-\frac{1}{2}}\big{)}\big{]},$
(140)

using the UV initial condition (139). We use linear extrapolation as spatial
boundary condition at $x_{\mathrm{max}}$. The UV cutoff is set to
$\Lambda=10^{6}$ at $t=0$. The integration is performed with Mathematica’s ODE
solver NDSolve [138] with a PrecisionGoal and AccuracyGoal of 10 and stopped
in the IR at $r(t_{\mathrm{IR}}=60)\approx 10^{-20}$ using the exponential
regulator shape function (8). We find that these parameters suffice to produce
decent results, as discussed in the following subsubsections. There, we
quantitatively analyze sources and severity of possible errors.

We first provide qualitative results of our numerical methods in Figs. 5, 6,
and 7, where we show the FRG flow of the effective potential $U(t,\sigma)$ and
its derivative $u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$ from the UV (blue)
to the IR (red).

Figure 5: The FRG flow of the effective potential $U(t,\sigma)$ (upper panel)
and its derivative $u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$ (lower panel)
for the zero-dimensional $O(3)$ model with initial condition Eq. (139)
evaluated at $t=0,\,2,\,4,\,\ldots,\,60$ (integer values for $t$ were chosen
for convenience and readability). The blue curves correspond to the UV and the
red curves to the IR. We used the exponential regulator Eq. (8) with UV cutoff
$\Lambda=10^{6}$. For the sake of readability, the plot does not show the
region $x=5$ to $x=10$, because the tiny differences between $u(t,\sigma)$ and
$u(0,\sigma)$ are not visible in this region and vanish for large $x=\sigma$
anyhow.

In the beginning, i.e., in the UV, the system is in the broken phase. This
changes only marginally until $t\approx 7$, which indicates that the UV cutoff
is chosen sufficiently large. When $r(t)$ reaches the scales of the model at
$t\gtrsim 11$ most of the dynamics takes place (symmetry restoration) and
$u(t,\sigma)$ changes rapidly and drastically until it freezes out in the IR.
In the IR the system is in the restored phase. After $t\approx 26$ the
potential does not change anymore, which indicates that one has reached a
sufficiently small IR scale to stop the integration. We render this statement
more precise in the following subsubsections. Note that the evolution in $t$
is logarithmic and corresponds to a change in scale over $25$ orders of
magnitude. At first glance this range sounds excessive, but its necessity is
explained in detail below.

During the FRG evolution one observes that diffusive processes smear out the
discontinuities at the non-analytic points at $\sigma=2$ and $\sigma=3$. We
also find that the diffusion acts in both directions – towards larger and
smaller values of $\sigma$ – as expected from our discussion in Sub.Sec. IV.1.
Nevertheless, the diffusion effects do not reach the computational boundary,
which is outside the plot range at $\sigma_{\mathrm{max}}=10$. This suggests
that $\sigma_{\mathrm{max}}=10$ is sufficiently large. Additionally, we
observe a propagation of the conserved quantity $u$ towards smaller values of
$\sigma$ via advection. Close to the initial discontinuities these advective
processes can be interpreted as rarefaction waves. In a more pictorial
language, the advection and diffusion “fill up the pit” in $u(t,\sigma)$ at
small values of $\sigma$ by moving more and more of the quantity $u$ from
larger values of $\sigma$ to smaller $\sigma$ (via advection and diffusion) as
well as from small negative $\sigma$ to small positive sigma (via diffusion).
Eventually the symmetry is restored and $u(t,\sigma)$ is smoothed towards the
IR by diffusion. Furthermore, the conserved quantity $u$ does not “pile up” at
$\sigma=0$ after symmetry restoration, because at negative $\sigma$ exactly
the opposite dynamics happens, due to the $\mathbb{Z}_{2}$ antisymmetry of
$u(t,\sigma)$, which is encoded in the boundary condition at $\sigma=0$, see
Sub.Sec. IV.4. The differences between advective and diffusive contributions
become apparent when comparing the same system for different $N$, see below.

Figure 6: The RG flow of the derivative of the effective potential
$u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$ for the zero-dimensional $O(3)$
model with initial condition Eq. (139). This is a 3-dimensional rendering of
the flow displayed Fig. 5 (lower panel).

From a numerical perspective, the KT scheme is able to handle the highly non-
linear dynamics, including the non-analyticities in $u(t,\sigma)$, without any
spurious oscillatory behavior (under-/over-shooting) and allows for a stable
$t$ integration to extremely small IR scales.

For the sake of completeness and illustrative purposes, we also provide the RG
flow of the integral of $u(t,\sigma)$, i.e., the effective potential
$U(t,\sigma)$, in Figs. 5 and 7. Here, the integration constant was set to
zero323232$U(t,0)=0$ is dictated by our choice of normalization for the zero-
point function(s), see Eq. (4).  and the integration was performed via Riemann
summation333333At this point we should mention that numerical results for
$U(t,\sigma)$ via Riemann summation should be treated with great caution:
Numerical errors in the cell averages $\bar{u}(t,x_{j})$ which arise from the
numerical FRG flow can accumulate during integration (here summation) along
$\sigma=x$. More precise quadrature methods should be used if precise,
quantitative values for $U(t,\sigma)$ are needed. This will be discussed
elsewhere [153].  of the discrete cell averages,

$\displaystyle U(t,x_{i})=\Delta
x\sum_{j=0}^{i}\frac{\bar{u}(t,x_{j})}{(1+\delta_{j0}+\delta_{ji})}\,.$ (141)

Figure 7 perfectly illustrates the restoration of the $O(3)$ symmetry of the
potential $U(t,\sigma)$ during the FRG flow via “vaporization” of the
infinitely many minima. Nevertheless, we find that it is hardly possible to
intuitively understand the contributions of the radial $\sigma$ mode and the
pions to the FRG flow on the level of $U(t,\sigma)$ only. This complements the
discussion in Sub.Sec. IV.1 and lends support to our claim that the fluid-
dynamical interpretation of the FRG flow in terms of $u(t,\sigma)$ is superior
to the canonical formulation in terms of $U(t,\sigma)$ commonly used in the
FRG literature.

Figure 7: The RG flow of the effective potential $U(t,\sigma)$ for the zero-
dimensional $O(3)$ model with initial condition Eq. (139). This is a
3-dimensional rendering of the flow displayed Fig. 5 (upper panel).

Before discussing quantitative results and sources of (numerical) errors in RG
flows, we briefly return to the interpretation of the radial $\sigma$ mode as
diffusion and the interpretation of the pions as advection in the FRG flow
along the field space direction. To this end, we discuss FRG flows for the
same initial UV potential (139) as before, but for different $N$. This
corresponds to a different number of pions in the flow and different advection
velocities (113). All other parameters remain unchanged. In addition to the
$N=3$ case in Fig. 5, we provide the RG flows of $u(t,\sigma)$ for
$N=1,10,100$ in Fig. 8.

Figure 8: The FRG flow of the derivative of the effective potential
$u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$ for the zero-dimensional $O(N)$
model for $N=1,\,10,\,100$ with initial condition Eq. (139). All parameters
are identical to those in Fig. 5.

The figure again demonstrates on a qualitative level that the sigma mode
enters the FRG as diffusion, while pions enter as advection: Increasing the
number of pions the problem becomes more advection-driven exhibiting more
pronounced waves and faster propagation. This can be seen by comparing the
plots at equal RG times. For $N=1$, the problem reduces to the pure diffusion
equation (38), where the dynamics is slowest and the diffusion propagates the
fluid from negative $\sigma$ to small positive $\sigma$ close to $\sigma=0$.
Furthermore, one observes stronger smearing of the discontinuities at
$\sigma=2$ and $\sigma=3$. The $N=100$ case is extremely advection-
dominated343434We will discuss the qualitative and quantitative differences
between RG flows at large and infinite $N$ in part III of this series of
publications [31]., similar to the situation in large-$N$ limit of the $O(N)$
model, where the $\sigma$ mode can be completely ignored. In the fluid-
dynamical language, this corresponds to a complete suppression of diffusion,
which is clearly observed from the fast propagation of the rarefaction waves
and almost negligible smoothing of the discontinuities at $\sigma=2$ and
$\sigma=3$.

####  V.1.2 Tests of the spatial resolution $\Delta x$

This subsubsection is dedicated to quantifying numerical errors in the FRG
flow arising from the finite spatial resolution $\Delta x$ of the cells in the
KT scheme. Any spatial discretization comes with a discretization error. The
KT scheme, which is used throughout this paper, is of second-order accuracy
$\Delta x$. Therefore, the numerical errors arising from the spatial
discretization should scale with $\Delta x^{2}$ when $\Delta x$ is decreased.

As test values (observables) we use the modulus of the relative errors of the
1PI $n$-point vertex functions $\Gamma^{(2n)}$ for $n=1,2,3$,

$\displaystyle\bigg{|}\frac{\Gamma^{(2n)}_{\mathrm{KT}}}{\Gamma^{(2n)}}-1\bigg{|}\,,$
(142)

where $\Gamma^{(2n)}_{\mathrm{KT}}$ is calculated from the FRG (via the KT
scheme) and $\Gamma^{(2n)}$ from the (functional) integral, see Table 1. In
order to determine how much of the relative numerical error (142) is directly
related to the spatial discretization, we have to optimize all other
parameters of the computation in order to reduce other sources of errors. We
basically choose the same parameter set and UV initial condition (139) which
was also used for the qualitative analysis in the previous subsection and only
vary the number of cells $n$ to change the resolution $\Delta x$. We keep
$\Lambda=10^{6}$, $x_{\mathrm{max}}=10$ and $t_{\mathrm{IR}}=60$, which turn
out to be decent choices as will become clear in the following. Without loss
of generality we keep $N=3$.

Before we embark on our discussion, we remark that the spatial-discretization
error enters the relative errors (142) in a twofold way:

  1. 1.

There is the discretization error which comes from the KT scheme during the
FRG time integration. This error shows up directly in the IR values
$u(t_{\mathrm{IR}},x_{i})$ and should scale according to $\Delta x^{2}$ for
the chosen second-order KT scheme.

  2. 2.

There is a discretization error which is related to the extraction of the 1PI
$n$-point vertex functions from the discrete $\bar{u}(t_{\mathrm{IR}},x_{i})$.
They have to be calculated at the physical point (the minimum at $x=\sigma=0$)
via numerical differentiation, which also comes with a discretization error.
The latter is also related to the spatial resolution $\Delta x$. In general
(especially in higher-dimensional models) it is not clear whether these
numerical derivatives at the physical point are always well-defined. We argued
before that $u(t,\sigma)$ might exhibit non-analytical behavior also at the
physical point in the IR, cf. Refs. [26, 27, 137, 32], such that a naive
numerical differentiation is not allowed in general. This will be discussed
and analyzed in detail elsewhere [153]. In the special case of zero-
dimensional QFTs, we prove in App. B that the IR effective action and the IR
potential $U(t\rightarrow\infty,\vec{\varphi}\,)$ are smooth, which also
translates to $u(t\rightarrow\infty,\sigma)$, such that numerical
differentiation (e.g., via finite-difference approximations) is well-defined.

However, even though finite-difference formulas are reliable approximations
for derivatives of a smooth function and have a well-defined truncation-error
scaling in powers of $\Delta x$, there remains a well-known subtlety: While
decreasing the resolution $\Delta x$, one eventually will reach a point where
the error starts increasing again contrary to the formal truncation-error
scaling. This is related to the ill-conditioned nature of finite-difference
formulas and to explicit rounding errors in floating-point arithmetic, which
increase the error of the numerical derivative below a certain $\Delta x$,
see, e.g., Chap. 5.7 of Ref. [43]. Related spurious cancellations occur if the
discrete data of a smooth function include some sort of noise. In our case
this “noise” is directly related to the spatial-discretization errors from the
KT scheme and the errors from RG time integration using numerical ODE solvers.
These errors might be tiny, but can easily inflate the errors of the numerical
derivatives, especially for higher-order derivatives.

In conclusion, while decreasing $\Delta x$ we expect that long before the
relative errors from the KT scheme start increasing again (because the KT
scheme begins operating close to machine precision or because the error of the
time stepping becomes relevant) the relative errors (142) start increasing due
to the numerical differentiation of slightly “noisy data”. This phenomenon is
expected to set in at larger $\Delta x$ for approximations for higher-order
derivatives and long before the formal error scaling of the KT scheme is no
longer valid.

Our explicit results for the scaling of the relative errors with decreasing
spatial resolution are presented in Fig. 9, where we show the relative errors
(142) for the two-, four-, and six-point functions in a double-logarithmic
plot over $\Delta x$.

Figure 9: The relative error of the numerical results (blue dots) from the KT
scheme for the 1PI $n$-point vertex functions $\Gamma^{(2n)}$ for $n=1,2,3$ as
a function of $\Delta x$ with initial potential (139). The numerical
derivatives at $\sigma=0$ of $u(t_{\mathrm{IR}}=60,\sigma)$ were calculated
via the second-order accurate central schemes (151), (154), and (157). The
plot was produced with $x_{\mathrm{max}}=10$, but could have been calculated
for any sufficiently large $x_{\mathrm{max}}$. We used the exponential
regulator (8) with UV cutoff $\Lambda=10^{6}$. The yellow straight lines are
for optical guidance.

For $\Gamma^{(2)}$ and $\Gamma^{(4)}$ we find that the corresponding relative
errors scale with $\Delta x^{2}$ (or even slightly better) as expected from
the error scaling of the KT scheme as well as the error scaling of the finite-
difference stencils (151) and (154). We observe two groups of points for
$\Gamma^{(2)}$ (upper panel of Fig. 9), which are related to discretization
errors of the discontinuous initial condition (139) at $x=2$ and $x=3$. The
error scaling of $0.02\,\Delta x^{2}$ for $\Gamma^{(2)}$ is a conservative
estimate for the observed errors, which are actually lower for $\Delta
x>0.005$. For $\Delta x<0.005$ we observe deviations from the conservative
estimate for the error scaling of $\Gamma^{(2)}$ related to other error
sources. In the middle panel of Fig. 9, we clearly see that there is an
optimal minimal $\Delta x\approx 0.02$ where the correct formal scaling of the
numerical derivative breaks down and the relative error of $\Gamma^{(4)}$
increases again for smaller $\Delta x$. We can be sure that this breakdown of
the error scaling is related to the numerical differentiation and not the KT
scheme because we observe scaling with at least $\Delta x^{2}$ for
$\Gamma^{(2)}$ in the upper panel of Fig. 9 well below $\Delta x\approx 0.02$.
This is expected for lower-order numerical derivatives. Furthermore, we find
that for $\Gamma^{(6)}$ (lower panel of Fig. 9) the order of the numerical
derivative is already too large, such that the theoretical error scaling of
the KT scheme cannot be seen at all and is completely obscured by the errors
from the numerical differentiation of $\bar{u}(t_{\mathrm{IR}},x_{i})$.

We conclude that the KT scheme is perfectly suited for the spatial
discretization of the RG flow equation for $u(t,\sigma)$ and shows correct
scaling with decreasing spatial resolution $\Delta x$. This is also confirmed
by tests with different initial conditions, see below.

In addition, we actually found that a more severe problem is the correct
extraction of physical observables from the IR values
$\bar{u}(t_{\mathrm{IR}},x_{i})$, which are usually related to derivatives of
$u(t_{\mathrm{IR}},\sigma)$. We further conclude that this problem is expected
to be worse in higher dimensions, were the IR potential is no longer
guaranteed to be smooth. We therefore suggest to search for better ways of
calculating those derivatives as well as for careful analysis tools for
numerically calculated 1PI $n$-point vertex functions in the vicinity of non-
analyticities in general. However, this is beyond the scope of the present
work.

We remark that our numerical findings indicate that – independent of the
specific numerical discretization scheme – the number of grid points or
expansion coefficients etc. might have been chosen too small in previous
studies to obtain a decent resolution. However, other works, cf. Refs. [227,
129, 130, 131, 180, 137], which also discuss the limitations of their
numerical schemes in detail, have used a rather large number of discretization
points – in some cases to compensate the demand for continuity of the specific
scheme.

In the next subsections we will mostly use a spatial resolution of

$\displaystyle\Delta x=\frac{x_{\mathrm{max}}}{n-1}\simeq 0.025\,,$ (143)

where we can trust the results for the two- and four-point functions. The
relative errors for the six-point function will only be plotted for the sake
of completeness, but cannot be included in any reasonable quantitative
analysis of other sources of (numerical) errors in RG flow equations although
they are still at an acceptably small level.

####  V.1.3 Tests of the size of the computational domain

In this subsection, we discuss the influence of the size of the computational
domain $[0,\sigma_{\mathrm{max}}]$ on the relative errors of the IR
observables (142). As discussed in Sub.Sec. IV.4, we expect that, if the
spatial boundary conditions are not implemented with great caution and the
computational domain is too small, one cannot trust the results from the
numerical integration of the RG flow. If the computational domain is too
small, we expect large errors, because the boundary conditions at
$\sigma_{\mathrm{max}}$ are no longer valid due to wrong extrapolation to the
ghost cells and consequently wrongly estimated in-flux.

In the case with UV initial condition (139), the boundary condition at
$\sigma_{\mathrm{max}}$ is implemented as a linear extrapolation of
$u(t,\sigma)$ to the two ghost cells of the KT scheme to mimic the asymptotic
behavior of $u(t,\sigma)$. As long as $\sigma_{\mathrm{max}}$ is sufficiently
large, we expect only tiny deviations of $u(t,\sigma)$ from its initial UV
value $u(t_{\mathrm{UV}}=0,\sigma)$ next to $\sigma_{\mathrm{max}}$. However,
if $\sigma_{\mathrm{max}}$ is too small and approaches to the model scales, we
expect the diffusion effects to reach the boundary of the computational
domain, such that a linear extrapolation is no longer a good approximation in
order to determine the spatial boundary condition.

To this end, we test the scaling of the relative errors (142) with decreasing
computational domain size $x_{\mathrm{max}}=\sigma_{\mathrm{max}}$ for $N=1$
(purely diffusive) and $N=3$. The results and (numerical) parameters are shown
in Figs. 10 and 11.

Figure 10: The relative error for $\Gamma^{(2m)}$ for $m=1,2,3$ for the UV
potential (139) of the $O(1)$ model as a function of $x_{\mathrm{max}}$, while
keeping the cell size constant, $\Delta x=0.025$. $\Gamma^{(2m)}$ are computed
from the discrete values of the derivative of the IR potential
$u(t_{\mathrm{IR}}=60,\sigma)$ via the second-order accurate central finite-
difference stencils (151), (154), and (157) at $\sigma=0$. We use the
exponential regulator (8) with UV cutoff $\Lambda=10^{6}$. The yellow straight
line is for optical guidance. Figure 11: The relative error for
$\Gamma^{(2m)}$ for $m=1,2,3$ for the UV potential (139) of the $O(3)$ model
as a function of $x_{\mathrm{max}}$. All parameters are identical to Fig. 10.

In both cases we find that the relative errors are independent of
$\sigma_{\mathrm{max}}$ for sufficiently large $\sigma_{\mathrm{max}}$.
However, if the spatial cutoff $\sigma_{\mathrm{max}}$ is approaching the
model scales (here the discontinuity in $u(t_{\mathrm{UV}}=0,\sigma)$ at
$\sigma=3$, see Fig. 4) the relative errors for $\Gamma^{(2)}$ and
$\Gamma^{(4)}$ start rising exponentially.

Contrary to our expectations, the results for $N=1$ and $N=3$ are very similar
and the exponential rise of the relative errors sets in at a similar
$\sigma_{\mathrm{max}}$. We expected that for the purely diffusive scenario
with $N=1$, the diffusion effects arising from the large gradients at
$\sigma=3$ might have more time to reach and influence the shape of
$u(t,\sigma)$ at larger values of $\sigma$, which does not seem to be the
case. Our employed monitors for numerical errors – the 1PI $n$-point vertex
functions in the IR computed at $\sigma=0$ and $t=0$ – are rather intensive to
such changes. A possible explanation is the fact that errors from the boundary
at $\sigma_{\mathrm{max}}$ propagate into the computational domain at a finite
speed, which is rather low in the purely diffusive case and in general small
at large $\sigma$, and thus do not influence the physical point at $t=0$ and
$\sigma=0$.

Nevertheless, we conclude from Figs. 10 and 11 that it is extremely important
to use sufficiently large computational domains to minimize numerical errors
in field-dependent FRG flows. This implies that $\sigma_{\mathrm{max}}$ should
be chosen much larger than all relevant scales of the model.

From our findings, it is therefore expected that choosing a large
$\sigma_{\mathrm{max}}$ might even gain in importance in higher-dimensional
models, where the physical point may be located at a non-trivial minimum in
the IR (e.g., the quark-meson model in its broken phase [141]): The 1PI
$n$-point vertex functions are calculated at this physical point. If the
physical point is closer to the boundary of the computational domain the
relative errors for observables might even be larger than for our zero-
dimensional model where the physical point moves towards $\sigma=0$ during the
FRG flow. In terms of errors originating from the boundary at
$\sigma_{\mathrm{max}}$, the physical point at $\sigma=0$ is ideal since it
has the largest spatial and – in a sense causal, due to the finite speed of
propagation – distance to $\sigma_{\mathrm{max}}$.

Lastly, we have to warn that there is no panacea for the construction of a
sufficiently large computational domain and the choice of
$\sigma_{\mathrm{max}}$ has to be adjusted to the specific model and specific
initial condition under consideration. For some problems even more involved
approaches (like using several computational grids of different resolution
$\Delta x$) might be needed or are at least highly advantageous [26, 27]. In
any case one has to check that the IR results do not depend on the size of the
computational domain (even if exact reference values for observables are
unknown), cf. Refs. [130, 129, 32]. This can be done by fixing appropriate
values for the spatial resolution $\Delta x$ as well as for all other
(numerical) parameters and successively increasing $\sigma_{\mathrm{max}}$
until the IR observables do not change anymore.

####  V.1.4 Tests of the UV and IR scales

We now turn to a long-standing discussion in the FRG community, namely the
question: How do we have to choose the initial UV and numerical IR cutoff
scale for the calculation of the RG flow for a specific model?

A common argument is based on the energy scales of a given model. The UV
initial condition is fixed at UV cutoffs $\Lambda$ that are close to the
largest energy scale of the model. Higher $\Lambda$ are excluded by arguing
that at higher energy scales other physical degrees of freedom (e.g., other
interaction channels or new particles) are relevant and the model at hand is
only valid within a certain energy regime. On the other hand, the IR cutoff
$k_{\mathrm{IR}}$ scale is oftentimes fixed by arguing that if it decreases
below the lowest energy scale of the model, the FRG flow is effectively
“frozen in” and the effective potential no longer changes anyway. A relatively
low UV initial scale and a high IR cutoff lead to rather short flow times of
only $t_{\mathrm{UV}}-t_{\mathrm{IR}}\approx 3-4$.

Another approach, which is sometimes employed in conjunction with the first
strategy, is guided by the principle of “numerical stability” of the RG flow,
where cutoffs are chosen in a certain way to “improve performance and
stability” during the numerical RG time integration. In turn, in Refs. [227,
129, 130, 131, 137] relatively small IR cutoff scales are reached due to the
use of numerical stable schemes or the control of stability. Careful
extrapolations into the deep IR like the ones discussed in, e.g., Refs. [227,
26, 27] are another possibility to achieve low IR cutoffs. Note that, for
theories in $d>0$ dimensions, numerical integration into the (deep) IR becomes
very demanding due to multiple reasons, see also Refs. [227, 26, 27, 32]. This
is probably the main reason why often too large numerical IR cutoffs are used.

In general, however, there is a well-defined strategy for the choice of the UV
cutoff scale, which is known as RG consistency. It states that the full
effective action $\Gamma[\Phi]$ in the IR must be independent of the UV
initial scale [96],

$\displaystyle\Lambda\,\frac{\mathrm{d}\Gamma[\Phi]}{\mathrm{d}\Lambda}=0\,.$
(144)

In this framework, one should consider effective models as “stand-alone
theories”, which have their own functional integral and their own well-defined
expectation values. This framework is consistent with our introduction to
zero-dimensional QFTs and FRG in Sec. II and implies that the UV cutoff scale
$\Lambda$ has to be much larger than all scales in the model. (In higher-
dimensional QFTs the “scales of the model” are defined via couplings of
dimension energy as well as external parameters like temperature or chemical
potential.) In this sense, a high initial UV scale is necessary to include all
fluctuations. It was already demonstrated in Ref. [96] that if the initial UV
scale $\Lambda$ is chosen too small and too close to the model scales or
external scales, physical results are spoiled drastically by slightly varying
$\Lambda$ and Eq. (144) is not fulfilled anymore, cf. Refs. [97, 98, 99, 100]
for related discussions in the context of low-energy effective theories of
QCD.

A lower limit for $\Lambda$ arises from the fact that for a given initial
condition $U(t=0,\sigma)$ and at all $\sigma$

$\displaystyle\Lambda+\tfrac{1}{\sigma}\,\partial_{\sigma}U(t=0,\sigma)>\,$
$\displaystyle 0\,,\vphantom{\bigg{(}\bigg{)}}$ (145)
$\displaystyle\Lambda+\partial_{\sigma}^{2}U(t=0,\sigma)>\,$ $\displaystyle
0\,,\vphantom{\bigg{(}\bigg{)}}$ (146)

must hold to have a non-singular flow equation (100). This is discussed, e.g.,
in Refs. [227, 228] and represents a minimal requirement for $\Lambda$ when
considering a given initial condition $U(t=0,\sigma)$. However, guaranteeing
the inequalities (145) and (146) does by itself not guarantee RG consistency
in the sense of Eq. (144).

For higher-dimensional QFTs it is actually complicated to quantify the
relative error of observables from violations of Eq. (144), because “exact”
reference values, e.g., by numerical calculation of expectation values from
the functional integral, are rarely known, especially for low-energy effective
models. In higher dimensions one can only quantify the changes of observables
while varying $\Lambda$ and search for a sufficiently large $\Lambda$, where
Eq. (144) holds. In zero-dimensional QFT this is different, because we can
directly calculate the relative errors for observables like 1PI $n$-point
vertex functions, cf. Eq. (142), for different values of $\Lambda$.

Similar arguments apply to the IR cutoff, where the numerical integration of
the RG flow is stopped. Here, one must clearly state that the full effective
average action $\Gamma[\Phi]$ in the IR is unambiguously defined via the limit
${t\rightarrow\infty\Longleftrightarrow r(t)\rightarrow 0}$ of
$\bar{\Gamma}_{t}[\Phi]$, cf. Eq. (27). In practice, a direct integration to
$t\rightarrow\infty$ is numerically impossible, which implies that one has at
least to make sure that the numerical RG time integration is stopped no
earlier than when all observables of interest do not change anymore, or one
has to systematically extrapolate to $t\rightarrow\infty$, see, e.g., Refs.
[26, 27]. It is worth mentioning that, depending on the specific observable,
these “freeze-out scales” can be extremely different, see Fig. 12.

Figure 12: The RG flow of the minimum $\sigma_{\mathrm{min}}(t)$ (blue) of the
effective potential $U(t,\sigma)$ as well as the FRG flow of the curvature
mass $m_{\sigma}^{2}(t)$ of the $\sigma$ mode (red-dashed) evaluated on the
equations of motion (23) (at the flowing minimum) during the RG flow. The blue
curve sets in after a unique minimum at $\pm\sigma_{\mathrm{min}}(t)$ has
formed. As UV initial condition we use Eq. (139). We used the exponential
regulator (8) with UV cutoff $\Lambda=10^{6}$. The curvature mass
$m_{\sigma}^{2}(t)$ was extracted from $u(t,\sigma)$ via Eq. (153) at the
moving $\sigma_{\mathrm{min}}(t)$. The horizontal (yellow) line denotes the
exact IR result for $\Gamma^{(2)}$ at $\sigma=0$, which must agree with
$m_{\sigma}^{2}$ in the IR, where $\sigma_{\mathrm{min}}(t)=0$.

In the following, we will therefore explicitly explore the influence of UV and
IR cutoff scales on the relative errors (142) for the $\Gamma^{(2n)}$. We
start our discussion by providing results for the relative errors (142)
depending on the RG time $t$ for different $N$ of $O(N)$ and UV initial
condition (139). In Figs. 13, 14, and 15 we plot the relative errors of
$\Gamma^{(2n)}$ for $n=1,2,3$ for $N=1,3,10$, which are all extracted via
various finite-difference stencils from $u(t,\sigma)$ at the physical point
$\sigma=0$ and different $t$ during the FRG flow. All (numerical) parameters
are mentioned in the figures or the respective captions.

Figure 13: The relative error for $\Gamma^{(2m)}$, for $m=1,2,3$, calculated
with the KT scheme as a function of the RG time $t$ for the $O(1)$ model. The
initial UV potential is given by Eq. (139). We use the exponential regulator
(8) with UV cutoff $\Lambda=10^{6}$. The computational grid has 400 cells and
$\sigma_{\mathrm{max}}=x_{\mathrm{max}}=10$. $\Gamma^{(2m)}$ are extracted
from $u(t_{\mathrm{IR}}=60,\sigma)$ via the finite-difference stencils (151) –
(158). Figure 14: The relative error for $\Gamma^{(2m)}$, for $m=1,2,3$,
calculated with the KT scheme as a function of the RG time $t$ for the $O(3)$
model. The initial condition and all other parameters are identical to those
of Fig. 13. Figure 15: The relative error for $\Gamma^{(2m)}$, for $m=1,2,3$,
calculated with the KT scheme as a function of the RG time $t$ for the $O(10)$
model. The initial condition and all other parameters are identical to those
of Fig. 13.

For all three figures (i.e., for $N=1,3,10$) and independent of the choice of
discretization of the numerical derivatives, we observe plateaus for the
relative errors for $\Gamma^{(2n)}$ at the beginning and the end of the FRG
time evolution. The plateau at small $t$ corresponds to the UV regime and
indicates that the UV cutoff is chosen sufficiently large because no
fluctuations are present at the IR physical point until $r(t)$ reaches the
scales of the model. RG consistency (144), hence UV-cutoff independence should
therefore be fulfilled, as long as we initialize our RG flow at some RG scale
which is at the far left of this plateau. Such a plateau at small $t$ is a
sufficient condition for RG consistency but not a necessary one, because
quantum fluctuations could already work at positions in field space away from
the IR physical point and only influence higher-order correlation functions.
We will quantify this within the next paragraphs. In the plots various finite-
difference stencils with distinct error scaling in $\Delta x$ are used to
demonstrate that the plateaus are independent of other sources of errors, like
spatial discretization errors353535Incidentally, Figs. 13, 14, and 15 also
underline our statement that the spatial discretization errors deriving from
the numerical differentiation of $u(t,\sigma)$ are much more severe than the
discretization errors of the KT scheme. Otherwise, the curves for the various
finite-difference stencils would coincide in the IR..

For intermediate $t$, we observe strong dynamics and fast changes in the
relative errors for the $\Gamma^{(2n)}$. The actual values of the relative
errors for intermediate $t$ is irrelevant for the current discussion on UV and
IR scales.

The plateau at late RG times $t$ corresponds to the IR scale of the theory and
indicates that the physical observables are frozen and do not change anymore,
such that the numerical time integration can be stopped. As expected, we find
that the explicit IR scale strongly depends on the choice of $N$, thus the
number of pions and the amount of advection. The smaller $N$ and the more
diffusive the system, the longer it takes to reach the IR363636This is a well-
known observation from all kinds of fluid-dynamical systems. It takes much
longer to reach a (thermal) equilibrium state via diffusion only than by
including advective processes. : For $N=10$ the freeze-out already sets in at
$t\approx 26$, while for $N=1$ one has to wait until $t\approx 30$ to find
that the dynamics ends. This is a difference of two orders of magnitude in the
RG scale. In general, our toy-model tests indicate that rather small IR scales
are needed to actually reach the regime where the observables are frozen.
Still, for $N=10$, $r(t\approx 26)\approx 5\cdot 10^{-6}$, i.e., the IR regime
begins six orders of magnitude below the model scales.

This observation might also partially translate to higher-dimensional models,
meaning that commonly used IR cutoffs might be systematically chosen too
large, such that predictive power is lost. Nevertheless, we expect this
problem to be the less severe the higher the space-time dimensionality of a
model under consideration, because of the larger phase-space (momentum
suppression)373737The smaller the space-time dimension of a model, the more
important are long-range interactions – quantum fluctuations at small RG
scales $k$ – for the macroscopic observables, which is of course most extreme
for $d=0$. Furthermore, field space interactions tend to become more important
for a small number of space-time dimensions..

Furthermore, we observe from Fig. 12 as well as Figs. 13, 14, and 15 that the
freeze-out scale is slightly different for different observables, because
higher 1PI $n$-point vertex functions seem to be more sensitive to tiny
changes in $u(t,\sigma)$. In particular, we observe from Fig. 12 that the
minimum $\sigma_{\mathrm{min}}$ is already frozen at $t\approx 14$, while the
curvature mass $m_{\sigma}^{2}$ still changes drastically after $t\approx 14$
over several orders of magnitude in RG scale. This is especially interesting
for higher-dimensional models: Oftentimes the freeze-out of the minimum is
considered a suitable IR scale to stop the FRG flow, which is definitely not
justified, since the derivatives of the potential at the physical point are
usually still changing. Using the changing rates of the curvature mass instead
of the position of the minimum as a monitor for the dynamic range – viable
numerical IR cutoffs – has proven crucial in the FRG study [32] of the Gross-
Neveu(-Yukawa) model by several authors and their collaborators.

Next, we explicitly quantify the relative errors of $\Gamma^{(2n)}$, which
derive from too small UV cutoffs $\Lambda$ and the violation of RG consistency
(144). To this end, we plot the relative errors (142) as a function of the UV
cutoff $\Lambda$, while keeping the IR cutoff scale fixed at extremely small
$r(t_{\mathrm{IR}})=10^{-20}$.

Figure 16: The relative error for $\Gamma^{(2m)}$ for $m=1,2,3$ from the KT
scheme as a function of the UV cutoff scale $\Lambda$ for the initial
potential (139). We use the exponential regulator (8) and keep the IR cutoff
scale constant at $r(t_{\mathrm{IR}})=10^{-20}$ for all runs. Furthermore, for
all data points the computational grid size is fixed at
$\sigma_{\mathrm{max}}=x_{\mathrm{max}}=10$ and the number of volume cells is
set to $n=400$. $\Gamma^{(2m)}$ are calculated from
$u(t_{\mathrm{IR}}=60,\sigma)$ via the approximations (151), (154), and (157)
for the numerical derivative. The yellow straight line is for optical
guidance.

In Fig. 16 we observe that the IR observables become independent of $\Lambda$
at rather large $\Lambda\approx 10^{6}$. This is several orders of magnitude
above the model scales, contrary to what is often used in FRG studies in
higher dimensions. If the UV cutoff is chosen too small, we find that the
relative errors of $\Gamma^{(2n)}$ grow proportional to $\frac{1}{\Lambda}$,
as estimated in Eq. (33). Surprisingly, it turns out that the RG-consistency
condition (144) is already violated at rather large UV cutoff scales
$\Lambda\approx 10^{5}$ and is only fulfilled for $\Lambda\gtrsim 10^{5}$. We
conclude that it requires great care when specifying the UV scale in a FRG
calculation.

Before we close this discussion, we provide a natural measure to estimate the
correct UV and IR scales of a model or theory, even if there are no exact
reference values for observables that can be used for comparison with the FRG
results. To this end, we plot in Fig. 17 the shifted logarithm of the changing
rates of the $\bar{\Gamma}^{(2n)}(t)$ at the IR minimum $\sigma=0$ over RG
time $t$. These quantities have to vanish in the UV and the IR, when the
relative errors (142) are not changing.

Figure 17: The rate of change in $t$ of $\bar{\Gamma}^{(2m)}(t)$ at the IR
minimum $\sigma=0$ for $n=1,2,3$ during the RG flow. This rate of change is
defined as the numerical RG time derivative
$\partial_{t}\bar{\Gamma}^{(2m)}(t)$ over the RG time.
$\partial_{t}\bar{\Gamma}^{(2m)}(t)$ are calculated via a finite-difference
approximation ${[\bar{\Gamma}^{(2m)}(t)-\bar{\Gamma}^{(2m)}(t-\Delta
t)]/\Delta t}$, where $\Delta t=0.2$. $\bar{\Gamma}^{(2m)}(t)$ are obtained
via numerical derivatives (151), (154), and (157) of $u(t,\sigma)$ at
$x=\sigma=0$. For convenience, we added $1$ and took the logarithm to
pronounce the regions of high changing rates of the observables
$\bar{\Gamma}^{(2m)}(t)$ and to identify the freeze-out plateau, where they
have to vanish. We used the exponential regulator (8) with UV cutoff
$\Lambda=10^{6}$.

A similar investigation can be done for any other model or theory and can be
used as an indication to ensure sufficiently large UV and sufficiently small
IR cutoffs: A first estimate may be obtained by choosing $\Lambda$ and
$t_{\mathrm{IR}}$ in a way that the plateaus (or scaling regimes) in figures
similar to Fig. 17 are of approximately equal RG time duration than the time
interval in which the actual dynamics takes place. In the absence of an
explicit and accessible error estimate rates of change are a cheap and simple
tool to study the UV and IR limits of RG time evolution, cf. Ref. [32].

###  V.2 Test case II: $\phi^{4}$ theory

The second test case is a zero-dimensional version of $\phi^{4}$ theory with a
“Mexican hat”-type potential well-known from standard textbook discussions of
spontaneous symmetry breaking [117, 112]. The initial UV potential reads

$\displaystyle
U(\vec{\varphi}\,)=-\tfrac{1}{2}\,\vec{\varphi}^{\,2}+\tfrac{1}{4!}\,(\vec{\varphi}^{\,2})^{2}\,.$
(147)

The corresponding initial condition for the FRG flow is illustrated in Fig.
18. The reference values for the exact IR 1PI vertex functions $\Gamma^{(2n)}$
of the $O(N)$ model (73) – (75) are calculated numerically from the UV
potential (147) and are listed for selected values of $N$ in Table 2, which
also includes reference values at $N=4$ for the analogous potential with
positive mass term $+\tfrac{1}{2}\vec{\varphi}^{\,2}$.

In the following subsubsections, we discuss results obtained with the KT
scheme and numerical results computed using the Taylor (vertex) expansion
outlined in Sub.Sub.Sec. III.3.2.

Figure 18: The UV potential $U(\sigma)$ (red-dashed) and its first derivative $u(\sigma)=\partial_{\sigma}U(\sigma)$ (blue-solid) of our test case (147) evaluated on the constant background field configuration.  Table 2: The exact results for $\Gamma^{(2n)}$ for the $O(N)$ model with the initial UV potential (147) for $N=1$ and $N=4$ as well as $N=4$ with a positive mass term. The results are obtained by a high-precision one-dimensional numerical integration of the expectation values $\langle(\vec{\phi}^{\,2})^{n}\rangle$ using Mathematica’s numerical integration routine NIntegrate [138] with a PrecisionGoal and AccuracyGoal of $10$. Here, we present the first ten digits only.  $N$ | $\Gamma^{(2)}$ | $\Gamma^{(4)}$ | $\Gamma^{(6)}$
---|---|---|---
$1$ | $0.1995098930$ | $0.0622583604$ | $0.1077442107$
$4$ | $0.5064440744$ | $0.1824153741$ | $0.2802880035$
$4$383838Using $U(\vec{\varphi}\,)=+\tfrac{1}{2}\,\vec{\varphi}^{\,2}+\tfrac{1}{4!}\,(\vec{\varphi}^{\,2})^{2}$, thus Eq. (147) with a positive mass term. | $1.5809246562$ | $0.6118483669$ | $0.5686310448$

####  V.2.1 Results obtained using the KT scheme

In this subsubsection we will discuss selected numerical results of the
application of the KT scheme for the analytic initial condition (147). We have
performed the full set of numerical tests discussed in Sub.Sec. V.1 for this
initial condition and found results supporting the general statements made
there. For the sake of brevity, we will therefore not repeat the complete
discussion of that section.

##### UV and IR scales:

In Fig. 19 we present the RG flow of the derivative of the effective potential
$u(t,\sigma)$ from the UV (blue) to the IR (red). For the smooth initial
condition – in the absence of large gradients – the highly non-linear
advection and diffusion contribute almost an equal amount to the dynamics.
Between $t\approx 25$ and $t\approx 30$ we observe significant changes in the
shape of the potential: the non-trivial minimum moves towards $\sigma=0$ and
vanishes at $t\approx 28$ resulting in a convex potential with a trivial
minimum at $\sigma=0$ as expected and required. At small and large $t$ outside
the apparent dynamic range between $t\approx 25$ and $t\approx 30$ we observe
only very marginal changes in Fig. 19.

Figure 19: The FRG flow of the effective potential $U(t,\sigma)$ (upper panel)
and its derivative $u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$ (lower panel)
for the zero-dimensional $O(4)$ model with initial condition Eq. (147),
evaluated at $t=0,\,2,\,4,\,\ldots,\,60$ (integer values for $t$ were chosen
for convenience and readability). The (overlapping) blue and violet curves
correspond to the UV and the red curves to the IR. We used the exponential
regulator (8) with UV cutoff $\Lambda=10^{12}$. The plot does not show the
region $x=5$ to $x=10$, because the tiny differences between $u(t,\sigma)$ and
$u(t_{\mathrm{UV}},\sigma)$ are not visible in this region and vanish for
large $x=\sigma$ anyhow.

A close inspection of the relative errors for the first three non-vanishing
$n$-point vertex functions in Fig. 20 reveals that actually the relevant
dynamics sets in much earlier at $t\approx 10$ and the values for the
$n$-point vertex functions freeze out at late times around $t\approx 40$,
which is due to the diffusion close to $\sigma=0$. On the level of
$u(t,\sigma)$ these subtle changes in the $n$-point vertex functions cannot be
observed by a simple visual inspection of Fig. 19.

The plateaus in the UV (at small $t$) and the IR (at large $t$) support the
choice of $\Lambda=10^{12}$ and $t_{\mathrm{IR}}=60$ to be valid initial UV
and IR cutoff scales in terms of RG consistency. The present initial UV scale
is larger when compared to $\Lambda=10^{6}$ used for most computations
involving the non-analytic potential considered in the previous section.
($\Lambda=10^{6}$ corresponds to $t\approx 14$ in the present case, which is
already in the dynamic region in Fig. 20.) Hence, the inclusion of a quartic
interaction term in Eq. (147) seems to require higher initial UV scales to
ensure RG consistency. This supports the statements made in Sub.Sub.Sec.
V.1.4: RG consistency and UV and IR scales have to be re-evaluated when
changing the initial condition in the UV or the model under consideration,
since characteristic internal scales then also change.

Figure 20: The relative error for $\Gamma^{(2m)}$, for $m=1,2,3$, calculated
with the KT scheme as a function of the RG time $t$ for the $O(4)$ model. The
initial UV potential is given by Eq. (147). We use the exponential regulator
(8) with UV cutoff $\Lambda=10^{12}$. The computational grid has 400 cells and
$\sigma_{\mathrm{max}}=x_{\mathrm{max}}=10$. $\Gamma^{(2m)}$ are extracted
from $u(t_{\mathrm{IR}}=60,\sigma)$ via the finite-difference stencils (151),
(154), and (157).

##### Computational grid and domain size:

We conclude this subsubsection with a brief discussion regarding the
computational grid and domain size. The relative error for the first three
non-vanishing $n$-point vertex functions is shown as a function of the cell
size $\Delta x$ in Fig. 21. For the two-point function we recover a perfect
error scaling with $\Delta x^{2}$ down to extremely small $\Delta x$. The last
data point in Fig. 21 is at ${\Delta x\approx 3.3\cdot 10^{-3}}$ corresponding
to $n=3000$ cells. For the two-point function the rounding errors of the
employed finite-difference extraction (151) for $\Gamma^{(2)}$ and the finite
precision of the ODE integrator (NDSolve from Mathematica [138] with a
PrecisionGoal and AccuracyGoal of 10) seem to be small for all depicted
$\Delta x$ in this scenario. A comparison with the present perfect error
scaling for $\Gamma^{(2)}$ supports the comments made about discretization
errors for the discontinuous initial condition (139) in Fig. 9. For the
higher-order $n$-point vertex functions $\Gamma^{(4)}$ and $\Gamma^{(6)}$,
however, we find that rounding errors related to the finite-difference
extractions (154) and (157) limit the achievable precision. Again, we identify
$\Delta x\approx 0.025$ as an optimal cell size for the extraction of
$\Gamma^{(4)}$ and $\Gamma^{(6)}$ but note that typical relative errors for
$\Gamma^{(6)}$ are rather high at $\approx 4\%$ around $\Delta x\approx
0.025$.

Figure 21: The relative error as a function of the cell size $\Delta x$ for
the numerical results (blue dots) from the KT scheme for the coefficients
$\Gamma^{(2n)}$ for $n=1,2,3$ with initial potential (147). The numerical
derivatives at $\sigma=0$ of $u(t_{\mathrm{IR}}=60,\sigma)$ were calculated
via the second-order accurate central schemes (151), (154), and (157). Here,
$x_{\mathrm{max}}=10$, but we could have used any sufficiently large
$x_{\mathrm{max}}$. We used the exponential regulator (8) with UV cutoff
$\Lambda=10^{12}$. The yellow straight lines are for optical guidance.

In Fig. 22, we study the effect of the size of the computational domain
$x_{\mathrm{max}}$ on the achievable relative errors for $\Gamma^{(2)}$,
$\Gamma^{(4)}$, and $\Gamma^{(6)}$ at a constant $\Delta x=0.025$. One major
difference between the $\phi^{4}$ potential (147) studied in this section and
the non-analytic potential (139) of the previous subsection is their
asymptotic behavior for large $\sigma$. For large $\sigma$ the leading-order
term of the $\phi^{4}$ potential is – as the name suggests – quartic while the
non-analytic potential of the previous section grows only $\sim\sigma^{2}$. In
terms of the conserved quantity $u=\partial_{\sigma}U$ one might expect
problems when using a linear extrapolation for the ghost cells at large
$\sigma$ as discussed in Sub.Sub.Sec. IV.4.2 with a potential where $u$ grows
$\sim\sigma^{3}$ for large $\sigma$. For the non-analytic initial condition
(147) we avoided this possible source of error by construction. However,
considering the results plotted in Fig. 22 together with the perfect error
scaling displayed in the previous Fig. 21, we conclude that a linear