You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: 'Within the framework of BCS theory and Skyrme-Hartree-Fock model, we employ various microscopic pairing gaps and effective pairing interactions to study pairing effects on the equation of state (EOS) of neutron matter and the symmetry energy at subsaturation densities. We find pairing effects may have considerable contributions to the EOS of neutron matter at very low densities ($\lesssim 0.02~\rm{fm}^{-3}$), while only have a small impact on the symmetry energy at subsatruation densities. In addition, the reliability of the parabolic approximation for the isospin asymmetry dependence of nuclear matter EOS with pairing correlations included is also discussed.'
author:
- 'Zhen Zhang[^1]'
- 'Lie-Wen Chen[^2]'
title: Pairing effects on neutron matter equation of state and symmetry energy at subsaturation densities
---

Introduction
============

Nuclear symmetry energy $E_{\mathrm{sym}}(\rho)$ and the equation of state (EOS) of pure neutron matter (PNM) $E_{\mathrm{PNM}}(\rho)$ have profound impact on many important physics problems in nuclear physics and astrophysics[@Lat04; @Ste05; @Bar05; @LCK08; @Oer17] as well as some issues in new physics beyond the standard model [@Hor01; @Sil05; @WLC09; @ZZC14; @ZSC15]. For instance, the density dependence of the symmetry energy or neutron matter EOS at subsaturation densities is intimately related to the neutron skin thickness of finite nuclei [@Bro00; @Bro13; @Typ01; @Hor01a; @Fur02; @Yos04; @Tod05; @CLK05; @CKLX10; @Cen09; @Roc11; @Agr12; @Zha13], the properties of neutron star crust [@Hor01a; @Xu09; @Duc11; @Zha14], the cluster formations in nuclear matter at low densities [@Kow07; @Nat10; @Wad12; @Typ10; @ZhangChen17], and the isospin diffusion in heavy ion collisions at Fermi energies [@Tsa04; @Tsa09; @CKL05; @LC05].

In last few decades, a lot of efforts based on various phenomenological models or microscopic theories have been devoted to exploring $E_{\mathrm{sym}}(\rho)$ and $E_{\mathrm{PNM}}(\rho)$ at subsaturation densities. From analyses of various experimental observables using phenomenological models, our knowledge on nuclear matter EOS at subsaturation densities have been significantly improved. For example, in the density region of about $\rho_0/3\sim\rho_0$ ($\rho_0\approx0.16~\mathrm{fm}^{-3}$ is the nuclear saturation density), constraints on nuclear symmetry energies have been extracted from analyses of nuclear masses [@Zha13; @Bro13] as well as isobaric analog states and neutron skin data [@Dan14] using Skyrme-Hartree-Fock (SHF) model, and also from transport model analyses of experimental observables in mid-peripheral heavy ion collisions of Sn isotopes [@Tsa09]; around $\rho_0/3$, both $E_{\mathrm{sym}}(\rho)$ and $E_{\mathrm{PNM}}(\rho)$ have been well determined by the electric dipole polarizability in $^{208}$Pb [@Zha15]. These constraints from phenomenological models are essentially consistent and the uncertainty can be as small as $\approx 1~\mathrm{MeV}$ (see, e.g., Ref. [@Zha15]).

In addition, great progress in constraining $E_{\mathrm{sym}}(\rho)$ or $E_{\mathrm{PNM}}(\rho)$ at subsaturation densities has also been made by ab initio calculations, especially with the development of modern chiral nuclear force. For examples, the QCD sum rules have been shown to provide important information on EOS of neutron matter at subsaturation densities [@Cai18]. Based on modern chiral nuclear forces, various many-body theories such as chiral effective field theory [@Mac11; @Tew13], self-consistent Green function method [@Car14], coupled-cluster method [@Hag14] and Quantum Monte Carlo technique [@Car15] have also provide important information from first principle on nuclear matter EOS at subsaturation densities.

It should be pointed out that the constraints on $E_{\mathrm{sym}}(\rho)$ or $E_{\mathrm{PNM}}(\rho)$ from phenomenological models (e.g., SHF model and relativistic mean field theory) are usually obtained without including nucleon pairing correlation, while many ab initio calculations include it. It is well known that the $^1\mathrm{S}_0$ pairing gap in nuclear matter can play an essential role at low densities (e.g., $\rho_0/10$ ) [@Lom99; @Dea03]. Therefore, it is important and interesting to examine the pairing effects on the EOS of nuclear matter at subsaturation densities. This is especially the case when the phenomenological models give more and more stringent constraints on the $E_{\mathrm{sym}}(\rho)$ and $E_{\mathrm{PNM}}$ at subsaturation densities (see, e.g., Ref. [@Zha15]), which provides the main motivation of the present work.

Theoretically great efforts have been made to determine the density dependence of pairing gap in neutron matter using various microscopic many-body approaches, e.g., Bardeen-Cooper-Schrieffer (BCS) theory [@Heb07; @Zha10], correlated-basis-function calculations [@Che93; @Fab05], renormalization group [@Sch03], Brueckner theory [@Cao06] and Quantum Monte Carlo [@Gez10; @Gan08; @Gan09]. In particular, at extremely low densities, or equivalently low $k_{\mathrm{F}}a$ with $a\approx-18.5$ fm being neutron scattering length and $k_{\mathrm{F}}$ neutron Fermi momentum, there exists a well-known analytical pairing gap of $\Delta(k_F)=\frac{1}{(4e)^{1/3}}\frac{8}{e^2}\frac{\hbar^2k_F^2}{2M}\mathrm{exp}\left( \frac{\pi}{2ak_F}\right)$ with $e$ being Euler’s number [@Gor61]. Thanks to the analytical limit and microscopic calculations, the pairing gap in neutron matter at densities $< \rho_0/10$ is under good control. But its higher density behaviors are still largely uncertain [@Gez14]. Moreover, the isospin-dependence of the pairing gap is even more poorly known and microscopic calculations are still inadequate for the study of pairing effects on the symmetry energy. An alternative perspective is starting from effective pairing interactions which are usually constructed by fitting properties of finite nuclei 

[@Kha09; @Yam09; @Yam12] or microscopic pairing gaps [@Mar07] with a hypothetical functional form. Using the effective pairing interactions, one can study the effects of pairing correlation on nuclear symmetry energy within the framework of BCS theory.

In this work, within the framework of BCS theory together with the SHF model, we study pairing effects on the EOS of neutron matter at subsaturation densities by invoking various microscopic pairing gaps or effective pairing interactions. Subsequently, we use several effective pairing interactions to study pairing effects on the symmetry energy at subsaturation densities. The reliability of the parabolic approximation for the isospin asymmetry dependence of nuclear matter EOS is further verified in the case of including the effect of pairing correlations.

The paper is organized as follows. In Sec. \[Sec:PF\] we introduce the effective pairing interactions used in this work and the calculation of pairing energy density. Sec. \[Sec:Results\] presents the results for pairing effects on neutron matter EOS and the symmetry energy at subsaturation densities. We then end the paper with a summary and conclusions.

Effective pairing interaction and pairing energy density {#Sec:PF}
========================================================

In the framework of BCS theory, the equation for the pairing gap $\Delta_q$ ($q=n, p$) in nuclear matter is given by $$\label{Eq:gap}
\Delta_q(k)=-\int \frac{d^3k'}{(2\pi)^3} \frac{v_q(k,k')\Delta_q(k')}{2\sqrt{[\epsilon_q(k')-\lambda_q]^2+\Delta_q^2}},$$ where $v_q$ is the pairing strength in momentum space, $\lambda_q$ is the effective chemical potential, and $\epsilon_q(k)=\hbar^2k^2/2m_q^{\ast}$ is the single-particle kinetic energy with $m_q^*$ being the nucleon effective mass. It should be pointed out that in the case of contact pairing interaction, both $v_q$ and $\Delta_q$ are momentum independent. Note that the momentum independent mean-field potentials have been absorbed into the effective chemical potential, which can be determined by the nucleon density $$\label{Eq:den}
\rho_q = \int \frac{d^3k}{(2\pi)^3}\left[ 1-\frac{\epsilon_q(k)-\lambda_q}{E_q}\right],$$ where $E_q = \sqrt{\left[\epsilon_q(k)-\lambda_q\right]^2+\Delta_q^2}$ is the quasi-particle energy. For contact pairing interactions, once given the pairing strength $v_q(\rho_n, \rho_p)$ \[the pairing gap $\Delta_q(\rho_n, \rho_p)$\] and the nucleon effective mass $m_q^{\ast}(\rho_n, \rho_p)$ at neutron density $\rho_n$ and proton density $\rho_p$, the $\Delta_q(\rho_n, \rho_p)$ \[$v_q(\rho_n, \rho_p)$\] and effective chemical potential $\lambda_q(\rho_n, \rho_p)$ can be easily obtained by solving Eqs. (\[Eq:gap\]) and (\[Eq:den\]). The pairing energy density in asymmetric nuclear matter then can be expressed as $$\begin{aligned}
\label{Eq:enden}
\varepsilon_{\rm{pair}}&=&\sum_{q=n,p}\int \frac{d^3k}{(2\pi)^3}
\left\lbrace \epsilon_q(k)\left[1-\frac{\epsilon_q(k)-\lambda_q}{E_q} \right]
-\frac{1}{2}\frac{\Delta_q^2}{E_q} \right\rbrace  \notag\\
&&-\sum_{q=n,p} \frac{3}{5}\frac{\hbar^2}{2m_q^{\ast}}\rho_q k_{F,q}^2 ,\end{aligned}$$ where $k_{F,q}=(3\pi^2\rho_q)^{1/3}$ is the Fermi momentum. In the weak coupling approximation ($\Delta_q\ll \hbar^2k^2_{F,q}/2m_q^{\ast}$) [@Cha10], the pairing energy density can be approximated as $$\label{Eq:ana}
\varepsilon_{\rm{pair}}=-\frac{1}{2}\left( N_n\Delta_n^2+N_p\Delta_p^2\right),$$ with $N_q = m_q^*k_{F,q}/2\pi^2\hbar^2$ being the density of states.

While the pairing effect on neutron matter EOS can be studied by directly invoking microscopic pairing gaps, for the symmetry energy, due to our poor knowledge on the isospin dependence of pairing gaps, one has to introduce effective pairing interactions. In this work, following Ref. [@Kha10; @Mar14], we use several different effective contact pairing interactions to study the pairing effects on nuclear symmetry energy at subsaturation densities. Note that the integral in Eq.(\[Eq:gap\]) is divergent for contact pairing interactions. Therefore, cutoff momenta are usually introduced in the effective contact pairing interactions in different prescriptions corresponding to different physics problems (see, e.g., Ref. [@Mar07]). In the construction of pairing interactions used in this work, the cutoff is defined with respect to the quasiparticle energy $\sqrt{(\epsilon_q(k)-\lambda_q)^2
+\Delta_q^2}<E_C$, where $E_C$ is the cutoff energy and is taken to be $60$ MeV [@Yam12; @Mar07].

The commonly used effective contact pairing interaction has the form of $$v_q(\bm{r},\bm{r}^{\prime})
=v_0\left[1-\eta\left( \frac{\rho}{\rho_c}\right)^{\alpha} \right]
\delta(\bm{r}-\bm{r}^{\prime}),$$ with $\rho_c =0.16~\rm{fm}^{-3}$ and $v_0$ being the strength parameter. The parameter $\alpha$ is usually taken to be one for simplicity, and then the parameter $\eta$ determines the density dependence of the pairing interaction with $\eta=1$ for surface pairing and $\eta=0$ for volume pairing. In recent studies, $\eta= 0.5$ (mix pairing) is preferred as it can well reproduce the mass dependence of the odd-even mass staggering parameter [@Dob02].

Considering the isospin dependence of the pairing gaps, an extended pairing interaction with the inclusion of isovector density $\rho_1=\rho_n-\rho_p$ has been introduced in Ref. [@Yam09; @Yam12] as $$\begin{aligned}
v_q(\bm{r},\bm{r}^{\prime})&=&
v_0\left[ 1-\eta_0\frac{\rho}{\rho_0}-\eta_1\tau_3\frac{\rho_1}{\rho_0}-\eta_2\left(\frac{\rho_1}{\rho_0} \right)^2 \right] \notag \\
& &\times \delta(\bm{r}-\bm{r}^{\prime}),\end{aligned}$$ where $\tau_3 = 1(-1)$ for $q=n(p)$. The additional $\rho_1$ terms are important to describe the isospin dependence of experimental pairing gaps. In this work, we consider three parametrizations from Ref. [@Yam12]: 1) $v_0 = -370.8~\mathrm{MeV\ fm}^3$, $\eta_0=0.75$ and $\eta_1=\eta_2=0$ (SLy4+IS); 2) $v_0 = -396.47~\rm{MeV\ fm}^3$, $\eta_0=0.75 $, $\eta_1 = 0.270$ and $\eta_2=2.5$ (SLy4+IV); 3)$v_0 = -388.60~\rm{MeV\ fm}^3$, $\eta_0=0.75 $, $\eta_1 = 0.4$ and $\eta_2=2.5$ (LNS+IV

). Here we denote them as ‘IS’ and ‘IV’ since they are isospin-independent (scalar) and isospin-dependent (vector), respectively. All the three parametrization forms are optimized to fit experimental pairing gaps in Hartree-Fock-Bogoliubov calculations using the SLy4 [@Cha98] or the LNS [@Cao06a] Skyrme interaction. It should be noticed that the SLy4 predicts negative neutron-proton effective mass splitting, i.e., $m_n^*<m_p^*$, in neutron-rich matter, while the LNS interaction, which is obtained by fitting predictions of Brueckner-Hartree-Fock calculations, has $m_n^*>m_p^*$. As can be seen from Eq. (\[Eq:ana\]), the isospin-dependence of the nucleon effective mass is related to the isospin behavior of the pairing energy. We will also consider the isospin dependent SLy4+MSH [@Mar07] pairing interaction which is given as [@Mar07] $$\begin{aligned}
v_q^{\mathrm{MSH}}(\bm{r},\bm{r}^{\prime})=v_0
\left[ 1-(1-\tau_3\delta)\eta_s\left( \frac{\rho}{\rho_0} \right)^{\alpha_s}\right. \notag \\
\left.-\tau_3\delta\eta_n\left( \frac{\rho}{\rho_0} \right)^{\alpha_n}\right]
\delta(\bm{r}-\bm{r}^{\prime}),\end{aligned}$$ with $v_0=-448~\mathrm{MeV~fm}^{3}$, $\eta_s=0.598$, $\alpha_s=0.551$, $\eta_n=0.947$, $\alpha_n=0.554$, and $\delta=(\rho_n-\rho_p)/\rho$ being isospin asymmetry. The SLy4+MSH parameter set is determined by fitting neutron pairing gaps in both symmetric nuclear matter and pure neutron matter predicted by microscopic Brueckner calculations [@Cao06] with the effective mass deduced from the SLy4 interaction.

Results and Discussion
======================

 \[Sec:Results\]

![Density dependence of $^1$S$_0$ neutron pairing gap in symmetric nuclear matter (a) and pure neutron matter (b) from BCS calculations using various effective pairing interactions and using the microscopic Argonne V18 two-body force with the LNS mean-field [@Zha10]. Panel (b) also includes corresponding results from calculations using microscopic nuclear forces (see the text for details).[]{data-label="Fig:gap"}](Gap.eps){width="1.\linewidth"}

In Fig. \[Fig:gap\], we show the density dependence of neutron $^1$S$_0$ pairing gaps in symmetric nuclear matter (SNM) (panel (a)) and PNM (panel (b)) obtained from the BCS calculations using SLy4+IS, SLy4+IV, SLy4+MSH and LNS+IV, as well as using the microscopic Argonne V18 two-body force [@Zha10] (denoted as AV$_{18}$ (BCS)). Also included in Fig. \[Fig:gap\] (b) are the corresponding results from three different types of pairing gaps in pure neutron matter predicted by different methods: ‘N3LO Chiral (BCS)’ is calculated using chiral N3LO two-body interaction at BCS level [@Heb07]; ‘Brueckner’ is the prediction of Brueckner theory using the Argonne V18 two-body force and a three-body force [@Cao06]; ‘AFDMC’ is obtained from the auxiliary field diffusion Monte Carlo (AFDMC) method calculation [@Gan08] using realistic two- and three-body nuclear force, i.e., the Argonne $v'_8$ and the Urbana IX. It is seen from Fig. \[Fig:gap\] that pairing gaps from the various effective pairing interactions or approaches have quite different density dependence, especially in neutron matter. Note that the pairing gaps in neutron matter from the BCS calculations with realistic two-body interactions \[i.e.,N3LO Chiral (BCS) and AV$_{18}$ (BCS) \] are in very good agreement with each other and roughly reflect the upper limit of microscopic gaps (see e.g. Ref. [@Gez14]) at low densities. As shown in Fig.\[Fig:gap\] (b), the pairing gaps in PNM from the effective contact pairing interactions and microscopic nuclear forces are rather different, especially at higher densities.

![ Binding energy per neutron in pure neutron matter as a function of density without and with different density-dependent pairing gaps (see the text for details). The inset shows the difference between the neutron matter EOSs with and without pairing correlation.[]{data-label="Fig:Epnm"}](Epnm.eps){width="\linewidth"}

As introduced in Sec. \[Sec:PF\], for contact pairing interaction, the pairing strength can be exactly determined by the pairing gap and the nucleon effective masse. Therefore, one can directly construct contact pairing interactions from microscopic pairing gaps and further study pairing effects on properties of nuclear matter and finite nuclei [@Cha10; @Gor09]. Combining the LNS Skyrme interaction with the four neutron pairing gaps in PNM from microscopic nuclear interactions, we calculate the EOS of pure neutron matter at low densities up to $0.04~\mathrm{fm}^{-3}$ and the results are shown in Fig. \[Fig:Epnm\]. For comparison, Fig. \[Fig:Epnm\] also includes the prediction of Hartree-Fock calculation using LNS interaction without neutron pairing correlation. It is seen that the neutron pairing has negative contribution to the EOS of PNM, especially around the density of $ \rho =0.005~\mathrm{fm}^{-3}$. To more clearly clarify the effect of pairing correlation, we define $E_{\mathrm{pair}}$ as the difference between the EOSs with and without neutron pairing, and exhibit $E_{\mathrm{pair}}$ as a function of neutron density in the inset of Fig. \[Fig:Epnm\]. One sees that the neutron pairing may play a considerable role at very low densities (e.g., in the case of ‘N3LO Chiral (BCS)’ and ‘AV$_{18}$(BCS)’, the neutron matter EOS is reduced by about $20\%$ at densities around $\rho= 0.005~\mathrm{fm}^{-3}$), but its effect turns to be negligible towards higher densities (e.g., $E_{\mathrm{pair}}<0.25~\mathrm{MeV}$ above about $\rho= 0.02~\mathrm{fm}^{-3}$). We would like to point out that the conclusion is essentially independent on the choice of Skyrme interaction, since at such low densities the effective masses of neutron are approximately equal to the neutron bare mass.

Furthermore, we study the pairing effects on the symmetry energy using the isospin-dependent pairing gap ‘AV$_{18}$ (BCS)’ [@Zha10] and the four effective pairing interactions, i.e., SLy4+IS, SLy4+IV, LNS+IV and SLy4+MSH interactions introduced in Sec. \[Sec:PF\]. In Fig. \[Fig:Esym\] (a), we show the density dependence of symmetry energy at subsaturation densities below $0.12~\mathrm{fm}^{-3}$ obtained from Hartree-Fock calculations using SLy4 and LNS interactions without including pairing correlation. Fig. \[Fig:Esym\] (b) shows contributions of pairing correlation on the symmetry energy, $E_{\mathrm{sym}}^{\mathrm{pair}}(\rho)$, for the ‘AV$_{18}$ (BCS)’ and the four effective pairing interactions. Here $E_{\mathrm{sym}}(\rho)$ is numerically calculated according to its definition $E_{\mathrm{sym}}(\rho)=\left.\frac{1}{2!}\frac{\partial ^2E(\rho,\delta)}{\partial \delta^2}\right\vert_{\delta=0}$ with $E(\rho,\delta)$ being the EOS of asymmetric nuclear matter. It can be seen that the ‘AV$_{18}$ (BCS)’ pairing gaps lead to rather small and thus negligible effects on the symmetry energy. For effective pairing interactions, while the SLy4+IS interaction always reduces the symmetry energy, SLy4+MSH, SLy4+IV and LNS+IV interactions can provide either positive or negative contributions depending on the density. For the three effective pairing interactions combined with SLy4 interaction, the magnitude of $E_{\mathrm{sym}}^{\mathrm{pair}}(\rho)$ is less than about $0.5$ MeV, while in the case of LNS+IV, the $E_{\mathrm{sym}}^{\mathrm{pair}}(\rho)$ reaches a maximum of about $0.9$ MeV at the density of $\rho=0.08~\mathrm{fm}^{-3}$. Although the SLy4+IV and LNS+IV pairing interactions are constructed by fitting the same experimental data with the same form of parametrization, the effect of LNS+IV on the symmetry energy is much larger than that of SLy4+IV, which could be related to their different isospin dependence of nucleon effective masses, i.e.,the SLy4 and LNS interactions predict respectively negative and postitive neutron-proton effective mass splittings in neutron-rich nuclear matter. It is seen that the maximum of $E_{\mathrm{sym}}^{\mathrm{pair}}(\rho)$ is only about $4.5\% $ of the $E_{\mathrm{sym}}(\rho)$ (about $20$ MeV) and thus usually can be neglected. Nevertheless, with the more in-depth research and more precise constraints on the symmetry energy at low densities, the pairing effects would become nonnegligible. At this point, we would also like to point out that the cluster formations in nuclear matter at subsaturation densities have a substantial impact on nuclear symmetry energy [@Nat10; @Typ10; @ZhangChen17]. It would be nice to investigate simultaneously the effects of pairing correlations and cluster ormation in a self-consistent framework, and this may be pursued in future.

It is also interesting to check the pairing effects on the applicability of the widely used parabolic approximation for the isospin asymmetry dependence of nuclear matter EOS, in which the symmetry energy is approximated to be $E_{\mathrm{sym}}^{\mathrm{PA}}=E_{\mathrm{PNM}}(\rho)-E_0(\rho)$. To this end, we present in Fig. \[Fig:Esym\](c) the deviation of the symmetry energy in parabolic approximation from its exact value, namely, $\Delta E(\rho)=E_{\mathrm{sym}}^{\mathrm{PA}}(\rho)- E_{\mathrm{sym}}(\rho)$, for different pairing interactions. It is seen that the $\Delta E(\rho)$ is negligible compared with the magnitude of $E_{\mathrm{sym}}(\rho)$ (the relative deviation is about $2\%$). One can thus conclude that the parabolic approximation for the isospin asymmetry dependence of nuclear matter EOS is still valid for the EOS of asymmetric nuclear matter with contributions from pairing correlations.

Conclusion
==========

Within the framework of BCS theory together with the SHF model, we have investigated the pairing effects on the equation of state of pure neutron matter and the symmetry energy at subsaturation densities. For neutron matter, invoking microscopic pairing gaps together with the LNS Skyme interaction, we have found the pairing correlations may have essential impact on the EOS of neutron matter at very low densities below about $0.02~\mathrm{fm}^{-3}$, while they turn to be negligible at higher densities. For symmetry energy, using various effective pairing interactions and the ‘AV$_{18}$ (BCS)’, we have found the inclusion of pairing correlations only slightly affect the magnitude of the symmetry energy and usually can be safely neglected. In addition, the parabolic approximation for the isospin asymmetry dependence of nuclear matter EOS has been proved to be a reasonable approximation in the case of including the pairing correlations.

We thank Shi-Sheng Zhang for providing us the isospin-dependent pairing gap from BCS calculation using the Argone V$_{18}$ two-body nuclear force. This work was supported in part by the National Natural Science Foundation of China under Grant No. 11625521 and 11905302, the Major State Basic Research Development Program (973 Program) in China under Contract No. 2015CB856904, the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education, China, and the Science and Technology Commission of Shanghai Municipality (11DZ2260700).

[99]{}

J.M. Lattimer and M. Prakash, Science **304**, 536 (2004); Phys. Rep. **442**, 109 (2007). A.W. Steiner, M. Prakash, J.M. Lattimer, and P. J. Ellis, Phys. Rep. **411**, 325 (2005). V. Baran, M. Colonna, V. Greco, and M. Di Toro, Phys. Rep. **410**, 335 (2005). B.A. Li, L.W. Chen, and C.M. Ko, Phys. Rep. **464**, 113 (2008). M. Oertel, M. Hempel, T. Klahn, and S. Typel, Rev. Mod. Phys. **89**, 015007 (2017). C.J. Horowitz, S.J. Pollock, P.A. Souder, and R. Michaels, Phys. Rev. C **63**, 025501 (2001); T. Sil, M. Centelles, X. Vi$\tilde{\text{n}}$as, and J. Piekarewicz, Phys. Rev. C **71**, 045502 (2005). D.H. Wen, B.A. Li, and L.W. Chen, Phys. Rev. Lett. **103**, 211102 (2009). H. Zheng, Z. Zhang, and L.W. Chen, J. Cosmo. Astropart. Phys. **08**, 011 (2014). H. Zheng, K.J. Sun, and L.W. Chen, Astrophys. J. **800**, 141 (2015).

B.A. Brown, Phys. Rev. Lett. **85**, 5296 (2000). B.A. Brown, Phys. Rev. Lett. **111**, 232502 (2013). S. Typel and B.A. Brown, Phys. Rev. C **64**, 027302 (2001).

R.J. Furnstahl, Nucl. Phys. A **706**, 85 (2002). S. Yoshida and H. Sagawa, Phys. Rev. C **69**, 024318 (2004). B.G. Todd-Rutel and J. Piekarewicz, Phys. Rev. Lett. **95**, 122501 (2005). L.W. Chen, C.M. Ko, and B.A. Li Phys. Rev. C **72**, 064309 (2005). M. Centelles, X. Roca-Maza, X. Vin$\tilde{\text{a}}$s, and M. Warda, Phys. Rev. Lett. **102**, 122502 (2009). X. Roca-Maza, M. Centelles, X. Vin$\tilde{\text{a}}$s, and M. Warda, Phys. Rev. Lett. **106**, 252501 (2011). L.W. Chen, C.M. Ko, B.A. Li, and J. Xu, Phys. Rev. C **82**, 024321 (2010).

B.K. Agrawal, J.N. De, and S.K. Samaddar, Phys. Rev. Lett. **109**, 262501 (2012).

Z. Zhang and L.W. Chen, Phys. Lett. B **726**, 234 (2013). C.J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. **86**, 5647 (2001).

J. Xu, L. W. Chen, B. A. Li, and H. R. Ma, Phys. Rev. C **79**, 035802 (2009); Astrophys. J. **697**, 1549 (2009). C. Ducoin, J. Margueron, C. Providencia, and I. Vidana, Phys. Rev. C **83**, 045810 (2011) Z. Zhang and L.W. Chen, Phys. Rev. C **90**, 064317 (2014)

S. Kowalski, J.B. Natowitz, S. Shlomo, R. Wada, K. Hagel, J. Wang, T. Materna, Z. Chen, Y.G. Ma, L. Qin, A.S. Botvina, D. Fabris, M. Lunardon, S. Moretto, G. Nebbia, S. Pesente, V. Rizzi, G. Viesti, M. Cinausero, G. Prete, T. Keutgen, Y. El Masri, Z. Majka, and A. Ono, Phys. Rev. C **75**, 014601 (2007).

J. B. Natowitz, G. Röpke, S. Typel, D. Blaschke, A. Bonasera, K. Hagel, T. Klähn, S. Kowalski, L. Qin, S. Shlomo, R. Wada, and H.H. Wolter, Phys. Rev. Lett. **104**, 202501 (2010).

R. Wada, K. Hagel, L. Qin, J. B. Natowitz, Y.G. Ma, G. Röpke, S. Shlomo, A. Bonasera, S. Typel, Z. Chen, M. Huang, J. Wang, H. Zheng, S. Kowalski, C. Bottosso, M. Barbui, M. R. D. Rodrigues, K. Schmidt, D. Fabris, M. Lunardon, S. Moretto, G. Nebbia, S. Pesente, V. Rizzi, G. Viesti, M. Cinausero, G. Prete,T. Keutgen, Y. El Masri, and Z. Majka, Phys. Rev. C **85**, 064618 (2012). S. Typel, G. Röpke, T. Klähn, D. Blaschke, and H. H. Wolter, Phys. Rev. C **81**, 015803 (2010). Z. W. Zhang and L. W. Chen, Phys. Rev. C **95**, 064330 (2017).

M.B. Tsang, T.X. Liu, L. Shi, P. Danielewicz, C.K. Gelbke, X.D. Liu, W.G. Lynch, W.P. Tan, G. Verde, A. Wagner, H.S. Xu, W.A. Friedman, L. Beaulieu, B. Davin, R.T. de Souza, Y. Larochelle, T. Lefort, R. Yanez, V.E. Viola, R.J. Charity, and L.G. Sobotka, Phys. Rev. Lett. **92**, 062701 (2004). M. B. Tsang, Y. Zhang, P. Danielewicz, M. Famiano, Z. Li, W. G. Lynch, and A. W. Steiner, Phys. Rev. Lett. **102**, 122701 (2009) L.W. Chen, C.M. Ko, and B.A. Li, Phys. Rev. Lett. **94**, 032701 (2005). B.A. Li and L.W. Chen, Phys. Rev. C **72**, 064611 (2005).

P. Danielewicz and J. Lee, Nucl. Phys. A **922**, 1 (2014).

Z. Zhang and L.W. Chen, Phys. Rev. C **92**, 031301(R), (2015).

B.J. Cai and L.W. Chen, Phys. Rev. C **97**, 054322 (2018); Phys. Rev. C [**100**]{}, 024303 (2019).

R. Machleidt and D. R. Entem, Phys. Rep. **503**, 1 (2011). I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. **110**, 032504 (2013).

A. Carbone, A. Rios, and A. Polls, Phys. Rev. C **90**, 054322 (2014). G. Hagen, T. Papenbrock, A. Ekström, K. A. Wendt, G. Baardsen, S. Gandolfi, M. Hjorth-Jensen, and C. J. Horowitz, Phys. Rev. C **89**, 014319 (2014). J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. R. Schiavilla, K. E. Schmidt, and R. B. Wiringa, Rev. Mod. Phys. **87**, 1067 (2015).

U. Lombardo, in *Nuclear Methods and the Nuclear Equation of State*, Int. Rev. Nucl. Physics, Vol. 9, edited by M. Baldo (World Scientific, Singapore, 1999).

D. Dean and M. Hjorth-Jensen, Rev. Mod. Phys. **75** , 607 (2003).

K. Hebeler, A. Schwenk, and B. Friman, Phys. Lett. B **648**, 176 (2007). S. S. Zhang, L. G. Cao, U. Lombardo, E. G. Zhao, and S. G. Zhou, Phys. Rev. C **81**, 044313 (2010).

J. Chen, J. Clark, R. Dave, and V. Khodel, Nucl. Phys. A **555**, 59 (1993)

A. Fabrocini, S. Fantoni, A. Y. Illarionov, and K. E. Schmidt, Phys. Rev. Lett. **95**, 192501 (2005).

A. Schwenk, B. Friman, and G. E. Brown, Nuclear Physics A **713**, 191 (2003).

L. G. Cao, U. Lombardo, and P. Schuck, Phys. Rev. C **74**, 064301 (2006).

A. Gezerlis and J. Carlson, Phys. Rev. C **81**, 025803 (2010).

S. Gandolfi, A.Yu. Illarionov, S. Fantoni, F. Pederiva, and K.E. Schmidt, Phys. Rev. Lett. **101**, 132501 (2008).

S. Gandolfi, A.Yu. Illarionov, F. Pederiva, K.E. Schmidt, and S. Fantoni Phys. Rev. C **80**, 045802 (2009). L.P. Gorkov and T.K. Melik-Barkhudarov, Sov. Phys. JETP **13**, 1018 (1961). A. Gezerlis, C.J. Pethick, and A. Schwenk, arXiv:1406.6109 \[nucl-th\] (2014).

E. Khan, M. Grasso, and J. Margueron, Phys. Rev. C **80**, 044328 (2009).

M. Yamagami, Y. R. Shimizu, and T. Nakatsukasa, Phys. Rev. C **80**, 064301 (2009).

M. Yamagami, J. Margueron, H. Sagawa, and K. Hagino, Phys. Rev. C **86**, 034333 (2012).

J. Margueron, H. Sagawa, and K. Hagino, Phys. Rev. C **76**, 064316 (2007).

N. Chamel, Phys. Rev. C **82**, 014313 (2010).

E. Khan, J. Margueron, G. Colò, K. Hagino, and H. Sagawa, Phys. Rev. C **82**, 024322 (2010). J. Margueron, E. Khan, G. Colò, K. Hagino, and H. Sagawa, Eur. Phys. J. A **50**, 18 (2014).

J. Dobaczewski, W. Nazarewicz, and M.V. Stoitsov, Eur. Phys. J. A **15**, 21 (2002). E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A **635**, 231 (1998). L. G. Cao, U. Lombardo, C. W. Shen, and N. V. Giai, Phys. Rev. C 73, 014313 (2006).

S. Goriely, N. Chamel, and J. M. Pearson, Phys. Rev. Lett. **102**, 152503 (2009).

[^1]: <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>

[^2]: lwchen$@$sjtu.edu.cn