Supernova constraints on an axion-photon-dark photon interaction

We present the supernova constraints on an axion-photon-dark photon coupling, which can be the leading coupling to dark sector models and can also lead to dramatic changes to axion cosmology. We show that the supernova bound on this coupling has two unusual features. One occurs because the scattering that leads to the trapping regime converts axions and dark photons into each other. Thus, if one of the two new particles is sufficiently massive, both production and scattering become suppressed and the bounds from bulk emission and trapped (area) emission both weaken exponentially and do not intersect. The other unusual feature occurs because for light dark photons, longitudinal modes couple more weakly than transverse modes do. Since the longitudinal mode is more weakly coupled, it can still cause excessive cooling even if the transverse mode is trapped. Thus, the supernova constraints for massive dark photons look like two independent supernova bounds super-imposed on top of each other.


Introduction
One of the most extreme environments in the universe are the centers of core-collapse supernovae. For many seconds right after the collapse, supernovae reach nuclear densities and temperatures of order 10s of MeV before they are cooled via neutrino emission. This makes them a potential source for particles beyond the Standard Model (BSM) if their masses are below about 100 MeV. The detection of neutrinos from the supernova SN1987A [1][2][3] confirmed the theoretical picture that the proto-neutron stars formed during core collapse supernova cool via neutrino emission. Thus, one can place bounds on BSM particles by requiring that they cannot be emitted so much as to cool the supernovae faster than neutrinos [4][5][6]. This allows one to constrain particles that are so weakly coupled that they can escape the core without interacting, which results in a bulk emission constraint. For larger couplings, the new particles become trapped and the cooling is effectively only done by particles emitted at the last scattering surface, which acts like a blackbody of the trapped particles. The trapping regime leads to an upper bound on the couplings that can be constrained by cooling while bulk emission leads to a lower bound. Generically, the lower bound curve from bulk emission and the upper bound curve from surface emission in the mass vs coupling plane intersect at large masses, leading to a closed region that can be excluded by cooling considerations.
In this paper, we study the supernovae cooling constraint on the axion-photon-darkphoton coupling, and show that this model leads to two new qualitative features in supernovae that can play a role in non-minimal scenarios where there are more than one new degree of freedom. The first feature we discuss is that different polarizations of particles can behave differently, and thus their production, and in particular their trapping, needs to be treated separately rather than averaged over 1 . An example of where something like this occurs is when discussing transverse and longitudinal modes of a dark photon. The longitudinal mode couples more weakly than the transverse mode and needs its own cooling and trapping constraints.
The second feature arises when the new particles must be pair produced and the two species have different masses. For fixed coupling, as the mass of one particle becomes larger, their joint production is suppressed because the threshold energy for pair production is increased. However, because the same interaction vertex appears in scattering processes, the scattering of the lighter particle becomes suppressed as well, as it does not have enough energy to up-scatter into the heavier particle. As a result, both production and scattering become suppressed so that the associated bounds both become weaker exponentially, and the trapping limit and bulk emission limit curves do not intersect.
In the axion-photon-dark-photon model we study in this paper, there are two new particles, the axion and the dark photon, and their leading coupling is given by where a is the axion and F (F D ) is the field strength tensor of the photon (dark photon). The axion and the dark photon are two of the most motivated candidates for physics beyond the standard model. Axions are pseudo-Goldstone bosons that appear very naturally in string theory [23] and can solve a variety of problems such as the strong CP problem [24][25][26][27], or dark matter [28][29][30]. Dark photons appear in a wide range of motivated models of dark matter [31][32][33] or as a solution to various anomalies [34][35][36][37][38]. While in many models the leading coupling of the dark photon to our sector is through kinetic mixing, this coupling can be highly suppressed making the interaction in Eq. 1.1 the most relevant one for phenomenology. In particular, if the dark sector has a charge conjugation symmetry, C D , then kinetic mixing is removed. In this case, as long as the axion is odd under both C D and CP , Eq. 1.1 is the leading coupling between the dark sector and the visible sector. A coupling of the form shown in Eq. 1.1 has been considered before in many different contexts, e.g. Ref. [39][40][41][42][43][44][45][46][47][48]. Additional motivation for studying this coupling comes from dark photon dark matter. Dark photons are a compelling dark matter candidate and so production mechanisms for it are interesting. Above a keV, dark photon dark matter can be produced thermally, either through freeze-out or freeze-in. As it becomes lighter than keV, one needs to consider non-thermal production mechanisms [49][50][51][52][53][54]. Dark photon production from axions using the coupling in Eq. 1.1 is one of the few ways to realize dark photon dark matter lighter than a µeV [55]. This coupling also leads to the two new aforementioned features on the cooling bound, it has a dark photon whose longitudinal mode and transverse modes couple to matter with different strengths. Additionally, because the coupling Eq. 1.1 contains both the dark photon and the axion, production and scattering necessarily involve two separate species, which must be taken into account in cooling constraints. The bounds on this coupling are shown in Fig. 8 where both of these features can be seen. At low mass, there are two separate cooling and trapping constraints while at high mass, both the trapping and cooling bounds become exponentially weaker.
In Sec. 2, we review how core-collapse supernovae can be used to place constraints on new physics and derive the cooling constraints on our model. In Sec. 3, we derive the trapping constraints. We discuss in detail the longitudinal mode of the dark photon and its constraints in Sec. 4. Non-supernova constraints are discussed in Sec. 5. Finally, in Sec. 6 we conclude.

The Energy Loss Argument and Emission Constraints
The core collapse supernova SN1987A has provided, and continues to provide, valuable constraints on new physics. Measurements of the energy emitted into neutrinos and the duration of the emission were used to infer the cooling rate of the supernova. This cooling rate was found to be consistent with Standard Model predictions. Any new particle weakly coupled to the Standard Model provides a new channel for energy loss in the supernova. This would increase the cooling rate, pushing it out of agreement with Standard Model predictions, and provides a constraint on the coupling of new light particles to the Standard Model. This constraint can be expressed in a simple form with the Raffelt criterion [6], which states that the luminosity of new particles, L x , can be no larger than the luminosity of neutrinos, L ν , L x ≤ L ν = 3 · 10 52 erg·s −1 . (2.1) Depending on how strongly the new particles interact, cooling can occur in two different regimes. If the new particles interact very weakly, after being produced they carry away energy without further interactions inside the proto-neutron star. This corresponds to the bulk emission regime and is used to place a lower bound on the coupling, dependent only on the amount of energy going into the new particles. Once the interactions of the new particles are sufficiently strong, they become trapped in the proto-neutron star and the cooling is only done through emission at their last scattering surface (trapping surface), i.e. as in a blackbody. As the coupling becomes stronger, this trapping surface moves to larger radii and thus to lower temperatures, which results in a smaller luminosity. This leads to an upper bound on the couplings that can be excluded by Raffelt's criterion. In this section we will examine the constraint arising from bulk emission and discuss the trapping constraint in Sec. 3. In our model the axion and dark photon play the roles of these weakly interacting particles. To find the luminosity of axions and dark photons we must find the emissivity,˙ , or energy emission rate per unit volume. Integrating this over the volume of the supernova gives the luminosity of the bulk emission, Next, we must compute the emissivity for each process that can produce axions and dark photons in the supernova. To find these emissivites, we simply take the interaction rate per unit volume for that process and weigh it by the energy carried away by the axion and dark photon in the final state [15], (2.2) Here, momenta indexed with i represent initial state particles while momenta indexed with f represent final state particles. f i represent the thermal weights of the initial state particles. The matrix element is averaged over initial spins and summed over final spins. E a and E D are the energies of the produced axion and dark photon. For computational ease, we neglect Pauli blocking and Bose enhancement of the final states. The latter is never a large effect in the processes under consideration, while the former can significantly reduce the rates for processes involving scattering with electrons at the core. As we will show later, those processes are subdominant even without taking Pauli blocking into account, and thus we will not include them when deriving constraints. When Eq. 2.2 is simplified for processes involving two initial state particles (A and B), the emissivity can be nicely written as the product of initial densities multiplied by a thermally averaged, energy weighted cross sectioṅ where σ E is the ordinary cross section weighted by the sum of energies of the axion and dark A similar simplification can be done for the emissivity for any process with a single initial state particle (A). The result is a product of the initial densities and the average decay rate weighted by energy is˙ where Γ is the usual decay rate. There are 4 relevant processes for the production of axions and dark photons in supernova. 3 of them are collision processes involving 2 initial state particles: the annihilation of an electron and positron, Compton scattering of an electron and photon, and nuclear bremsstrahlung of a proton and neutron. An example diagram for each of these processes is shown in Fig. 1. The emissivites for these processes are computed in App. A and the results are given in where E tot is the average of the sum of the initial state particle's energies. The Q factors in each expression, defined in Eqs. A.8, A.15, A.20 in App. A, contain all of the dependence on the masses of the axion and dark photon. They are normalized so that they are O(1) when m a = m D = 0 and then decrease as the axion or dark photon gain mass.
The annihilation and Compton emissivities have a simple structure. They include the couplings involved in the process, the initial densities and the average of the total initial energy. The bremsstrahlung emissivity has a different structure due to how it was computed, which we do following Ref. [16]. In order to deal with the nuclear scattering, which does not allow for a simple perturbative calculation in the energy regime of interest, we assume that the process is dominated by scatterings where the radiated particles (and hence the virtual photon in the diagram) have energies much smaller than the kinetic energy of the nucleons 2 . This approximation, called the soft radiation approximation, allows us to factorize the emissivity into two factors: one factor accounting for the nuclear scattering and another capturing the decay of a virtual photon into the axion and dark photon. The factor Σ pn is the piece that accounts for the dynamics of the nuclear scattering and is proportional to the momentum transfer cross section of proton-neutron scattering. In detail, it is given by In order to compute the momentum transfer cross section we use experimental data from Ref. [14] and for more details see App. A. Now we consider the decay processes that could produce axions and dark photons. To leading order, there is only one such process: the decay of plasmons into dark photons and axions. The diagram of this process is shown in Fig. 2. In order to compute the emissivity, we've used the polarizations and dispersion relations for the plasmon given in Ref. [56]. Again, the details of the computation are given in App. A, but the result is˙ where Q p defined in Eq. A.24 contains all of the dependence on the mass of the axion and dark photon. It is O(1) when m a = m D = 0 and decreases for increasing dark photon and axion masses. ω p is the plasmon frequency defined in Eq. A.22. Fig. 3 shows a plot of the emissivities for m a = m D = 0 as a function of supernova radius. These emissivities are computed using the temperature and density profiles given in App. B.  Near the center of the core, Compton scattering and plasmon decay dominate the energy emission while at larger radii, annihilation dominates. Note that the Compton emissivity does not take Pauli blocking into account, and thus the true rate would have a large suppression compared to what is shown in Fig. 3. We will not include the Compton contribution to the luminosity, which is a conservative approach. Since the Compton contribution to the total luminosity is subdominant, neglecting it only leads to a small effect. Upon integration of each emissivity over the volume of the supernova, we find that the annihilation is the dominant process, with plasmon decay being a significant contribution for lower masses. Finally, these emissivities are summed into a total emissivity and integrated over the volume of the supernova to give the total luminosity of axions and dark photons. This total luminosity is then applied to the constraint in Eq. 2.1 for various axion and dark photon masses and is shown in the lower line in Fig. 5.

Trapping Constraints
As discussed in Sec. 2, at stronger couplings cooling is done through surface emission due to particles becoming trapped. Using this surface emission, one can derive an upper bound on the couplings from the Raffelt criterion in Eq. 2.1. Once the coupling is large enough, any new particle, x, will be trapped inside the supernova below some radius r x . This radius, r x , is defined as the point where the optical depth τ is equal to 2/3 The optical depth is defined as where λ x is the mean free path of the new particle and is related to the scattering process for the new particle x with Standard Model particles. Because of rapid absorption and reemission, the emitted energy due to this new particle is no longer the bulk emission seen in Sec. 2 but can be approximated as blackbody radiation from this radius 3 r x given by Here m x is the mass of the particle and g ,x is its number of degrees of freedom. This resembles the familiar expression for blackbody radiation with g ,xπ 2 120 being the Steffan Boltzmann constant in natural units and an additional factor h(m x /T ). This factor is to account for the mass of the radiated particle and is defined by This integral is such that when the particle is massless, h(0) = 1 and we recover the familiar blackbody formula. In the limit m x T , h ∼ e −mx/T meaning that for very massive particles, blackbody radiation is highly suppressed.
The luminosity of this blackbody radiation has a different dependence on the coupling than the bulk emission considered in Sec. 2. Since for radii larger than the core radius, ∼ 10 km, T falls off faster than r −1/2 (see App. B), the total luminosity decreases with increasing r. This means that as the new particle becomes more and more weakly coupled, it can travel through the star easier and so the radius of blackbody emission shrinks. As the blackbody radius shrinks, the emitted luminosity increases (assuming it is outside the proto-neutron star's core) and so in the trapping limit, the Raffelt criterion excludes all couplings weaker than a certain value. Now we will apply this trapping constraint to our model. Both the axion and dark photon have their own mean free paths λ a and λ D from their decay and scattering with SM particles. From these mean free paths we define the two radii r a and r D via Eq. 3.1. Effectively, these radii can be thought of as the radius below which the corresponding particle is trapped and outside of which it free-streams.
The difference between r a and r D introduces extra complications that precludes a simple treatment of the luminosity of both dark photons and axions as arising from a blackbody. This arises from the fact that in the region between the two radii neither of the particles are in thermal equilibrium with the rest of the proto-neutron star, and so the assumptions that go into treating the outer radius as a blackbody fail. To illustrate these complications, let us consider the case where r D < r a , i.e. the dark photon becomes free-streaming at a smaller radii than the axion. While the dark photon can be treated as a blackbody emitting from r D , the axion emission is not so simple. In the region between r D and r a , the axion scatters or decays before traveling a large distance, but the dark photon can easily escape the region without subsequent interactions. Because any interaction/decay of the axion converts it into a dark photon, even though from Eq. 3.1 we find that the axion is trapped, after a single interaction a portion of its energy goes into dark photons that can escape the region. Therefore, we can treat the smaller radii as a blackbody for both particles, but we must take into account the transmission probability for the energy that is radiated into axions.
Once produced at the inner trapping radius, the particle associated with the larger trapping radius travels a finite distance before having an interaction (scattering or decay) that converts it into the free-streaming particle. The luminosity associated with the trapped particle depends on the probability that the free-streaming particle created in that interaction travels outward (so that it does not return to the inner radius), P θ , and the fraction of the initial particle energy carried by the free-streaming particle, . Calculating those factors precisely is computationally intensive and amounts to small changes in the final result. To simplify, we take P θ and to have the same value for all scatterings. So that the total luminosity is L tot (r in ) = L Blackbody (r in , m in , g ,in ) + P θ L Blackbody (r in , m out , g ,out ), (3.4) where m in , m out , g ,in , and g ,out are the masses and degrees of freedom of the particle corresponding to the inner and outer radius respectively. In reality, P θ depends on how far the initial particle makes it from the blackbody before scattering since the further it makes it out, the more solid angle it can scatter into without being directed back into the blackbody. Thus we can say that 1/2 ≤ P θ ≤ 1. To get an estimate for we can easily see that depending on the mass of the initial particle and whether the event is a scattering or a decay that in the center of mass frame the final particle takes between 1/2 and all of the initial particle's energy. Thus we have 1/2 ≤ ≤ 1. In order to get a conservative estimate on our coupling we note that if we underestimate our total luminosity, we will also underestimate the radius of the blackbody once we solve L tot (r in ) = L ν . This means we will be underestimating the constraint on the coupling thus giving us a conservative upper bound. We therefore take the conservative estimate of = P θ = 1/2. In order to get a sense of the maximum value our constraint could be off by, we look at the other extreme = P θ = 1. This changes the value of the constraint curve by at most around 20%, so we use the value P θ = 1/4 as it provides a good conservative estimate of the constraint. Now we compute the mean free paths of the axion and dark photons. The mean free path for both particles is determined by their scattering cross section with the three electrically charged fermions in the supernova: electrons, positrons, and protons as well as their decay width (where only the heavier between the axion and dark photon decays). The diagrams for these scatterings and decays are shown in Fig. 4 with ψ designating any of the three fermions. In order to compute the mean free path, we must compute the cross sections and decay widths from these interactions. However, rather than using the ordinary cross sections to determine the mean free paths, we will use the momentum transfer cross section defined as Using the momentum transfer cross section means that we are favoring scatterings with a more severe angular deflection because they are more effective at trapping particles than those with a small angular deflection. The momentum transfer cross section also helps to regulate the t channel divergence that appears in our diagrams when m a = m D . The computation of these cross sections is detailed in App. A, where we show the cross-sections take the form where m ψ is the mass of the fermion the axion or dark photon is scattering with. The definition for I is given Eq. A.27 in App. A, and is exactly 1 for m a = m D = 0. The decay widths are computed in App. A and are where E is the energy of the initial particle. The mean free path λ for the axion and dark photon are defined as where in the sum ψ represents: electrons, positrons and protons. σ a,ψ is a short hand notation for the momentum transfer cross section of an axion with fermion ψ and likewise σ D,ψ for the dark photon. Additionally, the subscript ψ in the thermal average in Eq. 3.8 is meant to convey that we are taking the average of initial fermion states only. For the initial boson (axion or dark photon) energy we simply input the average energy of a boson at radius r i . It is important to highlight why inserting the average energy for the initial boson (call it X), rather than averaging over initial energies is the more physically accurate thing to do. This stems from the fact that particles of different energies decouple at different radii, and thus treating this radiation as a blackbody emission is an approximation. This approximation can become very inaccurate when applied to particles that have inelastic scatterings. Consider the scenario where the final state boson has a very large mass, M . In order to have enough energy to produce the much more massive final state boson, the particle ψ that X scatters off of needs to have very high energy. The lower the energy of X, the higher the energy of ψ needs to be. This implies that X with a lower energy will scatter less frequently than X of higher energy. In the limit that M T , most of the initial state bosons are effectively free while only those with very high energy have a chance at scattering. Averaging over initial state X energies fails to capture this behavior effectively because it will assign every X the same mean free path regardless of energy. This would lead us to conclude that the typical thermal bosons in the supernova scatter much more effectively than they actually do which in turn would lead us to overestimate the extent to which these particles are trapped in the supernova. In order to capture the dynamics, which is dominated by particles with energies close to the average energy, it is more physically accurate to use the mean free path for particles with energy E b . So, the thermal average used in Eq. 3.8 is defined as Now we can use these cross sections and the corresponding mean free paths to place

Bounds for Massive Dark Photon
Bulk Emission Trapping Figure 5. These plots show the lower constraint from bulk emission and the upper constraint from trapping as a function of particle mass. The first plot shows these constraints for a massive axion and massless dark photon and the second plot shows them for a massive dark photon and a massless axion. Note that the upper trapping limit only contains the contributions from trapped axions and transverse dark photons since we must treat the longitudinal dark photons separately (See Sec. 4).
the constraint on the coupling. The radii r a and r D , are implicitly defined as functions of the coupling using Eq. 3.1 and using the mean free paths defined in Eq. 3.8. Then we can numerically solve for the value of f a such that the emitted luminosity given in Eq. 3.4 is equal Since we expect both radii to decrease with increasing f a , and that the blackbody luminosity should increase as the radii decrease, we can exclude couplings below 1/f a (but larger than the constraint from bulk emission). The region of parameter space excluded by supernova is shown in Fig. 5. The first plot considers a massive axion and massless dark photon while the second considers a massless axion and massive dark photon. The high mass behavior of these constraints is quite different from the usual supernova constraint behavior. Ordinarily, the trapping constraint curve comes down to intersect the bulk emission curve to form a closed constrained region. However, we see that in our constraints, the region does not close off, but is an exponentially increasing strip. We will now look into the limiting behavior of the constraints in the high mass limit, M T , where M is the mass of the massive final state boson to determine the origin of this unique behavior.
Origin of the high mass scaling: The high mass behavior of the bulk emission curve is straightforward to understand. The initial particles must have enough energy to produce this very massive particle. This leads to the luminosity being Boltzmann suppressed for M T . The constraint curve is found by demanding L bulk (f ) = L ν so in order to compensate for the exponentially decreasing luminosity, the coupling must increase in exactly the same way. This leads to the exponentially increasing constraint curve for the bulk emission shown in Fig. 5. This type of scaling is expected in supernova constraints for any novel massive particle.
The trapping constraint curve is a bit more difficult to understand. To understand the scaling at high mass, first consider a more typical model that contains only one novel particle radiating like a blackbody. For this type of model, the constraint decreases with mass M for M T due to the fact that the radiation of a blackbody decreases as the particles mass increases. In order to keep L trapping (r, M ) = L ν the radius of the blackbody must decrease to compensate for this mass suppression. This smaller radius then implies that the constrained coupling is weaker and gives a decreasing coupling constraint in the high mass limit.
Contrast this behavior with a model, like ours, that contains two novel particles. As one of the particles gains mass, M , it is still true that the blackbody emission of that particle gets suppressed. However, the blackbody emission of the massless particle does not get suppressed. So in the high mass limit, the total luminosity is only due to blackbody radiation of this massless particle. Demanding that L tot = L ν then fixes the radius of emission to a radius r 0 , which is independent of the heavy mass M . Now, the trapping constraint curve is found by demanding that τ (r 0 ) = 2/3. This optical depth is determined by the strength of the scattering of the massless particle into the massive particle (with mass M ). This up-scattering process will be Boltzmann suppressed for M T and so to compensate, the coupling must increase in exactly the same way. This gives the exponentially increasing constraint curve that we observe in Fig. 5.
We will now investigate in more detail the source of these two exponentially increasing curves to determine the shape of their fall off. In particular, we will show that the bulk emission constraint scales exponentially in M/T while the trapping constraint scales exponentially in M 2 /T 2 . To do this we can look at the scaling of each constraint in the limit M T ignoring O(1) factors. For both curves, we found that this scaling is the inverse of the Boltzmann suppression of the initial particles in each interaction. Let us start with the bulk emission constraint. The dominant process here is the annihilation process. The electrons and positrons are largely relativistic so that we can ignore their mass m e . A simple computation gives, For the trapping constraint, the dominant scattering is the scattering off of electrons. In order for this scattering to be kinematically allowed, s > (M + m e ) 2 . Again, we ignore m e to find, Where E e − is the energy of the initial electron and E b is the energy of the initial, massless boson. Since E b is fixed to be the thermal average energy, E b ∼ T , we conclude that This leads to a Boltzmann suppression of the scattering that scales as e −M 2 /T 2 . The scaling of the trapping constraint curve is thus Since the trapping curve scales as e M 2 /T 2 and the bulk emission curve scales as e M/T the two curves never meet and diverge from one another. Note that if E b is not fixed to be the averaged energy, then we are free to let both E b ∼ M and E e − ∼ M so that our Boltzmann suppression will look like e −M/T . This will lead to a trapping constraint that scales as e M/T for M T . This very different asymptotic scaling of the trapping constraint highlights how important it is to use the averaged initial energy rather than averaging over the initial energies of the boson.

Longitudinal Dark Photons
A massive vector boson like the dark photon has a longitudinal polarization which behaves quite differently from its transverse counterparts. In particular, as we will show, for small dark photon mass, the coupling of the longitudinal dark photon is highly suppressed. This difference in coupling motivates us to consider the transverse and longitudinal dark photons as two separate particles. This allows us to get two constrained regions by applying the constraints from Sec. 2 and Sec. 3 to each polarization separately.
Let us now consider longitudinal dark photons in the small m D limit. We will show that the coupling of the longitudinal dark photons gets suppressed by a factor of its mass, m D . This is a result of the Goldstone Boson Equivalence Theorem. To see this, consider leaving the unitary gauge of our model through a gauge transformation The mass term for the dark photon gives rise to a kinetic term for the Goldstone boson, ζ, and a 2-point interaction/mixing term with the dark photon The Goldstone Boson Equivalence Theorem states that in the low mass/high energy limit, the amplitude for a matrix element containing a longitudinal dark photon is equal to the amplitude for the same process with the longitudinal dark photon replaced with the Goldstone boson. In the low mass/high energy limit the mixing term m D A µ D ∂ µ ζ can be treated as an interaction term with coupling m D . Since this is the only interaction the Goldstone boson has, we know that any matrix element with an external Goldstone boson must scale as m D . So for masses with m D T , we should find M L ∼ m D where M L is any matrix element involving a longitudinal dark photon in the initial or final state. This argument is illustrated pictorially in Fig. 6. Figure 6. A pictorial representation of the suppression of matrix elements containing longitudinal dark photons in the low mass/high energy limit. The first equality is the Goldstone boson equivalence theorem stating that in the low mass/high energy limit any matrix element with an external longitudinal photon is equal to the same matrix element with the longitudinal dark photon replaced with the Goldstone boson. The second equality makes use of the 2 point Goldstone bosons-dark photon interaction term that results from leaving unitarity gauge. This coupling scales as m D which gives the desired result. Now let us compute the constraints coming from the longitudinal dark photon. We start with the bulk emission constraint from Sec. 2. In order to compute this bulk emission constraint we need to compute the emissivity of longitudinal dark photonṡ This is exactly the same as Eq. 2.2, but rather than summing over all dark photons spins in the final state, we simply insert the longitudinal polarization for the dark photon (where longitudinal is defined in the frame of the star). We also are only accounting for the energy emitted by the longitudinally polarized dark photons, since we will assume that the axion is trapped in the supernova for the couplings of interest. We will consider two process which produce longitudinal dark photons. First, we will consider the annihilation process, since it was the dominant process in the full emissivity. Second, we consider the plasmon decay process because the mass suppression has a smaller impact on this process compared to the others. A general computation given in App. A shows that any production process will have a matrix element that scales as m 2 D /s in accordance with the Goldstone boson equivalence theorem. For the plasmon decay, the center of mass energy in this process is the plasmon mass which is of the same order as the plasmon frequency, ω p . The plasmon frequency is around a full order of magnitude smaller than the typical center of mass energy for annihilation so the m 2 D /s suppression factor leads to the plasmon decay emissivity being significantly less suppressed by the mass of the dark photon than the other processes. In fact after computation, we find that plasmon decay is the dominant process in the production of longitudinal dark photons  Just as before Q a L and Q p L are O(1) numbers for low m D and m a and are defined in App. A. Comparing Eq. 4.4 to the emissivities given in Eq. 2.5 and Eq. 2.7 for the annihilation and plasmon decay emissivities, we can see that they have the same form but with an additional m 2 D /s-type suppression factor in both.
For the longitudinal trapping bounds we will need to compute the cross section for longitudinal dark photons scattering with fermions, as well as their decay width. The decay width for the longitudinal dark photon can be seen to be the same as the decay for the transverse dark photon by considering the decay in the dark photon's rest frame. Thus the decay width is the same as that given in Eq. 3.7.
Next we turn to the scattering cross section for longitudinal dark photons. The matrix element for this process is the same as the second diagram in Fig. 1, with the polarization of the initial dark photon chosen to be longitudinal. The computation is very similar to that of the full scattering cross section and is detailed in App. A. The result is .

Other Astrophysical Constraints
As a final exercise, we can improve on our constraints, especially in the low mass regions by using our cross section computations in other astrophysical situations. We will examine white dwarf cooling bounds and BBN N ef f bounds.
White Dwarfs First we consider white dwarfs. A cooling bound on the coupling, similar to that found for SN1987A, can be derived from white dwarf data. For low temperature white dwarfs, the dominant source of energy emission is the emission of surface photons via blackbody radiation. This, along with the degeneracy of the white dwarf, allows one to derive a simple power law relationship between the luminosity of a white dwarf and the relative abundance of white dwarfs with that luminosity in the universe known as "Mestel's Law". For higher luminosity white dwarfs, bulk neutrino emission begins to play a significant role in cooling the white dwarfs resulting in a dip in the white dwarf abundance at higher luminosity. Observations of white dwarf abundance agree well with Mestel's Law and the socalled neutrino dip. An additional channel of significant bulk emission would increase this dip spoiling this agreement. It has been found that any additional mode of bulk particle emission must have an emissivity less than that of the neutrinos in order to preserve experimental agreement [59]. Therefore, the cooling constraint for white dwarfs can be summarized aṡ For neutrino emission and our axion/dark photon emission, plasmon decay is the dominant source of particle production in white dwarfs and thus will be the only process considered.
Ref. [6] has computed a nice expression for the emissivity of Standard Model neutrinos This shares a similar form to our expression for the emissivity of axions and dark photons, given in Eq. 2.7. In our constraints we take the white dwarf to have a constant temperature (T = 2 keV) and constant plasmon frequency (ω p = 40 keV) as is typical for white dwarf constraints [6]. Applying Eq. 2.7 and Eq. 5.2 to the constraint given in Eq. 5.1 we get a constraint on the coupling in agreement with previous considerations of this bound [41,48]. At low masses Q 3 ≈ Q p , but as mass of the axion or dark photon increases Q p decreases. An upper trapping bound on the excluded region from white dwarfs could be obtained but is rendered irrelevant by even the most conservative BBN ∆N ef f bound.
BBN Bounds Now we consider a very simple, conservative BBN ∆N ef f bound. In our simple constraint we demand that axions and dark photons were not in equilibrium at nucleosynthesis as the addition of 3 (or 4 for massive dark photons) bosonic degrees freedom in the radiation density would induce too large of an increase in the number of effective neutrinos. We take the temperature at nucleosynthesis to be T BBN = 2 MeV, both m a and m D to be smaller than T BBN , and consider pair production of axions and dark photons through electron positron annihilation. To implement ∆N ef f bounds, we demand that where n e is the number density of electrons and n a,γ D is meant to convey the equilibrium density of whichever particle is massless. Neglecting the electron mass since T BBN > m e , we find v mol σ e + e − →aγ D = α e 24f 2 a .

(5.5)
Inserting this with the densities and Hubble at 2 MeV, we get that the allowed couplings must satisfy, 1/f a < 4.5 × 10 −7 GeV −1 for massless axions and 1/f a < 6.5 × 10 −7 GeV −1 for massless dark photons. This constraint applies until the mass of the axion or dark photon is Ref. [60]. Note that these constraints are only shown in the massive dark photon plot since Ref. [60] only considers the m D m a limit. Similar collider constraints exist for the opposite limit but have not yet been computed. Also note that slightly stronger bounds than the white dwarf cooling bounds can be obtained by using the recent analysis of the tip of the red giant branch [43,61]. around 1 MeV at which point they no longer contribute to the radiation density.
These two constraints as well as the supernova and collider constraints taken from Ref. [60] are shown in Fig. 8. Note that only the constraints for the massive dark photon show collider constraints since Ref. [60] derive their constraints explicitly in the limit m D m a . Similar constraints could be derived in the limit m a m D although this has not been done. We expect that collider constraints for m a m D will exclude a similar region of the parameter space as those for m D m a , but have refrained from including it in Fig. 8.

Conclusion
In this article we considered the supernova bounds on the well motivated topological axionphoton-dark photon coupling. This coupling is interesting because it doesn't require or generate kinetic mixing between the vectors, while at the same time it can be responsible for generating the abundance of dark matter. The supernova bound on this coupling highlighted two generic features that appear in supernova constraints. The first feature occurs when the interaction being probed at the supernova involves multiple new particles with different masses, an example familiar from similar to the dynamics in dark matter models involving co-annihilation [62] or co-scattering [63]. When this situation occurs, the production and scattering cross section decreases whenever one of the particles obtains a mass so that the bulk emission and trapping bounds never meet and instead both decrease exponentially.
The second feature is a bit more ubiquitous and involves the importance of treating longitudinal and transverse polarizations of particles independently, with their own cooling and trapping constraints. Because the longitudinal mode coupling is m D /E suppressed in the low mass limit, it acts as its own weakly coupled particle that can cool supernova giving a cooling and trapping bound separate from the transverse modes.
Because these properties are generic, they will appear in all manners of motivated scenarios. The most obvious model where a new feature occurs is in dark photon models with kinetic mixing. If the dark photon is heavier than twice the electron mass, trapping is dictated by decay rather than scatter. Since decay treats all polarizations the same, the differences between transverse and longitudinal modes doesn't become relevant until the dark photon mass is less than an MeV, where supernova bounds are superseded by other constraints.

A Computational Details
In this appendix, we provide details on how to calculate the averaged cross-sections used in the constraints. First, we define some notation to simplify the calculations. We will be frequently integrating over final states, for which we will use the shortened notation We use the bracket notation to refer to a thermal average over initial states A and B. For some general function g(p A , p B ) of these momenta, this average is defined as In order to capture kinematic features of most of our processes we define two functions α and β as follows, Physically, in the center of mass frame of two particles with center of mass energy √ s, the energy and momentum of the first particle are E 1 = √ s 2 α(s, m 1 , m 2 ) and |p| = √ s 2 β(s, m 1 , m 2 ), with similar expressions for the second particle.

A.1 Dark Photon and Axion Production Emissivities
In this section we give the details of the three emissivities in Eq. 2.5 and the emissivity in Eq. 2.7. We start with Eq. 2.2. As a convention, in all processes, the initial momenta will be indexed by capital letters while the final momenta will be indexed with a number. p 1 and p 2 represent the axion and dark photon momentum respectively. k will be the momentum of the photon that produces the axion and dark photon and q will be the total momentum. We will see that the matrix element for all processes can be factorized as The label X on M X µ will label the process, taking values A for annihilation, B for bremsstrahlung, C for Compton and P for plasmon decay. The factor A µ contains all of the matrix element's dependence on the axion and dark photon final states. Writing our matrix element in this way allows us to factor out our final state phase space integral into two and write Eq. 2.2 as two separate integrals, We now have one factor that contains all of the process-dependent information contained in M X µ and a second piece that contains all of the dependence on the axion and dark photon momentum. This second piece can be computed immediately, We know that for every process (except plasmon decay) k µ M X µ = 0. So the first term doesn't contribute and we get a simple expression for the emissivity of 2 → N processeṡ In each process, we must identify the correct portion of the matrix element M X µ and insert it into this equation.
Annihilation Following the first diagram from Fig. 1, the matrix element for the annihilation process is Comparing this with Eq. A.4, we can identify, Finally we define the factor so that we may simply write˙ a = α e 24f 2 a n e + n e − E tot Q a , It is worth nothing that in the limit m D = m a = m e = 0, Q a = 1. So, given that electrons and positrons are dominantly relativistic in the supernova, Q a is O(1) for small m a and m D .
Nuclear Bremsstrahlung In the soft radiation approximation, there are two diagrams which contribute to the amplitude for bremsstrahlung production: In order to write down the matrix element for this process, we define a Dirac matrix Γ pn by so that the amplitude for proton-neutron scattering can be written as We can then write down the matrix element for our bremsstrahlung process So, we can identify the portion of the matrix element M µ B as Because we are working in the soft radiation approximation we take |k| |p|, where p is any of the nucleon momenta in the process. With that we get a simplified expression for M µ B , We can then easily square this and sum over the remaining final spins and average over the initial spins. Using the fact that the nuclei are non-relativistic and expanding in their velocities we find Next, we insert this into Eq. A.6, The remaining phase space integral will contain a delta function δ 4 (p A + p B − p 3 − p 4 − k). Once again, we use the soft photon approximation to simplify and we ignore the photon momentum k in the delta function. However, we must ensure that the photon energy is less than the center of mass energy of the nuclei, and so we impose, This condition could have been achieved with a hard cutoff on the photon energy [14], however the exponential simplifies the computation and leads to similar results (see e.g. [16,58]). This substitution also allows us to factor the emissivity into two separate factors, (A.13) The first factor can be seen to be proportional to the momentum transfer cross section for pn → pn scattering, where, The second factor in Eq. A.14 can be simplified significantly by performing the angular integrals and changing integration variables from k 2 and k to k 0 and v = |k|/k 0 . The result is, This integral can be computed exactly in the massless dark photon and axion limit, With this in mind we can write this second factor as αeT 3 Compton There are two relevant diagrams for the Compton scattering process, The matrix element is So we can identify the piece M C µ as Using this matrix element, and after some algebra we finḋ where t = l 2 , and After doing the integral over final states dΠ 2 , the end result iṡ where, The integral over k 2 runs from (m a + m D ) 2 to ( √ s − m e ) 2 . Unlike the other two processes, the m a = m D = 0 limit does not lead to a simple expression. However, we may define Q c is still roughly O(1) in the limit m a = m D = 0 since E A ≈ E tot /2.
Plasmon Decay From the diagram in Fig. 2, we can easily write the matrix element for plasmon decay as: Where P µ is the plasmon polarization. If we square this and sum/average over final/initial spins and integrate over final states using the results of Eq. A.5, we can get We will compute the matrix element of longitudinally polarized plasmons and transversely polarized plasmons separately since their dispersion relations and field normalizations are different. The polarizations and dispersion relations are taken from Ref. [56]. The polarizations are Where T · k = 0. The dispersion relations are found by numerically solving, Where, These analytic equations for the dispersion relations are derived by assuming that the velocity of the particles in the plasma are dominated by single velocity, v * given by v * This also defines the plasma frequency ω p . Finally, the field renormalizations for the plasmon are given by .
Using all of this, our longitudinal and transverse matrix elements are Then for each polarization We can write the emissivity using Eq. 2.4 pπ t , m a , m D )

A.2 Scattering Cross Section
The matrix elements for each process are shown in Fig. 4 and their amplitudes are simple to write, using notation developed in the previous sections.
Q ψ is the electric charge of the fermion, ψ. We can easily square this matrix element, sum over final spins and average over initial state spins: g ,k represents the degrees of freedom in the initial particle that is either the axion or dark photon and t is the Mandelstam variable equal to the square of the transferred momentum. We will also use k to denote the initial particle four momentum and k to denote the final particles four momentum regardless of if that particle is an axion or dark photon. Framing the computation this way allows us to compute both matrix elements at once. As discussed in Sec. 3, we are interested in computing the momentum transfer cross section defined in Eq. 3.5. For computational ease we consider a relativistically invariant generalization of the momentum transfer cross section by replacing the 1 − cos(θ) with 2t/(sA). Where and It is easy to see that in the limit m ψ /s 1 (valid for electrons and positrons which are relativistic in the supernova) and the limit m 2 ψ k 2 , k 2 (valid for protons) This property allows us to compute the momentum transfer cross section without sacrificing relativistic invariance.
We have checked that the final results are insensitive to this particular deviation from standard convention.
In order to compute the cross section will need to compute the final state integral This integral can be computed by expanding in terms of Lorentz covariant tensors involving the initial state momenta and metric. We can then compute the necessary coefficients in this expansion by picking specific components for the indices in the center of mass frame or contracting. When all is done, the result is, Finally, plugging this into the cross section, we find Inserting m a and m D into the appropriate spots for each cross section gives us Eq 3.6.

A.3 Longitudinal Scattering Cross Section
In order to compute the longitudinal bounds, we need to compute the cross section for γ L D e − → ae − . The method described above transfers nicely to this specific case. Starting with Eq. A.25, we replace the sum over dark photon polarizations with just the longitudinal polarizations. Since the dark photon is in the initial state here, these polarization don't interfere with the final state integrals which are expanded and computed just as before. The only remaining difference is the contraction with the epsilon tensors which now must be contracted with the explicit longitudinal polarization rather than the metric. After some simplification, the result is

A.4 Decay Widths
Let us now compute Eq. 3.7 From the diagram given in Fig. 4 it is easy to see that the amplitude for scattering can be written as iM = i µ γ A µ (p 1 , p 2 ).
Here, p 1 will be the photon momentum and p 2 will be the other final state particle, axion or dark photon. Squaring this, summing over final spins, averaging over initial spins, and integrating over final states is straight forward. To find the decay width we simply divide by 2E where E is the intial energy and we recover Eq. 3.7.

A.5 Production of Longitudinal Dark Photons
The emissivity for a given process of longitudinal dark photons is given by Eq. 4.3. Just as before, we separate the matrix element as For the two processes we are considering here, the final state is only the axion and dark photon so we can writė Where i n i is the product of the densities of the initial state(s). The final state integral can be computed by considering the Lorentz covariant tensor L µνα ≡ dΠ 2 (p 1 , p 2 )A µ L A ν L p α 2 of which we wish to take the α = 0 component. This integral can be computed by contracting the µ and ν indices, computing in the supernova frame where the longitudinal polarization takes the simple form L = 1 m D |k|,kk 0 , and then properly relating the supernova frame momenta to the center of mass frame momenta so that we may do the final state integral in the center of mass frame. The result is that,

B Supernova Profiles
In this section we discuss the profile used in placing our constraints and the variation in the constraints placed that arise when using a different profile. The temperature (T ), mass density (ρ), and electron fraction (Y e ) profiles used in our computations use analytical fits to simulation used in Ref. [58]. These profiles are shown in Fig. 9.  Our expressions for the emissivities and scattering cross sections involve number densities and chemical potentials. From these profiles we can calculate the number densities for the electron, proton, neutron, photon and positron. From the number densities, we can compute the chemical potentials of the electron, proton and neutron (µ e , µ p and µ n respectively).
In order to access the dependence of our results on the profile used, we redid our calculations using the 18 M profile found in Ref. [64]. To compare the results visually, we show the total luminosity as a function of axion mass for the two different profiles in Fig. 10. In both profiles, annihilation is the dominant source of cooling. The difference between the bulk cooling bound placed on f a when using the two different profiles is about a factor of 1.5, which is within the expected uncertainty of cooling bounds from supernovae.  Figure 10. The total luminosity of various channels as a function of axion mass. The solid lines show the profile in Ref. [58] while the dashed lines show the profile in Ref. [64]. The difference in the bounds on f a is about a factor of 1.5.