Vector and Axial-Vector Mesons in Nuclear Matter

As a first step towards a realistic phenomenological description of vector and axial-vector mesons in nuclear matter, we calculate the spectral functions of the $\rho$ and the $a_1$ meson in a chiral baryon-meson model as a low-energy effective realization of QCD, taking into account the effects of fluctuations from scalar mesons, nucleons, and vector mesons within the Functional Renormalization Group (FRG) approach. The phase diagram of the effective hadronic theory exhibits a nuclear liquid-gas phase transition as well as a chiral phase transition at a higher baryon-chemical potential. The in-medium $\rho$ and $a_1$ spectral functions are calculated by using the previously introduced analytically-continued FRG (aFRG) method. Our results show strong modifications of the spectral functions in particular near the critical endpoints of both phase transitions which may well be of relevance for electromagnetic rates in heavy-ion collisions or neutrino emissivities in neutron-star merger events.


I. INTRODUCTION
The properties of matter under extreme conditions in temperature and/or density, as encountered for instance in the early universe, the core of neutron stars, and binary neutron star mergers are in the focus of ongoing theoretical as well as experimental and observational efforts. Hot and dense strong-interaction matter is created and studied in relativistic heavy-ion collisions at the world's most powerful accelerator facilities while the properties of neutron stars and their dynamics are inferred from observations of electroweak signals and more recently from gravitational radiation in merger events. In both cases, key challenges include the investigation of the in-medium modifications of hadrons and their connection to changes in the underlying symmetries, the equation of state, transport properties, and the phase structure of strong-interaction matter. Here, the basic features of Quantum Chromodynamics (QCD), dynamical chiral symmetry breaking and confinement, play a decisive role, as it is expected that at high temperatures and baryon densities chiral symmetry gets restored, and that quarks and gluons are liberated as confinement disappears. For overviews on heavy-ion measurements, astrophysical observations and theoretical studies, see e.g. the reviews [1][2][3][4][5][6][7][8][9][10][11][12] and references therein.
To determine the electroweak response of compressed and hot nuclear matter, a realistic theoretical description of the in-medium ρ and a 1 spectral functions is required. Both are parity as well as chiral partners of the global chiral SU (2) L × SU (2) R symmetry of QCD for N f = 2 light quark flavors. Its realization therefore plays an important role for their in-medium properties. However, a consistent description that incorporates this (approximate) chiral symmetry and its spontaneous breaking as well as the effects from critical fluctuations, most notably from a critical endpoint of a possible phase transition inside a dense nuclear environment, is still missing.
In this work we therefore present a new setup for the calculation of vector and axial-vector meson spectral functions in dense nuclear matter. Our approach is based on the Functional Renormalization Group (FRG) which represents a non-perturbative framework that is capable of including both quantum and thermal fluctuations [13][14][15][16][17][18][19][20]. Most notably this includes order-parameter fluctuations due to the dynamics of collective excitations which are not accessible in mean-field approximations. This makes the FRG particularly well suited to study the critical behavior of the corresponding correlations. Moreover, as other functional methods, the FRG is not hampered by the fermion sign problem encountered in lattice QCD [21] at finite baryon density. While the region of finite temperature at low net-baryon density, for baryon chemical potentials µ B less than about twice the temperature T where lattice QCD results are available, is important to benchmark effective theories and approximations, including the unavoidable truncations in functional methods, the FRG can therefore then be applied with some confidence also in the dense region of the phase diagram of strong-interaction matter, with µ B of the order of the nucleon mass at temperatures that are at least an order of magnitude lower than that. In addition, the necessary truncations can be made to preserve the global symmetry structure and its breaking patterns as here described by the underlying effective theory.
To construct an effective low-energy description for nuclear matter that is consistent with chiral symmetry and its breaking pattern, the notion of parity-partners in the bosonic sector has to be extended to massive fermions. This is accomplished in the parity-doublet model (PDM), or mirror baryon model [22][23][24]. The PDM describes nucleons along with their parity partners and can account arXiv:2105.00861v2 [hep-ph] 20 Aug 2021 for a finite nucleon mass in a chirally-invariant fashion. It is motivated by the assumption that a large fraction of mass of the nucleon in QCD is generated by the gluonic contribution to the scale anomaly [25], and not through dynamical chiral symmetry breaking. The PDM also provides a natural description for the parity-doubling structure of the low-lying baryons observed in recent lattice-QCD calculations [26,27]. The mean-field phase diagram [28][29][30][31][32][33][34][35][36][37][38][39] is known to consist of two distinct first-order phase transitions: the usual nuclear liquid-gas transition together with a second, chiral transition at higher chemical potentials. In particular the existence of this chiral transition inside the dense nuclear-matter phase was confirmed including fluctuations within the FRG in [40] to be a robust prediction of the PDM. Just like ρ and a 1 , the nucleon N (938) and its parity partner, commonly assigned to the 1/2 − N * (1535) resonance become (almost) degenerate at this transition, with a common and chirally invariant finite baryon mass m 0,N from the scale anomaly. This model thus serves as suitable effective theory to describe a chiral phase transition inside nuclear matter entirely in terms of hadronic degrees of freedom.
Experimentally testable predictions from this scenario could range from an enhanced production of η mesons in heavy-ion collisions at low beam energies, when a population imbalance between N and N * is created in such a transition, to an enhanced dilepton signal from critical fluctuations when the trajectory in the phase diagram of the expanding system comes close to the associated chiral critical endpoint. We will focus on the latter here, and calculate as a first step towards obtaining electromagnetic spectral function and thermal dilepton rates [41], the in-medium spectral function of the ρ vector meson and its chiral partner, the a 1 axial-vector meson in nuclear matter within this PDM setting.
For the calculation of the spectral functions we use the analytically-continued FRG (aFRG) framework developed in [42][43][44]. The aFRG method avoids the need for numerical reconstruction schemes, see for example [45][46][47][48][49][50][51][52][53], and it is thermodynamically consistent in that the thermodynamic grand potential and the spectral functions are calculated on the same footing. The aFRG method has been successfully applied in different situations, for example to calculate in-medium spectral functions of pions and the scalar σ meson [43,44,54], the quark spectral function [55,56] as well as vector-and axial-vector meson spectral functions at finite temperature and density in extended linear-sigma models with quarks [57,58], together with the corresponding electromagnetic spectral function and thermal dilepton rates [41] inside quark matter.
In this paper we use the PDM as our effective theory describing the two iso-doublets of N and N * as parity partners with chiral representations in mirror assignment and chirally invariant mass term. These interact via chirally invariant Yukawa couplings with pions and the scalar σ meson as well as ρ and a 1 mesons as the re-spective chiral partners in the (pseudo-)scalar and the (axial-)vector meson channels. The formalism to describe the fluctuations of massive vector and axial-vector mesons in an effective theory based on (anti-)selfdual field strengths is taken from [58]. Compared to the precursor study in [57], the inclusion of fluctuating vector and axial-vector mesons in the aFRG flows allows, in particular, to account for important additional contributions to their spectral functions such as, e.g., the three-body resonance decay a 1 → ρπ → 3π.
Our work represents a first step towards a realistic description of vector and axial-vector mesons in nuclear matter. While the former are particularly relevant for the interpretation of dilepton spectra in heavy-ion collisions, both are needed to determine the neutrino emissivities in proto-neutron stars and binary neutron-star mergers. The in-medium spectral functions provide access to realtime quantities such as pole masses and decay widths but also to other observables such as pertinent transport coefficients. Moreover, since our FRG treatment is thermodynamically consistent and symmetry preserving, the in-medium modifications of the spectral functions can be stringently connected to the restoration of chiral symmetry.
The remainder of this paper is organized as follows. In Sec. II the theoretical setup as well as results for the phase diagram and the Euclidean (screening) masses are presented. In Sec. III we discuss results for the real-time two-point functions and the in-medium spectral functions of the ρ and the a 1 meson. We close with a summary and outlook in Sec. IV. Further details are deferred to an appendix.

A. The Effective Average Action
The question of how to describe baryons within effective models incorporating the principle of chiral symmetry has a long history. Here, important work was done by Walecka who introduced a hadronic model consisting of nucleons, scalar mesons and vector mesons [59] and by Lee and Wick who reformulated the model of Walecka as a chirally invariant version [60]. The latter model is often called the chiral Walecka model and basically corresponds to the quark-meson model with nucleon instead of quark fields. The problem here is that the baryonic degrees of freedom become essentially massless in the chirally restored phase (Lee-Wick matter) as in these models their mass, in the chiral limit, is entirely generated by spontaneous chiral symmetry breaking. The occurrence of a Lee-Wick phase is circumvented by including the parity partners of the nucleons in a chirally invariant way, leading to parity-doublet models.
For an FRG treatment of the PDM we need an ansatz for the corresponding effective average action Γ k , the central object in the FRG approach formulated by Wetterich [61], where k is the renormalization-group scale. In this work we will use the following ansatz for the effective average action of the PDM, extended by vector and axialvector mesons, thus combining the FRG framework for the PDM presented in [40] and the strategy to include massive spin-1 (axial-)vector mesons presented in [58], The nucleon fields N 1 and N 2 are defined to have opposite parity and respectively represent the iso-doublet of nucleons, (p, n), and their parity partners, to which we assign the N * (1535). The chirally-invariant bare nucleon mass is given by m 0,N ; µ B denotes the baryon chemical potential, and the h's label the various Yukawa couplings between mesons and baryons. In this work we choose the scalar and vector-couplings to be the same, i.e. h s,1 = h v,1 and h s,2 = h v,2 , as also done in [57,58]. The scalar and pseudo-scalar meson fields are combined in φ 2 = σ 2 + π 2 with φ = (σ, π) T . U k (φ 2 ) is the O(4) symmetric effective potential, and the term cσ provides the explicit chiral-symmetry breaking that arises from the small but finite current masses of the light quarks in perturbative QCD. In principle, the effective potential and hence the thermodynamic grand potential can depend on all field combinations allowed by symmetry. The guiding principle here is to include fluctuations due to collective excitations such as those of order-parameter fields as in Landau-Ginzburg-Wilson effective theories. Because neither the ρ nor the a 1 meson are expected to develop non-vanishing expectation values in symmetric nuclear matter, their fluctuations are not included in the effective potential, at this level. This ansatz represents the leading order in a derivative expansion, also called local potential approximation (LPA) [62,63].
To describe the dynamics of massive vector and axialvector fields and their couplings in an effective theory [58], right and left-handed vector mesons are first introduced as (anti-)selfdual field strengthsρ ± µν = ±ρ ± µν which transform according to the (1, 0) and (0, 1) representations of the Euclidean O(4) replacing the proper orthochronous Lorentz group for massive spin-1 particles, with (1, 0) ↔ (0, 1) under parity, Here, T R and T L denote the so(4) Lie algebra matrices, see App. A for explicit expressions and conventions, of the generators of the chiral SU (2) L × SU (2) R in the ad-joint representation. 1 The iso-triplet vector ρ µ and axial-vector a 1µ fields with common mass m v are then obtained from these field strengths as where T V = T R + T L and T A = T R − T L . These represent conserved four-vector fields by construction, and, in particular, no π − a 1 mixing arises because the direct derivative coupling of the pion ∝ ∂ µ π with the conserved a 1µ field vanishes [64]. The interactions of the (axial-)vector fields, combined in the so(4) matrix V µ = ρ µ T V + a 1µ T A , with the SO(4) vector of scalar and pseudo-scalar σ and π mesons are determined from minimal coupling with D µ = ∂ µ + igV µ , see App. A.

B. Flow of the Effective Potential and Numerical Implementation
The ansatz for the effective average action Γ k formulated in Eq. (1) is now used in the Wetterich equation [61] which defines the 'flow' of Γ k and is given by where R k is a regulator function that suppresses momentum modes with momenta smaller than k, Γ second functional derivative with respect to the fields, and the supertrace runs over all internal indices, over bosonic and fermionic field space, in momentum space including an integration over internal momenta or thermal Matsubara sums as well as the fermionic minus signs and factors of two. At the ultraviolet (UV) scale k = Λ, Γ k is essentially given by the bare action. By solving the Wetterich equation and lowering the scale k the effects of quantum and thermal fluctuations are gradually included until the full effective action Γ = Γ k=0 is obtained in the limit k → 0.
The regulator function R k has to be chosen appropriately for different types of fields [65]. In this work we use three-dimensional regulator functions that only regulate spatial momenta but not the energy components, at the expense of slightly breaking the Euclidean O(4) symmetry [42]. While in principle also four-dimensional regulator functions can be used [66,67], the three-dimensional regulators allow to analytically perform the integration over the internal energy component, or the corresponding Matsubara sum at finite temperature, as included in the supertrace of Eq. (6). This in turn allows to apply the aFRG analytic continuation procedure as necessary for the calculation of the real-time two-point functions and spectral functions in the following. Explicit expressions for the different regulator functions are given in App. A.
When inserting the ansatz (1) into the Wetterich equation (6), one obtains the flow equation for the effective potential, Therein, we introduced the number of flavors N f = 2, the bosonic and fermionic occupation numbers, n B and n F , as given explicitly in App. A, and the scale-dependent particle energies for the sigma meson, the pion as well as for the nucleon and its party partner. The effective quasi-particle energies are defined as with the effective masses of the pion an the sigma meson given by where primes denote derivatives with respect to the chirally invariant φ 2 ≡ σ 2 + π 2 , and φ 2 0 = σ 2 0 is the global minimum of the effective potential in the IR. The masses of the nucleon and its party partner are defined by the eigenvalues of the mass matrix and are given explicitly by As expected, for σ 0 → 0 the two masses become degenerate and reduce to the bare mass m 0,N , while for m 0,N = 0 the model reduces to a sum of two independent fermion-meson models with masses for the nucleon and its parity partner given by h s,1 σ and h s,2 σ. In this case the fermion masses are generated by spontaneous chiral symmetry breaking as in the Lee-Wick model [60]. Note that the flow equation for the effective potential, Eq. (7), does not depend on quantities related to vector mesons. This is because, for isospin-symmetric matter with an equal number of protons and neutrons, the isovector ρ and a 1 mesons are not expected to develop non-vanishing expectation values and are therefore not included in the Ansatz for the effective potential. The isoscalar ω vector meson, in principle, contributes to the effective potential, see e.g. [40,68], but has not been included here. The inclusion of the ω meson as a dynamical field with finite expectation value inside nuclear matter, and fluctuations contributing to effective potential and thermodynamics will be deferred to future work.
It remains to specify the UV initial conditions and the parameters used to solve the flow equation for the effective potential. At the UV scale Λ we choose the effective potential to be of the form and then solve the flow equation numerically using the socalled "grid method", see for example [69]. This method is based on a discretization of the field φ in σ-direction while the pion field is set to its expectation value π = 0. Derivatives of the potential in field direction, as needed in the flow equation, are then obtained by a finite-difference scheme. We have checked explicitly that this numerical setup gives the same results as other approaches like finite-element or finite-volume methods, see also [70,71].
The numerical values for the parameters used for the effective potential, the bare nucleon mass, and the Yukawa couplings are summarized in Tab I. These are chosen such as to reproduce phenomenologically reasonable values for the pion decay constant and the particle masses in the vacuum at k → 0, where the pion decay constant f π is identified with the value of the sigma field at the global minimum of the effective potential. The resulting values for f π , the pion mass m π , the sigma mass m σ , the nucleon mass m N1 , and the mass of its parity partner m N2 are also given in Tab I. For the UV cutoff we use Λ = 1000 MeV and for the IR scale k IR = 40 MeV.
As usual in O(N )-Yukawa models, the symmetry breaking is generated from the fermionic fluctuations, where the fermionic minus sign acts like a negative index of refraction to drive the expectation value of the scalar order-parameter field away from its symmetric minimum. Although the fermionic fluctuations arise here from baryons with a sizable bare mass of m 0,N = 800 MeV, a UV cutoff of Λ = 1 GeV turns out to be just large enough to generate the right amount of sym- RG-scale dependence of the Euclidean particle masses in the vacuum. We observe strong effects from spontaneous chiral symmetry breaking between RG scales of k ≈ 900 MeV and k ≈ 600 MeV. The scale-dependent masses shown here serve as input for the calculation of the ρ and a1 spectral functions and determine the locations of the thresholds corresponding to the various decay channels. metry breaking starting from an effective potential with only the symmetric minimum for our choice of UV parameters. This mechanism of dynamical chiral symmetry breaking by the baryonic fluctuations is demonstrated in Fig. 1. Starting from the UV cutoff Λ the mass of the nucleons first decreases while that of the 1/2 − baryons increases. Incidentally, the chiral symmetry breaking scale at k χ ∼ 850 MeV here coincides with the scale m N2 ∼ k at which the heavier 1/2 − baryons decouple. The nucleon mass starts increasing at this point so that the fermionic fluctuations eventually cease to dominate the flow. A second scale around k ∼ 600 MeV then emerges below which the mesonic fluctuations eventually dominate. They tend to act symmetry restoring and this thus explains why the breaking pattern is not monotonously increasing with the flow towards the infrared where it levels at the desired physical values as it would in a purely mesonic model.
The fact that this can be achieved in this way, with no clear separation of scales between the initial fermion mass of m 0,N = 0.8 GeV and the UV cutoff scale Λ = 1 GeV might be surprising at first. It is reassuring for our effective hadronic theory which would otherwise start to lose credibility, if either considerably higher UV cutoff scales where needed or the fermionic fluctuations were irrelevant in the first place.

C. Phase diagram
In order to obtain the phase diagram of the PDM for the specified parameters we solve the flow equation for the effective potential at different combinations of temperature and baryon chemical potential, and plot the chiral order parameter, i.e. the value of σ 0 (µ B , T ) at the global minimum of the effective potential in the IR. The resulting phase diagram in the regime of high chemical potentials and comparatively low temperatures, as relevant for nuclear matter, is shown in Fig. 2.
As also found in [40], we observe two distinct phase transitions. The phase transition at lower chemical potentials represents the liquid-gas transition of nuclear matter while the second phase transition at higher chemical potentials inside dense nuclear matter can in the chiral limit be identified as the chiral phase transition. Both phase transitions consist of a first-order line at low temperatures connected to a critical endpoint (CEP). With our current parameters these CEPs are located at (µ B ≈ 896 MeV, T ≈ 10 MeV) and (µ B ≈ 925 MeV, T ≈ 33 MeV), respectively.
The position and the strength of the liquid-gas tran-sition strongly depend on the bare nucleon mass m 0,N . The larger the value of m 0,N the more does the location of the discontinuity move towards larger µ B while at the same time the strength of the transition gets weaker resulting in a larger in-medium condensate and hence a smaller nucleon sigma term. Obtaining phenomenologically acceptable values for the binding energy per nucleon, the nuclear saturation density, and the correct inmedium condensate all at the same time is known not to be possible within the present FRG setup [40]. This will require the proper inclusion of fluctuating ω mesons which we defer to a future study. For our first qualitative study here, we chose m 0,N = 800 MeV as a reasonable compromise as concluded in [40].
At the mean-field level, the inclusion of the ω meson in the effective action is known to result in a simple shift of the chemical potential, and is hence effective in adjusting the binding energy per nucleon essentially without influence on the strength of the nuclear liquid-gas transition. That the same mechanism, from a meanfield gap equation for the ω meson, does not work within the FRG framework for the order parameter fluctuations used here, was shown in [40]. The present effective theory framework for fluctuating vector mesons from [58] can in principle be extended to include also the repulsive contributions from fluctuating ω mesons in the flow for the effective potential and hence the thermodynamic grand potential to improve this situation. This requires further technical developments, however, and is therefore left for future work. Another issue is the slope of the first-order lines at low temperatures. As pointed out in [72], from a Clausius-Clapeyron relation, a positive slope dT c /dµ c of the first-order line implies a negative jump in the entropy density when going from the gaseous to the liquid phase. While this by itself is not necessarily unphysical at finite temperature, when the magnitude of the jump gets too large it leads to negative entropy densities on the liquid side of the transition line which then certainly is unphysical. The inclusion of a scale-dependent gap equation for a mean-field description of the ω meson was recently shown to be able to remedy the analogous unphysical effect in the quark-meson model [68]. It might therefore be reasonable to expect that it will also be affected when the ω-meson fluctuations are properly included within our FRG framework for a more realistic description of the thermodynamics of nuclear matter from the PDM in the future, as discussed above.

D. Vector and Axial-Vector Meson Propagators and Masses
We now turn to the calculation of the Euclidean vector meson masses, which will use the scale-dependent effective potential as input. In Ref. [58] it was shown that the vector-meson part of the effective action in Eq. (1), with Eqs. (3) and (4) corresponds to tree-level two-point functions for free (axial-)vector mesons with mass m v of the form, It was furthermore explicitly verified that this form, upon analytic continuation, in the interacting theory correctly describes that of the corresponding single-particle contributions to the spectral representations of the propagators of massive (axial-)vector fields which are related to the analogous current-current correlation functions by current-field identities [58]. For the inclusion of (axial-)vector fluctuations of this form within the FRG we also follow the strategy of Ref. [58] and temporarily add artificial longitudinal terms in order to be able to invert the correlation functions. This then leads to an ansatz for the scale-dependent (axial-)vector propagators which reads as follows, with a regulator shape function r(y) and y = p 2 /k 2 as defined in App. A, the transverse and longitudinal projection operators, and a scale-dependent mass parameter which differs from the running vector-meson pole mass m v,k in that it includes an equally scale-dependent wavefunction renormalization factor Z k . At the UV cutoff, we start with Z Λ = 1 and typically m 0,Λ = m v,Λ ≈ Λ, i.e. the common pole mass of the transverse vector and axial-vector fluctuations starts out at an initial value of the same order as Λ in the UV. The UV mass of the corresponding longitudinal fluctuations, ξm 2 v,Λ , is of the same order. The dimensionless parameter ξ can be introduced to further suppress these initial longitudinal fluctuations, it is here chosen as ξ = 10. Because of the additional factor of Λ 2 /k 2 this longitudinal mass then quickly increases with lowering k and the unphysical longitudinal fluctuations decouple. Varying the parameter ξ, one furthermore verifies that the results are in fact completely independent of these minute longitudinal fluctuations which strictly speaking violate the current conservation laws and require modified Ward identities at finite k. In the limit k → 0 the propagators become purely transverse again, however, and thus fulfill the usual Ward identity For further details, see Ref. [58]. Explicit expressions for the transverse and longitudinal regulators are given in App. A. In order to calculate the masses of the (axial-)vector mesons we first need to solve the flow equation for the scale-dependent vector mass m v,k , from which we obtain the ρ and a 1 masses as The flow of the vector mass m v,k can be obtained from the flow equation for m 0,k by observing that the flow of the product of m 2 v,k and m 2 0,k vanishes, We thus get with the flow of m 0,k given by The required flow equations for the two-point functions can be obtained by taking two functional derivatives of the Wetterich equation with respect to the appropriate fields which leads to the general structure where explicit expression for the three-and four-point vertices are given in App. A. The resulting flow equation for the ρ and the a 1 two-point function is represented diagrammatically in Fig. 3. Here, we take all loops into account that have up to one internal vector meson, as also done in [58]. This in particular gives rise to the process FIG. 4. Euclidean particle masses as a function of baryon-chemical potential, µB, at fixed temperatures of T = 10 MeV (left) and T = 33 MeV (right). Left: We observe the effects of the liquid-gas CEP at µB ≈ 896 MeV, where the sigma mass decreases rapidly, followed by the discontinuities generated by the first-order phase transition at higher chemical potentials. Right: The liquid-gas CEP still affects the behavior of the masses in a crossover from low to high density at µB ≈ 890 MeV while the chiral CEP gives rise to strong modifications at µB ≈ 925 MeV. At both temperatures chiral symmetry is restored to a large extent for chemical potentials beyond the second phase transition, giving rise to degenerate masses of chiral partners. a 1 → ρ+π which is important to describe the a 1 spectral function. We also note that in the following we will only be dealing with the transverse two-point function which can be obtained from the corresponding flow equation, The flow equations for m 0,k and m v,k are then solved using the scale-dependent effective potential and its derivatives at the IR minimum, σ 0,IR , as input. As initial condition we use m v,Λ = m 0,Λ = 1081 MeV and g = 8.78 for the dimensionless coupling of (axial-)vector to (pseudo-)scalar mesons (in their covariant derivative). These parameters are chosen such as to obtain phenomenological values for the ρ and the a 1 pole mass, as discussed in the following.
Our results for the flow of the Euclidean masses of the pion, the sigma meson, the ρ and the a 1 , as well as for the nucleon and its parity partner in the vacuum are shown in Fig. 1. Here we evaluate the masses at the scale-dependent minimum of the potential and not at the fixed IR minimum in order to show the effects of chiral symmetry breaking more clearly. Starting at the UV scale where chiral symmetry is restored, the masses of the chiral partners π − σ, ρ − a 1 , and N 1 − N 2 are degenerate. Taking fluctuations into account by lowering the scale k, we observe the effects of chiral symmetry breaking with the masses splitting up. At the IR scale these Euclidean mass parameters then arrive at the values for (pseudo)scalar mesons and nucleons listed in Tab. I, together with m ρ ≈ 840 MeV and m a1 ≈ 1170 MeV for the vector and axial-vector mesons. Note that these mass parameters do not represent the physical masses of the ρ and the a 1 which are in turn given by the pole masses obtained from their aFRG flows and result with these parameters to be m p ρ ≈ 775 MeV and m p a1 ≈ 1230 MeV, see below. We now turn to the dependence of the Euclidean particle masses on temperature and chemical potential. In Fig. 4 we show these masses at temperatures of T = 10 MeV and T = 33 MeV, i.e. at the temperatures of the two CEPs, as a function of baryon chemical potential µ B . We find that the masses are almost constant for a wide range of chemical potentials at these low temperatures, as expected from the Silver Blaze property [73]. Very close to the CEP at µ B ≈ 896 MeV and T ≈ 10 MeV, however, we observe drastic changes in the masses. In particular the sigma mass strongly decreases at this second-order phase transition. In principle, the sigma meson should become exactly massless here, since it is connected to the critical long-range correlations in the density fluctuations.
At T = 10 MeV and chemical potentials larger than the critical value of the liquid-gas CEP, one then encounters the discontinuous behavior generated by a firstorder transition. At even lower temperatures this second, chiral, phase transition is stronger than the nuclear liquid-gas transition, in that the chiral condensate and the masses change by a larger amount, cf. also Fig. 2. We also note that chiral symmetry becomes almost completely restored here, as evident from the fact that the masses of the chiral partners become almost degenerate. At T ≈ 33 MeV we only see smooth changes of the masses at chemical potentials near µ B ≈ 896 MeV while near the second CEP at µ B ≈ 925 MeV we again observe strong but continuous changes, as expected. These masses and their scale-dependence will serve as input for the calculation of the vector-meson spectral functions discussed in the following.

A. Analytic continuation and real-time two-point functions
So far, we have been working in Euclidean spacetime. In order to calculate real-time quantities like retarded two-point functions and spectral functions we will now perform an analytical continuation of the flow equations for the two-point functions from imaginary to real energies using the previously introduced analyticallycontinued FRG (aFRG) technique [42][43][44]. This technique utilizes the fact that the FRG flow equations always have a one-loop structure and that therefore the well-known analytic continuation procedure for one-loop calculations can be applied, see e.g. [74,75].
To be more specific, the flow equations for the twopoint functions, cf. Eq. (27) are analytically continued from imaginary to real energies by first using the periodicity of the bosonic and fermionic occupation numbers, which result from the Matsubara summation over the loop energy, w.r.t. the discrete external Euclidean energy p 0 , i.e.
In a second step, the Euclidean energy p 0 is replaced by a continuous real frequency ω in the usual way, where we use a small but finite value of = 0.1 MeV in our numerical implementation. The resulting flow equations for the retarded two-point functions are then solved using the scale-dependent effective potential as well as the flow of the masses m v,k and m 0,k evaluated at the IR minimum σ 0,IR as input.
The initial values of the retarded two-point functions are given by This initial shape will change as fluctuations are included by solving the flow equations. In particular, the twopoint functions will show sudden changes and thresholds induced by the different processes possible in our setup. The ones representing a decay of an off-shell meson into two on-shell particles can occur even in the vacuum while the processes involving an additional particle in the initial state are only possible at finite temperature and/or density when a thermal medium of such particles is available. We note that the usual kinematic constraint for a decay process like ρ * → π + π is given by ω ≥ 2m π while for a thermal capture process like a * 1 + π → σ we must have ω ≤ m σ − m π .
For an off-shell rho-meson ρ * with energies ω up to 2 GeV the relevant processes here are ρ * → π + π, (33) ρ * → a 1 + π, ρ * + π → a 1 , Note that ρ * + a 1 → π is in principle also possible, at very large values of the baryon chemical potential, when eventually m π > m a1 , cf. Fig. 4. Other processes either require higher energies or involve anti-baryons in the heat bath which are exponentially suppressed. For the a 1meson we have, analogously, where again also the more exotic processes a * 1 + ρ → π and a * 1 + a 1 → σ are possible in those regions of the phase diagram at very large µ B where the (pseudo-)scalar masses increase beyond those of the (axial-)vectors. More importantly, however, when the sigma mass drops below the pion mass, we can also have This occurs only in the critical regions close to both CEPs, and it might hence serve as a potential signature of the existence of the chiral CEP, in particular. In Fig. 5 we show the IR result for the real parts of the ρ and the a 1 two-point functions in the vacuum. We use the zero-crossings of these functions as an approximation for the pole masses of the respective resonances, which we fix to the phenomenological values of m p ρ ≈ 775 MeV and m p a1 ≈ 1230 MeV, cf. [76], by adjusting the parameters to the values listed in Tab. I and in Sec. II D. We also observe the effects from different vacuum decay processes which give rise to sudden changes in the two-point functions.
In Fig. 6 we show the imaginary parts of the ρ and the a 1 two-point functions in the vacuum, as well as at T = 33 MeV and µ B = 890 MeV as a representative example for medium effects in a region where both endpoints might have an influence at the same time, that of the nuclear liquid-gas transition, here as a crossover from low to high density, as well as the chiral CEP approached at this temperature for higher µ B .
The imaginary part of the ρ-meson two-point function in the vacuum is practically zero below the lowest decay threshold which is determined by the process ρ * → π + π and is located at ω ≈ 280 MeV.
In particular for the nucleon-(anti)nucleon threshold, but also in ρ * → a 1 + π, we observe the effects of the finite value for which produces small imaginary parts and hence results in non-zero values of the spectral function even below the thresholds. We note that the flow equations for the imaginary parts of the retarded twopoint functions can in principle also be solved for = 0 exactly, see e.g. [57]. Here, however, we have used a finite value of in order to keep the numerical results consistent with the ones for the real parts, where the limit → 0 is not so straightforward to take numerically with the principal value prescriptions that are involved there. Moreover, a small value for has the additional benefit of mimicking a phenomenological two-particle-into-twoparticle scattering continuum at arbitrarily low energies. At T = 33 MeV and µ B = 890 MeV, the imaginary part shows additional changes due to capture processes involving particles from the thermal medium, in particular from ρ * + π → a 1 and ρ * + N 1 → N 2 .
For the a 1 two-point function we observe in principle the same behavior as for the ρ two-point function, especially in the vacuum where the lowest threshold starts at a * 1 → π + σ ≈ 610 MeV. This here represents the σ-resonance contribution to the three-particle decay a 1 → 3π which should of course start somewhat lower, at about 420 MeV, reflecting the smoother onset of the 3π continuum. To achieve this, one needs to feed the result for σ two-point function as the broad 2π resonance back into the aFRG flow equations in a fully selfconsistent calculation in the future.
At finite T and µ B an increased number of processes and hence more complicated structures arise. In particular, we also observe a van Hove-like peak at ω ≈ 35 MeV in the contribution from a * 1 + π → σ which originates from an approximate saddle point in the difference of the corresponding scale-dependent quasi-particle energies E σ,k − E π,k , see also [57]. Another interesting effect is visible in the contribution from the a * 1 → a 1 + σ process where a peak near the threshold is forming at ω ≈ 1370 MeV. This enhancement is due to the dropping sigma-meson mass in the crossover region above the nuclear liquid-gas transition. These structures, as all the effects visible in the imaginary parts of the two-point functions, directly translate into the shape of the spectral functions discussed in the following.

B. Spectral functions
In this section we present our results for the in-medium ρ and a 1 spectral functions in nuclear matter. The spec- The ρ spectral function shows a prominent peak at its pole mass of m p ρ ≈ 775 MeV while the a1 spectral function exhibits a rather broad maximum which is strongly influenced by the decay a1 → ρ + π, cf. Fig. 6. tral function is generally defined as the discontinuity at the cut in the propagator along the timelike invariantmomentum axis and hence given by the imaginary part of the retarded Greens function G R , which can be expressed in terms of the retarded two-point function as (37) In this work, we set the external spatial momentum p = 0 to zero which makes an additional splitting of the spectral functions into a part transverse and longitudinal to the medium unnecessary.
In Fig. 7 we show the ρ and the a 1 spectral functions in the vacuum. The ρ spectral function shows a prominent peak at its pole mass of m p ρ ≈ 775 MeV and a full width of Γ ≈ 100 MeV. The only process contributing in this energy regime is the decay into two pions, ρ * → π + π. Comparing to the experimental width of 147.5 MeV for the charged ρ into two pions [76] our decay width is somewhat small but of the right order. At higher energies the decay channels ρ * → a 1 + π and ρ * → N 1 +N 1 give rise to additional thresholds at around 1300 MeV and 1880 MeV, cf. Fig. 6. The a 1 spectral function shows a broad maximum between ω ≈ 1000−1500 MeV where the width is due to the processes a * 1 → σ + π and a * 1 → ρ + π. At higher energies we observe the a * 1 → a 1 + σ threshold while the a * 1 → N 1 +N 1 contribution is very small below ω ≈ 2 GeV, cf. Fig. 6. We note in particular that this is the first time that the ρ and a 1 spectral functions have been obtained within an aFRG setting without suffering from unphysical decay thresholds into quark-antiquark pairs. This is one of the reasons why we are using the hadronic effective theory which contains nucleons and their parity partners in the place of the quarks in chiral quark models such as the Nambu-Jona-Lasinio or quarkmeson models, no matter whether these are enhanced by Polyakov-loop variables to model confinement or not.
In Fig. 8 we show the ρ and a 1 spectral functions at different temperatures and baryon chemical potentials. As also observed for the two-point functions, the spectral functions are essentially independent of µ B at low temperatures until very close to the phase transition. The size of this 'critical regime' is indeed very small, as evident from the large changes between µ B = 895 MeV and the CEP at µ B = 896 MeV at T = 10 MeV. Here, we observe some effects from the capture processes a * 1 + π → ρ and a * 1 + N 1 → N 2 at lower energies as well as an additional peak structure forming in the a 1 spectral function at ω ≈ 1380 MeV. This effect is due to the process a * 1 → a 1 + σ which is strongly affected by the dropping of the sigma mass at this second-order phase transition.
At T = 33 MeV the effects from the capture processes can be observed more clearly, both at µ B = 890 MeV MeV one sees the appearance of an additional peak structure in the a1 spectral function at ω ≈ 1380 MeV which stems from the a * 1 → a1 + σ process, with the σ meson encoding the critical behavior. At higher temperatures of T = 33 MeV the effects from capture processes become more pronounced, e.g. from the a * 1 + N1 → N2 process at ω ≈ 500 MeV (bottom left). At higher chemical potential (bottom right) we observe additional peak structures arising, this time also in the ρ spectral function, as well as a progressing degeneration of the spectral functions due to the restoration of chiral symmetry. See text for details. and µ B = 915 MeV, e.g. near ω ≈ 500 MeV, mainly due to the baryonic capture processes ρ * /a * 1 + N 1 → N 2 . Closer to the chiral CEP, at µ B = 915 MeV, we again observe the formation of an additional peak structure in the a 1 spectral function due to the critical effects entering via the process a * 1 → a 1 + σ. A similar peak structure is also visible in the ρ spectral function which, however, originates from an additional zero-crossing forming in the real part of the two-point function which is connected to the process ρ * → π + π. We also observe that the ρ and the a 1 spectral functions become increasingly degenerate due to the progressing restoration of chiral symmetry. In fact, one can show analytically that the flow equations of the ρ and the a 1 two-point functions become degenerate in the limit σ 0 → 0.
The most relevant low-energy processes at T = 33 MeV and µ B = 924 MeV, near the chiral CEP, are shown in Fig. 9 where the different contributions to the imaginary parts of the ρ and a 1 two-point functions are plotted up to 1 GeV (higher energies become increasingly difficult to compute as manifest in unphysical sign changes that can occur for higher energies in this critical region). Although we can clearly identify the potential signature of criticality in the process a * 1 + σ → π as mentioned above, the strength of this signal in the a 1 two-point function below 100 MeV turns out to be very weak. The by far dominant low energy features in both two-point functions here are the contributions from the nucleon capture processes ρ * + N 1 → N 2 and a * 1 + N 1 → N 2 . These baryonresonance formation processes give rise to rather strong ] ω [MeV] a 1 +σ → π, a 1 → σ+π a 1 +π → ρ, a 1 → ρ+π a 1 +σ→ a 1 , a 1 → a 1 +σ a 1 +N 1 → N 2 , a 1 → N 1 +N 1 FIG. 9. Imaginary part of the ρ (left) and the a1 (right) two-point functions at T = 33 MeV and µB = 924 MeV, close to the chiral CEP. Here, a particularly small value of = 0.01 MeV was needed in order to be able to resolve weak low-energy contributions from capture processes such as the critical a * 1 + σ → π. As before, the separate components are extracted from the different loops shown in Fig. 3. Of all contributions to the imaginary parts discussed above, the most prominent medium modifications of the critical spectral functions are the baryon-resonance formation processes ρ * + N 1 → N 2 and a * 1 + N 1 → N 2 which give rise to pronounced low-energy peaks around ω ≈ 240 MeV, below all other thresholds. The occurrence of these peaks is a unique prediction of the baryonic mirror assignment and its observation through enhanced dilepton pair production in the vicinity the chiral CEP would be an important confirmation of this picture of mass generation in QCD. The critical capture process a * 1 + σ → π at even lower energies turns out by far too weak to be potentially significant. It is about six orders of magnitude lower than the baryonic capture process in the a 1 spectral function in Fig. 10. Other than that we observe only a small mass shift and broadening in the ρ spectral function with considerably stronger medium modifications near the quasi-particle peak in the a 1 , indicating the emerging restoration of chiral symmetry on the level of the eventually complete degeneration of the spectral functions of the chiral partners ρ and a 1 at high density.

IV. SUMMARY AND OUTLOOK
In the work presented here we discuss results on vector and axial-vector meson spectral functions at finite temperature and baryon-chemical potential, in order to functions near criticality, as mentioned above, lead to positivity violations in the spectral functions at higher energies, above 1 GeV. Whether these are related to thermodynamic instabilities observed in dense regions of the phase diagram with fluctuations [72] remains to be investigated. assess the impact of chiral symmetry restoration in dense, low-temperature nuclear matter on the redistribution of spectral strength in both channels. As low-energy effective theory we use a chiral baryon-meson model, namely a parity-doublet model, which contains pions, sigma mesons, ρ and a 1 mesons as well as nucleons and their parity partners chosen to be the N * (1535). Choosing hadronic degrees of freedom avoids unphysical quarkantiquark thresholds in the spectral functions in the confined phase. The vector and axial-vector mesons are introduced using a novel FRG formulation for massive vector fields based on (anti-)selfdual field strengths [58]. In our opinion, this extended parity-doublet model captures the essential features of mass generation in QCD, in that hadron masses only partially result from the spontaneous breaking of chiral symmetry. On the other hand, the degeneracy in the spectral functions of parity partners in the restored phase is entirely driven by the evolution of the chiral condensate. The effects of thermal and quantum fluctuations are taken into account by using the FRG approach, which is known to describe phase transitions and critical phenomena in a way that is superior to a thermodynamic mean-field description mainly by including the dynamics of order-parameter fluctuations due to collective excitations.
Within this theoretical setup we have calculated the phase diagram for isospin-symmetric nuclear matter as a function of T and µ B as well as the screening masses of the various hadrons involved. The distinctive feature of the model is that it exhibits a nuclear liquid-gas phase transition as well as a chiral phase transition at a higher chemical potential where the nucleons and their parity partners become approximately degenerate but remain massive. Similar to the chiral partners ρ and a 1 in the vector-meson channel, the splitting between nucleons N and their parity partners N * gradually disappears as chiral symmetry gets restored at finite density by predominantly the resonance mass dropping down to their common chirally invariant mass from the scale anomaly. Near the two first-order phase transitions as well as at their respective CEPs the Euclidean mass parameters of mesons and baryons all show the expected behavior and serve as input for the evaluation of the spectral functions.
For the calculation of the real-time two-point functions and the spectral functions we used the so-called aFRG method, of solving analytically continued FRG flow equations. Within this method one performs the analytic continuation from imaginary to real frequencies directly on the level of the FRG flow equations for the two-point functions and thus avoids the need for any numerical reconstruction. Moreover, it is thermodynamically consistent in that the effective potential and the spectral functions are calculated on the same footing. Using this approach, we have calculated the ρ and the a 1 spectral functions at different temperatures and baryon-chemical potentials.
In the vacuum, the ρ spectral function shows a prominent peak whose width is solely determined by the de-cay into two pions, as expected. The a 1 spectral function, in contrast, exhibits a very broad peak at higher energies which is determined by the decay into a pion and a sigma meson as well as into a rho meson and a pion, representing the σ and ρ-meson resonance contributions to its three-pion decay width. For small temperatures and chemical potentials, below the liquid-gas phase transition, the spectral functions essentially coincide with those in the vacuum, as expected from the Silver-Blaze property. In the vicinity of the two critical endpoints, however, we observe significant changes with additional peaks emerging. Most strikingly are the modifications near the chiral CEP, were a prominent low-energy peak around 240 MeV shows up. Its origin can be traced back to the resonance excitation of the in-medium N * (1535) in both the vector and axial-vector channel. In fact, the excitation strength is nearly identical, reflecting the signature of parity doubling. Preliminary estimates indicate that this effect, which is strongest in the vicinity of the chiral CEP and possibly near the first-order boundary of the chiral phase transition might be observed experimentally in the vector channel through an increased dilepton yield at correspondingly low invariant masses measured in heavy-ion collisions at a few GeV/nucleon with high statistics. Its detection would yield strong evidence in support of the parity-doubling scenario as providing the mechanism for chiral symmetry restoration inside dense nuclear matter.
To arrive at a satisfactory description of dense and warm nuclear matter and its spectral properties further improvements are called for. One relates to a quantitative description of the nuclear matter. For a phenomenologically acceptable description of the binding-energy per nucleon, the nuclear saturation density and the equation of state (EoS) of symmetric nuclear matter we will have to include the ω vector meson as a dynamical field in the calculation of the thermodynamic grand potential. From Walecka-type mean-field studies it is known that the repulsive nature of the ω meson is essential for a realistic nuclear matter EoS. Possible further improvements include taking into account higher truncation orders in the effective average action or using a self-consistent setup where the spectral functions are back-coupled into flows of effective potential and two-point functions. In addition, a calculation of the expected dilepton yields in heavy-ion collisions at a few GeV/nucleon is left for future work. Other phenomenologically important extensions will include the EoS of highly isospin-asymmetric nuclear matter and determination of neutrino emissivities relevant for binary neutron-star merger events from the corresponding weak (axial-)vector spectral functions in warm and dense neutron matter. work was supported by the Deutsche Forschungsgemeinschaft (DFG) through the grant CRC-TR 211 "Strong-interaction matter under extreme conditions" and the German Federal Ministry of Education and Research (BMBF) through grants No. 05P18RFFCA and No. 05P18RGFCA.

Appendix A: Explicit Expressions
In this appendix we summarize explicit expressions for various quantities used in this work. We begin with the regulator functions, for which we use three-dimensional Litim-type regulators [77], which allow for an analytic evaluation of Matsubara sums. More explicitly, we employ the following regulator functions for (pseudo-)scalar mesons, (axial-)vector mesons, and nucleons, R σ/π,k (p) = (k 2 − p 2 )Θ(k 2 − p 2 ), (A1) The corresponding regulator shape function is given by with y = p 2 /k 2 . Upon evaluating the (unregulated) Matsubara sums over internal energies in the various flow equations, we encounter the bosonic and fermionic occupation number factors which are given by The minimal coupling between (pseudo-)scalar and (axial-)vector mesons is then defined by the covariant derivative D µ = ∂ µ + igV µ with the matrix valued vector field V µ = ρ µ T V + a 1µ T A . Explicitly, this leads to The three-point interactions of order g then yield, ( π × ρ µ )·∂ µ π − σ a 1µ ·∂ µ π + a 1µ · π∂ µ σ , and the quartic interactions of order g 2 become Finally, we list the explicit expressions for the three-and four-point vertices which can be obtained by taking three and four functional derivatives of the ansatz for the effective average action with respect to the various fields, cf. Eq. (1). From the relevant contributions to the effective average action listed explicitly above, we obtain the three-point vertex functions involving vector mesons used in this work as where all momentum arguments denote the incoming momenta of the (pseudo-)scalar mesons, q µ for the sigma meson and q i µ for the isovector pions π i . The four-point vertices are analogously obtained as Γ (4) ρ i µ ρ j ν π k π l = g 2 δ µν (2δ ij δ kl − δ ik δ jl − δ il δ jk ) , (A16) Γ (4) a i 1µ a j 1ν π k π l = g 2 δ µν (δ ik δ jl + δ il δ jk ) , The three-point couplings of the (axial-)vectors to the two baryons doublets N d , d = 1, 2, are readily obtained from the ansatz for the effective average action and are given by