Critical parameters of liquid-gas phase transition in thermal symmetric and asymmetric nuclear matter

TL;DR

This study investigates the critical parameters and phase diagram structure of liquid-gas phase transition in thermal symmetric and asymmetric nuclear matter using covariant density functional theory, highlighting the role of symmetry energy and isospin asymmetry.

Contribution

It explores correlations between critical parameters and symmetry energy, revealing how isospin asymmetry influences the phase transition properties in nuclear matter.

Findings

Critical temperature correlates with isospin asymmetry and symmetry energy slope.

Symmetry energy significantly affects the size of the liquid-gas coexistence region.

Phase diagram patterns are determined by critical temperature at non-zero isospin asymmetry.

Abstract

The properties of critical parameters and phase diagram structure of liquid-gas phase transition are investigated in thermal symmetric and asymmetric nuclear matter with the covariant density functional (CDF) theory. Although uncertainty remains in predicting the critical parameters such as the critical temperature and pressure from various CDF functionals, several correlations are explored numerically and verified to be approximately linear between them. These correlations become worse when nuclear matter is more isospin asymmetric, resulting mainly from the effects induced by symmetry energy. By looking over the isospin dependence of the critical temperature, the role of the symmetry energy in LG transition properties of asymmetric matter is realized. The change of critical temperature with isospin asymmetry is found to be correlated well with and as a consequence could be constrained by the density slope of symmetry energy at saturation density. Then, the structure of phase diagram of thermal nuclear matter is analyzed carefully. It is revealed that the contribution from symmetry energy dominates the size of liquid-gas phase coexistence area. Moreover, the specific pattern of the phase diagram could be determined by the critical temperature at non-zero isospin asymmetry, illustrated from the correlations of the temperature with pressures at several characteristic points, paving the possible way to further explore the structure of liquid-gas phase diagram of thermal nuclear matter.

Tables

Table 1. Table 1: Critical parameters of LG phase transition for symmetric nuclear matter, i.e., the critical temperature T C subscript 𝑇 𝐶 T_{C} (MeV), critical density ρ C subscript 𝜌 𝐶 \rho_{C} (fm -3 ), critical pressure P C subscript 𝑃 𝐶 P_{C} (MeV ⋅ ⋅ \cdot fm -3 ), compressibility factor Z C ≡ P C / ( ρ C ​ T C ) subscript 𝑍 𝐶 subscript 𝑃 𝐶 subscript 𝜌 𝐶 subscript 𝑇 𝐶 Z_{C}\equiv P_{C}/(\rho_{C}T_{C}) , and critical incompressibility K C subscript 𝐾 𝐶 K_{C} (MeV), as well as the incompressibility K 0 subscript 𝐾 0 K_{0} (MeV) at saturation density of symmetric nuclear matter at zero temperature.

	$T_C$	${\rho }_C$	$P_C$	$Z_C$	$K_C$	$K_0$
PKO1	14.53	0.045	0.191	0.286	-77.29	250.24
PKO2	15.76	0.042	0.220	0.332	-93.83	249.60
PKO3	14.57	0.048	0.198	0.283	-75.03	262.47
PKDD	14.91	0.049	0.225	0.308	-81.67	262.18
NL3	14.60	0.046	0.200	0.297	-77.90	271.73
PK1	15.11	0.049	0.223	0.305	-82.83	282.68

Table 2. Table 2: Similar as Table 1 , but for the asymmetric nuclear matter at isospin asymmetry δ = 0.5 𝛿 0.5 \delta=0.5 .

	PKO1	PKO2	PKO3	PKDD	NL3	PK1
$T_C^{\delta }$	12.56	14.13	12.85	12.79	12.31	12.78
${\rho }_C^{\delta }$	0.048	0.045	0.051	0.054	0.051	0.053
$P_C^{\delta }$	0.198	0.239	0.221	0.252	0.208	0.233
$Z_C^{\delta }$	0.328	0.376	0.337	0.365	0.331	0.344

Table 3. Table 3: The pressure of EC point P EC subscript 𝑃 EC P_{\rm{EC}} and its difference P CP − P EC subscript 𝑃 CP subscript 𝑃 EC P_{\rm{CP}}-P_{\rm{EC}} from the critical pressure for LG phase diagram at temperature T = T C δ = 0.5 𝑇 superscript subscript 𝑇 𝐶 𝛿 0.5 T=T_{C}^{\delta=0.5} , taken from Fig. 7 . Correspondingly the components in P CP − P EC subscript 𝑃 CP subscript 𝑃 EC P_{\rm{CP}}-P_{\rm{EC}} are given according to Eq. ( 31 ). The values are in unit of MeV ⋅ ⋅ \cdot fm -3 .

	$P_{\mathrm{EC}}$	$P_{\mathrm{CP}}-P_{\mathrm{EC}}$	$\Delta P_{E_0}$	$\Delta P_S$	$\Delta P_{E_0}+\Delta P_S$	$\Delta P_{E_S}$
PKO1	0.112	0.086	-0.442	0.407	-0.035	0.121
PKO2	0.154	0.085	-0.403	0.377	-0.026	0.113
PKO3	0.124	0.097	-0.469	0.433	-0.036	0.135
PKDD	0.125	0.127	-0.589	0.533	-0.056	0.180
NL3	0.106	0.102	-0.530	0.480	-0.050	0.150
PK1	0.120	0.113	-0.571	0.517	-0.054	0.168

Table 4. Table 4: The Pearson’s linear correlation coefficients r G subscript 𝑟 𝐺 r_{G} ( r L subscript 𝑟 𝐿 r_{L} ) between T C δ = 0.5 superscript subscript 𝑇 𝐶 𝛿 0.5 T_{C}^{\delta=0.5} and P Gas subscript 𝑃 Gas P_{\mathrm{Gas}} ( P Liquid subscript 𝑃 Liquid P_{\mathrm{Liquid}} ) at various isospin asymmetry δ 𝛿 \delta shown in Fig. 7 , see text for details.

$\delta$	$0$	$0.1$	$0.2$	$0.3$	$0.4$	$0.5$
$r_G$	0.978	0.978	0.977	0.974	0.969	0.957
$r_L$	0.978	0.976	0.964	0.933	0.895	0.880

Equations

$$E\equiv$$

$$\Omega =$$

$${\rho }_b=$$

$$S=$$

$$\displaystyle\big{[}\bm{\gamma}\cdot\bm{p}^{*}+M^{*}\big{]}$$

$$n_{\tau }(p)=$$

$$E^{\mathrm{kin}}=$$

$$E_{\sigma }^D=$$

$$E_{\varphi }^E=$$

$$\displaystyle\hskip 10.00002pt\times\Big{[}A_{\phi}+\hat{M}_{\tau}(p)\hat{M}_{\tau^{\prime}}(p^{\prime})B_{\phi}+\hat{P}_{\tau}(p)\hat{P}_{\tau^{\prime}}(p^{\prime})C_{\phi}\Big{]},$$

$${\rho }_s=$$

$${\rho }_b=$$

$${\rho }_{b3}=$$

$$E_{N.L.}=-\frac{1}{6}{g}_2{\sigma }^3-\frac{1}{4}{g}_3{\sigma }^4+\frac{1}{4}{c}_3{\omega }_0^4.$$

$$\Sigma (p)u(p,s,\tau )=\frac{\delta }{\delta \bar{u}(p,s,\tau )}\sum _{\varphi }\left[E_{\varphi }^D+E_{\varphi }^E\right].$$

$$P=$$

$$P=P_{E_0}+P_{E_S}+P_S$$

$$P_{E_S}={\rho }_b^2\frac{\partial }{\partial {\rho }_b}\left[{\delta }^2{E}_S({\rho }_b)\right].$$

$$\displaystyle L=3\rho_{b}\frac{\partial E_{S}(\rho_{b})}{\partial\rho_{b}}\Big{|}_{\rho_{b}=\rho_{0}},$$

$${\mu }_{\tau }^L(T,{\rho }_b^L,{\delta }^L)=$$

$$P^L(T,{\rho }_b^L,{\delta }^L)=$$

$${\rho }_b{\left(\frac{\partial P}{\partial {\rho }_b}\right)}_{T,\delta }>$$

$$\displaystyle\frac{\partial P}{\partial\rho_{b}}\Bigg{|}_{T=T_{C}}=$$

$$\displaystyle\frac{\partial\mu_{\tau}}{\partial\delta}\Bigg{|}_{T=T_{C}}=$$

$$\displaystyle K_{C}=9\rho_{b}^{2}\frac{\partial^{2}}{\partial\rho_{b}^{2}}\frac{F}{\rho_{b}}\Big{|}_{\rho_{C}}.$$

$$\displaystyle\frac{\partial}{\partial\rho_{b}}\left(\rho_{b}^{2}\frac{\partial}{\partial\rho_{b}}\frac{F}{\rho_{b}}\right)\Bigg{|}_{\rho_{C}}=0,$$

$$K_C+18\frac{P_C}{{\rho }_C}=0,$$

$$T_C=a{P}_C+b,$$

$$Z_C=\frac{P_C}{{\rho }_C{T}_C}$$

$$T_C^{\delta }=c(P_C^{\delta }-P_{C,E_S}^{\delta })+d,$$

$$\Delta T_C^{\delta }\equiv T_C^{\delta =0}-T_C^{\delta }.$$

$$\Delta T_C^{\delta }=c{P}_{C,E_S}^{\delta }+\left[a{P}_C^{\delta =0}-c{P}_C^{\delta }+b-d\right].$$

$$P_{\mathrm{CP}}-P_{\mathrm{EC}}=\Delta P_{E_0}+\Delta P_{E_S}+\Delta P_S,$$

Full text

I Introduction

Features of nuclear matter at finite temperature are of fundamental importance in nuclear physics. Among them, the so-called liquid-gas (LG) phase transition in sub-saturated thermal nuclear matter has drew a lot of attention, which occurs due to the van der Waals behavior of the nucleon-nucleon interaction Siemens (1983). The LG phase transition of thermal nuclear matter have been studied both experimentally and theoretically in a variety of works over the past several decadesFinn et al. (1982); Siemens (1983); Panagiotou et al. (1984); Jaqaman et al. (1984); Kapusta (1984); Su et al. (1987); Müller and Serot (1995); Baldo and Ferreira (1999); Ma (1999); Natowitz et al. (2002); Gobet et al. (2002); Karnaukhov et al. (2003); Chomaz et al. (2004); Das et al. (2005); Pichon et al. (2006); Brown et al. (2007); Li et al. (2008); Elliott et al. (2013); Borderie and Frankland (2019), and its important impact has been illustrated on many aspects of nuclear physics, such as heavy ion collisionsChomaz et al. (2004); Das et al. (2005); Brown et al. (2007); Li et al. (2008) and nuclear astrophysicsPethick (1992); Prakash et al. (1997); Lattimer and Prakash (2004, 2016).

The occurrence of LG phase transitions has been confirmed in both symmetric and asymmetric nuclear matter, being recognized from the survey of nuclear caloric curve and multifragment distribution in heavy-ion collisions Finn et al. (1982); Panagiotou et al. (1984); Natowitz et al. (2002); Gobet et al. (2002); Le Fèvre et al. (2005); Sfienti et al. (2009); Elliott et al. (2013). In the works for symmetric nuclear matter, one usually concerned the critical temperature $T_{C}$ of LG phase transition as an important and characteristic quantity Su et al. (1987); Baldo and Ferreira (1999); Rios et al. (2008); Rios (2010); Lourenço et al. (2016). In general, $T_{C}$ is predicted with a large uncertainty, around $10\thicksim 20$ MeV, from several theoretical models of thermal nuclear matterSu et al. (1987); Müller and Serot (1995); Zuo et al. (2004); Rios et al. (2008); Xu et al. (2007); Rios (2010); Vovchenko et al. (2015); Vovchenko (2017). On the other hand, its precise value is hardly constrained experimentally as well Panagiotou et al. (1984); Natowitz et al. (2002); Karnaukhov et al. (2003); Elliott et al. (2013). First, the uncertainty comes from the limitation that the available experiments only extrapolate $T_{C}$ to infinite matter from fragmentation reactions on finite nuclei Natowitz et al. (2002); Das et al. (2005); Sfienti et al. (2009). Moreover, the finite-size effects are also found to be influential in estimating the critical parameters for experiments, where a model dependence is involved inevitably Moretto et al. (2003, 2005); Elliott et al. (2013). The knowledge about the dependence of $T_{C}$ on various bulk properties of nuclear matter then would make sense on eliminating such an uncertainty. It was revealed that $T_{C}$ could be correlated with the incompressibility of symmetric nuclear matter at zero temperature and the nucleon effective mass at nuclear saturation density $\rho_{0}$, paving the way to deduce $T_{C}$ by constrain related quantities Lourenço et al. (2016). Recently, one also concerned the correlation among critical parameters themselves, i.e., critical temperature $T_{C}$, critical density $\rho_{C}$ and critical pressure $P_{C}$, to understand the behaviors of LG phase transition Rios (2010); Elliott et al. (2013).

Recently, a lot of experimental efforts have been made on nuclear LG phase transition with extreme interest in its isospin asymmetry dependence Le Fèvre et al. (2005); Sfienti et al. (2009); McIntosh et al. (2013a, b). With various kinetic thermometer approaches, the dependence of the nuclear caloric curve on the neutron-proton asymmetry can be extracted, correspondingly providing experimental information on the limiting temperatures of finite nuclei which is correlated with the critical temperature of nuclear matter Natowitz et al. (2002). For asymmetric nuclear matter, various theoretical works predict a phase diagram structure of LG phase transition in a wide range of isospin asymmetry and pressure Müller and Serot (1995); Xu et al. (2007); Dutra et al. (2008); Sharma and Pal (2010); Zhang and Jiang (2013); Fedoseew and Lenske (2015). Some of them argued that the phase diagram properties of LG phase transition could be correlated and affected by the bulk properties of nuclear matter, especially the symmetry energy Müller and Serot (1995); Xu et al. (2007); Sharma and Pal (2010); Zhang and Jiang (2013). However, the quantitative evidence is still not ample to deduce a clear correlation among each characteristic quantities of the LG phase diagram.

To investigate the LG phase transition of thermal nuclear matter, the covariant density functional (CDF) theory Serot and Walecka (1986); Reinhard (1989); Ring (1996); Bender and Heenen (2003); Vretenar et al. (2005); Meng et al. (2006); Nikšić et al. (2011); Meng and Zhou (2015); Meng2016 has been extended to the case of finite temperature, with many important achievements in describing the EoS of thermal nuclear matterWaldhauser et al. (1987, 1988); Zimanyi and Moszkowski (1990); Delfino et al. (1995); Avancini et al. (2004); Typel et al. (2010); Shen et al. (2010), the physics of supernova and proto-neutron starShen et al. (1995, 1998); Broderick et al. (2000); Menezes and Providência (2003); Avancini et al. (2008, 2010); Shen et al. (2011); Hempel et al. (2012), and properties of excited hot nucleiGambhir et al. (2000); Niu et al. (2009, 2011, 2013a); Liu et al. (2015); Zhang and Niu (2017, 2018), etc.. The critical temperatures of LG phase transition within the CDF calculations, in general locating around $T_{C}\thickapprox 15$ MeV, still keep model dependence Müller and Serot (1995); Sharma and Pal (2010); Zhang and Jiang (2013); Fedoseew and Lenske (2015).

In recent years, the CDF approach with Fock terms, namely the relativistic Hartree-Fock (RHF) theory, was also developed in terms of the density dependent meson-nucleon couplingLong et al. (2006, 2007, 2010). Significant improvements were obtained by the RHF theory, with the involvement of the exchange diagram and the self-consistent tensor force effects Jiang et al. (2015a); Zong and Sun (2018); Wang et al. (2018), in describing not only the properties of nuclear ground stateLong et al. (2009); Lu et al. (2013); Wang et al. (2013); Li et al. (2014, 2016) but the excitation and decay modesLiang et al. (2008); Niu et al. (2013b). Besides, several topics on the isospin properties of nuclear matter were studied as well, demonstrating the important role of Fock terms in the nuclear symmetry energy and neutron star properties Sun et al. (2008); Long et al. (2012); Jiang et al. (2015b); Zhao et al. (2015); Liu et al. (2018). As an application of the RHF theory to hot nuclei, the pairing transition in both stable and weakly bound nuclei has already been studiedLi et al. (2015). A series of novel phenomenon could occur when the contribution of continuum states is dressed by a finite temperature. For instance, a pairing re-entrance is predicted for drip-line nucleus 48Ni Belabbas et al. (2017). However, further systematical study still need to verify the robustness of these predictions. Alternatively, it is interesting to investigate the properties of thermal nuclear matter such as the LG phase transition within these newly developed CDF approaches, especially their model dependence on the selection of effective interactions.

In this work, based on the finite temperature CDF theory with and without Fock terms, the properties of liquid-gas phase transition in thermal symmetric and asymmetric nuclear matter will be studied. The critical parameters of LG phase transition and their correlations with several bulk properties of nuclear matter will be analyzed in detail. The paper is organized as follows. The formalism of the CDF theory for thermal nuclear matter is briefly introduced in Section II. In Sec. III we present the results within CDF calculations and discussion, including the critical point properties of LG phase transition in nuclear matter in Sec. IIIA, the properties of LG phase diagrams in Sec. IIIB, and the correlations between the critical temperature and characteristic pressures in LG phase diagram in Sec. IIIC. Finally, a short summary is given in Sec. IV.

II Thermal nuclear matter under the CDF theory

In this section, the general formalism of the covariant density functional theory for thermal nuclear matter will be described briefly. In order to eliminate the model dependence of the analysis as soon as possible, we will utilize three different meson-exchange types of CDF theory, namely, the relativistic mean-field approach with the nonlinear self-coupling of mesons (denoted as NLRMF), the density-dependent relativistic mean-field (DDRMF) and the relativistic Hartree-Fock (DDRHF) approaches. The corresponding formalism at zero temperature has already been addressed in several referencesBan et al. (2004); Sun et al. (2008).

Based on the meson exchange diagrams of nuclear force, the theoretical starting point — Lagrangian density can be deduced associated with the degrees of freedom of nucleons ($\psi$), two isoscalar mesons ($\sigma$ and $\omega$), two isovector mesons ($\pi$ and $\rho$), and photons ($A$). For uniform nuclear matter systems, the photon field, describing the electro-magnetic interactions between protons, is ignored naturally. Following the standard procedure Bouyssy et al. (1987), the energy functional is then obtained by taking the expectation of the Hamiltonian operator $\mathcal{H}$ with respect to the ground state $|\Phi_{0}\rangle$,

[TABLE]

where $E^{\rm kin}$ and $E_{\phi}$ denote the kinetic and potential energy densities, respectively, and for the latter the Hartree-Fock approach leads to two types of contributions: the direct (Hartree) $E_{\phi}^{D}$ and exchange (Fock) terms $E_{\phi}^{E}$. According to the specific CDF functional $\phi$ could be $\sigma,\omega,\rho,\pi$ etc.. The further details can be found in Ref. Sun et al. (2008).

The CDF theory at finite temperature is then deduced by considering the grand canonical ensemble in quantum statistical mechanics, where the thermal equilibrium state for a statistical $N$-body system can be determined by the variation of grand canonical potential $\Omega$,

[TABLE]

here $F$, $E$, $S$ and $T$ are the free energy, the total energy, the entropy and the temperature, respectively. The associated Lagrange multiplier $\mu$, also referred as chemical potential, is introduced to preserve the particle number at average. Different from the standard CDF approach, the thermal excitation will lead to the spreading of valence particles over the states around the Fermi surface such that the occupation probability $n_{i}$ of the state $i$ is not 1 or 0 any more. Therefore, the nucleon density and particle number $N$ read as,

[TABLE]

where $u_{i}$ is the Dirac spinor for state $i=(p,s,\tau)$, which satisfies the normalization condition $u^{\dagger}_{i}u_{i}=1$. Correspondingly, the entropy $S$ is,

[TABLE]

In this work, the finite-size effects, for instance discussed in Refs. Qian et al. (2001); Pawłowski (2002), will not be considered for simplicity since the motivation here focuses mainly on systematical exploration of correlations among critical parameters of LG phase transition based on a series of CDF functionals.

The variation of grand canonical potential $\Omega$ shall be performed with respect to the Dirac spinor $u_{i}$ and the occupation probability $n_{i}$, respectively, which leads to the nucleon equation at finite temperature, namely the Dirac equation, and Fermi-Dirac distribution $n_{\tau}(p)$,

[TABLE]

where $\varepsilon(p,\tau)$ is the single-particle energy of the state $i=(p,s,\tau)$, and the spin index $s$ is omitted since the single particle states are degenerated for $s=\pm 1/2$. Notice that the Dirac equation (5) is formally unchanged as compared to the one in Ref. Sun et al. (2008), and the temperature effects lie implicitly in the starred quantities, $\bm{p}^{*}=\bm{p}+\hat{\bm{p}}\Sigma_{V}$, $M^{*}=M+\Sigma_{S}$, and $\varepsilon^{*}=\varepsilon(p)-\Sigma_{0}$, which are reflected by the occupation probability $n_{\tau}(p)$ in the self-energies $\Sigma_{S}$, $\Sigma_{V}$ and $\Sigma_{0}$.

After considering the occupation probability induced by the thermal excitation, the energy density functionals, i.e., the kinetic part $E^{\rm kin}$, the potential parts $E_{\phi}^{D}$ and $E_{\phi}^{E}$ in Eq. (1) can be obtained as,

[TABLE]

where $\mathcal{I}_{\phi}(\tau,\tau^{\prime})$ represents the isospin-related factor, $A_{\phi},B_{\phi},C_{\phi}$ are the angle integral coefficients, and $\hat{P},\hat{M}$ are the hatted quantities, see Ref.Sun et al. (2008) for details. The scalar density $\rho_{s}$, and the baryon density $\rho_{b}$ and its third component $\rho_{b3}$ read as,

[TABLE]

with $\tau=1$ for neutron and $-1$ for proton, respectively. For NLRMF models, an extra contribution from nonlinear self-coupling of mesons $E_{\rm{N.L.}}$ should be appended in the Hartree terms of potential energies $E_{\phi}^{D}$,

[TABLE]

After performing the variation to the potential energy densities, the nucleon self-energy $\Sigma(p)$ is obtained, namely

[TABLE]

Via a self-consistent procedure of self-energies, the properties of thermal nuclear matter can be determined at given density $\rho_{b}$, the isospin asymmetry $\delta=(N-Z)/(N+Z)$ and the temperature $T$. With the free energy $F$, the pressure of thermal nuclear matter is then derived from the thermodynamic relation,

[TABLE]

According to the definition of free energy, the pressure in Eq. (15) can be divided further into

[TABLE]

where the terms $P_{E_{0}}$ and $P_{E_{S}}$ are originated from the binding energy per nucleon $E/\rho_{b}-M$ which is divided further by the isospin symmetric part $E_{0}$ and the symmetry energy related one $\delta^{2}E_{S}$, and $P_{S}$ from the entropy. For instance, the symmetry energy related part $P_{E_{S}}$ is expressed as

[TABLE]

It is clear that the contribution of $P_{E_{S}}$ is discarded in symmetric nuclear matter as $\delta=0$. For asymmetric matter, $P_{E_{S}}$ plays a role in the total pressure, and from the definition its value is found to be ascribed qualitatively to the density slope $L$ of symmetry energy at saturation density $\rho_{0}$. Since $L$ is denoted as

[TABLE]

one then find $P_{E_{S}}\propto\rho_{b}L$ approximately at a given density $\rho_{b}$ (actually fulfilled strictly at $\rho_{0}$).

To reveal the liquid-gas (LG) phase transition in thermal nuclear matter, one needs to solve the phase coexistence equations,

[TABLE]

which correspond to the Gibbs conditions, i.e., the identical pressures and chemical potentials for both liquid (L) and gas (G) phases at given temperature $T$. When solving the phase coexistence equations, the stability condition shall be also satisfied as,

[TABLE]

At the critical points of LG phase transition, the temperature, density and pressure of nuclear matter is denoted as $T_{C}$, $\rho_{C}$ and $P_{C}$, respectively. For symmetric nuclear matter, the critical point is determined by the inflection point of pressure curve with respect to the baryon density, which is,

[TABLE]

while for asymmetric nuclear matter, instead the critical parameters should be solved by the inflection point of chemical potential isobars, namely,

[TABLE]

Moreover, one can also introduce the critical incompressibility $K_{C}$, which is defined as the second derivative of free energy $F$ with respect to the baryon density $\rho_{b}$ at finite temperature,

[TABLE]

For symmetric nuclear matter, the first condition at critical point in Eq. (21) can be expressed further as

[TABLE]

according to the definition of the pressure Eq. (15). One then readily find a relation between the critical parameters $K_{C}$ and $P_{C}$,

[TABLE]

which makes an alternative way to determine the critical point of LG phase transition in symmetric nuclear matter.

III Results and discussion

In this work, the analysis based on the CDF theory will be carried out by three kinds of meson-exchange types of CDF functionals: (1) NLRMF functionals NL1Reinhard (1989), NLZBender et al. (1999), NLZ2 Bender et al. (1999), NL3Lalazissis et al. (1997), NL3∗ Lalazissis et al. (2009), NL-SHLalazissis et al. (1997), NL$\rho$Liu et al. (2002), TM1Sugahara and Toki (1994), TM2Sugahara and Toki (1994), TMAToki et al. (1995), GL-97Glendenning (2000), PK1Long et al. (2004) and PK1RLong et al. (2004); (2) DDRMF functionals TW99Type and Wolter (1999), DD-ME1Niks̆ić et al. (2002), DD-ME2Lalazissis et al. (2005) and PKDDLong et al. (2004); (3) DDRHF functionals PKO1Long et al. (2006), PKO2Long et al. (2008) and PKO3Long et al. (2008).

In ordinary nuclear matter, the integration over momentum $p$ is carried from zero to Fermi momentum $p_{F}$. For the nuclear matter at finite temperature, the thermal excitation will lead to the spreading of the valence nucleons over the states nearby the Fermi surface, such that the integration over $p$ shall be done from zero to infinity. Several numerical techniques to this kind of integration have been discussed such as in Ref. Gong et al. (2001). However, for the cases concerned in this work where the temperature is lower than 20 MeV, the diffusion of Fermi surface is somewhat weak so that the occupation probability $n_{i}$ drops down to zero promptly. It has been checked that a Gauss-Legendre integration up to about $5p_{F}$, as momentum cutoff condition adopted here, has guaranteed the convergence numerically in momentum space. Moreover, the phase coexistence equations (19) as a set of non-linear equations are solved numerically with the Powell hybrid method Powell (1970), which overcomes the deficiency of possible divergence compared to the classical Newton-Raphson method by introducing a ”hybrid algorithm” in the iteration of Jacobian matrix.

III.1 Critical point properties of LG phase transition in thermal nuclear matter

Critical parameters are very important characteristic quantities in determining properties of liquid-gas phase transition, among which the critical temperature $T_{C}$ is especially concerned. For symmetric nuclear matter, $T_{C}$ is estimated in the range of $10\thicksim 20$ MeV in previous studies. To reduce its predicted uncertainty theoretically, the correlations among various critical parameters of LG phase transition account for and need to be discussed not only in symmetric but in asymmetric nuclear matter.

For symmetric nuclear matter, the critical point of LG phase transition is determined according to Eq. (21), which is relevant to the inflection point on $T_{C}$ isotherm. Taking the RHF functional PKO1 as an example, the calculated pressure curves of thermal symmetric nuclear matter with the baryon density $\rho_{b}$ are shown in Fig. 1. At finite temperature, it is seen that the pressure curves behave a characteristic $S$ shape of van der Waals-like isothermGoodman et al. (1984); Silva et al. (2008); Rios et al. (2008); Rios (2010); Vovchenko et al. (2015); Lourenço et al. (2016); Vovchenko (2017). When the temperature is lower than a certain value which defines the critical temperature $T_{C}$, the pressure curve presents a non-monotonic trend with increasing density. Accordingly, the spinodal instability would occur in the density range between two extreme points (points in which $\partial P/\partial\rho_{b}=0$), leading to the LG phase transition. For PKO1, $T_{C}$ is found to be $14.53$ MeV, and the critical pressure $P_{C}$ of LG phase transition is 0.191 MeV$\cdot$fm*-3* (see Table 1 for others). For classical van der Waals (VDW) gas, it has been deduced that a linear relation between $T_{C}$ and $P_{C}$ exists as $T_{C}/P_{C}=8b$, where $b$ is the VDW parameter that describes repulsive interaction Vovchenko et al. (2015). After considering Fermi statistics, the VDW-like equation of state could be established analytically for thermal nuclear matter Goodman et al. (1984); Lourenço et al. (2016), and it is found the $T_{C}-P_{C}$ linear relation is still preserved under several approximations such as those to the effective mass and equation of motion of nucleons. In the following, we will check such a relation numerically within CDF functionals.

Figure 2 shows the critical parameters of LG phase transition for symmetric nuclear matter given by the selected CDF functionals, which are determined based on Eq. (21). It is seen that $T_{C}$ given by these CDF functionals vary from $13\thicksim 18$ MeV, the range of which is too large to constrain $T_{C}$ or compare directly with the experimental data. Hence, some model-independent relations or correlations of the critical parameters within themselves or with bulk properties of nuclear matter, once verified, would be very helpful to a better constraint of $T_{C}$. In the left panel of Fig. 2, the critical temperature $T_{C}$ is plotted versus the corresponding critical pressure $P_{C}$ given by the selected CDF functionals. A linear correlation between $T_{C}$ and $P_{C}$, which has been claimed well in ideal VDW gas, is achieved approximately in present numerical studies. Then the results can be fitted in terms of

[TABLE]

where $a=24.55~{}{\rm fm^{3}}$, $b=9.67~{}{\rm MeV}$ for symmetric nuclear matter, and the Pearson’s coefficient is $r=0.967$ which indicates notably the robustness of such a $T_{C}-P_{C}$ linear correlation to the choice of models.

In addition, from the linear relationship between the critical incompressibility $K_{C}$ and the ratio $P_{C}/\rho_{C}$ illustrated in Eq. (25) for symmetric nuclear matter, it is natural and readily to establish a $T_{C}-\rho_{C}K_{C}$ correlation via $T_{C}-P_{C}$ one. Since the values of $\rho_{C}$ are close to each other for the CDF functionals, as seen in Table 1, one would then expect a possible $T_{C}-K_{C}$ correlation. As shown in the right panel of Fig. 2, the linear correlation between $T_{C}$ and $K_{C}$ is verified numerically in CDF approaches, with the Pearson’s correlation coefficient $r=-0.969$. For convenience, one usually introduce a dimensionless parameter to describe such a correlation, namely, the compressibility factor $Z_{C}$ at critical point of LG phase transition which is defined as

[TABLE]

It has been checked that the values of $Z_{C}$, with samples listed in Table 1 (and also in Table 2 for asymmetric matter), are in general located around $0.3$ from present CDF calculations, in consistence with the previous results analyzed by various density functional approaches Goodman et al. (1984); Rios (2010); Lourenço et al. (2017). Furthermore, one notice that these values are also compatible to (although always smaller than) those from standard VDW gas which is known as 3/8 Vovchenko et al. (2015); Lourenço et al. (2017), indicating again the VDW gas-like nature of thermal nuclear matter in CDF approaches.

To better constrain the critical temperature, its dependence on a series of bulk quantities of cold nuclear matter should be investigated as well. In previous works Rios et al. (2008); Rios (2010); Lourenço et al. (2016), it is suggested that the critical temperature $T_{C}$ of thermal nuclear matter could be correlated with the properties of symmetric nuclear matter at zero temperature such as the incompressibility at saturation density $K_{0}$. It is essential to confirm the conclusion within various different nuclear models. If a distinct correlation do exists between $T_{C}$ and $K_{0}$, the constraint on $K_{0}$ from a lot of experiments Garg and Colò (2018), for example the giant monopole resonance Stone et al. (2014), could be utilized to get more strict value of $T_{C}$. In Table 1, in addition to the critical parameters of LG phase transition, we also list the incompressibility coefficients $K_{0}$ from six characteristic CDF functionals. However, by checking the Pearson’s correlation coefficient (here as $r=0.623$), the quantity $T_{C}$ is only weakly dependent on $K_{0}$ in CDF calculations (20 selected functionals). Thus, careful analysis from different parts of the free energy need further so as to find the influence of finite temperature on the incompressibility.

To investigate the change of the feature of LG phase transition with isospin asymmetry $\delta$, it is valuable to look for the possible correlations of LG critical parameters in asymmetric nuclear matter as well. When $\delta\neq 0$, the critical point should be determined by the condition in Eq. (22). The $T_{C}-P_{C}$ correlation is checked again but for the case of $\delta=0.5$, as shown in the left panel of Fig. 3. In comparison with the symmetric one shown in Fig. 2, the critical temperature $T_{C}^{\delta=0.5}$ is no longer linearly correlated well with the critical pressure, as $r=0.880$, while its value locates in the range of $11\thicksim 15.5$ MeV. From Eq. (16), the contribution to $P_{C}^{\delta}$ from different components could be quantified and be helpful to clarify the physical origin of such a destruction of correlation. As compared to the symmetric part $P_{C,E_{0}}^{\delta}$ and entropy part $P_{C,S}^{\delta}$, it is found that the symmetry energy related part $P_{C,E_{S}}^{\delta}$ actually has a larger model dependence. Therefore, it is rational that the exclusion of $P_{C,E_{S}}^{\delta}$ from $P_{C}^{\delta}$, namely $P_{C}^{\delta}-P_{C,E_{S}}^{\delta}$, exhibits a partly recovered correlation ($r=0.955$) with $T_{C}^{\delta}$ for asymmetric nuclear matter, as seen in the right panel of Fig. 3. Therefore, one could introduce a possible linear relation as

[TABLE]

where $c=29.10~{}{\rm fm^{3}}$, $d=10.57~{}{\rm MeV}$ for the case of $\delta=0.5$. Besides, the Pearson’s correlation coefficient between $T_{C}^{\delta=0.5}$ and $K_{0.1}$ (or $K_{0}$) is calculated. The smaller value of $r=0.571$ ($r=0.585$) than that in $\delta=0$ case indicates that the correlation between the critical temperature and the incompressibility become worse due to the isospin asymmetry.

To compare further with the results of symmetric nuclear matter in table 1, Table 2 shows the critical parameters of LG phase transition at isospin asymmetry $\delta=0.5$. It is found that the critical temperatures $T_{C}^{\delta=0.5}$ are smaller than those of symmetric matter, while the critical density and pressure increase slightly as compared to $\delta=0$ case. Recently, one has paid considerable attention to the dependence of LG critical parameters on the isospin asymmetry from experiments of the nuclear caloric curves, where the evolution of the limiting temperature of finite nuclei with mass number and isospin is illustrated Le Fèvre et al. (2005); Sfienti et al. (2009); McIntosh et al. (2013a, b). Here, the isospin dependence of the critical temperature is demonstrated as well within the CDF theory, as seen in Fig. 4 with several CDF functionals. It is revealed that the critical temperature $T_{C}^{\delta}$ goes down monotonically with increasing isospin asymmetry $\delta$, in agreement with previous analysis adopting CDF approaches Sharma and Pal (2010); Zhang and Jiang (2013). At small isospin asymmetry, the change of $T_{C}^{\delta}$ with $\delta$ is moderate, while the value is suppressed drastically after $\delta\gtrsim 0.5$, despite a slight model dependence.

As has been discussed around Eq.(17), $P_{C,E_{S}}$ could be somewhat related to the symmetry energy, correspondingly being the critical temperature $T_{C}^{\delta}$ at finite isospin asymmetry. It is helpful to elucidate such a possible relation numerically with present CDF calculations. Assuming the evolution of $T_{C}^{\delta}$ with $\delta$ is controlled by the symmetry energy (some discussion see Refs. Sharma and Pal (2010); Zhang and Jiang (2013)), it is worthwhile to define a scaled critical temperature $\Delta T_{C}^{\delta}$ for a certain $\delta$ as

[TABLE]

Based on Eqs. (26) and (28), $\Delta T_{C}^{\delta}$ is then expressed as

[TABLE]

The terms inside square brackets contribute an intercept $\backsimeq[-3.3,-1.6]$ MeV within the selected CDF functionals, showing a weak model dependence. In combination with the relation $P_{E_{S}}\propto\rho_{b}L$, then it is deduced that $\Delta T_{C}^{\delta}$ is proportional to $\rho_{C}^{\delta}L$ roughly, as the ratio $\rho_{C}^{\delta}/\rho_{0}$ can be treated as an constant approximately in CDF calculations Lourenço et al. (2017). Utilizing the selected 20 CDF functionals, such a $\Delta T_{C}^{\delta}-\rho_{C}^{\delta}L$ relation or $\Delta T_{C}^{\delta}/\rho_{C}^{\delta}-L$ relation equivalently is verified numerically, as shown in Fig. 5 for the case of $\delta=0.5$. The Pearson’s coefficient $r=0.932$ indicates a good linear correlation between the scaled critical temperature and the density slope parameter of symmetry energy. With the constraint on the density slope $L=58.7\pm 28.1$ MeV taken from Ref. Oertel et al. (2017), it is then proposed from Fig. 5 that the value of $\Delta T_{C}^{\delta}/\rho_{C}^{\delta}$ at $\delta=0.5$ is about $31.2\pm 6.6$ MeV$\cdot$fm3.

Hence, it is seen in CDF cases that although the critical temperature of LG phase transition has a clear model dependence, both in symmetric and asymmetric nuclear matter, several linear and model-independent correlations between the critical temperature and other LG critical parameters or bulk properties of nuclear matter could exist. These correlations become worse when nuclear matter is more asymmetric, resulting mainly from the uncertainty of symmetry energy related contributions. However, a linear correlation between the critical temperature at isospin asymmetric case and the density slope of symmetry energy is unveiled, which paves a possible way to constrain the critical parameters of LG phase transition.

III.2 Properties of LG phase diagram in thermal nuclear matter

Phase diagram provides essential information about matter structure at a certain circumstance. Specifically, the liquid-gas phase diagram for thermal nuclear matter is substantial to understand several aspects in heavy-ion collision and nuclear astrophysics Chomaz et al. (2004); Lattimer and Prakash (2016). Following the above discussion, it is convenient to study LG phase diagrams within CDF functionals, which for the case of symmetric nuclear matter are given in Fig. 6. The boundary between two phases can be fixed by solving Eqs. (19). It is found that the phase diagram is divided into three regions: (I) the gas phase at low density; (II) the mixed phase and (III) the liquid phase at high density. Because of the deviation of $T_{C}$ as listed in Table 1, there exist an obvious model dependence (particularly the temperature) of the critical points of LG phase transition (filled circles) for the selected CDF functionals, where the RHF functional PKO2 gives the highest $T_{C}$.

The LG phase diagram for symmetric nuclear matter has also been discussed in many works quantitativelyRios et al. (2008); Rios (2010); Fedoseew and Lenske (2015); Lattimer and Prakash (2016). For the case of asymmetric nuclear matter, the influence of symmetry energy is supposed to be important in deciding the pattern of the phase diagram Müller and Serot (1995); Xu et al. (2007); Sharma and Pal (2010); Zhang and Jiang (2013). It is argued that the behavior of liquid-gas phase coexistence could be correlated with the symmetry energy at saturation density Müller and Serot (1995) or just its density slope Xu et al. (2007); Sharma and Pal (2010); Zhang and Jiang (2013). To investigate such a topic more quantitatively, here we plot the LG phase diagrams with the selected CDF functionals in Fig. 7 by fixing $T$ for each functional at its own critical temperature $T_{C}^{\delta=0.5}$ given in Table 2, which is different from the common treatment of exploring at a constant temperature.

In order to clarify the structure of phase diagram, it is salutary to define three characteristic points: the critical pressure (CP) point (filled circles) determining the maximum pressure $P_{\rm{CP}}$ that the LG phase transition could occur, the maximum asymmetry (MA) point which is given by $\delta=\delta_{\max}$ of the gas phase during the phase transition, and the equal concentration (EC) point of the phase diagram at $\delta=0$. When the phase diagram is plotted in manner of fixing $T=T_{C}^{\delta}$, the pressure $P_{\rm{CP}}$ at CP point is just $P_{C}^{\delta}$ mentioned in subsection III.1. Correspondingly, the phase diagram is divided into two branches by the CP and EC points, namely the high-density liquid phase line (left branch) and the low-density gas phase line (right branch), and the region surrounded by two lines is the phase coexistence area. When $\delta$ is larger than one at CP point, namely $\delta>\delta_{\rm{CP}}$ ($\delta_{\rm{CP}}=0.5$ in the case of Fig. 7), the system will not change completely into the liquid phase Sharma and Pal (2010). The positions of characteristic points then determine more or less the size of coexistence area, namely, the lower(higher) $P_{\rm{EC}}$ ($P_{\rm{CP}}$) is, the larger the phase coexistence area becomes. Since the small divergence of $\delta_{\max}$ for the selected models as seen in Fig. 7, one can adopt the pressure difference between CP and EC points, i.e., $P_{\rm{CP}}-P_{\rm{EC}}$, to indicate the size of phase coexistence area of LG phase diagrams. It is seen that $P_{\rm{CP}}$ and $P_{\rm{EC}}$ is clearly model dependent from the picture, leading to the uncertainty of diagram pattern. For instance, a remarkable enhancement of $P_{\rm{EC}}$ is given by PKO2 functional, while its $P_{\rm{CP}}$ is generally comparable with other model predictions, so that a relatively smaller LG phase coexistence area appears in PKO2 case.

From Fig. 7, it is necessary to extract the pressure values at various characteristic points, as listed in Table 3, so as to explore how the bulk properties of nuclear matter affect the size of phase coexistence area. For the difference $P_{\rm{CP}}-P_{\rm{EC}}$, PKDD gives the largest value among all functionals, corresponding to the most extensive area of phase coexistence. With the help of Eq. (16), the contribution of $P_{\rm{CP}}-P_{\rm{EC}}$ can be separated into

[TABLE]

where $\Delta P_{E_{0}},\Delta P_{E_{S}}$ and $\Delta P_{S}$ represent isospin symmetric, isospin asymmetric (symmetry energy) and entropy part, respectively. Since $P_{E_{S}}=0$ at EC point, $\Delta P_{E_{S}}=P_{C,E_{S}}^{\delta}$ which is just the value at CP point. As revealed in Table 3, the value of $P_{\rm{CP}}-P_{\rm{EC}}$ is mainly ascribed to the contribution of symmetry energy part $P_{C,E_{S}}^{\delta}$, while $\Delta P_{E_{0}}$ and $\Delta P_{S}$ almost cancel each other although their respective contributions are relatively large. From Eq. (31), one can expect directly a linear correlation between $P_{\rm{CP}}-P_{\rm{EC}}$ and $\Delta P_{E_{S}}$ as well, which is drawn in Fig. 8. The Pearson’s correlation coefficient is obtained as good as $r=0.968$, indicating the significant role of the symmetry energy in the size of phase coexistence area in LG phase diagram, in agreement with the conclusion in previous works Müller and Serot (1995); Xu et al. (2007); Sharma and Pal (2010); Zhang and Jiang (2013). As an alternative case, the linear correlation between $P_{\rm{CP}}-P_{\rm{EC}}$ and $\Delta P_{E_{S}}$ for LG phase diagram at temperature $T=10$ MeV is done as well, and the above conclusion is confirmed with a correlation coefficient $r=0.898$.

III.3 Correlations of critical parameters in LG phase diagram

In subsection III.1, the correlations among the critical parameters of LG phase transition in thermal nuclear matter, in particular the critical temperature $T_{C}$, have been discussed, which is supposed to be a possible way to constrain the critical parameters from bulk properties of nuclear matter. The structure of phase diagram could be reflected in a similar way, if the properties at characteristic points are confirmed to be associated with the critical parameters of LG phase transition as well.

In the left panel of Fig. 9, the pressures at EC point $P_{\rm{EC}}$ in phase diagrams of thermal nuclear matter at temperature $T=T_{C}^{\delta=0.5}$ is treated to correlate with the critical parameter $T_{C}^{\delta}$ at $\delta=0.5$. Such correlation tend to be well linear with the Pearson’s coefficient $r=0.978$, which can be illustrated readily from the satisfied $T_{C}^{\delta}$ correlation with $P_{C}^{\delta}-P_{C,E_{S}}^{\delta}$ shown in Fig. 3 in combination with the relationship deduced from Eq. (31) where $P_{C}^{\delta}-P_{C,E_{S}}^{\delta}$ is equivalent to $P_{\rm{CP}}-\Delta P_{E_{S}}$. Besides, it is realized from the right panel of Fig. 9 that the pressure at maximum asymmetry (MA) point $P_{\mathrm{MA}}$ is also relevant to $T_{C}^{\delta=0.5}$ although a relatively smaller $r=0.928$ than one in $P_{\mathrm{EC}}$ case.

At a certain isospin asymmetry $\delta$, two characteristic pressures associated with the gas phase and liquid phase lines can be defined further as $P_{\mathrm{Gas}}$ and $P_{\mathrm{Liquid}}$, which are extracted from Fig. 7. As plotted in Fig. 10 for an example of $\delta=0.3$, both these two quantities are demonstrated to correlate with $T_{C}^{\delta=0.5}$. Furthermore, it is unveiled that $T_{C}^{\delta=0.5}$ and $P_{\mathrm{Gas}}$ address a correlation with $r=0.974$, while for $P_{\mathrm{Liquid}}$ case it has a relatively smaller Pearson’s coefficient of $r=0.933$, suggesting a better linear correlation for gas phase line than that for liquid phase line. The rule is also proved to be satisfied at other isospin asymmetries, as shown in Table 4 for $0\leqslant\delta\leqslant\delta_{\rm{CP}}$. It is found that the Pearson’s correlation coefficients $r_{L}$ for the cases of liquid phase line are always smaller than $r_{G}$ for those of gas phase line when $0<\delta\leqslant\delta_{\rm{CP}}$, which could be interpreted by the fact of a larger CDF model dependence in describing the liquid phase than the gas one since the density of the former is larger and the interaction between nucleons stronger.

From the above discussion, several linear correlations are illustrated between the critical temperature of LG phase transition $T_{C}^{\delta}$ at $\delta\neq 0$ and the characteristic pressures of LG phase diagrams for asymmetric nuclear matter, including $P_{\rm{EC}}$, $P_{\rm{MA}}$, $P_{\rm{Gas}}$ and $P_{\rm{Liquid}}$, which are demonstrated to be better than $T_{C}^{\delta}-P_{C}^{\delta}$ correlation as revealed in Fig. 3 of subsection III.1. It is then expected that these correlations could be utilized to constrain the structure of LG phase diagram with the progress in determining $T_{C}^{\delta}$ at various isospin asymmetries.

IV Summary

In conclusion, by adopting the covariant density functional theory, namely, NLRMF, DDRMF and DDRHF approaches, the liquid-gas phase transition in thermal nuclear matter specifically its properties at critical point has been studied in this work. The thermal nuclear matter in CDF calculations behaves like van der Waals gas as illustrated in the shape of pressure isotherms and the compressibility factor. It is seen that the critical parameters, including the critical temperature $T_{C}$, critical density $\rho_{C}$, critical pressure $P_{C}$ and critical incompressibility $K_{C}$, are clearly model dependent in both symmetric and asymmetric nuclear matter. However, it is verified numerically within CDF functionals that there exist linear correlations approximately between critical parameters and bulk properties of nuclear matter, such as between $T_{C}$ and $P_{C}$ ($K_{C}$). These correlations become worse for larger isospin asymmetry, which can be attributed from the uncertainty of the contribution $P_{C,E_{S}}^{\delta}$ due to the symmetry energy. Correspondingly, the role of the symmetry energy in the isospin dependence of LG transition parameters is focused further. It is unveiled from the CDF calculations that the scaled quantity $\Delta T_{C}^{\delta}/\rho_{C}^{\delta}$ can be well determined by the density slope of symmetry energy $L$. Thus, more constraints on nuclear symmetry energy would be crucial and necessary to better understand the critical parameters of LG phase transition at various asymmetric isospin. With recent empirical value of $L$ Oertel et al. (2017), the value of $\Delta T_{C}^{\delta}/\rho_{C}^{\delta}$ at $\delta=0.5$ is suggested to be about $31.2\pm 6.6$ MeV$\cdot$fm3.

Then in the later parts of Sec. III, the structure of LG phase diagram of thermal nuclear matter is investigated, especially via analysis of the pressure associated with equation of state or entropy. It is found that the size of LG phase coexistence area, determined approximately by the pressure difference $P_{\rm{CP}}-P_{\rm{EC}}$, is well correlated with the pressure part $P_{C,E_{S}}^{\delta}$ due to symmetry energy, which is in agreement with the conclusion in previous studies. After extracting the pressure values at several characteristic points in LG phase diagrams, namely, $P_{\rm{EC}}$, $P_{\rm{MA}}$, $P_{\rm{Gas}}$ and $P_{\rm{Liquid}}$, their linear correlations with the critical temperature $T_{C}^{\delta}$ at non-zero isospin asymmetry are confirmed. Therefore, a possible way is established to depict the pattern of LG phase diagram directly from the critical temperature at virous isospin asymmetries. If $T_{C}^{\delta}$ can be well constrained such as by the density slope $L$ of symmetry energy, the uncertainty of theoretical prediction to the LG phase diagram will be diminished substantially owing to these correlations, and the physics of liquid-gas phase diagram of thermal nuclear matter will be clarified explicitly.

Acknowledgements.

The authors are grateful to Prof. Wen Hui Long and Dr. Jian Min Dong for helpful discussions. This work is partly supported by the National Natural Science Foundation of China (Grant Nos. 11675065 and 11875152).