{'fraction_non_alphanumeric': 0.04850633675316838, 'fraction_numerical': 0.023787970227318447, 'mean_word_length': 4.062889879057925, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Goldstone Boson Emission From Nucleon Cooper Pairing in Neutron Stars and Constraints on the Higgs Portal Models

When a neutron star cools to below the critical temperature for the onset of superfluidity, nucleon pair breaking and formation (PBF) processes become the dominant mechanism for neutrino emission, while the modified URCA and the nuclear bremsstrahlung processes are suppressed. The PBF processes in neutron stars have also been used to set upper limits on the properties of axions, which are comparable to those set by supernova SN 1987A. We apply this constraint on Weinberg's Higgs portal model, in which the dark radiation particles (the Goldstone bosons) and the dark matter candidate (a Majorana fermion) interact with the Standard Model (SM) fields solely through the mixing of the SM Higgs boson and a light Higgs boson. We compare the Goldstone boson emissivity with that of the neutrinos by considering several superfluid gap models for the neutron singlet-state pairing in the neutron star inner crust, as well as in the core region. We find that the PBF processes in the superfluid neutron star interior can indeed probe Weinberg's Higgs portal model in a new parameter space region. Together with our previous works on the constraints from supernovae and gamma-ray bursts, this study demonstrates further the competitiveness and complementarity of astrophysics to laboratory particle physics experiments.


Introduction
After its birth, a neutron star with an initial core temperature of ∼ 10 11 K cools via neutrino emission from the interior during the first 10 5 years, and subsequently by photon thermal emission from the surface [1,2,3,4,5,6,7]. The basic neutrino production mechanisms are the modified URCA and the nuclear bremsstrahlung processes. The direct URCA processes are possible only if the proton fraction in the degenerate nuclear matter exceeds a certain threshold value (see e.g. Ref. [8].) Analytical approximations of neutrino luminosities for those three processes have also been derived [9,10], which can be used to simplify the comparison of theoretical modelling with the observations of neutron star cooling. However, as the neutron star interior cools to below the critical temperature for the onset of nucleon superfluidity, theoretical calculation is much complicated. Superfluidity and superconductivity have two important effects on neutron star cooling via neutrino emission. On the one hand the modified URCA emission rate is reduced due to the appearance of an energy gap at the Fermi surface, which suppresses single particle excitations of the paired nucleons. On the other hand, the nucleon pair breaking and formation (PBF) processes are switched on and become the dominant neutrino emission mechanism at this stage [11,12,13,14,15,16,17,18,19,20,21]. It is generally assumed that neutrons undergo singlet-state Cooper pairing in the neutron star inner crust, and triplet-state pairing in the core [22]. Protons are expected to undergo singlet-state pairing in the core. The NASA Chandra X-ray observatory has recorded a steady decline of the effective surface temperature of the supernova remnant Cassiopeia A. This finding is interpreted as the evidence for the existence of superfluidity in its core [23,24]. Neutron star superfluidity may also be studied using pulsar glitches, as shown by Refs. [25,26]. For recent reviews on the superfluidity in neutron stars, we refer to Refs. [4,27,28,29,30]. Beyond the Standard Model (SM) of particle physics, axion emission from nucleon PBF processes in neutron stars has been studied in Refs. [30,31,32,33]. Simulation of neutron star cooling by axions in addition to neutrinos performed in Ref. [33] determined an upper bound on the axion mass of m a (0.06 − 0.12) eV, comparable to those set by supernovae (SN 1987A) and white dwarfs. See Ref. [30] for a most recent summary of this topic. In this work we show that another interesting example to study in the context of neutron cooling is provided by Weinberg's Higgs portal model [34], which was proposed to account for the dark radiation in the early universe [35] (see, however, also Ref. [36].) In this model, Weinberg considered a global U(1) continuous symmetry associated with the conservation of some quantum number, and introduced a complex scalar field to break it spontaneously. The radial field of the complex scalar field acquires a vacuum expectation value (vev), and mixes with the SM Higgs field. The Goldstone bosons arising from the symmetry breaking would be massless and very weakly-interacting. Thus they can decouple from the early universe thermal bath at the right moment, and be a good dark radiation candidate. Previously we have examined energy losses due to the resonant emission of Weinberg's Goldstone bosons in a post-collapse supernova core [37,38], as well as in the initial fireballs of gamma-ray bursts [39]. In this work we focus on the Goldstone boson resonance production when the neutron star cools to just below the critical temperature for the onset of the superfluidity. In section 2 we briefly summarise current attempts to study neutrino emission from the nucleon PBF processes using the Green function approach in the Nambu-Gor'kov formalism. Section 3 contains a short review on Weinberg's Higgs portal model. In Section 4 we describe our method for estimating the Goldstone boson emissivity from the nucleon PBF processes. In section 5 we first determine the profiles of the neutron superfluid gap energy and the critical temperature in the neutron star inner crust and core region. Based on this information we compare the Goldstone boson and neutrino emissivities by adopting different gap models, and at different radii inside the neutron star. Variations of our neutron PBF bounds on Weinberg's Higgs portal model are discussed, and the bounds are confronted with those set by laboratory experiments and by high-energy astrophysics. In section 6 we summarise this work.

Neutrino Emissivities from the Nucleon Cooper Pair Breaking and Formation Processes
First we summarise current knowledge of the neutrino emissivities from the nucleon Cooper pair breaking and formation (PBF) processes where {NN} represents a nucleon Cooper pair, with N = n (neutron) or p (proton). The low-energy effective Lagrangian for neutral weak current interactions can be written as where G F is the Fermi constant, and the neutrino weak current is l µ =νγ µ (1 −γ 5 )ν. The nucleon weak current is with c V and c A the vector and the axial-vector coupling, respectively. For the reactions with neutrons, c V = 1 and c A = g A , where g A ≈ 1.26 is the axial coupling constant. In the non-relativistic limit, the nucleon weak vector and axial-vector currents become respectively, with σ = (σ 1 , σ 2 , σ 3 ) the Pauli matrices acting on the nucleon spinors ψ N . At finite temperature T , neutrino emissivity per neutrino flavour is calculated by means of the imaginary part of the retarded weak polarisation tensor Π µν , as (for derivations see e.g. Refs. [40,41]) Here the four-momentum transfer is q = q 1 + q 2 = (ω, q), with q i = (ω i , q i ), i = 1, 2 the four-momenta of the two emitted neutrinos, and f B (ω) = [exp(ω/T ) − 1] −1 the Bose distribution function. By inserting d 4 q δ (4) (q − q 1 − q 2 ) = 1 into the above equation and making use of Lenard's integral, the leptonic trace Tr (l µ l * ν ) is integrated to be

Green functions and dressed vertices in neutron superfluid
The retarded polarisation tensor of the Z boson, Π µν (q), is to be calculated using the dressed weak verticesΓ µ and the nucleon quasi-particle propagatorsĜ (see e.g. Ref. [42]), whereÔ denotes operators in the Nambu-Gor'kov space. In the Nambu-Gor'kov formalism [43], a two-component spinor field operator is introduced to study the particle-hole dynamics. Here c † p ↑ and c p ↑ are the creation and annihilation operators for quasi-particle of four-momentum p, and spin direction ↑. The one-particle Green's function matrix is defined aŝ for α, β = 1, 2, where T is the time-ordering symbol. The ground state (vacuum) |0 is that of the modified zero-order Hamiltonian [44] for an average number of nucleons N 0 . Here H HF = p,s ε HF (p) n p,s is the Hartree-Fock potential, where p,s n p,s = N gives the total number of nucleons. The term H φ is added for describing the pairing correlations, and µ is the chemical potential. The Pauli matricesσ i , i = 1, 2, 3, operate in the Nambu-Gor'kov space. Choosing φ p 2 = 0 as permitted by the rotational invariance in theσ space, φ p 1 is identified as the energy-gap parameter ∆ p . The quasi-particle energy is then whereǭ p ≡ ǫ p + ε HF (p) − µ, with ǫ p = | p| 2 /(2m * N ), and m * N the nucleon effective mass.
Assuming that the nuclear medium is in thermal equilibrium at temperature T , one can adopt the Matsubara Green's function technique. All energy variables associated with fermion lines are replaced by p 0 → ip n = iπ (2n + 1) T , with n = 0, ±1, ±2, .... The boson energy variables are similarly replaced by q 0 → iω m = iπ 2m T , with m = 0, ±1, ±2, .... Consider the singlet-state ( 2S+1 L J = 1 S 0 ) pairing of neutrons first. In the spectroscopic notation specifying the nucleon-nucleon scattering, S is the total spin, L the total orbital angular momentum, with L = 0, 1, 2, 3, ... ⇒ S, P, D, F, ..., and J = L + S the total angular momentum. In the frequency-momentum space, the Green's function for nucleons iŝ The diagonal elements are the normal propagators for the particle and the hole, respectively. Their relation is determined by the space and spin inversion symmetry properties of the system, and in the above equation is that in the singlet-state pairing case. The off-diagonal elements are termed the anomalous propagators, which describe the transition of a particle into a hole and a condensate pair (and vice versa) in the 1 S 0 pairing case. The Z boson polarisation tensor or self-energy is then [44,45] Π µν (ω m , q) = T pn d 3 p (2π) 3 Tr [γ µ (p n , p; p n + ω m , p + q) ×Ĝ(p n + ω m , p + q)Γ ν (p n + ω m , p + q; p n , p)Ĝ(p n , p)] , (13) where p and p + q are the fermion four-momenta in the loop. Analytic continuation of the polarisation tensor with retarded boundary condition, iω m → ω + iη, gives rise to a delta function of δ(ω 2 − 4E 2 p ) for the imaginary part of the retarded weak polarisation tensor, ImΠ µν (q). In the Nambu-Gor'kov formalism, the bare vector vertex can be written aŝ whereÎ is the identity matrix in the Nambu-Gor'kov space. The dressed vector vertex is determined by the Bethe-Salpeter equation where V is the attractive two-body potential between the nucleons. Here the Green's function for the interacting system is approximated by using the Dyson equation where the self-energy Σ(p) is a functional of the Green's function, Σ = Σ[G]. Note that the relation between the dressed vertex function and the Green's function is governed by the generalised Ward identity [43,44] so that gauge invariance is ensured.

Neutrino emissivity: vector and axial-vector current contributions in different pairing states
The mixed terms between the vector and axial-vector terms are antisymmetric, therefore they vanish after the contraction with the symmetric leptonic current tensor (q µ q ν −q 2 g µν ) in Eq. (5). It has been pointed out that at densities higher than the nuclear saturation density, triplet-state (P -wave state 3 P 2 mixed with F -wave state 3 F 2 ) pairing is more favourable [46]. Following Refs. [47,48], one usually considers three types of BCS superfluidity: singlet-state 1 S 0 , triplet-state 3 P 2 with m J = 0, and triplet-state 3 P 2 with |m J | = 2. Vector-current contribution in the case of the 1 S 0 pairing for N ν neutrino species is found by various approaches to be [14,18,49] The neutron Fermi momentum is p F = (3π 2 n n ) 1/3 , with n n the neutron number density, and the neutron Fermi velocity is defined by v F = p F /m * N . In Refs. [11,13] first estimates for Q (s) νν,V were made by using the non-relativistic approximation In the calculations therein only the temporal component of the bare vector vertices was considered, in which case the generalised Ward identity, Eq. (17), could not be satisfied. The result was an overestimation of the neutrino emissivity. Up to now different approaches agree that the early estimates [11,13] are reduced by a factor of about 0.05 v 4 F . Axial-vector current contribution in the case of singlet-state pairing was first considered in Ref. [50]. Other estimates [18,51,52] all agree with it within 25%, and can be approximated by assuming an effective nucleon mass m * N = 0.7 m N . Being suppressed by a factor ∝ v 2 F relative to v 4 F , the axial-vector current contribution thus dominates over the vector one. Neutrino emissivity due to the triplet-state ( 3 P 2 ) Cooper pairing of neutrons was first calculated in Ref. [13], and then more recently in Refs. [21,53]. It is generally assumed that the vector-current contribution is negligible, and the axial-vector current contribution is that obtained in Ref. [13] suppressed by a factor of 1/4 (for |m J | = 0) or 1/8 (for |m J | = 2). In the case of triplet-state pairing, the gap energy ∆ p depends on the direction of the quasi-particle momentum p. The expression for the 3 P 2 pairing with |m J | = 0 is [14,54,55,56] (see also Refs. [17,33]) For the case of |m J | = 2, the neutrino emissivity is Q . In the non-relativistic limit the axial-vector current represents the nucleon spin density (cf. Eq. (3).) Since the nucleon spins can fluctuate freely in the pairing state where the total spin S = 0, the axial-vector current contributions in the triplet-state pairing are not suppressed by the Fermi velocity v F . Thus they are regarded as the dominant neutron star cooling channels when the core temperature drops below the critical temperature for the onset of triplet-state pairing of neutrons, T cnt .

Weinberg's Higgs Portal Model
In this section we briefly summarise Weinberg's model [34] following the convention of Refs. [37,57]. Consider the simplest possible broken continuous symmetry, a global U(1) symmetry associated with the conservation of some quantum number W . A single complex scalar field S(x) is introduced for breaking this symmetry spontaneously. With this field added to the Standard Model (SM), the Lagrangian is where Φ is the SM Higgs doublet, µ 2 , g H , and λ are real constants, and L SM is the usual SM Lagrangian. One separates a massless Goldstone boson field α(x) and a massive radial field r(x) in S(x) by defining where the fields α(x) and r(x) are real. In the unitary gauge one sets Φ T = (0, ϕ + ϕ(x)) / √ 2, where ϕ(x) is the physical Higgs field. The Lagrangian in Eq. (21) thus becomes where the replacement α(x) → α(x)/ (2 r ) was made in order to achieve a canonical kinetic term for the α(x) field. The two fields ϕ and r mix due to the g H (S † S)(Φ † Φ) term, with their mixing angle given by where m H and m h are the masses of the two resulting physical Higgs bosons H and h, respectively. The heavier one is identified with the SM Higgs boson with mass m H = 125 GeV, and vacuum expectation value (vev) of ϕ = 246 GeV. The lighter one is assumed to have a mass in the range of MeV to hundreds of MeV. The vev of the radial field r must satisfy the perturbativity condition on the quartic self-coupling of the S field, In this model, the interaction of the Goldstone bosons with the SM fields arises entirely through the SM Higgs boson in the mixing of the ϕ and r fields. The Lagrangian for the interaction of the Goldstone bosons with the SM fields contains the Yukawa terms and the terms which couple the Goldstone boson pair with the SM Higgs boson and the light Higgs boson This model can also be extended to include a dark matter candidate by adding one Dirac field and assigning a charge U(1) W = 1 for it. After the radial field r acquires a vev, diagonalising the ψ mass matrix generates two mass eigenvalues m ± = m ψ ± f χ r for the two mass eigenstates ψ ± , which are Majorana fermions. The Lagrangian contains the term for their interactions with the SM and the light Higgs bosons. The lighter fermion is stable due to unbroken reflection symmetry, thus can play the role of the dark matter, with mass m − ≡ M χ in the GeV to TeV range. Its relic density has been calculated in Ref. [58]. Model parameters in the minimal set-up are the Higgs portal coupling g H , mass of the light Higgs boson m h , and its vev r . In the extended version the WIMP dark matter mass M χ , and its coupling to the Higgs bosons f χ are included. From the SM Higgs invisible decay width, a collider bound on the Higgs portal coupling has been derived in Ref. [57] g H < 0.011 .
Laboratory experimental limits on meson invisible decay widths have also been turned into constraints on the ϕ-r mixing angle θ H in Ref. [59]. We have scrutinised the production and propagation of Weinberg's Goldstone bosons in the post-collapse supernova core [37,38], as well as in the initial fireballs of gamma-ray bursts [39]. We found that the supernova bound based on the SN 1987A neutrino burst observations surpasses those set by laboratory experiments for all light Higgs boson masses m h 500 MeV. Using generic values for the GRB initial fireball energy, temperature T 0 2 · 10 11 K, radius, expansion rate, and baryon number density, we found that the GRB bounds are competitive to current constraints set by B meson and radiative Upsilon invisible decays in the parameter range m h /T 0 10-15. In Weinberg's Higgs portal model including the dark matter candidate, exclusion limits on the WIMP-nucleon elastic cross section set by the null results of the direct search experiments have been found to put very strong bounds on the ϕ-r mixing angle [59].

Goldstone Boson Emissivity from the Nucleon Cooper Pair Breaking and Formation Processes
The effective Lagrangian for the interactions of the Goldstone boson with the nucleons is given by The Higgs effective coupling to nucleons, f N m N / ϕ ≡ g N N H , has been calculated for the purpose of investigating the sensitivities of the dark matter direct detection experiments. So far various estimates agree that f N ≃ 0.3 (see e.g. Ref. [38] for a brief summary.) For the Goldstone boson emissivity from the neutron PBF processes we calculate by means of the imaginary part of the Higgs boson self-energy Π, as It decays dominantly to a pair of Goldstone bosons, with the decay width given by When kinematically allowed, there is also a probability for h decaying into a pair of SM fermions, as well as a pair of pions.

Resonance effect in Goldstone boson production
We evaluate the integral over the Goldstone boson momenta first: The dimensionless integral is defined by withω ≡ ω 1 /ω, and θ α is the angle between the two emitted Goldstone bosons. In the resonance region, we make use of the limit of the Poisson kernel and obtain Since this approximation is valid when m 2 h /ω 2 ≈ 2ω (1 −ω), where the latter ≤ 1, it is only applicable for ω √ 2m h and Γ h ≪ ω. On the other hand, in the large m h limit, the dimensionless integral in Eq. (36) is simply As we will show in the next section, useful constraints on Weinberg's Higgs portal model can only be obtained from the nucleon PBF processes when the Goldstone boson emission is resonantly enhanced. The typical energy of the Goldstone boson pairs emitted from the nucleon PBF processes is ω = 2E p ∼ O(1)∆ p , therefore we expect that the production of a real light Higgs boson of mass m h 15 MeV should dominate.

Suppression of scalar particle production in neutron superfluid
Assuming that the nuclear medium is in thermal equilibrium at temperature T , the Higgs boson self-energy can be calculated by adopting the Matsubara Green's function technique, Tr [γĜ(p n + ω m , p + q)Γ(p n + ω m , p + q; p n , p)Ĝ(p n , p)] .
(39) Here all energy variables associated with fermion lines are replaced by p 0 → ip n = iπ (2n + 1) T , with n = 0, ±1, ±2, .... The boson energy variables are similarly replaced by q 0 → iω m = iπ 2m T , with m = 0, ±1, ±2, .... The bare scalar vertex isγ =σ 3 in the Nambu-Gor'kov space (cf. Eq. (14).) The quasi-particle Green's functionĜ is given in Eqs. (10), (11), and (12) in the case of singlet-state pairing. Analytic continuation of the polarisation tensor with retarded boundary condition, iω m → ω + iη, gives rise to a delta function of δ(ω 2 − 4E 2 p ) for the imaginary part of the retarded scalar polarisation tensor, ImΠ(q). The dressed scalar vertexΓ in dense nuclear medium is determined by the Bethe-Salpeter equation similarly to Eq. (15) for µ = 0. However, unlike the vector vertex, the scalar vertex is not protected by the generalised Ward identity given in Eq. (17). Therefore there is a finite shift, as pointed out in e.g. Ref. [60]. In Ref. [61] the form factors for the strong interaction vertices of the scalar and the vector meson with the nucleons in dense nuclear matter (in the normal state) are calculated in dependence of the momentum transfer Q 2 and nuclear density ρ B .
Here a = s, v for the scalar and the vector meson, respectively, and F vac a (Q 2 ) are the corresponding empirically determined form factors in vacuum. For the scalar form factor F s (Q 2 , ρ B ), the reduction factor R s (Q 2 , ρ B ) contains a shift which does not vanish for momentum transfer Q 2 → 0. It is found that R s ≈ R v at normal nuclear matter density ρ B ∼ ρ 0 = 0.17 fm −1 and at momentum transfer Q 2 1 GeV 2 . Ref. [15,45,49] provides the expressions for the temporal (∝Γ 0 ) and spatial (∝ Γ i , i = 1, 2, 3) components of the neutrino vector-current emissivity from the nucleon PBF processes (cf. Eq. (18)) where , ω)). Both the temporal and the spatial components are of the order v 4 F , On the other hand, by keeping only the first term containing the temporal component in the full K V (| q|, ω) expansion, which is (| q| 2 /ω 2 )K V,0 (| q|, ω), one obtains exactly the same value for Q (s) νν,V . In the non-relativistic limit, the temporal component of the vector-current interaction can represent scalar interaction. There should be finite shifts in the temporal and the spatial components of the vector-current contribution. However, the Ward identity guarantees that both shifts cancel each other when the temporal and the spatial components are added together. Since the shift in the case we are considering is unknown, we make the assumption that in the superfluid dense nuclear medium, the scalar vertex is modified by the same amount as the temporal component of the vector vertex. Goldstone boson emissivity in the singlet-state ( 1 S 0 ) pairing is then with the suppression factor F s ∼ 0.05 v 4 F from Eq. (18). In the case of triplet-state pairing, the Goldstone boson emissivity is also suppressed, similar to the vector-current contribution for the neutrino emissivity.

Constraints on Weinberg's Higgs Portal Model
In this section we evaluate the Goldstone boson emissivity from the nucleon PBF processes, and compare with those of the neutrinos. It is necessary to know in which regions of the neutron star interior and at what temperature those processes are most efficient.

Neutron superfluid gap energies and critical temperatures
The gap energy ∆ p depends on the temperature, the quasi-particle momentum, as well as the type of superfluidity. It is well understood at densities less than one tenth of the nuclear saturation density, while significant uncertainties at higher densities arise due to complicated nature of nucleon-nucleon interactions, and difficulties in solving many-body problem. Parametrisations for the critical temperatures T c and/or the gap energies ∆ p at the neutron Fermi momentum p F can be found in Refs. [62,63,64,65], which have the general form The fit parameters ∆ 1 , k 0 , k 1 , k 2 , and k 3 are determined for various superfluid gap models in the literature. For some representative nuclear equation of states (see e.g. Ref. [66] for a recent review), the phenomenological superfluid gap models are tested against the cooling behaviour of the supernova remnant Cassiopeia A in Ref. [65]. We use the C++ package LORENE [67] with various equations of state (EoS) supplements to obtain the density profile of a non-rotating neutron star of gravitational mass M g = 1.41M ⊙ . For the DD2 equation of state [68,69,70], the circumferential equatorial radius of this neutron star is determined to be R circ = 13.24 km. Ref. [71] defines the inner crust of a cold neutron star to be the region between the density where neutrons drip out of nuclei (ρ ≃ 4 · 10 11 g/cm 3 ), and the density for the transition to homogeneous nucleonic matter at about half of the nuclear saturation density (ρ ≃ 0.5 ρ 0 ≃ 1.5 · 10 14 g/cm 3 ). Following this definition, the inner crust of this neutron star configuration is the region at the radius 9.75 km r 10.8 km.
The proton fraction decreases from Y p ≃ 0.036 at r = 9.7 km to Y p ≃ 0.016 at r ≃ 10.1 km, and then increases to to Y p ≃ 0.32 at r = 10.7 km. Fig. 1 shows the profile of the neutron singlet-state ( 1 S 0 ) pairing gap energy ∆ p (r) for the nine superfluidity gap models summarised in Ref. [65]. For a collection of more recent gap models, see e.g. Ref. [30]. The local Fermi momentum is related to the nuclear density ρ(r) by p F (r) = (3π 2 Y n ρ(r)/m * N ) 1/3 , where Y n = 1 − Y p is the neutron fraction. Most of the neutron singlet-state pairing gap models are confined to the inner crust, with three models extending into the core region. As proposed by Ref. [25], these models are useful for explaining the observed pulsar glitches. The dependence of the gap energy on the temperature T and superfluidity type is parametrised in Ref. [13] as where θ is the angle between the quasi-particle momentum p and the quantisation axis in different pairing states. The dependences are F (θ) = 1 for the 1 S 0 pairing, (1 + 3 cos 2 θ) for the 3 P 2 pairing with m J = 0, and sin 2 θ for the 3 P 2 pairing with |m J | = 2. Ref. [30] provides a detailed discussion on the temperature dependence of the gap function. For the singlet-state pairing, the high-temperature asymptotic is somewhat larger than the earlier analytical fit given in Ref. [13]. Measurement of the rapid cooling of Cassiopeia A provides the first constraints on the critical temperature for the onset of superfluidity of core neutrons and protons [23,72]. It is found that T cnt ≈ (5-9) · 10 8 K for neutron triplet-state pairing, and T cps ≈ (2-3) · 10 9 K for proton singlet-state pairing. Using observations of 18 isolated neutron stars (including Cas A), Ref. [73] found the best-fit values T cnt = 2.09 +4.37 −1.41 · 10 8 K and T cps = 7.59 +2.48 −5.81 · 10 9 K. On the other hand, the effect of the neutron singlet-state pairing on the temperature evolution occurs early on, at the neutron star age 100 yrs. Therefore the Cas A data do not provide useful constraints on T cns .

Goldstone boson emissivity and the neutron star PBF bound
We numerically evaluate the Goldstone boson emissivity from neutron singlet-state pairing, Q p F = 0.77 fm −1 , and fix the nucleon effective mass at m * N = 0.7m N . We choose a gap energy value of ∆ p = 1 MeV, which corresponds roughly to the peak of the AWP3 [74] and WAP [75] gap models considered in Ref. [65] (cf. Fig. 1.) For the critical temperature for the neutron singlet-state pairing, we adopt the approximate relation given in Ref. [65]: T cns ≈ 0.5669 ∆ p ∼ 7 · 10 9 K. In Fig. 2 we plot the Goldstone boson emissivity scaled by the Higgs portal coupling g 2 H at different temperatures below T cns , for three representative cases of the light Higgs boson mass and the vacuum expectation value (vev) of the radial field r: (m h , r ) = (1 MeV, 100 MeV), (10 MeV, 1 GeV), and (100 MeV, 1 GeV). Numerically we found in order that the Poisson kernel approximation (Eq. (37)) works well, the integration over dE p in Eq. (43) should be carried out for E p 0.55 m h instead of E p ≥ ∆ p . More explicitly, we have in the resonance region, where y ≡ E p /T , and z ≡ ∆ p /T . For m h 15 MeV, the phase space distribution functions (e Ep/T + 1) −2 render the resonance enhancement negligible. Therefore the three solid curves in Fig. 2 represent three distinct cases: for m h = 1 MeV, resonance effect can enhance the Goldstone boson emissivity in proportion to r 2 /(1 MeV) 2 in the entire temperature range shown. On the m h = 10 MeV curve, the transition from the resonance region to the large m h limit (cf. Eq. (38)) is clearly seen around T ∼ 4 · 10 9 K. On the other hand, the temperature in the neutron star superfluid inner crust is too low for the production of a real light Higgs boson of m h = 100 MeV, so the corresponding curve is completely fixed in the large m h limit. Also shown in Fig. 2 are the vector (Q In the resonance region of producing a real light Higgs boson h, where the approximation with Poisson kernel limit is applicable, we can simply scale Eq. (47) by g 2 H r 2 to satisfy this criterion. In Fig. 3 we present our main result of this work: the neutron PBF bounds on g H r for various m h , determined with i) ∆ p = 1 MeV, p F = 0.77 fm −1 , T cns = 7 · 10 9 K, at T = 6.85 · 10 9 K; ii) ∆ p = 1.75 MeV, p F = 0.77 fm −1 , T cns = 10 10 K, at T = 6.85 · 10 9 K; iia) ∆ p = 1.75 MeV, p F = 0.77 fm −1 , T cns = 10 10 K, at T = 9.8 · 10 9 K; iii) ∆ p = 1 MeV, p F = 1.37 fm −1 , T cns = 7 · 10 9 K, at T = 6.85 · 10 9 K.
In all cases the nucleon effective mass is fixed at m * N = 0.7m N , and for simplicity we make the emissivity comparison of at only one temperature in each case. Here case ii) and iia) correspond roughly to the peak of the CLS [76] and the MSH [77] neutron 1 S 0 pairing gap models considered in Ref. [65] (cf. Fig. 1), at the neutron star radius r ≃ 10.2 km. Case iii) corresponds roughly to the peak of the SCLBL [78] gap model at r ≃ 9.7 km, located in the neutron star core region.

Discussions and comparisons with laboratory and other astrophysics bounds
Comparing among the four cases i), ii), iia), and iii) listed in the previous subsection, one sees that choosing a higher gap energy ∆ p weakens the PBF bound for m h 3 MeV by less than a factor of 2. This is because for 0.55 m h < ∆ p , the increase of the Goldstone boson emissivity due to larger ∆ p is less than the increase of the neutrino emissivity. At higher temperatures, the Goldstone boson resonance production is more probable due to the larger phase space factor, especially through heavier light Higgs boson h. Thus applying the emissivity criterion in Eq. (48) at higher T improves the PBF bound for m h 3 MeV, and extends its validity  48) in the neutron star core region where the Fermi momentum p F is larger, the PBF bound can be strengthened proportionally for all m h values. This is due to the fact that the Goldstone boson emissivity is suppressed more by the Fermi velocity than that of the neutrinos: Fig. 3 are the bounds inferred from the experimental results of observing meson invisible decays, first applied to Weinberg's Higgs portal model in Refs. [59,79]. Here the upper limits on the decay branching ratio: B(B + → K + + h) < 10 −5 , B(K + → π + + h) < 10 −10 , as well as B(Υ(nS) → γ + h) < 3 · 10 −6 are used. Not shown explicitly in this plot are the collider bound of g H < 0.011 in Eq. (30), and the perturbativity condition in Eq. (25). The latter can be recast as g H r ≥ 3.08 · 10 −6 GeV g H 0.011 and one can verify that both of them are satisfied. Fig. 3 demonstrates that the neutron PBF processes in neutron stars can indeed probe Weinberg's Higgs portal model in the parameter space of low mass (m h 10 MeV), small ϕ-r mixing angle (θ H ≈ 0.0157 g H ( r /1 GeV) 10 −4 GeV), which is unaccessible to current laboratory experiments. In Fig. 4 we present a comparison of our neutron star PBF bounds with other astrophysics bounds we derived previously. The supernova (SN 1987A) bounds [37,38] were obtained by invoking Raffelt's analytical criterion on the energy loss rate per unit mass due to the emission of an exotic particle species X: ǫ X 10 19 erg·g −1 ·s −1 , at the typical proto-neutron star core temperature T = 30 MeV, and nuclear density ρ 0 = 3·10 14 g/cm 3 . In Ref. [38] we adopted two distinct estimates for the amplitudes of the nuclear bremsstrahlung processes NN → NNαα, which are the one-pion exchange (OPE) approximation, and the global fits for the nucleon-nucleon elastic scattering cross section data. On the other hand, nuclear medium effects were not included. The gamma-ray burst (GRB) bounds [39] were obtained by invoking the energy loss criterion Q e + e − →αα ∆t ′ E/(Γ 0 V 0 ). Here ∆t ′ is the time duration in the fireball comoving frame for the GRB initial fireball to expand from the initial radius R 0 to R 0 + ∆R 0 , V 0 is the initial fireball volume, and Γ 0 the initial Lorentz factor of the expanding fireball. The largest uncertainty in this consideration is the unknown initial fireball temperature, therefore we assumed two generic values T 0 = 18 MeV and 8 MeV.
In the extended version of Weinberg's Higgs portal model, latest exclusion limits published by the dark matter direct search experiments LUX [80], PandaX-II [81], and XENON1T [82] are translated into constraints on the parameter combination f χ g H r /m 2 h for WIMP mass M χ ranging from 6 GeV to 1 TeV [59]. Shown in Fig. 4 are the DM constraints for M χ = 10 GeV and 100 GeV, where the WIMP coupling f χ is fixed by requiring the relic density to be Ω χ h 2 ≃ 0.11. We conclude that in most of the parameter space of Weinberg's Higgs portal model, the supernova constraints surpass those set by laboratory experiments or by energy loss arguments in other astrophysical objects, such as the gamma-ray bursts. Nevertheless, the neutron PBF processes in the superfluid neutron star interior provide the unique possibility to explore the low mass (m h ≃ O(1) MeV), small ϕ-r mixing angle (θ H ≈ g H ( r /1 GeV) 1.6 · 10 −5 ) region. DM bounds for M χ 100 GeV remain the strongest constraints among all on the extended version of Weinberg's Higgs portal model.

Summary
Weinberg's Higgs portal model is another good example to elucidate that high-energy astrophysical objects such as the supernovae and gamma-ray bursts are excellent laboratory for probing particle physics. In this model, massless Goldstone bosons arising from the spontaneous breaking of a U(1) symmetry play the role of the dark radiation in the early universe. They couple to the Standard Model fields solely through the mixing of the ϕ and r fields, which give rise to the SM Higgs boson and a light Higgs boson h. Goldstone boson production in the hot proto-neutron star core formed in stellar collapse is dominated by the emission of a real light Higgs boson in nuclear bremsstrahlung processes and its subsequent decay. After the neutron star cools to below the critical temperature for the onset of superfluidity, resonant production of Goldstone bosons becomes possible again through the neutron Cooper pair breaking and formation processes. Theoretical calculation of the neutrino emissivity due to neutron PBF processes is a difficult task. So far all estimates by various approaches agree well, although we note that this problem is not settled yet. In this work we assume that in the superfluid phase, the scalar vertex is modified to the same extent by the nuclear medium effects as the vector vertex, and neglect the unknown shift present in the dressed scalar vertex. We compare the Goldstone boson emissivity with that of the neutrinos by considering several superfluid gap models for the neutron singlet-state pairing in the neutron star inner crust, as well as in the core region. For a typical gap energy of 1 MeV, resonance production of Goldstone boson pairs is efficient when the light Higgs boson is lighter than about 15 MeV, in which case useful constraints can be obtained. Assuming a larger gap energy increases the Goldstone boson emissivity, but does not improve the PBF constraints. In gap models which predict higher critical temperatures, the PBF bound can be improved and its applicability covers larger light Higgs boson mass m h . In those neutron singlet-state pairing gap models which extend to the neutron star core, the Goldstone boson emissivity is less suppressed due to the larger neutron Fermi momentum therein, so the PBF bounds are strengthened. Overall we found that the neutron PBF processes in the superfluid neutron star interior offer the unique possibility to explore the low light Higgs boson mass (m h ≃ O(1) MeV), small ϕ-r mixing (θ H ≈ 0.0157 g H ( r /1 GeV) 1.6 · 10 −5 ) region in the parameter space of Weinberg's Higgs portal model. Together with our previous works on supernovae and gamma-ray bursts constraints, this study demonstrates further the competitiveness and complementarity of astrophysics to laboratory particle physics experiments. case i) case ii) case iia) case iii) Figure 3: Neutron star PBF upper limits on g H r , the product of the Higgs portal coupling with the vacuum expectation value of the radial field r, for various light Higgs boson mass m h . Several gap models of the neutron singlet-state ( 1 S 0 ) pairing superfluid are considered at two different radii of the neutron star interior, as described in section 5.2. Case i): bound derived by invoking Eq. (48) at the neutron star radius r ≃ 10.2 km where the Fermi momentum is p F = 0.77 fm −1 , and assuming the gap energy ∆ p = 1 MeV, so that the corresponding critical temperature is T cns ≈ 7 · 10 9 K. The emissivity comparison is made at the crust temperature T = 6.85 · 10 9 K (red solid line). Case ii): at p F = 0.77 fm −1 , assuming ∆ p = 1.75 MeV, T cns ≈ 7 · 10 9 K, at T = 6.85 · 10 9 K (blue dashed). Case iia): same as case ii), but the emissivity comparison is made at a higher crust temperature T = 9.8 · 10 9 K (pink dotted). Case iii): same as case i), but the comparison is made at the radius r ≃ 9.7 km where p F = 1.37 fm −1 (green dot-dashed). Also shown are the upper limits set by laboratory experiments (dash-double dotted lines, from top to bottom), such as radiative Upsilon decays Υ(nS) → γ + h, B meson invisible decay B + → K + + h, as well as K meson invisible decay K + → π + + h.
[2] D. G. Yakovlev Fig. 3, here including the upper limits set by considering other astrophysical objects, as well as by dark matter detection when the extended version of Weinberg's Higgs portal model is considered. Solid lines in red, blue, pink, and green are the neutron PBF bounds derived assuming case i), ii), iia), and iii), respectively. The two dashed brown lines below the "SN 1987A" label are those we derived in Ref. [38] by applying Raffelt's criterion on the energy loss rate of an exotic particle species in the proto-neutron star core. Two distinct estimates were adopted for the amplitudes of the nuclear bremsstrahlung processes N N → N N αα: the global fits for the nucleon-nucleon elastic scattering cross section data (upper), and the one-pion exchange (OPE) approximation (lower). The two dotted orange lines with the "GRB" label between them are those we derived in Ref. [39] by invoking energy loss argument on the initial fireballs of gamma-ray bursts. Two GRB initial fireball temperature values T 0 = 18 MeV (lower) and 8 MeV (upper) were assumed, and the Higgs portal coupling g H was taken to saturate the current collider bound in Eq. (30). The two dot-dashed yellow-green lines with the "DM" label between them are the upper limits set by the dark matter direct search experiment LUX, for WIMP mass M χ = 10 GeV (upper) and 100 GeV (lower), respectively.