Nonperturbative quark matter equations of state with vector interactions

Nonperturbative equations of state (EoSs) for two and three quark flavors are constructed with the functional renormalization group (FRG) within a quark-meson model truncation augmented by vector mesons for low temperature and high density. Based on previous FRG studies without repulsive vector meson interactions the influence of isoscalar vector $\omega$- and $\phi$-mesons on the dynamical fluctuations of quarks and (pseudo)scalar mesons is investigated. The grand potential as well as vector meson condensates are evaluated as a function of quark chemical potential and the quark matter EoS in $\beta$-equilibrium is applied to neutron star (NS) physics. The tidal deformability and mass-radius relations for hybrid stars from combined hadronic and quark matter EoSs are compared for different vector couplings. We observe a significant impact of the vector mesons on the quark matter EoS such that the resulting EoS is sufficiently stiff to support two-solar-mass neutron stars.


Introduction
With the recent dawn of multimessenger astrophysics new data will become available with the possibility to further scrutinize models of the structure of compact objects. In particular, several observational projects are underway or planned to pin down the compact star equation of state (EoS) not yet fully known to date. Let us mention for example precise mass determinations from pulsar timing with current [1][2][3][4] and future instruments such as the SKA [5], observations of binary neutron star (BNS) mergers by the LIGO/Virgo collaboration [6][7][8] introducing constraints on the EoS via the measurement of the tidal deformability of the inspiraling stars [7,6], and measurements of the neutron star's x-ray emission giving information on masses and radii [9,10], see also the recent NICER results [11,12]. Post-merger oscillations are very sensitive to the dense matter EoS, too [13], and can among others indicate the existence of a strong first-order phase transition in compact star matter [14], such that future observations by the LIGO/Virgo/Kagra collaboration or via various projects such as the Einstein Telescope or the Cosmic Explorer are promising tools to reduce uncertainties.
However, despite this bright future on the experimental side, a precise determination of the underlying EoS will not be fully conclusive because on the theoretical side the EoS alone cannot resolve for the detailed matter composition and interactions. Only a full theoretical understanding of dense QCD will enable a satisfactory final picture which entails the need for a quantitative EoS grounded in first-principle QCD.
This theoretical challenge is aggravated by the confinement property of QCD: due to the running of the QCD gauge coupling the fundamental degrees of freedom lose their significance at low temperatures and densities and are confined into colorless, composite states, the hadrons. After confinement, a relatively small residual nuclear interaction is left which binds the nucleons into atomic nuclei with a typical binding energy per nucleon of the order of 1-10 MeV. This energy scale is two orders of magnitude smaller than the confinement scale but still very strong. Decades of considerable effort together with constraints from experimental data on nuclei and theoretical calculations have led to reliable models for the NS crust and homogeneous nuclear matter up to roughly twice the saturation density, ρ 0 ∼ 0.16 fm −3 , see e.g. [15][16][17] and references therein for a discussion. Above this density, not only the models are less under control, but non-nucleonic degrees of freedom might appear and the situation becomes more complicated.
Presently, the composition and in particular the phase structure of the interior of compact stars is not known. In addition to the possible onset of non-nucleonic hadronic degrees of freedom, such as hyperons [18][19][20], mesons, or ∆-baryons [21], quark matter could exist in the cores of compact objects where due to extreme gravitational forces densities of up to ≈ 8-10ρ 0 are expected [22,23]. At extremely high densities ρ > 60ρ 0 perturbative QCD can be used to guide reliably the construction of the EoS. However, the density of matter relevant for the description of compact stars is not high enough for a perturbative treatment to be conducive. Alternative approaches for a proper theoretical description of this intermediate density regime from low-energy QCD are inevitable. Unfortunately, ambitious first-principle lattice QCD simulations are not applicable in this regime due to a generic sign problem at finite densities. Therefore, so far mainly either parameterised interpolations between the lowdensity nuclear part and the perturbative QCD regime have been used, e.g. [24,23], generic parameterisations of the quark matter part, e.g. [25][26][27][28] or phenomenological approaches, such as bag models, Nambu-Jona-Lasinio (NJL)-type models or the quark-meson coupling model in mean field approximation have been employed to describe quark matter in this intermediate density regime, e.g. [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45]. Some models include constraints from lattice calculations, see e.g. [46], and recently holographic approaches have been developed, e.g. [47][48][49]. It is obvious that still much effort is needed and it should in particular be stressed that the fundamental dynamics of fluctuations is mostly ignored in mean-field calculations or simple parameterizations, see e.g. [50,51].
An alternative approach for the elaboration of the EoS in the intermediate density regime is based on the functional renormalization group (FRG) method. It is a promising nonperturbative realization of Wilsonian renormalization group idea in the continuum and not limited to small couplings. The straight and controlled computation of the EoS from first-principle QCD quark and gluonic degrees of freedom becomes conceivable [52]. However, towards lower densities the quarks cluster into nucleons and the emergence of long-range correlations between nucleons will increasingly complicate the QCD-based FRG approach [53][54][55].
Within this framework, the impact of fluctuations in the (pseudo)scalar interaction channel on the compact star EoS has recently been studied [56]. It was found that the fluctuations decrease the sound speed of the quark matter even below the mean-field value of c 2 s = 1/3, leading to a rather soft EoS at high densities. This rather soft EoS does not allow for the construction of hybrid stars with a three-flavor quark core in agreement with present neutron star mass measurements [2][3][4]. Therefore, here we will extent the work of [56] and investigate the impact of additional vector meson interactions. While the (pseudo)scalars fields are integrated out within the functional renormalization group framework, the vector mesons are treated on a constant background level. This idea has already been employed to studies of the phase diagram at finite temperature, see e.g. [57][58][59][60].
Vector interactions are expected to add repulsion [61]. This can be seen from classical models of the nucleon-nucleon interaction, from phenomenological models allowing for hy-perons in the neutron star core, see e.g. [62,63], or from many phenomenological quark matter studies, see e.g. [64][65][66]. The inclusion of a repulsive vector interaction in the quarkmeson model should thus stiffen the quark matter EoS, leading to a higher speed of sound and allowing for constructing hybrid stars compatible with observations.
Within this first study, we neglect the possibility of diquark pairing and do not enter the discussion of the extremely rich phase structure of color superconducting matter in the density range of neutron stars [67,68]. Diquark pairing is important for transport properties, but, being a Fermi surface phenomenon, has only little influence on the equation of state we are focusing on here.
As already mentioned for very low baryon density a neutron star can be characterized by nonrelativistic nucleons via nuclear forces while at high densities quarks have more and more of an impact such that the EoS in the low density regime can be constructed with techniques of nuclear matter theory. At high baryon density, due to the failure of first-principle QCD approaches, usually phenomenological quark models such as NJL-type models are adopted mostly on a mean-field level to incorporate the dynamics of quarks. In order to gap the bridge between these two density regimes nonperturbative effects as the generation of constituent quark masses by a spontaneous chiral symmetry breaking must be taken into account.
The work is organized as follows: In the following Sec. 2 a chirally symmetric effective model with quarks and mesonic degrees of freedom is introduced with scalar-pseudoscalar as well as axialvector-vector channels to account for quark matter including chiral symmetry breaking. Later, the focus lies mainly on the two and three lightest quark flavors relevant for the strongly interacting intermediate transition regime. The introduced model setup serves then as truncation for a nonperturbative integration of the quark and meson dynamics with the FRG in Sec. 3. The augmentation of the FRG with vector mesons is presented in Sec. 3.2. Now equipped with a nonperturbative EoS neutron star properties can be investigated. The numerical results are presented in Sec. 4 wherein the β-equilibrated EoS, mass-radius relation and tidal deformability with and without strangeness for pure quark matter stars as well as hybrid stars for various vector couplings are discussed. In Sec. 5 a summary with a conclusion is given. A list of the used input parameters and some technical details for the numerical solution are collected in the appendices.

Quark-meson model with vector mesons
A chirally symmetric effective theory is considered with N f flavors of constituent quark fields q and N 2 f (pseudo)scalar and (axial)vector meson fields. With the usual U (N f ) flavor transformation generators T a with a = 0, . . . , N 2 f − 1 the meson matrices can be rewritten as for the (pseudo)scalar mesons and for the (axial)vector mesons.
The mesons interact via a scalar g s and a vector g v Yukawa coupling with the quarks which is encoded in the quark-meson Lagrangian in Euclidean space The anomalous breaking of the axial U (1) A -symmetry is realized via 't Hooft determinants ξ = [det Φ † + det Φ] with a constant parameter c A . The mesonic (self-)interactions are parametrized with the chiral symmetric potential and is in general a function of N f independent chiral invariants ρ n . The highest chiral invariant ρ N f is usually omitted in the chiral potential because it often depends on the other invariants.
A small explicit breaking of chiral symmetry is implemented with a linear term in the meson fields Tr[H(Φ † + Φ)] and a corresponding symmetry breaking parameter matrix H.
We assume isospin symmetric matter such that the two lightest quark flavors up and down are degenerated in the masses and only one index l = u = d is needed. Hence, for three quark flavors N f = 3 the explicit symmetry breaking term in (2) reduces to two different terms proportional to −c l σ l − c s σ s with the strange quark flavor index s. A constant rotation establishes the relation of the non-strange-strange basis (σ l , σ s ) and the singlet-octet basis (σ 0 , σ 8 ) in the scalar meson sector such that the isospin-symmetric scalar vacuum condensates are diagonal in flavor space For N f = 2 quark flavors only one independent invariant ρ 1 in the chiral potential survives and hence one explicit symmetry breaking parameter is needed. Since the U (1) A -symmetry breaking term scales with the meson fields to the power N f the 't Hooft determinant corresponds to a mesonic mass term for N f = 2 and constitutes a mass splitting between the scalar and pseudoscalar meson multiplets {σ 0 , π 1 , π 2 , π 3 } ↔ {π 0 , σ 1 , σ 2 , σ 3 }. Hence, the parameter c A is dropped for N f = 2 and only the first multiplet corresponding to lighter mesons, i.e., the scalar resonance and the pseudoscalar Goldstone bosons ϕ = σ, π π π T T , is considered dynamically. In this way the axial U (1) A -symmetry breaking is assumed to be maximally broken. For more details see e.g. [69,70].
The generalization to finite temperature and baryonic densities is achieved within the Matsubara formalism, wherein the time-component is Wick-rotated t → −iτ and the imaginary time τ is compactified on a circle with radius equal to the inverse temperature β = 1/T . In general, N f independent quark chemical potentials µ f can be implemented in the quark part of the Euclidean Lagrangian, Eq. (2), by However, for cold neutron star matter a weak equilibrium with neutrinos that leave the star without further interactions is present such that not all chemical potentials are independent anymore. For three light quark flavors a common quark chemical potential µ can be introduced which is related to the baryon number via µ = µ B /3. Respecting the different electrical charges of the quarks by an additional electron chemical potential µ e which is the negative charge chemical potential one finds Taking electrical charge neutrality of the star into account yields for the quark and electrical densities such that only one independent chemical potential remains. Our choice for two and three quark flavors is the common quark chemical potential µ. Despite the fact that isospin symmetry is broken in the way the chemical potentials are introduced we still assume only one light condensate σ l as approximation. Consequently, even for isospin asymmetric matter the light quark masses are degenerated m u = m d = m l . For the following investigation of the equation of state (EoS) the total grand potential Ω is needed. It is given by the logarithm of the grand partition function, which in general is a path integral over all involved quantum fields and hence incorporates all quantum, thermal and density fluctuations of the studied system. In the literature usually the integration over the quark fields is performed whereas the dynamics of the remaining fields are drastically truncated and are taken into account on a mean-field level, see e.g. [35].
In this work we consider in addition the fluctuations of the remaining mesons via the FRG method which we will employ in the next Section to evaluate the grand potential at finite temperatures and densities and determine the EoS in a nonperturbative manner.

Functional Renormalization Group Approach
For a consistent implementation of successively integrating out quantum, density and temperature fluctuations from large to small energy scales we employ Wilson's functional renormalization group idea [71,72] in terms of the Wetterich equation [73] This flow equation is a functional differential equation for the evolution of the scale dependent effective action Γ k where the logarithmic scale derivative is denoted by ∂ t = k d dk . The effective average action Γ k interpolates between a microscopic or bare UV action S Λ = Γ k→Λ and the full quantum effective action Γ = Γ k→0 in the infrared. It thus governs the dynamics of the field expectation values after the integration of quantum fluctuations from the UV scale Λ down to the infrared scale k IR . The infrared regulator R k specifies the regularization of quantum fluctuations near an infrared momentum shell with momentum k. It is a diagonal matrix for mesons and symplectic for quarks. The scale-dependent IR regulator R k can be interpreted as momentum-dependent masses that suppress the infrared modes of the associated fields. The derivative term ∂ t R k in (10) ensures UV-regularity. The second functional derivative of the effective average action with respect to the fields of the given theory is generally denoted as Γ (2) k . The trace runs over all discrete and continuous indices, i.e., color, spinor and the loop momenta and/or frequencies. Since the full nonperturbative propagators enters in the flow equation the Wetterich equation is highly non-linear and includes higher loop contributions in perturbation theory despite its simple one-loop structure. The application of this approach is not restricted to the existence of a small expansion parameter and hence is applicable in any nonperturbative regime. For QCD-related applications of this approach see e.g. the reviews [74][75][76][77][78][79].
For the explicit solution of the functional equation some truncations are required that turn it into a system of finite-dimensional partial differential equations. Even though any truncation of a functional equation might induce a certain dependence of physical observables on the employed regulator the impact can be minimized by choosing optimized regulators or by implementing RG consistency [80]. Here, a modified three-dimensional flat regulator has been used [81,82].
For the solution of the flow equation (10) an initial condition at the UV scale k = Λ must be supplemented. A complete solution of the entire Lagrangian is presently beyond the scope of this work such that we simplify the system in the following way: The full dynamical fields in the flow equations are the quarks and the (pseudo)-scalar mesons which affect the effective potential and the size of the remaining vector mesons. For three quark flavors the UV initial condition for the flow equation (10) for Γ k=Λ reads wherein the effective potential U (2+1) k=Λ depends on the first two chiral invariants ρ 1 and Here ρ 2 has been shifted by ρ 2 1 /3 for computational simplicity. For two quark flavors, the effective potential U (2) k=Λ (ρ 1 ) depends only on one chiral invariant. This truncation for the effective action corresponds to a leading-order derivative expansion with standard kinetic terms for the (pseudo)scalar meson fields. In this local potential approximation (LPA) of the initial action no scalar wave function renormalizations and no scaledependence in the Yukawa couplings between quarks and mesons are taken into account. Note that in this simple truncation the important dynamical back-reaction of the mesons on the quark sector of the model is already included. For more details see the literature, e.g. [69,70].
Plugging the truncation (11) of the action into the Wetterich equation (10) yields finally an IR and UV finite flow equation for the effective potential where the flow of the (pseudo)scalar mesonic degrees of freedom is fully taken into account. The single-particle energies E b = k 2 + m 2 b include the RG scale dependent (pseudo)scalar meson masses m b which are obtained by diagonalizing the mass entries of the matrix Both the potential U (N f ) k and the mass matrix (13) are evaluated at the vacuum expectation value given in Eq. (5). Details and the lengthy explicit expressions of the eigenvalues can be found in the literature [69,70]. The corresponding single-particle energies for the quarks are E u = E d ≡ E l = k 2 + g 2 s σ 2 l /4 and E s = k 2 + g 2 s σ 2 s /2, respectively. The full thermodynamic potential evaluated at the solution of the gap equation, i.e., the minimum of the grand potential is obtained by evolving the system towards the infrared.
As initial UV condition for the flow (12) the effective potential is parameterized for which contains the scale-dependent chiral effective potential and the 't Hooft determinants ξ evaluated for the light and strange condensates For three quark flavors only three scale-dependent expansion coefficients a ij in the chiral potential U (2+1) χ,Λ are needed and all remaining parameters are kept constant. For two quark flavors, the number of coefficients reduces to two because the second chiral invariant is not independent anymore.
From the RG point of view at the initial scale it is sufficient to take only relevant and marginal operators into account since meson fluctuations are small at high energies (due to their larger masses) and irrelevant operators are in addition dimensionally suppressed. However, this does not mean that irrelevant operators can be ignored in general. They are generated by the RG flow at smaller scales and are of relevance, see e.g. [83].

Implementation of vector mesons
The vector mesons are implemented on a mean-field level in the quark-meson Lagrangian (2), i.e., as static background fields such that their kinetic terms are not of further relevance anymore. Due to rotational symmetry all components of the vector meson condensates except the temporal ones are assumed to vanish [84]. For three quark flavors the only non-vanishing meson fields are in principle the diagonal scalar fields and the isoscalar vector fields ω and φ and the third isovector vector field ρ 3 0 . For asymmetric isospin matter and two quark flavors two scalar condensates σ 0 and σ 3 emerge and three scalar condensates σ i , i = 0, 3, 8 for three flavors in general. Due to technical reasons 2 we slightly simplify our approximation further by dropping the isospin breaking condensate σ 3 in the flow equation. Hence, to be consistent within our approximation scheme we also omit the isovector vector condensate ρ 3 0 since this introduces a further isospin asymmetry.
For three quark flavors this yields finally a diagonal matrix in flavor space for the ω-and φ vacuum expectation values This also assumes an ideal quark mixing such that the quark content of the ω-meson consists purely of up and down quarks while the φ-meson is purely strange. This in turn leads to the additional Yukawa-type coupling in the quark sector of the Lagrangian which can be interpreted as a shift in the corresponding chemical potentials, giving now rise to modified effective chemical potentials with g φ := √ 2g v . The constant vector meson vacuum expectation values also contribute to the mean-field potential wherein the negative sign expresses the repulsive nature of the vector interactions. The masslike parameters m 2 ω and m 2 φ are basically unconstrained and not to be identified with the physical vector ω-and φ-meson masses. They just inherit their names due to their mass dimension. However, we fix them to the measured vector meson masses m ω = 782 MeV and m φ = 1020 MeV such that the vector meson coupling g v is the only remaining free parameter of the system and of the order of one. For two quark flavors, only the ω-meson is accounted for. However, in the following we will establish the N f = 2 + 1 flavor equations and suppress the flavor index. The two flavor results arise in an obvious way.
This approximation can now be augmented with the FRG by adding the vector meson potential U vec to the scale-dependent chiral effective potential. For an arbitrary renormalization scale k, the total effective potentialŨ k reads Since the masses associated with the ω and φ bosons are related to the inverse range of the isoscalar short-distance NN interactions and are large compared to the relevant low-energy scales the ω and φ fluctuations should be more suppressed. Hence as a sort of inert degrees of freedom they can be treated as background fields. Contrarily, the fluctuations in the pseudoscalar channel (e.g. the sigma and the pions for two flavors) and also the particle-hole excitations of the quarks around the Fermi surface are fully taken into account. The respective condensates are determined in the infrared by solving the gap equations whereŨ IR denotes the fully evolved effective IR potential including all dynamic quark and (pseudo)scalar meson fluctuations. The condensates, i.e. the fields at the minimum of the potential, depend on the temperature and chemical potentials. The gap equations for the ω-and φ-condensates (22) can be rewritten The subscript "gap" in (23) labels the gap equation solution meaning that the potential is evaluated at these field configurations which solve the corresponding gap equations. Thus, both gap equations for the vector condensates (23) are self-consistent and can be solved numerically by root finding. Finally, the infrared potentialŨ IR evaluated on the gap equations can now be identified as the grand potential yielding the equation of state in a standard thermodynamic manner. For vanishing temperature the normalized pressure and energy density are given by with the quark number densities The quark number densities are defined as the derivative of the grand potential Ω(T, {µ f }) with respect to the corresponding chemical potentials. Since the implicit dependence of the infrared potentialŨ IR on the chemical potentials through the condensates vanishes by virtue of the gap equation, i.e.
where the ellipses represent similar derivative terms for all other condensates, we can identify the derivative terms in (23) with the quark number densities: Note that the gap equation (27) is solved including the full underlying nonperturbative contributions from the FRG in the (pseudo)scalar channel. This is in contrast to a similar two quark flavor FRG study [59] where the gap parameter for the isoscalar ω-condensate is evaluated from a mean-field flow which ignores the back-coupling of the FRG flow.
Furthermore, the inclusion of the vector mesons into the (pseudo)scalar sector appears solely by the replacement of the chemical potentials in (12) with the effective chemical potentials given in (19).

Numerical Results
In the following we will present our findings obtained with the FRG quark-meson truncation including isoscalar vector-mesons for two and three quark flavors. All results are obtained for β-equilibrated and electrical charge neutral matter. Since we are primarily interested in the physics of older neutron stars for which temperature effects can be neglected, all flow equations are strictly solved for vanishing temperature.
We have employed two different and complementary numerical solution strategies for the flow equations as explained in App. B. The two quark flavor results are obtained with an upwind finite difference scheme while for the three flavor calculations a two-dimensional grid of the two scalar field variables has been used. In principle, this enables us to estimate possible numerical artefacts. We found excellent agreement in particular at low chemical potentials with no strange quarks populated, showing the robustness of the numerical scheme. This can be seen in Fig. 1 where the isoscalar ω-(dashed and solid lines) and φ-meson (dash-dotted lines) condensates for three different vector couplings as a function of the quark chemical potential are shown. The difference in the ω-condensate with and without strangeness is negligible. From the gap equations (23) it is clear that the vector condensates are proportional to the (respective) number densities. Hence, for T = 0 the ω-condensate vanishes in the chirally broken phase where all occupation numbers are zero, and the φ-condensate (dashdotted lines) is zero until µ ≈ 400 MeV when the strange quark states get populated. Since an increase in a vector condensate means a decrease in the respective effective chemical potential(s), cf. Eq. (19), which in turn decreases the number density from its initial value, there is always a unique solution to the vector meson gap equations (23). Note that the assumption of one light chiral condensate for both the up and down flavors breaks down close to the chiral phase transition under the assumptions of β-equilibrium and charge neutrality [56]. Hence, only data points above the chiral phase transition in the light scalar sector with µ > 310 MeV are considered. The effect of the background isoscalar vector mesons on the equation of state is displayed in Fig. 2. In this figure the (normalized) pressure obtained with the FRG quark-meson truncation including the vector meson condensates is displayed for different vector couplings g v as a function of the corresponding energy density. As already mentioned, β-equilibrium and charge neutrality have been implemented. The two quark flavor EoS (solid lines) is stiffer than the corresponding EoS with strangeness (dashed lines) for energy densities beyond the onset of strangeness. This is not astonishing since an additional degree of freedom generally reduces pressure and thus softens the EoS.
Vector mesons contribute to the EoS with two effects. Firstly, since the vector meson potential (20) gives a negative contribution to the grand potential, it follows from (24) that an increasing vector meson condensate leads to an increasing overall pressure. Secondly, at the same time, the effective chemical potentials are lowered which reduces the contributions to the pressure and also the energy density via the particle densities from the (pseudo)scalar and quark sectors. Altogether, a larger vector coupling leads to an increase of the EoS's stiffness both for N f = 2 and N f = 2 + 1 quark flavors as expected.
Knowing now the EoS for pure quark matter, the mass-radius relation of a compact star can be obtained as solution of the Tolman-Oppenheimer-Volkoff equations, assuming a perfect fluid and a non-rotating star. Such pure quark stars [85] could exist under the hypothesis  of absolutely stable strange quark matter [86,87]. Within our setup, quark matter is not absolutely stable, but it is nevertheless instructive to investigate the mass-radius relation of quark stars with different strengths of the vector coupling. The results are summarized in Fig. 3, in the left panel for two flavor quark matter and in the right panel including strangeness. The colored horizontal bands indicate the measured two-solar-mass pulsars [2][3][4]. Increasing the vector coupling shifts the masses to larger values and increases the radii. Maximum masses are all compatible with observed pulsar masses, and radii are generally larger than current neutron star observations suggest [1,6,11,12], see also [56].
As mentioned above, within our setup quark matter is not absolutely stable and we will now turn to the construction of a hybrid matter EoS. Please note that while generally an increasing vector coupling increases the EoS's stiffness the pressure decreases with increasing vector coupling for a given quark chemical potential. This can be explained by the aforementioned reduction of the effective chemical potential in the quark loop, leading to an overall pressure reduction for a given chemical potential. This has significant consequences for the possible occurrence of hybrid stars, i.e. neutron stars with a quark matter core. We construct a hybrid EoS by assuming two separate phases, a nucleonic phase described by the HS(DD2) hadronic equation of state [88,89] and a quark matter phase described by the present FRG quark-meson EoS with additional vector meson interactions. The results are shown in Fig. 4, where the pressure as a function of the energy density is displayed (left panel for two quark flavors and right panel for N f = 2 + 1). Both phases are separated by a clear boundary and individually fulfill the weak equilibrium and charge neutrality conditions. A first-order transition is obtained via a Maxwell construction (horizontal dotted lines in the figure). It can be characterized by an onset energy density ε trans in the hadronic phase and a gap ∆ε given by the difference between the energy density in the quark phase at the end of the transition and the onset ε trans . The onset energy density also defines the transition pressure p HS(DD2),trans ≡ p HS(DD2) (ε trans ).
Due to the decreasing quark matter pressure at increasing vector interaction strength, the phase transition gradually moves to higher quark chemical potentials, i.e. to a higher intersection pressure and a higher ε trans . ∆ε also increases for increasing g v . For N f = 2+1,     EoS (gray color) and the quark matter phase by the FRG quark-meson EoS with vector mesons (see Fig. 2). Both phases separately meet weak equilibrium and charge neutrality and are connected via a Maxwell construction. For N f = 2 + 1, energy densities larger than ε ∼ 1200 MeV/fm 3 corrresponding to µ > 500 MeV are dropped.
the transition generally occurs at lower pressures and with larger ∆ε than for N f = 2 due to the additional strange degree of freedom in the quark matter EoS. The size of the discontinuity of the energy density ∆ε determines, too, the stability of the hybrid star against gravitational collapse: a large discontinuity destabilize the star immediately at the transition point p = p trans whereas for a small discontinuity a small quark core forms and the star remains stable. This scenario can be summarized in terms of the Seidov limit [90] ∆ε crit ε trans = 1 2 + 3 2 p trans ε trans . (28) ∆ε crit denotes here the threshold value below which a stable hybrid star branch is connected to a hadronic star branch. Thus, above the Seidov limit the sequence of stars become unstable immediately. All EoS displayed in Fig. 4 remain well below this limit. If ∆ε > ∆ε crit , a so-called "third family" [91] stable sequence of hybrid stars may exist at higher central densities for certain conditions, leading eventually to twin or even triplet configurations [92,25]. The conditions for the existence of such twin or triplet configurations has been discussed in detail in [25], characterising the transition by the two parameters ε trans and ∆ε together with a constant speed of sound parameterisation of the quark phase. In Fig. 5 we show, together with the FRG hybrid EoS with g v = 0, such a paremetrisation , ε < ε trans p HS(DD2) (ε trans ) , ε trans < ε < ε trans + ∆ε p HS(DD2) (ε trans ) + s [ε − (ε trans + ∆ε)] , ε > ε trans + ∆ε . (29) ε trans has thereby been chosen very close to the transition density in the FRG two-flavor hybrid EoS, whereas ∆ε/ε trans = 0.6 is close to the Seidov limit. It is obvious, as already anticipated from the Seidov limit, that the FRG hybrid stars remain stable. Two different values for the speed of sound in the quark phase have been chosen. The first one, s = 1/3, is close to the FRG sound speed, whereas the second one s = 1 represents the causality limit, i.e. the stiffest possible quark matter EoS. As already discussed in [56], for s = 1/3, the pressure in the quark matter phase is not sufficient to counteract the strong gravitational pull due to the large energy density of the quark core, cf. Ref. [25], and thus does not support a stable hybrid star branch. For s = 1, a stable third family branch with twin configurations exists [56]. Since even with the inclusion of vector interactions ∆ε does not exceed the Seidov limit and the sound speed is only insufficiently increased in the quark phase, the findings discussed here for the case g v = 0 remain valid for nonzero vector coupling. We thus confirm the conclusion of [56] that the occurrence of twin stars in our model is ruled out due to the small energy gap at the phase transition from nuclear matter to quark matter and due to the small stiffness of the quark matter EoS.
The shift of the phase transition in a hybrid star to higher densities with increasing vector coupling can be seen in the mass-radius relations, too, shown in Fig. 6 from the combined HS(DD2) and the present FRG quark-matter EoSs (left panel two quark flavors and right with strangeness). An increase of the vector coupling leads to a continuously smaller quark matter core in the hybrid star, but an increasing maximum mass. Especially for N f = 2 + 1 (right panel), where without vector interactions the two solar mass limit cannot be satisfied, the maximum mass is in agreement with current observations for g v 1. However, a quark matter core is only found as a small, continuous branch. For example, for g v = 1 the N f = 2 + 1 hybrid star model yields for the heaviest stable star a quark matter core with radius around 3.2 km, constituting about 4% of the star's total gravitational mass. For g v = 2, the heaviest star's quark core radius is 1.6 km and makes up only 0.6% of its total mass. Concerning the possibility of twin stars, the mass-radius relations confirm the absence of twin configurations within our model.
Another interesting quantity that is experimentally accessible is the tidal deformability of neutron stars. For a static, spherically symmetric star, placed in a static external quadrupolar tidal field E ij , the tidal deformability λ can be defined to linear order as where Q ij represents the star's induced quadrupole moment. The tidal parameter λ can then be computed from a perturbation of the spherical TOV solution, see [93] for more details.     The results for the dimensionless tidal deformability Λ = λ/M 5 as a function of the star's gravitational mass are shown in Fig. 7 for hybrid and quark stars with two and three quark flavors, respectively, and different values of the vector coupling. The pure quark stars lead to tidal deformabilities which are significantly too large compared with the GW170817 observations [7,6]. This is also in line with findings from parameterized EoS's in [94]. As discussed before, we do not expect pure quark stars to exist within our setup since quark matter is not absolutely stable. The hybrid star tidal deformabil-ities only differ from the HS(DD2) ones close to their respective maximum masses, i.e., the quark cores are too small to have an impact on the tidal deformabilities for all inspiral stars of masses below ≈ 1.8 M or even higher, depending on g v . The HS(DD2) tidal deformability is in slight tension with the GW170817 observations. A hybrid construction with other nucleonic EoS leading to lower tidal deformabilities could thus be appropriate. However, in [56] we found no intersection of the quark matter EoS with other such nucleonic EoS, since the nucleonic pressure over the entire relevant range exceeded the quark matter one for given chemical potentials. Since with increasing vector coupling, the pressure is reduced for given chemical potential, we confirm that we do not find hybrid stars with lower tidal deformabilities within the present FRG approach to the quark matter EoS.

Summary and conclusions
Based on a previous work [56] a two-and three flavor quark-meson model that fully incorporates chiral symmetry breaking has been augmented with vector mesons to investigate the quark matter EoS relevant for neutron star physics. Quantum and density fluctuations of the quarks and the (pseudo)scalar meson channels beyond the mean-field approximations are treated using the nonperturbative functional renormalization group method wherein the coupling of the isoscalar vector mesons to the FRG flow is taken into account.
As an application pure quark stars as well as hybrid stars are analyzed for different vector meson couplings always in β-equilibrated neutral matter. Since quark matter within our setup is not absolutely stable we construct hybrid stars where the EoS obtained with the FRG are combined with a nuclear EoS via a Maxwell construction. In general, an increase of the vector interaction increases the EoS's stiffness while the pressure decreases for a given chemical potential due to a reduction of the effective chemical potential in the quark loop.
Similar to the previous findings of the FRG EoS without vector interactions the inclusion of strangeness softens the EoS at high densities. This is in contrast to the addition of repulsive vector interactions to the system which stiffens the EoS. However, because of the decreasing quark matter pressure at increasing vector interaction strength the phase transition is shifted to higher quark chemical potentials. Also the energy density gap at the transition increases with increasing vector couplings. Including strangeness, the transition generally occurs at lower pressures but with larger energy density gaps which is plausible due to additional strange degree of freedom in the EoS.
As a consequence pure quark matter as well as hybrid stars with vector interactions lead to larger maximum masses such that the stars are consistent with existing observed pulsar masses. For pure quark stars, the radii are generally larger than observations suggest. Note that for hybrid stars without vector interactions the two solar mass limit cannot be reached when strangeness is included.
However, even though the hybrid stars' maximum mass increases with vector interaction, the quark matter core becomes smaller. Despite the fact that a quark matter core is only found as a small continuous branch in the mass-radius relation the largest core radius is found for a vanishing vector coupling. Within our setup the possibility of twin stars can be excluded.
Furthermore, hybrid stars with lower tidal deformabilities can also be excluded with the present FRG setup. The tidal deformability is mainly determined by the nucleonic EoS and differs from it only in the vicinity of its maximum mass where a quark matter core forms. The nucleonic pressure required to produce smaller tidal deformabilities would be too large to allow for a hybrid construction over the relevant range of chemical potentials.
In total, regarding repulsive vector interactions the occurrence of hybrid stars with extensive quark cores seems to be disfavored within our setup, in particular in view of experimental restrictions on the masses, radii, and tidal deformabilities. Besides these findings, a number of open issues remain: A systematic improvement of the employed FRG truncation, specifically including the running of higher derivative couplings in the (pseudo)scalar channel, might lead to further insights on the impact of fluctuations and on the robustness of the approach. Finally, the role of many-body quark correlations relevant for a more realistic description of dense matter -such as residual six-quark or diquark correlations in the quark phase -needs to be considered in the future.

A Input Parameters
In this appendix all input parameters for the FRG calculations are given. For further information see also [69,70,56]. We start with N f = 2+1 quark flavors and then simplify to N f = 2 quark flavor. The finite pseudoscalar masses of the pions, m π = 138 MeV, and kaons, m K = 496 MeV, are fixed via the explicit chiral symmetry breaking terms c l = (120.73 MeV) 3 and c s = (336.41 MeV) 3 . Without an explicit symmetry breaking all pseudoscalar masses would vanish in a chirally invariant or spontaneously broken theory due to the Goldstone theorem. The summed squares of the η and η masses m 2 η + m 2 η = (1103.2 MeV) 2 are reproduced with the axial U (1) A symmetry breaking parameter c A = 4807.84 MeV. The remaining three parameters in the ultraviolet effective potential are fixed with the (broad) sigma meson resonance mass which we have chosen to be of the order of m σ = 560 MeV and the two vacuum condensates σ l,0 = 92.4 MeV and σ s,0 = 94.5 MeV that yield the pion and kaon decay constants, f π = 92.4 MeV and f K = 113 MeV. The constituent light and strange quark masses follow from a single Yukawa coupling g = 6.5, i.e. m l = gσ l,0 /2 ≈ 300 MeV and m s = gσ s,0 / √ 2 ≈ 434 MeV. For N f = 2, we proceed in an analogous fashion. Only the remaining non-strange condensate σ l,0 and m σ are set by the two free parameters in the chiral potential in the ultraviolet whereas m l and m π are fixed by g and c l as before.
The used input parameters of the chiral potential for an initial UV cutoff of Λ = 1 GeV and an infrared cutoff of k IR = 80 MeV are summarized in Tab Table 1: FRG input parameters for the chiral ultraviolet potential, Eq. (15), for N f = 2 and N f = 2 + 1 quark flavors.

B Numerical Solution Details
In this work, one of the numerical challenges consists of solving for multiple gap equations (23) simultaneously, whereas for each trial point a full FRG flow equation (12) in the (pseudo)scalar sector has to be solved. In order to retain reasonable computation times, we solve the vector meson gap equations by computing data points on a discrete set of vector meson condensates and interpolating key quantities like the number densities (25). This allows us to use the same basic set of points for all coupling strengths. The same process is used for the electron chemical potential to satisfy charge neutrality, see (8). Instead of keeping the interpolated values for the equation of state, we then calculate new data points with the appropriate shifts in the chemical potentials inserted. This way, we have a method to check to what extent the gap equations and charge neutrality condition are actually fulfilled and to thereby gauge the quality of the interpolation. The solution of the FRG flow equation (12) for the (pseudo)scalars is performed numerically by discretizing field space [95,96]. This leads to a set of coupled ordinary differential equations. Considering recent numerical advances in the field, see [97], for N f = 2 flavors we employ an upwind finite difference scheme for the determination of field derivatives while in parallel solving the flow for the potential and its derivative. For N f = 2 + 1 flavors, the two-dimensional grid of the two scalar field variables is linked via cubic splines as outlined in [69] and also applied in preceding works [56]. While we found the latter method to lead to small numerical inaccuracies for large vector couplings g v at chemical potentials higher than the onset of strangeness around µ ≈ 430 MeV, the close agreement with the vastly different method for N f = 2 flavors at smaller chemical potentials gives us confidence in the general validity of the employed approach.