Stellar probes of dark sector-photon interactions

Electromagnetically neutral dark sector particles may directly couple to the photon through higher dimensional effective operators. Considering electric and magnetic dipole moment, anapole moment, and charge radius interactions, we derive constraints from stellar energy loss in the Sun, horizontal branch and red giant stars, as well as from cooling of the proto-neutron star of SN1987A. We provide the exact formula for in-medium photon-mediated pair production to leading order in the dark coupling, and compute the energy loss rates explicitly for the most important processes, including a careful discussion on resonances and potential double counting between the processes. Stringent limits for dark states with masses below $3\,$keV ($40\,$MeV) arise from red giant stars (SN1987A), implying an effective lower mass-scale of approximately $10^9\,$GeV ($10^7\,$GeV) for mass-dimension five, and $100\,$GeV ($2.5\,$TeV) for mass-dimension six operators as long as dark states stream freely; for the proto-neutron star, the trapping of dark states is also evaluated. Together with direct limits previously derived by us in Chu et al. (2018), this provides the first comprehensive overview of the viability of effective electromagnetic dark-state interactions below the GeV mass-scale.


I. INTRODUCTION
The prospect that new physics might be hiding under our noses in form of light dark states that have been in kinematic reach for decades is most intriguing if not seemingly preposterous. In fact, cases exist where new interactions are of comparable strength to the ones encountered in the Standard Model (SM), while being compatible with all to-date searches. The direct test of such physics, i.e. new particles and interactions below the GeV-scale has become a major field in recent years [1,2], and provides a complementary direction to the beyond-SM searches at the energy frontier.
Whereas the GeV-mass scale might comprise somewhat of a "blind-spot" that allows for the existence of new physics with appreciable interactions to the SM, once the mass enters the keV-regime the landscape changes fundamentally. Astrophysical constraints on long-lived dark states that are derived from stellar cooling arguments [3] are typically so severe that the cases for laboratory detection drastically diminish. Of course, the observable signatures of dark states depend on the nature of the coupling to the SM. For example, a new force can be mediated by new scalar or vector particles. Benchmark models are then derived based on minimality of the SM extension and on the dimensionality of the interaction operator, and within this framework the interplay between cosmological and astrophysical implications and direct tests is fleshed out.
In this work we will consider-from the low-energy effective theory point of view-an even more minimal possibility than the dark photon, namely, that the SM photon is the new physics mediator. Beyond carrying a millicharge, DM may also interact directly with the photon through a number of higher dimensional operators that encapsulate magnetic or electric dipole moment interactions (MDM or EDM), an anapole moment (AM) or a charge radius interaction (CR). These possibilities were originally considered in [21][22][23] with further studies on the phenomenology found in [24][25][26][27][28][29]. Motivated by the intense efforts to search for sub-GeV dark sector states [1,2], the topic of form-factor interactions was recently revisited in detail by some of us [30].
In [30] we focused on prospects of detecting electromagnetic (EM) form factor interactions of a dark sector Dirac particle χ with mass at or below the GeV-scale. The direct production of pairs χχ was constrained with data from BaBar [31], NA64 [32] and mQ [33] and future improvements in sensitivity were derived for Belle-II [34], LDMX [35] and BDX [36]. The direct sensitivity was then compared with indirect probes such as electroweak precision tests, flavor physics constraints, as well as with results from LEP and LHC. It was found that, owing to the higher dimensionality of the operators, high energy probes provide superior sensitivity. The conclusions are independent of the lifetime of χ, as long as is stability is guaranteed while traversing terrestrial detectors.
In contrast, if χ is long-lived, additional constraints from cosmology, astrophysics, and direct DM searches apply, and in [30] we have considered the most impor-tant ones that are crucial in the MeV-GeV mass bracket of χ. However, once we allow the χ-mass to drop into the keV-region, additional constraints from the production of χχ-pairs in stars become important [3]. In this work we complement our previous results derived in [30] with astrophysical limits that apply once the dark state is stable on a macroscopic time scale, without necessarily demanding that that sub-MeV χ particles make up the DM. We derive the limits from stellar cooling that arise from red giant (RG), horizontal branch (HB) stars, and the Sun, and revisit our calculation of the supernova bound, taking into account all major production channels.
Concretely, we are considering the following fundamental dark state emission processes, highlighting in brackets the stellar system(s) for which the process is most relevant, Plasmon decay: Annihilation: Bremsstrahlung: e − N → e − N χχ (RG, HB, Sun), Compton scattering: e − γ T,L → e − χχ (all).
The respective processes are decay of in-medium longitudinal (L) and transverse (T) modes of thermal photons γ T,L which we will simply refer to as "plasmons", electron-positron annihilation, electron bremsstrahlung on protons and nuclei, nucleon-nucleon bremsstrahlung and Compton scattering with the emission of a χχ-pair. Exemplary respective diagrams are shown in Fig. 1. The paper is organized as follows: in Sec. II we first set the stage by listing the effective operators that mediate χ-photon interactions. Section III gives a brief account on stellar energy loss arguments. Our calculations on χ particle emission are presented in Sec. IV. The ensuing constraints are then collected in Sec. V before concluding in Sec. VI. Several appendices provide details on the calculations and are referenced in the main text.

II. ELECTROMAGNETIC FORM FACTOR INTERACTIONS
A Dirac fermion χ may have a number of interactions with the photon gauge field A µ or its field strength tensor F µν . At mass dimension-5 the interaction terms of the Lagrangian are given by where µ χ and d χ are the MDM and EDM coupling which may be measured in units of the Bohr magneton, µ B ≡ e/(2m e ) = 1.93 × 10 −11 e cm; m e is the mass of the electron and σ µν = i 2 [γ µ , γ ν ]. At mass dimension-6 we have where a χ and b χ are the AM and CR coefficients. All coupling strengths in Eqs. (5) and (6) are real. At massdimension-7 the interactions involve two photons at the vertex and hence require a dedicated treatment. For this reason we restrict our study to dim-5 and dim-6 operators.
The effective interactions in Eqs. (5) and (6) may, e.g., arise from the compositness of χ [37][38][39] or perturbatively, from a UV completion that contains electrically charged states [40]. In the latter case, MDM and EDM moments are e.g. generated by loop-induced axial or vector Yukawa interactions y A,V of χ with additional scalars and fermions. Parametrically, one expects where Q is the electric charge of the mediator and M is some common mass-scale of these new states. In turn, the strength of AM and CR interactions may be expected as a χ , b χ ∼ Q|y A,V | 2 /M 2 . It should be noted, however, that these estimates may be significantly enhanced by the lightness and/or mass-degeneracy of the spectrum of states [41]; a systematic study on EDMs induced by CP violation from light dark sectors was recently performed in [42]. In what follows, we treat the interactions (5) and (6) independent of their embedding.
For the Feynman-diagrammatic computation, one assembles the interactions into the matrix element of the effective EM current of χ, where p i,f and q = p i − p f are four-momenta. For a neutral particle χ the vertex functions reads, Here we regard the various moments as being generated at an energy scale well above the energies involved in the stellar production; they are hence q-independent.

III. STELLAR OBSERVABLES
In this section we review the arguments on stellar energy loss. Active stars such as RG, HB, or the Sun are systems of negative heat capacity: if energy is lost, either through photon emission or through new, anomalous processes, the decrease of total energy causes the gravitational energy to become more negative. By virtue of the virial theorem, the average kinetic energy and thereby the photon temperature increases. The system heats up leading to a faster consumption of its nuclear fuel while the overall stellar structure remains largely unchanged. In contrast, dead stars such as white dwarfs or the protoneutron star formed in core-collapse SN are supported by degeneracy pressure and stellar energy loss implies a cooling of the system. Constraints are then derived based on an observationally inferred cooling curve.
FIG. 1. Shown are the pair production processes of χχ that are calculated in this paper, namely, (a) plasmon decay, (b) e + e − annihilation, (c) 2 → 3 Compton scattering, (d) electron bremsstrahlung and (e) nucleon bremsstrahlung; for (c)-(d) we only show one of two relevant diagrams. The four momentum of the χχ-producing photon is denoted by k throughout the paper.

A. RG and HB stars
In globular clusters (GCs), the population of stars on the red giant branch vs. horizontal branch is directly related to the lifetime of stars in the respective phases. Their observationally inferred number ratio agrees with standard predictions to within 10%. Anomalous energy losses shorten the helium-burning lifetime in HB stars, creating an imbalance in the number of HB vs. RG stars. This constrains the luminosity in non-standard channels to be less than approximately 10% of the standard helium-burning luminosity of the HB core [3], core dVQ < 10% × L HB (HB) .
Following [3], L HB will be taken as 20 L for a 0.5M core below. The values of the Solar mass and luminosity are M = 1.99 × 10 33 g and L = 3.83 × 10 33 erg/s, respectively. The computation of the anomalous energy loss rate per unit volume and time,Q, will be the subject of the next section.
A constraint for RG stars may be derived from an agreement between predicted and observationally inferred core masses prior to helium ignition. Energy loss delays the latter and the core mass keeps increasing as the hydrogen burning "ashes" fall onto the degenerate He core. Preventing an increase in core mass by no more than 5% yields the constraint [3], Here,Q is to be evaluated at an average density of ρ = 2×10 5 g/cm 3 and a temperature of T = 10 8 K 8.6 keV, slightly higher than that of HB stars. The criterion (8) on energy loss can be improved utilizing high precision photometric observations of GCs. For example, considering the brightness of the tip of the RG branch, [43] has provided a detailed error budget and new limits on neutrino dipole moments from GC M5 were derived based on predictions of absolute brightness in the presence of anomalous energy loss that are obtained with dedicated stellar evolutionary codes. It was found, however, that previously derived limits based on (8) remain largely intact, as there appears to be a slight preference for anomalous energy loss channels [43]. In the following, for our purposes it will hence be entirely sufficient to employ the simple condition (8) to arrive at constraints on the EM form factors.

B. Sun
Solar neutrino fluxes are a direct measure of the nuclear fusion rates inside the Sun. For example, not only the 8 B neutrino flux is very well measured but also the sensitive dependence of the responsible reaction on temperature provides an excellent handle for constraining anomalous energy losses. The ensuing constraint is then phrased in terms of the total Solar photon luminosity [44,45], as It is important to note that (9) is basically insensitive to the long-standing "solar opacity problem": the measured 8 B neutrino flux is situated in the overlap region of the nominal error ranges between the discrepant highand low-metallicity determinations of the Solar chemical composition [8]; see the respective references [46] and [47]. Hence, (9) suffices as a criterion, awaiting further developments on Solar opacity determinations.

C. Supernova
New particles that are emitted from the proto-neutron star and that stream freely may quench the electroweak rates of neutrino emission during the cooling phase. The involved processes and their dynamics are highly complex. However, an approximate but very useful criterion to constrain additional energy loss is the condition that the total luminosity due to non-standard processes should not exceed the neutrino luminosity at one second after core bounce [3], core dVQ < L ν = 3 × 10 52 erg/s (SN). (10) The applicability of the bounds above are contingent on the produced particles being able to escape the dense environment of the (proto-)stellar core. Below, we will account for this so-called "trapping-limit" in the case of SN. For all other systems introduced above, trapping is either irrelevant, or happens in a parameter region that is excluded otherwise.

IV. PRODUCTION CROSS SECTIONS AND ENERGY LOSS RATES
In this section we first provide the general formula for χχ pair production in the thermal bath, before breaking it down into the most relevant pieces that dominate the in-medium production cross sections and, thereby, the stellar cooling rates.

A. Exact formula for χχ pair production
In thermal field theory, the production rate of a decoupled fermion per volume per time may be obtained from its relation to the imaginary part of its self-energy in medium [48] viȧ where Im Π χ (E χ , p χ ) =ū(p χ )Σ(E χ , p χ )u(p χ ) is the discontinuity of the thermal self-energy of χ, Σ(E χ , p χ ); u(p χ ) andū(p χ ) are free particle spinors with fourmomentum p χ = (E χ , p χ ). To lowest order in the dark coupling, Σ(E χ , p χ ) is found from the one-loop diagram with a dressed photon propagator attached to the χ fermion line. A general exposition on calculating discontinuities in the thermal plasma is found in [48,49]. Below, in Eq. (13), we are using a different formulation and the equivalence may be appreciated in the following way: when cutting the self-energy diagram for χ, the optical theorem implies that the production rate may also be obtained by computing all graphs where a photon γ * of four-momentum k = p χ + pχ emerges from a SM current and is being dotted into the dark current of the χχ pair. The SM-process that leads to the creation of γ * is in turn related to the imaginary part of the photon self-energy in the medium, Im Π µν , where Here T,L are the transverse and longitudinal photon polarization vectors and Π L,T (k) is thermal photon selfenergy for the respective polarization; explicit expres-sions are given in App. A. Identifying the leading contributions to Im Π L,T (k) in various medium then allows to account for the dominant χ pair production channels. The exact differential production rate per volume of χχ pairs via a photon of 4-momentum k = (ω, k) emerging from any SM process to lowest order in the dark current can be obtained by borrowing the results from dilepton production in hot matter, see e.g. [50,51]. Adopted to our purposes (see App. B) it reads, where s χχ = k 2 is the invariant mass of the χ-pair and the internal degrees of freedom of two polarization modes are g T = 2, g L = 1. The differences in the various interaction possibilities are entirely captured in a factor that will repeatedly appear and that was obtained in our preceding work [30], Equation (13) is the general expression of the weakly coupled χ pair-production rate from the thermal medium; details are found in App. B. The contribution to χχ production to leading order in α is given by the pole in (13), i.e. for s χχ = Re Π L,T . When this condition is met, (13) reduces to the decay rate of thermal photons γ L,T → χχ. Hence, resonant χχ production is fully accounted for by γ L,T decay. The decay itself becomes possible by virtue of the in-medium (squared) mass of γ L,T : it is given by Re Π L,T (ω L,T , k), where ω L,T denotes the solution of ω(| k|) of the corresponding longitudinal and transverse dispersion relations ω 2 − | k| 2 − Re Π L,T (ω, k) = 0. Plasmon decay is discussed in the following subsection, and explicitly calculated in (B7-B9) in App. B. The expressions for Re Π L,T and finite-temperature dispersion relations are found in (A5-A6). Production off-the-pole to Eq. (13) can be elucidated by studying the contributions to Im Π using the optical theorem, illustrated in Fig. 2. The left hand side shows the fully dressed vacuum polarization of an off-shell photon γ * , found by considering loop-diagrams of increasing order in α illustrated in the first equality. When those 2 Im = 2 Im Optical theorem relating the imaginary part of the photon self energy to the sum of all SM processes that create an offshell photon γ * . The first equality shows the leading individual contributions to the self-energy. When the latter loop-diagrams are cut, they correspond to the scattering processes shown in the second line, where dΠi symbolizes the phase space integral of all external particles, except γ * . When the scattering diagrams are deformed in a way such that two SM particles are in the initial state, the processes correspond to annihilation, Compton scattering and bremsstrahlung (from left to right). Any diagrams with χ particles involved yield contributions to the production rate (13) that are of higher order in the dark coupling.
loop diagrams are cut, their imaginary parts are given by the tree-level production processes for γ * shown in the last equality. The leading α contribution to Im Π is then given by the electron one-loop diagram. Although it is well known that on-shell plasmon decay γ L,T → e + e − remains forbidden at finite temperature [52], an electron loop still contributes to Im Π in the off-shell case. The associated process is then e + e − annihilation to χχ, i.e. process (2b). The second and third diagrams in the last line of Fig. 2 are related to χχ production in Compton scattering and bremsstrahlung. Here, it is important to note that s χχ = Re Π L,T can also be met in the photon propagator that produces the χ-pair with invariant squared mass s χχ . However, including such resonances would amount to double-counting. As we have seen above, the pole contributions are already captured by plasmon decay. 1 In our calculations, we explicitly avoid this situation by setting Π L,T → 0 in the propagator if the resonance is kinematically allowed for the photon that directly couples to the dark current. We have numerically verified that our results remain otherwise unaffected by neglecting the thermal shift in the photon propagator.
Finally, there is also a potential double counting between Compton scattering and bremsstrahlung processes, which happens when in the bremsstrahlung process the photon exchanged between two initial particles carries 4-momentum q (see Fig. 1d) that satisfies the dispersion relation q 2 − Re Π L (q 0 , q) = 0, leading to the exchange of an on-shell longitudinal plasmon. The process then becomes equivalent to Compton scattering e/N +γ L → e/N +χ+χ. This has been reported for axion production processes, where the contribution of latter is mostly covered by that of bremsstrahlung [55]. To avoid such double-counting, we take the static approximation (q 0 = 0) for the thermal mass of the photon exchanged in bremsstrahlung processes, which is a valid limit as the nucleon mass is large. As q 2 < 0 and Π L (q 0 = 0, q) is always positive, the exchanged photon can not become on-shell in bremsstrahlung processes, thus double counting is avoided (see Sec. IV D for more details).

B. γT,L decay
The on-shell process of photon decay to χχ (Fig. 1a) becomes possible in the medium and has an important analogy in the literature, the plasmon decay to neutrinos. Since the dispersion relation for transverse and longitudinal thermal photons are distinct, it is again helpful to separate the two polarizations in the calculation. Explicitly, we obtain for the decay rate per degree of freedom where ω T,L = ω T,L (| k|) for each polarization mode, as defined above. Details on the definition of the wave function renormalization factors Z T,L and the calculation are again given in App. B. In the limit of m χ → 0, the decay widths for MDM agree with the well-known formulae for plasmon decay to a neutrino pair [3]. For the plasmon decay processes, the energy loss rate can be expressed as [3], The expression for k max is given in Eq. (A7). For a non-relativistic medium (HB, RG, Sun), the dispersion relation crosses the light-cone at | k| = k max , signaling the damping of longitudinal modes (i.e. Landau damping); for a relativistic plasma (SN) k max → ∞. The relative factor of 2 between the expressions reflects the counting of polarization degrees of freedom. Finally, the last factor is a kinematic restriction on the phase space, ω 2 T,L −| k| 2 ≥ 4m 2 χ . For transverse mode thermal photons, the integral becomes bounded from below since ω 2 T − | k| 2 increases as | k| increases according to the dispersion relation. For the longitudinal case, the integral is additionally bounded from above since the trend in ω 2 L −| k| 2 with respect to | k| is reversed.
The degenerate plasma of the proto-neutron star with temperature T m e contains a population of e + , allowing for dark state pair-production through e + e − annihilation (Fig. 1b). The calculation for the pair production cross section is detailed in App. B.
In terms of the invariant s = (p 1 + p 2 ) 2 and the sum/difference of incoming e ∓ energies E 1,2 in the frame of the thermal bath, i.e., E ± ≡ E 1 ± E 2 , the corresponding cross section mediated by the transverse polarization part of the propagator reads, (17) For the longitudinal part we obtain (18) Note that there is no interference term between the two. Furthermore, the sum of both cross sections, σ T + σ L , becomes Lorentz invariant in the limit of Π T,L → 0. 2 Before using (17) and (18) in the calculation of the energy loss rate, a comment on the analytic structure is in order. Although it appears that the process may be significantly enhanced when s = Re Π T,L , this condition is never met: for the same reason that the decay of thermal photons into an electron-positron pair (γ T,L → e + e − ) is forbidden [52], the finite-temperature corrections to m e prevent the process (2) from going on-shell. It is for this reason that we have explicitly evaluated the thermal electron mass for the employed radial profile of the protoneutron star; see App. A. In other words, we use a thermal electron mass in SN, and use the zero-temperature electron mass in HB, RG and Sun, where e + e − annihilation is of little relevance. The values of chemical potential µ e are self-consistently adjusted to match the numerical SN profiles from the literature (see below).
The energy loss rate of e + e − annihilation is found by weighing the emission process by the total radiated final state energy E 3 + E 4 = E 1 + E 2 and by the probability of finding the initial states with the respective energies E 1 and E 2 , Here, f e ± are the phase-space distributions of e ± , with internal degrees of freedom g e ± = 2, and |M ann | 2 is the squared matrix element for e + e − annihilation into the dark state pair. It is further assumed that the population of dark states is small, The energy loss rate can be written in terms of the cross sections σ T,L . Borrowing from the discussion on phase space in [57], we find explicitlẏ The distribution functions f e − and f e + read Here, µ e is the chemical potential of electrons and T is the temperature. The threshold mass m th is equal to max{m e , m χ }.

D. Compton scattering
For 2 → 3 Compton scattering (e − /N + γ T,L → e − /N + χ +χ) with an initial γ T,L (Fig. 1c), we calculate the differential cross section via Here, σ 2→2 (s χχ ) is the cross section of the two-body Compton scattering with the final-state photon having a mass √ s χχ . We are only required to consider the process on electrons, e − + γ T,L → e − + χ +χ, as Compton scattering on protons is strongly suppressed. Following the treatment in [53] and our discussion above, we neglect the thermal mass of the final state photon in σ 2→2 to avoid any potential double counting with γ T,L decay. For the initial state photon in the integration of energy loss rate, Eq. (23) below, the thermal mass is properly taken in account through the dispersion relation (A6). Furthermore, note that there is no double counting between the Compton process and bremsstrahlung either in our treatment. A double counting would appear if the t-channel photon exchange in bremsstrahlung, with 4-momentum q (see Fig. 1d), goes on resonance. This is in principle possible for the longitudinal mode, since Π L in the propagator could become negative once the dispersion relation of γ L crosses the light cone. Nevertheless, in the electron bremsstrahlung process discussed belowmost relevant for RG, HB and the Sun-the proton recoil and hence the energy exchange are extremely small. Therefore, the propagator can be taken in the static limit (energy exchange q 0 → 0). This limit amounts to Debye screening, characterized by Π L (q 0 → 0, | q|). Since the screening scale is always positive, a resonance is never met. Therefore, we include the contribution from e − + γ L → e − + χ +χ to capture the t-channel resonance contribution of electron bremsstrahlung, although it is less important than plasmon decay.
The energy loss rate from Compton scattering is calculated in a similar way as (19), but here E loss is given by the energy carried by the virtual photon in the medium frame, where f 1,2,3 are the distribution functions of the incoming electron, γ T,L and outgoing electron, respectively, with g e − = g T = 2 and g L = 1 the internal degrees of freedom for the incoming electron and γ T,L . Pauli blocking is accounted for by including the factor (1 − f 3 ). The energy loss E loss = Eχ + E χ can be expressed in terms of variables defined in the medium frame. Moreover, for RG, HB and the Sun, the relativistic corrections induced by transforming from the center-of-mass (CM) frame to the medium frame are very small, and are neglected for simplicity.
In this subsection we consider dark state pair production from bremsstrahlung of electrons on protons or other nuclei (Fig. 1d). Similar to the Compton scattering above, we also relate the 2 → 4 cross section to a 2 → 3 process of eN → eN + γ * T,L in which the emitted photon, γ * T,L , has an invariant mass √ s χχ In the following we shall only consider photon-emission from electrons, as the emission from the nuclear leg is suppressed by a factor of (Zm e /m N ) 2 1 where Z and m N are the charge and mass of the nucleon/nucleus. Furthermore, ordinary electron-electron bremsstrahlung is a quadrupole emission process and correspondingly smaller in practice. We therefore also neglect such production channel.
The eN process is sensitive to the details of in-medium corrections. To this end, recall that the t-channel photon exchange in Fig. 1d has a well-known Coulomb divergence in the limit of vanishing momentum-transfer. This issue is mitigated by two factors: first, the divergence is not met kinematically as long as m χ = 0 since a minimum momentum transfer is necessary to create the final state pair. Second, the medium itself regulates the process through the Debye screening of bare charges characterized by a momentum scale k D . The latter appears as the static limit of Π L (q 0 → 0, q) and for a classical plasma reads, For the numerical results, we have calculated σ 2→3 in (24) using the propagator (A4), neglecting, for simplicity, ion contributions. We separate the squared amplitude into transverse and longitudinal parts and include the static limits of Π T,L in the respective propagators. For the longitudinal part, the zero-temperature propagator q −2 is replaced by (q 2 − k 2 D ) −1 . In contrast, there is no magnetic screening in the static limit (Π T (q 0 → 0, | q|) = 0), hence there is no thermal screening for the propagator of the transverse mode. We find that in the non-relativistic limit the contribution of the longitudinal mode dominates.
To avoid any double counting between this process and γ T,L decay, we need to subtract the contribution when the virtual photon that directly couples to χ goes onshell. As stated above, this is achieved by setting Π T,L in the corresponding propagator to zero. Since this should over-estimate the production rate at s χχ ≤ Π T,L , we have also tested an opposite option of choosing Π T,L → −Π T,L (E χ + Eχ) to avoid the singularity, which underestimates the production rate. We find that both prescriptions lead to same results at percent level, which justifies our simplification of taking Π T,L ≡ 0 for the photon that directly couples to χ.
For dark state pair production in e − bremsstrahlung on protons and nuclei, the energy loss rate is expressed aṡ where f 1,2,3 are the distribution functions of the incoming electron, proton/nucleus and outgoing electron, with g 1,2 the internal degrees of freedom for the incoming particles. We have neglected the Pauli blocking factor for finalstate protons/nuclei as it plays little role. The Møller velocity v M = F/(E 1 E 2 ) is given in terms of the flux factor F found in (B19). The energy carried-away by the dark states is E loss = Eχ + E χ and its expression in the medium frame is introduced in App. B.
Making the approximation that protons and other nuclei are at rest, their phase-space integral gives dΠ 2 f 2 = n N /(2m N g 2 ), where n N is the number density of the protons/nuclei. Hence we arrive aṫ with | p 1 | = E 2 1 − m 2 e and where σ 2→4 is obtained from integrating (24) over appropriate boundaries (see App. B 5). Generically, bremsstrahlung is less effective when pair annihilation or plasmon decay are open as production channels, but it can be dominant at low temperatures where the latter processes are kinematically suppressed.
Before ending this subsection, it is worth commenting on the so-called soft photon approximation, which states that in the limit that the emitted photon energy is small compared to the available kinetic energy (i.e., ω E kin ), the process of eN → eN +γ * T,L factorizes into a product of elastic scattering times a factor describing the additional emission of γ * T,L . While this approximation works well for the emission of a massless photon, it breaks down if the off-shell photon's effective mass is large, √ s χχ ∼ E kin . Overall, we find that the soft photon approximation describes the 2 → 4 process well for small m χ in the non-relativistic limit. However, for 2m χ ∼ E kin or for relativistic initial states the approximation fails, and it is ultimately related to the UV-sensitivity of the cross section (see App. C for details). Even though calculations simplify considerably in the soft photon limit, it cannot be applied for the whole m χ -range in electron bremsstrahlung and we therefore calculate σ 2→4 exactly, relegating details of the calculation to App. B. However, we will use the soft photon approximation in its region of validity to estimate the energy loss from nucleon bremsstrahlung in the next subsection.

F. Nucleon Bremsstrahlung
Proton and nuclear bremsstrahlungs are strongly suppressed in low-temperature environments due to the negligible thermal velocities of the initial states. However, in the interior of a proto-neutron star, the typical nucleon velocity is v 1/3 and N N -bremsstrahlung contributes to the total energy loss.
The photon that pair-creates the dark states is emitted from the proton-leg in proton-proton (pp) and neutronproton (np) scattering. Radiating off the neutron-leg through the neutron magnetic dipole moment is suppressed. In addition, since pp-scattering is associated with quadrupole radiation, it is suppressed with respect np scattering by a factor v 2 [58]. Therefore, we only consider np-scattering in the following. The interaction of protons and neutrons is mostly mediated by pions whose mass is of the order of the average momentum transfer in elastic collisions in SN, allowing for a separation of the phase space into an elastic and an emission piece (see App. C for details). The energy loss rate in the nonrelativistic, non-degenerate limit 3 is then given by [58] Q np = n n n p √ π (m n T ) where n n and n p are the neutron and proton number densities, m n is the average nucleon mass and E kin the available kinetic center of mass energy. For the elastic np-scattering transport cross section σ T np we use the numerical data from Fig. 3 in Ref. [58] for E kin 1 MeV and Fig. 2 in Ref. [59] for smaller energies; the emission piece I ω is given by Eq. (C8).
The energy loss rates from nucleon bremsstrahlung are comparable with the rates from plasmon decay and are depicted for m χ = 1 MeV as a function of the SN radius in Fig. 4.

V. CONSTRAINTS ON THE EFFECTIVE COUPLING
After calculating the energy loss rates induced by those relevant processes in each stellar environment (e.g. see Figs. 3-4), we apply the luminosity criteria introduced in Sec. III to obtain the upper bounds on the EM formfactors of light dark states.

A. Limits from RG, HB, and Sun
In this subsection, we derive the constraints coming from HB and RG stars utilizing the above calculated anomalous energy loss rates. For HB, we consider a representative star of 0.8 M and utilize the stellar profiles for density, temperature and chemical partition between hydrogen and helium from [3,60], reproduced in Fig. 7 in App. A. The luminosity of its helium-burning core is L HB = 20L to which (7) is then applied. For RG we use the prescription detailed below (8): a 0.5 M helium core with a constant density of ρ = 2 × 10 5 g/cm 3 and a temperature T = 10 8 K.
For the Sun, we use the standard Solar model BP05(OP) [61] to calculate the total power radiated into χχ which in turn is constrained from (9). For bremsstrahlung we take the contribution of electron scattering on H, 4 He and other less abundant nuclei ( 3 He, C,  N, and O). For simplicity, we assume all targets are in a fully ionized state. We find numerically that the contribution from the second class of elements contribute 10% of the total energy loss rate from bremsstrahlung, as the coherent enhancement from atomic charge number Z somewhat compensates for their scarcity in number.
The energy loss rates as a function of fractional stellar radius for HB (Sun) for all operators considered in (5) and (6) are shown in the left (right) panel of Fig. 3 for m χ = 10 eV and µ χ (or d χ ) = 10 −6 µ B and a χ (or b χ ) = 0.1 GeV −2 . MDM and EDM as well as AM and CR lines essentially yield identical results. This is owed to the fact that production proceeds in the kinematically unsuppressed region T m χ for which the energy loss rates match; the γ 5 factor discriminating the interactions of same mass-dimension only plays a role when χ particles become non-relativistic, hence close to kinematic endpoints. As can be seen, for dimension-5 operators the decay process (dotted lines) dominates over bremsstrahlung (dash-dotted lines) and Compton scattering (dashed lines) processes in both HB and Sun. For dimension-6 operators, the contribution of Compton scattering is comparable to that of decay processes in HB while in the Sun all three processes are of comparable importance.
Applying the criteria for the maximum allowable energy loss of Sec. III, we obtain the excluded shaded regions in Figs. 5 and 6 as labeled. The strongest limits are provided by RG stars. They have a higher core temperature, T = 8.6 keV, compared to HB stars or the Sun, favoring an emission process that is UV-biased because of the considered higher-dimensional operators. In the low mass region, for 2m χ < ω p , the limits are governed by γ T,L decays, and become independent of χ mass quickly. Once the decay process is kinematically forbidden, the limits become determined by the bremsstrahlung and Compton scattering processes. As can be seen, the critical values of m χ where this happens for RG, HB, and the Sun are reflective of the differing core-plasma frequencies (A3) of the respective systems. Furthermore, the mass-dimension 5 constraints on MDM and EDM are practically identical; differences only appear in the kinematic endpoint region.

B. Limits from SN
Limits on χ-photon interactions from SN have previously been estimated in our earlier paper [30], largely following the approach of [54], and considering e + e − annihilation but with plasmon decay neglected. Here we revisit these constraints in light of a more detailed calculation. Dark state pairs with mass m χ 400 MeV can be efficiently produced inside SN, predominantly through e + e − annihilation as positrons are thermally supported. Nevertheless, we will consider all processes in Fig. 1 except for electron bremsstrahlung as it is significantly weaker than the others; see Fig. 4 for one example with m χ = 1 MeV.
When the particles stream freely after production and are hence able to escape from the SN core, the limit (10) applies. We set the size of proto-neutron star to be r core = 15 km and model the proto-neutron star from which χχ pairs are emitted using the SN explosion simulation results of a 18M progenitor in [62] (see Fig. 7 in App. A). We adopt the total energy density ρ(r), temperature T (r) and electron abundance Y e (r) profiles at 1s after the core bounce. The number density of baryons can be computed as n b (r) ρ(r)/m p and the number density of electrons can be written as n e (r) n b (r)Y e (r). Other quantities such as chemical potential of electrons µ e (r), plasma frequency ω p (r) and effective electron mass m eff e (r) are derived from n e (r) and T (r). m eff e (r) is recursively solved at each radius using Eq. (A20); see App. A for details. In the calculations that relate to SN emission, m e is understood to be m eff e . The result is shown by the lower boundary of the region labeled SN1987A in Figs 5 and 6. Compared to our previous result in [30] where only e + e − annihilation was taken into account, the constraint for MDM and EDM is improved. This is traced back to the fact that the energy loss rate of γ T,L , nucleon-bremsstrahlung and Compton scattering for MDM and EDM are comparable to e + e − annihilation into a χχ pair. For AM and CR, however, the results from [30] remain largely unchanged, as γ T,L decay, nucleon-bremsstrahlung and Compton scattering are less efficient.
Once the effective coupling becomes large enough, the produced χ particles will eventually come into thermal equilibrium with SM particles. Here we follow [63,64] to divide the radial region into an inner "energy sphere" of radius r ES , within which light dark particles are thermalized, and an outer diffusion zone, where the χ luminosity gets attenuated by a transmission coefficient S ES . Both quantities are obtained from the energy-exchange and momentum-exchange mean-free-paths, i.e., λ E (r) and λ M (r), of the χ particle in the medium [63]: and where r inf is set to be 35 km, beyond which χ particles free-stream.
In the SN environment of interest here, the energyexchange mean-free-path is mostly governed by e-χ elastic scattering. We estimate it at each radius r through Here, ... denotes a thermal average over all participating particles; n e and n χ are the number densities of electrons and χ particles which are all assumed to be in equilibrium. dσ eχ /d cos θ is the differential cross section of e-χ elastic scattering. The energy-exchange in the scattering, |∆E|, is a function of the scattering angle θ in the rest frame of the medium. Although not explicitly written in the expression, we take the Pauli blocking factor for the final state electron into account. Meanwhile, the momentum-exchange mean-free-path λ M (r) turns out to be dominated by χ-nucleon scattering, where n i is the number density of neutrons (n) and protons (p), provided by the SN model [62]. The differential scattering cross section with χ is denoted as dσ iχ /d cos θ.
In obtaining our numerical results, we have neglected the nucleon velocities in the medium frame for simplicity. As we assume that χ stay in chemical equilibrium at r ES , its blackbody luminosity L χ (r ES ) is described by where g χ = 4. These χ particles are emitted towards the exterior from the energy sphere, but they continue to scatter elastically with the medium inside the diffusion zone. The χ flux-attenuation during its propagation from r ES to r inf can be estimated via in the diffusion limit [63]. That is, such attenuation only relies on the momentum-transfer cross section. Requiring L inf ≤ L ν then leads to upper boundaries of SN exclusion limits shown in Figs. 5 and 6. We note in passing, that the location of the upper boundaries is still conservative, as additional "freeze-in" χ-production for r > r ES has been neglected. With respect to our previous work [30], where the transmission ratio S ES was in practice taken as a simple unit-step function, the upper limits are improved by up to a factor of two.

C. Related works
Stellar bounds on the EM properties of light dark states have been studied in the literature, mostly in the context of EM properties of eV-scale (SM) neutrinos; see [65,66] and references therein. In these studies, the mass of neutrino is essentially zero. Therefore, in the limit m χ → 0 our results can be compared with previously derived constraints on neutrino EM interactions.
For higher-dimensional operators, [73] estimated the anomalous energy loss rate in SN through electron pair annihilation into light right-handed neutrinos, limiting its charge radius to be below 3.7 × 10 −34 cm 2 , that is 9.5×10 −7 GeV −2 , about seven times weaker than the one presented above. This is partially due to the fact that [73] assumed an one order of magnitude larger luminosity as the maximum permissible energy loss.

D. Cosmological constraints
Light dark states may lead to extra radiation in the early Universe, and thus its coupling to the SM bath is constrained by both the predictions from big-bang nucleosynthesis (BBN) and the observed cosmic microwave background (CMB). On the one hand, for the mass region considered here the CMB bounds depend on how it annihilates/decays. On the other hand, primordial abundance measurements of D and 4 He suggest that extra relativistic degrees of freedom need to be less than that of one chiral fermion during the nucleosynthesis (see e.g. [74][75][76]). Thus here we require that the Dirac fermion χ is thermally diminished at T ∼ 100 keV, either due to a feeble EM form-factor coupling or by a Boltzmannsuppression induced by its mass.
The relevant bounds are also given in Figs. 5 and 6. They only constrain the parameter region with m χ 1 MeV. In the same figures, we also show the line which corresponds to the thermal freeze-out scenario which generates the observed dark matter abundance, although such scenario has been excluded by various constraints for this model; see our previous work [30]. The dominant annihilation channel is into two photons at m χ < m e and into a pair of electrons at m χ ≥ m e , which explains the sharp decrease of the relic density curve at m χ ∼ m e seen in Fig. 5.

VI. CONCLUSIONS
In this paper we explore the sensitivity of stellar systems to neutral dark states that share higher-dimensional interactions with the SM photon. To this end we choose a Dirac fermion χ that is coupled to mass dimension 5 MDM and EDM operators with respective dimensionful coefficients µ χ and d χ and mass dimension 6 AM and CR operators with respective coefficients a χ and b χ . We consider anomalous energy losses from the interior of RG and HB stars, of the Sun, and of the proto-neutron star of SN1987A. Together with previously derived direct, indirect, and cosmological limits by us in [30], this work adds astrophysical constraints to draw a first comprehensive overview of light dark states with masses (well) below the GeV-scale and EM moment interactions.
The thermal environments of stellar interiors significantly affect (or enable) production processes of χχ pairs. Before breaking it down to individual contributions, we establish the exact formula, Eq. (13), for the pair-production rate in leading order of the dark coupling. The expression factorizes into a piece that represents the probability to produce an off-shell photon γ * from a SM current, and a piece that describes the production of the χχ pair from that photon. The former is proportional to the imaginary parts of the longitudinal and transverse thermal photon self energies Im Π L,T . The latter are model-dependent but otherwise universal factors that represent the choice of interaction, Eq. (14). The optical theorem then allows us to identify all major production processes by studying the contributions to Im Π L,T . The approach also allows us to clarify the role of thermal resonances in these processes, i.e., the kinematic situation when the pair-producing photon goes on-shell, k 2 = Re Π L,T . We find that resonant production is entirely captured by the decay of transverse and longitudinal thermal photons or "plasmons", γ T,L → χχ.
We compute the rates of χ-pair production and its ensuing energy loss from plasmon decay and Compton production for all systems. In addition, we evaluate eN bremsstrahlung for RG, HB and the Sun, and e + e − annihilation and N N bremsstrahlung for SN1978A. For MDM and EDM interactions, plasmon decay dominates in HB and RG stars and in the Sun. For the interactions of increased mass dimension, AM and CR, the Compton (bremsstrahlung) production dominates in HB and RG (Sun). In SN, e + e − annihilation dominates the anomalous energy loss for r 7 km. In the most inner region the population of positrons becomes extremely Boltzmann suppressed by a decrease in temperature, and plasmon decay and np bremsstrahlung take over as the most important production channels. For all processes we have taken into account all important finite-temperature effects. Furthermore, in the evaluation of rates, we explicitly avoid any double counting between plasmon decay and an on-shell emitted photon in bremsstrahlung and between Compton production and an on-shell exchanged t-channel photon in bremsstrahlung.
The rates when integrated over stellar radius then become subject to the observationally inferred limits on anomalous energy loss. The resulting restrictions on the parameter space are found in Figs. 5 and 6. In the kinematically unrestricted regime m χ 1 keV, the stellar limits are dominated by RG with µ χ , d χ ≤ 2 × 10 −12 µ B and a χ , b χ ≤ 6 × 10 −5 GeV −2 . All interactions are additionally constrained from SN1987A, in the windows 10 −10 µ B ≤ µ χ , d χ ≤ 10 −8 µ B and 10 −7 GeV −2 ≤ a χ , b χ ≤ 10 −3 GeV −2 for m χ 10 MeV. The SN constraining region is bounded from above by the trapping of χ particles, which we evaluate in some detail. The presented astrophysical constraints add to a program that we have started in [30] and that aims at charting out the experimental and observational sensitivity to effective dark state-photon interactions. The stellar constraints on anomalous energy loss derived in this work yield the most important limits on the existence of effective dark sectorphoton interactions for χ-particles below the MeV-scale.
Acknowledgments The authors are supported by the New Frontiers program of the Austrian Academy of Sciences. JLK and LS are supported by the Austrian Science Fund FWF under the Doctoral Program W1252-N27 Particles and Interactions. We acknowledge the use of computer packages for algebraic calculations [77,78].

Appendix A: Photons in a thermal medium
The processes depicted in Fig. 1 are fundamentally affected (or enabled) by the in-medium modified photon dispersion. Here we collect the central results that go into the computation of the energy loss rates (our convention largely follows [3]). The central quantity measuring the strength of the medium-effect is the plasma frequency ω p , obtained through where v = p/E is the velocity of electrons or positrons, and f e − and f e + are their respective Fermi-Dirac distributions, f e ± = [e (E±µe)/T + 1] −1 . Equation (A1) takes on the following analytic forms in the classical, degenerate and relativistic limit respectively, where α is the fine-structure constant, n e is the number density of electrons, p F = (3π 2 n e ) 1/3 is the Fermi momentum, E F = m 2 e + p 2 F is the Fermi energy and v F = p F /E F is the Fermi velocity. Here "classical" refers to a non-relativistic (T m e ) and non-degenerate (T µ e − m e ) plasma. The proto-neutron star of a SN is both in a relativistic and degenerate regime and we find that the relativistic limit above yields a better fit to the general form of ω p in (A1) than the degenerate limit; the latter exhibits a 10% deviation. The core of a RG star is non-relativistic but degenerate whereas HB stars and the Sun are well described by the classical limit. In our numerical calculations, we adopt ω p computed from Eq. (A1), avoiding any ambiguities of taking limiting cases. Representative values of ω p at the cores of all stellar objects are summarized as Sun's core 2.6 keV HB's core 8.6 keV RG's core 17.6 MeV SN's core . (A3) The computation of most of the processes requires the inmedium photon propagator. Picking Coulomb gauge, for a photon carrying 4-momentum k = (ω, k), the latter divides into longitudinal (L) and transverse (T) parts [79], where k i is the Cartesian component of the photon threemomentum (magnitude | k|). Using [79] and adopting the conventions of [3,5], the real part of the polarization functions Π T,L in the rest frame of the (isotropic) thermal bath reads, The full expressions for the dispersion relations k 2 − Π L,T = 0 then relate the energies of an on-shell photon, ω T,L , to its momentum k to order α [79], via Equations (A6) are also valid to order | k| 2 at small | k| for all temperatures and electron number densities. Throughout the paper, we always use ω T,L , as functions of | k|, to denote the energy of an on-shell thermal photon, which satisfies Equations (A6), and use ω for off-shell photons.
Longitudinal photons are populated up to a wavenumber k max , beyond which the longitudinal dispersion relation crosses the light cone and L-modes become damped, with and in the relativistic limit k max → ∞. In these equations, the mobility of charges is captured by the typical velocity of electrons, v * ≡ ω 1 /ω p , where In the three limits mentioned previously, v * can be ap- Finally, as alluded to in the main text, the processes we consider are non-resonant in the photon exchange and Im Π T,L can be neglected throughout.
Furthermore, the in-medium coupling of the photon to the EM current is modified by the vertex renormalization constants Z T,L ≡ (1 − ∂Π T,L /∂ω 2 T,L ) −1 . For the convention adopted here, they are equivalent to the ones given in [3] (A11) These factors are attached to each zero-temperature vertex factor involving an external photon state. For internal photons, this effect is already accounted for in the momentum-dependent self-energy Π T,L (ω, k).
For thermal corrections to the electron mass which is relevant for SN, we closely follow [80]. For an electron with a general 4-momentum p = (E, p) in a neutral medium where the positron number density is negligible, we first introduce the four functions below where q here is the absolute value of the 3-momentum of medium particles (electron, photon) that is integrated over and L 1,2,3,4 are functions of q in terms of Here m e is the zero-temperature mass of electron, 0.511 MeV, while f γ (p) and f e − (p) give the thermal momentum distribution functions of photon and electron (per degree of freedom). We have set f e + (p) = 0 in the above equations. In the end, we take the approximation made in [81] to obtain that We have neglected thermal corrections to χ states. In the phenomenologically relevant regime, their coupling to the thermal bath is very weak. Finally, we have reproduced the profiles of the HB model from [60], and SN model from [62], adopted in this work, as shown in Fig. 7, where neutrality and µ e − + µ e + = 0 at each radius have been taken for granted for the SN profile.
where the coefficients A and B are given by with s χχ = k 2 . In terms of the functions f (s χχ ) defined in (14), the factor I νσ is then explicitly given by, Putting all of the above together, we obtain the differential production rate per volume (13) found in the main text which we repeat here for convenience, where f B (ω) = (e ω/T −1) −1 is the Bose-Einstein momentum distribution of thermal bosons. During the derivation we have used that

Leading contributions to Im ΠT,L
We now demonstrate that (13) or, equivalently, (B3) contain the leading production mechanisms considered in this paper. In particular we clarify the role of resonances, and that they are accounted for by the process γ T,L → χχ. To this end, we isolate the pole contribution to the total production rate, i.e. the case s χχ = Re Π L,T . To this end, we adopt the narrow width approximation, where Im Π L,T < 0. Then noting that Re Π L,T is also a function of s χχ and writing d 4 k = d 3 kds χχ /(2ω) yielḋ where g T = 2 and g L = 1, counting the degrees of freedom of the photon modes. Now both ω T,L and Re Π L,T need to satisfy the photon dispersion relation with a 3momentum k due to the δ-function above. As will be calculated below and given explicitly in (15), the term in the parentheses is exactly the decay rate of γ L,T into χ pairs.
In a next step we further verify that the contribution of the one electron loop (OEL) to Im Π µρ induces the production rate of χ from electron pair annihilation 4 , process (2b). In this case, it is easier to start with (B3), where according to the in-medium optical theorem (see Fig. 2) we may write where T µρ e = M µ γ * →e + e − M ρ e + e − →γ * , and f e ∓ gives the momentum distribution function of e − (p 1 ), e + (p 2 ) per degree of freedom as defined above. Moreover, terms that are kinetically forbidden for k 2 > 0 have been neglected [48]. The presence of (1 − f e − − f e + ) is due to quantum statistics, and would disappear for classical particles. Substituting this expression into (B3) gives Then for the Fermi-Dirac distribution function f e ± and the Bose-Einstein distribution f B (ω) with the energy conservation E 1 + E 2 = ω, there exists the relation allowing us to rewrite the number production rate per volume above aṡ after integrating over d 4 k on both sides. The last expression transforms precisely to the corresponding energy loss rate (19), once both the energy-loss factor (E 1 + E 2 ) and fermionic degrees of freedom f e ± are taken in account.

γT,L decay to dark states
In the following appendices we calculate the leading processes in the usual Feynman-diagrammatic approach using tree-level perturbation theory augmented by the thermal corrections outlined in Appendix. A. The decay of a transverse or longitudinal photon of 4-momentum k to a pair of dark statesχ(pχ) + χ(p χ ) is described by the spin-summed squared matrix element where Z T,L is the vertex renormalization factor in (A11), µ is the photon polarization vector and T µν χ is given in (B2). The decay rate is given by the phase-space integral, where ω T,L is the energy of the external transverse or longitudinal photon. It is useful to employ (B5). In terms of I νσ defined in (B6), we can write the decay rate as The explicit expression, given by (15) in the main text, is then found by using the expressions (A10) for the polarization vectors when the initial state propagates in positive z-direction, i.e. for k µ = (ω, 0, 0, k); note that the term proportional to k µ k ν in I µν does not contribute due to the Ward identity.

e + e − annihilation to dark states
Here we consider the process e − (p 1 ) + e + (p 2 ) → χ(p χ ) +χ(pχ). By setting p i = (E i , p i ) and k = p 1 + p 2 , one can define the cross sections 5 in terms of the squared matrix element for annihilation |M ann | 2 , where T µρ e reads T µρ e = −2(sg µρ − 2p µ 1 p ρ 2 − 2p ρ 1 p µ 2 ), and I νσ χ is given in (B6); here s = k 2 . Furthermore, v M is the Møller velocity defined as v M = F/(E 1 E 2 ) and the flux factor F is given by Contracting the Lorentz indices then yields where the interference term vanishes in the Coulomb gauge, as can also be seen from (13), and the cross section for each polarization mode is given as equations (17) and (18) in the main text. 5 We emphasize that this leads to a Lorentz-invariant total cross section up to the thermal mass of photons, which is convenient for the phase space integral.
5. e − N bremsstrahlung production of dark states The 2 → 3 amplitude squared |M 2→3 | 2 can be split into three parts as spins |M 2→3 | 2 (4πα) 2 g 1 g 2 = D ρβ (q)D σγ * (q)W ρσ L µν βγ * µ (k) ν (k), (B21) where q = p 2 − p 4 is the momentum transfer between the initial states, L µν βγ stands for the leptonic part, W ρσ is the hadronic tensor and ν (k) is the polarization vector of the emitted photon of virtual mass s χχ = k 2 . Detailed forms for L µν βγ and W ρσ are given in the App. A of our previous work [30].
The integration boundary of s 4 is given by, The physical region for p 1k is expressed by n × n asymmetric and symmetric Gram determinants, G n and ∆ n .
It reads, and the Jacobian |∂φ R4k Putting everything together, the full 2 → 4 cross section is given by The integration boundaries of s χχ are given by where µ is the reduced mass, v is the relative velocity and E (E ) is the initial (final) total CM energy of the colliding particles. Even though in deriving the expression in Eq. (C2) we have assumed √ s χχ E kin and ω E kin , integrating (C2) over ω in the region √ s χχ < ω < E kin (C3) still gives a very good approximation to the full cross section. Obviously, the approximation breaks down if √ s χχ ∼ E kin where the integration region of ω gets very small and the integral is dominated by large emission energies.
To obtain the 2 → 4 cross section from the 2 → 3 cross section in Eq. (24), σ soft 2→3 gets multiplied by the factors in Eq. (14) corresponding to the EM form factor interactions. This leads to the following parametric dependence on s χχ , that is, for dim-4 operators, like millicharged states, the 2 → 4 cross section is dominated by small s χχ , whereas for higher dimensional operators, the expression is UV biased. Hence, the main contribution to the integral comes from s χχ -values for which the soft approximation breaks down as one probes the kinematic endpoint region. It turns out that in the non-relativistic regime and for m χ + mχ E kin , using Eq. (C2) reproduces the exact 2 → 4 cross section up to a factor 2 or 3. However, for relativistic particles, the error at √ s χχ ∼ E kin gets larger. Due to the s χχ -dependence in Eq. (C4), this still results in a decent description of millicharged χχ emission, but produces errors of several orders of magnitude in the relativistic regime for the EM form factors considered in this paper. Equation (C1) can be further simplified by separating the phase space. This is possible, if the elastic scattering cross section is insensitive to an angular cut-off in the forward or backward direction, e.g. if the interaction is mediated by a massive particle such as the pion in np scattering. 7 Then, where σ T 2→2 is the transport cross section and the emission piece I(s χχ ) is obtained by executing the integral in Eq. (C1) In the first line, we have divided by (1 − cos θ) to cancel the θ-dependent part in the emission piece, which we have absorbed into the elastic cross section. In Sec. IV F we make use of this factorization in calculating the energy loss rate for neutron-proton scattering in SN. For that, the integral in Eq. (C7) is weighted with ω to obtain with x = √ s χχ /E kin , which is in agreement with the findings of Ref. [58].