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 temperatureTplane, known as the $ \mu_{B} $ -Tplane [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 scalekobeys 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 scalek. 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 scalekin 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 potentialUdepends 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 scalek. 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 scalek[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 thekdependence in $ \mu_{eff} $ , and carry out the integration overk. 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 theIRcutoff, 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. 1illustrates 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 inFig. 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, seeFig. 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 inFigs. 1and2, 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].

InFig. 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. 5depict 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 ofFig. 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]. InFigs. 4,5and6, 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. 7illustrates 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 inFigs 4,5and6. As illustrated inFig. 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.

InFig. 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) inFig 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) seeFig 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. 9shows 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 inFig. 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 lowT. 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.