Stellar limits on scalars from electron-nucleus bremsstrahlung

We revisit stellar energy-loss bounds on the Yukawa couplings $g_{\rm B,L}$ of baryophilic and leptophilic scalars $\phi$. The white-dwarf luminosity function yields $g_{\rm B}\lesssim 7 \times 10^{-13}$ and $g_{\rm L}\lesssim 4 \times 10^{-16}$, based on bremsstrahlung from ${}^{12}{\rm C}$ and ${}^{16}{\rm O}$ collisions with electrons. In models with a Higgs portal, this also implies a bound on the scalar-Higgs mixing angle $\sin \theta \lesssim 2 \times 10^{-10}$. Our new bounds apply for $m_\phi\lesssim {\rm 1~keV}$ and are among the most restrictive ones, whereas for $m_\phi\lesssim 0.5\,{\rm eV}$ long-range force measurements dominate. Besides a detailed calculation of the bremsstrahlung rate for degenerate and semi-relativistic electrons, we prove with a simple argument that non-relativistic bremsstrahlung by the heavy partner is suppressed relative to that by the light one by their squared-mass ratio. This large reduction was overlooked in previous much stronger bounds on $g_{\rm B}$. In an Appendix, we provide fitting formulas (few percent precision) for the bremsstrahlung emission of baryophilic and leptophilic scalars as well as axions for white-dwarf conditions, i.e., degenerate, semi-relativistic electrons and ion-ion correlations in the ``liquid'' phase.


Introduction
The emission of radiation through bremsstrahlung, due to the deceleration of a particle when deflected by another one, is a staple of classical and quantum field theory. The radiated power was derived by Larmor in 1897 [1], the quantum-mechanical problem was first solved by Sommerfeld in 1931 [2], and numerical solutions were analyzed in a seminal paper by Karzas and Latter in 1961 [3]. While the problem of bremsstrahlung emission could appear quaint, the topic has been revisited several times over the years. The full relativistic cross section was found only in 1969 [4]. In recent years, the numerical results have been updated [5,6], and new approximate formulae were obtained [7]. Particles other than photons can be emitted as well. In an astrophysical plasma, electron-proton bremsstrahlung can copiously produce neutrino pairs, as first proposed in Refs. [8,9]. Neutrino bremsstrahlung is the dominant neutrino emission process in stars with low temperature and high electron density [10,11]. The emission from a non-degenerate, non-relativistic plasma was obtained again in Ref. [12]. Electron-proton bremsstrahlung in the Sun provides the largest keV-range neutrino flux at Earth [13].
The neutrino-proton bremsstrahlung process ν + p → ν + p + γ has been proposed as an ambitious approach to measure neutrino masses and to distinguish their Dirac vs. Majorana nature by looking at the kinematic endpoint [14]. Moreover, putative particles beyond the standard model could interact with electrons and nucleons. The emission of photons in dark matter-proton bremsstrahlung processes has been proposed as a detection channel for sub-GeV dark matter searches [15,16]. These examples are cases of "inverted kinematics" in the sense that the collision is between a heavy and a light particle, where the energy comes from the light partner, whereas the radiation is emitted by the heavy one. We will pay careful attention to answer precisely this question: how is bremsstrahlung modified by inverting the roles of who provides the energy and who radiates in a collision.
In stellar plasmas, feebly interacting bosons instead of neutrino pairs can be produced through bremsstrahlung as well [17]. For example, axions coupling to electrons [18][19][20] are produced in ee or eN collisions. The emission rates of axions and photons are related to each other, with the former being suppressed compared to the latter by a factor O(ω 2 /m 2 e ) in the non-relativistic limit where ω ≪ m e [20]. Perhaps surprisingly, the vector and axial currents contribute at the same level to neutrino bremsstrahlung [12], so one should be careful in the parametric rescaling of the results concerning different couplings. Moreover, non-relativistic expansions should be handled with care (compare e.g. the results of Refs. [9,12] in the context of neutrinos, and [21,22] in the context of axions).
Novel CP-even bosons that couple to ordinary matter can be copiously produced in stars. Scalar production through bremsstrahlung in a non-relativistic plasma was considered in Ref. [23] and revisited in Ref. [24], where it was shown that they can be produced through resonant conversion of longitudinal plasmons. Electrophilic scalars with a coupling g ϕ ϕēe are mostly produced in this way. Nucleophilic scalars g ϕ ϕN N can be emitted by this process as well with a proton thermal loop (spectator). Parametric estimates suggested that the strongest constraint on g ϕ came from the evolution of red giants [24]. Particles coupling to both nucleons and electrons in a Higgs-portal fashion, (m f /v) sin θ ϕf f with v = 246 GeV, can emerge from this process with both protons and electrons as bystanders. The largest contribution to the scalar production would come from the electron coupling, providing a bound sin θ ≲ 3 × 10 −10 for m ϕ ≲ 1 keV based on the evolution of red giants [24].
Subsequently, other authors found that Higgs-portal scalars were mostly produced in electron-nucleus bremsstrahlung, with the nucleus radiating the particle [25,26]. The kinematics is somewhat peculiar: a light particle shakes a heavy one, which in turn emits the radiation. White-dwarf cooling then implied a bound sin θ ≲ 10 −17 [25,26], many orders of magnitude more stringent than the previous one. If true, this powerful constraint would strongly impact the parameter space of scalars mixing with the Higgs [27], such as in relaxion models [28][29][30], and suggests that white dwarfs can probe the CP-violating scalar coupling of QCD axions originating in the weak sector of the standard model [31].
Unfortunately, this amazing constraint is suspicious because the alleged bremsstrahlung rate is not suppressed by the factor (m e /m N ) 2 that one would expect for dipole radiation and that was explicitly found in the classical Larmor formula for scalars [32] 1 .
Motivated by the important consequences of this question, we revisit the emission of baryophilic and leptophilic scalars in electron-nucleus bremsstrahlung. We obtain explicit expressions for the energy-loss rate and confirm the expected (m e /m N ) 2 suppression relative to the results of Refs. [25,26]. We obtain analytical results for a non-relativistic, nondegenerate plasma and for the degenerate case, showing that bremsstrahlung is of course strongly suppressed by Pauli blocking, contrary to previous findings. Moreover, for the conditions of a white-dwarf interior, we properly treat the strong ion-ion correlations and provide simple and accurate fitting formulas.
The rest of the paper is organized as follows. In section 2 we show with simple arguments the expected scaling of baryophilic vs. leptophilic scalar bremsstrahlung radiation, i.e. how the mass of the radiating particle in a collision affects the rate. In section 3 we find the emission rates, including screening effects for different conditions. In section 4 we derive bounds on novel scalars from the white-dwarf luminosity function in analogy to earlier studies for axions. Finally, section 5 is dedicated to a summary and discussion. Several technical issues are relegated to appendices. In particular, we develop fitting formulas for the emission of baryophilic and leptophilic scalars as well as axions for white-dwarf conditions. The formulas for axions are somewhat more precise than earlier ones in the literature.

Scalar bremsstrahlung: Quantum mechanics and classical limit
In this section, we consider non-relativistic bremsstrahlung in electron-nucleus collisions using time-dependent perturbation theory in quantum mechanics and also using the classical limit. We show that the emission of scalars is perfectly analogous to that of vectors (photons) except for the different number of polarization states. In both cases, in a collision or in an atomic transition, the radiation emitted by the heavy partner (e.g. the proton in a hydrogen atom), is suppressed by an approximate factor (m e /m p ) 2 that was overlooked in previous studies of baryophilic scalar bremsstrahlung [25,26]. Our elementary reasoning supports the same finding in a detailed quantum-field theory calculation in Section 3, although in the semi-relativistic plasma of a white dwarf, there are small corrections to this simple scaling.

Photon radiation in quantum mechanics
Let us start considering the emission of photons by two interacting charged particles. This includes bremsstrahlung (free-free emission), but also free-bound or bound-bound processes, the latter equivalent to atomic transitions. Following textbook discussions (e.g. Weinberg [33]) and also a recent detailed study of quadrupole radiation [34], we consider two particles with masses m n , n = 1 or 2, and electric charges e n = Z n e with e the positive unit of electric charge, defining the fine-structure constant as α = e 2 /4π. The initial and final two-particle states |i⟩ and |f⟩ are assumed to be eigenstates with energies E i,f of the Hamiltonian where r = r 1 −r 2 is the relative coordinate between the particles and V (r) a central potential, here essentially a screened Coulomb potential. The non-relativistic interaction Hamiltonian of charged particles with photons in Coulomb gauge is −(e n /m n ) p n · A(r n ), where A(r n ) is the photon vector potential in the interaction picture at location r n of the particle n = 1 or 2. The matrix element of the interacting-particle states for the emission of a photon is therefore where ϵ is the real photon polarization vector (i.e. describing linear polarization states) and q its momentum.
Only the relative motion of the two particles, not the CM motion, can lead to radiation and so one uses the CM coordinates P = p 1 + p 2 , R = (m 1 r 1 + m 2 r 2 )/M , r = r 1 − r 2 , and p = (m 2 p 1 −m 1 p 2 )/M with M = m 1 +m 2 the total mass. One easily confirms that R and P as well as r and p fulfill canonical commutation relations, whereas r and P as well as R and p commute. The reverse mapping is p 1 = m 1 P/M + p, p 2 = m 2 P/M − p, r 1 = R + m 2 r/M , and r 2 = R − m 1 r/M . In the new canonical variables, the particle Hamiltonian is is the reduced mass. The operator sandwiched between |i⟩ and |f⟩ in Eq. (2.2) reads in the new variables Following the textbook literature (e.g. Weinberg Sec. 11.7 [33]) we notice that the factor e −iq·R introduces a recoil on the radiating system by the momentum q of the emitted radiation, which we neglect in the "long wavelength approximation" where q is much smaller than the momenta of the radiating particles. In a bremsstrahlung process, this means to neglect q in the momentum-conserving δ function. In a free-bound or bound-bound transition, it means to ignore the recoil of the final-state bound object. In the CM frame, we may also ignore the term proportional to P.
In the remaining term, we expand the exponentials up to first order, finally leading to the matrix element in the CM frame The first term represents dipole (E1) radiation, whereas the second one corresponds to quadrupole (E2) radiation as well as magnetic dipole (M1) radiation, the latter related to orbital angular momentum. We have not included possible magnetic dipoles of the charged particles that would also contribute on that order if the magnetic moment is roughly that of a Dirac fermion. Actually, a spin-flip transition can be the dominant effect as e.g. in the 21 cm hyperfine transition in hydrogen or the 14.4 keV nuclear transition in 57 Fe that has been used in solar axion searches [35][36][37][38]. However, we are primarily interested in the emission of scalars where such effects do not occur, in contrast to pseudoscalars such as axions.
For normal atomic transitions or electrons colliding with charged particles (Z, A) in a stellar plasma, we have m 1 = m e , e 1 = e, m 2 = Am N (nucleon mass m N ), and e 2 = Ze, implying e 1 /m 1 ≫ e 2 /m 2 and the reduced mass m ≃ m e so that |A 1 | ≃ |A 2 | ≃ e. For electron-electron collisions, the dipole term vanishes and the quadrupole term dominates [34]. In this case one may say that the center-of-mass and the center-of-charge coincide so that there is no time-changing electric dipole moment that is needed to emit radiation. However, the dipole term dominates unless it cancels for particles with equal e/m or unless in an atomic transition it is forbidden by the quantum numbers of the participating atomic states.
More specifically, the quadrupole term is the next order in an expansion in q·r n . That is to say, the quadrupole operator is suppressed with respect to the dipole by a factor q/p rel , where p rel is the relative momentum between the particles. In a thermal medium this ratio is of the order of T /( √ m e T ) = T /m e , where T is the temperature of the plasma. This means that for our temperatures of interest (around 1-10 keV), quadrupole processes will be suppressed by a factor T /m e ∼ 10 −2 -10 −3 . A further relative suppression derives from the ratio of A 2 /A 1 discussed below.
Concerning the dipole term, from the commutation relation between the reduced Hamiltonian H 0 = p 2 /2m + V (r) with r one finds that ⟨f|p|i⟩ = i(E f − E i )mr fi with r fi = ⟨f|r|i⟩. Therefore, we find for the dipole term where the emitted photon energy is ω = E i − E f . In this way it is obvious that the matrix element is the same independently of the magnitude of the dipole moment A 1 itself or how the two interaction partners contribute. Beyond electromagnetism, we may imagine that the electron alone carries a "leptonic charge" g L or the nucleons alone a baryonic one g B . In this case, the corresponding dipole moments are A L 1 = g L /m e or A B 1 = g B /m N with m N the nucleon mass. The ratio is A B 1 /A L 1 = (g B /g L )(m e /m N ). Therefore, apart from the obvious ratio of squared coupling constants, the baryonic emission rate is suppressed by the squared-mass ratio (m e /m N ) 2 . This insight is the main result of this discussion.
Quadrupole radiation is generically suppressed relative to dipole radiation as explained earlier. Moreover, in our exotic example, there is a factor A B 2 /A L 2 = (g B /g L )(m e /m N ) 2 so that besides the squared ratio of coupling constants, the baryonic emission rate is suppressed by the ratio (m e /m N ) 4 . Therefore, the relative suppression is even larger and we may safely neglect quadrupole radiation in all cases of interest.

Scalar radiation in analogy to photons
We now pass to consider the scalar case, which is the focus of this work. For scalars ϕ interacting with electrons or nuclei, the potential created by the radiation is simply −g ϕ. Therefore, the matrix element in the CM frame is (2.7) The first term 1 in the expansion of the exponentials does not lead to radiation because the initial and final states are orthogonal. Expanding up to second order provides For electron-nucleus collisions and for a baryonic interaction, the same hierarchy A 1 ≫ A 2 arises as in the baryonic photon case discussed earlier. The dominant dipole term can be written in the form whereq is a unit vector in the direction of the emitted radiation and ω its frequency. This is precisely the same form as for photon emission Eq. (2.6) with the replacement ϵ →q. The angular integration of the squared matrix element leads to a factor 1/3 in both cases, but in the emission of photons a factor of 2 appears for two polarization states. Otherwise the emission of a scalar or vector is the same. Therefore, considering scalar leptonic vs. baryonic emission we conclude, in analogy to the vector case, that the baryonic dipole emission rate is relatively suppressed by the factor (m e /m N ) 2 on top of the squared coupling-constant ratio.
Similar techniques have been used to relate low-energy emission processes of radiation with different spin parities. For example, the spectral axion emission from the Sun was found by similar scaling laws from the tabulated optical opacity [20], and the same can be achieved for keV-range neutrino emission [12]. Likewise, incomplete axion free-bound transition rates in the Sun or Earth were corrected using systematically the methods of non-relativistic quantum mechanics [22].

Classical limit and Larmor formula
The question of scalar radiation by the heavy partner in a binary collision has two aspects. One is to compare scalar with vector radiation and the other is the unusual kinematics, where the light partner provides the energy that can be radiated, whereas the heavy partner is the one doing the radiation. Scalar bremsstrahlung in the classical limit was discussed, for example, by Ren and Weinberg [32] who showed that the emitted power by an accelerated scalar charge moving on a prescribed trajectory is half that of the corresponding electromagnetic case, a point that also follows from our quantum-mechanical discussion in the non-relativistic limit, while Ren and Weinberg used general kinematics.
The electromagnetic power radiated by an accelerated charge in the non-relativistic limit is given by the Larmor formula [1]. For scalar emission it is then in natural units (2.10) For two particles interacting by a central potential, the mutual force is opposite equal. Because "force = mass × acceleration," the acceleration for the two partners is inversely proportional to their respective mass. If only one of them radiates, assuming we have either a leptonic charge or a baryonic one, the radiation power is inversely proportional to the squared mass of who is radiating. This is the same conclusion that we reached earlier in the quantum mechanical discussion. Based on the Larmor formula we can also estimate the energy-loss rate of a plasma. In fact, let us consider a gas of non-relativistic electrons and ions, these latter with electric charge Z e and mass m i . Let us also focus on the emission from the ions only. These experience an acceleration simply due to Coulomb interaction with electrons where b is the impact parameter. The power emitted during a single "collision" is therefore (2.12) In a plasma, rather then considering single particles, we must consider clouds of particles with number densities n i and n e , respectively for ions and electrons. If the relative velocity between the particles is v rel , then the total number of collisions per unit volume for the ions during an interaction interval ∆t will be n e n i v rel ∆t b 2π db. Given the impact parameter and the relative velocity, the typical time scale is simply set by ∆t ∼ b/v rel . The total emitted power is then where in the last step we considered b min ∼ 1/(m e v rel ) as set by the uncertainty principle. (In this "classical" argument we actually do need a vestige of quantum mechanics.) In a thermal plasma with temperature T one has v rel ∼ T /m e . Therefore we are left with the energy loss rate per unit volume i n e n i m e T . (2.14) In the next section we will see that this simple scaling with the temperature and the particle masses is indeed obtained by a rigorous quantum field theory computation.

Summary
We have studied radiation from two interacting non-relativistic particles, notably with very different masses such as electron-proton interaction. The main concern was to understand the modification between the usual situation when the light particle (the electron) carries the radiating charge (e.g. a leptonic charge) and when the heavy particle carries the radiating charge (e.g. a baryonic charge). Of course, for the usual electromagnetic case, it is the electron which mostly radiates, although the proton contributes subdominantly to the radiating dipole moment. In the relevant dipole approximation, we have found that the only modification is a factor (m e /m p ) 2 that can be easily gleaned from the classical Larmor formula or from quantum-mechanical perturbation theory. This factor was unfortunately missed in recent discussions of stellar energy losses, leading to excessively restrictive bounds on the coupling constant of new baryophilic scalars [25,26]. Still, it remains somewhat surprising that the factor (m e /m p ) 2 is the only modification, applying to free-free, free-bound or bound-bound transitions. In the latter (bremsstrahlung), the spectrum of the emitted radiation is the same in both cases because the phase-space factors are the same, including Pauli blocking in a degenerate stellar medium. It is only when relativistic modifications come in, for us in the semi-relativistic plasma of a white dwarf, that this simple scaling receives corrections as we will see in the following Section.
Other unusual cases of bremsstrahlung, where the heavy partner radiates, were studied in the recent literature, including photon emission by nuclei that are hit by a small-mass darkmatter particle [15]. The crucial point was that the bremsstrahlung spectrum extends to the maximum available energy, i.e., the kinetic energy carried by the light particle. Notice that in this situation, the center-of-mass frame of the colliding particles is nearly identical with the rest frame of the heavy one and it is at first surprising that the full kinetic energy stored in the light particle can be emitted by the heavy one. Another example is coherent neutrino scattering on nuclei, where the photon endpoint carries information of the neutrino mass because the bremsstrahlung photon can take up all the energy of the incoming neutrino [14].

Electron-proton bremsstrahlung
We now turn to a quantum-field theory calculation of scalar bremsstrahlung emission from a stellar plasma. In particular, we consider a novel CP-even scalar ϕ that interacts with protons and electrons according to L ⊃ g p ϕpp and g e ϕēe. (3.1) For the moment, we leave open if g e and g p are universal leptonic or baryonic "charges" that we used in the previous section, or if they are related to each other, for example, by a Higgs-portal interaction or if protons and neutrons carry different coupling constants. We focus on the bremsstrahlung process where p and e are respectively protons and electrons in the star of interest, with their appropriate thermal distributions. We will find that the energy-loss rate per unit volume caused by proton or electron bremsstrahlung scale as Q ϕ p /Q ϕ e = (m e /m p ) 2 as anticipated with our more elementary arguments in the previous Section, except for small corrections in a white dwarf, where the electrons are semi-relativistic.

Emission rate for general electrons conditions
The squared amplitude for the bremsstrahlung from electron-proton collisions, the latter non-relativistic, and quasi-massless scalars (m ϕ ≪ T ), is depends on the emitting particle i = e or p and the "spectator" particle j = p or e , (3.3) where Q = (ω, q) is the four-momentum of the emitted scalar, M i,1 (M j,1 ) and M i,2 (M j,2 ) are respectively the initial and final four-momenta of the particle i(j) = e or p emitting (not emitting) the scalar, either the electrons with four momenta K 1,2 = (E 1,2 , k 1,2 ) or the protons with four momenta P 1,2 . (An extension to scalars with larger masses is provided in Appendix B.) One can see from Eq. (3.3) that in the limit of non-relativistic electrons, the two amplitudes are the same up to a factor g 2 i m 2 i . This expression diverges for (M j,1 − M j,2 ) 2 → 0, but in a plasma, charged particles are subject screening effects. The most naive inclusion of this effect would be to assume an effective in-medium mass for photons, but there is no simple fundamental method for treating screening effects to a consistent order of perturbation theory. We will return to this subject later and in Appendix E, whereas for the moment we simply augment the squared matrix element with a screening factor S(∆M j ) that will be made more precise later in the context of specific assumptions about the medium.
Armed with the amplitude squared, we can compute the energy loss per unit volume due to the production of scalars from the species i in a generic stellar plasma, where m p is the proton mass, and we already integrated over the momentum of the final protons using the delta function for momentum conservation, which in the long-wavelength approximation reads p 1 + k 1 ≃ p 2 + k 2 . Moreover, we already simplified the energy delta function. In fact, energy conservation imposes where in the last step we neglected the proton kinetic energies. This approximation applies to bremsstrahlung from either protons or electrons, i.e., it is always the electron providing the emitted energy as stressed earlier in Section 2.
At this point, one should explicitly write down the squared amplitudes, but it becomes difficult to treat the emission from electrons and protons on the same footing. Here we provide our final master formulae for the two cases separately, which can be obtained after some tedious algebra. For the proton, the energy-loss rate per unit volume becomes where x 12 ≡k 1 ·k 2 is the cosine of the angle between the initial and final electrons, and we introduced the adimensional variables Notice that the argument in the structure function is simply (M j,1 − M j,2 ) 2 /m 2 e , i.e., the squared quadrimomenta exchange in the Coulomb propagator, normalized by m 2 e . We also introduced the electromagnetic fine-structure constant α = e 2 /4π and the analogous one for the scalar interaction, α p ≡ g 2 p /4π. If the scalar is emitted from the electron, we find a similar, but more cumbersome, expression, where α e ≡ g 2 e /4π, x 1 ≡k 1 ·q is the angle between the emitted scalar and the incoming electron, and x 2 ≡k 2 ·q = cos ϕ 1 − x 2 1 1 − x 2 12 + x 1 x 12 is the angle between the scalar and the outogoing electron.
Our results are easily generalized to the case in which electrons scatter off non-relativistic ions with atomic number Z and mass number A. Let us consider the more generic Lagrangian L ⊃ g p ϕpp + g n ϕnn and g e ϕēe (3.8) and the scattering process where (Z, A) is an ion with atomic number Z and mass number A. In this case Eq. (3.6) and Eq. (3.7) are easily modified introducing the "effective" couplings where the factor Z 2 comes from the ion electric charge, [g p Z + g n (A − Z)] 2 is a coherence factor to be included in the nucleophilic case because in the long wavelength approximation all the nucleons emit radiation coherently, and the factor 1/A 2 takes into account the fact that now the entire ion needs to be accelerated to emit radiation. We stress that, as long as nucleons are non-relativistic, the results in this section are exact and they can be computed for any electron chemical potential and temperature of interest. In the following we provide compact formulae which apply for different limiting cases for the electrons conditions.

Non-relativistic, non-degenerate electrons
Let us now consider the case in which electrons are non-relativistic and non-degenerate. The simple results derived here apply with good precision to the Sun and horizontal-branch stars. In such a weakly correlated plasma, Coulomb screening is well approximated by Debye screening, resulting in the static structure function where k = k 1 − k 2 and k s is the screening wave number. The screening scale receives contributions from both free electrons and ions. In the non-degenerate limit they are where n e is the electron number density and n j the number densities of ions with electric charge Z j e. A degenerate electron gas is much more "stiff" with regard to electric polarization and so they contribute much less to screening. With this screening prescription, the squared amplitude in the non-relativistic nondegenerate limit reduces to 13) and the energy loss rate per unit volume is (3.14) If electrons are non-relativistic, the occupation number is given by a Maxwell-Boltzmann distribution 15) and the energy-loss rate is We therefore find a compact expression for the scalar emission rate per unit volume, Bremsstrahlung production rate for non-relativistic and non-degenerate electrons where the coefficient is We can see that Eq. (3.17) indeed agrees with the scaling of the classical result of Eq. (2.14). One can perform a similar computation for the vector case. We find the same result of the scalar, multiplied by a factor of 2 coming from the sum over the polarizations, as expected from the quantum mechanics computation in Section 2. In Appendix D we also sketch the analogous quantum-mechanical calculation of the emission rate at second-order in perturbation theory, also in Born approximation. This calculation proceeds along the lines of Section 2, but assuming the initial and final states to be plane waves, with the Coulomb interaction included as a perturbation on the free-particle Hamiltonian.

Degenerate electrons
We now turn to the case of degenerate electrons, which will be relevant for scalar production in RG cores and WDs. Our calculations closely follow those in Ref. [39], where one of us computed the production rate of electrophilic pseudoscalars from bremsstrahlung. For completeness and comparison, we also include this case here. Assuming the electrons are degenerate we can take their momenta to be k 1,2 ∼ k F , where k F is the Fermi momentum. Furthermore, the momentum transfer between electrons can be approximated as |k 1 . With these approximations, Eqs. (3.6) and (3.7) reduce to Bremsstrahlung production rate for degenerate electrons where in each case HereŜ(x 12 ) is a function that takes care of screening as a function of scattering angle and would be unity without screening. In the non-relativistic limit (β F = 0), the integral kernels are G i = 1, whereas in general they are where we recall that x 2 = cos ϕ 1 − x 2 1 1 − x 2 12 + x 1 x 12 . These expressions are even functions in β F because the variables x 1 and x 2 vary homogeneously on the interval −1 to +1, so after integration, odd terms in β F must disappear. The full analytic expressions are given in Appendix C, but they are too complicated to be illuminating. A low-order expansion of the G-functions in β F is For a RG near helium ignition or WDs with masses of around 0.6 M ⊙ , a typical average density is 10 6 g cm −3 and the composition is either of 4 He or 12 C and 16 O, in all cases with Y e = Z/A = 1/2 electrons per baryon. In this case the Fermi momentum is k F = 409 keV and the velocity at the Fermi surface is β F = k F /(k 2 F + m 2 e ) 1/2 = 0.625 and thus β 2 F = 0.39 is not very small so that relativistic corrections are not completely negligible.
The scaling of axion emission in Eq. (3.19) agrees with the general finding that in the non-relativistic limit, the bremsstrahlung axion emission rate is 1 2 (ω/m e ) 2 that of the photon one [20], whereas we found that the scalar emission rate is 1 2 times that of photons. In other words, the non-relativistic scalar and pseudoscalar ones are the same up to a factor (ω/m e ) 2 in the latter. The overall factor (π 2 /5)(T /m e ) 2 ≃ 2.0 (T /m e ) 2 in the integrated rate reflects that in bremsstrahlung ⟨ω⟩ ≃ T is rather soft, but harder for axions than for scalars.

Screening effects for red-giant conditions
The actual emission rate strongly depends on how we deal with screening effects that are discussed in more detail in Appendix E. The degenerate electrons are "difficult to polarize" and essentially form a neutralizing homogeneous background charge density in which the nuclei (or ions) are immersed. Neglecting screening by electrons and treating the nuclei as essentially static, screening is governed by the Debye limit of the static ion-ion structure function S i (k) given in Eq. (3.11). For a single species of nuclei with charge Ze, the ionion screening scale is given by k 2 i = 4παZ 2 n i /T . Expressing the screening function in our degenerate limit it terms of the scattering angle provideŝ where ρ 6 = ρ/10 6 g cm −3 and T 8 = T /10 8 K. For a species with atomic weight A, the number density is n i = ρ/Am u with m u the atomic mass unit and we have assumed that Z/A = 1/2. In a RG core before helium ignition, T ≃ 0.7 × 10 8 K [40] and for helium Z = 2 so that κ 2 = 0.21 ≪ 1. In the non-relativistic limit (β F = 0), where all integral kernels are G i = 1, the dimensionless emission rates for all processes given by the Debye expression are

Screening effects for white-dwarf conditions
However, our main interest are WDs where the Debye screening prescription is no longer appropriate. The degree of correlation among the ions is measured by the plasma parameter Γ = Z 2 α/a i T , which is the ratio of the ion-ion Coulomb interaction energy over their thermal kinetic energy. Here a i is the ion-sphere radius given by n −1 where we have used A = 2 Z. The plasma parameter is also connected to our parameter κ 2 through Γ κ 2 = 3π 2 2 Strong correlations begin for Γ ≳ 1, corresponding to κ 2 ≳ 0.12 for 12 C, which we may call the liquid phase. For Γ ≳ 178, the ions begin to crystallize in a lattice.
For the luminosity function of low-mass WDs, we are primarily interested in the liquid phase. The static structure function must be determined numerically as discussed in more detail in Appendix E. In all cases (axions and scalars) the F i functions in Eq. (3.19) can be written as 28) and the coefficient functions are found to be well fitted by the functional form where x = log 10 (ϱ) with ρ in units of g/cm 3 . The numerical coefficients differ for different atomic charge Z and different bosons. Their values, and an extended discussion, can be found in Appendix. E.3 (see in particular Table 1).
4 Bounds from the white-dwarf luminosity function

Introduction
New low-mass particles can be systematically constrained by their emission from hot stellar plasmas, leading to observable consequences for the evolution of various well-observed stars or classes of stars. For scalar particles, the intriguing phenomena of resonant conversion from longitudinal plasmons was proposed some years ago, leading to very restrictive bounds g L < 0.7×10 −15 and g B < 1.1×10 −12 based on the brightness of the tip of the red giant (RG) branch [24]. 2 The difference between the leptonic and baryonic bounds actually represents the ratio of g L /m e and 4g B /m4 He that is now familiar from our bremsstrahlung argument. Far more restrictive bounds were derived in Refs. [25,26] using inter alia bremsstrahlung emission in WDs. Unfortunately, for baryonic scalars they did not include the generic m e /m p factor that we have argued in Sections 2 and 3. In addition, they did not use degeneracy effects correctly, another motivation to revisit the WD argument because it continues to provide one of the most restrictive limits even after these corrections.
WDs often provide very restrictive limits because, while they are about as hot as the Sun inside (around 1 keV) and have perhaps a half solar mass, they are only about the size of the Earth and therefore very dim because of their small surface, despite of being very hot ("white"). Therefore, volume particle emission competes only with a small photon surface luminosity. Moreover, they no longer burn nuclear fuel so that their evolution is a benign cooling process. Probably these points were made for the first time nearly 40 years ago in the 2 We observe that a small error has crept into the resonant emission rate Eq. (2.31) of Ref. [24] where the middle factors should read k 3 ωp ωp instead of k 2 ωp ω 2 p . We thank E. Hardy for confirming this point. However, it makes no difference for scalar bounds in the massless limit.
PhD work of one of us [41] in the context of axion emission. Since that time, many authors have studied specifically axion emission from WDs, sometimes even observing a tentative excess cooling that can manifest itself in a drift of the oscillation frequency of variable WDs and in some cases this drift can be amazingly well measured. For more details and references to the original literature we refer to a recent review by some of the original authors [42].
In our argument, we will primarily use the WD luminosity function (WDLF) in the galactic disk [43,44], i.e., the distribution of WDs as a function of luminosity. WDs are the compact remnants of low-mass stars after undergoing the red-giant and asymptotic-giant phase, after which they tend to shed their envelope in the form of a beautiful planetary nebula and remain as a glowing ember. Assuming an approximately constant birth rate, the number of WDs in each brightness interval is a direct measure of the cooling speed that can be enhanced by scalar emission. Moreover, the slope of the luminosity function would be very different compared to the case when standard surface photon emission dominates. Because scalar emission scales with 2 fewer powers of T compared with axions, the situation is now very different because axions change the shape of the luminosity function in much more subtle ways. We now venture to elaborate these general arguments in some detail.

Rescaling of axion bounds
However, before turning to a detailed analysis, we can get a rough idea of what to expect using the results obtained for axions. In fact, we have seen above that the functional forms of the emission rates are very similar, Eq. (3.19), and therefore easy to rescale.
The most stringent axion bound of g ae < 1.6 × 10 −13 at a nominal 95% C.L. was derived from the brightness of the tip of the RG branch [45]. With the scaling of Eq. (3.19) and using F e ≃ 1 and F ae ≃ 1 and T = 10 8 K = 8.6 keV, the RG axion bound translates to a bound on the scalar electron coupling of g e < 4 × 10 −15 , to be compared with a more stringent RG bound from resonant plasmon conversion of g e < 0.7 × 10 −15 [24].
For axions, a comparable bound of g ae < 2.8 × 10 −13 at a nominal 99% C.L. derives from the WDLF [46]. From the scaling observed in the previous section, we can foretell the expected sensitivity for scalars. In fact, from Figs. 4 or 6 of Ref. [46] one gleas that axion cooling gets constrained mainly by WDs with bolometric brightness M bol ∼ 7-9, where M bol = 4.74 − 2.5 log 10 (L/L ⊙ ). (4.1) This range corresponds to internal T ∼ 2-3 keV. Therefore, using again the β F → 0 limit, the bound on scalars would be g e < g ae × √ 2T /m e ∼ 10 −15 . This is only a rough estimate because scalar emission affects the WDLF in different ways compared to pseudoscalars. In particular, scalars are more important for colder and thus older WDs. In any case, this simple scaling suggests that the WDLF can provide a limit comparable to that from resonant plasmon conversion in RGs.
In both RGs and WDs we can also rescale the axion bounds to the scalar baryon case with the scaling factor (m e /m p ) 2 . For a scalar coupling to baryon number, the radiating "charge" of a nucleus is enhanced by a factor of its atomic number, but its mass receives the same factor, so indeed the scaling is (m e /m p ) 2 , independently of the chemical composition of the RGs (mostly helium) and WDs (mostly carbon and oxygen). Therefore, the estimated bremsstrahlung bounds from RGs would be g B ≲ 7 × 10 −12 , to be compared with the more stringent g B < 1.1 × 10 −12 based on resonant plasmon conversion [24]. Our naive scaled WD bound is g B ≲ 0.6 × 10 −12 , again suggesting that the WDLF can give stringent constraints.

Energy-loss rates
Here we provide for clarity the explicit energy-loss rates per unit mass for the WD case ϵ ϕ i ≡ Q ϕ i /ρ. We have checked that the degenerate approximation works extremely well and one can use directly Eq. (3.19), instead of the more cumbersome general equations. We consider an equal mixture of carbon and oxygen. Assuming g n = g p ≡ g B for the baryonic case, one has where X j is the mass fraction of the element j, and the function F i,j are given for carbon and oxygen by the fitting formulae of Eq. (E.22) and we made explicit the dependence on temperature through Eq. (E.10). The energy-loss rate per unit mass can be integrated over the entire stellar profile to yield the total scalar luminosity, The WD core can be considered isothermal. The density profile can be obtained enforcing hydrostatic equilibrium and mass conservation (see e.g. Ref. [47,48] and Section 3.5 of Ref. [49]), dP (r) where M (r) is the enclosed mass, P (r) is the pressure, and ρ(r) is the density at radius r. The two boundary conditions are used to fix e.g. the pressure at the core boundary and the central density. Finally, the equation of state can be described assuming that electrons form an ideal Fermi gas with approximately zero temperature that prevents the star from its gravitational collapse. Defining the dimensionless "relativity parameter" x ≡ k F /m e , one finds In the zero-temperature limit, k 3 F = 3π 2 Y e ρ/m u with Y e = 0.5, and the system can be solved. We will make the crude assumption that all WDs have a central density ρ c = 3.5×10 6 g/cm 3 , which corresponds to M WD = 0.607 M ⊙ , approximately the average mass of DA WDs [50][51][52], that constitute the largest population of WDs. The inclusion of General Relativity and Coulomb corrections can be neglected at our level of accuracy. We show in Fig. 1 the profile in mass coordinates, that one can compare with Fig. 1 of Ref. [53], the density profile of a 0.602 M ⊙ WD. We conclude that our approximations should give a representative profile up to perhaps a few tens percent.

WD cooling and luminosity function
WDs have no nuclear energy sources and their evolution is basically a cooling process, based on the emission of photons and neutrinos, and potentially of new particles X. The number density of WDs in a given magnitude interval is (see e.g. Eq. 2.9 of Ref. [17], and Ref. [47]) where B 3 is the (constant) birthrate normalised to 10 −3 pc −3 Gyr −1 . There are a number of assumptions needed to obtain Eq. (4.6). The relationship between the surface luminosity and the internal temperature is obtained assuming the Kramer's opacity, so that L γ = C γ L ⊙ T 7/2 , where T is the temperature in the core and C γ is determined by fitting the data of the WDLF corresponding to photon-dominated cooling. Moreover, the thermal energy is considered to be stored in the nuclei. We neglect additional effects such as physical separation processes, convection, the contribution of electrons to the specific heat, and magnetic fields [54,55]. The cooling of hot WDs is dominated by neutrino emission through plasmon decay [56]. However, once the WD is cool enough, plasmons get suppressed and cooling is dominated by photon emission from the surface. Neglecting for the time being also L X one finds Assuming an equal mixture of carbon and oxygen, and taking M = 0.6 M ⊙ and B 3 = 1, we obtain Mestel's cooling law [57] log 10 (dN/dM bol ) = 2 7 M bol − 6.84 + log 10 (B 3 ), (4.8) which provides a very good fit to data for intermediate luminosity. In this region, the luminosity can be written as L = C γ L ⊙ T 7/2 , with C γ = 8.5 × 10 −4 [58]. One can therefore obtain the core temperature by inverting this relationship, which in principle is valid only as far as cooling is dominated by photon emission. The emission of novel feebly interacting particles modifies Mestel's cooling law. As we assumed the WD core to be isothermal, the effect of scalars can be parametrized by the T dependence of the energy-loss rate, and the coupling of scalars to electrons or protons. In Fig. 2, we show the WDLF data with 3 σ error bars from Ref. [59], together with Mestel's law (thin blue), and two curves corresponding to scalar emission (thick purple) and axion emission (thick orange), parametrized as L X = L X,0 T n keV . Scalar and pseudoscalar emission rates can be parametrized respectively with n = 2 and n = 4. From Fig. 2 we see that the effect of axions is particularly pronounced at smaller M bol (hotter WDs), while scalars kick in at lower internal T (larger M bol ).  Figure 2. Comparison between the WDLF data (red dots) from Ref. [59], the simple Mestel's law (thin blue) and a WDLF in the presence of an extra cooling due to axions (orange) or scalars (purple). The scalar case departs from Mestel's law at large M bol , while the pseudoscalar case at small ones. For all curves the WD birth rate was fixed to B 3 = 1. Figure 2 gives already an idea of the maximum extra cooling allowed by data, depending on temperature dependence. In order to be more quantitative, however, we run a simple statistical test. We consider the theoretical ("th") expression for the WDLF in the case of scalar emission from the species "i" (electrons or baryons) (4.9) where we highlighted in red the two free parameters of the model. Given a WD model and an array of measured M bol , which in turn determines a temperature array, the scalar cooling is entirely determined up to a g 2 i rescaling. We then consider the experimental data from Ref. [59], shown in Fig. 2, with their associated errors bars σ(M bol ) and build a twoparameters χ 2 statistic where we use only data in the bolometric magnitude range 7.75 < M bol < 12.75. Taking this subset of data is justified for two reasons. On the one hand, at low magnitudes neutrino cooling cannot be neglected. On the other hand, for large magnitudes (very cold WDs) crystallization effects become relevant. In fact, for very old and cold WDs, the ions begin to freeze into a regular lattice structure [60]. Nevertheless, for bolometric luminosity M bol < 12.25, crystallization should not be relevant yet [44] and therefore our scalar emission rate is precise. Of course, a truly self-consistent treatment should closely follow the procedure of Ref. [59], and one should evolve WDs models which include the extra cooling process of scalar emission ab-initio. We shall perform a dedicated study in a future work, nevertheless, the present procedure should provide the correct ballpark for the excluded values. We therefore minimize Eq. (4.10) and find the exclusion limits for g i . We assumed an equal mixture of carbon and oxygen. We checked that a one-zone model with constant density ρ = 1.3×10 6 g/cm 3 and total mass M = 0.6 M ⊙ gives similar results. For a baryophilic scalar, the best fits are B 3 = 1.02 and α B = 1.15 × 10 −26 , with a reduced chi squared χ 2 red = 2.04. For a leptophilic scalar, we find a best fit for B 3 = 1.04 and α e = 5.02×10 −33 , with a reduced chi squared χ 2 red = 2.03. We find the nominal 95% C.L. limits For the corresponding Higgs portal case, we use g e = (m e /v) sin θ and find sin θ ≲ 1.9 × 10 −10 . (4.12) In Fig. 3 we show this bound in the context of that from red giants and long-range force experiments. Our WD bounds are still somewhat more restrictive than those derived from plasmon resonant conversion [24], but many orders of magnitudes weaker than those of Refs. [25,26] because of an incorrect emission rate. More precisely, their bound on sin θ is 7 orders of magnitude more restrictive, i.e., a difference of 14 orders of magnitude in the emission rate.
The main sources of discrepancy in the emission rate are the factor (m e /m p ) 2 ∼ 10 −6 and their use of nondegenerate approximations for the WD environment, which leads to another missing factor ∼ (T /E F ) 2 ∼ 10 −6 . These two factors only led to an overestimate of the scalar flux by roughly 12 orders of magnitude. Furthermore, the WD temperature and luminosity assumed for their one-zone model bound do not match at any point of Mestel's cooling law. Finally, their limit on sin θ was derived from baryonic emission, while emission from electrons should prevail. The bounds derived in this work are the strongest ones for m ϕ ≳ eV, while for smaller masses fifth-force experiments prevail [28]. Nevertheless, a truly precise comparison between RG bounds and our new WD constraints is beyond the scope of this work and requires a dedicated effort. Neither here nor in Ref. [24] a robust statistical and astrophysical analysis has been undertaken. Furthermore, Ref. [24] assumed non-degenerate and non-relativistic  Fifth-force experiment bounds (gray region) are taken from Ref. [28], while RG bounds from resonant conversion (red) from Ref. [24].
electrons. This approximation allowed them to write the scalar self-energy in a very simple form, analogous to the photon self-energy. However, a RG core at helium ignition is degenerate and semi-relativistic like our WDs, only somewhat hotter. We expect this approximation to affect the RG bound at most by a factor of a few in coupling.

Conclusions and outlook
Recent interest in the existence of putative scalar particles coupling to ordinary matter (electrons and nucleons) has prompted us to revisit stellar bounds, focusing on white dwarfs (WDs), which often offer competitive bounds on novel particles. The main production mechanism for both baryophilic and leptophilic scalars in WDs relies on electron-nucleus bremsstrahlung, motivating an explicit evaluation of the energy-loss rate due to this process. We have found the energy-loss rate for any plasma condition, and we have obtained compact expressions for the emission of scalars by a non-relativistic, non-degenerate plasma, as well as by a degenerate plasma for any degree of relativistic electron motion. While the emission of scalars from electrons is conceptually very similar to the emission of electromagnetic radiation caused by the acceleration of the light charged particle in the collision, the emission of baryophilic scalars is less trivial, and has generated some confusion in the recent literature. We have shown that, while somewhat surprising, the only modification to obtain the baryonic emission rate from the electron emission rate (in turn related to the photon emission rate) is the inclusion of the factor (m e /m p ) 2 . This result applies to free-free, free-bound or bound-bound transitions and to any degree of electron degeneracy. This simple scaling applies only in the non-relativistic limit, whereas for the semi-relativistic conditions in WDs small corrections (tens of percent) arise.
Following earlier studies of axion emission, we have obtained novel bounds from the effect that the emission of scalars has on the WD luminosity function (WDLF) in the galactic disk. We have found that the recent evaluation of these bounds from WDs were overly stringent by several orders of magnitude, the difference arising from erroneous bremsstrahlung rates. Despite this reduction of sensitivity, the WDLF continues to provide one of the most restrictive limits, slightly more restrictive than the estimated bounds from resonant conversion of longitudinal plasmons in red giants at helium ignition existing in the literature.
Besides the specific constraints derived in our paper, we identify several directions for future work. For baryophilic scalars, it appears that the WDLF provides the most restrictive limits. To substantiate them, one should perform a self-consistent evolution of WD models, including the emission of scalars, that can have non-trivial effects on the WDLF that could not be captured by our simple treatment. To take advantage of the data from the dimmer end of the WDLF, both of common WDs as well as heavier WDs that show crystallization, one needs to compute the energy-loss rate of scalars in a strongly coupled plasma. We shall perform these computation in a future work.
For leptophilic scalars, the resonant conversion of longitudinal plasmons in the core of red giants at helium ignition looks like the most powerful argument. To substantiate these results, one needs to include degeneracy effects and semi-relativistic electrons in the plasmon conversion rate. It also would be interesting to evaluate directly the impact of this emission rate on the brightness of the tip of the red giant branch as it has been done for axions and neutrino dipole moments.
Note Added: Shortly after our paper had appeared on arXiv, an independent study appeared that found similar conclusions with regard to scalar bremsstrahlung as well as to the relevance of white-dwarf cooling [61]. electrons. We can therefore define a small parameter, ϵ ∝ p 2 i /m N ≪ 1, and perform a perturbative expansion in this parameter for all the scalar products of interest. Let us also define the following quadri-momenta: Q = (ω, q) for the emitted scalar, K 1,2 = (E 1,2 , k 1,2 ) for the incoming and outgoing electron, P 1,2 = (E N 1,2 , k 1,2 ) for the nucleons. For the computation of the squared amplitude, few products will be needed and it is important to keep track of their order in the ϵ expansion. One has To obtain our main results it is enough to keep only the leading ϵ terms in these scalar products.
Let us now evaluate the relevant Feynman diagrams. Only two diagrams contribute to the simple process under consideration, one with the scalar attached to the outgoing nucleon and one with the scalar attached to the incoming nucleon. The sum of the two amplitudes reads where we have factorized the amplitude into a piece concerning electrons and one concerning the nucleons as follows: M µ e =ū e (k 2 )γ µ u e (k 1 ) , The amplitude squared and summed over spins is therefore The electron part is easily evaluated as spin M µ e M * ν e = 4(K µ the expression for the nucleon part is more cumbersome, but when contracted with the electron piece -using the scalar products defined above and keeping the lowest order in the ϵ-expansion -we obtain Eq. (3.3).

B Master formula for bosons with mass
In this appendix we generalize our results to the case of a scalar with mass m ϕ . If the radiated scalar is massless, then its energy is of the order of T , and therefore kinematically small compared to the masses and momenta of the other particles. If the scalar mass was much larger than T , then this assumption need not be true, but of course the emission would be exponentially suppressed. Therefore we still assume that the energy and momentum of the scalar are small compared to the other energies. Under this assumption, we find that the energy loss rate in Eq. (3.4) remains unchanged and therefore is We now parameterize the scalar 4-momentum as Q = ω (1, β ϕ ), with β ϕ the scalar speed, so that the master formula Eq. (3.6) for the emission of baryophilic scalars thus becomes where we recall that z i = y 2 i − 1. Likewise, the master formula Eq. (3.7) for the emission of leptophilic scalars is of similar form where in red we highlighted the difference as compared to the massless case. In terms of the variables E i ≡ y i m e we find for the scalar velocity These differences can be understood as follows. The integration lower limit on y 1 comes from the fact that E 1 = ω+E 2 ≥ m ϕ +m e , while the upper limit on y 2 from E 2 = E 1 −ω ≤ E 1 −m ϕ . Finally, from the parametrization of Q we get the scaling: where in the non-relativistic limit F(β ϕ ) ∝ β 2 ϕ , as can be seen neglecting the y 1 − y 2 term in the last factor of Eq. (B.3), thus explaining the extra factor of β 2 ϕ in Eq. (B.2). In Fig. 4 we show the ratio of the emission rates for the massive and massless cases as a function of the scalar m ϕ . For this plot we fixed the density of the medium to be ρ = 10 6 g/cm 3 , the temperature to T = 1 keV and the screening scale to k s = 1200 keV, which are typical values for the inner parts of WDs. It is evident that the production rate gets heavily suppressed as soon as m ϕ ≳ T . This is also why in the main text we limit our WDLF analysis up to masses m ϕ ∼ keV, given that T ∼ keV is the typical temperature in the WD core.

C Explicit integral kernels
The integral kernels Eqs. (3.21b) and (3.21c) can be worked out explicitly. Setting x = x 12 and β = β F for compactness of typography, they are We recall that is an odd function of its argument. The use of these explicit kernels makes the numerical evaluation of the emission rates much faster.

D Bremsstrahlung in quantum mechanics using the Born approximation
We now connect the quantum-field theoretical calculation of the emission rate in Section 3 (that uses the Born approximation) with the general quantum-mechanical argument about the mass scaling in Section 2. To this end, we here sketch the quantum-mechanical calculation of the emission rate also in Born approximation. In Section 2 we assumed the initial and final wave functions of the interacting particles to be exact solutions of the interacting system before the interaction with scalars was included. So these could have been atomic wave functions or, in the free-free case, scattering states involving Coulomb wave functions. Now, on the other hand, we assume the initial and final states to be plane waves, whereas the Coulomb interaction itself is included as a perturbation on the free-particle Hamiltonian, implying that we need to go to second-order perturbation theory. Analogous computations for photon emission are found in Ref. [62] for example. Specifically we consider electron-proton collisions and assume that the new scalar ϕ only couples to protons with a Yukawa strength g p . The emission is therefore described by the total Hamiltonian where the Coulomb and scalar-field potentials are We use rationalized units, where α = e 2 /4π. The amplitude for the process e(k 1 ) + p(p 1 ) → e(k 2 ) + p(p 2 ) + ϕ(q) arises at second order in the Born approximation and reads where the sum is over intermediate states |a⟩ with energy E a . As usual, the sum includes on matrix element for the radiation emitted before and one after the Coulomb interaction. We will now provide the different matrix elements without writing, for simplicity, the delta functions that enforce momentum and energy conservation. They will be reintroduced in the final result. The matrix elements of V C are given by the Fourier transform of the Coulomb potential and are because the exchanged momentum is always k 1 − k 2 in the long-wavelength approximation where the momentum q carried by the emitted radiation is ignored. Up to O(ω) corrections, following the computations in Section 2, the matrix elements of V ϕ are whereβ ϕ is a unit vector in the direction of motion of the emitted radiation. We now compute the propagator (E i − E a ) −1 by enforcing momentum conservation. When the Coulomb scattering occurs after the emission of the scalar as in M C,ϕ , the scalar must be included in the intermediate state energy E a , so that On the contrary, in M ϕ,C , the scattering occurs before emission, and the propagator is Putting everything together, we find Using non-relativistic phase space factors, the energy loss rate is then given by: which exactly matches the limit of non-relativistic electrons of Eq. (3.4) except for the screening correction in the Coulomb propagator.
E Screening prescription in the degenerate limit

E.1 General formulation
Particle emission from a medium is only approximately represented by individual processes among particles that interact as if they were in a vacuum. Even for simple examples such as Thomson scattering of photons on electrons in the Sun, one needs to include correlations to go beyond a rough estimate. Such correlations arise from the Pauli exclusion principle (an electron is less likely than average in the same location as another electron), but also from their Coulomb repulsion [63], an effect that is easily overlooked in this context. For the bremsstrahlung processes discussed in Section 3, correlation effects are more dramatic because without them, the rate would diverge because of the infinite-range Coulomb interaction. On the other hand, it is clear that in an electrically neutral medium, the forward-scattering rate must vanish instead of diverge. For axion (pseudoscalar) emission, this question was explicitly addressed in Refs. [17,39,[64][65][66][67] for the environments relevant in a RG core near helium ignition or in WDs, which are both electron-degenerate environments with 10 6 g cm −3 range densities and temperatures in the 10 6.5 -10 8 K range, corresponding to T ≃ few keV.
The following synopsis of this subject borrows heavily from these papers. In our cases of interest, the nuclei are heavy compared with the electron mass or energy and compared with the emitted radiation. Therefore, we can think of the nuclei as static and we are essentially considering electrons scattering on the static Coulomb field of nuclei fixed in space. In this idealized situation, there are two sources of modification of the vacuum Coulomb field of a single nucleus. One is the screening provided by the electrons themselves which are highly degenerate and therefore "difficult to polarize," whereas the other is the spatial correlation of the nuclei caused by their Coulomb repulsion. At low T , they actually arrange themselves in a body-centered cubic lattice and then are strongly correlated, not located independently at random relative positions. Strong correlations are more important in heavier WDs that have larger densities. The crystallization process was recently observed in the WD luminosity function of 0.9-1.1 M ⊙ [68], but plays no strong role for 0.6 M ⊙ WDs, let alone in RG cores. We will not have to worry about outright crystallization, yet we will have to worry about going to the intermediate correlation regime ("liquid phase") beyond Debye screening.
In this Appendix we use the notation that the electron momenta in the scattering process are k 1 and k 2 , whereas their momentum transfer is k = k 1 − k 2 . For very degenerate electrons, those able to scatter are at the Fermi surface with |k 1 | = |k 2 | = k F , the latter being the Fermi momentum. Therefore, as in the main text, k 2 = 2k 2 F (1 − x 12 ), where x 12 is the cosine of the angle between the in-and outgoing electron.
The deformation of a fully degenerate and homogeneous electron gas by an external test charge is governed by the Thomas-Fermi (TF) wave number that is 6 , (E.1) where ρ 6 is ρ in units of 10 6 g cm −3 and we have assumed that Y e = Z/a = 1 2 electrons per baryon as will be the case in a medium consisting of 4 He, 12 C, or 16 O. We also recall that The TF scale provides a Yukawa modification of the Coulomb potential or equivalently, in Fourier space, the squared Coulomb propagator gets modified as |k| −4 → (k 2 + k 2 TF ) −2 . Therefore, in a Coulomb integral, we should include if degenerate electrons were the only source of screening. The quasi-static nuclei (usually called ions in this context) require a different treatment. The electron scattering amplitudes from the ensemble of nuclei interfere coherently, where the interference term would average to zero if the targets were at random locations. Otherwise, the scattering process requires a "static structure factor" S i (k) of the momentum transfer k, where the index i stands for "ion." Therefore, in a Coulomb integral, we should include Averaging over directions of k, the structure factor is only a function of |k|. The limiting behavior of S i (|k|) is 1 for large |k| and it behaves as |k| 2 for small |k|.
Notice that the squared matrix elements in Section 3 involve a squared Coulomb propagator 1/|k| 4 . However, the phase-space integration over momentum transfers d 3 k = 4π d|k| k 2 introduces a factor k 2 in the numerator. So without screening, the rates would have a simple 1/k 2 divergence that is logarithmic in the integrated rate, even though this may not be directly apparent from the expressions in Eqs. (3.21a)-(3.21c). Therefore, the k 2 scaling at low |k| of the structure function is enough to moderate the divergence. As expected in a neutral medium, Coulomb scattering processes vanish in the forward direction (i.e. for vanishing momentum transfer). Unlike the TF prescription, that behaves as |k| 4 , one here does not modify the Coulomb field with a Yukawa factor and the resulting rates are not those that one would obtain from a screened Coulomb field, but we still refer to this modification as a screening effect.
In a weakly correlated medium (sufficiently large T ), the ion structure factor is given by the Debye formula The mass density is ρ = n i Am u with A the atomic mass number (assuming only a single species) and m u = 0.931 GeV the atomic mass unit. Therefore, Z 2 n i = Z(Z/A) ρ/m u = (Z/2) ρ/m u because in our media of interest, Y e = Z/A = 1/2. Therefore, numerically where T 8 = T /10 8 K and Z 2 = Z/2, corresponding to 4 He, i.e., the reference conditions roughly correspond to a RG near helium ignition. (Recall that 10 8 K = 8.6 keV.) A dimensionless parameter that we often use in the main text is again for Z/A = 1/2 and the numerical factor corresponds roughly to a RG core. As anticipated, the screening scale from degenerate electrons is much smaller compared to ion correlations. One way of including both effects is the prescription [39] 1 In Ref. [65] we can see for example in their Eqs. (12) or (21) that they include electron screening with Jancovici's static dielectric function ϵ(k, 0) that is given in their Eq. (3).
(Notice that they define the Thomas-Fermi scale with the non-relativistic formula where E F = m e .) So they use [1/ϵ(k, 0)] 2 where we use S TF (k). We have checked that the two expressions differ from each other only by a k-dependent factor of the order of α/π and therefore agree on the relevant level of perturbation theory.

E.2 Strongly correlated plasma
Beyond the Debye approximation, our picture is that of mobile ions immersed in an "infinitely stiff" electron background, i.e., a homogeneous neutralizing charge density. This is the traditional picture of a correlated one-component plasma. The usual "plasma parameter" to measure the strength of the correlations is the ratio of the ion-ion Coulomb interaction energy over their thermal kinetic energy, Γ = Z 2 α/a i T , where a i is the ion-sphere radius given by n −1 Therefore, up to a factor, both quantities convey the same information. It is also worth noting that k 2 i a 2 i = 3Γ and therefore the Debye structure factor can be expressed as .
Overall we conclude that the Debye prescription is certainly good enough for RGs near helium ignition, whereas in WDs, especially toward the colder end of the luminosity function, the approximation is not necessarily sufficient. In the numerical simulations of axion emission in WDs [46], the used emission rate actually included the screening prescription of Ref. [65] and therefore took account of strong correlations for Γ > 1.
In this regime one may use tabulated values for S i (ak), where here ak = |a i k| is a dimensionless momentum transfer in terms of the ion-sphere radius. Tabulations for a onecomponent plasma are found in Ref. [64] for Γ = 1, 3, 6, 10, 20, 40, 80, 100, 125, and 160, whereas in Ref. [67] for Γ = 2, 5, 10, 20, 40, 80, 125, and 160, both going back to Ichimaru and collaborators. In the regions of overlap, both sets are practically identical. In Fig. 5 we show the numerical results (solid lines) from these data, i.e., S i (ak) and compare them with the Debye approximation (dashed lines). We show the results for Γ = 1, 2, 5, 10, 20 and 40, from upper left to lower right, with colors blue, orange, green, and so forth.
To estimate quantitatively the impact of ion-ion correlations, we may imagine that the rates are expanded in powers of β F , where actually only even powers of β F appear. In the absence of screening effects, we need integrals of the form of Eq. (3.20), i.e., we need . (E.13) Here, n = 0 is the simple Coulomb integral that we need for the β F = 0 limit, wheres we need also n = 1 if we go to the next order β 2 F . To include the static structure function, we observe that the momentum transfer |k| ranges from 0 to 2k F and that S i (ak) is provided in terms of the dimensionless momentum transfer ak = |a i k|. Therefore, under the Coulomb integrals we must include the factor S i a i k F 2(1 − x 12 ) , (E.14) where a i k F = 9πZ 4  [64,67]. From upper left to lower right, following the colors from blue, orange, green, and so forth, these correspond to plasma parameters Γ = 1, 2, 5, 10, 20, and 40. The dashed lines are the corresponding Debye structure functions of Eq. (E.12) which asymptotically approach 1 for large ak, where screening is irrelevant, and agree with the full structure functions at small ak. The dimensionless momentum transfer is defined as ak = |a i k| in terms of the ion-sphere radius a i defined in Eq. (E.9).
where we have chosen the ion charge 6 as a reference value for 12 C in a WD. To include Thomas-Fermi screening, we need the further factor Therefore, overall the Coulomb integrals are . (E.18) In the Debye limit (Γ ≪ 1), and ignoring the TF term (κ 2 TF = 0), the first two integrals are The integrals can also be done including the TF term, but the expressions are too complicated to be illuminating.  In Fig. 6 we show the Coulomb integrals F 0 (Γ) and F 1 (Γ), assuming Z = 6 (carbon) and a density of ϱ = 10 6 g cm −3 . The blue dots are the numerical integrals for the tabulated static structure functions and also include the TF screening by degenerate electrons. We also show the Debye result, and in the upper panel also the result in the absence of TF screening (dashed lines). We also show the fit functions (blue lines) Itoh et al. [64] have provided analytic fit functions, where our F 0 is what they call 2⟨S −1 ⟩ and F 1 is what they call 4⟨S +1 ⟩. For F 0 , their fit function virtually overlays with ours and the agreement is very good. For F 1 we show their fit function as a green line. The agreement is slightly worse, but still, as F 1 would appear together with a factor β 2 F , the overall error would be small.

E.3 Full numerical integration for carbon and oxygen
Clearly the emission rates will be significantly larger in a WD than predicted by a naive application of Debye screening. Moreover, the non-relativistic expansion using only the β F = 0 limit is somewhat rough for WD conditions. Therefore, we consider the full expressions for baryophilic scalar, leptophilic scalar, and axion emission through their electron coupling. Our data can be represented by a fit function of the form F fit (ρ, Γ) = A(ρ) Γ −0.37 + B(ρ) Γ +0.03 .
(E. 22) For all cases, the coefficient functions are found to be well fitted by the functional form where x = log 10 (ϱ) with ρ in units of g/cm 3 . The numerical coefficients are different for different atomic charge Z and different bosons ( Table 1). The fit applies to the range 4 ≤ x ≤ 7 and 1 ≤ Γ ≤ 160, where it is typically good at the few % level. For x ≳ 5, the fit works better than 1%. In order to be concrete and to illustrate the quality of our fitting formulas, we consider once more WD conditions with the density ρ = 10 6 g cm −3 and Z = 6 (carbon). In this case, the various physical parameters are: Fermi momentum k F = 409 keV. Thomas-Fermi wave number: k TF = 51.0 keV. Velocity at Fermi surface: β F = 0.625. Screening scale from ions: k i = 1216 keV/T 7 , where T 7 = T /10 7 K. Ion-sphere radius: a −1 i = 117 keV. Plasma parameter: Γ = 35.8/T 7 . In Fig. 7 we show the emission rates for our three generic scalar boson models.  Figure 7. Energy-loss rate for ρ = 10 6 g cm −3 and Z = 6 (carbon) for different bosons. The curve "Debye" given by Eq. (E.19) is the same in all panels and, apart from global factors, is the emission rate for β F = 0 and includes only ion-ion correlations in the Debye limit. The data points come from numerical integrations with the full ion-ion correlation function and include Thomas-Fermi screening by the electrons. The blue lines are the fit function of Eq. (E.23). For axions, the green curve is the fit function of Nakagawa et al. [66] that is claimed to be accurate to better than 20%, but much better for this example.
In every panel, we show the "naive Debye" rate as an orange line. This is simply the Coulomb integral F 0 with the inclusion of only the static ion-ion correlation. The full numerical integration, including the Thomas-Fermi screening, for the available tabulations of the ion-ion correlation are shown as blue dots. The results of the fitting formulas Eq. (E.22) are shown as blue lines. For axions, analytic fitting formulas were already provided by Nakagawa et al. [66] (green line), which agree with our results at their claimed level of accuracy of better than 20%.