Tetraneutron condensation in neutron rich matter

In this work we investigate the possible condensation of tetraneutron resonant states in the lower density neutron rich gas regions inside Neutron Stars (NSs). Using a relativistic density functional approach we characterize the system containing different hadronic species including, besides tetraneutrons, nucleons and a set of light clusters ($^3$He, $\alpha$ particles, deuterium and tritium). $\sigma,\omega$ and $\rho$ mesonic fields provide the interaction in the nuclear system. We study how the tetraneutron presence could significantly impact the nucleon pairing fractions and the distribution of baryonic charge among species. For this we assume that they can be thermodynamically produced in an equilibrated medium and scan a range of coupling strengths to the mesonic fields from prescriptions based on isospin symmetry arguments. We find that tetraneutrons may appear over a range of densities belonging to the outer NS crust carrying a sizable amount of baryonic charge thus depleting the nucleon pairing fractions.


Introduction
The essential identity of nuclear forces acting between nucleons lies behind the idea of the possible existence of resonant tetraneutron states ( 4 n). Long standing searches of this four neutron state were initiated by studying its possible formation in the break-up reaction 14 Be → 10 Be+ 4 n [1]. Despite the fact that the existence of bound tetraneutrons was theoretically excluded with high confidence [2,3,4,5,6,7,8], no arguments able to rule out the resonant 4 n state were presented. A new splash of enthusiasm to studies of the four neutron resonance is motivated by a recent experimental claim that such a state with an excitation energy E 4n = 0.83 MeV and width Γ 4n = 2.6 MeV was detected in the exchange reaction α+ 8 He → 8 Be+ 4 n [9]. Intensive theoretical studies confirmed that tetraneutrons can exist in the form of resonances. In Ref. [10] and in what can be assumed to be vacuum conditions its excitation energy and width were predicted to be E 4n = 0.8 MeV and Γ 4n = 1.4 MeV, respectively. Larger values of these quantities of about E 4n =7.3 MeV and Γ 4n = 3.7 MeV were also reported in [11]. Other robust theoretical approaches predict even larger widths of tetraneutrons [2,7,12].
Such large widths seem to indicate that tetraneutrons are not likely to exist under the form of stable aggregates although this may be possible if confined by strong external fields [2,10]. The simplest system exhibiting such a tetraneutron confinement could be associated to 8 He. In such a system, nuclear forces create a significant potential barrier which prevents the break up of tetraneutrons and confines them around α-cores. Along the same line, it may be of interest to investigate whether in an extended system, such as an strongly interacting nuclear matter medium, the existence of these resonances could be due to the presence of strong mesonic fields mediating nucleonnucleon forces. Such an effect could be crucial for the description of neutron rich matter, where the probability of creation of a 4 n resonance could be, in principle, statistically favoured. If the fraction of such resonant states was high enough they could significantly affect the properties of β-equilibrated nuclear matter even despite their associated short lifetime. NSs seem to be suitable places to help shed light on the possible existence of tetraneutrons. As we will argue later in the manuscript, these particles can populate the stellar medium with an onset density that can be well below nuclear saturation density. The main question we address in this contribution is, thus, twofold. First, whether tetraneutrons could exist in a neutron rich medium and second, how they could affect its microscopic properties.
The possible existence of tetraneutrons in a nuclear medium is linked to their resonant character [3,4,5,6,7,8,12,13]. An effective treatment based on non-fundamental hadronic degrees of freedom has proven to be successful when describing the interior of NSs at intermediate energies. Therefore, in this work we will consider the appearance of tetraneutrons in β-equilibrated nuclear matter by using a relativistic field theoretical model based on, besides this new species, baryons, leptons, light clusters and mesonic degrees of freedom in a similar fashion to that done in previous works [15,16,17,18]. Further, such a treatment is also analogous to that followed for dibaryons in [23,24]. Additional work on the presence of other resonances has also been performed, for example in [25] where the effects of the mass distribution and the associated mass-dependent lifetimes of ∆(1232) resonances were partially considered, although a full in-medium treatment is still missing. Further works on this include [26,27,28]. It is important to remark that for low densities further corrections from correlations beyond a relativistic mean field (RMF) treatment are needed to precisely characterize the system. The inhomogeneous phases yield an optimized and lower free energy when compared to that of homogeneous matter. This determines what heavier clusters (nuclei) appear in the system at fixed temperature, density and proton fraction. Some approaches obtain the actual matter configuration analyzing energy functionals with contributions including that from the bulk, computed in the RMF theory [29] or Brueckner-Hartree-Fock theory [30] for uniform nuclear matter from the neutron and proton densities, surface, exchange (if needed), electron and Coulomb terms. Microscopic simulations can also provide a deeper knowledge, see [31,32]. As a first approach, the present contribution can provide a preliminary view on the species composition in a β−equilibrated gaseous system using an effective description. In this spirit we will consider a prescribed minimal coupling of these resonances to mesonic fields and solve the set of self-consistent equations arising. Within this framework, tetraneutrons, being scalar entities with non zero baryon charge, can be represented by a single-component complex field.
By using an approach based on specifically targeting resonant states [33] and taking into account the finite quantum mechanical width of these particles one could explicitly introduce this ingredient to the underlying RMF model. In particular and as a result, the calculation of the contribution of tetraneutrons to thermodynamic quantities, such as pressure or energy density, is performed in a different fashion than for other non-resonant states due to the inclusion of a mass distribution function ρ 4n (m). This approach has been applied in various versions of the hadron resonance gas model for the description of experimental hadron multiplicities produced in heavy ion collisions [34,35]. Technically, it involves an integration over masses exceeding the dominant decay channel threshold, m th 4n , while ρ 4n (m) gives an integration weight. In this contribution we set m th 4n = 4m n , with m n being the neutron mass, and use a relativistic Breit-Wigner distribution centered at m = m 4n ≡ 4m n + E 4n with an associated width Γ 4n to describe the spread of the tetraneutron mass as we discuss later.
The structure of this contribution is as follows. In Sec. 2 we present the Lagrangian model used and describe the hadronic content including bosonic and fermionic species. We introduce the formalism to treat the mass spread of the tetraneutron resonances under the approximation taken. Later, in Sec. 3, we focus on the effect of the tetraneutron condensate and discuss its impact on the pressure, density and nucleon pairing. We also explain other thermodynamical considerations when solving our set of self-consistent equations. In Sec. 4 we present our results and discuss the baryonic charge fraction carried by the tetraneutron resonance in the range of densities allowing their conden-sation. We discuss the constraints to bear in mind when considering our results with non-zero fraction of tetraneutrons, putting especial emphasis in the decay width values that can accommodate such solutions. Finally, in Sec. 5, we give our conclusions.

Lagrangian model and resonances
The model Lagrangian considered in our work includes contribution of nucleons i.e. protons (p) and neutrons (n), electrons (e), light nuclear clusters (deuterons, tritiums, 3 He nuclei and α-particles) as well as tetraneutrons in a zero temperature, electrically neutral and β−equilibrated medium. Note that the ∆(1232) resonance will be likely produced [25,28] at much higher densities (several times nuclear saturation density) than those of interest to this work and will not be considered. Thus we will restrict to the exploration of the finite-width effects only regarding the condensation of tetraneutrons.
In our modelling nuclear interaction is mediated by vector-isoscalar ω, vector-isovector ρ and scalar-isoscalar σ mesons. Baryons and stable nuclear clusters are coupled to mesons within a minimal coupling scheme which is exhaustively described in Refs. [15,16,17,18]. Due to the lack of knowledge of the tetraneutron couplings we also choose these resonances to be coupled in the same fashion. In addition, to include further corrections due to finite size of tetraneutrons we consider a spatial radial extent for these resonances up to R ∼ 5 fm as it is quite close to their reported size [13]. This value is significantly larger than ∼ 0.4 fm found for nucleons from analysis of hadron multiplicities measured in heavy ion collisions [34,35]. Therefore, although we consider the present value of R as an overestimated approximation of the in-medium spatial extent of tetraneutrons, we include it in order to account for effects of their finite size. Coupling strengths are parametrized by constants g jω , g jρ and g jσ in terms of the nucleon ones (g ω , g ρ and g σ ) as g jω = x jω g ω , g jρ = x jρ g ρ and g jσ = x jσ g σ . Let us remind here that these coupling ratios are largely uncertain in our calculation as we will discuss later. Here the index j runs over the species considered. Values of g ω = 9.479, g ρ = 8.424 and g σ = 8.487 are chosen in order to reproduce zero pressure, binding energy per nucleon equal to 16.3 MeV for the symmetric nuclear matter ground state at saturation density n 0 = 0.153 fm −3 as well as the symmetry energy coefficient a sym = 32.5 MeV. Vacuum rest masses of nucleons and stable nuclear clusters are defined as in Ref. [17] while for mesonic masses we used [37]. The mass of tetraneutrons m 4n = 4m n +E 4n is given in terms of m n and their excitation energy [9]. For electrons we set m e = 0.511 MeV. In addition, we assume no magnetic field is present. If that was the case a more detailed analysis including polarization effects [41,42] should be included.
Let us analyze the different contributions in the model Lagrangian used. First, it includes fermionic terms of nucleons, electrons, tritiums and 3 He nuclei. Here f stands for the Dirac bispinor. Contribution of deuterons is described through the Lorentz vector field d µ (we sum over repeated indexes) The Lagrangian term for α-particles is written as Finally, the contribution of tetraneutrons is described with a Lorentz scalar field φ averaged over the Breit-Wigner mass distribution with N a normalization constant. The expression for the resonant tetraneutron Lagrangian is In our prescription, covariant derivatives are defined as where I j is the isospin vector of the corresponding particle species. The isospin third component I 3 j = Q j − Bj 2 is defined through the electric Q j and baryonic B j charges. In the case of nucleons medium masses are m * j = m j − g jσ σ while for stable nuclear clusters, instead, m * j = m j + δm j − g jσ σ. In addition, in our setting the medium mass of tetraneutrons is defined as Note, that for Γ 4n = 0 the expression of m * 4n recovers the form of stable species since the mass distribution ρ 4n (m) acts as Dirac δ-function of argument m − m 4n .
It is already known that effective medium masses of stable nuclear clusters suffer a shift of their binding and excitation energies from the Pauli blocking corrections [15]. Following the formalism of Ref. [18] we write the corresponding mass shift (including tetraneutrons) as where the index k labels the aforementioned species, Z k and N k are their proton and neutron numbers, while n N and N with N = n, p represent the particle and energy densities of gas nucleons of a given sort. It is worth mentioning at this point that this correction is included in [18] for light systems starting from A = 2 (deuteron) even if it seems to be more justified for medium to heavy nuclei. However, it has proven to give a reasonable effective energy shift for light clusters and, in the same spirit, we introduce this Pauli blocking shift for the tetraneutron resonance. We will nevertheless later discuss the validity of our ansatz. Further, a few-body treatment would involve more refined calculations such as those claimed in [19] where they use a method based on an exact integral version of the Faddeev-Yakubovsky equations governing the 4-fermion system proposed by Alt, Grassberger, and Sandhas (AGS) [20] or in [21] for α-particles, however the high complexity of our many-body system prevents its application in this context. Mesonic contributions to the Lagrangian are the ones of free vector (ω and ρ) and scalar (σ) fields in addition to the well known self-interaction of σ-field [16,17,18,24,25,28,37,43] We set b = 6.284 × 10 −3 and c = −3.492 × 10 −3 so the present model fits the incompressibility factor K 0 = 250 MeV and the effective nucleon mass m * N = 0.75m N at saturation density of symmetric nuclear matter. Note that, although allowed, we do not include other possible nonlinear terms involving ω and ρ mesons in order to keep our modelling simple. Regarding the density dependence of the symmetry energy slope we find it to be L = 91.2 MeV in our model. This is somewhat larger than the global average value around L = 61 MeV and more in line with values arising from studies of nuclear masses [38] and heavy-ion collisions [39], see Fig. 24 in [40].
Isoscalar couplings of stable clusters are set as in Ref. [17]. This set up is consistent with the values of their binding and dissociation energies [15] as well as with experimental predictions of the Mott transition densities [36].
A more careful discussion involves the resonant tetraneutron states. At the moment there is neither theoretical nor experimental information about their coupling strengths. Therefore, the existing degree of uncertainty only permits performing a study in an exploratory fashion. In this spirit and following [18] we will consider the two sets shown in Table 1 where tetraneutron-meson couplings are set x 4nω /A = 1 and x 4nσ /A = 0.85 or 1 with A = 4. As we will later discuss the x 4nσ parameter is the most influential in our solution and if taken larger than the quoted values it could overpredict tetraneutron onset densities. Isovector couplings are set the same for all particle species, i.e. x jρ = 1 except electrons (x eω = x eρ = x eσ = 0). This simple parametrization in terms of fractional ratios x ij is well tested and proven to be successful in many previous works [15,16,17,18,25].
The spatial extent of tetraneutrons is probably larger than the one of nucleons and bound light nuclear clusters [13]. This gives rise to additional repulsion between these resonances, which is accounted for by the vector meson driven repulsion and Pauli blocking only in part. In addition, the finite size of tetraneutrons can lead, in principle, to their overlapping. Including an excluded volume correction can effectively prevent this and take into account meaningfully the spatial extent of tetraneutrons as we will see later. The order of magnitude estimate of the baryonic density at which tetraneutrons overlap can be obtained with their typical estimated size R ∼ 5 fm [13] and corresponding eigenvolume V 4n = 4π 3 R 3 . Using a typical fraction of baryonic charge carried by tetraneutrons 4n 4n /n B about 0.1 (see Fig. 3) and n 4n = V −1 4n we obtain n B 30 πR 3 0.08 fm −3 , which is well above the range of densities where tetraneutron may condense as we actually find in this work. Thus since these resonances are expected to exist only at small densities, the Van der Waals approximation results sufficiently accurate. A corresponding correction can be introduced to the present model within the prescription of [14]. According to it, the effective chemical potential of tetraneutrons is reduced by 4V 4n p, where p is the total pressure.

Effect of Bose condensation of Tetraneutrons
In our model, p can be written as a sum of partial pressures of non-interacting quasiparticles, with the medium masses m * j and effective chemical potentials defined through baryonic µ B and electric µ Q chemical potentials, mean values of the scalar field σ and temporal components of the ω and ρ the vector fields. In the case of tetraneutrons while for other pointlike species An excluded volume of tetraneutrons in Eq. (10) is introduced within the Van der Waals approximation. In the grand canonical ensemble this framework leads to the pressure p V dW (µ) = p 0 (µ − 4V p V dW ) defined through the eigenvolume of particles V and the pointlike pressure p 0 (µ). Formally, the excluded volume of tetraneutrons appears due to the repulsion between the hard nuclear cores. Therefore, a corresponding negative correction −4V 4n p arises as for any repulsive term. Besides, physical chemical potentials are µ j = µ B B j + µ Q Q j . This form of µ j automatically insures β-equilibrium since in this case µ n − µ p = µ e by construction. Values of µ B and µ Q are defined jointly by the value of the baryonic density n B and the requirement of electrical neutrality. Finally, the expression for the total pressure at zero temperature can be written in the form where d f = 2 is the spin degeneracy factor and the Fermi Note that formally Eq. (12) is an equation to be solved self-consistently in order to find the pressure since it enters on the righthand side indirectly through µ * 4n . The first sum runs over all fermionic species. The second one accounts for deuterons, α-particles and tetraneutrons which at zero temperature exist as the Bose-Einstein condensate (BEC) [44]. Note that this BEC is significantly delocalized since it exists in the lowest energy state and its wave function is characterized by a spatial spread that is much larger than the typical cluster size of a few fm. The real numbers ζ b represent amplitudes of zero modes in the field operators of these bosons. Physical values of ζ b are obtained by maximization of pressure, i.e. from the condition This yields either ζ b = 0 or Thus, the BEC contribution to Eq. (12) is always zero. At the same time the density of condensed bosons has a finite value if the previous equality is fulfilled. Therefore, Eq. (14) is a condition for the BEC existence. Note, that the denominator of this expression is caused by the fact that part of the system volume is excluded by non pointlike tetraneutrons. Note that the same denominator appears in expressions for densities of the rest of all species. For fermions it reads In the case of tetraneutrons it ensures that they never overlap. Indeed, even at the maximal packing of tetraneutrons reached at ζ 2 4n µ * 4n → ∞ each of them occupies a cell of volume V cell = n −1 4n = 4V 4n , which is four times larger than their eigenvolume.
The energy density can be obtained from the previous expressions considering contributions of both bosonic and fermionic species, ε = f,b µ j n j − p. Mean mesonic fields can be self-consistently found from conditions of maximal pressure ∂p ∂ω = 0, ∂p ∂ρ = 0 and ∂p ∂σ = 0. At given baryonic density they define all termodynamic quantities of electrically neutral nuclear matter at β-equilibrium. We paid special attention to the analysis of the possibility of the tetraneutron BEC existence. As mentioned before, there is a large uncertainty regarding Γ 4n values. Therefore and in order to be practical we first set this parameter equal to an average value Γ 4n 9 MeV which is close to the inverse vacuum lifetime of the tetraneutron τ ∼ 10 −22 s as reported in Ref. [9]. However, when considering the isospin asymmetric nuclear medium, τ could be, in principle, modified by the interaction with neutrons [22] whereas rare collisions with other particles could be neglected. Thus τ and, accordingly, Γ 4n are expected to suffer a possible in-medium widening with typical values belonging to the interval of Γ 4n 10 − 30 MeV at large densities. In our case, however, the reported densities for condensation are well below nuclear saturation density, see below.
On the other hand, the fraction of baryonic charge carried by the tetraneutron BEC crucially depends on the tetraneutron width. We found that for any realistic value of x 4nσ /4 ≤ 1 tetraneutrons can exist only in some limited range of densities. At the same time, for any set of coupling constants and value of the hard-core radius, the topology of the solution can be different regarding the value of Γ 4n . For example, Fig. 1 shows this fraction as a function of baryonic density calculated for set B i.e. x 4nσ /4 = 0.85, x 4nω /4 = x 4nρ = 1 and R 4n = 0 considering several values of the tetraneutron width. As it is seen from the figure, at small Γ 4n (solid curve) some baryonic densities around n B ∼ 0.0275 fm −3 support the simultaneous existence of two different values of the tetraneutron fractions, which is unphysical. This situation happens only for the tetraneutron widths smaller than some critical value Γ min .
In order to understand this behaviour let us discuss now the possible appearance of two solutions for a given density of tetraneutrons. The problem of finding n 4n can be solved in two steps. The first one corresponds to solving the field equations for σ, ω and ρ for some arbitrary amplitude of the zero mode of tetraneutrons ζ 4n . With these mean mesonic fields one can find the effective mass m * 4n and chemical potential µ * 4n of tetraneutrons. Since the amplitude of zero mode and density of tetraneutrons are related by Eq. (15), then µ * 4n and m * 4n can be considered functions of n 4n . On the second step n 4n can be found by requiring µ * 4n = m * 4n , which agrees with the condition (14) for the existence of the tetraneutron BEC. Fig.  2 illustrates the situation when at small Γ 4n the condition µ * 4n = m * 4n can be fulfilled in different manners. This shows how the two solutions for a condensate of tetraneutrons with Γ 4n = 7 MeV at n B = 0.027 fm −3 depicted in Fig. 1 arise. At the same time, the larger values of the tetraneutron width suppress these aggregates making their condensation impossible since in this case the condition µ * 4n = m * 4n can not be fulfilled by any value of n 4n . Therefore, we conclude that a physically meaningful Γ 4n should be larger than Γ min . Furthermore, in Fig. 1 large values of Γ 4n (dotted curve) make the fraction of tetraneutrons tiny, while at Γ 4n exceeding some critical value Γ max , they totally disappear. This allows us to conclude that for a single and physically correct solution with non zero fraction of tetraneutrons (dashed and dotted curves on Fig. 1) we must have Γ min ≤ Γ 4n ≤ Γ max . We show in Fig. 3 the coloured blue (pink) bands of Γ 4n providing existence of a condensate of pointlike (finite size) tetraneutrons as a function of x 4nσ /4 at x 4nω /4 = 1. This phenomenological constraint indicates the need for a care- ful determination of the in-medium width Γ 4n and its dependence on the coupling x 4nσ .

Results
We now comment on the results found in our calculation. Let us first remind that, in what follows, our analysis and calculations are performed for the two selected values of coupling parameter sets, A and B, and two values of spatial tetraneutron extent, R, and decay width, Γ 4n [45]. Since set A predicts larger baryonic densities where these resonances may exist, increased Γ 4n values are allowed (see Fig. 1). It is important to notice that in what follows we have chosen values of Γ 4n allowing a physical BEC for both sets A and B but if values are outside the interval, tetraneutrons would not be able to exist as predicted in our scenario. Values of x 4nω and x 4nσ which correspond to different onset, n os 4n , and dissolution densities, n dis 4n , of tetraneutrons are listed in Table 1.
In Fig. 4 we show the total energy per nucleon /n B − m n for set A (solid line) and set B (dashed line) as a function of baryonic density. We also depict the energy of a tetraneutron free system for reference. A pointlike treatment for tetraneutrons has been used. It can be clearly seen that for the two parameter sets used, matter with tetraneutrons is energetically favoured over that where tetraneutrons are not allowed. In other words, the presence of tetraneutrons in neutron rich matter is energetically favoured even despite the short lifetime of these resonances.
In order to better understand the role of tetraneutrons in the nuclear system under study we plot in Fig. 5 the difference of effective chemical potentials and masses of nuclear clusters in the case of pointlike (R = 0) tetraneutrons as a function of baryonic density for sets A and B. It is convenient to analyse this quantity since a nuclear species j can exist only if µ * j − m * j ≥ 0 in the case of fermions or µ * j − m * j = 0 in the case of bosons. The region of the tetraneutron BEC existence is defined by the condition µ * 4n − m * 4n = 0 appearing as horizontal line segments of dashed (set A) and solid (set B) green curves. It is clearly seen that in the case of set A all nuclear clusters have µ * j − m * j < 0 in the region of the tetraneutron BEC. In other words, tetraneutrons do not coexist with these clusters. The same situation happens in the case of set B for all clusters except α-particles, which can exist simultaneously with tetraneutrons in a narrow range of n B = 0.007 − 0.01 fm −3 . At the same time, fractions of α and 4 n in the overlap region are so small that their impact on each other is almost absent. This explains why for deuterons, tritiums, 3 He and α-particles µ * j − m * j is the same for sets A and B. As we have verified, accounting for the finite size of tetraneutrons leads to the decoupling of regions where stable nuclear clusters and tetraneutrons can coexist. This happens due to the increase of the tetraneutron onset density. As a generic conclusion we can say that the present couplings of stable nuclear clusters disfavour their coexistence with tetraneutrons.
In Fig. 6 we show the baryonic charge fraction for set A (upper panel) and set B (lower panel). We depict the different components, i.e. n (blue), p (red) and tetraneutron (green) species as a function of baryon density n B calculated for R = 0 (dashed line) and R = 5 fm (solid line). We also include the 4 n free solution (thick solid line). We have scaled n curves by a factor 1/5 and p curves by a factor of 10 in order to facilitate the reading. In our β−equilibrated system the fraction of protons remains tiny at all densities The regions of the tetraneutron BEC is shown by the horizontal line sigments of dashed (set A) and solid (set B) green curves signaling the vanishing value for the corresponding density range. Note, that 3 He, deuterium and tritium do not condensate nor coexist with tetraneutrons, while for α-particles in the case of set B there is an narrow overlapping region at nB = 0.007 − 0.01 fm −3 . and the appearance of light clusters is suppressed due to the combined effect of the negative contribution of electric chemical potential and the selected set of coupling parameters. We can see that onset (dissolution) densities of tetraneutrons are larger (smaller) when including finite volume corrections.
Note that the particle density fraction of 4 n is obtained dividing each value on the tetraneutron curves by B 4n = 4 since their particle number density is B 4n times smaller than the corresponding baryonic charge density. We can see that for the two parameter sets used in this work x 4nσ /4 1 and 4 n are restricted to the lower densites in a small fraction up to ∼ 4% of the baryon density. In the case of the depicted set A (upper panel) 4 n exist in the range n B (0.02 − 0.031) fm −3 , in line with typical Mott densities. For densities larger than those it is not so energetically favourable to gather charge into these resonances due to Pauli blocking.
To illustrate this we show in Fig. 7 the mass shift induced by Pauli blocking for the case of two sets of couplings, A and B, for pointlike tetraneutrons and for αparticles. The effect of this term in Eq. (8) clearly induces an extra energetic cost for composite species since δm > 0. It is also seen that in the case of tetraneutrons it is about two times stronger than for the case of α-particles. This is caused by the fact that most of the Pauli blocking induced shift of mass comes through neutrons, while the contribution of protons can be neglected due to their small density. Consequently, with a good accuracy we can conclude that In order to further explore the microscopic consequences of the presence of tetraneutrons we have studied the nucleon pairing into spin-zero Cooper pairs. We have selected the most attractive 1 S 0 channel using for this purpose the strategy of Ref. [46]. Note that we consider the BCS approximation although more refined treatments are indeed possible [47,48] quoting in particular those includ- ing short-range and long-range correlations to account for medium effects and polarization [49]. The dependence of the nucleon pairing gap ∆ N on its momentum k is defined by the gap equation , (17) where N (k , m * N ) = k 2 + m * 2 N and the matrix elements of the two-nucleon interaction potential V are V (k, k ) = ∞ 0 dr r 2 j 0 (kr)V (r)j 0 (k r).
Here j 0 denotes the first kind Bessel function of order zero. The consistency with the present model Lagrangian is provided by the Yukawa parametrization of V , which includes repulsive contributions from the ω and ρ mesons as well as an attractive one from the σ meson. Thus, in the coordinate space A phenomenological parameter A = 0.0435 is chosen in order to provide reasonable characteristics of this potential. Its minimum is located at r = 0.61 fm and has a depth of 50 MeV. The factor 1 2 in the ρ-meson term comes from the nucleon isospin. Note, that such a potential can be derived within the one boson exchange (OBE) approximation [50]. Safely, as it has been shown in [51] for standard BCS calculations, gap energies are quite similar for different realistic interactions and we choose this parametrization for the sake of simplicity. We believe this treatment captures the essence of the nucleon pairing in our diverse population scenario. The behaviour of the neutron pairing gap as a function of momentum is shown in Fig. 8 for model set A (upper panel) and set B (lower panel). The baryonic density is fixed at n B = 0.0275 fm −3 for set A and n B = 0.0150 fm −3 for set B, close to the maximum of the baryonic charge fraction for tetraneutrons depicted in Fig. 6. We find that the corresponding gap for protons is much smaller, in agreement with standard calculations [52] and is not shown. We depict the cases of pointlike (thin dashed line) and finite-size tetraneutrons with R = 5 fm (thin solid line). The case with no tetraneutrons is shown with thick solid line and for that case neutron Fermi momentum is k F = 0.9323 fm −1 (set A) and k F = 0.7625 fm −1 (set B). When tetraneutrons are present set B yields more pronounced differences at low momentum. This, in turn, translates into the pairing fractions as we will later see in the manuscript. We consider Γ 4n values as they appear in Table 1. We can see that the main difference arises for the low momentum when the BEC manifests more clearly the difference among paired and unpaired neutrons. When tetraneutrons are present the amplitude is somewhat decreased with respect to the case without them. In addition, the dependence on the coupling ratio x 4nσ along with the decay width is nevertheless weak leading to slight changes in the gap profile. Note that for values of Γ 4n outside the allowed bands (see Fig. 3) no solution with tetraneutrons would exist.
Pairing of nucleons is controlled by the gap ∆ N , which is strongly influenced by the tetraneutron BEC. Therefore, this condensate also significantly changes the density of paired particles. For nucleons it reads . (20) We show in Fig. 9 the behaviour of the ratios of densities of paired nucleons in presence and in absence of tetraneutrons n pair /n pair no 4n as a function of the baryonic density for model sets A (upper panel) and B (lower panel). The cases of finite-size (solid line) and pointlike tetraneutrons (dashed line) are also shown. Clearly, before the tetraneutron onset densities and beyond their dissolution ones this fraction has a unit value. The situation changes once the onset density for formation of the tetraneutron condensate is reached (see actual values in Table 1) as there is a pronounced decrease for both nucleon types. As it is seen, this reduction is larger for protons at densities belonging to the region of the outer NS crust. Although our model only considers light clusters it is worth noting that the addtional pressence of a fraction of heavier species, nuclei, is to be considered in future works as it could lead to a supression of the light bound clusters as found in [53].

Conclusions
We have explored the possible condensation of tetraneutrons in neutron rich matter inside Neutron Stars. We assumed that they can be produced in a thermodynamically equilibrated medium whose properties are controlled by the corresponding chemical potentials. As a first step, neglecting higher order correlations, we started by using a relativistic density functional approach and we find that scanning a prescribed range of couplings of tetraneutrons to the σ, ω, ρ fields based on arguments of isospin symmetry, similar to those used in the literature for other clusters, their decay width largely determines the actual presence of a tetraneutron condensate. If that was the case it can lead to a partial suppression of the S-wave nucleon pairing manifested through a reduction of the fraction of paired protons and neutrons in the system. This happens due to a more energetically favourable combination of neutrons to tetraneutrons condensing to the lowest energy state. Pauli blocking effects have been partially included using an effective treatment in the same fashion already used for stable nuclear clusters. In our model, the fraction of tetraneutrons depends on their actual decay width and, if allowed, it remains small (up to 4% of baryonic density) restricted to densities about one tenth of nuclear saturation density, thus we expect that they will have a very mild impact on the equation of state of dense nuclear matter or NS masses. We expect, however, that it could most likely affect the microscopic behaviour of the neutron rich matter in the crust. Further work is needed to clarify this latter aspect and it is left for future contributions.