Influence of broken flavor and C and P symmetry on the quark propagator

Embedding QCD into the standard model breaks various symmetries of QCD explicitly, especially C and P . While these effects are usually perturbatively small, they can be amplified in extreme environments like merging neutron stars or by the interplay with new physics. To correctly treat these cases requires fully backcoupled calculations. To pave the way for later investigations of hadronic physics, we study the QCD quark propagator coupled to an explicit breaking. This substantially increases the tensor structure even for this simplest correlation function. To cope with the symmetry structure, and covering all possible quark masses, from the top quark mass to the chiral limit, we employ Dyson-Schwinger equations. While at weak breaking the qualitative effects have similar trends as in perturbation theory, even moderately strong breakings lead to qualitatively different effects, non-linearly amplified by the strong interactions.


Introduction
With the detection of gravitational waves [1] a whole new era of astronomy has begun. Eventually, this will allow to investigate neutron star mergers. In such dense environments, weak interactions become so prevalent that the dynamical backcoupling between the weak and the strong interaction becomes relevant, see, e.g., [2][3][4][5][6][7][8][9]. This requires therefore a fully coupled, and necessarily nonperturbative, description. This is so far only possible at the level of comparatively simple effective models, but not yet in an ab initio calculation. This is not the only reason to consider this problem. Beyond the known standard model (hidden) sectors may exist in which strongly interacting parity-conserving and parity-breaking interactions both exist.
Both these considerations motivate to understand such backcouplings better. The main hallmark of both scenarios is the appearance of explicit C and P symmetry breaking, as well as, to a lesser extent, flavor breaking. The presence of additional, or non-negligible, symmetry breaking effects implies always a more involved tensor structure of correlation functions. Therefore, we focus here on the simplest object, which exhibits the full additional complexity, the quark propagator (QP). Since we are mainly interested in how the strong interaction may amplify, or modify, the symmetry breaking effects, we will consider a e-mail: axel.maas@uni-graz.at b e-mail: walid.mian@uni-graz.at only an explicit source of the symmetry breaking, as will be discussed in greater detail in sect. 2. For the weak interactions, which, due to explicit masses of the W and Z bosons do not show a strong momentum dependence at low energies, this is a sufficient approximation.
However, the inclusion of such a breaking, and the large differences in relevant energies, already limit the possible choices of methods. Especially, lattice gauge theory is not suitable. This is, on the one hand, due to the immense computational costs for the vastly different energy levels involved. On the other hand, there is no fully proven way yet to upgrade the static breaking considered here to the full weak interactions using lattice methods [10].
As an alternative, here functional methods in the form of Dyson-Schwinger equations (DSEs) will be employed. These have been used very successfully to determine the quark propagator in QCD in various levels of sophistication [11][12][13][14][15]. In the more exploratory investigations here, the most important features are the dynamical mass generation as well as the correct implementation of chiral symmetry. These features are particularly well implemented in the so-called rainbow-ladder truncation [11,12,14]. This completes our setup, which we describe in much more detail in sect. 2 as well as in appendices A and B. Especially in appendix A, we will discuss the tensor structure of the quark propagator, which has now four matrix-valued, rather than two real-valued, dressing functions, demonstrating the much higher complexity compared to QCD alone.
The most important results will be discussed in sect. 3. Probably the most relevant insight gained, aside from the necessary technology to deal with the increase of complexity, is that the fully coupled system behaves as expected from perturbation theory only at very weak breaking. Already at moderately small breaking, the amplification by the strong interaction can lead to qualitatively different behaviors for various dressing functions of the quark propagator than perturbatively expected. Of course, only future investigations of observable quantities will tell what the implications for physics are, but the present results mandates caution with extrapolation of perturbative notions. In this section, we will also present information on the analytic structure of the quark propagator using its Schwinger function, an issue which is already in QCD alone highly non-trivial [11,12,14].
Finally, we list many further results in appendix C, providing a complete picture of the quark propagator in this setup for a wide range of parameters and quark masses. All of the insights and results are finally summarized in sect. 4.
Some preliminary results have already been reported in [16].

Ansatz
The breaking of the symmetries will be realized by including an explicitly symmetry-violating term in the Lagrangian. The breaking is thus generated already at tree level. Since the weak interactions motivate our study, the term will break C, P , and flavor within generations. Thus, our QP becomes matrix-valued in flavor space, with the off-diagonal elements mediating flavor changes.
Flavor violation within a generation is actually not possible without involving further particles, due to electric charge conservation. In the standard model, this is ensured by the emission of a lepton and a (anti)neutrino. To avoid this additional complexity, here the leptons are modeled as an external background field. Given that our ultimate interest is in neutron star mergers, where such a reservoir is readily available, this appears like a reasonable approximation.
In the following the subscript u and d denotes up-like quarks and down-like quarks and the superscript L and R left-handed and right-handed quarks, respectively. The strength of the weak interaction, or more precisely the coupling to the symmetry breaking external field, will be denoted by g w , the effective weak strength. We will vary the value of g w from small values to large values to turn on the effects smoothly. Finally, the quark fields are denoted by ψ.
All together leads to the Lagrangian where A i are the gluon fields and T i are the generators of SU (3). m u and m d are the masses for up-like and down-like quarks. L Rest is the remainder of the QCD Lagrangian, which includes the gluon self-interaction and the gluon-ghost part, and does not play an explicit role in the following. It also contains the gauge-fixing terms of our choice of (minimal) Landau gauge. For brevity, we also suppressed the renormalization constants.

Quark propagator
In the following, P AB is the propagator from flavor A to B with A, B ∈ {u, d} for up-like and down-like quarks. Its inverse will be denoted by P −1 AB . Because of parity violation the QPs have in addition to the usual vector and scalar channels also non-vanishing axial and pseudo-scalar channels. The standard notation for the vector-and scalarchannel dressing functions of the inverse Propagator in the literature is A and B. We will keep this and denote the dressing functions for the axial and pseudo-scalar channels of the inverse propagator with C and D. The corresponding dressing functions of the propagator are denoted by a tilde.
This leads to the form The dressing functions of the propagator and its inverse are related with each other. In general the dressing functions of the propagator depends on all the dressing functions of the inverse propagator in a complicated way, see for details appendix A, and generically denoted as Instead of splitting the Lorenz channels into the vector and axial channels, it can also be split into left-handed L(p 2 ) and right-handedR(p 2 ) components, leading to with the relations At tree level the propagator reads where the common denominator is given by .
The quantity σ ab is 1 for a = u and b = d and −1 in the other case. The tree-level propagator already reveals the major contributions for the QPs. They are separated in the different channels at tree level, but will mix in the full case. Consider the denominator. In the second line of eq. (7) we have factorized the denominator to see both poles of the tree-level propagator. For g w → 0 M l goes to the mass of the lighter quark and M h to the mass of the heavier quark. By increasing g w , the value of M l is decreased and M h is increased and thus increases the effect of the mass splitting.
Note that M l and M h diverge at g w = 0.5. At this point the poles turn imaginary, indicating a breakdown of the trivial vacuum around which the perturbative expansion is performed. This feature will actually not be lifted in the full non-perturbative treatment, and we are only able to find solutions as long as g w 0.4. Since this is a very large breaking, probably far too strong for the setting of neutron star mergers guiding this work, we did not endeavor to find out what happens beyond this point, and restrict ourselves to breaking strengths below this value.
The explicit chiral symmetry breaking by the tree-level masses manifests itself in the scalar channel of the pure and mixed flavors. In the chiral limit the scalar channel vanishes at tree level, but due to dynamical chiral breaking the scalar channel does not vanish for the full propagator.
The vector channel of the mixed propagator is proportional to g w (m u m d − p 2 ), which changes its sign for p 2 > m u m d . The influence of this is a contribution in different direction for large and low momenta. For heavier bare quark masses this contribution is shifted to higher momenta.
The most remarkable contribution, a difference of both quark masses, appears in the pseudo-scalar channel of the mixed propagator. This will have a significant impact on the full propagator. It is remarkable that this contribution appears with opposite sign for the propagator from up-like quarks to down-like quarks and the other way around: Although we have taken the same strength for the propagation of both mixed QP in our ansatz, we get a difference, if the quarks have different masses.

DSEs
We can write the QP in a matrix form, where the diagonal elements are the QPs for pure flavor and the off-diagonal elements are the QPs for mixed flavors, see for details appendix A. A graphical representation of the DSEs is given in fig. 1. The equations look similar to those of QCD, where the QP can also be given in a matrix form, but with vanishing off-diagonal elements. It is only the appearance of the off-diagonal tree-level elements, which gives rise to all differences.
The full quark-gluon-vertex appears in the self-energy graph, which is determined by a separate DSEs involving even higher-order correlation functions.
To avoid this complication, the DSEs are truncated at this level, which is known as the rainbow truncation [11,12,14]. The required inputs of the full gluon propagator and the full quark-gluon-vertex are then replaced by a bare gluon propagator and a quark-gluon vertex given by the tree-level tensor structure, but dressed with an effective running coupling α. Denoting the tree-level inverse propagator with P −1 0,AB , the final DSE reads where Z 2,A and Z 2,B are the quark wave functionrenormalization constants for flavor A and B. Λ represents a translationally invariant regularization with the UV-cutoff Λ. μ is the renormalization point and k = q − p.
For the dressing function α we choose the Maris-Tandy coupling [17] α(q 2 ) = π ω 6 Dq 4 e − q 2 Here the parameters are adapted to describe pions in the vacuum adequately. In the literature these parameters are fitted for degenerate masses of up and down quarks [18]. For a detailed analyses of different parameter sets see, e.g., [19,20].  [21]. In addition, we will also consider the chiral limit as well degenerate cases. Let us note here that the PDG values are in a different renormalization scheme, actually even different schemes for different flavors, and at different renormalization points for the different flavors. Here, our aim is to compare different flavors directly. We therefore chose as a uniform renormalization scheme that our quark masses should equal the values listed above at μ = 10 6 GeV, and use a miniMOMtype scheme. At the perturbative resumed one-loop level in our miniMOM scheme, this amounts to a difference of about 25-35% larger masses, depending on flavor, compared to the PDG renormalization scheme. Note that the small difference is due to the appearance of the running coupling at the renormalization point in this scheme. At the qualitative level of our study, this will make little difference. Ultimately, if we would be able to solve the DSEs exactly, the scheme would anyhow not matter for physical observables.
Further details can be found in appendix B.

Schwinger function and masses
To obtain information on the analytic structure, we also determine the Schwinger function [22,23]. It is defined as σ AB (p 2 4 ) is one of the dressing functions from the propagator evaluated at zero spatial momenta ( p = 0).
The actual analytic form of the propagator is yet unknown, but poles and/or cuts appear likely [11,14,22]. If there were only an ordinary mass pole, the Schwinger function would show an exponential decay [23] Δ(t) ∼ e −mt (11) and m would be the mass. However, investigation of pure QCD in the rainbow truncation yielded rather a structure with complex conjugated poles [11,14,22], which is, e.g., expected in the Gribov-Stingl scenario [24][25][26]. Note, however, that this may be a truncation artifact. In this case the Schwinger function is roughly given by The decay rate is given by the real part and the oscillation frequency by the imaginary part of the mass pole.

Results
For each quark generation, we have 2 dressing functions for pure flavor and 2 dressing functions for mixed flavor and each has 4 Lorentz channels, resulting in 16 dressing functions. We will consider 6 cases in total, the chiral limit and three physical quark generations and two cases of degenerate masses. Therefore we have numerical results for 192 dressing functions and each of them as a function of the weak strength. In addition to that, we also have the Schwinger function for each dressing function. To avoid cluttering up the main text, most of these results are relegated to appendix C. Here, only the qualitatively most remarkable results will be analyzed. The results in the appendix do not add any new concepts to this section. The result for the QP in QCD are usually [11,14,22] given in terms of the wave function renormalization Z(p 2 ) = 1/A(p 2 ) and mass function M (p 2 ) = B(p 2 )/A(p 2 ). For ease of comparison, the case with g w = 0 will serve as reference. Therefore in sect. 3.2 the results for Z andÃ AA will be explored. Afterwards, the mass function and the relatedB AA will be discussed in sect. 3.3. In sect. 3.4 we will analyze the results for the axial channel C AA , and study the impact of parity violation. The results for the Schwinger function will be discussed in sect. 3.5. But first, it is necessary to discuss the involved scales, as the problem is now a multi-scale one. Note that in the chiral limit there is no difference between the up-like quark and down-like quark, and thus the flavor-diagonal elements coincide. In these cases, always the up-type one will be shown.

Relative scales
g w is dimensionless. Thus, a possibility for the comparison of the breaking strength with the strong interaction is  up-down charm-strange top-bottom given by g w and the average strength of the Maris-Tandy coupling α in the IR region. Neglecting the perturbative part of the Maris-Tandy coupling, we get where x = q 2 /ω 2 and the integration region for the average is taken from q = 0 to q = 1 GeV. The weak coupling is constant throughout. The smallest values of the weak strength g w = 10 −6 is many order of magnitude smaller than α and thus represent the natural case. The largest value of g w = 0.4, which is only one order of magnitude smaller than α, substantially deviates from nature. This also justifies our choice to restrict to not too large values of g w . However, such effects may play a role in theories with strongly interacting chiral sectors. The interaction also manifests itself into a momentum scale in particular appearing in the scaler channel. The weak strength is transmuted into a momentum scale in the flavor off-diagonal element, which is zero without breaking. It is shown for various quark masses in fig. 2.
It is seen that in the IR this dressing function is approximately proportional to the weak strength for g w ≤ 0.1. Note that the generated scale depends on the quark masses. For both light generations the difference to the chiral limit is small. For top and bottom quark, the generated scale is one order of magnitude bigger at the same g w . Thus, there is a linking of the different involved scales.

Wave function renormalization
From eq. (B.3) follows that A uu andÃ uu are directly linked with each other, and thus Z = 1/A uu is also directly related toÃ uu . These dressing functions are shown in fig. 3 for different values of g w in the chiral limit. For values of g w 0.01 no appreciable effect is seen. At larger valuesÃ uu slightly increases in the UV when increasing g w . In the mid momenta regime it is slightly decreased and in the IR it is significantly increased. A consequence of this is that Z also increases in the IR. This can be understood from (B.3), as A uu is obtained from integrating A uu multiplied with a kernel over all momenta. Because   A uu is increased in the UV very little and more decreased in the mid range, Z is slightly decreased in the UV range due to the integration. For the same reason Z is not increased as much asÃ uu is increased in the IR. For other quark masses the same behavior is seen, as shown in fig. 4. The graph shows thatÃ is increased in the IR for all quark flavors and thus Z is also increased.
The effect comes from different sources. One is from g w and the other from a combination of g w and the bare quark masses. Especially, A and Z are increased for the up quarks more than in the chiral limit. Also, in general the value for up-like quarks is increased more than for downlike quarks. This is also seen in fig. 15 in appendix C.1 in more detail.
This can be understood from eq. (6) for the tree-level case. One of the contributions arises from the bare quark masses and another from mass splitting with different signs. This creates the cross-talks leading to the observed effects. We will return to this later in sect. 3.4. In addition, the absolute value is decreased for higher bare quark masses, as anticipated because the masses of the quarks enter in the denominator of the QP. This can be seen already for the tree-level propagator in eq. (7).

Mass function
The relations between the dressing functions of the QP and its inverse are more involved as in QCD, see also appendix A. In QCD, the relation is given bỹ To be able to compare, we therefore choose to define a (pseudo) mass function as which by construction coincides with the usual one in the QCD case. Of course, neither in QCD nor here this function needs to coincide with the actual mass. Any such statement requires the Schwinger function in sect. 3.5. Nonetheless, we will stick here with the usual convention and call this quantity mass function. The dependence on g w of this mass function is shown in fig. 5. As in QCD, the mass function is non-zero, indicative of chiral symmetry breaking. The mass function starts to change appreciably for g w 0.01, like the wave function renormalization. The same is true forB uu , which is also shown in fig. 5. Since the connection between B uu and B uu , due to eq. (B.3), is similar as for A uu andÃ uu , the same analysis as before applies, and the response to g w follows the same pattern.
At non-zero masses, the picture changes. This is shown in fig. 6 for up and down quarks. The mass function M of the up quark is decreased by increasing g w for g w 0.3. For larger values it increases again, but here our approximations start to break down. This replicates the result of the tree-level propagator in sect. 2.2: The mass of the heavier quark in a generation is increased and the mass of the lighter quark is decreased.
The other flavors are shown in fig. 7. The second generation follows the pattern of the first, but not the third. In the latter case the mass function of both quarks increases. This implies different contributions to the mass functions. One increases the mass function of the heavier quark and decreases the one of the lighter quark. The other contribution increases with the mass of both quarks. An indication of this is already seen in the second generation, albeit not creating a qualitative change.

Parity violation
In the following the handiness of the quarks is investigated, using the definition (5). This requires the axial channel, shown in fig. 8 for the chiral limit. The corresponding dressing functionC is for the flavor-diagonal elements found to be positive for higher momenta and negative in the IR. At the same time, for increasing g w the absolute value ofC also increases. Of course, at large momenta the dressing function goes to its tree-level part, which from eq. (6) isC which is positive, actually for all momenta. Therefore the backcoupling to QCD forces it to be negative in the IR. This happens at a transition scale of approximately 1 GeV 2 , which is the typical QCD scale.
To assess the consequences of this for the left-handed and right-handed contributions, it is useful to define their relative ratio as SinceÃ is always positive, the sign ofr is given by the sign ofC. This already entails a change of sign, and that the left-handed part is larger in the infrared. This is also shown in fig. 8. The effect increases non-linearly with g w : For g w 0.1 the absolute value |r| in the UV and IR is increased by two order of magnitudes, when g w is increased by one order of magnitude. The flavor-off-diagonalC is always negative in the chiral limit, butÃ changes its sign, see figs. 17 and 22 in appendix C. At the same transition scale of approximately 1 GeV 2 as for the flavor-diagonal case, and in the chiral limit,Ã has a zero crossing andC not. This leads to a divergingr ud at this scale. Thus also in this case there is a transition from right-handed in the UV to left-handed in the IR. This is shown in fig. 9. Increasing the mass, the situation for up and down quarks is shown in fig. 10. The absolute value ofC is different for up quark and down quark, but for g w = 10 −5 the behavior for up and down quark is as in the chiral case. Slightly increasing g w to 5 · 10 −5 entails a drastic qualitative change. For the up quarkC is still positive in the UV and negative in the IR, but for the down quark it remains positive for all momenta. Therefore the up quark still flips its chirality at long distances, but the down quark does not do so.
To understand the origin of this effect, it is helpful to study the degenerate mass case, also shown in fig. 10. There is no (numerically detectable) difference between up and down quark 1 , andC changes again sign, as in the chiral limit. For higher values of g w the absolute value of C just increases in the IR. This implies that the different behavior ofC for the non-degenerate case is due to the mass splitting of the quarks. Since at tree level only in the pseudo-scalar channel a contribution proportional to the mass splitting occurs, this effect must have been propagated by the QCD interaction to the axial channel. Moreover, the effect becomes already important for a very small mass splitting and a very small breaking scale of g w ≈ 5 · 10 −5 , compared to QCD. This can only happen if there is a strong non-linear amplification mechanism at work. This implies that the QCD medium strongly affects the helicity at long ranges, but only for non-degenerate quark masses. That was certainly not expected.
The same is true for the other quark generations, see for the second generation fig. 9. The effect is still there for the third generation, with its very large mass splitting, see fig. 11, and there occurs already at an even smaller breaking strength of g w = 10 −6 . Thus the absolute value of the involved mass scales amplifies the non-linear backcoupling, such that it occurs at weaker breaking strength.
The corresponding relative ratios are shown in fig. 25 in appendix C.5. These graphs support the existence of a transition scale, where the left-handed and right-handed contribution change their relative contribution.
For the mixed flavor case the effect is solely driven by the absolute value of the mass splitting, and the breaking strength plays only a minor role. This is shown in fig. 12. Whether there is a mass splitting or not plays    only a role if the mass splitting is large enough, i.e. in the third generation. Only then the behavior with or without mass splitting differ qualitatively in the infrared. In fact, already the second generation is sufficient for this, as can be seen in fig. 9.

Schwinger function
As noted, the analytic structure is accessible through the Schwinger function (10). Let us now consider the Schwinger function in the chiral limit for the flavordiagonal dressing function B. The other flavor-diagonal dressing functions do not lead to qualitatively new results, and are even quantitatively similar, so these will be skipped here.
The results are shown in fig. 13. The Schwinger function shows an oscillatory behavior, consistent with the form (12), and thus complex conjugate poles, as in pure QCD for the rainbow-ladder truncation [22]. In fact, a fit using a more detailed ansatz, see appendix D, of this type works very well, as is also shown in fig. 13. The values of the fit parameters are listed, for completeness, in appendix D.
However, the oscillation period starts to substantially increase for g w > 0.01, up to a point where at large g w the first zero crossing has moved to a time which we can no longer numerically resolve reliably. Thus, the imaginary part shrinks with increasing breaking. The curvature at short distances is still not quite right for a physical particle. A similar behavior, though with a suppression of oscillations for decreasing interaction strength, has already been observed for adjoint scalar particles [27,28]. This strongly suggests that the interaction strength plays a crucial role for the scale at which negative norm contributions become relevant, even though the coupling does not differentiate between positive-norm and negative-norm states.
At the same time the steepness decreases, making the real part smaller. Thus, the increase in g w moves in total the poles closer to the origin, as both real and imaginary part decrease.
The flavor-off-diagonal Schwinger function, again only the scalar part as the others are very similar, is shown in fig. 14. In principle, it shows a very similar behavior as for the flavor-diagonal part, except that it always retains a first zero crossing, relatively independent of g w , at very short times. It is still possible to fit it using the same fit form. The results for the fit are also listed in appendix D. It is found that the zero crossing at small momenta comes from the phase shift δ in (12). It is very close to π/2 and causes the sign change for the Schwinger function at small t. Still, the position of the pole also moves towards the origin with increasing breaking strength.
The Schwinger function for the first two generations show the same behavior as in the chiral limit, see fig. 26 in appendix C.6. For the third generation the fall-off was too fast, due to the large real part, as that any unambiguous statements could be drawn before numerical noise drowns out the signal.
It is a quite interesting result that the breaking pushes the poles closer to the origin. As we expect a change of physics when crossing the threshold g w 0.4, this could be a first indication of a drastic change at strong breaking. However, this is probably not of relevance to neutron star physics. On the downside, the decreasing distance to the origin will create additional problems in any mesonic correlators in rainbow-ladder calculations [11,12,14]. In these cases more elaborate schemes will be necessary than a tree-level breaking, which we are currently developing.

Conclusion
We have calculated the quark propagator in the presence of explicit flavor, C and P symmetry breaking. Moreover, we took into account the non-linear back-coupling from QCD in the rainbow-ladder truncation. The latter leads to qualitative effects, even for relatively small explicit breaking strengths, at long (hadronic) distances. They also couple in a highly non-trivial way the various dressing functions to each other. This was particularly visible in the way how effects from mass splitting and mass averages surfaced in various dressing functions. The non-linear amplification also surfaced in other ways. This is a very important insight: External perturbations, even in rainbowladder truncation, can be substantially amplified by the strong interactions. This must be regarded as a warning that even small effects can play a non-perturbatively large role when QCD is involved.
From the point of view of physics, another interesting insight is obtained when considering how left-handed and right-handed particle propagations change. Under particular conditions, flips between handedness can be amplified at long distances by the strong interaction. This can deplete or enlarge the available particles in some handedness. As the weak interactions only couple to a particular handedness, this can increases or decrease the reservoir of particles which are weakly interacting in a system. If this pertains to the full system, this can influence the dynamics in forming or merging neutron stars, as this could alter, e.g., the opacity for neutrinos. This is even more important as the typical range where this occurs is only of the size of a hadron.
Concluding, this investigations showed that weak interactions effects, even if themselves small, can be amplified by the strong interactions, and this backcoupling can have qualitative impact. Keeping this in mind will be important in the next step, when relaxing the assumption of a reservoir, and taking the weak interactions explicitly into account, including the emitted neutrinos and electrons. This will require to work on a hadronic level, which is our next aim.
Open access funding provided by Karl-Franzens-University Graz. We are grateful to Helios Sanchis-Alepuz, Jordi Paris-Lopez and Adrian Lorenz Blum for helpful discussions. WM has been supported by the FWF doctoral school W1203-N16.

Appendix A. Structure of the quark propagator
In the following a number of useful relations between the dressing functions of the quark propagator and its inverse will be collected. Combining the flavor elements in a matrix, e.g. for the vectorial channel as allows for a compact notation. Writing the propagator and its inverse in terms of these matrices yields P (p 2 ) =Ã(p 2 )i / p +B(p 2 )1 1 +C(p 2 )i / pγ 5 +D(p 2 )γ 5 , The propagator satisfies the condition yielding relations between the matrix-valued dressing functions While this system is a system of linear equations for either the matrix elements of the propagator or its inverse, an explicit solution is of little use. The expressions become extremely lengthy, and therefore prohibitively expensive to evaluate during numerical calculations. Therefore, in our investigations we always solved such equations numerically at double precision.

Appendix B. DSEs for quark propagators
To derive the DSEs for the different dressing functions, we insert in eq. (8) the quark propagator of eq. (2) and project out the different channels, by taking suitable traces. We define the following two kernels: We further define the functional Π i for i = 1, 2 as This yields the DSEs of the different dressing functions for flavor-diagonal elements A ∈ {u, d} and for mixed flavors, A, B ∈ {u, d}, A = B, This system of equation is technically very similar to the ordinary rainbow-ladder truncation. Therefore a numerical solution using standard fixed-point iteration schemes is possible, and was done here. Only the quark propagator was numerically inverted at every step, due to the involved structure, as discussed in appendix A. There are, however, a few more subtle numerical issues to be mentioned. To perform the integral for Π i,A the dressing functionsÃ,B,C andD are evaluated for various momenta using interpolation. We performed our calculations using linear and cubic interpolation. If we use a precision of 5 × 10 −5 , which is considered to be sufficient in the standard fixed-point iteration scheme, for A, B, C and D then we get different numerical solutions for linear and cubic interpolation. By increasing the precision the solution from the linear interpolation approaches the solution from the cubic interpolation for all dressing functions except for D AA andD AA . Especially, using linear interpolation we get different solutions for D AA andD AA by using different start values for the iteration. In contrast to this, we get the same solutions for different start values using the cubic interpolation. Thus, we consider the solutions from the cubic interpolation as stable, and used them throughout this work.
In addition, we used a precision of 5 × 10 −7 , instead of the standard value, for A, B, C and D and 2 8 = 256 grid points. In this case the difference for the solutions from the linear and cubic interpolation were at most in the third significant digit, and thus lead essentially to the same results.
Let us note that our precision for the iteration procedure is for the dressing functions of the inverse propagator (A, B, C and D). The dressing functions of the propagator (Ã,B,C andD) are calculated by a numerical inversion, with a precision of roughly 10 −20 .

Appendix C.1. Vector channel
The flavor-diagonal vector dressing functions for the different generations and values of the breaking strength are shown in fig. 15. A detailed discussion is given in sect. 3.2. A comparison between the different generations, the chiral limit, and tree level for the flavor-off-diagonal vector dressing function is shown in fig. 16 for a fixed value g w = 0.2. The tree-level value is given by, see eq. (6), This demonstrates that at tree level the dressing function is negative for momenta p 2 ≥ m u m d and positive for lower momenta. The full dressing function exhibits the same behavior, but the absolute value is substantially larger. Note especially that in the chiral limit the tree-level propagator is proportional to −p 2 and thus negative for all momenta. The full dressing function is, however, positive in the IR. As the masses of up and down quark are comparatively small, the result for them is close to the one in the chiral limit.
The dependence on the generation and g w is shown in fig. 17. In every generationÃ ud increases with g w . The zero ofÃ ud shifts for higher momenta for heavier mass. This is already the case for the tree-level propagator, where the zero is determined by the condition p 2 = m u m d .

Appendix C.2. Scalar channel
Complementing the results in sect. 3.3 the flavor-diagonal scalar dressing functions for the different quarks are shown in fig. 18.
A comparison to the tree-level case of the third generation flavor-off-diagonal scalar dressing function is shown in fig. 19. Note that the tree-level result is, see eq. (6),    Thus, at tree levelB 0,ud is negative for all momenta, in particular in the UV. The latter is also seen in the full case. But it switches to a positive value in the IR. This qualitative behavior is also seen for the other generations, as is also plotted in fig. 19, but at differing absolute values. This persists even in the chiral limit. The value of this quantity is found to also increase when increasing g w , which is shown in fig. 20.