Functional renormalization group study of ρ meson condensate at a finite isospin chemical potential in the quark meson model

Quantum chromodynamics (QCD) describes the strong interactions between quarks and gluons. Its phase structure, dependent on temperature and chemical potential, reveals various phases of matter [1−6]. The phase structure is often depicted in the baryon chemical potential $ \mu_{B} $ versus temperature T plane, known as the $ \mu_{B} $-T plane [7−11]. This diagram illustrates the various phases of QCD matter under extreme conditions [3, 12−15], such as those in high-energy collisions or dense astrophysical objects like neutron stars.

The isospin chemical potential ($ \mu_{I} $) significantly influences the structure of the QCD phase diagram by shifting the boundaries of critical transitions such as chiral symmetry restoration and quark-gluon plasma (QGP) formation to higher values of baryon chemical potential ($ \mu_{B} $) [16−19]. This is due to the isospin asymmetry between up and down quarks [20]. When $ \mu_{I} $ exceeds the pion mass, a pion condensation phase occurs, resulting in spontaneous symmetry breaking [21−24]. Furthermore, large values of $ \mu_{I} $ alter the formation of the QGP and modify the behavior of color superconducting phases, especially in high-density environments like neutron stars [25]. In this connection, lattice (QCD) is free from the sign problems, thereby permitting the utilization of standard Monte Carlo theory to calculate thermodynamic properties and to draw the phase diagram as a function of temperature (T) and isospin ($ \mu_{I} $). This benefit allows confrontation with low-energy effective theories, such as chiral perturbation theory [26−28], as well as with models like the Nambu-Jona-Lasinio [29−38] and the quark-meson model [39−47].

In Ref. [48], holographic QCD was employed to investigate ρ-meson condensation, providing valuable insights from the strong-coupling regime that extend our understanding of this phenomenon beyond conventional approaches. In the present work, we address a related but simpler problem: ρ condensation in the presence of a finite isospin chemical potential $ \mu_{I} $. Using a two-flavor quark–meson model, we explore the formation and behavior of the ρ condensate under varying $ \mu_{I} $, temperature (T), and coupling constant ($ g_{\rho} $). Our objective is to advance the understanding of QCD phase structure and the properties of dense hadronic matter. Calculations within such models are typically performed at the mean-field level; going beyond this approximation is an essential step toward a more complete and realistic description. To achieve this, we employ the functional renormalization group (FRG) method [49−55], which systematically incorporates quantum and thermal fluctuations beyond mean-field theory. The Functional Renormalization Group is a potent non-perturbative method that allows quantum and thermal fluctuations to be incorporated into a field theory. The FRG has been utilized extensively to investigate the QCD phase diagram with chiral effective models beyond the MF, such as the NJL and the QM models.

This paper is organized as follows: The 2-flavour quark-meson model, including the ρ vector meson, and the FRG method, is presented in Sect. 2. In Sect. 3 the results are discussed. Finally, in Sect. 4 the conclusions is presented.

The Lagrangian of the two-flavour Quark-Meson model with ρ vector meson in Euclidean space is

The field strength tensors of the vector bosons $ {\boldsymbol{\rho_\mu}} $ is generally given as: $ {\boldsymbol{R}}^{(\rho)}_{\mu\nu}=\partial_\mu {\boldsymbol{\rho}}_\nu-\partial_\nu {\boldsymbol{\rho}}_\mu-g_\rho {\boldsymbol{\rho}}_\mu \times {\boldsymbol{\rho}}_\nu $. A field ψ is the light two-flavor quark field $ =(u,d)^T $, and coupled to scalar (σ) and isoscalar (π) fields transforming as a four-component field $ (\sigma,{\boldsymbol{\pi}})^T $ under the chiral group. A bold symbol stands for a vector, and $ {\boldsymbol{\tau}}=(\tau_1,\tau_2,\tau_3) $ are the Pauli matrices in isospin space, and introducing isospin chemical potential $ \mu_I=\mu_u-\mu_d $. Treated as background fields, only $ \rho^3_0 $ is non-vanishing, and the non-abelian of $ {\boldsymbol{R}}^{(\rho)}_{\mu\nu} $ does not contribute in practice. The spatial components vanish because of the assumption of homogeneous isotropic matter. A Hubbard-Stratonovich transformation bosonizes these interactions, introducing effective vector-isovector fields, $ {\boldsymbol{\rho}}_\mu $. While the fluctuations of the π and σ fields will be included non-perturbatively, $ {\boldsymbol{\rho}}_\mu $ will be treated as mean fields. These vector bosons conveniently parametrize unresolved short-distance physics. The potential for $ \sigma,{\boldsymbol{\pi}} $ and $ {\boldsymbol{\rho}}_{\mu} $ is

where $ f_\pi $ is the pion decay constant. We will use the value $ f_\pi=93 $ MeV, and $ m_\rho\sim 1 $ GeV. The parameters in our model are $ g_s $, $ g_\rho $ and λ. The values of these parameters for the FRG calculations are seeted to reproduce the same value for quantities such as the constituent quark mass of $ \sim 300 $ MeV. As for the value of $ g_\rho $, $ m_\rho $, in our calculation, they always appear in the form of $ g_\rho / m_\rho $, so we will not discuss their values independently.

Functional Renormalization Group (FRG), a powerful non-perturbative method that allows incorporating quantum and thermal fluctuations in a field theory [52, 56] and has been extensively applied to effective QCD models [47, 57, 58]. The effective average action $ \Gamma_k $ with a scale k obeys the exact functional flow equation

Where $ \Gamma^{(2)}_k $ is the second functional derivative of the effective average action with respect to the fields, the trace includes momentum integration as well as traces of overall inner indices. An infrared regulator $ R_k $ was introduced to suppress fluctuations at momenta below the scale k. Following the RG scale, the regulator may assume a functional form [59]. In this investigation, quarks serve as the dynamical fields in the flow equation, σ, and π, and they affect the effective potential and the $ \rho_{0}^3 $ field. Contrary to vector fields' spatial components, because it is not coupled to the time derivative, the $ \rho_{0} $ field is not dynamical. Therefore, the value of $ \rho_{0}^3 $ is completely fixed by specifying the values of other fields. At each scale k in the flow equation, we determine the value of $ \rho_{0}^3 $ field by solving the consistency equation for given σ, and π. These $ \rho^3_{0} $ field in turn appear in the effective chemical potential for quarks, affecting the dynamical fluctuations in the flow equations. Throughout our study, we neglected the flow of all wave-function renormalization factors in the so-called Local Potential Approximation (LPA)

The scale-dependent effective potential can be expressed by replacing the potential U with the scale-dependent one $ U_{k} $,

with the Euclidean Lagrangian from Eq.1, Finite temperatures are treated within the Matsubara formalism. The time-component is Wick-rotated, $ {t}\rightarrow -i\tau $ and the imaginary time τ is compactified on a circle with radius $ \beta=\dfrac{1}{T} $, where T is the temperature, for which after introduced $ \int d^4 x \equiv \int^{1/T}_0 d x_0 \int_V d^3 x $. Due to the chiral symmetry, the potential U depends on σ and π only through the chiral invariant

As mentioned, the vector field $ \rho^3_0 $ appears here only as mean fields. The complete k-dependence is in the effective potential $ U_k $. In analogy to the mean-field potential, the effective potential has a chirally symmetric piece, $ U_k^\phi $, the explicit chiral symmetry breaking term, and the mass terms of the vector bosons:

Starting with some ultraviolet (UV) potentials $ U_{\Lambda} $ as our initial conditions, we integrate fluctuations and obtain the scale-dependent $ U_k $. The form of $ U_k^\phi $ will be determined without assuming any specific forms, while for the potential of the ρ field:

To use Wetterich’s equation, a regulator function, that respects the interpolating limits of the effective average action, has to be chosen. We employ the so-called optimized or Litim regulator function [60], for bosons and fermions, respectively given by:

because of the structure of the regulators, the dependence on three-momenta is eliminated, and only integral over theta function remains. The flow equation for the potential $ U_k^\phi $ can be obtained as:

Here, $ \nu_{q} = 2\,({\rm{spin}}) \times 2\,({\rm{flavor}}) \times 3\,({\rm{color}}) = 12 $ is the quark degeneracy factor, and $ E_{q} = \sqrt{p^{2} + m_{{\rm{eff}}}^{2}} $ with the effective quark mass $ m_{{\rm{eff}}} = g_{s}\sigma $. The effective energies are given by:

for pion, sigma-meson, and quark, respectively. And the scale-dependent particle masses are:

and we also define $ U'_k=\dfrac{\partial U_k}{\partial \phi^2} $.

The effective chemical potential,

depends also on the field $ \rho_{0,k} $, and depend on the scale k. The extended occupation numbers simplify to the usual Fermi-Dirac distribution functions for boson and fermion occupation numbers:

According to the flow equation for the effective potential, the ρ field can be calculated self-consistently. Therefore, at each momentum scale k [61] we solve the equation for $ \rho^{3}_{0,k} $.

The dependence on $ \rho^{3}_{0,k} $ appears in the mass term and the fermion loop. The flow equation reads:

This equation constitutes our flow equations for the ρ field as functions of ϕ. Note that the flow equations for $ \rho^3_{0} $ can be solved for a given ϕ, independently of the potential $ U^{k}_{\phi} $ (which only tells us where the minimum of ϕ is). To understand the behavior of $ \rho^3_{0} $, for the moment we ignore the k dependence in $ \mu_{eff} $, and carry out the integration over k. Finally, the initial conditions for the flow equations must be set up. The UV scale Λ should be sufficiently large in order to take into account the relevant fluctuation effects and small enough to render the description in terms of the model degrees of freedom realistic [62]. In our calculation we follow the choice of Ref. [63], Λ = 500 MeV. The initial for the potential is

and set the parameters $ g_{s} $ = 3.2, λ = 8 with the vacuum effective potential from the FRG computation having the minimum at $ \sigma_{vac}\simeq 93 $ MeV, which is regarded as $ f_{\pi} $. We used the value of λ, which enforces ϕ to stay near $ f_{\pi} $, is $ (\lambda \sim 8) $. If we start with another initial condition with an additional $ \phi^{2} $ term to give the mass, we need to readjust λ but obtain qualitatively similar results; in fact, starting with the condition Eq.(21), the scale evolution first generates the $ \phi^{2} $ terms, reflecting the universality. The initial condition for the ρ field has not been examined in detail, and we simply try.

We also examined a different initial condition but do not present it here, as it does not affect our main results.

Assembling all these elements, we calculate the effective potential with the fluctuations integrated to $ \Gamma_{IR} $ = 0. The final step is to find Φ = $ \sigma^{*} $, which minimizes the effective potential. At the location of minimum, the effective potential is identified as the thermodynamic potential,

In practice, it is numerically expensive to reduce the IR cutoff, and we typically stop the integration under $ k_{IR} \simeq20 MeV $.

In this section, we investigate the effects of varying isospin chemical potentials and coupling constants on the chiral phase transition in the chiral limit. The analysis is performed using the FRG flow equation for the ρ meson, given in Eq. (20), over a range of temperatures and chemical potentials. Additionally, we determine the ρ meson condensate by solving the FRG effective potential equation, Eq. (20), for different parameter sets to facilitate comparison.

Fig. 1 illustrates the impact of the vector coupling constants on the chiral phase boundary. We observed that, at fixed temperatures, the chiral phase transition occurs at higher chemical potentials. Although ρ vector mesons contribute to the system's increased stability, their influence on the phase boundary becomes progressively limited as the coupling strength increases. The boundary of the first-order phase transition shifts to lower temperatures.

Figure 1.  (color online) The FRG $ {T-\mu} $ chiral phase diagram with different vector couplings. The solid lines show the first-order phase transition, and the dashed lines show the second-order phase transition. The stars show the TCPs. The parameters are set as: $ f_{\pi} $ = 93 MeV, $ g_{s} $ = 3.2, λ = 8, the ultraviolet cutoff $ \Lambda_{FRG} $ =500 MeV, $ {\mu_I} $ = 200 MeV.

The chiral phase diagram derived from the functional renormalization group equation at various isospin chemical potentials is presented in Fig. 2. As the isospin chemical potential increases, the chemical potential required for the chiral phase transition decreases, leading to a leftward shift of the critical endpoint toward lower temperature and chemical potential regions, where RG flow evolution is significantly influenced by low-temperature fluctuations. This shift of the chiral boundary with increasing $ \mu_I $ is consistent with both analytical calculations and lattice QCD results [22]. While Ref. [64]reported comparable behavior under specific conditions, the present study extends these findings and provides a more comprehensive demonstration.

Figure 2.  (color online) The FRG $ {T-\mu} $ chiral phase diagram with different $ \mu_{I} $. The solid lines represent the first-order phase transition, and the dashed lines represent the second-order phase transition. The stars show the TCPs. The parameters are set as: $ f_{\pi} $ = 93 MeV, $ g_{s} $ = 3.2, λ = 8, the ultraviolet cutoff $ \Lambda_{FRG} $ =500 MeV, the coupling constant $ {g_{\rho}} $ $ m^{-1}_{\rho} = 0.006 \ {\rm{MeV}}^{-1} $.

In order to understand the precise process of phase transition under the influence of isospin chemical potential, see Fig. 3. The change of order parameters in the low-temperature region of the phase diagram with chemical potential can be observed. Evidently, a first-order phase transition is occurring in this region. As the chemical potential increases, the chiral condensate changes and drops to zero at the phase transition. On the other hand, we note that with the increase of the isospin chemical potential, the initial value of the chiral condensate gradually decreases smoothly. It can be imagined that when the isospin chemical potential gradually increases to a certain value, the first-order phase transition region will eventually disappear, which is also expressed in the phase diagram.

Figure 3.  (color online) Chiral condensates as a function of quark chemical potential under different isospin chemical potentials, calculated at $ T=10 $ MeV, $ {g_{\rho}}/m_{\rho} = 0.006\; {\rm{MeV}}^{-1} $. The different colored lines correspond to different isospin chemical potentials.

Based on the analyses in Figs. 1 and 2, we selected optimal parameter values to compute the ρ meson condensation, considering different couplings $ {g}_{\rho}/m_{\rho} $. These choices were guided by the observed shifts in the phase boundary and the influence of varying isospin chemical potentials on the system’s stability. The data for $ {g}_{\rho} \rho^{3}_{0} $ were obtained by directly solving Eq. 20. The figures illustrate how combinations of temperature, chemical potential, and coupling strength affect meson condensation. By carefully tuning these variables, we accurately determined the condensation values, providing deeper insights into chiral phase transitions under varying conditions. Near chiral restoration, vector and axial-vector modes broaden and nearly degenerate. The ρ channel competes with the chiral condensate and reshapes $ U_k $, producing the mutual "pull" observed in the curves. This reflects a coupled dynamical interplay rather than two separable phenomena, arising from their shared origin in the QCD Lagrangian, where both structures emerge from the same fermionic degrees of freedom and their interactions. These results highlight the role of the ρ meson in stabilizing the system, particularly in regions of the phase diagram where fluctuations are pronounced, and are essential for understanding the behavior of hadronic matter in extreme environments such as heavy-ion collisions or astrophysical settings [65, 66].

In Fig. 4, we compare the results for the ρ condensate as a function of isospin chemical potential, for μ = 100 MeV, and different ρ couplings panel (a) $ {g_{\rho}} $ $ m^{-1}_{\rho} $= 0.002 MeV⁻¹, panel (b) $ {g_{\rho}} $ $ m^{-1}_{\rho} $= 0.006 MeV⁻¹, panel (c) $ {g_{\rho}} $ $ m^{-1}_{\rho} $= 0.008 MeV⁻¹, and panel (d) $ {g_{\rho}} $ $ m^{-1}_{\rho} $= 0.010 MeV⁻¹, respectively. The temperature for this calculation is T = 10 MeV. In this scenario, the chiral condensate decreases with increasing $ \mu_{I} $ while the ρ condensate starts to grow for $ \mu_{I} $ greater than the critical value. But at even higher isospin chemical potential, a new first-order transition occurs, and the chiral condensate drops to zero. In this new region, the ρ condensate becomes dominant when the coupling value for the ρ meson increases.

Figure 4.  (color online) The chiral and ρ condensates as a function of isospin chemical potential, shown for (a) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.002 $ MeV⁻¹, (b) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.006 $ MeV⁻¹, (c) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.008 $ MeV⁻¹ and (d) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.010 $ MeV⁻¹. The chemical potential is $ \mu = 100 $ MeV, temperature is $ T = 10 $ MeV and the ultraviolet cutoff $ \Lambda_{FRG} $ =500 MeV

Fig. 5 depict the chiral and ρ condensates as functions of the isospin chemical potential for $ \mu = 200 $ MeV, respectively. Across all panels, the behavior differs from that at lower μ. The ρ condensate begins to grow when $ \mu_I $ exceeds the critical value for a second-order phase transition, while the chiral condensate simultaneously drops to zero. In this regime, the ρ condensate remains dominant for larger ρ meson couplings. The chiral condensate has vanished in all panels of Fig. 6, corresponding to $ T = 10\; {\rm{MeV}} $ and $ \mu = 300\; {\rm{MeV}} $. Comparing $ g_\rho / m_\rho = [0.002,\, 0.010]\; {\rm{MeV}}^{-1} $, one finds that larger $ g_\rho $ values enhance isovector repulsion, causing the ρ condensate to dominate the order-parameter dynamics while the chiral condensate becomes almost completely suppressed. At strong coupling, the system is increasingly governed by the vector channel, reinforcing the conclusion that isovector repulsion reshapes the low-temperature phase structure [57]. In Figs. 4, 5 and 6, we observe the influence of the chemical potential on the chiral and ρ condensates; the isospin chemical potential takes different critical values when the transition of the ρ condensate occurred from zero to non-zero.

Figure 5.  (color online) The chiral and ρ condensates ρ as a function of isospin chemical potential, shown for (a) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.002 $ MeV⁻¹, (b) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.006 $ MeV⁻¹, (c) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.008 $ MeV⁻¹ and (d) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.010 $ MeV⁻¹. The chemical potential is $ \mu = 200 $ MeV, temperature is $ T = 10 $ MeV and the ultraviolet cutoff $ \Lambda_{FRG} $ =500 MeV

Figure 6.  The chiral and ρ condensates ρ as a function of isospin chemical potential, shown for (a) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.002 $ MeV⁻¹, (b) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.006 $ MeV⁻¹, (c) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.008 $ MeV⁻¹ and (d) $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.010 $ MeV⁻¹. The chemical potential is $ \mu = 300 $ MeV, temperature is $ T = 10 $ MeV and the ultraviolet cutoff $ \Lambda_{FRG} $ =500 MeV

Fig. 7 illustrates the evolution of the ρ condensate from microscopic to macroscopic scales with different isospin chemical potentials. In the mean-field approximation, the condensation of the ρ meson occurs only when the isospin chemical potential exceeds its vacuum mass, i.e., $ \mu_I > m_\rho $. However, beyond the mean-field level, such as in the Random Phase Approximation (RPA) and Chiral Perturbation Theory (CPT) [22], fluctuation effects significantly lower this critical chemical potential to $ \mu_I = m_\pi $. This value agrees well with our own calculations in Figs 4, 5 and 6. As illustrated in Fig. 7, our results further confirm that when fluctuations are included, ρ condensation appears at the IR scale.

Figure 7.  (color online) ρ condensate as a function of the field ϕ at IR cutoff = 20 MeV and the ultraviolet cutoff $ \Lambda_{FRG} $ =500 MeV, shown for (a) $ \mu = 100 $ MeV, (b) $ \mu = 200 $ MeV, (c) $ \mu = 240 $ MeV, and (d) $ \mu = 300 $ MeV. At the coupling $ \dfrac{{g_{\rho}}}{m_{\rho}} = 0.006 $ MeV⁻¹ and temperature is $ T = 10 $ MeV.

In Fig. 8, we show the ρ condensate as a function of the chemical potential for different isospin chemical potentials $ \mu_{I} = 200 $ MeV, $ \mu_{I} = 300 $ MeV,$ \mu_{I}= 400 $ MeV and $ \mu_{I}= 500 $ MeV, panel (a, b, c), and (d), respectively. It can be read that when the chemical potential increases, the ρ condensate increases first and then drops to zero. Also the critical chemical potential shifts left when the isospin chemical potential increases at a fixed temperature. Furthermore, it is seen that the ρ condensate increases with increasing couplings at a fixed temperature. It is worthy of mentioning that ρ meson becomes condensate when the isospin chemical potential is larger then 210 MeV for fixed chemical potential $ \mu = 200 $ MeV, see the panel (a) in Fig 5. Also it's more clear at isospin chemical potential $ \mu_{I} $ = 400 MeV panel(c) when the chemical potential is below than μ = 100 MeV, which is almost the critical chemical potential for the ρ meson condensate (in addition to the chiral condensate) see Fig 4. This happened for all the isospin chemical potential relation with chemical potential for the chiral and ρ condensates.

Figure 8.  (color online) ρ meson condensate as a function of chemical potential for different values of isospin chemical potential: (a) $ \mu_{I} = 200 $ MeV, (b) $ \mu_{I} = 300 $ MeV, (c) $ \mu_{I} = 400 $ MeV, and (d) $ \mu_{I} = 500 $ MeV. Here, $ T = 10 $ MeV.

Finally, Fig. 9 shows the ρ meson condensate as a function of temperature for different isospin chemical potentials: $ \mu_{I} = 200 $ MeV (panel (a)), $ \mu_{I} = 400 $ MeV (panel (b)), and $ \mu_{I} = 500 $ MeV (panel (c)), at a fixed $ \mu = 200 $ MeV and $ {g}_{\rho}/m_{\rho} = 0.006\; {\rm{MeV}}^{-1} $. The condensation value of the ρ meson for $ \mu_{I} = 200 $ MeV (panel (a)) and $ \mu = 200 $ MeV is consistent with the results obtained in Fig. 5, confirming the reliability of our calculations. The coupling dependence shows that a larger $ g_\rho $ enhances the magnitude of the ρ condensate without altering its qualitative temperature dependence. At finite temperature, the isovector order tends to melt, whereas stronger repulsion reinforces it at low T. This interplay explains why ρ effects are most pronounced in cold, dense matter. As the temperature increases, the ρ mode softens and approaches ρ–$ a_1 $ degeneracy, with the decreasing condensate following the chiral trend—signaling proximity to (partial) chiral restoration [41].

Figure 9.  (color online) ρ meson condensate as a function of temperature for different values of isospin chemical potential: (a) $ \mu_{I} = 200 $ MeV, (b) $ \mu_{I} = 300 $ MeV, (c) $ \mu_{I} = 400 $ MeV, and (d) $ \mu_{I} = 500 $ MeV. Here, $ \mu = 200 $ MeV.

We investigated the ρ meson condensation in isospin chemical potential by applying the FRG, using a two-flavor quark-meson model with the ρ meson. We have also investigated the impact of vector mesons and isospin chemical potential on the phase structure of the chiral phase transition. The primary conclusions are categorized into two facets: the impact of vector couplings and the isospin chemical potential on the phase structure. The phase boundary moves as a unit to the low temperature and low density region as the isospin chemical potential increases, with the vector coupling strength remaining constant. Consequently, the temperature of TCP gradually decreases. In contrast to changing only the vector coupling, as in [16], the isospin chemical potential also very slightly reduces the temperature at which the phase transition occurs at low chemical potential, and this effect is very slight, with no observed back bending behavior, consistent with the explanation provided in [64]. On the other hand, we investigated the ρ condensation as a function of the isospin chemical potential with different vector coupling constants. Beyond the mean-field approximation, fluctuation effects lower the critical isospin chemical potential for ρ meson condensation from $ \mu_I > m_\rho $ to $ \mu_I = m_\pi $, at which point the ρ meson condenses alongside the chiral condensate, consistent with the results from RPA and CPT [22], and confirmed by our FRG calculations. We noted that at large chemical potential and large isospin chemical potential, the ρ meson condensate dominates. Increasing the coupling constant for the ρ meson enhances the condensate value, though the critical isospin chemical potential remains relatively stable around 200 MeV. Increasing ρ meson coupling slightly shifts the boundary of the phase transition.

We thank Hai-cang Ren and Moran Jia for useful discussions.