Shocks and quark-meson scatterings at large density

We discuss the phase structure of the two-flavour quark-meson model including quantum, thermal, density and critical fluctuations with the functional renormalisation group. This study combines two technical advances in the literature, that are also chiefly important for the quantitative access of the phase boundary of QCD at large density or baryon chemical potential. Specifically we allow for the formation and propagation of shocks as well as a fully self-consistent computation of the order parameter potential for chiral symmetry breaking.

We discuss the phase structure of the two-flavour quark-meson model including quantum, thermal, density and critical fluctuations with the functional renormalisation group. This study combines two technical advances in the literature, that are also chiefly important for the quantitative access of the phase boundary of QCD at large density or baryon chemical potential. Specifically we allow for the formation and propagation of shocks as well as a fully self-consistent computation of the order parameter potential for chiral symmetry breaking.
Chiefly important are the introduction of a Fierzcomplete basis of four-quark scattering vertices as well as quantitative access to order parameter potentials for homogeneous and inhomogeneous condensates. The latter allows us to discuss the eminently important question of the location of phase transition lines, that of the symmetry breaking pattern and the order of the phase transitions. It has been shown in the past decade that functional QCD flows towards QCD-assisted low energy effective models for energy scales below 1 GeV, for a detailed discussion see in particular the recent works [12,17]. With dynamical hadronisation [33][34][35][36] the LEFT is the (Polyakov-loop enhanced) quark-meson model (QM-model), or more generally the quark-hadron model. For recent fRG-works with the (P)QM model on the phase structure of QCD revelant for the present work see e.g. [16,, for a recent overview see [17]. This emergence of LEFTs from first principle QCD flows is well understood and quantitatively explored in the vacuum, see [34,[62][63][64][65][66]. It entails that the infrared critical dynamics is dominated by the low energy fluctuations of quarks and hadrons. For small baryonchemical potentials, µ B /T < ∼ 4, the relevant hadronic degrees of freedom are simply the pseudoscalar pions and the sigma mode, see [9,12,14,23,24]. In turn, for baryon-chemical potentials µ B < ∼ 4/T the situation is less clear, but we expect sizable diquark contributions, see [14].
In the present work we make significant steps towards such a quantitative control of the phase structure of high density QCD within functional methods. It combines two systematic advances in the past years: The first one was the development of self-consistent approximations for the computation of order parameter potentials, [47]. The second one was the development of a numerical approach for solving flow equations that also enables us to discuss discontinuities in the flows such as shocks that are potentially relevant for the correct description of first and second order phase transitions, [67]. Within this approach we compute the phase structure of the quark-meson model (QM-model) at finite temperature and density.
An important benchmark is already provided in the large N f -limit with an infinite number of flavours. It is argued that within an 't Hooft-type limit we can mimic the two-flavour QM-model well (or any other flavour), and in particular reproduce well it's non-universal properties such as the location of the phase boundary. Moreover, in this limit the numerical approach with discontinuous Galerkin set-up in [67] is fully developed and we have a quantitative access to the shock-development and propagation. The respective results are compare with the currently most advanced approximation (including shocks) to the self-consistent approximation including the order parameter potentials in [47] for the N f = 2-flavour quark-meson model. The results include also the regime µ B /T > ∼ 4. In this regime the current model has to be augmented with a diquark channel which is done in a forthcoming work. Still, the present work is a necessary and important study also in this regime.

II. QUARK MESON MODEL
The quark-meson model describes the dynamics of quarks and mesons at low energy. Within functional QCD this low energy effective theory (LEFT) emerges arXiv:2102.01602v1 [hep-ph] 2 Feb 2021 naturally from the momentum scale flow of the theory at momentum or cutoff scales k < ∼ 1 GeV, [12,34,[62][63][64][65][66]68]. In this regime the gluonic degrees of freedom decouple from the dynamics due to the gluonic mass gap in QCD, for a detailed discussion see [12,17,65], for a discussion of the emergent LEFT see [66].

A. Emergent LEFTs and their range of validity
The key ingredient for this emergence is the scaledependent four-quark scattering, whose dynamics at large momentum scales is driven by a box diagram with a two-gluon exchange between quark currents. For the discussion of its low-momentum behaviour we restrict ourselves to the momentum-independent tensor structures, that is 10 tensor structures in two-flavour QCD and 28 (32) tensor structures in three-flavour QCD, the relevant cases for the discussion of the phase structure of QCD. It has been shown in [64,65] that in the vacuum the scalar-pseudoscalar channel is dominating the dynamics by far, both above and below the chiral symmetry breaking scale of k ≈ 500 MeV: switching of all other channels leads to negligible effects for most physical observables. Moreover, in [14] is has been shown for N f = 2-flavour QCD in the chiral limit, that qualitatively this dominance persists up to large densities or chemical potentials, µ B /T c (0) ≈ 6, where T c (µ B ) is the chiral crossover or phase transition temperature at a given baryon chemical potential µ B . This highly interesting first dominance study in QCD is based on qualitative approximations, and a conservative error estimate leads us to µ B /T c (0) < ∼ 4 − 8 for the (total) dominance regime of the scalar-pseudoscalar channel.
This supports the computations in [12], where the phase structure of 2-and 2+1-flavour QCD was computed within a one-channel approximation (scalarpseudoscalar) to the Fierz-complete tensor structure for µ B /T (µ B ) < ∼ 6 or µ B /T (µ B ) < ∼ 4. Then, dynamical hadronisation takes into account multi-scattering events of the resonant channels (multi-scatterings of pions and the scalar σ-mode) that are relevant for the critical dynamics in a regime with second or first order transitions. In summary we estimate the reliability regime of the present approximations in functional QCD (see also respective considerations in DSEs ( [9,23,24]) to be The critical end point (CEP) computed both within the most recent fRG-computations, µ B /T (µ B ) = 5.59 from [12] and DSE, µ B /T (µ B ) = 5.54 from [23], for the physical case of 2+1-flavour QCD agree well, which sustains the respective reliability of these estimates. Still it is not within the regime of quantitative reliability of the current approximation. Consequently, (1) entails that for a quantitatively sound prediction of the CEP the current approximation to the full first principle QCD-flow has to be improved systematically in two directions for chemical potentials µ B /T (µ B ) > ∼ 4: First we need to include at least the dominant tensor structure at large densities, the cscor diquark-channel. This extension will be considered elsewhere. Second the self-consistent computation of the order parameter potential set-up in [7] is required. This is done in the present work within a recently developed numerical approach that also allows the inclusion of the formation and propagation of shocks, [67]. In this section we briefly recapitulate the fRGapproach to the (Polyakov-enhanced) Quark-Meson model (QM-model). The inclusion of the dynamical mesons as low energy effective degrees of freedom has to be seen as an efficient and convenient book-keeping device for the respective resonant interaction channels. In particular, this substitutes the rather tedious inclusion of the resonant parts of the higher-order scattering processes of quarks. Still, if used on a quantitative level, even for large UV-cutoff scales its effective action does not reduce to a simple local classical actions. For more details and in particular its quantitative properties as an emergent low energy theory in QCD we refer to [12,[63][64][65]. Validity checks, benchmarks and bounds in comparison to QCD have been provided in [66].
As discussed before, in the present work we restrict ourselves to a globally rather qualitative approximation to the effective action. Here, we are predominantly concerned with the quantitative access to the effective potential of the chiral order parameter. The systematic inclusion of the present quantitative setup within functional QCD flows is straightforward due to the modular nature of the fRG-approach and will be considered elsewhere.
The scales of the present LEFT are gauged by the pion decay constant in the chiral limit. We use f π,χ = 88 MeV and measure all other scales with these units.
We choose the UV-cutoff scale of the QM-model as Λ = 650 MeV. We consider this to be a good compromise between integrating-out as many momentum-fluctuations as possible and stretching the validity-bound of the LEFT. The momentum fluctuations with momentum scales k ≤ Λ are included with the functional renormalisation group (fRG). This approach has been used intensively in the past 25 years for the inclusion of low energy dynamics of the QM-model. For the setup of the flow equation for the effective action, and the derivation of the respective flow equations for (fielddependent) couplings we refer to the fRG-reviews see [17,[69][70][71][72]. Applications relevant for the present work can be found in [14,47,73], the derivations and flows for the present approximation can be found in [47].
The effective action Γ k [q,q, φ] of the N f -flavour QM-model is used in the following approximation, with τ µ being related to the Pauli matrices, τ = 1/2(1, iγ 5 σ), and the quark-meson coupling incorporates the SU(2) ∼ = SO(3) symmetry of the pseudoscalar subgroup. The O(4)-scalar field φ and the respective O(4)-invariant ρ are given by In (2) we have also x = 1/T 0 dx 0 d 3 x as an abbreviation for the finite temperature spatial integration.
We emphasize that the the Yuakwa coupling h k (ρ) is considered fully field-dependent. It multiplies the O(4)-invariant operatorq τ φ q, hence it only depends on the O(4)-invariant ρ.
The field-dependence of the Yukawa-coupling takes into account higher-order point-like scatterings of the resonant scalar-pseudo-scalar channels with the quark-anti-quark pair. The inclusion of these processes is necessary for a fully consistent zeroth order derivative expansion, and has been introduced in [47]. For further works in Yukawa models with fielddependent Yukawa coupling see [57,[74][75][76][77]. This is easily seen by performing a perturbative one-loop computation within the QM-model. Then, the quark loop with h(ρ) contributes to the full effective potential. Of course higher terms in the derivative expansion also contribute to the effective potential, but the Yukawa-term contains no derivatives. Accordingly, its full field-dependence should be accounted for in a consistent lowest order derivative expansion.
Finally, the scalar effective potential V k (ρ) takes into account the remaining part of the higher orders scatterings of the mesons. The linear term introduces explicit chiral symmetry breaking (finite current quark masses). Evidently, it drops out on the right hand side of the flow equation and does not run. Consequently, the full flow and hence the full effective potential V k does not know anything about explicit chiral symmetry breaking, and we do not consider it any further.
The next systematic step beyond the zeroth order derivative expansion would be the inclusion of wavefunction renormalisations Z q (φ), Z φ (φ) for quarks and mesons. This can be done either fully field-dependent (1st order derivative expansion) or field-independent (usually called LPA'). The latter approximation has been used in [47] together with a field-dependent Yukawacoupling. While technically in reach, we have chosen to drop these terms for the sake of concentrating on the quantitative discussion of the full effective potential. Hence, with (2) we assume implicitly, The flow equation for the complete set of couplings, h k (ρ), V k (ρ), and wave function renormalisations, can be found in [47]. We use the same setup here, including the choice of regulators, three-dimensional flat or Litim regulators, [78].
For the effective potential we simply evaluate the flow for Γ k [q,q, φ] for constant scalar fields φ and vanishing quark fields, q,q = 0. This leads us to with ρ-dependent quark-and meson-dispersion relations, and the RG-time t = ln k/Λ. The RG-time involves a reference scale in the logarithm, which we have set to be the initial scale. The masses m q , m φ are obtained by evaluating the respective two-point functions at constant fields. Note that m q , m φ are the curvature and not the pole or screening masses of quarks and mesons, for respective definitions and discussions see [46]. The meson curvature masses are defined with and hence are curvature-coefficients of the effective potentials. In turn, the quark mass is proportional to the Yukawa-coupling, It is left to discuss the flow equation for the fielddependent Yukawa coupling, for details we again refer to [47]. We can project the flow for Γ k onto the Yukawa coupling h(ρ) by evaluating the quark two-point function at vanishing quark and pion fields, q,q, π = 0, and constant σ. With (2) we arrive at where we have dropped the momentum conservation δ p,p in the last line with δ p,p = (2π) 4 δ(p − p ) in the vacuum. (9) reflects the fact that the Yukawa term simply is the ρ-dependent mass quark term, m q,k (σ) = h k (ρ)σ/2. Accordingly, the flow of the Yukawa coupling h k (ρ) can be derived from that of the scalar part of the quark twopoint function: it is simply σ/2 ∂ t h k (ρ) as ∂ t c σ = 0 by definition. Thus we get, In (10) we have used that ρ = σ 2 /2 for π = 0. The flow (10) is depicted in Figure 1. From (10) and the approximation (2) we finally get, (1,1) (m 2 q,k , m 2 π,k ; T, µ) with v d−1 = 1 The threshold functions l The the first and second lines from (11) are contributions of the first two diagrams with mixed fermionic and bosonic loops in Figure 1, whereas the third and fourth lines are that from the bosonic loop with the four-vertex.

C. The Large-N Limit
Most of the numerical results in the present work are obtained in the large-N limit of these equations, as it simplifies the numerical treatment significantly: it eliminates the σ-loop in the flow equation, and hence the second derivative terms in the sigma meson mass term. We are left with only the pion loops as well as the quark loop. The pion loops constitute the flow of a pure O(N )-theory in the large N -limit as considered in [67] with discontinuous Galerkin methods. Such a non-linear first order system is solved using approximate Riemann solvers. These solvers rely on the assumption that the solution is dominated by one strong wave, for more details see Appendix A. This assumption holds if the flow is dominated by contributions of the pion and quark loops, which is always ensured in the large-N limit.
This simplification is also helpful when considering systems of multiple differential equations. However, it will also be instructive to make a comparison between both the finite N f case and the large N f limit in the case with constant Yukawa coupling. Moreover, we can simulate the physics case, N f = 2 and N c = 3, with a suitable chosen large N -limit: To begin with we keep the ratio of color and flavour fixed to that of the N f = 2 quark-meson model in QCD, With (13a) we keep the flavour-colour balance of QCD intact. This ensures that the contributions of the quarkloop are not suppressed by 1/N f . Moreover, the flavourcolour ratio is certainly of crucial importance for e.g. the question of the existence and size of a quarkyonic phase. Finally, we consider the following generic rescaling of ρ, V k (ρ) and h k (ρ), The factor N π in (13b) is introduced to simulate the flows of a quark-meson model with N π pions instead of one sigma meson and three pions. Both cases are relevant for the physics of two-flavour QCD or the two-flavour QM-model. In the chirally symmetric phase for large temperatures and cutoff scales, the pions and the sigma are degenerate on-shell at ρ = 0. The second derivative term vanishes 2ρV (ρ)| ρ=0 = 0, and the on-shell σ-propagator agrees with the pion one, and the (on-shell) flow equation resembles that with four pions.
In turn, in the broken phase, the σ-mode develops a mass and quickly decouples from the dynamics of the system. Then, the (on-shell) dynamics of the theory is driven by the three (massless or light) pions. From previous fRG investigations of the quark-meson model as well as QCD we know that the mesonic dynamics in the broken phase is of sub-dominant importance for not too large chemical potential. This suggests that the N = 4 case should mimic the two-flavour case best. A full discussion of the comparison is provided in Section IV B 1 and Section IV B 2.
With the limit N f → ∞ and (13) we derive the flow equations for large-N f Yukawa coupling, h lN k (ρ), and effective potential, V lN k (ρ), and This concludes our derivation of the set of flow equations solved in the present work: the system of flows at finite N is given by (5) and (11), those in the large N -limit are given in (14) and (15). Numerical results for both systems will be presented in Section IV, the discontinuous Galerkin setup, with which the numerical results are achieved, are discussed in the next section.

III. DISCONTINUOUS GALERKIN METHODS IN THE FRG
Most of the flow equations introduced in the previous Section II can only be solved numerically. In the present work we use Discontinuous Galerkin methods (DG-methods), which have been introduced to the fRG in [67] by the example of the large-N limit in an O(N )model. In contrast to the set of flow equations discussed in the present work for the QM-model, the flow equation for the effective potential in the large N -limit of the O(N )-model is given by an hyperbolic equation of first order that can be written in a conservative form. This type of equation is well studied and understood. Many different numerical schemes were developed to obtain a stable solution in a weak sense see e.g. [79].
In the present case, the flow equations (14) and (15) in the large-N limit are not conservative anymore. Indeed, for a constant Yukawa coupling the QM-model can be understood as a driven O(N )-model, where the driving force is provided by the quark loop. If this approximation to the effective action is upgraded to one with a cutoffdependent quark two-point function, there are backcoupling effective from the meson loop into the quarkloop, and the driving-force is not (fully) independent anymore.
In any case the system of differential equations ceases to be conservative. For the non-conservative hyperbolic problem, like the large N f equation with running Yukawa coupling, the concept of a weak solution was introduced relatively recent in [80,81] and applied in a context of Finite Volume and Discontinuous Galerkin schemes [82][83][84][85][86][87][88][89] for multiple physical systems. Hyperbolic equations in non-conservative form occur rather frequently in modeling physical system, as example viscous relativistic hydrodynamic equations are of this type [90][91][92][93][94][95][96][97] and recently also a formulation of general relativity has been solved in this formulation [98] highlighting more advanced stability properties. Moreover, in case of finite N f it also contains a diffusion term that originates in the σ-loop. Thus, on the technical level, the present work represents a non-trivial extension of [67]. The different extensions are discussed in Section III A (nonconservative systems) and Section III B (diffusion terms). With respect to these extensions the present work should be considered as a first step to the full implementation of DG-methods for non-conservative systems, more details and further extensions will be considered elsewhere. For a more detailed introduction to DG-methods see also [99].
In the context of the FRG, further work concerning the inclusion of higher order derivatives and non-conservative formulations has been achieved in [100], which will be published soon.
Pseudo-spectral methods are an integral part of DGmethods. They are applicable to FRG equations in the absence of shocks, and have been used successfully in e.g. [101][102][103][104][105].

A. Non-conservative flux equations
In this section the extension of DG-methods for the fRG to non-conservative flow equations is set up. To this end we consider a system of differential equations of the form, (16) where u = (u 1 , u 2 ) T and i, j ∈ {1, 2}. The s i are source terms and f i conservative fluxes. In (16) we also allow for non-conservative terms a i . In the full quark-meson model the flux is additionally separated into a convective and a diffusive contribution depending also on ∂ ρ u i . We note that the splitting into conservative and non-conservative terms is not unique in these equations. More details on the numerical treatment can be found in Appendix A, whereas an evaluation of the convergence properties of the scheme is performed in Appendix B. In the following (5), (14) and (15) are reformulated to fit (16).

Flow of the effective potential
Equation (5) and (14) are rearranged to ensure that they have the conservative form required by (16) such that DG methods are applicable in a purely conservative formulation. Similarly to [67] we observe a non-linear dependence of the potential V k (ρ) on its derivative with respect to the field expectation value ∂ ρ V k (ρ) and in case of (5) also on its second derivative. The dependency on the first derivative can be eliminated by introducing it as a new variable, which coincides with the pion mass squared: Since (14) is not dependent on itself we can simply take a ρ derivative, which turns u k (ρ) into a conserved quantity which is fit for DG schemes. This procedure is also applied to the QM-model. In this case we need an additional expression for the second derivative of the potential in the sigma mass in (8). We obtain this expression by taking another ρ derivative of the polynomial basis functions φ n introduced in Appendix A:

Flow of the field-dependent Yukawa coupling
The flow of the Yukawa coupling at finite N is given by a highly non-linear equation of second order. Since it can not be made to fit the form given in (16), it is not solved within the introduced framework. However, the expression simplifies significantly in the large-N limit and (15) can be written to suit the formalism. Equation (15) is rewritten in terms of the quark mass squared m 2 q (ρ), as we are primarily interested in physical observables. Thus, a new variable is introduced, We refer to Appendix D for the explicit calculation. Introducing this new variable proves to be very helpful for the computation. Appendix E explains how the ambiguity in splitting the conservative and nonconservative contributions to the flux are used to accommodate boundary conditions. For completeness the final form of the equations is stated, This version of the equation has the advantage that the non-conservative product is rather small in comparison to the conservative part.

B. Finite N equations
For finite N, the equation is parabolic. Apart from the convection term (Goldstones), there is also a diffusion term that arises form the σ-loop. The equations for finite N are highly non-linear and of second order, therefore we refrain from considering the field dependency of the Yukawa coupling. Schematically the flow equation of the potential is written as The weak formulation and the stability of this type of equation has not been fully understood until now. The presence of the diffusion modifies the numerical flux significantly. However, in the convection dominated regime, and in the absence of a discontinuity, it is possible to neglect this diffusion numerical fluxes and formulate the Discontinuous Galerkin method for this equation as follows, where a i,h , s i,h and f i,h are computed form the field u i and their local approximation of the derivative, no other numerical fluxes are introduced into the numerical scheme. The absence of numerical fluxes for the extra derivative present in the equation correspond to the assumption of continuity of this field and the DG scheme somehow reduces to a pseudo spectral method. This approximation is acceptable, whenever the flow is rather smooth and no shock or rarefaction wave are generated during the simulation. In turn, this scheme will fail in the vicinity of a first order phase transition. There we expect shock-formation and propagation in the flow equation.
In conclusion, for the rest of the phase diagram the present approximation can be considered as a sufficiently accurate solution of the flow equation due to the local high order accuracy of the DG scheme.

IV. RESULTS
In this section we present and discuss our results for the phase structure of the QM-model in the different approximations. We chose our initial conditions such that the quark mass and the pion decay constant reproduce physical values in the vacuum, this is discussed in detail in Section IV A. Note however, that its not the main objective of the present work to produce quantitatively reliable results, the vacuum bench marking simply facilitates the understanding of our results. In the present work we rather focus on the qualitative behavior of the matter sector of QCD at large densities, quantitatively reliable results require full QCD flows and will be considered elsewhere. Such a set-up entails, that while we compute and present a phase structure at large densities, our present low energy effective theory gradually looses predictability for larger µ B /T . Such an estimate in functional QCD leads to a predictability bound of µ B /T < ∼ 4, if the currently existing state of the art computations are combined, [12,23,24] and estimates for missing channels and effects are considered as well, [12,14,106]. In the present class of low energy effective theories (QM, NJL-type, PQM, PNJL), a respective estimate leads to µ B /T < ∼ 2. We first present results within an approximation where only the effective potential depends on the cutoff  The scales are fixed with the pion decay constant in the chiral limit f π,χ = 88 MeV, that is m σ /f π,χ ≈ 3.603 and m q /f π,χ ≈ 3.532. In the present approximation we have f π,χ = σ 0 , and in the model the dimensionless ratios are simply m σ /σ 0 and m q /σ 0 . In the chiral limit we also have m π = 0.
scale, the local potential approximation (LPA), for both, the finite N f and the large-N f limit in Section IV B. They serve as a benchmark for the more advanced approximations discussed in Section IV C. Additionally, the results in Section IV B also serve as benchmark for results in the literature within the QM-and Polyakov loop-enhanced QM models, in particular at large density, where DG-methods or similar numerical approaches are mandatory for reliable results.
In Section IV C we present results for the coupled system of effective potential V k (ρ) and Yukawa coupling h k (ρ) in the large-N f limit. This investigation allows us to solidify the results in [47] concerning the flattening of the quark mass m q (φ).
Lastly, the technical advances made here readily carry over to first principle QCD within the fRG, as discussed in the introduction, they are one of two missing ingredients for reliable predictions of the QCD phase structure for µ B /T > ∼ 4.

A. Initial Conditions
We initiate the flow at a cutoff scale k = Λ ≈ 0.650 GeV with the classical action of the QM-model. Then, the parameter in the initial effective action Γ Λ is the φ 4 -coupling in the classical potential, as well as the Yukawa coupling h Λ . For the sake of simplicity we use a initial meson quark mass, m 2 φ = 0. All our scales are measured in the pion decay constant in the chiral limit f π,χ = 88 MeV. Within the present approximation of the QM-model we have f π ≈ σ 0 , and hence we define σ 0 = 88 MeV. Then, the two model parameters λ Λ , h Λ are fixed such that they lead to a 'physical' constituent quark mass 1/ √ 2hσ 0 , and a 'physical' mass of the sigma resonance, m σ . The parameters for the couplings of the effective theory and their relation to physical observables are summarised in  Table I. The dimensionless ratios in the models at k = 0 are given by and follow with the initial parameters in Table I.

B. Results for the effective potential with constant Yukawa coupling
In this section we compare the numerical results of the physical case with N f = 2, and in the large-N f limit with three and four degrees of freedom. This is done in LPA, where we solve the flow equation for the effective potential, (14). We first discuss the asymptotic regimes: vacuum, large temperatures, and large chemical potential, Section IV B 1. Then we show that the chiral phase transition, or rather its non-universal properties, agree quantitatively for all models, Section IV B 2. The shock-development at large chemical potential is discussed in Section IV B 3. Finally, we compare the phase structure for all three cases in Section IV B 4.

Asymptotic regimes
We begin with an evaluation of the asymptotic regimes: the vacuum with µ q , T = 0, large temperatures with µ q = 0, and large chemical potentials with T = 0. For these cases we show the field-dependence of the pion mass m π,k (ρ) = ∂ ρ V k (ρ), see (17) for different cutoff scales. This resolves the effective potential, obtained from an integration over ρ, in terms of a physical observable.
For the numerical solution of (14) we use an interval of ρ ∈ [0, 0.02] GeV, which is expanded in K = 80 cells with polynomials of order N p = 2. The length of the interval is The corresponding fit parameters to (25) can be found in Table II. chosen such that it includes all relevant phenomena: the flux at the outer boundary is very small. This ensures the numerical convergence for the entire phase diagram.
The benchmark case is the vacuum, where the present models are anchored, see Section IV A. The fielddependence of the pion mass is shown in Figure 2. For the initial cutoff k = Λ, the pion mass is simply a straight line, m π,Λ (ρ) = λ Λ ρ, where the slope is the initial mesonic self-coupling, λ Λ . With decreasing cutoff scale, the pion mass develops a flat regime, which is related to the emergence of non-trivial minima ρ 0 = σ 0 /2 in the potential, indicating chiral symmetry breaking.
We also conclude, from the comparison of the pion masses in the different models, that the cutoffdependence of the physical two-flavour case is best mimicked by the large-N f limit with four degrees of freedom: most of the dynamics takes place in the symmetric regime or close to it. Technically, this regime is described with 1 + m 2 σ (ρ 0 )/k 2 ≈ 1 + m 2 π (ρ 0 )/k 2 , owing to the fact that the total mass of the respective modes is k 2 + m 2 π/σ . Finally, we determine the relative values of the pion decay constants in the chiral limit, f π,χ = σ 0 , in the different models. All scales are measured in the pion decay constant f π = 88 in the two-flavour case. The expectation value σ 0 or ρ 0 = σ 2 0 /2 follows from (7) as the maximal field value with m π (ρ 0 ) = 0. However, since we stop the numerics at a small but finite cutoff value, k min = 50 MeV, we extrapolate the expectation value σ 0 (k min ) to σ 0 (0) within a linear fit: we use data from 10 equally spaced RG-scales from k = 90 MeV to k = 50 MeV, and fit The self-consistency of this linear fit is checked by the perfect agreement of the linear fit with the data, see Figure 3. The respective values for σ 0 are given in Table II. The mass m σ of the scalar mode is extrapolated to k = 0 from the same data. Note however, that once   (25) to the zero point at 5 equally spaced RG-scales k = 90 MeV to k = 50 MeV. The error is computed from the grid resolution and the error to the fit parameters. The mass m σ is extrapolated from the same data points as the fit. The error of m σ is obtained analogously to the one of σ 0 , it might be underestimated due to the kink developing at σ 0 .
the kink enters the cell in which σ 0 is located, the precise determination of the derivative is difficult. Therefore, the flattening of the potential most likely causes an underestimation of m σ .
For large temperatures we safely stay in the symmetric regime and the mesons simply acquire a thermal Debye mass. This is seen in Figure 4a. In the symmetric regime we have four mesonic degrees of freedom in the twoflavour case. Consistent with our expectations, that the cutoff-dependence of the pion mass in the model with four degrees of freedom in the large-N f limit has the best agreement with the two-flavour case.
We close with a discussion of the large chemical potential asymptotics.
The respective pion mass (squared) is depicted in Figure 4b. The sudden increase of the pion mass for k < ∼ k on with k on ≈ µ in Figure 4b is related to the Silver Blaze property, [107], for the discussion in the fRG-approach see [50,108,109]. This property entails that correlation functions below the density onset are simply functions of p 0 ∓ iµ q for quark and anti-quark frequencies respectively. Accordingly, observables do not depend on the chemical potential for µ q < µ q,on , where µ q,on is the onset chemical potential. For µ q > µ q,on , the medium leads to deformations of the quark-meson scattering processes, comprised in medium meson-dispersions. In the presence of thermal fluctuations this onset is washed out with increasing temperature.
where for N F → ∞ we have β = 1/2, the mean field critical exponent. In turn, for the present LPA-study of the Yuakwa model with the O(4)-universality class we have used the three-dimensional spatial flat or Litim regulator, [78], for both quarks and mesons. See also Appendix H, (H1), (H2). This leads us to β ≈ 0.40, see [110] with [111,112] and in particular the recent work in the QM-model, [19]. Note, that more advanced approximations of the fRG provide β ≈ 0.39 consistent with conformal bootstrap and Monte-Carlo results. For a recent compilation see [113], for the QM-model see [19] that also includes an investigation of the Z 2 -universality class.
The respective scaling regimes are already very small in the O(4)-model and even shrink in the presence of the (driving) fermion loop, see the discussion in [19]. While possible, we do not aim at a precision estimate of critical exponents here, as we focus on the location of the phase boundaries. Accordingly, we have simply checked the consistency of the scaling law (26) with β ≈ 0.4 (N f = 2) and β = 1/2 (N f → ∞) for small reduced temperatures 1−T /T c → 0 − . This also allows us to determine the respective scaling regimes. Consistent with the observation above that they should be even smaller as the already small scaling regime in O(N )models we find scaling for Moreover, a scaling fit with (26) in the regime (27) allows us to determine T c as well as the prefactor c cr . We see from Table III, that the two-flavour critical temperature agrees well with the large-N f limit with four degrees of freedom. This is expected from the theoretical analysis and our results on the asymptotics in Section IV B 1. This good agreement extends to the full temperature dependence, as can be seen from Figure 5a. In turn, the order parameter from the large-N f limit   (26) to the mean field expectation values in Figure 5a which are underlined by the transition scaling. The error in the data is expected to be higher, since the numerical precision is limited by the grid resolution. An exact reconstruction of the zero-crossing is not possible.   with three degrees of freedom seemingly shows a slightly different behaviour.
However, the two large-N f models are obtained by a simple rescaling of the fields and hence are identical to each other. They can be mapped onto each other by the relative rescaling. Put differently, the temperature dependence of the order parameters should agree if plotted in dimensionless units, σ(T )/σ(0) and T /T c . This comparison is shown in Figure 5b: as expected, the temperature-dependence of the order parameter of the large N f models agree. More importantly, also the two-flavour case agrees quantitatively, though with small deviations. Trivially, the non-trivial critical scaling of the two-flavour case does not agree with the trivial meanfield scaling for N f → ∞, but the scaling regimes are very small, see (27).
In summary, the thermal properties of the models in the large N f -limit and the physics case N f = 2 agree very impressively.

Shock Development and First Order Phase Transition at High Densities
In this section we discuss the shock development and propagation at intermediate densities and very low temperatures. This is also used to discuss the first order regime.
In these computations we use a grid ρ = [0, 0.2] GeV and expand in K = 200 elements with a local approximation order of N p = 2. This finer grid is required for the shock resolution at low temperatures and chemical potentials close to the onset chemical potential. Indeed, the full resolution of some of the features in this regime (e.g. the precise location of the transition line in the absence of shocks) requires an even higher resolution. While technically possible, we have refrained from doing so, as the related aspects have been not in the main focus of the present work.
We first note, that the running of the pion mass stops quickly below the onset RG-time t on = ln Λ µ : the RG-flow is proportional to the Fermi-distribution. Hence it stops at k on at T = 0, for finite T the Fermi-distribution is softer but for small temperatures there is still a strong exponential suppression for k ≤ k on . This suppression leads to two competing effects at finite densities: • The onset amplitude, and thus m 2 π (ρ 0 ) in the symmetric phase, are linked to the suppression of the quark-contribution. This contribution dominates initially, but due to the constant Yukawa coupling it is quickly suppressed with k 5 .
• For field values with positive meson masses m 2 π,k (ρ) the meson loop in the flow is suppressed with k 5 . In turn, for negative meson masses the meson loop is suppressed with k 4 . Note also, that the flow increases with decreasing values of m 2 π,k , which is closely linked to the restoration of convexity. The mesonic flow contribution is reminiscent of the spreading of waves in hydrodynamics, where its value corresponds to the wave velocity: if we consider the solution m 2 π (ρ) as a wave packet, it flows with the RG-time in the direction of smaller field values with a ρ-dependent propagation velocity. The velocity of the solution and its effect on convexity is inspected more closely in Appendix C. GeV is plotted at different densities. We find a second order phase transition of the order parameter. The parameters of a scaling fit, see (28), are given in Table IV. of shocks and a first order phase transition at low temperatures.

The interplay of both effects leads to the creation
The increased propagation speed of negative modes is blocked by the slowed propagation of positive modes. The shock travels towards smaller field values during the RG-time evolution, but eventually freezes when the shock amplitude is too high. An illustration of this process can be found in Appendix F.
Naturally, the occurrence of shock development depends on the choice of initial conditions, specifically those that a trigger stronger dynamics of the system, for a respective discussion in the O(1)-model see [67]. With physical initial conditions we find shock development in the large N f limit with 3 DoF, whereas the dynamics for 4 DoFs and finite N f are not strong enough to generate a shock at finite temperature T > 10 MeV. This is an important observation: we have used the same initial conditions for all models, fixed within the twoflavour case. As discussed before, the two models in the large N f -limit only differ by a rescaling of the fields and parameters. Accordingly, they can be interpreted as the same model with different initial conditions, as we do not apply any rescaling to the initial condition. However, these changes are marginal, as can be seen from the small variation of the pion decay constants and σ-masses in the vacuum, see Table II. In conclusion the physical case is very close to the situation where shocks may form during the RG-time evolution. Whether or not this also occurs in QCD requires further investigation: (i) The embedding of the present model as part of the matter sector will lead to additional driving forces in the flow. This may be mimicked with   (28) to the shock positions in Figure 6.
a T, µ-dependent change of the initial conditions here. naturally these changes can go either way, they may support the shock development or soften it.
(ii) The additional diffusion terms in the finite-N f case, see Section III B, may structurally soften the RGtime evolution and remove any shock development. It is also unclear whether shocks develop in the presence of diffusion terms allow, as the diffusive flux counteracts the formation of a discontinuity.
The resolution of these aspects is crucial for an access to the QCD phase boundary at large chemical potential and low temperatures. This goes far beyond the scopes of the present work, and is subject of ongoing work.
In the large N f limit with 3 DoFs we can use the shock development at low temperatures for an accurate determination of the phase transition line. The shock position ξ final at k = 0 is extrapolated by fit, utilizing the exponential decay of the flows, We use the shock position at 6 equally spaced time steps between RG-times t = 3 and t = 3.5. We expect the same power law behavior as in [67] for the final shock position as a function of chemical potential, As an explicit example we concentrate on T = 10 MeV. The coefficients of the χ 2 -fit are provided in Table IV. From this second order phase transition we obtain an accurate estimate for the critical chemical potential of µ crit = 0.30460 ± 0.00033 GeV. The phase transition and shock positions are depicted in Figure 6. Shock formation occurs only in a relatively small area of the (T, µ) plane, being confined to intermediate densities 290 MeV < 360 MeV and small temperatures up to T = 20 MeV.
In the absence of a shock a very fine grid has to be used to pin down the phase transition line for low temperatures. In the present work we have simply narrowed down the location of the phase transition line for small temperatures T < ∼ 30 MeV to a small interval µ crit ∈ [270, 290] MeV.   For the resolution of shock formation a finer grid is required. We observe a formation of shocks around the first order phase transition in the large-N case at densities around µ = 0.3 GeV. In this area we expand in K = 250 cells to reduce oscillations and ensure numerical convergence. The flow is evaluated up to t = 4, which corresponds to a momentum scale of k = 0.001 GeV. The field expectation value σ is chosen as order parameter and evaluated as demonstrated previously for the approximate vacuum.
The result for the large-N limit is depicted in Figure 7. The crossover region is discernible by the color gradient that smoothly transitions between both phases, whereas in the first order regime a jump is clearly visible. Figure 7b also illustrates how shock development shifts the critical chemical potential to higher values.
The phase diagram for the finite N f case is given in Figure 8, the computation did not converge at high densities, this is further discussed in Appendix B 1. As discussed before, it will be interesting to see, how this regime changes, if the present model is embedded as part of the matter sector of QCD in full QCD-flows. This is subject of ongoing work.

C. Quark-meson scatterings in the Large-N limit
The discussion of the results in LPA with a constant Yukawa coupling have revealed a very intricate structure at about onset chemical potentials and low temperatures. In particular the occurrence of shocks is very sensitive to  We now present our results for the QM-model with a field-dependent Yukawa coupling in the large-N f limit with four DoFs, based on the combined numerical solution of the flows (14) and (15), as formulated in (20).

Dynamics in the vacuum
For the discussion of the dynamics in the vacuum, (20) is solved on a grid with varying cell sizes. A local approximation order of N p = 3 is used with K = 100 cells in ρ ∈ [0, 0.02]. Figure 9 depicts the solutions of the pion and quark masses in approximate vacuum, the pion mass in comparison to the case with constant Yukawa coupling. We can see that in approximate vacuum the pion mass remains unchanged for both models. An exponential fit is performed on the position of the zero point of ∂ ρ V k (ρ) using 5 equidistant RG-scales from k = 65 MeV to k = 25 MeV and we obtain This is consistent with the results in the previous sections with a constant Yukawa coupling (LPA). Consequently, this confirms previous findings in [47], that LPA or rather higher orders in the derivative expansion are a good approximation for vacuum QCD. Figure 9b depicts the field-dependent quark mass m 2 q,k (ρ) for different RG-times. As argued in [47], the quark mass flattens for meson fields ρ ≤ ρ 0 . The computation in the present work puts these conceptual and preliminary numerical findings on a sound numerical footing. In summary, at vanishing cutoff scale k = 0, this leaves us with a field-dependent quark mass m 2 q (ρ) with Note that while conceptually the field value ρ 0,q , below which the mass function flattens, has to agree with ρ 0 , the solution of the EoM, numerically this is not fully guaranteed. Hence, this provides a further consistency check of the present scheme. For performing the respective reliability check, we have determined the position of the kink by subtracting −ρ from the solution and taking the minimum. An exponential fit gives, which coincides within its error σ 0 in (29). The error is computed from the grid resolution and the error to the fit parameters. The quark mass in the flattened area in Figure 9b is computed as follows: one computes the average value of the quark mass up to the kink for the previously mentioned 5 RG-scales and performs an exponential fit. This leads to the physical quark mass, The linear ρ-dependence of m 2 q (ρ) in Figure 9b for ρ > ρ 0 entails, that the Yukawa coupling is constant with ∂ ρ h(ρ) ≈ 0, already observed in [47] for finite N f . The constant approximation for these field values works so well, that we can use the input Yukawa coupling h Λ to confirm σ 0 with the consistency relation σ 0 = m q (σ 0 )/h Λ . This leads us to σ 0,q = 86.01(24) MeV, where we take the error from the fit and the mean deviation in the flattened region.
In summary the present numerical analysis confirms quantitatively the conceptual and preliminary numerical analysis in [47]: In the broken phase the flattening-out of the field-dependent quark-mass m 2 q (ρ) is triggered by the flattening of the effective potential. In turn, in the symmetric phase the field dependent quark mass does not flatten-out and the quark mass vanishes on the solution of the equation of motion, ρ 0 = 0. Respective plots of the field-dependent quark mass m 2 q (ρ) and the pion mass m 2 π (ρ) for high temperature and density values are discussed in Appendix G. Importantly, apart from the flattening of the quark mass, they do not deviate significantly from the results in LPA. In particular this applies to their values of the equations of motion. Note however, that this is bound to change for finite N f .

Phase structure
As in LPA with a constant Yukawa coupling, we finally present our results on the phase diagram of the QM-model in the large-N f limit including quark-meson scatterings via a field-dependent quark mass or Yukawa coupling.
The computation uses the same resolution as above: N p = 2 and K = 80 cells in ρ ∈ [0, 0.02] up to the RGtime t = 4, that is k = 0.001 GeV. The result is shown in Figure 10. The computations did not complete the time integration for µ ≥ 0.3 GeV, respective upgrades are currently investigated.
As already discussed in the last section, secsec:DynVac, in the phase with chiral symmetry breaking the fielddependent quark-mass m 2 q is necessarily flattened for ρ ≤ ρ 0 . In turn, The quark mass function does not flatten in the symmetric phase, and the quark mass is found to be zero in the symmetric phase. Plots of the field dependent quark mass m 2 q ang pion mass m 2 π for high values of temperature and density can be found in Appendix G and do not significantly deviate from the results with a constant Yukawa coupling.
We close this section with a comparison of the phase structure in Figure 10 with that in LPA, Figure 7b, in the same setting: large-N f limit with four DoFs. While the phase boundaries do not change significantly, the crossover gets softened, if quark-meson scatterings are taken into account. This is clearly visible in Figure 11, where we depict the chiral order parameter as a function of temperature for different densities with µ q = 0, 100, 150 MeV. This entails, that the quark-meson scatterings considered here give sizable contributions to important observables measured in heavy ion collisions. First of all, fluctuation observables will be sensitive to such a widening of the crossover. These effects may be even more prominent at large chemical potential where the freeze-out line most probably deviates from the chiral crossover line. Moreover, it can also be deduced from Figure 11, that the quark-meson scatterings may have a sizable delaying effect on a possible critical end point. We conclude, that these scatterings have to be taken into account for a quantitative prediction for the location of the CEP. : Chiral order parameter σ 0 (T ) in the large-N f limit with four DoFs as a function of temperature and quark chemical potential. We compare the LPA results with quark-meson scatterings encoded in m q (σ) for different densities. The data is interpolated from the phase diagrams in Figure 7 and Figure 10.

V. CONCLUSIONS
We have presented a study of the QM-model with an emphasis on a quantitative access to order parameter potentials at finite chemical potentials. This allows us to discuss the eminently important question of the location of phase transition lines, that of the symmetry breaking pattern and the order of the phase transitions.
The present study combines two systematic advances in the past years: The first one was the development of self-consistent approximations for the computation of order parameter potentials, [47]. The second one was the development of a numerical approach for solving flow equations that also enables us to discuss discontinuities in the flows such as shocks that are potentially relevant for the correct description of first and second order phase transitions, [67].
Within this approach we have computed the phase structure of the quark-meson model (QM-model) at finite temperature and density. An important benchmark is already provided in the large N f -limit with an infinite number of flavours. We have argued that within an 't Hooft-type limit we can mimic the two-flavour QM-model well (or any other flavour), and in particular reproduce well its non-universal properties such as the location of the phase boundary.
Moreover, in this limit the numerical approach within the discontinuous Galerkin set-up in [67] is fully developed and we have a quantitative access to the shockdevelopment and propagation, even in the presence of non-conservative forces.
The present approach works very well except of a small regime at low temperature and onset densities. This is a merely technical problem and related upgrades of the present schemes are in development. Moreover, already for smaller ratios µ q /T close to the crossover line we have to also improve the current approximation of the matter sector of QCD. This follows already from [12,14]. The results there indicate the potential relevance of nontrivial meson dispersions as well as the diquark channel at larger chemical potentials, µ q /T > ∼ 4/3. Moreover, in the vicinity of a potential critical end point we also have to take into account the density channel, that mixed with the critical σ-mode. This is work in progress and we hope to report on it in the near future.

ACKNOWLEDGMENTS
We thank L. Corell, L. Kades, A Koenigstein, F. Rennecke, M. Steil and J. Urban for discussions. This work is done within the fQCD collaboration [114], and is supported by EMMI, the BMBF grant 05P18VHFCA, and is part of and supported by the DFG Collaborative Research Centre SFB 1225 (ISOQUANT) as well as by the DFG under Germany's Excellence Strategy EXC -2181/1 -390900948 (the Heidelberg Excellence Cluster STRUCTURES).

Appendix A: Implementation of DG-Methods in Dune
This section gives an introduction to the numerical and computational framework used to solve the equations (16). We made use of the DUNE (Distributed and Unified Numerics Environment) library [115][116][117][118][119][120], which is a modular toolbox for solving partial differential equations with grid-based methods. The library is an open source initiative to create a common interface for many different numerical methods and supports high performance computing .

Weak formulation and discrete problem
The system of equations is solved on a computational domain Ω h , which is composed of K disjoint elements, called cells, D k such that: For the purpose of the calculations in this work we used the Dune-grid YaspGrid, which is contained in the module dune-grid and allows for n-dimensional cubic grids and parallelised computation. In this paper a one dimensional grid is used. It represents the computational domain Ω h , with the grid-cells being disjoint intervals D k of possibly differing lengths, as discussed in Section III. In a more general formulation the domain Ω h would be given as an n-dimensional rectangular grid and the elements D k would be implemented as cubic grid cells.
The solution in each cell D k is approximated by where u k h (t, x) is the local solution in each cell and the index h denotes the approximation. The local solution in turn is then approximated by a polynomial of degree N = N p − 1 such that in each element D k : The local approximation u k h (t, x) is given by a modal expansion, where {q n } is a local polynomial basis with time dependent expansion coefficientsû k n (t). Thus the global solution consists of K local polynomial solutions of order N . The local approximaltion was implemented using the dune-pdelab module, specifically using the class QkDGLocalFiniteElementMap. In the one-dimensional case basis functions q n are given by the Legendre-Polynomials up to order N p .
For the purpose of higher dimensional computations the basis functions are taken from the polynomial space Q k of the Legendre-Polynomials.
In our calculations we use a locally defined weak formulation of the convergence requirement: The right hand side of the equation contains the standard numerical flux f * as well as an additional nonconservative flux term D.n is the outward pointing normal vector. We chose to use the local Lax-Friedrichs flux for f * , which averages the flux on both sides of boundary and adds an additional diffusion term smoothing out jumps across the boundary: where and the indices + and -denote the outward (neighboring) and interior element at the boundary. The brackets denote a jump across the boundary, f max is the local maximal wave speed, which corresponds to the speed of the fastest propagating mode across the boundary. In one dimensions this corresponds to The local Lax-Friedrichs flux is the most natural extension from the analytic solution of linear conservation laws to the non-linear case. It relies on the so called Roe condition, which reflects the assumption that the system is dominated by one strong wave.

Non-conservative product
The additional non-conservative flux across a boundary is given by D. The theory of non-conservative fluxes was developed in [80,81] and is applied in the context of Finite Volume and Discontinuous Galerkin schemes [82][83][84][85][86][87][88][89]98]. To compute this quantity we need to consider the general form of a flux across an interface. For this purpose we consider a path φ i (s) along the solution u i , with start and endpoint u L i and u R i respectively and the parameter s ∈ [0, 1]. The formal definition of the flux along this path for a non-conservative flux contribution a i ∂ ρ u j (see (16)) is then given by: We remark that in the non-conservative case the flux is dependent on the chosen path. By choosing the right and left sides of a boundary u R = u + and u L = u − we are able to compute the flux from one cell to another. Similar to the numerical flux, D has to satisfy the jump property for consistency, The last remaining degree of freedom are the boundary conditions for the outer boundary of Ω h . In our case they are given by the in-/out-flowing flux, which is implemented by setting u + i,h = u − i,h at the outer boundaries, effectively adding an imaginary additional cell.
It follows that the non-conservative flux is not fit for flux-boundary conditions, since the nonconservative flux vanishes at the outer boundaries due to the jump property. Therefore the equations need to be reformulated such that the boundary-conditions can be met using the conservative flux. This is done in Appendix E.

Time stepping
The solution is computed by an explicit third-order time-stepping scheme from the dune-pdelab module, where we additionally implemented Courant-Friedrichs-Lewy (CFL-)conditions. The time step ∆t is thus limited by the propagation speed of the flow: where ∆x is the size of the grid cell, N the polynomial degree used in the computation, such that the denominator indicates the total amount of grid points within the respective cell. f max is the maximal propagation speed of the information and was defined in (A4). Additionally we use a minmod slope limiter before each computation step to suppress oscillations around kinks and jumps in the solution.  Table V

Appendix B: Convergence Properties
Results from the large-N computations with constant Yukawa coupling ( Section IV B) of different approximation orders N p are compared to benchmark the accuracy of the computations. Since there is no analytic solution available to compute the numerical error generated by the DG-scheme, we use the result of a computation with N p = 5 and K = 700 elements u ref as a reference. The results are compared at an energy scale of k = 140 MeV. The discrete time-stepping was adjusted such that the time-step computed by the CFL-condition is lowered to ensure k = 140 MeV is reached exactly. The discrete solution is used to generate an interpolating function u, from which the L 2 norm ||u − u ref || L 2 ,Ω h , is computed. This is done at temperature T=10 MeV for µ = 0 MeV and µ = 400 MeV. The results are depicted in Figure 12a and Figure 12b respectively. In vacuum we recover good convergence properties observed in [67] and we perform a fit of the first 10 data-points to the parametrization, The convergence behaves like a power law when increasing the number of elements K and shows   Figure 12b using Equation (B2). exponential convergence when increasing the local approximation order N p up to about K = 10 2.1 ≈ 125.
The fit parameters are given in Table V and the fit to the data points is included in Figure 12a. For higher K we the error decreases even faster which is due to the numerical error of u ref which e.g. also contains some of the resolution issues around the kink in the potential. When compared to [67] the relative error between u and the reference solution u ref significantly bigger. This is in part caused by the fact that the functions were interpolated and not reconstructed using the original basis polynomials, since the Dune output only contains the values at cell boundaries with a precision of eight digits. Adding a finite chemical potential increases the effect of the source term in Equation (D1). At low temperatures the silver-blaze property introduces a sharp onset of chemical potential that shows better convergence for lower and odd N p .
However, since the exact solution is not known it is difficult to judge which local approximation is the best choice as a reference. Figure 12b illustrates the convergence properties in the area in which the source term dominates the equation over the flow.

Convergence in Systems with Diffusion
In this section we will comment on the Convergence of the equations in the Finite N f case.
• We retain convergence properties similar to the previous section in regions of the phase diagram where the sigma mode is not critical. This is supported by the observation that the Finite-N f and the Large N f simulations behave very similarly (see eg. Figure 2, Figure 4a). The Courant number is a factor C chosen to ensure the inequality in (A9). In hydrodynamics C = 0.01 is a common choice for diffusive systems, which is appropriate if the flow is not diffusion-dominated.
• In the diffusion dominated scenario, the solutions are not convergent. Additionally the time stepping becomes effectively 0. This is due to shock development in the solution, which creates a steep negative slope, i.e. a very high (divergent) diffusion flux contribution.
The lack of convergence is explained by the fact that approximate Riemann solvers are only applicable in convection dominated systems. The convergence of diffusion dominated flows can only be ensured using a new formulation of the fluxes, such as the local-DG methods [121].

Convergence with a non Conservative Flux
In this section we will comment on the convergence of the system with non conservative flux, i.e. the computations for field dependent Yukawa coupling.
We can again distinguish two cases: • The area of the phase diagram where the pion mass m 2 π has a single minimum for the entirety of the flow: In this case we retain similar convergence properties as in section Appendix B.
• At high chemical potential ∂ ρ m 2 π begins to take on high negative values. This directly feeds back to the non-conservative flux. In this area the nonconservative flux dominates the flow, which is not contained in the CFL conditions.

Appendix C: Convexity restoration and time-stepping
In this section we are going to inspect the timestepping and its related problems more closely. The equations are solved by a numerical step-wise integration of the RG-time using the CFL-conditions introduced in Appendix A 3. The size of the integration step is dependent on the information flux between cells, the local wave speed f max , which is defined in (A4). The local wave speed is plotted in Figure 13 for approximate vacuum in the broken symmetry phase and for high temperatures and chemical potential in the symmetric phase.

Approaching Convexity in the Broken Symmetry Phase
It can be observed from Figure 13 that the broken symmetry phase has a steadily increasing maximum wave speed. This corresponds to steadily decreasing time steps and leads to long computation times. This behavior is caused by the time-steps inverse proportionality to the flux: ∆t ∝ (k 2 + u) n/2 , for some positive integer value of n. The two-point function Γ (2) can have negative eigenvalues during the RG-flow, which is what happens to the computed function u = m 2 π = Γ (2) ππ , the pion mass, in approximate vacuum. The flow is self-regulating, ensuring that the expression in the square root in the proportionality remains positive: The closer the root gets to becoming negative, the stronger the flow increases u, causing the modulo |u| to always be slightly smaller than k 2 . This must hold for k → 0 from which it follows that u → 0, such that Γ (2) = 0 at infinite RG-time and convexity is restored. |u| teeters on the edge of becoming bigger than k during the entire integration which results in a big flux between grid cells and very small time steps.

Convexity in the Symmetric Phase
The pion mass becomes positive at some point during the interpolation in the symmetric phase and convexity is restored before k = 0. Positive values of u also significantly decrease the flux between grid cells as can be seen from Figure 13 and increase the size of time steps.
At high temperatures the positiveness of u is caused by the fact that the quark contribution to the flow that initially decreased u is much smaller due to the dampening by the Fermi-Dirac Distribution. This translates into k 2 −u never being remotely close to |0| and therefore no increased convexity ensuring flux. Initially the maximum wave speed at high chemical potentials is similar to the approximate vacuum. However, the sudden onset of density at ln Γ µ drives u to positive values and the information flux decreases.
It can be seen that the introduction of the field dependent Yukawa coupling only slightly increases the wave speed in the broken symmetry phase and has no effect in the symmetric phase, which is to be expected from the observation that m 2 q barely changes during the RG-time evolution at high temperatures or densities made in Section IV C.

Problems and Challenges with Time-Stepping
A recurrent struggle while solving the non-linear RGequations using CFL-conditions is time stepping. The RG-flow continuously works to restore convexity, relying heavily on the fact that k 2 > |u| for negative u. For some cases where k 2 > |u| is very small, a numerical error, for example an oscillation around a shock, can cause the radicant to become negative, and therefore no longer well defined. Limiting the time step by using (A4) is spoiled further at high densities, where the sigma-mode becomes critical. The introduction of a steep slope in the potential introduces contributions to the flow that are not accounted for by time stepping.

Appendix D: Flow Equations of Pion and Quark Masses
In this section the equations are reformulated to simplify their numerical treatment. The flow equation of the pion mass is obtained by taking a ρ derivative of the effective potential. In case of the large-N model this is given by, The flow equation of the Yukawa coupling in (15) is rewritten in terms of the quark mass squared m 2 q (ρ). To this aim, we multiply the original flow equation by 4h(ρ)ρ, which gives where corresponds to the contribution of the pion tadpole diagram and (1,1) (m 2 q,k , m 2 π,k ; T, µ) , to the mixed contribution in Figure 1. The explicit form of the threshold functions is given in Appendix H.
Appendix E: Calculation of the Non-Conservative Numerical Flux The flow equation for the Yukawa coupling (15) was reformulated in Appendix D to suit the general form of the partial differential equations given in (16) and contains a non-conservative flux term and a source term s, The exact definition and derivation of the appearing terms is given in Appendix H.
The non-conservative flux is computed using the integral derived by the jump condition in Appendix A in (A8). We chose a straight path across an interface, w(s) = w − + s(w + − w − ) , u(s) = u − + s(u + − u − ) .
We note again that this is a path along the solutions u and w and not a path in the 'spatial' coordinate ρ. The straight path was chosen because it is often the simplest choice for the evaluation of the integral. In our case the expression simplifies so much that it can be evaluated analytically, due to the explicit form of the equations, where the non-conservative flux is given by A(u(s)) = ∂ u g (u(s)) = 1 u + − u − ∂ s g (u(s)) .
where we used in the last equality that ∂ s w is a constant expression. Instead, the constant C is simply the absolute value of the jacobian matrix. Often, it can be approximated as the maximal characteristic speed of the non conservative product. Note that for constant u across the interface g(u + )−g(u − ) u + −u − = A(u), such that we recover a conservative flux for constant u. There is a large set of paths across the interface that lead to the same value in the integral, due to the fact that A can be written as a derivative of u. This hints at the possibility that there might be a conservative formulation for the system of equations.
Since this formulation only allows flux boundary conditions for conservative fluxes a partial integration is performed, We now have a conservative flux A(u k )w k the proper in-/out-flow boundary conditions for w k and a very small non-conservative flux contribution D accounting for jumps in u k , This contribution is very small when u k is smooth and only contains small jumps across interfaces. It obviously vanishes at the outer boundary since there we have u + k = u − k . Thanks to this formulation of the equation the maximal wave speed of the non conservative product is rather small and can be neglected in practice. In the general case however the inclusion of this term is important especially if the non conservative product are the only convective term in the equation, since it introduce the necessary numerical dissipation to make the numerical scheme stable.