††thanks:  J. Stoll, N. Zorbach, A. Koenigstein, and M. J. Steil contributed
in equal shares to this work. S. Rechenberger was involved in early stages of
this project before leaving academia and contributed to the final draft.

#  Bosonic fluctuations in the $(1+1)$-dimensional Gross-Neveu(-Yukawa) model
at varying $\mu$ and $T$ and finite $N$

Jonas Stoll<EMAIL_ADDRESS>Institut für Theoretische Physik,
Johann Wolfgang Goethe-Universität,
Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany  Niklas Zorbach
<EMAIL_ADDRESS>Institut für Theoretische Physik, Johann
Wolfgang Goethe-Universität,
Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany  Adrian Koenigstein
<EMAIL_ADDRESS>Institut für Theoretische Physik, Johann
Wolfgang Goethe-Universität,
Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany  Martin J. Steil
<EMAIL_ADDRESS>Technische Universität Darmstadt,
Department of Physics, Institut für Kernphysik, Theoriezentrum,
Schlossgartenstraße 2, D-64289 Darmstadt, Germany  Stefan Rechenberger
<EMAIL_ADDRESS>Institut für Theoretische Physik, Johann Wolfgang
Goethe-Universität,
Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany

###### Abstract

Using analogies between flow equations from the Functional Renormalization
Group and flow equations from (numerical) fluid dynamics we investigate the
effects of bosonic fluctuations in a bosonized Gross-Neveu model – namely the
Gross-Neveu-Yukawa model. We study this model for finite numbers of fermions
at varying chemical potential and temperature in the local potential
approximation. Thereby we numerically demonstrate that for any finite number
of fermions and as long as the temperature is non-zero, there is no
$\mathbb{Z}_{2}$ symmetry breaking for arbitrary chemical potentials.

Gross-Neveu model, Gross-Neveu-Yukawa model, bosonic fluctuations, phase
diagram, $\mathbb{Z}_{2}$ symmetry, chiral symmetry, symmetry
breaking/restoration, finite $N$, FRG, numerical fluid dynamics

###### Contents

  1. I Introduction
    1. I.1 Phenomenology of the Gross-Neveu model
    2. I.2 To break or not to break – $\mathbb{Z}_{2}$ symmetry at finite $N$
    3. I.3 Main results
    4. I.4 Structure
  2. II The Gross-Neveu(-Yukawa) model
    1. II.1 In vacuum
    2. II.2 In medium
  3. III The Gross-Neveu-Yukawa model via the Functional Renormalization Group
    1. III.1 The Functional Renormalization Group
    2. III.2 Truncation and Renormalization Group flow equation(s)
    3. III.3 Comment on regulators
    4. III.4 Truncations and initial conditions – Differences and similarities of GN, bGN and GNY models
    5. III.5 Dynamical symmetry breaking and RG flows
  4. IV The FRG and (numerical) fluid dynamics
    1. IV.1 The LPA flow equation as a non-linear heat equation with sink/source term
    2. IV.2 Irreversibility and entropy production
    3. IV.3 Field space boundary conditions
    4. IV.4 The chemical potential as a “shock wave” in field space
    5. IV.5 Numerical implementation
  5. V Infinite-$N$ analysis within the FRG – Consistency checks and the UV initial condition
    1. V.1 Mean-field, infinite-$N$, and FRG
    2. V.2 The mean-field potential and asymptotic freedom
    3. V.3 Analytical results for the mean-field phase diagram in the renormalized limit $\Lambda\rightarrow\infty$
    4. V.4 UV initial condition for the effective potential for RG flows with and without bosonic fluctuations
    5. V.5 (Numeric) consistency check with mean-field results
  6. VI Bosonic quantum fluctuations in the Gross-Neveu-Yukawa model – Results
    1. VI.1 Symmetry breaking and restoration during the RG flow
    2. VI.2 The role of the chemical potential in the fluid dynamic setup
    3. VI.3 Varying the number of flavors $N$
    4. VI.4 Varying the temperature $T$
    5. VI.5 Varying the chemical potential $\mu$
    6. VI.6 Computations in vacuum
    7. VI.7 The phase diagram
  7. VII Conclusion and outlook
  8. A Conventions
    1. A.1 Metric and Clifford algebra
    2. A.2 Relation between the Minkowski functional integral and Euclidean partition function
    3. A.3 Fourier transformation
  9. B Hubbard-Stratonovich Transformation and bosonization of the GN model
  10. C Quantum fields on a cylinder – the thermal GNY model
  11. D The LPA flow equation at non-zero $\mu$ and $T$
  12. E Zero temperature calculations
  13. F Numerical tests
    1. F.1 Size of the computational domain
    2. F.2 Spatial resolution
    3. F.3 UV-cutoff independence
  14. G Discretization schemes for the source/sink term

##  I Introduction

Since its original publication in 1974, the Gross-Neveu (GN) model [1] – a
relativistic QFT of $N$ massless Dirac fermions that are self-interacting via
the scalar channel of the four-Fermi interaction – was subject of intensive
research w.r.t. various aspects of strongly interacting systems. In $1+1$
space-time dimensions it was studied as a (perturbatively) renormalizable
prototype model for asymptotic freedom [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
while for $2+1$ space-time dimensions it served as a toy model for
asymptotically safe QFTs [12] and is renormalizable with non-perurbative
methods [13, 14, 15, 16]. In $3+1$ dimensions it mimics the dynamics of the
important scalar channel of the Nambu-Jona-Lassinio model [17, 18] and is non-
renormalizable, but can be seen as being part of a low-energy effective model
of QCD, cf. the reviews [19, 20]. However, within this work and from now on,
we exclusively focus on its original formulation in $1+1$ space-time
dimensions, if not explicitly stated otherwise.

Furthermore, there are various connections to non-relativistic models in
solid-state physics. Indeed, GN-type models can arise in the continuum limit
(describing large distance physics) of one-dimensional solid-state systems.
Examples of such one-dimensional systems are quantum antiferromagnets
(described by spin-$\tfrac{1}{2}$ Heisenberg models [21]), interacting
electrons in a one-dimensional conductor (described by a Tomonaga–Luttinger
liquid [22, 23]), and polyacetylene polymer chains [24] –
${(\mathrm{C}\mathrm{H})_{x}}$ –(described in the limit $x\rightarrow\infty$
with the Su-Schrieffer-Heeger model [25] or the subsequent Takayama-Lin-Liu-
Maki model [26]). Recently, a one-dimensional probabilistic cellular
automaton, where classical bits can be interpreted as Ising spins, was shown
to be equivalent to a relativistic fermionic quantum field theory [27]. A
central concept behind the emergence of fermionic QFTs in two dimensions from
spin-systems in one dimension is the mapping of spin operators onto fermionic
creation and annihilation operators by means of the Jordan–Wigner
transformation [28], for a pedagogic discussion see Ref. [29]. GN type models
also arise naturally in the continuum limit of two-dimensional spin-systems
like the Ashkin-Teller model [30], which is related to the well-known Potts
model [31] and as such used to study various phenomena of solid-state physics,
see, e.g., Ref. [32] for a general overview. In the following paragraphs we
will list some explicit connections between the GN model and models used in
solid-state physics.

The lattice field theoretical formulation of the GN model in the limit
$N\rightarrow\infty$ is equivalent to an Ising model [33].

At finite $N$ the GN model can be considered as the continuum limit of the
$N$-color ($N$ Ising spin) Ashkin-Teller model [30, 34, 35], which describes
$N$ coupled Ising spins on a two-dimensional lattice [35].

For $N=1$ the GN model is equivalent to the Thirring model [36, 37] due to
Fierz identities, cf., Refs. [9, 10, 29]. In the continuum limit the one-
dimensional spin-$\tfrac{1}{2}$ Heisenberg model is equivalent to a $N=1$ GN
(Thirring) model [29]. The Thirring model also arises in the infinite volume
limit of the Luttinger model [23] with strictly local interactions [38], see
Ref. [29] for further details. The massive Thirring model is equivalent to the
sine-Gordon model [39, 40], which in turn is (among its other application in
mathematical physics) the continuum-limit of the Frenkel–Kontorova model [41].
The latter is a simple model of a harmonic chain in a periodic potential known
from solid-sate physics [42]. The equivalence of the Thirring model and sine-
Gordon model is based on an Abelian bosonization transformation, see Refs.
[29, 40] and references therein, which connects equivalent bosonic and
fermionic two-dimensional quantum field theories.

Apart form research efforts in high energy and solid-state physics the GN
model is also of interest in the context of holographic methods [43, 44]
especially in the study of the AdS/CFT correspondence involving higher spin
fields, see Ref. [45] for a recent review.

Despite this broad range of applications and interconnections as well as
intensive research, there are still open questions. Some of these are
addressed within this work.

###  I.1 Phenomenology of the Gross-Neveu model

A peculiar feature of the massless GN model is, that at leading order of an
$\tfrac{1}{N}$-expansion, thus in the infinite-$N$ limit (sometimes also ’t
Hooft limit), the GN model dynamically develops a mass gap for the fermions,
which is associated with an anti-fermion-fermion condensate
$\langle\bar{\psi}\,\psi\rangle\neq 0$. In turn, this results in the breakdown
of the discrete chiral symmetry of the initial microscopic ultraviolet (UV)
theory. The formation of a mass gap is a purely non-perturbative effect, see,
e.g., the Refs. [8, 10] for details, and is a prime example for dimensional
transmutation – the emergence of a dimensionful scale in a theory which has
only dimensionless couplings in its UV classical action, see, e.g., Refs. [46,
9, 10]. Hence, by summing up loop-contributions of all orders in the four-
Fermi coupling in a $\tfrac{1}{N}$-expansion the discrete chiral symmetry
spontaneously breaks down and is absent in the macroscopic theory in the
infrared (IR). In the partially bosonized version, this amounts to integrating
out the fermion-loop contribution to the bosonic effective potential, which
develops a non-trivial minimum in the IR – the condensate [8, 10, 3, 7, 11].

Shortly after D. J. Gross and A. Neveu had published their results [1], the
question came up, to what extent condensate formation is stable against
thermal effects due to non-zero temperature $T$ or an increase in density,
induced by a non-zero quark chemical potential $\mu$. Within the infinite-$N$
limit and allowing only for spatially homogeneous condensates, the answer to
these questions was quickly settled and is remarkable [47, 48, 49, 50, 51, 52,
53, 54, 55]: The phase diagram of the GN model in the $\mu$-$T$-plane is
actually unique. It consists of a region, where the discrete chiral symmetry
is broken and a region of vanishing chiral condensate, cf. Fig. 2. The phase-
transition line between these regions splits up into a second-order phase
transition (starting at $\mu=0$ and some critical temperature
$T_{\mathrm{C}}\neq 0$ and ending in a critical point $\mu_{\mathrm{CP}}\neq
0$ and $T_{\mathrm{CP}}\neq 0$) and a first order phase transition (starting
at the critical point and ending on the $T=0$ axis and some non-zero chemical
potential $\mu_{1}$). Here, being “unique” means that the entire phase diagram
solely depends on a single dimensionful parameter, which is related to a
renormalization condition, and can be chosen freely. All other dimensionful
quantities are fixed multiples of this parameter and choosing a different
renormalization condition (fixing some other parameter) corresponds to simple
rescalings of all dimensionful quantities, but does not change their ratios,
the phenomenology and the phase diagram.

Notwithstanding these successes, the discussion on the physics of the GN model
did not stop. One of the assumptions, which has lead to the above results, is
the assumption of a spatially homogeneous condensation of the fermions.
Relaxing this assumption, but retaining the $N\rightarrow\infty$ limit, it was
shown in Refs. [56, 57, 58] that there are regions in the $\mu$-$T$-plane,
where a spatially inhomogeneous but static condensate
$\langle\bar{\psi}\,\psi\rangle(x)\neq\mathrm{const}.$ is energetically
favored over homogeneous condensation including a vanishing of the condensate.
These so-called inhomogeneous phases attracted a lot of attention within the
last decades in the high energy physics community and seem to be a rather
robust feature of a lot of effective models of strongly correlated systems at
least in the infinite-$N$ limit [59, 60, 61]. We are not going to focus and
discuss the possibility for spatially inhomogeneous condensation in this work,
but instead concentrate on the other main assumption, which underlies the
above discussions – namely the infinite-$N$ limit. This limit for the GN model
leads to a complete suppression of bosonic quantum fluctuations.

###  I.2 To break or not to break – $\mathbb{Z}_{2}$ symmetry at finite $N$

The big question that immediately comes to mind if one relaxes the
infinite-$N$ limit and studies the GN model for finite $N$ is, if spontaneous
symmetry breaking (SSB) and condensation still takes place, especially when
medium effects (non-zero $\mu$ and/or $T$) are included.

At first sight, it seems as if already the Coleman-Mermin-Wagner-Hohenberg
theorem [62, 63, 64] forbids the formation of a condensate at non-zero
temperatures. Though, the theorem strongly relies on the presence of massless
Nambu-Goldstone bosons [65, 66, 67], which are only included in extensions of
the GN model, e.g., the chiral GN model [68, 69, 70] with continuous chiral
symmetry or other related models [37, 3]. We therefore believe that one should
exercise caution, when arguing directly with this theorem.

Nevertheless, a qualitative argumentation was put forward already by L. D.
Landau in 1950 [71], which should in principle forbid also discrete chiral
symmetry breaking in one-spatial dimension at non-zero temperature. L. D.
Landau and E. M. Lifshitz argue in Chap. 163 of Ref. [71] that for systems in
one dimension of infinite extend only one phase can exist at $T>0$, since
coexistence of more than one phase is energetically disfavored. To be
concrete, they considered a bistable system in one dimension of infinite
extend at non-zero temperature with $n$ interfaces between the two possible
phases per length $L$ and showed that the thermodynamic potential of the
system can be decreased by increasing the concentration of interfaces $n/L$,
which is directly related to the entropy, assuming a finite interface energy.
Thus the system breaks down into a macroscopic number of domains, which
rendered macroscopic phase coexistence at non-zero temperature impossible.
Landau’s argument can be applied to a broad range of effectively one-
dimensional systems including the Ising model in one spatial dimension [72] at
$T\neq 0$, which is always in its $\mathbb{Z}_{2}$ symmetric phase, cf. Refs.
[10, 3]. A dedicated discussion of Landau’s argument in the context of one-
dimensional systems can be found in Ref. [73].

In their study [50] of the GN model based on a large-$N$ expansion R. F.
Dashen, S. Ma, and R. Rajaraman were able to confirm Landau’s argument, see
also Ref. [74]. They found no SSB for any small but non-zero temperature and
finite $N$. Using in parts heuristic arguments, they showed that the entropic
gain of a field configuration of alternating kinks is large enough at finite
$N$ and $1/T$ to be energetically preferred over a homogeneous
configuration111Homogeneous and inhomogeneous field configurations in this
context refer to the configurations used to evaluate the partition function in
a saddle point approximation, see Ref. [50] for details on their computation,
and they are not to be confused with homogeneous and inhomogeneous classical
field configurations $\langle\bar{\psi}\psi\rangle(x)$ discussed in the
previous Sub.Sec. I.1.. Those field configurations alternating in kinks have a
vanishing chiral condensate $\langle\bar{\psi}\,\psi\rangle=0$. In the
infinite-$N$-limit (mean-field), the energy per kink becomes infinite and
consequently the density of kinks approaches zero realizing a homogeneous
field configuration compatible with $\langle\bar{\psi}\,\psi\rangle>0$
allowing for SSB. This is no contradiction to Landau’s argument since the
latter only holds assuming finite interface energies [71, 73]. Although R. F.
Dashen, S. Ma, and R. Rajaraman argue that $T_{\mathrm{C}}=0$ at finite $N$,
they do not discuss the situation for $T=0$ and $\mu\geq 0$ for finite $N$ in
Ref. [50].

Another discussion on the absence of symmetry breaking in one-spatial
dimension at non-zero $T$ can be found in Ref. [37] by E. Witten, which
discusses the absence of a phase with spontaneous symmetry breaking and long-
range order in accordance with the Coleman-Mermin-Wagner-Hohenberg theorem
[62, 63, 64] and the related possibility for a phase of Berezinski-Kosterlitz-
Thouless type [75, 76] with quasi long-range order for the $\mathrm{SU}(N)$
Thirring model [36, 37].

Other authors, e.g., U. Wolff [52], argue based on the duality of the spatial
and Euclidean time direction in $1+1$ dimensions for $T>0$ and $\mu=0$: The
thermal GN model in $1+1$ dimensions is equivalent to a QFT with a finite
spatial volume (but infinite Euclidean time direction) and hence no SSB takes
place since it is canonically considered as an effect only present in systems
with infinite volumes, see, e.g., Ref. [77]. While this reasoning – the
general absence of SSB in finite systems – might strictly speaking be sound,
sufficiently large volumes, the inclusion of small (possibly infinitesimal)
explicit symmetry breaking, or subtleties of the thermodynamic and/or infinite
volume limit, cf. Ref. [78], can lead to signatures reminiscent of SSB.

Indeed, only recently some of our colleagues and collaborators found some
indications via numerical lattice Monte-Carlo simulations [79, 80, 81, 82],
that some (inhomogeneous) condensation phenomena in the massless bosonized GN
model at finite $N$ and non-zero $\mu$ and $T$ still seem to be present.
Similar results were already found in earlier lattice Monte-Carlo studies [83,
84, 85] at finite $N$. However, in all of these works, the above arguments by
Landau et al. against condensation at non-zero $T$ could neither be completely
ruled out nor be confirmed. In fact, most of the results suffer from the facts
that proper continuum and infinite volume extrapolations were not performed.
Consequently finite volume effects and discretization artifacts limit the
predictive power of those results for the continuum theory in an infinite
volume. The finite sized spatial domain (and the related boundary conditions)
might have prevented a sufficient resolution of long-range fluctuations, which
are however of uttermost importance for condensation and in this context
especially vaporization phenomena in low-dimensional systems.

Recent lattice results presented in Refs. [79, 80, 81, 82] have sparked
further lattice studies of four-Fermi models in $1+1$ dimensions: We are aware
of these parallel developments and computations using lattice Monte Carlo
simulations in the GN model and related models (especially the chiral GN
model) in Refs. [86, 87], which are however not completed yet and therefore
omitted in the following discussion. For the chiral GN model, we expect some
interesting dynamics at finite $N$ due to competing effects from the Coleman-
Mermin-Wagner-Hohenberg theorem and the $U(1)_{\mathrm{A}}$ anomaly [70].

All of this lead us to the idea to study the phenomenon of $\mathbb{Z}_{2}$
symmetry breaking and/or restoration in the GN model at finite $N$, $T\geq 0$
and also $\mu\geq 0$ but in an infinite spatial volume within a different
framework – namely within the Functional Renormalization Group (FRG).222After
completing this work, we became aware of Ref. [88], where next-to-leading
order corrections of the $\tfrac{1}{N}$-expansion to the effective potential
of the GN model were calculated. Finding that this expansion breaks down in
the vicinity of mean-field critical temperature $T_{\mathrm{C}}$ the authors
also suggest to analyze the GN model within the FRG framework, which is the
main focus of this paper. We thank J. Braun for drawing our attention to this
interesting publication. We wanted to find out, if it is possible to
(numerically) confirm the arguments by Landau et al. against symmetry breaking
in the GN model or if there are some other competing effects, which are not
captured in the aforementioned mostly qualitative/heuristic discussions, that
allow for symmetry breaking or some long-range ordering.

###  I.3 Main results

Using the FRG method in the local potential approximation (LPA) we find no
spontaneous symmetry breaking in the Gross-Neveu-Yukawa (GNY) model, which is
a bosonized version of the GN model that is used for practical computations,
for any finite number of fermions, non-zero temperature and arbitrary chemical
potential. Our results therefore support the qualitative arguments by Landau
et al.. While fermionic interactions dominate the dynamics during the RG-flow
at high momenta and lead to the formation of a condensate – similar to the
situation at infinite-$N$ and finite $N$ lattice Monte-Carlo studies – long-
range bosonic quantum fluctuations restore the discrete chiral symmetry in the
deep infrared regime – at low momenta – of the RG flow independent of the
number of fermions and the value of the chemical potential for all non-zero
temperatures.333This dynamics is analogous to the phenomenon of
precondensation, which was first discussed in Refs. [89, 90, 91, 92, 93, 94]
in the context of RG flows. This is in accordance with the predictions of R.
F. Dashen, S. Ma, and R. Rajaraman put forward in Ref. [50]. The result
further supports the similarities between the one-dimensional Ising model [72]
and the GN model, because we find that the critical temperature for a phase
transition to the symmetry broken phase is $T_{\mathrm{C}}=0$ for all $\mu$
and finite $N$. The symmetry restoration deep in the IR might present a
substantial challenge when trying to resolve infinite volume physics with
simulations or computations in finite volumes.

Computations in the vacuum ($T=\mu=0$) limit indicate SSB even at finite $N$.
Direct computations at $T=0$ and $\mu>0$ could not be performed at finite $N$
within this work but computations at small $T$ and $\mu>0$ indicate, after
extrapolation to $T=0$ and under consideration of the vacuum results, a
possible quantum phase transition at non-zero chemical potential between a
phase of spontaneous broken symmetry at low chemical potentials and a restored
phase at high chemical potentials at $T=0$ and finite $N$.

Furthermore, as a consistency check, we show analytically that the mean-field
approximation of the LPA flow equation in the FRG framework, i.e. neglecting
bosonic fluctuations in the limit $N\rightarrow\infty$, is equivalent to the
canonical infinite-$N$ results in literature cf. Refs. [52, 95]. We use this
fact to check our numerics by reproducing the infinite-$N$ phase diagram.

To obtain numerical results at finite $N$, we use a novel reformulation of the
LPA flow equation in terms of a conservation equation and solve it by using a
finite-volume-method444We emphasize that finite-volume refers here to a
discretization scheme for partial differential equations and not to a finite
volume of the physical spatial direction of the GN model.. In this setup
fermions act like a sink/source term and introduce a moving discontinuity at
zero (low) temperatures and non-zero chemical potentials, which makes it
(practically) impossible to numerically calculate to arbitrary low RG-scales.
This might indicate a fundamental conceptual problem and prevents practical
numerical computations at $T=0$ and $\mu>0$ in this work.

###  I.4 Structure

We tried to keep the discussion in this work as self-contained as possible.
This should enable readers not familiar with all aspects of this work to
follow the entire line of argumentation without resorting to the literature.
This goal is the reason behind the substantial scope of this paper. In this
subsection we provide an overview of the structure of this work.

In Sec. II we introduce the GN model and its bosonized counterpart the GNY
model in vacuum and medium as well as their differences. Readers familiar with
the model may just skim over this section. In Sec. III we give an introduction
to and provide details on the FRG approach to the GN(Y) model. Readers
familiar with the FRG can skip Sub. Sec. III.1 and III.5 but might be
interested in the remainder of the section. Readers who are not interested in
the technical details of the FRG formalism in the context of this work may
omit Sub.Secs. III.3 and III.4 from their first reading. In Sec. IV we discuss
the relation between RG flow equations and (numerical) fluid dynamics and our
implementation. Readers who are not interested in technical details of the
implementation or who are familiar with computational fluid dynamics can skip
Sub.Sec. IV.3 and IV.5. Sub.Sec. IV.4 mainly addresses readers who work with
the FRG and can be bypassed by everyone else. The other parts of this section
are central for this work. In Sec. V we recapitulate some results for the GN
model in the infinite-$N$ limit. We further discuss initial conditions for the
RG flows at finite $N$ and show consistency checks of our numeric
implementation and choice of UV parameters. Readers familiar with the mean-
field phenomenology and phase structure of the GN model may just skim over
Sub.Sec. V.1 \- V.3. In Sec. VI we discuss the explicit RG flows and the phase
structure of the GNY model at finite $N$. This section comprises our main
results. We conclude our work in Sec. VII and provide an outlook to possible
future research. Several sections are accompanied by appendices, which are
supposed to facilitate a self-contained reading and give computational and
conceptual details omitted in the main part of this paper.

##  II The Gross-Neveu(-Yukawa) model

In this section, we introduce the Gross-Neveu (GN) model, its bosonized
counterpart and the Gross-Neveu-Yukawa (GNY) model in Euclidean space-time. We
comment on its symmetries and its in-medium generalization for non-zero
chemical potentials $\mu$ and temperatures $T$.

###  II.1 In vacuum

The Gross-Neveu model (GN model) in one spatial and one temporal dimension in
Euclidean space-time is defined via the classical action, cf. Ref. [1],

$\displaystyle\mathcal{S}_{\mathrm{GN}}[\bar{\psi},\psi]=\int\mathrm{d}^{2}x\,\big{[}\bar{\psi}\,\not{\partial}\,\psi-\tfrac{g^{2}}{2N}\,(\bar{\psi}\,\psi)^{2}\big{]}\,,$
(1)

where $\psi$ is an $N$-component object in flavor space ($f=1,2,\dots,N$ and
$N>1$)555We explicitly exclude $N=1$, where the GN model is identical to the
Thirring model [36, 37], which has a vanishing perturbative one-loop
$\beta$-function and different phenomenology than the GN model at $N>1$ [9,
10]. Some details and References concerning the Thirring model can be found in
Sec. I. and two-component spinor in Dirac space (see App. A for details on our
conventions including conventions for Euclidean space-time and the connection
to the original formulation [1] of the model in Minkowski space-time). The
action involves a kinetic term and a four-fermion-interaction term with the
dimensionless positive coupling constant $g^{2}$.

Apart from Euclidean space-time symmetries (translations and a rotation), the
action (1) is invariant under transformations of the symmetry group
${G=U(1)\times SU(N)\times\mathbb{Z}_{2}}$. The group acts on the fermion
fields as follows

$\displaystyle G\times\mathbb{C}^{2N}\to\,$
$\displaystyle\mathbb{C}^{2N}\,,\vphantom{\bigg{(}\bigg{)}}$ (2)
$\displaystyle((\alpha,\theta,n),\psi)\mapsto\,$
$\displaystyle\alpha\,(\theta\,\gamma_{\mathrm{R}}+n\,\theta\,\gamma_{\mathrm{L}})\,\psi\,,\vphantom{\bigg{(}\bigg{)}}$

where $n\in\mathbb{Z}_{2}\equiv\\{-1,+1\\}$,
$\gamma_{\mathrm{L}/\mathrm{R}}=\tfrac{1}{2}\,(\openone\mp\gamma_{\mathrm{ch}})$
are the left- and right-handed chiral projection operators, $\alpha\in U(1)$
and $\theta\in SU(N)$. The full group $G$ is defined as the direct product of
the groups $U(1)$, $SU(N)$ and $\mathbb{Z}_{2}$. Hereby, the $SU(N)$ symmetry
is usually called flavor or color symmetry and causes, according to Noether’s
theorem, the conservation of a vector current. The $\mathbb{Z}_{2}$ symmetry
is called discrete chiral symmetry. The $U(1)$ symmetry is called phase
symmetry and leads to a conserved Noether-charge density
$\bar{\psi}\,\gamma^{2}\,\psi$, which is usually called Baryon number density
[96, 97, 98, 82].666It is also worth mentioning, that the $SU(N)$ symmetry
group of the GN model in $1+1$ space-time dimensions is an $O(N)$ symmetry,
which is why it is sometimes also denoted as an $O(N)$ symmetric model [49,
10]. Furthermore the GN model in two dimensions has an additional hidden
$O(2N)$ symmetry between Majorana components of the fermion fields [51] which
prevents the appearance of different four-fermion interaction channels during
renormalization [7].

It can be shown, see App. B, that the partition function of the GN model is
equivalent to the partition function of the bosonized GN (bGN) model,

$\displaystyle\mathcal{S}_{\mathrm{bGN}}[\bar{\psi},\psi,\phi]=\int\mathrm{d}^{2}x\,\big{[}\bar{\psi}\,\big{(}\not{\partial}+\tfrac{h}{\sqrt{N}}\,\phi\big{)}\,\psi+\tfrac{h^{2}}{2g^{2}}\,\phi^{2}\big{]}\,,$
(3)

where $\phi$ is a (real) scalar “auxiliary” or “constraint” field [48, 49, 7,
99]. Being equivalent means that one obtains the same correlation functions as
for the purely fermionic action (1), cf. App. B, the discussion at the
beginning of Sub.Sec. III.4 or the textbooks [10, 9].

In the bosonized version (3), the four-Fermi interaction is replaced by a
Yukawa interaction term with coupling constant $h$ as well as a quadratic
(mass) term $h^{2}/g^{2}$ for the auxiliary field $\phi$. If we postulate

$\displaystyle G\times\mathbb{R}\to\,$
$\displaystyle\mathbb{R}\,,\vphantom{\bigg{(}\bigg{)}}$ (4)
$\displaystyle((\alpha,\theta,n),\phi)\mapsto\,$ $\displaystyle
n\phi\,,\vphantom{\bigg{(}\bigg{)}}$

then the action (3) is invariant under the same symmetry group $G$ as the
original action (1).

Within this work, we are especially interested in the discrete chiral symmetry
transformation, which we understand as the group element
$(\alpha,\theta,n)=(1,\mathbb{I},-1)\in G$, i.e.,

$\displaystyle\psi\mapsto\psi^{\prime}=\gamma_{\mathrm{ch}}\,\psi\,,$
$\displaystyle\phi\mapsto\phi^{\prime}=-\phi\,,\vphantom{\bigg{(}\bigg{)}}$
(5)
$\displaystyle\bar{\psi}\mapsto\bar{\psi}^{\prime}=-\bar{\psi}\,\gamma_{\mathrm{ch}}\,,$

It is this symmetry which prevents the GN model from perturbatively generating
a mass gap, see, e.g. Ref. [1, 10].

By including an additional kinetic term for the bosonic field we obtain the
Gross-Neveu-Yukawa model [8, 10, 15]

$\displaystyle\mathcal{S}[\bar{\psi},\psi,\phi]=\vphantom{\bigg{(}\bigg{)}}$
(6) $\displaystyle=\,$
$\displaystyle\int\mathrm{d}^{2}x\,\big{[}\bar{\psi}\,\big{(}\not{\partial}+\tfrac{h}{\sqrt{N}}\,\phi\big{)}\,\psi-\phi\,(\Box\phi)+\tfrac{h^{2}}{2g^{2}}\,\phi^{2}\big{]}\,,\vphantom{\bigg{(}\bigg{)}}$

which we will use for practical computations in the FRG framework. We
elaborate on the specific model choice and differences between the GN (bGN)
and GNY model in Sub.Sec. III.4.

###  II.2 In medium

In this work, we are mainly interested in the in-medium properties of the GNY
model. In order to work at non-zero Baryon density, we fix the net Baryon
number density by introducing a quark chemical potential $\mu$. Analogously to
standard statistical physics, see, e.g., Refs. [100, 101, 102, 50, 52] for
details, we subtract

$\displaystyle\mu N\equiv\mu\int\mathrm{d}^{2}x\,\bar{\psi}\,\gamma^{2}\,\psi$
(7)

from the classical (UV) action, respectively the probability distribution in
the grand canonical partition function777Note that $\gamma^{2}$ is
$\gamma$-matrix associated with the Euclidean temporal direction in our
conventions.,

$\displaystyle\mathcal{Z}\propto\int[\mathrm{d}\bar{\psi},\mathrm{d}\psi,\mathrm{d}\phi]\,\mathrm{e}^{-\mathcal{S}[\bar{\psi},\psi,\phi]+\mu
N}\,.$ (8)

Furthermore, we introduce non-zero temperature via a compactification of the
Euclidean time direction888Geometrically and technically the introduction of
non-zero temperature via compactification of the Euclidean temporal direction
is identical to defining the entire GNY model on a cylinder of infinite length
and finite circumference $\beta$ right from the beginning. For details, see
App. C., i.e., we replace
$\int\mathrm{d}^{2}x\mapsto\int\mathrm{d}x\int_{0}^{\beta}\mathrm{d}\tau$,
where $T=\tfrac{1}{\beta}$ is the temperature, and implement periodic and
anti-periodic boundary conditions in the Euclidean time direction for bosons
and fermions, respectively. In total, we obtain the in-medium or thermal GNY
model,

$\displaystyle\mathcal{S}[\bar{\psi},\psi,\phi]=\,$
$\displaystyle\int\mathrm{d}x\int_{0}^{\beta}\mathrm{d}\tau\,\big{[}\bar{\psi}\,\big{(}\not{\partial}-\mu\,\gamma^{2}+\tfrac{h}{\sqrt{N}}\,\phi\big{)}\,\psi\vphantom{\Bigg{(}\Bigg{)}}$
(9)
$\displaystyle-\phi\,(\Box\phi)+\tfrac{h^{2}}{2g^{2}}\,\phi^{2}\big{]}\,.\vphantom{\Bigg{(}\Bigg{)}}$

Note that the heat bath, introduced by the non-zero temperature, as well as
the chemical potential explicitly break (Euclidean) Poincaré invariance, i.e.,
invariance under rotations.

##  III The Gross-Neveu-Yukawa model via the Functional Renormalization Group

In this section, we introduce the general concept of the FRG as well as the
Exact Renormalization Group (ERG) equation as our setup of choice. The ERG
equation is used to investigate the research question outlined in the
introduction – the influence of (bosonic) quantum fluctuations on discrete
chiral $\mathbb{Z}_{2}$ symmetry breaking and restoration in the GNY model.
Furthermore, we present and discuss our truncation scheme (ansatz) for the ERG
equation and the resulting RG flow equation used for practical computations in
this work. Finally, we discuss how (dynamic) symmetry breaking realizes itself
during RG flows.

###  III.1 The Functional Renormalization Group

If summarized in a slightly oversimplified manner, the FRG maps the problem of
solving high-dimensional functional integrals like Eq. (8) to the problem of
solving a coupled system of highly non-linear partial integro differential
equations with the classical action $\mathcal{S}(\Phi)$ as initial condition.
For general reviews on the FRG, we refer to Refs. [103, 104, 105, 106, 107,
108, 109, 110, 107, 111].

At the same time, the FRG is much more than a practical tool for the
calculation of macroscopic observables, like expectation values, from
microscopic interactions. The FRG in fact mathematically formalizes the
abstract concept of Kadanoff’s block-spin transformations from small to large
scales (in position space) and Wilson’s corresponding idea of gradually
integrating over momentum shells from high energy scales (UV) to small energy
scales (IR) (in momentum space respectively) [112, 113, 114, 115].

In the modern approach the FRG is formulated in terms of the exact RG
equation999Being “exact” means that its solution for a specific model or
theory is completely equivalent to “solving the functional integral” or to
“knowing all correlation functions of the model or theory”.101010Earlier
attempts to formalize Wilson’s RG approach, using similar (functional)
differential equations, can be found in Refs. [116, 117, 118, 119, 120].
(oftentimes also referred to as Wetterich equation) [121, 122, 123, 124, 125,
126],

$\displaystyle\partial_{t}\bar{\Gamma}_{t}[\Phi]=\,$
$\displaystyle\mathrm{STr}\Big{[}\big{(}\tfrac{1}{2}\,\partial_{t}R_{t}\big{)}\,\big{(}\bar{\Gamma}^{(2)}_{t}[\Phi]+R_{t}\big{)}^{-1}\Big{]}=\vphantom{\Bigg{(}\Bigg{)}}$
(10) $\displaystyle=\,$
$\displaystyle\begin{gathered}\includegraphics{diagrams/wetterich_loop.pdf}\end{gathered}\,.\vphantom{\Bigg{(}\Bigg{)}}$
(12)

Here, $t$ is the RG time (note our sign convention – positive RG time),

$\displaystyle t=-\ln\big{(}\tfrac{k}{\Lambda}\big{)}\,,$ $\displaystyle
t\in[0,\infty)\,,$ (13)

which parameterizes the RG flow and $k(t)$ is the corresponding RG scale with
$\Lambda$ as a UV reference scale.

The ERG equation (10) is formulated in terms of the effective average action
$\bar{\Gamma}_{t}[\Phi]$, which is an RG-scale dependent version of the well-
known effective action $\Gamma[\Phi]$ – the generating functional of
$1$PI-$n$-point correlation functions [67, 127, 128, 129, 77, 9, 10, 130]. It
can be shown [122, 130, 104] that $\bar{\Gamma}_{t}[\Phi]$ approaches the
classical microscopic action $\mathcal{S}[\Phi]$ in the UV at $t=0$
($k=\Lambda$). In the IR, for $t\rightarrow\infty$ ($k\rightarrow 0$), the
effective average action $\bar{\Gamma}_{t}[\Phi]$ coincides with
$\Gamma[\Phi]$. Thus, the ERG equation (10) describes the RG-time evolution of
$\bar{\Gamma}_{t}[\Phi]$ from $\mathcal{S}[\Phi]$ to $\Gamma[\Phi]$ in terms
of a (functional) differential equation, with $\mathcal{S}[\Phi]$ being the
initial condition.

The field content of the theory is summarized in the “superfield” variable
$\Phi$, collecting all classical/mean fields under consideration. Thus the
“supertrace” $\mathrm{STr}$ stands for traces over all internal indices
(color, flavor etc.) as well as a loop integration in momentum/position space
[104]. The ERG equation (10) can be applied to a broad range of QFTs ranging
from gauge-theories, problems in solid-state theory to quantum gravity or
effective QFTs of strong interactions etc.. See, e.g., Refs. [123, 131, 132,
133, 134] for early works based on the ERG equation or Ref. [109] for a recent
review.

The regulator $R_{t}$, that enters into the full scale dependent propagator
(indicated by the solid black line in the Feynman graph) on the r.h.s. of Eq.
(10) acts as a scale- and momentum-dependent IR-mode suppressor. It suppresses
all quantum fluctuations at scales lower than the current RG scale
$k(t)=\Lambda\,\mathrm{e}^{-t}$ and can be interpreted as a scale dependent
mass. Hence, it ensures that only the quantum fluctuations from energy scales
greater than $k(t)$ enter the effective average average action
$\bar{\Gamma}_{t}[\Phi]$ at scale $k(t)$. Lowering $k(t)$, thus evolving in RG
time $t$, integrates out and includes more and more quantum fluctuations in
the effective average action $\bar{\Gamma}_{t}[\Phi]$, realizing an
implementation of Wilson’s concept of the RG. Considering the couplings
(moments of $\bar{\Gamma}_{t}[\Phi]$) we note that in principle all couplings
evolve as effective scale-dependent couplings and new couplings, which are in
accordance with the symmetries of the path integral (partition function), are
generated [117] in the RG flow.

On the other hand, also the scale-derivative of the regulator
$\partial_{t}R_{t}$ – the so called regulator insertion (indicated by the
$\otimes$ symbol in the Feynman graph) plays a crucial role in the ERG
equation (10). For “conventional” regulators it ensures momentum locality
during RG-time evolution, meaning that predominantly fluctuating modes of
similar energy to the energy of the RG scale $k(t)$ enter during an
integration step at the RG scale $k(t)$. Formally $\partial_{t}R_{t}$ peaks
the integrand of the momentum integrals in the supertrace $\mathrm{STr}$
around $p\approx k(t)$, where $p$ is the loop momentum. (For qualitative
sketches of the functional courses of these typical regulators, see, e.g.
Refs. [135, 107, 109].)

Although having a one-loop structure, the ERG equation (10) incorporates
effects from all loop orders from a perturbative point of view, cf. Refs.
[136, 137, 138]. Consequently the FRG is a non-perturbative method, which
allows for its application to systems including strong interactions – even at
non-vanishing chemical potentials111111At least in certain truncations the FRG
does not suffer from obvious complications at $\mu>0$ like the sign-problems
encountered in lattice Monte-Carlo simulations. Computations with involved
truncations for fermionic couplings are however not free of conceptional or at
the least significant practical/numerical challenges at $\mu>0$, see, e.g.,
Refs. [139, 140]. Within this work we encountered practical complications at
low temperatures (including $T=0$) and non-zero quark chemical potentials
$\mu>0$ within the LPA truncation with one-dimensional, optimized Litim
regulators see Sub.Sec. IV.4 and Sec. VI for conceptual considerations and
numerical results. We thank J. Braun and N. Wink for various discussions on
this issue. and temperatures.

However, the ERG equation (10) cannot be used directly for practical
calculations – and there is the rub. For any application so called truncations
(approximations) are needed121212There are some exceptional systems, where the
ERG equation can be solved without any truncations, see e.g. Refs. [141, 142,
110, 143, 144, 145]., which artificially restricts the entire theory space, in
which $\bar{\Gamma}_{t}[\Phi]$ can evolve, to a limited subspace. This is
necessary in order to convert the functional differential equation (10) to a
finite set of ordinary and/or partial differential equations. Choosing
appropriate and reasonable truncations for a specific problem under
consideration, which do not alter the physics and do not ignore significant
effects of the systems under consideration, is a key – but also a very
involved – task, when performing practical computations with the FRG.
Nevertheless, within the last three decades of usage of the ERG equation (10)
specific truncation schemes were established for various applications, see
Ref. [109] and references therein. While the effective (average) action in the
IR ${\bar{\Gamma}_{t\rightarrow\infty}[\Phi]=\Gamma[\Phi]}$ is independent of
the specific choice of regulator $R_{t}$ when considering the untruncated ERG
equation (10), IR results using truncations are no longer guaranteed to be
independent of the regulator choice, see, e.g., Ref. [146]. Checking for
independence of physical results from truncation schemes and regulator choice
is a tedious but in principle necessary task. We come back to this issue at
several points during our work.

We are aware that the FRG (at least in simple truncations like ours, see
below) is usually not the method of choice for quantitative high precision
predictions. But it has several advantages when it comes to the research
questions addressed in this work. The FRG naturally resolves fluctuation
effects – including fermionic and bosonic quantum and thermal fluctuations –
at different energy scales and provides direct access to our observables of
interest (the condensate, the curvature masses and the effective potential) at
all scales $k(t)$. The FRG as a continuum method can be used for direct
computations in infinite volumes, incorporating quantum and thermal
fluctuations non-perturbatively over a huge range of energies (wavelengths).
It therefore complements lattice Monte-Carlo simulations which are based on
finite sized space-time boxes.

Computations in finite volumes are also possible within the FRG framework.
Actually, we plan to repeat the analysis of this work for the GNY in a finite
spatial volume along the lines of Refs. [147, 148, 149, 150, 151] elsewhere,
in order to directly analyze the effects of a finite sized spatial volume and
to compare our results to the ones obtained with lattice Monte-Carlo
simulations [83, 84, 85, 81, 82, 79, 80].

The advantages and shortcomings of our FRG setup will become clear within the
next (sub)sections, where we introduce our truncation scheme and the explicit
RG flow equation, elaborate on the limitations of our approach, and discuss
symmetry restoration/breaking within the FRG approach.

###  III.2 Truncation and Renormalization Group flow equation(s)

In the context of this work, we use the local potential approximation (LPA) as
a truncation for Eq. (10). For the GNY model this means that the effective
average action is approximated as

$\displaystyle\bar{\Gamma}_{t}[\bar{\psi},\psi,\varphi]=\,$
$\displaystyle\int\mathrm{d}x\int_{0}^{\beta}\mathrm{d}\tau\,\big{[}\bar{\psi}\,\big{(}\not{\partial}-\mu\gamma^{2}+\tfrac{h}{\sqrt{N}}\,\varphi\big{)}\,\psi-\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle-\tfrac{1}{2}\,\varphi\,(\Box\varphi)+U(t,\varphi)\big{]}\
.\vphantom{\Bigg{(}\Bigg{)}}$ (14)

In this ansatz exclusively the so-called scale-dependent effective potential
$U(t,\varphi)$ is evolving with RG time $t$ (along RG scale $k(t)$), while all
other couplings (e.g., the Yukawa coupling $h$) are kept constant.

Although being a rather simplistic ansatz, the LPA has turned out as a common
and powerful truncation scheme and is widely – arguably sometimes wildly –
used in the FRG community, especially in the context of strongly interacting
systems, e.g. low-energy effective models, cf. Ref. [109] and references
therein. In the absence of fermions, the LPA can be viewed as the lowest order
contribution of a derivative expansion [152, 153, 103, 109]. The LPA is
assumed to be (and is for certain setups proven to be [141, 142, 154, 110,
143, 144, 155]) a good option for systems that are strongly coupled in field
space and systems, where interactions with low-momentum transfer dominate the
dynamics. For systems, where a high resolution in momentum space at all scales
$k(t)$ becomes relevant, a vertex expansion in powers of fields might be a
better choice, cf. Ref. [109] and references therein. For further discussions
on the quality or the comparison of truncation schemes, see, e.g. Refs. [156,
157, 158, 159, 160] and references therein.

When dealing with the effects of long-range interactions (low momenta) in a
low-dimensional model at non-zero $\mu$ and $T$ (strong coupling and
complicated dynamics in field space), that we are heading at, the LPA presents
as a natural starting point for an analysis beyond the mean-field
approximations.

We are aware, that the inclusion of scale-dependent and potentially even
field-dependent wave function renormalizations and couplings (especially a
scale- and potentially field-dependent Yukawa coupling $h(t,\varphi)$) would
result in a significantly improved truncation. We will come back to this in
Sub.Sec. III.4. For the moment we will just start off with the LPA, which is
already an improvement compared to the commonly used mean-field approximations
for the GN(Y) model, where the effects of bosonic quantum fluctuations are
usually completely ignored or compared to “improved” mean-field
approximations, see e.g., Ref. [50], which include only specific (effects of)
bosonic modes.

The RG flow equation for the effective potential in the LPA is obtained, by
inserting the ansatz for $\bar{\Gamma}_{t}[\bar{\psi},\psi,\varphi]$, Eq.
(14), into the ERG equation (10) followed by a projection onto a suited (here
constant) background field configuration ($\varphi(\tau,x)=\sigma$,
$\bar{\psi}(\tau,x)=0$, $\psi(\tau,x)=0$), cf. Chap. 23 of Ref. [46].
Additionally one has to specify proper, explicit regulators. For a discussion
on suitable choice and influences of different regulators on the RG flow and
the IR results within a truncation, we refer to Refs. [161, 162, 146, 140] and
the discussion in the next Sub.Sec. III.3. In the context of our work, we use
so-called one-dimensional LPA-optimized momentum-space regulators for fermions
and bosons [161, 162], see App. D.

Additionally and analogously to mean-field studies, we perform a rescaling of
the bosonic (background) field and the scale-dependent effective potential,

$\displaystyle\varphi\mapsto\,$
$\displaystyle\tilde{\varphi}=\tfrac{1}{\sqrt{N}}\,\varphi\,,\vphantom{\bigg{(}\bigg{)}}$
(15) $\displaystyle U(t,\varphi)\mapsto\,$
$\displaystyle\tilde{U}(t,\tilde{\varphi})=\tfrac{1}{N}\,U(t,\varphi)\,.\vphantom{\bigg{(}\bigg{)}}$
(16)

This rescaling allows for a comparison of calculations at different finite
values of $N$ and the infinite-$N$ limit. In the following we are exclusively
working in rescaled quantities $\tilde{\varphi}$ and
$\tilde{U}(t,\tilde{\varphi})$, such that we do not maintain the “tilde” in
our notation.

We obtain the RG flow equation for the scale dependent effective potential
$U(t,\sigma)$ in rescaled quantities,

$\displaystyle\partial_{t}U(t,\sigma)=\frac{1}{N}\,\Bigg{(}\begin{gathered}\includegraphics{diagrams/potential_02.pdf}\end{gathered}+\begin{gathered}\includegraphics{diagrams/potential_03.pdf}\end{gathered}\Bigg{)}=\,$
$\displaystyle-\frac{1}{\pi
N}\,\frac{k^{3}(t)}{2E_{\textrm{b}}(t,\partial^{2}_{\sigma}U)}\,\big{[}1+2\,n_{\mathrm{b}}(\beta
E_{\textrm{b}}(t,\partial^{2}_{\sigma}U))\big{]}+\vphantom{\Bigg{(}\Bigg{)}}$
(19)
$\displaystyle+\frac{d_{\gamma}}{\pi}\,\frac{k^{3}(t)}{2E_{\textrm{f}}(t,\sigma)}\,\big{[}1-n_{\mathrm{f}}(\beta[E_{\textrm{f}}(t,\sigma)+\mu])-n_{\mathrm{f}}(\beta[E_{\textrm{f}}(t,\sigma)-\mu])\big{]}\,,\vphantom{\Bigg{(}\Bigg{)}}$

where we introduced the abbreviations

$\displaystyle E_{\mathrm{b}}(t,\partial^{2}_{\sigma}U)\equiv\,$
$\displaystyle\sqrt{k^{2}(t)+\partial_{\sigma}^{2}U(t,\sigma)}\,,\vphantom{\bigg{(}\bigg{)}}$
(20) $\displaystyle E_{\mathrm{f}}(t,\sigma)\equiv\,$
$\displaystyle\sqrt{k^{2}(t)+(h\,\sigma)^{2}}\,,\vphantom{\bigg{(}\bigg{)}}$
(21)

for the Euclidean bosonic and fermionic “energies” (dispersion relations) and
used

$\displaystyle n_{\mathrm{b}}(x)\equiv\frac{1}{\mathrm{e}^{x}-1}\,,$
$\displaystyle n_{\mathrm{f}}(x)\equiv\frac{1}{\mathrm{e}^{x}+1}\,.$ (22)

to denote the Bose-Einstein [163] and Fermi-Dirac [164, 165] distribution
functions. Note the prefactor $\frac{1}{N}$ of the bosonic contributions (the
first term on the r.h.s. of Eq. (19)) realizing the aforementioned suppression
of bosonic fluctuations in the large-$N$ limit.

For the sake of completeness, all further details on the derivation of the
flow equation are provided in App. D.

Although looking complicated at first sight, the RG flow equation (19) for the
effective potential $U(t,\sigma)$ has a rather simple overall structure. It is
a non-linear partial differential equation in only two variables $\sigma$ and
$t$. Similar equations also appear in hydrodynamic problems. In the
hydrodynamic language $\sigma$ plays the role of a spatial coordinate and $t$
is a temporal coordinate (details are provided in Sec. IV). In the following
we adapt this language. The temporal and spatial domain of the PDE should not
be confused with the Euclidean space-time coordinates. The spatial domain of
the PDE is $\mathbb{R}$, hence $\sigma\in(-\infty,+\infty)$, while the
temporal domain is $\mathbb{R}^{+}$, hence $t\in[0,\infty)$. The PDE is first
order in temporal derivatives and second order in spatial derivatives.
Additionally, there are explicit time-($t$)-dependent factors and
position-($\sigma$)-dependent contributions. We will continue the systematic
analysis of this PDE in Sec. IV, when we draw some connections between RG flow
equations and flow equations in their literal sense.

Before we continue our main discussion, we present the zero-temperature limit
of the RG flow equation (19),

$\displaystyle\partial_{t}U(t,\sigma)=\vphantom{\bigg{(}\bigg{)}}$ (23)
$\displaystyle=\,$ $\displaystyle-\frac{1}{\pi
N}\,\frac{k^{3}(t)}{2E_{\mathrm{b}}(t,\sigma)}+\frac{d_{\gamma}}{\pi}\,\frac{k^{3}(t)}{2E_{\mathrm{f}}(t,\sigma)}\,\Theta\Big{(}1-\tfrac{\mu^{2}}{E^{2}_{\psi}(t,\sigma)}\Big{)}\vphantom{\bigg{(}\bigg{)}}$

as well as the vacuum limit for $T\rightarrow 0$ and $\mu\rightarrow 0$,

$\displaystyle\partial_{t}U(t,\sigma)=\,$ $\displaystyle-\frac{1}{\pi
N}\,\frac{k^{3}(t)}{2E_{\mathrm{b}}(t,\sigma)}+\frac{d_{\gamma}}{\pi}\,\frac{k^{3}(t)}{2E_{\mathrm{f}}(t,\sigma)}\,.$
(24)

The RG flow equation in vacuum is needed to fix the initial condition for all
RG flows – with and without the effects of a medium.

Note that Eq. (24) differs from the popular LPA flow equation of the effective
potential in vacuum, cf. Ref. [12, 15],

$\displaystyle\partial_{t}U(t,\sigma)=\,$ $\displaystyle-\frac{1}{4\pi
N}\,\frac{k^{4}(t)}{2E_{\mathrm{b}}^{2}(t,\sigma)}+\frac{d_{\gamma}}{4\pi}\,\frac{k^{4}(t)}{2E_{\mathrm{f}}^{2}(t,\sigma)}\,,$
(25)

due to the fact that we are using one-dimensional purely spatial LPA-optimized
regulators in Eq. (24) that do not regulate the Matsubara summation (before
taking the limit $T\rightarrow 0$), instead of using two-dimensional LPA-
optimized regulators employed in the derivation of Eq. (25). We comment on
this issue and its consequences for our work in the next Sub.Sec. III.3,
Sub.Sec. VI.6, and App. E.

###  III.3 Comment on regulators

At this point we ought to comment on the choice of our regulators and the
caveats that go hand in hand with this choice. The regulator shape functions
(120) and (121) that specify the regulators (112) and (113) only regulate the
spatial momentum (direction). Hence, if the spatial (loop) momenta are smaller
than the RG scale $k$, they are suppressed. However, the Matsubara sums are
not affected by the regulator at all and frequencies of all orders of
magnitude enter the RG flow at all scales $k(t)$. This has several direct
consequences. Although regulators should in principle be in accordance with
the symmetries of the theory, the regulators (112) and (113) explicitly break
(Euclidean) Poincaré invariance. Hence, we cannot expect to recover full
(Euclidean) Poincaré symmetry in the limits $\mu\rightarrow 0$ and
$T\rightarrow 0$, without the introduction of counter-terms (Ward identities)
to account for this discrepancy, e.g., Refs. [166, 167, 168]. This means that
the IR results of the RG-flow equation (24) do not necessarily coincide with
IR results of the vacuum LPA flow equation (25) that can be derived with two-
dimensional LPA optimized regulators, if one uses exactly the same UV initial
condition. Whether or not these differences manifest in physical observables
in vacuum or even medium, since computations at non-zero $T$ and $\mu$ usually
use an UV initial condition fixated in vacuum, depends on the observable,
model, and truncation under consideration, see, e.g., Refs. [166, 167, 168].
Within the scope of this work we performed some test comparing results
obtained with a two-dimensional LPA-optimized regulator to results obtained
using the one-dimensional LPA-optimized regulator in vacuum using identical
initial conditions. For a brief discussion see App. E. The situation for the
GN(Y) model might be discussed elsewhere in more detail especially regarding
regulator dependencies at non-zero temperature [169].

One might ask now, why we are – regardless of these facts – using one-
dimensional LPA-optimized regulators? The answer to this question has several
aspects:

A first drawback of using two-dimensional regulators is that large classes of
those regulators cause problems in the presence of chemical potentials and
violate the so called Silver-Blaze property [170, 171, 94]. However, we are
especially interested in calculations at non-zero $\mu$. Coping with this
challenge is part of state of the art research, see, e.g., Refs. [166, 172,
140], and we do not want to enter this discussion within this work.

Secondly, the analytic evaluation of the Matsubara sums or the loop-momentum
integrals in Eq. (119) might become impossible or at least extremely
challenging, which drastically complicates numerical computations. The
presence of numerical sums and integrals in the flow equation would
significantly increase computation time and hinder an in depth discussion at
variable $T$, $\mu$ and $N$.

Note that for non-zero temperature and chemical potential (Euclidean) Poincaré
invariance is broken anyhow, such that explicitly breaking this symmetry via
the regulators might not spoil the results too drastically.

Mainly to facilitate and speed up numerical computations and to avoid any
conceptual issues at $\mu>0$ we decided to use one-dimensional, spatial Litim
regulators within this work as a first significant step beyond mean-field
computations.

Interestingly, the approach of exclusively regulating spatial momenta is
similar to the common strategy employed in conventional mean-field studies at
non-zero temperature including the mean-field computations for the GN(Y) model
[47, 48, 49, 50, 51, 52]. There Matsubara summations are usually executed
analytically before momentum integrals are regulated at all. Divergent
contributions (usually associated with vacuum quantum fluctuations) to
expressions are separated from convergent (usually thermal) contributions and
only the divergent parts are regulated. Both approaches include all
thermal/Matsubara modes independent of RG scale or in case of mean-field
studies of the chosen regularization scheme.

###  III.4 Truncations and initial conditions – Differences and similarities
of GN, bGN and GNY models

In Sec. II we introduced three distinct models by specifying their actions:
the GN model defined by $\mathcal{S}_{\mathrm{GN}}$ in Eq. (1), the bGN model
defined by $\mathcal{S}_{\mathrm{bGN}}$ in Eq. (3) and ultimately the GNY
model defined by $\mathcal{S}$ in Eq. (6). In this subsection we will
elaborate on their differences and similarities in the FRG framework with
special focus on the in this context central issues of truncations and initial
conditions for their respective effective average actions
$\bar{\Gamma}_{t}[\bar{\psi},\psi,\varphi]$. As for all RG flows (and PDEs in
general) their respective solutions depend on the corresponding initial and
possible boundary condition(s), see, e.g., Refs. [173, 174]. We discuss the
(numerical) boundary conditions in Sub.Sec. IV.3 and solely focus on the
initial condition for the moment.

The GN model and its bosonized version – the bGN model – have an identical
partition function and are thus physically equivalent. The bosonization
procedure of App. B is a merely technical reformulation with the goal to
eliminate the four-Fermi coupling term $\propto(\bar{\psi}\,\psi)^{2}$ of the
GN model by introducing the auxiliary field $\phi$ with a mass term
$\propto\phi^{2}$ and a Yukawa-Coupling term $\propto\phi\,\bar{\psi}\,\psi$.
This reformulation facilitates computations especially for non-vanishing
condensates $\langle\phi\rangle\propto\langle\bar{\psi}\,\psi\rangle$, cf.
Ref. [1, 10, 9, 175].

As stated already in the introduction to the FRG method in Sub.Sec. III.1, the
initial condition for $\bar{\Gamma}_{t}[\Phi]$ for the ERG equation (10) is
the classical action $\mathcal{S}[\Phi]$. Within the LPA truncation the only
scale-dependent quantity is the potential $U(t,\sigma)$ for which an
appropriate initial condition at $k(t=0)=\Lambda$ can be read off directly
form the classical action,

$\displaystyle U(0,\sigma)=\,$
$\displaystyle\tfrac{h^{2}}{2g^{2}}\,\sigma^{2}\,,$ (26)

while the Yukawa coupling $h$ keeps its initial value throughout the RG flow
in the LPA. (Note the two rescalings in Eq. (15), which cancel for the
rescaled initial potential (26)). This initial condition $U(0,\sigma)$ is
valid both for the bGN and GNY model and will be discussed further in Sub.Sec.
V.4.

The only difference between the bGN (and by proxy the GN) model and the GNY
model is the kinetic term $-\phi\,(\Box\phi)$ in the GNY model. A
corresponding contribution $-\varphi\,(\Box\varphi)$ (in terms of
classical/mean fields $\varphi=\langle\phi\rangle$) is of course also present
in the LPA ansatz (14) for $\bar{\Gamma}_{t}[\bar{\psi},\psi,\varphi]$ of the
GNY model to establish the proper initial condition
$\bar{\Gamma}_{t=0}=\mathcal{S}$. This kinetic term presents an obvious
mismatch between the LPA formulation of the GNY model and a possible LPA
formulation of the bGN model. Nevertheless, such a kinetic term is needed in
the FRG framework to study the effect of bosonic quantum fluctuations of
$\phi$ at finite $N$. In consequence, the LPA of the GNY model for
$\bar{\Gamma}_{t}[\bar{\psi},\psi,\varphi]$ seems not capable of resembling
the bGN model (GN model by proxy) in the UV. One might argue that the results,
which are obtained from the RG flow equation (19) with initial condition (26),
are consequently not directly transferable to the GN model. It seems, as if
bosonic quantum fluctuations, which are linked to the kinetic term, are
artificially enhanced already at the beginning of the RG flow, while they
should actually be strongly suppressed (strictly speaking vanishing directly
in the UV when considering the limit $\Lambda\rightarrow\infty$) in the bGN
model (3), which starts without the kinetic term in the UV and generates the
term and bosonic fluctuations dynamically during the RG flow.

From the FRG perspective, the intuitive way to cure this problem is the
introduction of a bosonic wave function renormalization in the kinetic term,

$\displaystyle-\varphi\,(\Box\varphi)\rightarrow-
Z_{\varphi}(t)\,\varphi\,(\Box\varphi)\,,$ (27)

and initializing ${1\ggg Z_{\varphi}(0)>0}$ in the UV131313Initializing
$Z_{\varphi}(0)$ at exactly zero in the UV leads to complications in practical
computation since, e.g., the renormalized mass term
$\tfrac{1}{Z_{\varphi}(t)}\,\partial_{\sigma}^{2}U(t,\sigma)$ in the
dispersion relation (20) would diverge. Besides technical problems
$Z_{\varphi}(0)=0$ at a finite UV initial scale $\Lambda$ would weakly violate
RG-consistency [176] since fluctuations (especially fermionic ones) at scales
$k>\Lambda$ would have already generated a small wave function
renormalization. to make direct contact with the action (3) of the bGN model.
Initializing the wave function renormalization with ${1\ggg Z_{\varphi}(0)>0}$
in the UV at a sufficiently large UV initial scale $\Lambda$ – large enough at
the initial scale and in doing so realizing RG consistency [176] – leads to a
direct suppression of bosonic fluctuations at the beginning of the RG flow.
This is directly seen on the level of
$\bar{\Gamma}_{t}[\bar{\psi},\psi,\varphi]$, where the bosonic kinetic term is
suppressed by its small prefactor $Z_{\varphi}(t=0)$ in the UV. Most likely,
during the RG flow $Z_{\varphi}(t)$ will turn non-zero already by its
fermionic loop contribution, cf. Ref. [12, 99], and we expect to recover the
situation, where $Z_{\varphi}(t)$ is non-zero and facilitates bosonic
fluctuations in the IR.

Hence, we already conclude at this point that a natural generalization of the
LPA truncation in this publication is the inclusion of a bosonic wave function
renormalization, such that the UV initial condition for
$\bar{\Gamma}_{t}[\bar{\psi},\psi,\varphi]$ in the GNY model really resembles
the bGN and GN. This will be discussed elsewhere. (Of course, the additional
inclusion of further scale- and potentially field-dependent couplings and wave
function renormalizations, e.g., $h(t,\varphi)$, $Z_{\varphi}(t,\varphi)$,
$Z_{\psi}(t,\varphi)$, etc., is the natural step beyond the inclusion of just
$Z_{\varphi}(t)$.)

However, for the sake of this work, we note that the effect of a strong (or
even total) suppression of bosonic fluctuations in the UV (at fixed $N$) is
present already in the LPA truncation anyhow due to peculiarities of the RG
flow equation and the quadratic shape of the UV initial potential (26). The
first effect relates to the choice of $h$ and $g^{2}$ at $k(t=0)=\Lambda$. We
know that the GN model is asymptotically free, meaning
$\tfrac{1}{g^{2}}\sim\ln\big{(}\tfrac{\Lambda}{h}\big{)}$. Hence, for large
$\Lambda$ and finite $h$ the bosonic mass stemming from Eq. (26) is large
compared to the fermionic mass and bosonic fluctuations are expected to be of
minor importance – even if the bosonic wave function renormalization is fixed
$Z_{\mathrm{\varphi}}(t)=1$. The second effect relates to the quadratic shape
of the initial potential (26) and the (diffusive) character of the bosonic
contribution to the RG flow equation. A purely quadratic self-interaction
potential is the manifestation of a free, non-interacting bosonic sector,
which decouples completely from the fermionic one and has no dynamics during
RG flow. This holds true for the bosonic sector as long as the fermionic
contributions do not significantly alter the self-interaction potential away
from its initial quadratic shape. Both effects – the large mass and the
retention of a quadratic self-interaction potential in the UV – are discussed
in detail in Sub.Sec. V.4, after we firstly introduced the fluid-dynamic
reformulation of the flow equation (19) in Sec. IV and secondly discussed
asymptotic freedom and RG consistency in the context of the infinite-$N$ limit
(mean-field) in Sec. (V).

For the moment, we note that the LPA seems to turn out as a decent
approximation to capture the qualitative dynamics, without spoiling the actual
RG flow too much. This statement is refined in Sub.Sec. V.4 and supported by
our explicit numerical results in Sec. VI, where it is visible that the
bosonic contributions to the RG flow are relevant at much later RG times
(lower scales) than those of the fermions. Nevertheless, we plan to study
higher order truncation effects elsewhere.

Finally, we have to briefly comment on another issue, which also stems from
the truncation scheme and the mismatch of the initial conditions for the bGN
and GNY models. The original GN action (1) as well as its bosonized
counterpart in Eq. (3) without the kinetic term $-\phi\,(\Box\phi)$ have only
a single coupling constant $g^{2}$. The artificial Yukawa coupling $h$
introduced during the bosonization procedure can actually be absorbed in the
bosonic field. In mean-field calculations, see below, this is evident, because
all observables can be fixed via a single dimensionful parameter (e.g., the IR
bosonic curvature mass, the IR fermion mass, the critical temperature etc.)
for $\Lambda\rightarrow\infty$ and can be mapped into each other, if different
renormalization schemes were chosen. Including bosonic fluctuations, i.e. the
kinetic term for the bosons, the Yukawa coupling $h$ can no longer be
completely absorbed via appropriate substitutions. Hence, the model involves
an additional parameter that needs to be fixed via some renormalization
condition [10] – even if we stick to quadratic initial potentials like (26).

We also come back to this aspect in Sec. V, where we motivate our explicit
initial conditions for the RG flows. Though, we have to recapitulate some
mean-field results in advance.

###  III.5 Dynamical symmetry breaking and RG flows

Since the FRG is formulated in terms of the effective (average) action, we set
the following criteria for (spontaneous) symmetry breaking.

The vacuum of a QFT is defined as the field configuration with least energy,
which has to be a field configuration that minimizes the IR effective action
$\Gamma[\Phi]$, cf. Refs. [130, 77]. Hence it has to be a solution to the
quantum equations of motion for the mean fields $\bar{\psi}$, $\psi$, and
$\varphi$,

$\displaystyle
0\overset{!}{=}\frac{\delta\Gamma[\bar{\psi},\psi,\varphi]}{\delta\bar{\psi}}\,,$
$\displaystyle
0\overset{!}{=}\frac{\delta\Gamma[\bar{\psi},\psi,\varphi]}{\delta\varphi}\,.$
(28)

Within our truncation (14), these equations reduce to the following two
classical PDEs

$\displaystyle 0=\,$
$\displaystyle(\not{\partial}-\mu\gamma^{2}+h\,\varphi)\,\psi\,,\vphantom{\Bigg{(}\Bigg{)}}$
(29) $\displaystyle 0=\,$
$\displaystyle\Box\,\varphi-\partial_{\varphi}U(t_{\mathrm{IR}},\varphi)-\tfrac{h}{N}\,\bar{\psi}\,\psi\,,\vphantom{\Bigg{(}\Bigg{)}}$
(30)

where we did not separately list the corresponding equations for the
associated spinor $\bar{\psi}$.

At this point, it is natural to consider vanishing fermion mean fields
$\bar{\psi}(\tau,x)=0$, $\psi(\tau,x)=0$. Furthermore, if we assume the field
configuration to be invariant under any kind of space-time translations, thus
$\varphi(\tau,x)=\mathrm{const}.\equiv\sigma$, we arrive at

$\displaystyle 0\overset{!}{=}\partial_{\sigma}U(t_{\mathrm{IR}},\sigma)\,.$
(31)

In order to ensure that this extremum condition defines a minimum and not a
maximum, we need the additional sufficient condition
$\partial_{\sigma}^{2}U(t_{\mathrm{IR}},\sigma)>0$ at the extremum, viz. a
positive IR curvature mass for $\sigma$.

Consequently, in the FRG framework a minimum $\sigma_{\mathrm{min}}$ of the
scale dependent effective potential $U(t,\sigma)$ in the IR (for
$t\rightarrow\infty$) is considered to be a (translational invariant) ground
state of our system. The system is said to be in the symmetry broken phase, if
$\sigma_{\mathrm{min}}\neq 0$, and in the symmetric phase if
$\sigma_{\mathrm{min}}=0$, because non-trivial $\sigma_{\mathrm{min}}$ break
the $\mathbb{Z}_{2}$ symmetry of the vacuum.

The breaking or restoration of the $\mathbb{Z}_{2}$ symmetry takes place
dynamically. By integrating out quantum fluctuations from short wavelength
(UV) to long wavelength (IR), the fermion interactions may dynamically form a
condensate, which manifests as a non-trivial minimum
$\sigma_{\mathrm{min}}(t)$ of $U(t,\sigma)$ during the RG-flow at finite $t$.
On the other hand – as discussed in the introduction – there are arguments
[71, 72, 50, 14, 10] for the dynamical restoration of $\mathbb{Z}_{2}$
symmetry in the IR, which can be interpreted as a realization of the
precondensation phenomenon found by Refs. [89, 90, 91, 92, 93, 94] in the
context of RG flows. The underlying dynamics can be attributed to long-range
bosonic fluctuations in the FRG framework. Within this work, we can explicitly
resolve this complicated dynamics during the RG flows, compare Sec. VI.

We are aware, that there might be field configurations, which are not
invariant under translations, solve the quantum equations of motion (29) &
(30), and which might have an even smaller energy than the homogeneous
solution $\sigma_{\mathrm{min}}$. However, these so called “(spatially)
inhomogeneous (chiral) condensates”141414For further literature on
“(spatially) inhomogeneous (chiral) condensation” we refer to the Review [59]
and references therein and to the comprehensive work of M. Thies and
collaborators in the context of $(1+1)$-dimensional QFTs, see Tab. 1 of Ref.
[177] and references therein. are not part of this work and might be
considered elsewhere.

##  IV The FRG and (numerical) fluid dynamics

This section is dedicated to the formalities of the (numerical) solution of
the RG flow equation (19). Our discussion is twofold. Firstly, we discuss
analogies of the RG flow equation of the effective potential (19) and common
PDEs from the field of (numerical) fluid dynamics. This leads to an intuitive
understanding of the dynamics of the scale dependent effective potential
$U(t,\sigma)$ (or rather its derivative
$u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$) during the RG flow in terms of
fluid dynamic notions like diffusion or sink/source terms. Secondly, we
introduce the explicit implementation of our numerical scheme, which
originally stems from the research field of computational fluid dynamics.

###  IV.1 The LPA flow equation as a non-linear heat equation with sink/source
term

Within the last decades a huge variety of different numerical schemes was used
to solve field-dependent RG flow equations of similar type than Eq. (19).
However, only recently it was found by some of the authors and their
collaborators [142, 154, 110, 143, 144, 145] that the RG flow equation for the
scale dependent effective potential for a large class of models from QFT can
be recast as a conservation law in terms of a (advection-)diffusion-
source/sink equation and be entirely understood in terms of fluid dynamic
notions151515Earlier works, which already revealed some analogies between RG
flows and fluid dynamic PDEs are Refs. [178, 106, 179, 180, 120].. A fluid
dynamic reinterpretation is as advantageous for several reasons. On the one
hand, the different contributions to the RG flow can easily be understood via
their fluid dynamic analogues. On the other hand, this formulation provides
access to the highly developed toolbox of computational fluid dynamics [173,
174, 181, 182].

Hence, the starting point of our discussion is the conversion of the PDE (19)
into a non-linear diffusion equation. To this end, we act with
$\frac{\mathrm{d}}{\mathrm{d}\sigma}$ on both sides of Eq. (19) and obtain a
flow equation in conservative form161616Some authors, see, e.g. Refs. [183,
184, 185, 186], already studied flow equations of this form in a weak
formulation without the explicit identification and interpretation as a
conservation law.

$\displaystyle\partial_{t}u=\tfrac{\mathrm{d}}{\mathrm{d}\sigma}\,Q(t,\partial_{\sigma}u)+S(t,\sigma)\,,$
(32)

where $u=u(t,\sigma)\equiv\partial_{\sigma}U(t,\sigma)$ acts as the conserved
quantity. Henceforth, we refer to

$\displaystyle Q(t,\partial_{\sigma}u)\equiv\vphantom{\bigg{(}\bigg{)}}$ (33)
$\displaystyle\equiv\,$ $\displaystyle-\frac{1}{\pi
N}\,\frac{k^{3}(t)}{2E_{\textrm{b}}(t,\partial_{\sigma}u)}\,\big{[}1+2\,n_{\mathrm{b}}(\beta
E_{\textrm{b}}(t,\partial_{\sigma}u))\big{]}\,,\vphantom{\bigg{(}\bigg{)}}$

as a non-linear diffusion flux and

$\displaystyle S(t,\sigma)\equiv\vphantom{\bigg{(}\bigg{)}}$ (34)
$\displaystyle\equiv\,$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\sigma}\bigg{(}\frac{d_{\gamma}}{\pi}\,\frac{k^{3}(t)}{2E_{\textrm{f}}(t,\sigma)}\,\big{[}1-n_{\mathrm{f}}(\beta[E_{\textrm{f}}(t,\sigma)+\mu])-\vphantom{\bigg{(}\bigg{)}}$
$\displaystyle-
n_{\mathrm{f}}(\beta[E_{\textrm{f}}(t,\sigma)-\mu])\big{]}\bigg{)}\,,\vphantom{\bigg{(}\bigg{)}}$

as a local, non-linear, internal sink/source term, cf. Ref. [182]. On the
level of the PDE, especially from a fluid dynamic point of view, the RG time
$t\in[0,\infty)$ naturally corresponds to an effective temporal variable
(parameter), while the field space, spanned by $\sigma\in(-\infty,\infty)$,
presents as an effective (infinite) spatial domain. The derivative of the
effective potential $u(t,\sigma)$ plays the role of a fluid. The
interpretation of $Q(t,\partial_{\sigma}u)$ as a non-linear diffusion flux
becomes apparent when the derivative in Eq. (32) is performed,

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\sigma}\,Q(t,\partial_{\sigma}u)=D(t,\partial_{\sigma}u)\,\partial_{\sigma}^{2}u(t,\sigma)\,.$
(35)

Without the source term, we find that the problem of solving the RG flow
equation of the effective potential (19) has reduced to solving a heat
equation, cf. Ref. [187], with a time- and gradient-dependent highly non-
linear diffusion coefficient $D(t,\partial_{\sigma}u)$.

The interpretation of the purely fermionic contribution $S(t,\sigma)$ in terms
of a time dependent source/sink term is also rather natural. From Eq. (34) we
find that the fluid $u$ does not explicitly enter $S(t,\sigma)$ and the
contribution to the flow of the fluid $u$ solely depends on the time $t$ and
position in the spatial direction $\sigma$. Actually it presents on a formal
level in our flow equation as a classical potential (similar to a Newtonian
gravitational potential) of a conservative (external) force field acting on
the fluid, cf. Ref. [182]. Executing the $\sigma$-derivative in Eq. (34) one
finds that the term enters the flow equation (32) for $u$ with negative sign
for $\mu=0$, which explains the denotation as a “sink” for $\sigma>0$ and
“source” for $\sigma<0$ and is responsible for symmetry breaking during the
flow at $\mu=0$171717Note that the dynamics of $u(t,\sigma)$ is anti-symmetric
in $\sigma$. Due to this anti-symmetry and w.l.o.g. see Sub.Sec. IV.3, we
focus on $\sigma\in[0,\infty)$, where the “sink” interpretation is more
natural. For negative $\sigma$ the dynamics has of course opposite sign.. The
bosonic contributions in terms of non-linear diffusion tend to fill the “sink”
or distribute the “source” that is caused by the fermions and work against
symmetry breaking. At non-zero $\mu$ the sink/source term (34) is not always
manifest negative/positive: it has a rather intricate dynamic (especially at
low and zero temperature) and manifests as either source or sink depending on
the value for $\sigma$ at a given $\mu$ and $t$ for both – negative and
positive $\sigma$. This will be discussed at length in Sub.Sec. IV.4.

###  IV.2 Irreversibility and entropy production

Another direct consequence of the reinterpretation of the RG flow equation
(19) in terms of a diffusion-type equation (32) is that the irreversible
character of the RG as a semi group is directly hard-coded and understood on
the level of the PDE: Diffusion is a dissipative and irreversible process,
that introduces a “thermodynamic arrow of time” [188] into the problem. This
arrow of time clearly singles out the RG time (scale) as a temporal
coordinate. Additionally, the dissipative character in field space (space of
couplings) can be seen as a direct realization of Kadanoff’s irreversible
block spin-transformations (in position space) [120, 178, 115]. It follows
that the reverse process of integrating from the IR towards the UV is
generically impossible – also directly on the level of the PDE.

Though, the mean-field truncation (infinite-$N$ limit) seems to violate this
fact, especially concerning the discussion on asymptotic freedom, see Sec. V.
However, this has to be understood as an artifact of the oversimplified
restrictions on the evolution of $\bar{\Gamma}_{t}[\bar{\psi},\psi,\varphi]$
in theory space. For mean-field approximations the theory-space is entirely
described via a single field-independent coupling constant $g^{2}$, which is
evolving with $t$. Beyond the infinite-$N$ limit, this is an
oversimplification of the real dynamics – especially in the IR. However, via
the inclusion of bosonic fluctuations, the RG flow can retain its correct
irreversible diffusive character. Nevertheless, this is only possible if we do
not restrict the bosonic effective potential to a finite set of couplings.
Instead we have to resolve the dissipative dynamics in field space correctly –
also by using adequate numerical methods. For further discussions on these and
related issues, we refer to Refs. [120, 178, 110, 143, 144, 145].

###  IV.3 Field space boundary conditions

In order to make the solution of a PDE a well-posed problem, one needs to
specify initial and boundary conditions. Formally, our PDE (32) presents as a
pure initial value problem on $\mathbb{R}$, because the initial condition for
the PDE (the UV potential) $u(0,\sigma)$ is specified for
$\sigma\in(-\infty,\infty)$, cf. Eq. (26). Thus, on a formal mathematical
level, boundary conditions are not needed, because the problem is completely
defined via the asymptotics of the flow equation (32) and the initial
condition, the $\sigma$-derivative of Eq. (26), for
$|\sigma|\rightarrow\infty$. Asymptotically, the $\sigma$-derivative of the
diffusion flux (33) vanishes. The same holds true for the sink term, which is
suppressed for large $|\sigma|$. In total $u(t,\sigma)$ does not change for
$|\sigma|\rightarrow\infty$ and keeps its initial shape. But still, there is a
lot of dynamics, which however takes place only within a small subdomain of
the field space. Additionally, direct practical calculations are impossible on
an infinite interval. To facilitate practical (numerical) computations, we
artificially introduce boundary conditions to the problem. These boundary
conditions have to be constructed in a way, that does not spoil the original
dynamics of the system.

In a first step, we take advantage of the $\mathbb{Z}_{2}$ symmetry
$U(t,\sigma)=U(t,-\sigma)$, which implies a $\mathbb{Z}_{2}$ anti-symmetry of
$u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$, i.e., $u(t,\sigma)=-u(t,-\sigma)$
and therefore restrict our analysis to $\sigma\in[0,\infty)$ w.l.o.g.. This
amounts to the introduction of anti-symmetric reflective boundary conditions
at $\sigma=0$ on the level of the PDE (32). On the other hand, we argued that
the diffusion flux and sink-contribution fall off drastically for large
$|\sigma|$. Therefore, it is feasible to truncate the spatial domain at some
finite but large $\sigma_{\mathrm{max}}$, beyond which the $t$-evolution
$u(t,\sigma)$ is negligible. At this second artificial spatial boundary, we
can use linear extrapolating boundary conditions for $u(t,\sigma)$ due to the
linear asymptotics of the initial condition for $u(t,\sigma)$, the derivative
of the UV potential (26), which does not change during the flow. Obviously, we
have to test, which numerical values for $\sigma_{\mathrm{max}}$ are
sufficiently large. This is done in App. F.1.

For a detailed discussion on field space boundary conditions for RG flow
equations (and their explicit numeric implementation), we refer to Refs. [110,
144] and Sub.Sec. IV.5.

###  IV.4 The chemical potential as a “shock wave” in field space

Formulating the RG flow equation on the level of the derivative of the scale-
dependent effective potential $u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$ also
allows for a better understanding of the dynamics induced by a non-vanishing
chemical potential $\mu$. It turns out that the chemical potential induces a
kind of shock wave in $u(t,\sigma)$ in field space. This is understood as
follows.

Figure 1:  Source/sink term of Eq. (34) with Yukawa coupling $h=1$ at zero
(solid lines) and small temperature $T=0.01$ (dashed lines) at various RG
scales $k$ (close to and below the model scales) between $k=1.25$ in blue and
$k=0$ in red at $\mu=0$ (upper panel), $\mu=0.5$ (middle panel) and $\mu=1.0$
(lower panel).

For the sake of simplicity, we analyze the RG flow equation (19) at vanishing
(and small) temperature, hence Eqs. (23) and (24), but on the level of
$u(t,\sigma)$. We find that the difference between the $T=0$ flow equations
for $u(t,\sigma)$ at vanishing and non-vanishing chemical potential is the
contributions stemming from the derivative of the Heaviside function
$\Theta(1-\mu^{2}/E_{\mathrm{f}}^{2}(t,\sigma))$ in the source/sink term. The
sink/source term at different RG scales and chemical potentials is plotted in
Fig. 1 for $T=0$ and also small $T=0.01$ to support and illustrate the
following discussion. For small $T$ the Heaviside function gets smeared out.
We will discuss the dynamics induced by the sink/source term in the following
paragraphs. The approximate signs in the rest of this subsection hold at small
temperatures and become exact for $T=0$.

As long as $E^{2}_{\psi}\gtrsim\mu^{2}$, effects due the chemical potential
are not visible in the RG flow, because $\mu$ does not show up in the
source/sink term, see Fig. 1. However, when
$E_{\mathrm{f}}^{2}\approx\mu^{2}$, the chemical potential becomes relevant
and distinguishes RG flows with $\mu=0$ and $\mu\neq 0$. Therefore, it is
important to understand when (in terms of RG time $t$), where (in field space
$\sigma$), and how this is going to happen. This can be done by analyzing the
estimate

$\displaystyle E_{\mathrm{f}}^{2}=k^{2}+(h\sigma)^{2}\approx\mu^{2}$ (36)

for different positions in field space.

The first time during the RG flow, when $\mu$ is going to influence the
dynamics is, when $k^{2}(t)\approx\mu^{2}$, cf. Fig. 1 at $k=0.99$ and
$k=0.49$ in the middle and lower panel. At this RG time (scale), we find that
the chemical potential will influence the dynamics at positions in field space
close to $\sigma=0$. At slightly later RG times, when $t$ is larger and the RG
scale $k(t)$ is lowered a little, also field space positions at slightly
larger $|\sigma|$ will be influenced by $\mu$. The largest $|\sigma|$ that are
directly affected by $\mu$ in the PDE are of the order
$h^{2}\sigma^{2}\approx\mu^{2}$. This is the case when $t\rightarrow\infty$
and $k(t)\rightarrow 0$ in Eq. (36), cf. Fig. 1 at $k=0.0$ in the middle and
lower panel. The peaks induced by the chemical potential are located at
${\sigma\approx\pm\tfrac{1}{h}\sqrt{\mu^{2}-k(t)^{2}}}$ for $\mu\gtrsim k(t)$
and change the character of the sink/source term around their respective
position. For $\sigma>0$ ($\sigma<0$) the peak presents as a source (sink).

Of course, the diffusive character of the $\sigma$-loop contribution will
transport the effects of and the information about the chemical potential via
diffusion also to larger $|\sigma|$. The same holds true for the non-zero
temperature version of the RG flow equation, where the Heaviside function is
smeared out in terms of a Fermi-Dirac distribution function (22). However,
overall the RG time $t$ and scale $k(t)$ as well as the relevant positions in
field space $\sigma$ for the dynamics induced by the chemical potential are
very similar at low temperatures where the dynamics is dominated by the
chemical potential, cf. Fig. 1 for $T=0.01$.

The remaining question, that needs to be answered is, how the chemical
potential $\mu$ influences the dynamics of the RG flow.

We have argued before, that the fermionic contributions to the RG flow of
$u(t,\sigma)$ enter as local time-dependent sink terms. Now we found that the
sinking stops at certain positions in field space and certain RG times
suddenly during the RG flow due to the chemical potential. Hence, we expect
that the chemical potential introduces a sharp edge in $u(t,\sigma)$ at
$\sigma\approx 0$, when $k^{2}(t)\approx\mu^{2}$. This edge will “move” in
field space like a shock wave towards larger $|\sigma|$ until it reaches
$|\sigma|\approx|\mu/h|$, where it comes to a halt. This dynamics is imprinted
by the underlying dynamics of the sink/source term (34) visualized and
discussed earlier.

As already said before, this discontinuity is smeared by the diffusion an
washed out right from the beginning at non-zero temperatures. However, its
dynamics is clearly underlying all diffusive processes. In order to visualize
the drastic effects of this discontinuity, we plot and contrast mean-field RG
flows with RG flows involving bosonic fluctuations in Sub.Sec. VI.2. In mean-
field, only the fermions are active (the sink term), such that the induced
shock wave is not subject to diffusion effects and thus clearly visible.

In fact, it is the chemical potential that introduces a jump (mean-field) or
cusp (finite $N$) discontinuity in $u(t,\sigma)$, which seems to render
practical calculations at zero temperature involving bosonic fluctuations of
the $\sigma$-mode impossible: The discontinuity in $u(t,\sigma)$ corresponds
to a large jump in the derivative $\partial_{\sigma}u(t,\sigma)$ with negative
sign. However, this happens at rather small $k^{2}(t)\approx\mu^{2}$, such
that the bosonic energy function,

$\displaystyle
E_{\mathrm{b}}^{2}(t,\partial_{\sigma}u(t,\sigma))=k^{2}(t)+\partial_{\sigma}u(t,\sigma)\,,$
(37)

turns negative and one overshoots the pole of the bosonic propagator
$\tfrac{1}{E_{\mathrm{b}}}$ in the RG flow equations. Preventing further
numerical evolution towards the IR.

When considering a model involving also Goldstone modes, which enter the RG
flow as advection, cf. Refs. [142, 154, 110, 144, 145], the huge gradients
introduced via the chemical potential induce a shock wave that propagates
towards smaller $|\sigma|$.

A possible reason, why this subtle dynamics in field space was not discussed
in detail in the past within FRG studies might be that it is hardly visible
and understandable on the level of the scale dependent effective potential
$U(t,\sigma)$ itself, because it is the integral of $u(t,\sigma)$, where jumps
might only show up as tiny cusps, compare Sub.Sec. VI.2. Additionally, a lot
of previous studies did not perform calculations at sufficiently low
temperatures, where this effect is not completely washed out by the thermal
distribution functions (22). At this point the underlying reason for and
resolution of this practical problem are not clear to the authors. The
disconnection between the fermionic sector and the dynamics of $u(t,\sigma)$
in LPA, the regulator choice or the formulation181818A formal classification
and proper weak formulation of the PDEs (19) and (32) might be paramount to
understand the nature of the arising discontinuities and weak/physical
solutions in their presence. of the flow equation could be possible reasons
for this intricate problem. Solving this problem is beyond the scope of the
current work. Consequently, we claim that understanding and/or capturing this
effect correctly will be one of the central challenges in FRG at non-zero
$\mu$ within the next years – independent of the specific models.

For related discussions on these novel findings, we also refer to Refs. [154,
145, 189].

###  IV.5 Numerical implementation

In this subsection we present the explicit numeric implementation of the RG
flow equation (32). Thereby we make use of our findings of the previous
sections and related publications [110, 143, 144, 145] on the relation between
RG flow equations and flow equations in the true sense of the word.

To adequately capture highly non-linear diffusive effects as well as the
position-dependent source/sink terms, we use the so called Kurganov-Tadmor
(KT) scheme [182]. The KT scheme is an established numerical scheme, that can
cope with the aforementioned challenges. Furthermore, the KT scheme will allow
for a straightforward generalization of our numerics, if one extends the GNY
model to a model with continuous chiral symmetry group and related Goldstone
modes, which enter the LPA RG flow equation as non-linear advective
contributions [142, 154, 110, 144, 145]. The original KT scheme, which was
presented by A. Kurganov and E. Tadmor in Ref. [182], was in fact developed
primarily for non-linear hyperbolic PDEs. However, Ref. [182] also presented
generalizations of their scheme for PDEs that involve source/sink and
diffusion terms. For the scope of this work, we only use these extensions.

The scheme is formulated in a semi-discrete manner. Hence, position (field)
space is discretized via a finite volume discretization, while the temporal
coordinate is formally continuous. This enables the use of up-to-date generic
ODE time-steppers to solve the discretized equations.191919Within this work,
we use the time-stepper solve_ivp with the LSODA option using an Adams/BDF
method with automatic stiffness detection and switching from the SciPy 1.0
library [190]. For further details, see App. F. We also crosschecked some of
our results with the code-basis of Refs. [110, 143, 144, 145], which is
implemented with Mathematica [191] using NDSolve.

The explicit discretization of Eqs. (32), i.e. the diffusion flux (33) and
sink term (34) is as follows. The compact computational domain
$[0,\sigma_{\mathrm{max}}]$, which was introduced in Sub.Sec. IV.3, is
subdivided into $n\in\mathbb{N}$ equally sized (finite) volume cells of size
$\Delta\sigma$ centered at positions $\sigma_{i}$, $i=0,1,\ldots,n-1$. The
zeroth volume cell is centered at ${\sigma=\sigma_{0}=0}$, while the last cell
is centered at $\sigma_{n-1}=\sigma_{\mathrm{max}}$. Within a single volume
cell $\sigma_{i}$, the cell-average of the “fluid” $u(t,\sigma)$ is denoted as
$\bar{u}_{i}(t)$. The actual computation and scheme is entirely formulated via
these cell-averages that evolve with time.

For the implementation of boundary conditions two (one) so-called ghost cells
are required at each interval boundary when considering a problem involving
advective (solely diffusive) contributions in the semi-discrete KT-scheme. For
the GNY model which has no advective contribution in its LPA flow equation one
ghost cell at each interval boundary is sufficient. The ghost cells are of
size $\Delta\sigma$ and centered at $\sigma_{-1}=-\Delta\sigma$ and
$\sigma_{n}=\sigma_{\mathrm{max}}+\Delta\sigma$. As described in Sub.Sec. IV.3
we use the anti-symmetry of $u(t,\sigma)$ as boundary condition to fix the
cell-averages of the zeroth cell and first ghost cell to $\bar{u}_{0}(t)=0$
and $\bar{u}_{-1}(t)=-\bar{u}_{1}(t)$ for all times $t$ respectively. For the
ghost cell at $\sigma_{n}$ we use linear extrapolation, thus for the cell
average $\bar{u}_{n}(t)=2\,\bar{u}_{n-1}(t)-\bar{u}_{n-2}(t)$ for all $t$.

Additionally, we introduce the grid of cell interfaces, which are positioned
at $\sigma_{i+\frac{1}{2}}\equiv\sigma_{i}+\tfrac{\Delta\sigma}{2}$.

Within this setup, the semi-discrete scheme for the non-linear diffusion-sink
equation (32) and the evolution of the cell-averages $\bar{u}_{i}(t)$ reads
[182]

$\displaystyle\partial_{t}\bar{u}_{i}=\tfrac{1}{\Delta\sigma}\,\big{(}P_{i+\frac{1}{2}}-P_{i-\frac{1}{2}}\big{)}+S_{i}\,.$
(38)

Here, $S_{i}=S(t,\sigma_{i})$, where $S(t,\sigma)$ is the sink term (34)
(after analytic evaluation of the $\sigma$-derivative). For a detailed comment
on the implementation of the sink term, we refer to App. G.

Furthermore, according to Ref. [182] we specify

$\displaystyle P_{i+\frac{1}{2}}=\,$
$\displaystyle\frac{Q\big{(}t,\bar{u}_{i},\frac{\bar{u}_{i+1}-\bar{u}_{i}}{\Delta\sigma}\big{)}+Q\big{(}t,\bar{u}_{i+1},\frac{\bar{u}_{i+1}-\bar{u}_{i}}{\Delta\sigma}\big{)}}{2}\,,$
(39)

where for our purpose $Q(t,u,\partial_{\sigma}u)=Q(t,\partial_{\sigma}u)$ is
the diffusion flux (33), In fact, the independence of $u$ of our special
diffusion flux leads to further simplification of the original
$P_{i+\frac{1}{2}}$ from Ref. [182]. For an alternative implementation of the
diffusion flux in the KT scheme, see, e.g., Ref. [192].

For more details on the KT-scheme and its original design for advection driven
problems, we refer to the original publication Ref. [182]. For a comprehensive
discussion of the KT-scheme in the context of RG flow equations, we refer to
Refs. [110, 143, 144, 145], where also references for further reading and
extensions of the scheme are listed.

##  V Infinite-$N$ analysis within the FRG –
Consistency checks and the UV initial condition

In this section, we recap some of the well-known (analytical) results for the
infinite-$N$ limit (mean-field) of the GNY model. Furthermore, we demonstrate
that the FRG is capable of reproducing the latter results analytically and
numerically. Additionally, these mean-field calculations are used to motivate
a proper UV initial condition for the flow equation (19), but also serve as a
consistency check of our numerical implementation in the limit
$N\rightarrow\infty$, thus in the limit of vanishing bosonic fluctuations.

###  V.1 Mean-field, infinite-$N$, and FRG

Within the FRG framework (arguably even in general), the term “mean-field
approximation” has no universal, agreed upon formal definition. Usually
performing calculations on “mean-field level” in the context of fermionic
models refers to computations disregarding bosonic quantum and thermal
fluctuations and including only fermionic ones. Whether this includes
fermionic vacuum fluctuations and/or fermionic contributions beyond the
effective potential usually depends on the work under consideration. In the
FRG one way to define a mean-field approximation is the usage of a LPA
truncation disregarding the bosonic fluctuations. On a practical level, this
results in simply ignoring the bosonic loop contribution(s). On LPA level this
can be formally achieved by taking the infinite-$N$ limit for the LPA flow
equation under consideration, in this work Eq. (19), after appropriate
rescalings, like the ones introduced in Eqs. (15) and (16).

However, it can be shown that in general the infinite-$N$ limit and “ignoring
the bosonic loop in a LPA truncation” is not the same procedure. In fact, if
starting with a more advanced FRG truncation scheme, like LPA′, which also
includes wave function renormalizations, one finds, that even in the
infinite-$N$ limit, there are fermionic loop contributions to the bosonic
wave-function renormalization, e.g. Ref. [12, 99]. Hence, in general the order
of “limits” plays a crucial role and “choosing a truncation scheme in FRG” and
“taking the infinite-$N$ limit” do not commute. Furthermore in models
including Goldstone bosons a large-$N$ limit for fermions, e.g., in the
large-$N_{f}$ limit for chiral fermion flavors, does not lead to the desired
suppression of bosonic modes, because the Goldstone modes do not form as
flavor singlets. Another degree of freedom, e.g., the number of colors
$N_{c}$, has to be used to facilitate the large-$N$ limit and the desired
suppression of bosonic fluctuations.

For our purpose and the sake of simplicity, we simply start by definition on
the level of the LPA and take all limits like the infinite-$N$ limit or the
zero-$T$ and zero-$\mu$ limit afterward. Hence, within our truncation, the
mean-field limit and the infinite-$N$ limit are considered to be identical,
which simplifies the discussion and allows to make direct contact with
established conventional mean-field computations for the GN model, cf. Refs.
[52, 95], which consider only fermionic fluctuations on the level of the
effective potential.

Performing the large-$N$ limit For the flow equation (19), yields

$\displaystyle\lim\limits_{N\rightarrow\infty}\partial_{t}U(t,\sigma)=\vphantom{\Bigg{(}\Bigg{)}}$
(40) $\displaystyle=\,$
$\displaystyle\frac{d_{\gamma}}{\pi}\,\frac{k^{3}(t)}{2E_{\textrm{f}}(t,\sigma)}\,\big{(}1-n_{\mathrm{f}}[\beta(E_{\textrm{f}}(t,\sigma)+\mu)]-\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle-
n_{\mathrm{f}}[\beta(E_{\textrm{f}}(t,\sigma)-\mu)]\big{)}\,.\vphantom{\Bigg{(}\Bigg{)}}$

We find that the former PDE decouples in field space and reduces to a first
order ordinary differential equation in $t$ for each $\sigma$ separately. In
the fluid dynamic picture, on the level of
$u(t,\sigma)=\partial_{\sigma}U(t,\sigma)$, this is rather intuitive, since
the fermionic contribution to the flow equation (32) presents as a local time-
dependent source/sink term (34) and the spatial movement of the fluid (via
diffusion in field space) is totally suppressed.

When disregarding bosonic fluctuations completely in mean-field approximation
all three model variants – GN, bGN and GNY – discussed Sec. II and Sub.Sec.
III.4 are equivalent. The bGN and GNY models are identical in mean-field
approximation and the bGN model as the bosonized version of the GN Model is in
general physically equivalent to the GN model as discussed in Sec. II.

###  V.2 The mean-field potential and asymptotic freedom

Due to the decoupling in field space, the differential equation (40) can be
integrated analytically in $k(t)$. Using the definition of the RG time (13),
which implies $\partial_{t}=-k\,\partial_{k}$, we find

$\displaystyle U_{k=0}(\sigma)=\,$ $\displaystyle
U_{k=\Lambda}(\sigma)+\frac{d_{\gamma}}{\pi}\int_{0}^{\Lambda}\mathrm{d}k\,\frac{k^{2}}{2E_{\mathrm{f}}}\,\big{(}1-n_{\mathrm{f}}[\beta(E_{\mathrm{f}}+\mu)]-n_{\mathrm{f}}[\beta(E_{\mathrm{f}}-\mu)]\big{)}=\vphantom{\Bigg{(}\Bigg{)}}$
(41) $\displaystyle U_{0}(\sigma)=\,$ $\displaystyle
U_{\Lambda}(\sigma)+\frac{d_{\gamma}}{\pi}\,\bigg{[}\frac{k}{2}\,\big{(}E_{\mathrm{f}}+\tfrac{1}{\beta}\ln\big{[}1+\mathrm{e}^{-\beta(E_{\mathrm{f}}+\mu)}\big{]}+\tfrac{1}{\beta}\ln\big{[}1+\mathrm{e}^{-\beta(E_{\mathrm{f}}-\mu)}\big{]}\big{)}\bigg{]}_{0}^{\Lambda}\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle-\frac{d_{\gamma}}{2\pi}\int_{0}^{\Lambda}\mathrm{d}k\,\big{(}E_{\mathrm{f}}+\tfrac{1}{\beta}\,\ln\big{[}1+\mathrm{e}^{-\beta(E_{\mathrm{f}}+\mu)}\big{]}+\tfrac{1}{\beta}\,\ln\big{[}1+\mathrm{e}^{-\beta(E_{\mathrm{f}}-\mu)}\big{]}\big{)}\,,$

where the UV-initial condition for this trivial integrable “RG-flow” is given
by the classical UV-potential, see Sub.Sec. III.4,

$\displaystyle U_{\Lambda}\equiv
U(t=0,\sigma)=\tfrac{1}{2g^{2}}\,(h\sigma)^{2}$ (42)

and where we introduced the notation

$\displaystyle U_{k}(\sigma)\equiv U_{k(t)}(\sigma)\equiv U(t,\sigma)$ (43)

for the potential at a given RG scale $k$.

In the second line of Eq. (41) we integrated by parts in order to recover the
usual expression for the grand canonical potential density for
$\Lambda\rightarrow\infty$ (up to an infinite constant $\sim
k\,E_{\mathrm{f}}\big{|}_{k=\Lambda}$), cf. Refs. [48, 49, 52, 95]. The first
terms of the integrands in the first and third line of Eq. (41) lead to
divergences. These divergences have to be canceled by “renormalizing” the
coupling $g^{2}$ such that the IR observables are finite and could in
principle be matched with experimental observations. The renormalized version
of this so called vacuum contribution (the only contribution, that does not
depend on $\mu$ and $T$) is directly linked to the initial condition of our
RG-flows, when solving the PDE (19) and ODE (40) (numerically), see below.

To this end, we turn to the $\mu=0$ and $\beta\rightarrow\infty\Leftrightarrow
T\rightarrow 0$ limit of Eq. (41),

$\displaystyle U_{0;\mathrm{vac}}(\sigma)\equiv\,$
$\displaystyle\lim\limits_{\mu,T\rightarrow
0}U_{0}(\sigma)=\vphantom{\Bigg{(}\Bigg{)}}$ (44) $\displaystyle=\,$
$\displaystyle\frac{h^{2}\sigma^{2}}{2g^{2}}+\frac{d_{\gamma}}{2\pi}\int_{0}^{\Lambda}\mathrm{d}k\,\frac{k^{2}}{\sqrt{k^{2}+h^{2}\sigma^{2}}}\,.\vphantom{\Bigg{(}\Bigg{)}}$

and study the corresponding gap equation (45) at possible non-trivial minima
$\sigma_{0}\neq 0$, cf. Refs. [3, 52, 95],

$\displaystyle 0\overset{!}{=}\,$
$\displaystyle\frac{1}{h^{2}\sigma_{0}}\,\partial_{\sigma}U_{0;\mathrm{vac}}(\sigma)\big{|}_{\sigma=\sigma_{0}}=\vphantom{\Bigg{(}\Bigg{)}}$
(45) $\displaystyle=\,$
$\displaystyle\frac{1}{g^{2}}-\frac{d_{\gamma}}{2\pi}\int_{0}^{\Lambda}\mathrm{d}k\,\frac{k^{2}}{(k^{2}+h^{2}\sigma_{0}^{2})^{\frac{3}{2}}}=\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle=\,$
$\displaystyle\frac{1}{g^{2}}+\frac{d_{\gamma}}{2\pi}\,\bigg{[}\frac{k}{\sqrt{k^{2}+h^{2}\sigma_{0}^{2}}}\bigg{|}_{0}^{\Lambda}-\int_{0}^{\Lambda}\mathrm{d}k\,\frac{1}{\sqrt{k^{2}+h^{2}\sigma_{0}^{2}}}\bigg{]}=\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle=\,$
$\displaystyle\frac{1}{g^{2}}+\frac{d_{\gamma}}{2\pi}\,\bigg{[}\Big{[}1+\big{(}\tfrac{h\sigma_{0}}{\Lambda}\big{)}^{2}\Big{]}^{-\frac{1}{2}}-\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle-\mathrm{artanh}\bigg{(}\Big{[}1+\big{(}\tfrac{h\sigma_{0}}{\Lambda}\big{)}^{2}\Big{]}^{-\frac{1}{2}}\bigg{)}\bigg{]}\,.\vphantom{\Bigg{(}\Bigg{)}}$

Hence, as a first result, assuming that $\sigma_{0}$ is non-zero and finite,
we can study the asymptotic behavior of the equation for $\Lambda\ggg h$,

$\displaystyle\tfrac{1}{g^{2}}=\,$
$\displaystyle\tfrac{d_{\gamma}}{2\pi}\,\big{[}-1+\tfrac{1}{2}\ln\big{(}\big{(}\tfrac{2\Lambda}{h\sigma_{0}}\big{)}^{2}\big{)}\big{]}+\mathcal{O}\big{(}\big{(}\tfrac{h\sigma_{0}}{\Lambda}\big{)}^{2}\big{)}\,.$
(46)

This reflects the asymptotically free behavior of the four-Fermi coupling [1,
10, 9] of the original Gross-Neveu model (1), since

$\displaystyle\lim\limits_{\frac{\Lambda}{h}\rightarrow\infty}g^{2}=0\,.$ (47)

Furthermore, we can use (46) and solve for $\sigma_{0}$,

$\displaystyle\sigma_{0}=\pm\tfrac{2\Lambda}{h}\,\mathrm{e}^{-\big{(}\frac{4\pi}{d_{\gamma}g^{2}}+\frac{1}{2}\big{)}}\,.$
(48)

Hence, due to the asymptotic free behavior of $g^{2}$, there is
$\mathbb{Z}_{2}$ (discrete chiral) symmetry breaking via two non-trivial
minima for all non-zero $g^{2}$ in vacuum and for infinite-$N$ – a central
feature of the GN model [1].

Furthermore, we can insert the results from the gap equation (45) in Eq. (42)
and solve for the UV potential

$\displaystyle U_{\Lambda}(\sigma)=\,$
$\displaystyle\tfrac{d_{\gamma}}{4\pi}\,(h\sigma)^{2}\,\bigg{[}\mathrm{artanh}\bigg{(}\Big{[}1+\big{(}\tfrac{h\sigma_{0}}{\Lambda}\big{)}^{2}\Big{]}^{-\frac{1}{2}}\bigg{)}-\vphantom{\Bigg{(}\Bigg{)}}$
(49)
$\displaystyle-\Big{[}1+\big{(}\tfrac{h\sigma_{0}}{\Lambda}\big{)}^{2}\Big{]}^{-\frac{1}{2}}\bigg{]}\,.\vphantom{\Bigg{(}\Bigg{)}}$

In the limit $\tfrac{\Lambda}{h}\to\infty$ we approach the Gaussian fixed-
point (UV fixed-point) for the bosonized GN model (3), which becomes clear
when considering dimensionless quantities

$\displaystyle\tilde{h}=\tfrac{1}{\Lambda}\,h\,,$
$\displaystyle\tilde{U}_{\Lambda}(\sigma)=\tfrac{1}{\Lambda^{2}}\,U_{\Lambda}(\sigma)\,.$
(50)

Both, $\tilde{h}$ and $\tilde{U}_{\Lambda}(\sigma)$, vanish in the limit
$\tfrac{\Lambda}{h}\to\infty$. This implies – as expected – that the bGN model
(3) in vacuum in the infinite-$N$ limit, is also asymptotically free.

Turning to the IR mean-field potential (44) and using the previous results, we
find

$\displaystyle U_{0;\mathrm{vac}}(\sigma)=\vphantom{\Bigg{(}\Bigg{)}}$ (51)
$\displaystyle=\,$
$\displaystyle\tfrac{d_{\gamma}}{4\pi}\,(h\sigma)^{2}\,\bigg{[}\mathrm{artanh}\bigg{(}\Big{[}1+\big{(}\tfrac{h\sigma_{0}}{\Lambda}\big{)}^{2}\Big{]}^{-\frac{1}{2}}\bigg{)}-\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle-\Big{[}1+\big{(}\tfrac{h\sigma_{0}}{\Lambda}\big{)}^{2}\Big{]}^{-\frac{1}{2}}-\mathrm{artanh}\bigg{(}\Big{[}1+\big{(}\tfrac{h\sigma}{\Lambda}\big{)}^{2}\Big{]}^{-\frac{1}{2}}\bigg{)}\bigg{]}+\vphantom{\Bigg{(}\Bigg{)}}$
$\displaystyle+\tfrac{d_{\gamma}}{4\pi}\,\Lambda^{2}\,\sqrt{1+\big{(}\tfrac{h\sigma}{\Lambda}\big{)}^{2}}\,.\vphantom{\Bigg{(}\Bigg{)}}$

Considering the first derivative of Eq. (51) one can verify
${0=\partial_{\sigma}U_{0;\mathrm{vac}}(\sigma)|_{\sigma_{0}}}$ which has to
hold by construction and we note that

$\displaystyle\partial_{\sigma}^{2}U_{0;\mathrm{vac}}(\sigma)|_{\sigma_{0}}=\frac{d_{\gamma}}{2\pi}\,\frac{\Lambda^{3}}{[\Lambda^{2}+(h\sigma_{0})^{2}]^{3/2}}\,h^{2},$
(52)

which is manifest positive – again consistent with the notion of $\sigma_{0}$
as a non-trivial minimum by construction – and corresponds to the squared
curvature mass $m_{\sigma}^{2}$ of the sigma mode in vacuum.

Considering the limit $\tfrac{\Lambda}{h}\rightarrow\infty$, the divergent
contributions of the two $\mathrm{artanh}$ cancel exactly. For the last term,
we use

$\displaystyle\Lambda^{2}\,\sqrt{1+\big{(}\tfrac{h\sigma}{\Lambda}\big{)}^{2}}=\Lambda^{2}+\tfrac{1}{2}\,(h\sigma)^{2}+\mathcal{O}\big{(}\big{(}\tfrac{h\sigma}{\Lambda}\big{)}^{2}\big{)}\,,$
(53)

which results in an unobservable infinite constant and a finite contribution.
In total we find the well-known renormalized vacuum IR effective potential,
cf. Refs. [52, 95],

$\displaystyle U_{0;\mathrm{vac}}(\sigma)=\vphantom{\bigg{(}\bigg{)}}$ (54)
$\displaystyle=\,$
$\displaystyle\tfrac{d_{\gamma}}{8\pi}\,(h\sigma)^{2}\,\big{(}\big{[}\ln\big{(}\tfrac{h\sigma}{h\sigma_{0}}\big{)}^{2}\big{]}-1\big{)}-\tfrac{d_{\gamma}}{4\pi}\,\Lambda^{2}+\mathcal{O}\big{(}\big{(}\tfrac{h\sigma}{\Lambda}\big{)}^{2}\big{)}\,,\vphantom{\bigg{(}\bigg{)}}$

with its global minimum at $\pm\sigma_{0}$ and with a corresponding squared
curvature mass of

$\displaystyle m_{\sigma}^{2}=\tfrac{d_{\gamma}}{2\pi}\,h^{2}\,,$ (55)

cf. Ref. [1]. Finally, the result (54) can be used to also simplify the full
IR potential in medium (41) by replacing the vacuum contributions,

$\displaystyle U_{0}(\sigma)=\vphantom{\Bigg{(}\Bigg{)}}$ (56)
$\displaystyle=\,$
$\displaystyle\tfrac{d_{\gamma}}{8\pi}\,(h\sigma)^{2}\,\big{(}\big{[}\ln\big{(}\tfrac{h\sigma}{h\sigma_{0}}\big{)}^{2}\big{]}-1\big{)}-\tfrac{d_{\gamma}}{4\pi}\,\Lambda^{2}+\mathcal{O}\big{(}\big{(}\tfrac{h\sigma}{\Lambda}\big{)}^{2}\big{)}-\tfrac{d_{\gamma}}{2\pi}\int_{0}^{\Lambda}\mathrm{d}p\,\big{(}\tfrac{1}{\beta}\,\ln\big{[}1+\mathrm{e}^{-\beta(E_{\mathrm{f}}+\mu)}\big{]}+\tfrac{1}{\beta}\,\ln\big{[}1+\mathrm{e}^{-\beta(E_{\mathrm{f}}-\mu)}\big{]}\big{)}\,.\vphantom{\Bigg{(}\Bigg{)}}$

Let us again point out that mean-field IR potentials (54) and (56) of the
bGN/GNY model (3)/(6) are unique. For $\Lambda\rightarrow\infty$ we can read
off from Eq. (54) that the model contains only a single internal dimensionful
parameter, for instance $h$, because the bosonic field $\sigma$ is
dimensionless and $h$ and $\sigma$ only appear in combination. All other
dimensionful quantities ($\mu$, $T$, $\Lambda$, $U$) can be expressed via this
single parameter and results for different specifications of the reference
parameter can be mapped into each other via simple rescaling. A direct
consequence is that also the mean-field phase diagram is unique, as is
discussed in Sub.Sec. V.3.

###  V.3 Analytical results for the mean-field phase diagram in the
renormalized limit $\Lambda\rightarrow\infty$

In this subsection we will discuss the mean-field potential $U_{0}(\sigma)$ of
Eq. (56) in the limit $\Lambda\rightarrow\infty$. We will compare our
analytical results of this subsection with the existing renormalized mean
field results of specifically Refs. [52, 95] for the Gross-Neveu model at
vanishing bare fermion mass. Of special interest is the phase-diagram (see
Figure 1 of Ref. [95] or Fig. 3 of the preceding publication [52]) and the
underlying renormalized grand canonical potential density, see, e.g. Eq. (4)
of Ref. [95]. The results of this section serve as a reference in the
consistence check of our numerical implementation in Sub.Sec. V.5.

Consider the mean-field potential $U_{0}(\sigma)$ of Eq. (56) in the
renormalized limit $\Lambda\rightarrow\infty$,

$\displaystyle\lim_{\Lambda\rightarrow\infty}U_{0}(\sigma)-\tfrac{d_{\gamma}}{4\pi}\,\Lambda^{2}=\vphantom{\Bigg{(}\Bigg{)}}$
(57) $\displaystyle=\,$
$\displaystyle\tfrac{d_{\gamma}}{8\pi}\,(h\sigma)^{2}\,\big{(}\ln\big{[}\big{(}\tfrac{h\sigma}{h\sigma_{0}}\big{)}^{2}\big{]}-1\big{)}-\tfrac{d_{\gamma}}{4\pi}\int_{-\infty}^{\infty}\mathrm{d}k\,\big{(}\tfrac{1}{\beta}\,\ln\big{[}1+\mathrm{e}^{-\beta(E_{\mathrm{f}}+\mu)}\big{]}+\tfrac{1}{\beta}\,\ln\big{[}1+\mathrm{e}^{-\beta(E_{\mathrm{f}}-\mu)}\big{]}\big{)}\,,\vphantom{\Bigg{(}\Bigg{)}}$

where we used the symmetry ($k\rightarrow-k$) of the remaining integral. We
identify several noteworthy terms in Eq. (57): the remaining integral is the
canonical one-loop, convergent medium contribution of fermions in $1+1$
dimensions, the vacuum contribution carries a, for $1+1$ dimensions typical
[193, 194, 102], term proportional to $(h\sigma)^{2}\,\log(h\sigma)$ and we
subtracted a (in $\sigma$, $\mu$ and $\beta$) constant but otherwise divergent
contribution proportional $\Lambda^{2}$ to perform the limit
$\Lambda\rightarrow\infty$. In observables this divergent but constant
contribution cancels (since $U_{0}(\sigma)$ is only defined up to a constant)
and therefore it can be subtracted from $U_{0}(\sigma)$ without further
consequences. For the following we therefore consider

$\displaystyle V(\sigma)\equiv\,$
$\displaystyle\lim_{\Lambda\rightarrow\infty}U_{0}(\sigma)-\tfrac{d_{\gamma}}{4\pi}\,\Lambda^{2}=\vphantom{\bigg{(}\bigg{)}}$
(58) $\displaystyle=\,$
$\displaystyle\tfrac{1}{4\pi}\,\sigma^{2}\,[\log(\sigma^{2})-1]+I_{V}(T,\mu)\,,\vphantom{\bigg{(}\bigg{)}}$

where we used $d_{\gamma}=2$ as well as w.l.o.g. $h=1$ and $\sigma_{0}=1$ and
abbreviate the medium contribution with

$\displaystyle I_{V}(T,\mu)\equiv\,$
$\displaystyle-\tfrac{T}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}k\,\ln\big{[}1+\mathrm{e}^{-\tfrac{1}{T}(E_{\mathrm{f}}+\mu)}\big{]}-\vphantom{\Bigg{(}\Bigg{)}}$
(59)
$\displaystyle-\tfrac{T}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}k\,\ln\big{[}1+\mathrm{e}^{-\tfrac{1}{T}(E_{\mathrm{f}}-\mu)}\big{]}\,.\vphantom{\Bigg{(}\Bigg{)}}$

For a given $\mu$ and $T$ the physical value of the condensate $\sigma$ is
found by minimization of $V(\sigma)$, cf. Sub.Sec. III.5. When evaluated at
its global minimum $V(\sigma)$ is the renormalized grand canonical potential
per spatial volume. The renormalized (FRG) mean-field result of Eq. (58) for
$V(\sigma)$ obtained in the infinite-$N$ limit of the LPA flow equation with a
one-dimensional LPA-optimized regulator (40) coincides with the renormalized
mean-field results in literature, see, e.g., Refs. [52, 95], and explicitly
Eq. (4) of Ref. [95]. In the rest of this subsection we derive and
recapitulate known MF results for the homogeneous phase diagram of the GN
model, see, e.g., [52, 95]. This discussion may be skipped by readers familiar
with these results.

Figure 2:  Phase diagram of the renormalized GN model in the infinite-$N$
limit. The blue ($\alpha_{4}=0$) and red ($\alpha_{2}=0$) lines are obtained
from the Ginzburg-Landau type expansion and Eq. (64) and Eq. (65) for $n=2$
respectively. The critical point denoted by the red star separates the second
order phase boundary (red-dashed line) from the first order phase boundary
(solid green line). The red-dotted line and the green-dotted line are the left
and right spinodal lines. The first order phase boundary and the right
spinodal have been obtained by explicit numerical integration [195] and
subsequent repeated, local, numerical minimization [196] of the renormalized
grand canonical potential per spatial volume (58). The critical temperature
$T_{\mathrm{C}}$ is marked on the $T$-axis with a red disk, while the end
points at $\mu_{\mathrm{L}}$ and $\mu_{\mathrm{R}}$ of the spinodal lines are
marked on the $\mu$-axis with red and green disks respectively. The first-
order phase transition at zero temperature is marked with a green square at
$\mu_{1}$.

The phase diagram of Fig. 2 identical to the one presented in Fig. 1 of Ref.
[95]. By construction – due to the renormalization condition of Eq. (45) –
discrete chiral symmetry is broken with $\sigma=\sigma_{0}=1$ in the vacuum at
$T=\mu=0$. At high chemical potentials and temperatures the discrete chiral
symmetries is restored and therefore $\sigma=0$. At intermediate temperatures
and chemical potentials one observes a first order phase transition line
starting at zero temperature and non-zero chemical potential $\mu_{1}$ and
ending in a critical point at $(\mu_{\mathrm{CP}},T_{\mathrm{CP}})$. Above
$T_{\mathrm{CP}}$ discrete chiral symmetry is restored across a second order
phase transition line starting at the critical point and ending on the
$(\mu=0)$-axis at a non-zero temperature $T_{\mathrm{C}}$.

The position of the critical point as well as $\mu_{1}$ and $T_{\mathrm{C}}$
can be computed with the help of known functions without the need of numerical
minimization of the potential $V(\sigma)$. We devote the remainder of this
subsection to the derivation of the values for $T_{\mathrm{C}}$,
$T_{\mathrm{CP}}$, $\mu_{\mathrm{CP}}$ and $\mu_{1}$ given in Refs. [95, 52,
53] as well as expressions for specific lines in Fig. 2 .

For $T\rightarrow 0$ the integral $I_{V}(0,\mu)$ can be performed
analytically,

$\displaystyle I_{V}(0,\mu)=\vphantom{\bigg{(}\bigg{)}}$ (60)
$\displaystyle=\,$
$\displaystyle\tfrac{1}{2\pi}\,\Big{(}\mathrm{arsinh}\Big{[}\sqrt{\big{(}\tfrac{\mu}{\sigma}\big{)}^{2}-1}\Big{]}-\mu\,\sqrt{\mu^{2}-\sigma^{2}}\Big{)}\,\Theta\big{[}\big{(}\tfrac{\mu}{\sigma}\big{)}^{2}-1\big{]}\,,\vphantom{\bigg{(}\bigg{)}}$

and we note $I_{V}(0,0)=0$ as all vacuum contributions are already integrated
out and included in Eq. (58). An analysis of Eq. (60) reveals, that the
extremum at $\sigma=0$ becomes a local minimum for
$\mu>\mu_{\mathrm{L}}=\tfrac{1}{2}$. For
${\mu\in[\mu_{\mathrm{L}},\mu_{\mathrm{R}}]}$ the potential has three minima
(one at $\sigma=0$ and two at $\sigma=\pm 1$) and at
$\mu_{1}=\tfrac{1}{\sqrt{2}}\simeq 0.707107$ all local minima become global
minima signaling a first order phase transition at $\mu_{1}$. The notable
chemical potentials $\mu_{\mathrm{L}}$, $\mu_{1}$ and $\mu_{\mathrm{R}}$ are
marked on the $(T=0)$-axis in Fig. 2.

For $T>0$ the integral $I_{V}(T,\mu)$ can be performed numerically or it can
be rewritten in terms of known functions, see, e.g., Ref. [102] Eqs. (3.11)
ff. and Refs. [193, 194],

$\displaystyle I_{V}(T,\mu)=\vphantom{\Bigg{(}\Bigg{)}}$ (61)
$\displaystyle=\,$
$\displaystyle-\tfrac{2}{\pi}\,T^{2}\sum_{m=1}^{\infty}\tfrac{(-1)^{m+1}}{m}\,K_{1}\big{(}m\,\tfrac{\sigma}{T}\big{)}\,\cosh\big{(}m\,\tfrac{\mu}{T}\big{)}\,,\vphantom{\Bigg{(}\Bigg{)}}$

with the modified Bessel function of second kind $K_{n}(x)$ at $n=1$. This
representation is particularly useful when one is interested in a Ginzburg-
Landau-type series expansion of $V(\sigma)$, cf. Refs. [19, 59, 197], around
$\sigma=0$ at $T>0$, e.g., to study second order phase transitions within a
mean-field approach. Using Eqs. (58), (60) and (61) one can show202020The
computation at $T>0$ is based on the ascending series for $K_{1}(x)$, see,
e.g., Ref. [198] Eq. (9.6.11), and various identities for the polylogarithm
and Riemann zeta function, see e.g., Refs. [198, 199]. For $T=0$ the
coefficients $\alpha_{2n}$ can be computed directly from Eq. (60) with the
help of the series expansion for the $\mathrm{arcsinh}(x)$ for $|x|>1$, see,
e.g., Ref. [198] Eq. (4.6.31).,

$\displaystyle V(\sigma)=\,$
$\displaystyle\sum_{n=0}^{\infty}\alpha_{2n}\,\sigma^{2n}\,,$ (62)

where the zeroth order coefficient reads

$\displaystyle\alpha_{0}=\,$
$\displaystyle-\tfrac{\pi}{6}\,T^{2}-\tfrac{1}{2\pi}\,\mu^{2}\,,$ (63)

and the quadratic order coefficient is given by

$\displaystyle\alpha_{2}=\begin{cases}\tfrac{1}{2\pi}\,\ln(2\mu)\,,&\text{for}\quad
T=0\,,\mu>0\,,\vphantom{\bigg{(}\bigg{)}}\\\ -\tfrac{1}{2\pi}\,[\gamma-\ln(\pi
T)]\,,&\text{for}\quad T>0\,,\mu=0\,,\vphantom{\bigg{(}\bigg{)}}\\\
-\tfrac{1}{2\pi}\,\big{[}\gamma-\ln(2T)+\mathrm{DLi}_{0}\big{(}\tfrac{\mu}{T}\big{)}\big{]}\,,&\text{for}\quad
T>0\,,\mu>0\,.\vphantom{\bigg{(}\bigg{)}}\end{cases}$ (64)

The higher order coefficients $n\geq 2$ can also be expressed in terms of
known functions,

$\displaystyle\alpha_{2n}=\,$
$\displaystyle\begin{cases}-\tfrac{1}{\pi}\,\tfrac{(2n-3)!!}{2^{n+1}(n-1)!n!}\,\mu^{2-2n}\,,&\text{for}\quad
T=0\,,\mu>0\,,\vphantom{\bigg{(}\bigg{)}}\\\
\tfrac{(-1)^{n}}{\pi^{2n-1}}\,\tfrac{4^{1-2n}(4^{n}-2)(2n-2)!}{(n-1)!n!}\,\zeta(2n-1)\,T^{2-2n}\,,&\text{for}\quad
T>0\,,\mu=0\,,\vphantom{\bigg{(}\bigg{)}}\\\
-\tfrac{1}{\pi}\,\tfrac{2^{1-2n}}{n!(n-1)!}\,\mathrm{DLi}_{2-2n}\big{(}\tfrac{\mu}{T}\big{)}\,T^{2-2n}\,,&\text{for}\quad
T>0\,,\mu>0\,.\vphantom{\bigg{(}\bigg{)}}\end{cases}$ (65)

Here, we defined

$\displaystyle\mathrm{DLi}_{n}(z)\equiv\big{[}\tfrac{\partial}{\partial
s}\mathrm{Li}_{s}(-\mathrm{e}^{z})+\tfrac{\partial}{\partial
s}\mathrm{Li}_{s}(-\mathrm{e}^{-z})\big{]}_{s=n}\,,$ (66)

with the polylogarithm $\mathrm{Li}_{m}(x)$212121The numerical evaluation of
$\mathrm{DLi}_{n}(z)$ involves derivatives of the the polylogarithm
$\mathrm{Li}_{m}(x)$ w.r.t. to $m$ evaluated at $m=n\leq 0$. For the scope of
this work we compute $\mathrm{DLi}_{n}(z)$ numerically to 32 significant
digits using Maple [200].. We also used the Riemann zeta function $\zeta(x)$
and the Euler-Mascheroni constant $\gamma\simeq 0.577216$.

The expression (63) for $-\alpha_{0}$ is the Stefan-Boltzmann pressure of a
massless free Fermi gas in $1+1$ dimensions, cf. Ref. [46]. The logarithmic
divergent contribution $\log(\sigma^{2})$ for $\sigma\rightarrow 0$ to
$\alpha_{2}$ stemming from the vacuum term in Eq. (58) gets canceled exactly
in Eq. (64) by a corresponding contribution from the medium term.

The expansion (62) around $\sigma=0$ can be used to compute the second-order
phase-boundary in terms of known functions including the critical end point.
The second order phase transition between the restored and a broken phase with
small $\sigma>0$ occurs at $\alpha_{2}=0$ while $\alpha_{4}>0$. In the
vicinity of a second order phase transition one finds $\alpha_{2}>0$ and
$\alpha_{4}>0$ in the restored phase while $\alpha_{2}<0$ and $\alpha_{4}>0$
holds in the broken phase in this context. Using Eq. (64) we find the
transition temperature at $\mu=0$

$\displaystyle T_{\mathrm{C}}=\tfrac{\mathrm{e}^{\gamma}}{\pi}\simeq
0.566933\,.$ (67)

The curvature $\kappa$ of the second order phase boundary
$T_{\mathrm{C}}\big{(}\tfrac{\mu}{T}\big{)}=T_{\mathrm{C}}\,\big{[}1-\kappa\,\big{(}\tfrac{\mu}{T}\big{)}^{2}+\ldots\big{]}$
at $\mu=0$ can be computed using Eq. (64) and is given by $\kappa\approx
0.241671$.

The critical end point is located at the intersection of the $\alpha_{2}=0$
and $\alpha_{4}=0$ lines, cf. Refs. [59, 197]. Using Eq. (65) for $n=2$ we
determine $\tfrac{\mu_{\mathrm{CP}}}{T_{\mathrm{CP}}}\simeq 1.910669$ from the
only root of $\mathrm{DLi}_{-2}$, which is
$\mathrm{DLi}_{-2}(\tfrac{\mu_{\mathrm{CP}}}{T_{\mathrm{CP}}})=0$. Having the
ratio $\tfrac{\mu_{\mathrm{CP}}}{T_{\mathrm{CP}}}$, we determine

$\displaystyle(\mu_{\mathrm{CP}},T_{\mathrm{CP}})\simeq(0.608221,0.318329)\,,$
(68)

using Eq. (64). Above $T_{\mathrm{CP}}$ we have the second order phase
transition with $\alpha_{4}>0$, while below $T_{\mathrm{CP}}$ we have the
spinodal region with $\alpha_{4}<0$ and $\alpha_{6}>0$ in the vicinity of the
critical point. $\alpha_{2}>0$, $\alpha_{4}<0$ and $\alpha_{6}>0$ allows for
three local minima in the Ginzburg-Landau potential (62) when considered up to
$n=3$ which indicates a first order phase transition below $T_{\mathrm{CP}}$.
At the critical point ($\alpha_{2}=\alpha_{4}=0$ and $\alpha_{6}>0$) the non-
trivial minima from the ($\alpha_{2}>0$, $\alpha_{4}<0$)-scenario and
($\alpha_{2}<0$, $\alpha_{4}>0$)-scenario merge in $\sigma=0$. The
$\alpha_{2}=0$ line below $T_{\mathrm{CP}}$ signals the appearance of a local
minimum at $\sigma=0$ and is called the left spinodal line ending in
$\mu_{\mathrm{L}}$ at $T=0$, in accordance with the results at $T=0$ discussed
earlier in this subsection. The location of the first order phase transition
and the right spinodal line, indicating the disappearance of multiple local
minima, are not accessible with Ginzburg-Landau expansion around $\sigma=0$
and have to be determined by numerical minimization of renormalized grand
canonical potential per spatial volume (58). Phase boundaries, the
$\alpha_{2}=0$ and $\alpha_{4}=0$ lines and the spinodal region are displayed
in Fig. 2.

###  V.4 UV initial condition for the effective potential for RG flows with
and without bosonic fluctuations

Next, we again turn to the UV initial condition for the effective potential
and continue the discussion of Sub.Sec. III.4. For practical calculations
within the FRG framework, where the flow is not integrable analytically, we
can not initialize the RG flow directly at $\Lambda=\infty$. But rather we
have to choose some very large and finite $\Lambda$222222$\Lambda$ should be
significantly larger than the internal scales (here mainly $h\sigma_{0}$ and
$\partial_{\sigma}^{2}U(\sigma_{0})$) and external scales (here $\mu$ and
$T$). For a detailed discussion, we refer to Refs. [176, 110]., where we
specify the initial values for the flow via
$\bar{\Gamma}_{\Lambda}[\Phi]=\mathcal{S}[\Phi]$.

From the previous discussion, it is obvious (at least in the mean-field
approximation), how to specify the initial condition for a numeric solution of
Eq. (40) or rather Eq. (32) with $Q(t,\partial_{\sigma}u)=0$. Via our
discussion on asymptotic freedom, we were able to eliminate $g^{2}$ from the
initial condition (42) in favor of the UV cutoff $\Lambda$ and the combination
$h\sigma_{0}$. Hence, we can simply use Eq. (49) as the initial potential at
some large scale $\Lambda$. Initializing the (numeric) mean-field version of
the RG flow (32) with the $\sigma$-derivative of (49) and an arbitrary value
for $h$ at some scale $\Lambda\ggg h$, one always finds that the IR minimum in
vacuum is located at $\sigma_{0}$. Consequently and w.l.o.g. we can rescale
all dimensionful quantities in terms of $h$ and express the dimensionless
field space variable $\sigma$ in multiples of $\sigma_{0}$. On the level of
the equations, this amounts to setting $\sigma_{0}=1$ and $h=1$. Other choices
for $h$ and $\sigma_{0}$ correspond to different renormalization conditions,
but all results are unique and can be transformed into each other via simple
rescaling – as already mentioned. However, we still have to ensure that other
IR observables do not depend on the UV cutoff $\Lambda$, which is realized by
choosing $\Lambda$ much larger than all internal and external model scales.
This shown numerically in Sec. V.5 and discussed in the context of RG
consistency.

Including bosonic quantum fluctuations, it is less obvious, how to choose the
initial condition for the RG flow, that means how to choose a meaningful value
for $\frac{h^{2}}{g^{2}}$ in Eq. (26). When performing calculations at finite
$N$, each individual choice of $N$ represents a single model on its own.
Hence, even if there is symmetry breaking for the vacuum flow equation (24)
including bosonic quantum fluctuations, the IR physics for different $N$ is
not necessarily directly comparable. Thus, setting a unique renormalization
condition for all $N$ in the IR, like fixing the position of the IR minimum
and/or the IR curvature mass by tuning the various UV-initial conditions, is –
to the best of our knowledge – not useful.

We think that it is natural to use exactly the same initial condition for RG
flows, namely Eq. (49), with and without bosons and to fix the renormalization
condition for the bosonic RG flows in the UV. This might be counter-intuitive,
because in a lot of RG studies for effective models of strongly correlated
systems, the physics is fixed in the IR and the UV initial condition is tuned,
in such a way that the RG flow ends up with an IR effective action having
observables compatible with desired numerical values. Actually, our situation
is much closer to what is done in FRG studies for QCD, see, e.g., Refs. [201,
202, 203, 204, 205, 206, 207, 208, 209, 210] as well as the lecture notes
[111] and the recent review [109] – another asymptotically free QFT. Also
there, the asymptotic freedom, which is derived from the one-loop beta
function, is used to fix the gauge-coupling at some extremely large $\Lambda$
in the UV.

Overall our formal arguments for choosing Eq. (49) are as follows. Using the
same UV initial condition for all finite and infinite $N$ allows for a direct
comparison of calculations at different $N$. This is the case, because the
$\tfrac{1}{N}$-rescaled UV potential (49) always describes a theory of $N$
asymptotically free fermions and a $\sigma$ mode that decouples from the
system at the UV initial scale. Choosing $\Lambda\ggg h$ naturally leads to a
large curvature mass $\partial_{\sigma}^{2}U(t,\sigma)$ and suppression of
fluctuations of the $\sigma$ mode in the UV, which can be directly seen in Eq.
(49). This argument was already brought up in Sub.Sec. III.4. On a formal
level this can be seen by inspecting the fluid dynamic formulation of the RG
flow equation (32) and especially its bosonic contribution in terms of a
highly non-linear diffusion equation (35). One finds that the large curvature
mass $\partial_{\sigma}^{2}U(t,\sigma)=\partial_{\sigma}u(t,\sigma)$ in the
propagators $\tfrac{1}{E_{\sigma}}$ yields a small diffusion coefficient
$D(t,\partial_{\sigma}u)$ and therefore a suppression of the diffusion along
field space – the bosonic contribution to the RG flow.

Though, the more drastic argument, why bosonic fluctuations are actually
totally absent in the UV, is a fundamental property of all diffusion equations
of type $\partial_{t}u(t,\sigma)\propto\partial_{\sigma}^{2}u(t,\sigma)$.
Independent of the finite diffusion coefficient the term
$\partial_{\sigma}^{2}u$ vanishes exactly for spatially linear
$u(t,\sigma)\propto\sigma$ and the diffusion and dynamics stops. This is
rather natural, if one considers diffusive processes from our everyday life,
e.g., heat conduction. In the context of this work, it follows from the
quadratic UV potential (26) that $u(0,\sigma)\propto\sigma$ and consequently
the contribution from the $\sigma$ mode to the RG flow vanishes exactly.
Bosonic fluctuations will be suppressed as long as the fermionic source/sink
contributions to the RG flow do not alter the linear shape of $u(t,\sigma)$,
which is approximately the case until $k^{2}(t)\approx(h\sigma)^{2}$ for small
$\sigma$. We conclude that the RG trajectories in theory space for RG flows
including fermions and bosons at finite $N$ will approximately follow the
mean-field RG trajectories for infinite $N$, as long as the UV initial
potential is quadratic in $\sigma$. This behavior is indeed observed in our