Radiative transfer in stars by feebly interacting bosons

Starting from first principles, we study radiative transfer by new feebly-interacting bosons (FIBs) such as axions, axion-like particles (ALPs), dark photons, and others. Our key simplification is to include only boson emission or absorption (including decay), but not scattering between different modes of the radiation field. Based on a given distribution of temperature and FIB absorption rate in a star, we derive explicit volume-integral expressions for the boson luminosity, reaching from the free-streaming to the strong-trapping limit. The latter is seen explicitly to correspond to quasi-thermal emission from a"FIB sphere"according to the Stefan-Boltzmann law. Our results supersede expressions and approximations found in the recent literature on FIB emission from a supernova core and, for radiatively unstable FIBs, provide explicit expressions for the nonlocal ("ballistic") transfer of energy recently discussed in horizontal-branch stars.


Introduction
Dark sectors arising from physics beyond the standard model could provide explanations for various shortcomings of the standard model itself, including dark matter, neutrino masses, the baryon asymmetry, and the strong CP problem. One typical phenomenological consequence is the appearance of new, feebly-interacting bosons (FIBs) that can be experimentally searched and astrophysically or cosmologically constrained. One class of traditional arguments uses observational consequences of FIB emission from stars, an idea independently advanced by several groups in 1978 [1][2][3][4] when the Weinberg-Wilczek axion had been recognized as a consequence of the Peccei-Quinn solution of the strong CP problem. Ever since, the impact of many types of bosons in various astrophysical systems has been studied [5], sometimes posing interesting conceptual questions about FIB production or propagation in stars.
We here follow up one such case that has emerged in several recent studies of FIB production in supernova (SN) cores [6][7][8][9][10][11][12]. Actually the feeble interaction was taken strong enough to prevent free escape after production. In analogy to the SN "neutrino sphere," the FIBs emerge from a decoupling region that is traditionally pictured approximately as a black surface for thermal FIB radiation according to the Stefan-Boltzmann (SB) law [13]. The relevant temperature T SB is taken to be that of the SN medium at a radius R SB where the FIB optical depth is 2/3, and the radiating surface is 4πR 2 SB . We will see that this prescription is rather accurate, as physically it should be, but has evoked some doubts because clearly there is no hard surface of emission-the radiation must come from a shell with a geometric thickness corresponding to optical depth of around one.
Motivated by this question and doubts in the recent literature we take a fresh look at radiative transfer by FIBs that may or may not have a significant mass. In the diffusion limit, this problem was formulated a long time ago [14], following the standard theory of radiative transfer by photons. 1 Our focus here is to study explicitly the transition between the freestreaming and trapping (diffusion) limits, both in plane-parallel and spherical geometry. The latter is particularly interesting in a situation when the FIBs are unstable and deposit energy in regions far away from the compact emission volume, i.e., in a situation where the geometric extension of the "stellar atmosphere" is not much smaller, or even much larger, than the core radius of a SN or a horizontal-branch or red-giant star [20].
The main simplifying assumption, motivated by the boson interaction being "feeble", is to include only FIB absorption and emission from a medium in local thermal equilibrium, but not scattering between different FIB momenta or annihilation. In this case the only particlephysics ingredient is the "reduced absorption rate" Γ ω as a function of FIB energy ω, where Γ ω is equivalent to the imaginary part of the FIB self-energy, that also depends on the local conditions of the medium such as temperature, density, and chemical composition. In the absence of scattering, the stationary FIB occupation number on a given ray, corresponding to a given mode k of the FIB radiation field, can be expressed as an integral along this ray. Global solutions for plane-parallel or spherical geometries then follow as suitable superpositions of such single-ray solutions. In other words, for a given stationary stellar background model, the FIB radiation field is found from a quadrature. Explicit volume-integral expressions, notably in spherical geometry, are the main technical results of our paper. Based on Γ ω (r) and T (r) as functions of stellar radius, we thus provide integral expressions for the FIB luminosity L ω (r). Taken at spatial infinity, dω L ω (∞) provides the total FIB luminosity, e.g., of a SN core. Moreover, one can find the energy loss or deposition at a given radius through the radial variation dL ω (r)/dr.
Solutions derived from a prescribed and stationary background model are only useful, of course, in a physical situation when the thermal timescale exceeds the dynamical one. If this is not the case, and if the diffusion limit does not apply, the full Boltzmann collision equation needs to be solved, a task that is of course the main numerical effort in core-collapse SN simulations concerning neutrino transport.
Radiative transport by neutrinos, despite their weak interaction, is a much more complicated task than our FIB treatment. Neutrinos and antineutrinos of the electron and muon flavor can be absorbed and emitted by the medium through charged-current interactions, but neutral-current scatterings as well as annihilation and pair emission and absorption through bremsstrahlung and other processes occur on the same order of the coupling constant G 2 F . Moreover, besides energy also lepton number of different flavors is being transported.
In principle, our exercises are straightforward, but the devil is in the details, even for the much simpler problem of FIB transport. The correct expressions are apparently not available in the literature (and incorrect expressions or approximations have been floated), justifying our derivations, at the risk of being seen as a pedagogical exercise of standard radiativetransfer theory. In the same vein we also show explicitly the transition between a volume integral and a quasi-thermal surface integral in the strong-trapping limit. We believe that deriving these results from first principles, starting with the Boltzmann collision equation, is an instructive exercise that offers many interesting insights that may be useful for future studies of astrophysical particle bounds.

Radiative transfer by feebly interacting bosons
We begin with the Boltzmann collision equation (BCE) for new bosons a (reminiscent of "axion") that can be produced, for example, by processes of the type γ + B → B + a, that is to say axion-photon conversion by interaction with fermions (for example semi-Compton scattering on electrons or muons) or other charged particles as in the Primakoff case, but photon coalescence 2γ → a is also conceivable. On the other hand, scattering of the type a+B → B +a plays no role because the interaction is much more feeble than that of photons.

Freeze out from first principles
Ignoring FIB scattering from one momentum mode to another, we can focus on the evolution of a single mode with energy ω along some ray with spatial coordinate x. The BCE for the occupation number f is in this case where v is the particle velocity. Here Γ E is the spontaneous emission rate that appears multiplied with the boson stimulation factor 1 + f , whereas Γ A is the absorption rate, and in general both depend on ω and x. In the second expression, the terms proportional to f were consolidated and are proportional to the "reduced absorption rate" Γ * A = Γ A − Γ E that includes the effect of stimulated emission as a negative absorption rate.
If the medium is in local thermal equilibrium, detailed balance implies that locally Γ E = e −ω/T Γ A so that the reduced absorption rate is which we use as the absorption rate and which is the quantity that defines the optical depth. The spontaneous emission rate is then expressed as a relation between emission and absorption corresponding to Kirchhoff's Law.
In a stationary and homogeneous situation, the left-hand side (LHS) of Eq. (2.1) vanishes and the equation is solved by a thermal Bose-Einstein distribution f eq = (e ω/T − 1) −1 . So we may write the BCE instead for the deviation from equilibrium ∆f = f − f eq in the form So it is the reduced absorption rate Γ which damps the deviation of f from equilibrium, explaining its central importance for radiative transfer.
In the context of thermal field theory, the boson propagation properties are encoded in its self-energy Π within the medium. The imaginary part provides the rate-of-approach to thermal equilibrium as Im Π = −ωΓ [21], once more highlighting the role of the reduced absorption rate as the central interaction parameter.

Stationary state
Our main interest is a stationary situation, so only the gradient term on the LHS of the BCE survives and we need to solve v df dx = Γ E − Γf, (2.5) where the spontaneous emission rate Γ E is given in Eq. (2.3) in terms of the reduced absorption rate Γ under the assumption of local thermal equilibrium.
To solve this equation we notice that T and Γ are functions of x and we define the optical depth as where λ is the mean free path (MFP) for a FIB with velocity v. For massless FIBs we should use v = c = 1 everywhere. The optical depth τ (x) is measured relative to a distant observer at x = +∞ where τ (∞) = 0. So finally one finds the solution This is the intuitive answer that the occupation number at x is filled by spontaneous production up to this point, reduced by the absorption along the path from production to detection.
Here it was assumed that no radiation enters at the boundary at x = −∞, i.e., all radiation is generated by emission within the realm of integration.
Instead of x we may use τ (x) itself as a coordinate along the beam. Notice that this is a monotonically decreasing function of x and thus uniquely invertible to provide x(τ ). The limiting values are τ (∞) = 0 and reaches a maximum value τ max = τ (−∞). Notice also that dτ (x)/dx = −Γ(x)/v = −1/λ(x) and we introduce the blackbody occupation number at τ for the local temperature T (τ ) We see that the solution depends only on the temperature profile T (τ ) along the ray. The velocity v no longer appears explicitly because the optical depth is based on λ and not on Γ. If the medium is very opaque so that we cannot see through the star to the other side we may use τ max = ∞. For the occupation number at spatial infinity, corresponding to τ = 0, one finds in this opaque limit In the special case when the medium has the same T everywhere this is simply f (0) = f eq , the Bose-Einstein occupation number. So an optically thick object at temperature T radiates bosons with a thermal Bose-Einstein distribution. However, even if the radiating body has a hard material surface, the photons do not emerge from that surface but from a layer with thickness of a few MFPs.
We also consider the occupation number f − (τ ) of the "backward mode" moving in the opposite direction, toward the star, where it was assumed that at spatial infinity (τ = 0) the backward mode is not occupied. In this notation, the occupation f (τ ) of the outgoing mode is termed f + (τ ).
In this discussion we have implicitly assumed that the FIB absorption rate Γ depends on the background medium which is geometrically bounded so that it makes sense to use a distant observer as a point of reference when using the optical depth as a measure of distance. However, when FIB decay of the form a → 2γ is important, this approach is not justified. We will return to this question in the context of our spherically symmetric solution.

Particle flux
The radiation emerging from a source is usually not described in terms of the occupation numbers of the modes of the radiation field but rather by the corresponding particle or energy flux. The net particle flux in the outgoing direction is whereas the energy flux sports an additional factor ω. Assuming no backward occupation at spatial infinity, the outgoing flux for a distant observer is simply φ(0) = vf (0) given in Eq. (2.10). At intermediate positions, the flux can be expressed as So we find the intuitive result that the flux along some ray is driven by the temperature profile a few MFPs up-and downstream from the point of interest. Formally the function φ(τ ) on the interval 0 ≤ τ < ∞ is a certain linear transformation of the function f eq (τ ) on that same interval.

Example with power-law profile
We can illustrate FIB freeze out with a T profile inspired by a realistic Proto Neutron Star (PNS) profile of the form T (τ ) = T 1 τ p , (2.14) where 0 < p 1 is a small number for which we use p = 1/5 and T 1 is the temperature at unit optical depth. Moreover, we assume the FIB to be massless so that v = 1. For a typical boson energy of ω = 3T 1 we show the solutions f ± as well as the flux φ = f + − f − in Fig. 1. We see that the flux escaping from the star corresponds approximately to f eq at τ 0.8, but at this location the actual solution is far away from this value. The approach to the asymptotic solution is slow in the decoupling region.

Radius [km]
Occupation Number Figure 1. Solutions for the occupation numbers f ± and the flux φ for a massless boson and using our power-law profile for the temperature and using a typical energy ω = 3T 1 . The horizontal black line projects the escaping flux (the occupation number at the stellar surface) with the equilibrium one and marks the Stefan-Boltzmann optical depth for this energy, here approximately at τ 0.8 For illustration we can also go back to coordinate space and show these results as a function of geometric radius. We find it useful to take inspiration from a realistic model of a SN core, following in particular the Garching group's muonic model SFHo-18.8 [9,22] that we employed earlier for other studies [10,11]. In this case, one can see that the temperature varies approximately as T = T 1 (r 1 /r) 4 which is equivalent to T = T 1 τ 1/5 , implying τ = (r 1 /r) 20 where we take r 1 = 17 km. In this representation, the approach to the asymptotic solution looks more intuitive, but it remains true that the approach to the asymptotic solution does not happen at the nominal decoupling radius, but is considerably smeared out even though here we have a fixed energy and no energy dependence of the cross section.
So the picture that a Stefan-Boltzmann flux emerges from some narrow geometric range like "surface emission" is clearly not accurate. The bosons reaching infinity derive from a broad radial range, equivalent to a broad range of optical depth.

Strong trapping regime and plane-parallel atmosphere
The single-ray solutions of the previous section provide the full answer to the question of the stationary FIB radiation field based on a source distribution with prescribed properties (no feedback effects by particle emission on the medium). It remains to cast this result into a more explicit form for relevant overall geometries. To discuss more explicitly radiation decoupling in the strong trapping limit we use a plane-parallel atmosphere, where the temperature and optical depth are only functions of a cartesian coordinate z perpendicular to the atmospheric layering. In the example shown in Fig. 1, inspired by a realistic SN core model, the decoupling radius is some 17 km and the relevant shell has a thickness of a few km, so the plane-parallel approximation should provide a reasonable first description.

Intensities vs. occupation numbers
Solving the Boltzmann collision equation was most transparent using occupation numbers which appear directly in Bose stimulation factors. However, in the end we ask for the energy flux at some radial position. In this spirit we turn from occupation numbers to radiation intensities for a mode k of the radiation field where ω = (m 2 a + k 2 ) 1/2 . Notice that |k| = vω where v is the boson velocity. We have normalized the intensity such that the integral over energy and directions I k dω dΩ/4π is the local energy density. Whether or not to include the factor of 4π in the definitions of I k and the blackbody intensity B ω in Eq. (3.2) is a matter of convenience.
When the FIBs are in thermal equilibrium, the occupation number is f k = 1/(e ω/T − 1) and the equilibrium (blackbody) intensity is denoted as Here B ω is the blackbody intensity for one massless degree of freedom and v ω = 1 − m 2 a /ω 2 is the velocity for a boson with mass m a . For the massless case, the total energy density is For a nonvanishing mass, no simple expression exists.

Angular moments
The previous single-ray solution applies to a mode propagating in the radial direction, but now we consider one that is inclined by µ = cos θ such that µ = +1 is the outward direction and µ = −1 the inward one. We begin with Eq. (2.7) for the occupation at position x = z/ cos θ along the ray. As variable of integration we may use z, so we use the vertical depth as a measure of propagation distance, or equivalently, the optical depth τ in the vertical direction. Following the previous steps we find for the outgoing and incoming intensities where µ > 0 by definition, i.e., we treat inward-moving modes explicitly as backward moving ones with positive µ.
We are mostly interested in the energy flux, but in general one defines angular moments of the type Traditionally the zeroth moment (the energy density) is called J ω , the first moment (the energy flux) H ω , and the second moment K ω is related to pressure. For photons v = 1 and B ω acquires a factor of 2 for the two polarization states. The factors of v are understood in the sense that a flux (of energy or particles) requires one power of v compared with the massless case, whereas the pressure, being essentially a flux of momenta, requires one more v factor. Indeed, the spatial part of the stress-energy tensor dimensionally involves (momentum) 2 . The angle integrations in Eq. (3.4) can be performed explicitly. For the n th moment and using w = 1/µ one finds We use E m (t) only for positive arguments where it is positive and real. To consolidate the ± cases in Eq. (3.5) in a single expression it is convenient to define integral kernels of the form shown in Fig. 2. These are even functions of t for even n and odd functions of t for odd n.
With this notation, the moments of Eq. (3.5) are Notice that one factor of v ω comes from B v ω for particles with mass, whereas B ω is the massless intensity and thus only a property of the medium profile, not the particle mass. In the massless case, these are the Schwarzschild-Milne equations, providing us with the moments of the radiation field as linear transformations of the blackbody intensity on the interval 0 ≤ τ < ∞. The 0 th -order case, providing the local energy density, is known as the Λ-transformation. Figure 2. The n th -order integral kernels E n (s) defined in Eq. (3.7).

Diffusion regime
Asymptotically E n (t) = t −1 e −t for t → ∞ independently of n. Among other consequences, this implies that integrals over any power t n weighted with such kernels converge. It also implies that the local values of the moments only depend on the thermal radiation field a few MFPs up-and downstream. In particular, we consider a general function b(t) that we expand as a Taylor series Then we find which for the first three momenta gives explicitly Of course, this representation makes only sense at large optical depth where the lower limit of integration can be extended to −∞ and the Taylor expansion is really around a point τ 1. In this case we see that the kernels for the first two momenta at leading order effectively act as So deeply in the trapped regime, many MFPs away from the surface, the net diffusive flux is which is driven by the temperature gradient. (We prefer the letter F to H that is traditional in the theory of radiative transfer.) Recall that z is the coordinate perpendicular to the plane-parallel atmospheric layering, that the MFP is λ ω = v ω /Γ ω with the particle velocity v ω , that the Jacobian is dz/dτ = −λ, and that a factor v 2 ω comes from the first factor in Eq. (3.8). The diffusive flux is a good representation of the true flux when the MFP is short compared with the length scale of temperature variation. However, we can formally define F diff (τ ) everywhere, whether or not it is a good approximation of the true flux.
Finally we can define the Stefan-Boltzmann (SB) flux, given by the equilibrium intensity at a given radial position times an average angular flux factor 1/2 and times another factor 1/2 to count only the outward going modes. F SB (τ ) is the hypothetical FIB flux produced by a black surface at the radial position τ with the local T (τ ). Of course, the SB flux is simply another way of expressing the local temperature. Overall we define three different fluxes For the diffusive flux, we have rediscovered the usual factor 1/3 following directly from the properties of the exponential integrals. We recall that B ω is the intensity for one massless boson degree of freedom.

Integration over energies for a grey atmosphere
We are usually not interested in the detailed energy dependence unless there are resonant effects. So we may integrate over energies, but this requires to specify the energy dependence of the FIB interaction rate. The assumption that the reduced absorption rate Γ does not depend on energy is called the "grey-atmosphere approximation" in the theory of radiative transfer. Moreover, we now consider massless particles with v = c = 1. The grey-atmosphere approximation is surprisingly well motivated by FIBs absorbed by the Primakoff process as detailed in Sec. 5.1. Here we simply use this approach for the purpose of illustration. The integrated blackbody energy density for a single massless boson degree of freedom was given in Eq. (3.3) as B(τ ) = (π 2 /30) T (τ ) 4 , where the optical depth does not depend on energy by assumption. Then our three fluxes are explicitly Besides overall coefficients, the SB flux is a purely local quantity, the diffusive flux a spatial derivative, whereas the true flux involves a nonlocal operator, a convolution over all space, in practice a few units of optical depth upstream and downstream. So these three fluxes are nicely systematic about the FIB flux in the trapping limit. For a distant observer (τ = 0) and inserting the definition of E 1 (t), the true flux is found to be where E 2 (t) is the second exponential integral. The interpretation is that of every boson launched isotropically at optical depth τ , the probability to escape to infinity is given by the transmittance T(τ ) = 1 2 E 2 (τ ). The factor 1/2 accounts for all bosons launched away from the surface cannot escape, whereas the others have a chance of escape of E 2 (τ ). If all bosons were emitted either exactly toward or exactly away from the surface, the transmittance would be 1 2 e −τ . Due to the angular average e −τ → E 2 (τ ). We recall that τ here means the optical depth counted directly inward from the surface. The functional form of 1 2 E 2 (τ ) is, for positive τ , the orange curve in Fig. 2 marked E 1 . For large arguments it is 1 2 e −τ /τ .

Example with power-law profile
For illustration we return to the power-law profile of Eq. (2.14). Apart from a global factor that we now leave out, the three fluxes are For a typical case p = 0.2 we show the fluxes as a function of optical depth and radius in Fig. 3. We see that the diffusive and true fluxes become asymptotically close at large optical depth and then separate in the freeze-out region. This is most intuitively clear in the radial plot. The required optical depth for the SB flux to match the escaping true flux is τ SB 0.60. Notice that p = 1/4 is a special value where F diff = 1/3 is a constant and F SB = τ /4 increases linearly. We have not used this value to avoid an overly special case. In general, the true flux at the surface (τ = 0) is explicitly where we have used an expansion around the special value of p = 1/4 where this flux is near to a minimum. The condition where the approximation is good on the few-percent level in the entire range 0 < p < 1. For our special value τ SB (1/4) = 2/3 is exact.
In the neutrino decoupling region of a SN core, when diffusive transport is still appropriate, the neutrino flux itself, driven by the temperature gradient, is approximately constant. Therefore, the radiation density of neutrinos scales roughly linearly with neutrino optical depth. As the neutrino scattering rate is proportional to the density as assumed for our FIBs, the power-law index p 1/4 for the temperature as a function of optical depth is well motivated in the neutrino decoupling region and borne out from numerical models.
We may also ask where the emitted radiation reaching a distant observer is actually emitted. Equation (3.16) implies a distribution proportional to E 2 (τ )B(τ ). For p = 1/5 and thus B ∝ τ 4/5 , the normalized source distribution is (3.20) shown in the left panel of Fig. 4. As a function of geometric radius once more we assume τ = (r 0 /r) 20 with r 0 = 17 km, leading to the normalized source distribution shown in the right panel of Fig. 4. The vertical dashed lines show the location of the Stefan-Boltzmann radius.  We learn from this figure that the bosons reaching infinity originate from a fairly thick shell, not a sharp "boson sphere." The grey-atmosphere model, without any energy dependence of the cross section, provides the "sharpest" conceivable emission sphere. For neutrinos, the cross section varies with the square of energy and the "neutrino sphere" is much more smeared out and energy dependent.

Constant plus linear profile for the radiation density
The special power-law profile T (τ ) ∝ τ 1/4 corresponds to a linear profile for the radiation density B(τ ) ∝ τ . The next simple profile derives from adding an arbitrary constant where the letter q is traditionally used. The true flux is found through the convolution of Eq. (3.15a), leading to a complicated expression in terms of exponential integrals. At the surface (τ = 0) one finds the following true flux to be compared with the SB flux Thus the true flux at the surface is the same as the SB flux at the optical depth τ SB = 2/3, independently of the constant q. This is the formal derivation of where this particular reference number comes from that floats around in the literature. For other temperature profiles and for non-grey atmospheres, τ SB = 2/3 is only an estimate.
3.7 Self-consistent temperature profile and Eddington case In this paper we are considering FIB emission from a star or SN core with prescribed properties. On the other hand, in the trapping regime the FIB transfer of energy is not a perturbative effect, especially when they decouple at a radius larger than the neutrino sphere.
In this case, the atmospheric run of temperature is determined by FIB energy transport and, in a stationary state, is determined by the condition F true (τ ) = constant. Finding the corresponding B(τ ) is a formidable mathematical challenge that was solved in different ways as detailed, for example, in the book [23]. Expressing the solution in the form of Eq. (3.22), the solution q(τ ) is called the Hopf function that can be explicitly expressed, for example, as an integral that can be evaluated numerically.
We mention in passing that there is a surprisingly accurate approximation credited to Milne and Eddington that is given by the constant q = 2/3. From Eq.

True-flux convolution in geometric variables
Estimating the FIB flux with the SB approach is a good approximation that can be done for the energy-integrated flux or, if the monochromatic reduced absorption rate Γ ω strongly depends on energy, for every ω separately. However, many of the recent papers that have motivated our study used a time series of numerical SN models that were post-processed to obtain the FIB luminosity in the trapping limit. So if one anyway performs a numerical study of that type, one may as well compute directly the true flux for each energy ω based on the convolution integral Eq. (3.14a). However, while the optical depth τ as a measure of distance is very useful for conceptual discussions, it is somewhat abstract for practical implementation. More importantly, it has the disadvantage that Γ ω is assumed to decrease with increasing radius like the medium of a star so that spatial infinity corresponds to vanishing optical depth. However, for massive FIBs that can decay, for example by a → 2γ, the concept of optical depth is not directly appropriate and the FIB flux at infinity vanishes irrespective of the details of the source.
Both issues are resolved by returning to an integral over a geometric variable z which here is the coordinate perpendicular to the plane-parallel atmosphere. The convolution integral of Eq. (3.14a) for the monochromatic true flux becomes explicitly .
Here ω is the energy of a FIB with mass m a and velocity v ω = (1 − m 2 a /ω 2 ) 1/2 . The reduced absorption rate Γ ω (z) can also include free decay far away from the source. The thermal intensity B ω (z) = (π 2 /30)T (z) 4 is the one for a massless boson. Notice that one velocity factor in front of Eq. (3.14a) has cancelled against v −1 ω appearing in the Jacobian through |dτ ω /dz| = 1/λ ω = Γ ω /v ω . Moreover, λ ω (z) is the local MFP, based on the reduced absorption rate.
We have also introduced the volume energy loss rate, differential with regard to its variable ω, where the thermal FIB energy density was defined in Eq. (3.2). Recall that the spontaneous emission rate is Γ E,ω = Γ ω /(e ω/T −1), to be multiplied with the phase-space factor v ω ω 3 /(2π 2 ) to obtain the energy emission rate per energy interval dω. Together this implies Eq. (3.25) as a product of the reduced absorption rate times the blackbody FIB intensity. Notice that this applies to any process that absorbs the FIBs, including inverse bremsstrahlung or two-photon decay. The overall normalization (including a factor of 4π in B ω ) is such that is the local energy loss rate per unit volume, for example in units of erg cm −3 s −1 .
The non-appearance of a velocity factor in the flux expression of Eq. (3.24) is slightly confusing. Therefore, as a sanity check, we consider a uniform plane-parallel atmosphere at a constant temperature. The atmosphere ends at a surface at z = 0. So B ω (z) = 0 for z > 0 and is constant for z < 0. Likewise, Γ ω is constant for z < 0 and vanishes otherwise. The convolution integral can be solved analytically, however requiring many cases depending on the values of z, z and z . We find explicitly where E 3 is the third exponential integral discussed around Eq. (3.6). We show this solution in Fig. 6 where we see that F ω (z) develops a few MFPs below the surface and emerges with the Stefan-Boltzmann value v 2 ω B ω /4, including a factor v 2 ω in front of the massless-boson intensity. For a massive particle, the flux is reduced in two ways, the explicit v ω coming from the flux and one factor from the phase-space density of modes within an energy interval dω, not from the velocities of individual particles.

Including boson decay
We briefly illustrate the case where the FIBs can decay after emerging from the surface of an isothermal medium. So we assume that in the medium the (reduced) MFP is λ ω , caused by all kinds of processes, including photon coalescence. In vacuum, only free decay is possible for which we take schematically an MFP of 3λ ω . Performing the convolution, in analogy to Eq. (3.27) we find as shown in Fig. 6. The behavior in the medium depends only on the reduced interaction rate, not the individual contributions from different processes. In vacuum, where the source B ω = 0, only vacuum decay is relevant. Notice that the variation with distance is not exponential because the large-argument limit is E 3 (s) → e −s /s. The particles still decay exponentially on their trajectories, but the angle average implies that the overall flux decreases faster with distance. This behavior is a consequence of the plane-parallel model because at a large distance from a star, many stellar radii away, the flux decreases exponentially because the trajectories become more and more collinear with distance.

Rosseland average interaction rate for the diffusive flux
If the reduced MFP depends on energy, the energy-integrated true flux is given by Eq. (3.24) after performing the dω integral. The diffusive flux, on the other hand, is given by the dω integral of Eq. (3.14b). In geometric variables, one finds where B ω , given in Eq. (3.2), is the spectral blackbody density for a massless boson so that For a massless boson with energy-independent MFP, the flux expression is Therefore, if we wish to express the general diffusive flux in terms of an equivalent average MFP, comparing the two expressions yields In radiative transport, this effective MFP corresponds to the Rosseland average of the interaction rate.

Boson luminosity in spherical geometry
Our study is motivated by several recent papers concerning the FIB luminosity of a SN core and the associated energy loss. As we have argued in the previous section, the energy loss in the trapping limit can be estimated very well by quasi-thermal emission from a blackbody surface according to the Stefan-Boltzmann law. On the other hand, if one performs a numerical integration over an externally prescribed background model, one may as well use the exact expressions. Going beyond energy loss and asking for the nonlocal mode of energy transfer carried by FIBs, especially if these are radiatively unstable and can deposit energy far away from the point of production, a geometrically correct treatment is more important because a plane-parallel approximation is not appropriate if the energy is deposited far away from the compact stellar core. A similar question arises in the context of FIB energy loss and transfer in Horizontal Branch (HB) stars where the nonlocal transfer of energy was described as "ballistic" in contrast to that by diffusion [20]. Therefore, we now turn to formulating the FIB flux in spherical geometry.

Solution on a ray in geometric variables
The solution for the stationary flux in any geometry derives from the stationary solution on a given ray of the radiation field that was discussed in Sec. 2.2. Because FIBs are only absorbed or emitted, but not scattered, different momentum modes of the FIB radiation field are decoupled and so a single ray provides the mother of all solutions. Following Sec. 2.2 we thus consider a ray along some chosen FIB momentum direction, use a geometric coordinate s, and express the solutions in terms of intensities instead of occupation numbers, where ± refers to the intensities of the FIB modes at the point s along or opposite to the ray which has the direction of increasing s. The local FIB energy production rate Q ω (s) = Γ ω (s)B ω (s) = v ω B ω (s)/λ ω (s) was defined earlier in Eq. (3.25) and we assume that the blackbody intensity B ω (s) and the reduced MFP λ ω (s) are externally prescribed. The intensity at s is simply the integral over the emission from downstream of the respective direction, modified by exponential damping along the way.

Spherical volume integration: Observer perspective
As a first of two ways to calculate the boson flux at radius r in a spherically symmetric star we consider an observer at that radius and ask for the contribution of a given source to the outward or inward energy flux at that location and then integrate over all sources. Using the geometric setup shown in Fig. 7, we consider a ray with a coordinate s which is zero at r and positive in the inward direction for later convenience. The ray is tilted relative to the radial direction by an angle θ with µ = cos θ. Both Q ω (r) and λ ω (r) are assumed to be given as a function of stellar radius r. The impact parameter of this ray is b = r sin θ and half the secant line is a = r cos θ. At point s on this ray, the distance to the center of the star is given by The intensities for the two directions on this ray at s = 0 follow directly from Eq.  However, we are interested in the energy flux in the radial direction, not the intensity, and so we need another factor v ω µ, yielding at s = 0 the energy fluxes .
Notice that by our choice of direction of the s variable, it is I − ω that contributes to the radially outward flux F + ω and the other way around. Thus defined, both F ± ω,µ are positive if µ > 0 and the total flux is F ω,µ = F + ω,µ − F − ω,µ . To obtain the total flux, this expression must be integrated 1 0 dµ if we define the angle as the one between the ray and the radial direction as in Fig. 7. However, as F ω,µ consists of a piece with positive and one with negative µ, we may instead use only the first piece and integrate over all µ so that . (

4.4)
The second integral is the local intensity I ω,µ (r) times the particle velocity v ω as a function of direction. The monochromatic luminosity is L ω (r) = 4πr 2 F ω (r). As a cross check we consider the strong trapping limit, where the particle MFP is short compared with r and it makes sense to worry only about the region around r with a few MFPs upstream and downstream. In this case, the general volume integral should approach the plane-parallel result. After integrating over emission angles, Eqs. The variable s along the beam now always appears multiplied with µ and we introduce z = −µs, which is the radial distance to the position r. Notice that for positive µ a position at a radius larger than r has negative s by our convention for the direction of the considered ray. So the positive z-direction is the outward radiation direction and R s,µ → r + z. Because of strong trapping, regions a few MFPs upstream or downstream from r are suppressed by the exponential damping factor, we may nominally extend the dz integral to infinity. Therefore, the fluxes are The dµ-integrals can now be expressed in terms of exponential integrals as explained around Eq. (3.6). The total flux F ω = F + ω − F − ω can then be pieced together and reproduces a convolution integral in analogy to Eq. (3.24). The analogy becomes perfect if we reinterpret r as our z-variable and shift the integration variables accordingly.

Spherical volume integration: Source perspective
We next turn to the second picture, sketched in Fig. 8, where we consider emission from a given source and ask for its contribution to the outward and inward luminosities at some radius r. We begin with sources inside of r, all of which contribute to L + ω (r). We place the source at position a < r on the y-axis, radiating isotropically in all directions characterized by the angle β. Following a ray in the direction β with coordinate s (origin at the source), the corresponding radius vector in the x-y-plane is R < a,β,s = (s sin β, a + s cos β) so that R < a,β,s = a 2 + s 2 + 2sa cos β. . (4.10) The ray punches through the r-sphere with an angle θ and so the outward flux requires a factor cos θ. On the other hand, if we think of the ray as having a small cross section, the area of intersection with the r-sphere is increased by 1/ cos θ and so these two factors cancel. Actually if there were no damping of the emitted radiation, all bosons emitted from the source per unit time must pass the r-sphere and so the source contribution to the flux is the same, independently of the source position within the sphere. We only need to calculate the position-dependent average damping factor. Integrating over all source points within radius r, we find for the outward luminosity provided by sources inside the sphere r. Notice that in this case, contrary to the observer perspective, it is not possible to uniquely define the flux and we therefore work with the luminosity. The two quantities are however easy related F +,< ω (r) = L +,< ω (r)/4πr 2 . For a source outside the r-sphere (right panel in Fig. 8), a ray passes through the sphere if the angle β is constrained by sin β < r/a so that c min r,a = cos β max = 1 − r 2 /a 2 . (4.12) The distance from the stellar center of a point s on the ray is The length on the ray until the first and second points of intersection is r,a,β = a cos β ± r 2 − a 2 sin 2 β. (4.14) For the flux contribution of a ray intersecting a tilted surface, the same remarks apply as earlier.
The first intersection point contributes to the inward flux, the second one to the outward flux. Collecting everything, we find for the total flux where for a < r and E ω (r, a) = 0 for r = a. While this result looks far more complicated than the one found from the observer perspective, its structure is more reminiscent of the plane-parallel case in that we convolve the radial source distribution Q ω (r) with the kernel E ω (r, a) which is positive for a < r and negative for a > r and the local flux is determined by the sources a few MFPs inside and outside the considered radius.
The formal transition to the plane-parallel case is made by assuming the MFP is so small that the a-integration contributes only in a thin shell around r and we set a = r + z with |z| r. The angle integration for z < 0 covers only the range 0 < cos β < 1 and so for z < 0 we find . (4.18) Substituting z = −s cos β, this becomes The derivation for z > 0 is analogous if we notice that the contribution from the second intersection point can be dropped in the present limit. Therefore, E ω (r, a) in the small-λ limit is pieced together from exponential integral functions as in the plane-parallel case.

Luminosity at infinity
Another important limit is the FIB luminosity at infinity, corresponding to the total energy loss of the star or SN core in the form of FIBs. This quantity is only useful when the FIBs are essentially stable, otherwise the flux at infinity always vanishes. In practice, we take r to be much larger than the geometric size of the production region, but much smaller than the MFP against decay. As a consequence, there are no sources outside the very large R, so the large-R limit can be taken on the basis of the outward luminosity. Beginning with the "observer perspective," the starting point is the outgoing flux of Eq. (4.3a), implying a luminosity . (4.20) By assumption, regions far away from the star do not contribute to FIB production or decay, so the range of angles θ that contribute become infinitesimally small for r → ∞. This singularity is avoided by using instead the impact parameter b = r sin θ as integration variable. Moreover, the integration along the ray shown in Fig. 7 is performed in a shifted variable z = s − r, which amounts in the limit r → ∞ to putting the zero-point of this variable at the point of intersection of the impact line b with the ray. Using these coordinates, the radial position is R r,s,µ → R z,b = (b 2 + z 2 ) 1/2 . The angle integral thus becomes 4πr r contribute by assumption. Notice that 2πb is the circumference of a circle with radius b, so the b-integration is simply one over the stellar disk in terms of the radius (or impact parameter) on the disk. The volume integration has become one over the stellar disk and, transverse to it in the observer direction, over the new variable z. With L ω = L + ω (∞), collecting everything, and re-naming the variable of integration b → r we find The total energy loss finally requires an integration over dω. We next turn to the "source perspective" and note that for r → ∞ all sources are within the r-sphere, so as a starting point we may use Eq. (4.11) for the outward flux caused by sources within the r-sphere. In the limit r → ∞ the upper limit of integration becomes s max R,a,β → ∞. Collecting everything and re-naming the variable of integration a → r T ω (r) = e −τω,µ(r) angles .

(4.22)
The intuitive meaning is that we perform a volume integral over the radius-dependent energyloss rate, reduced by the angle-averaged transmittance T ω (r), where τ ω,µ (r) is the optical depth of the source point in a specific direction of emission. The expressions for L ω from the observer perspective Eq. (4.21) and from the source perspective Eq. (4.22) are both intuitive, yet look very different. However, one can show with a direct transformation of the integral expressions that they are indeed the same.

Transmittance in the strong-trapping limit
To calculate the FIB flux at infinity in spherical symmetry, the crucial geometric information in Eq. (4.22) is encoded in the angle-averaged transmittance where we have renamed cos β → µ. The integral in the exponential is the optical depth τ r,µ for a FIB emitted at radius r with direction µ = cos β relative to the radial direction.
In some recent papers [9,12], the transmittance was estimated as e −τ (r) , where τ (r) is the optical depth in the outward-radial direction (β = 0), i.e., the shortest way out. This approximation overestimates the transmittance except at the center of the star because otherwise the optical depth is larger in all directions compared with the radial-outward one. The prescription of Ref. [6] and re-used in Ref. [8] effectively employs an even larger transmittance. For a specific SN model, the difference between the naive transmittance e −τ (r) and that of Eq. (4.23) is shown in Fig. 10 below (dashed vs. solid blue line).
While there is no general answer concerning the difference, it is easy to estimate in the strong-trapping limit where the FIBs essentially emerge from a Stefan-Boltzmann sphere at τ (r SB ) = 2/3 and if we assume that the absorption rate decreases fast with radius at and beyond r SB . So if the geometric atmospheric height of the decoupling region is small relative to the decoupling radius we are back to the plane-parallel atmosphere approximation. In Eq. (4.23) this implies that s r for the contributing range, implying that r 2 + s 2 + 2rsµ → r + sµ. As a variable of integration we choose the vertical depth z = µs and we also note that in the strong-trapping limit the transmittance for inward-bound directions vanishes, so the dµ integral is only over positive µ. Collecting everything, we find in the plane-parallel approximation where E 2 (τ ) is the second exponential integral defined in Eq. (3.6). In the plane-parallel case, the transmittance only depends on optical depth, not geometric radial position, where here τ stands for the "outward radial" optical depth of a given source. We have also dropped the index ω for convenience. We may compare T(τ ) with the naive value e −τ in various cases. For τ = 0 we have e −τ = 1 and E 2 (τ ) = 1 and so T(0) = 1/2, which is 1/2 times the naive value of 1. The reason is that of all bosons launched at the surface, only the outward-moving ones escape. For very large τ , E 2 (τ )/ exp(−τ ) → τ −1 , so besides the previous factor 1/2 concerning the inwardbound bosons, the naive transmittance is τ times the true one and thus a vast overestimate because any trajectory that deviates only mildly from the exact radial direction implies much larger absorption. Finally, for the Stefan-Boltzmann value τ = 2/3, the ratio is 0.4968. In absolute terms, T(2/3) = 0.1239, meaning that around 1 in 8 FIBs produced at τ = 2/3 makes it to infinity. Counting only the outward-bound ones (µ > 0), almost exactly one in four escapes.
5 Explicit example I: Supernova energy loss through Primakoff production As a first example we consider axion-like particles (ALPs) with a generic two-photon interaction encoded in the coupling strength G aγγ and with a mass so small that decays are irrelevant and that we can use ultrarelativistic kinematics. In this case they are absorbed or produced only by the Primakoff process γ + Ze ↔ Ze + a on charged particles. It is reasonable to approximate the reduced absorption rate Γ ω as independent of energy [10], so this case comes close to the "grey atmosphere" approximation of radiative transfer theory. We use our expressions to calculate the total energy-loss rate L a for a prescribed numerical SN model as a function of G aγγ , compare it with the neutrino luminosity L ν , and find the two solutions for G aγγ where L a = L ν . The trapping solution is found to agree very well with the one from the Stefan-Boltzmann argument as anticipated.

Interaction model
We now consider massless ALPs that are assumed to interact with the electromagnetic field through the Lagrangian where G aγγ is a coupling constant with dimension (energy) −1 . It is the only particle-physics parameter entering our discussion. ALPs are dominantly absorbed by the Primakoff process a + Ze → Ze + γ on charged particles with a rate where n Z is the number density of targets, f S a screening factor, and f B a Bose stimulation factor for the final-state photon. The rate has been summed over final-state photon polarizations. An exact evaluation of this rate for the conditions of a SN core is not available because there are many complications as detailed in Sec. II.E of Ref. [10]. Electrons as targets are relativistic and degenerate and will be neglected. Charged nuclear targets are not only protons (as had often been assumed), but in the hottest and most important regions also small nuclear clusters. Neglecting electrons one can use Debye-Hückel screening [24], but here as well as in the target phase space we neglect degeneracy effects, probably not a bad approximation in the relevant hottest regions. As suggested in Ref. [10] we finally set where Y n is the neutron abundance (number of neutrons per baryon), keeping in mind that in general 1 − Y n is not the same as the proton abundance, although we call it the effective proton abundance. The screening factor varies only slowly in the range of energy relative to the screening scale and, given the relatively rough approximations used, we may as well set it to unity. Finally, the Bose stimulation factor is f B = (1 + f γ ) and because the targets do not recoil much, the photon energy is nearly the same as the ALP energy ω, so f γ = 1/(e ω/T − 1) and we note that f B = 1 + f γ = 1/(1 − e −ω/T ). Multiplication with (1 − e −ω/T ) to obtain the reduced absorption rate and collecting everything yields for the latter Numerically, the cross section is In this way we are naturally led to the simple case of a grey-atmosphere model which is defined by the reduced absorption rate not to depend on energy. In this case the energy integral in Eq. (4.22) can be done explicitly and the ALP luminosity at infinity is where the angle-averaged transmittance T(r) following from Eq. (4.23) is All we need to evaluate L a is a profile of (1 − Y n )ρ B and of the temperature.

Supernova model and its ALP flux
To illustrate these results we evaluate them explicitly for a numerical SN model. We use the Garching muonic SN model SFHo-18.8 at t pb = 1 s that was used in several recent studies of SN particle bounds [9][10][11]. These SN models include muons, which is a generic physical effect, although not crucial for our discussion. The models are spherically symmetric, but include convection in the form of a mixing-length treatment. (The fixed T gradient in the approximate range 8-15 km seen in the left-middle panel of Fig. 9 reflects convection.) The final neutron-star baryonic mass is 1.351 M , the final gravitational mass is 1.241 M , so the total amount of liberated gravitational binding energy is the difference which is 1.98×10 53 erg. Therefore, the released neutrino energy is at the lower end of the typical range, whereas the duration of neutrino emission is relatively short (due to convection) and the maximum temperature of around 40 MeV reached in the core is relatively small. More details are shown in Refs. [9,10], whereas the parameters relevant for us are plotted in Fig. 9.
To calculate the ALP luminosity, in principle one should include gravitational effects that are also included in numerical SN models, notably gravitational redshift as outlined in Refs. [10,11]. On the other hand, our entire treatment of radiative transfer has ignored such effects and in particular redshift and bending of trajectories. Here we are not performing a precision analysis of particle bounds but rather illustrate the relationship between volumeemission and boson-sphere Stefan-Boltzmann emission. Therefore we continue to ignore gravitational effects.
We express the ALP interaction strength in terms of the cross section Eq. (5.4) that we parameterized in terms of σ 41 = σ a /10 −41 cm 2 . The scale is chosen such that for σ 41 1 the ALP sphere will be close to the neutrino sphere at a radius of around 17 km. In the right-top panel of Fig. 9 we show the ALP production rate L a (r) for σ 41 = 1 defined in  Eq. (5.5). In normalized form and on a linear vertical scale it is the same as the red curve in the right-bottom panel. The maximum of emission is near the T maximum. In addition, the central stellar region is geometrically suppressed by the 4πr 2 factor.
In the right-middle panel, we show the transmittance of Eq. (5.6) for the indicated values of σ 41 , whereas in the right-bottom panel we show the product L a (r)T(r) in normalized form, i.e., the source distribution of the escaping ALPs. We see that the ALPs always originate from a shell of thickness of a few km. In the free-streaming limit (unit transmittance) this shell is simply given by the product of the T 4 andn profiles together with the geometric 4πr 2 factor. For larger coupling strengths, the emitting shell moves outward, driven by the transmittance that steeply falls for smaller radius, and the production rate, that steeply falls for larger r. However, the resulting shell is never very thin. The variation of widths of these normalized curves is also represented by their variation in height and we glean from the plot that the radial region of emission becomes less than a factor of 2 sharper for "surface emission" instead of free-streaming volume emission. This conclusion agrees with the schematic plane-parallel atmosphere model shown in Fig. 4.
Next we show in Fig. 10 the ALP luminosity thus derived as a function of the Primakoff cross section. We compare it with the neutrino luminosity L ν = 5.68 × 10 52 erg/s of this model. This value corresponds approximately to the neutrino-sphere region around 17 km, whereas after taking redshift effects into account it is L ν = 4.4 × 10 52 erg/s for an observer at infinity. However, as we do not include redshift effects in our ALP luminosity calculation, we compare the luminosities roughly in the local environment. On the free-streaming side, the two luminosities are equal for σ a = 1.0 × 10 −46 cm 2 , corresponding to G aγγ = 0.59 × 10 −10 GeV −1 . On the trapping side, they are equal for σ a = 1.27 × 10 −41 cm 2 , corresponding to G aγγ = 2.1 × 10 −6 GeV −1 . In which sense these G aγγ values should be seen as constraints has been discussed elsewhere [10]. Here we simply take them as the values where the ALP luminosity, calculated on an unperturbed SN model, equals L ν of that model.
In the trapping limit, we may compare the ALP flux with the one found from the Stefan-Boltzmann argument. In our model, the SB flux 4πr 2 SB (π 2 /120)T 4 SB equals L ν for r SB = 16.99 km. The cross section required to achieve τ = 2/3 at this radius is σ a = 1.03×10 −41 cm 2 , corresponding to G aγγ = 1.9 × 10 −6 GeV −1 . Therefore, within 10% one finds the same coupling strength as one found with the full transmittance-modified volume integration. The errors incurred by all other approximations, for example concerning the Primakoff cross section and concerning the impact of gravity, are of similar magnitude. Therefore, on this level of precision there is no particular benefit in performing the full volume integration that can be numerically cumbersome.
Notice that using the transmittance e −τ (r) with the optical depth only in the outwardradial direction (dashed line in Fig. 10) would lead, for the trapping regime, to the bound G aγγ = 6.1 × 10 −6 GeV −1 , a factor 3 more stringent than the correct one. This further stresses the importance of considering the correct angle-averaged transmittance as already discussed around Eq. (4.24).

Energy transfer by ALPs
Besides the SN energy loss (the luminosity seen by a distant observer) we may also ask for the ALP flux L a (r) as a function of radius in and near the SN. Its radial variation reveals the energy gain or loss by the local medium caused by ALP emission and absorption. In the source-perspective expression of Eq. (4.15) follows that the kernel E ω (r, a), in our present case, does not depend on ω and only on the radial variation of the MFP that here does not depend on temperature, so the kernel depends only onn(r) and the chosen value of σ a .
For illustration we use the trapping limit and specifically σ 41 = 1.27, where the ALP flux at infinity matches L ν = 5.68 × 10 52 erg/s. In Fig. 11 we show as a blue line the radial flux variation based on the diffusion approximation. As an orange line we show the true flux based on Eq. (4.15). The two curves separate in the decoupling region around 17 km where τ = 2/3. Beyond this region, the true ALP flux is constant. Deeper inside, it agrees with the diffusive result. We see that for radii smaller than the decoupling region, ALPs carry a significant energy flux and so would play a significant role for energy transfer within the star. 6 Explicit example II: Two-photon decay and photon coalescence As a second explicit case we consider ALPs with a mass m a so large that photon coalescence 2γ → a is the main production process, not Primakoff production which we now ignore. In a SN core, this situation pertains for m a 60 MeV [8] or in the core of horizontal-branch stars for m a 50 keV [20]. In this situation, the only information from the stellar model is the temperature profile, whereas for the ALP both the coupling strength and the mass enter.

Interaction model
Once more we consider generic ALPs with a two-photon coupling discussed in Sec. 5.1. In the ALP rest frame, the two-photon decay rate is which we use as our primary parameter to quantify the interaction strength.
The "absorption" rate caused by the decay a → 2γ for pseudoscalar FIBs was explicitly provided in the Supplementary Material of Ref. [11]. Starting from their Eqs. (S10) and (S11), the reduced absorption rate is and v ω = (1 − m 2 a /ω 2 ) 1/2 is the boson velocity. Here g B accounts for final-state Bose stimulation in the decay. Compared with the factor f B of Ref. [11], g B includes (1 − e −ω/T ) for the reduced absorption rate. In the limit T → 0 it is g B → 1 and we are back to the vacuum decay rate. The boson flux arising in Eq. (3.24) is here physically produced by photon coalescence 2γ → a, a process encoded in the reduced absorption rate of Eq. (6.2). In particular, the temperature of the background medium enters only through g B .
The local energy production rate in the form of ALPs is B ω v ω Γ ω or explicitly for example in units of erg cm −3 s −1 MeV −1 .

Diffusive energy transfer
To calculate the luminosity L ω (r) in Eq. (4.15) we need the MFP, which in our case is explicitly We show this result as a function of T /m a in Fig. 12. For T m a the effective MFP is exponentially suppressed. The interpretation is that we have defined it to describe energy transport relative to a massless boson and for large m a relative to T , the production of thermal bosons is suppressed.
To estimate the scale for the MFP required to have a significant impact on SN physics, we consider the temperature profile of our numerical SN model shown in Fig. 9. Around a radius of 10 km the temperature is around 30 MeV and the temperature gradient 4 MeV/km, then Eq. (3.31) implies a luminosity carried by ALPs of L a (λ eff /km) 66 L ν , where L ν = 5.68 × 10 52 erg/s is the neutrino luminosity of this model. In other words, unless λ eff 1 km, ALPs dominate the energy transport within the SN core. On the other hand, for a sufficiently large ALP mass, the effect is much smaller near the PNS surface where temperatures are much smaller.
We illustrate this point in Fig. 13, where we show the diffusive ALP flux for Γ −1 a = 1 km for the indicated range of masses. For small radii, where the T gradient is inward, the negative fluxes are shown as dashed lines. Taking the neutrino decoupling region to be around 17 km, we see that for m a 30 MeV, the ALP flux near the surface is smaller than L ν , whereas inside it is much larger. To avoid ALPs to dominate energy transfer within the entire SN core, and taking m a = 100 MeV, would require Γ −1 a 0.01 km and thus G aγγ 2 × 10 −6 GeV −1 .  Figure 13. Diffusive ALP energy flux carried within our numerical SN model shown in Fig. 9, based on λ eff given in Eq. (6.5) with Γ a = 1 km −1 and for the masses m a = 10, 30, 100, and 300 MeV (top to bottom). Negative fluxes (inward bound) shown as dashed lines. We also show the neutrino luminosity L ν (r) that reaches its final value between 15 and 16 km. Notice that the luminosities are in local variables, not for a distant observer, and they are in units of 5.68 × 10 52 erg/s, the local L ν near the decoupling radius.
In Ref. [8] an ALP exclusion plot is shown in the plane of m a and G aγγ , where our region of parameters is allowed. Therefore, there is a range of nominally allowed ALP parameters where they would contribute dominantly to energy transfer within the SN core, but not to energy loss. If this effect would actually make an observational difference is another question, but probably it would modify the appearance of convection in the PNS as well as the duration of the neutrino burst. In any event, we here have an explicit example of a particle that is too heavy and too short-lived to provide a SN energy-loss channel, yet has a significant effect for the energy transfer within the SN core.

Conclusions
Motivated by several recent studies about the role of feebly-interacting bosons in stars, notably in supernova cores, we have derived the equations for radiative transfer from first principles for such particles. The main simplification compared with photons is motivated by the feebleness of the interaction and leads us to neglect scattering. So we only consider boson emission and absorption by the background medium. We include systematically the effect of the boson mass that may be comparable to the local temperature or even much larger. After solving the Boltzmann collision equation for a single ray of the boson radiation field, solutions for plane-parallel and spherical geometry follow essentially from phase-space integrations, although these are not entirely trivial.
For the case of spherical geometry, the monochromatic boson luminosity at a radius r from the center of the star is expressed in the form L ω (r) = ∞ 0 dr 4πr 2 Q ω (r ) E ω (r, r ), (7.1) where Q ω (r) is the monochromatic energy-loss rate for the medium conditions at radius r and E ω (r, r ) is an integral kernel that depends on the reduced boson absorption rate as a function of r or equivalently, the corresponding MFP λ ω (r). One of our main technical results is to provide the integral kernel explicitly. The luminosity at a given radius depends on Q ω a few MFPs upstream and downstream. If this distance is short compared with the radius itself, the energy flux can be understood in the plane-parallel approximation. In this case the integral kernel simplifies considerably and corresponds to standard results in the literature. Moreover, when the MFP is small compared with the scale height of the temperature variation, one obtains the usual diffusion-limit result, where the energy flux is proportional to the MFP and the temperature gradient.
Our discussion applies to a stationary situation, when dynamical time scales are long compared with the time it takes for the FIB flux to relax to a stationary solution. In other words, we assumed the FIB flux could be calculated on the basis of a prescribed stellar model without feedback effects. Calculating the FIB flux is then a matter of integrating Eq. (7.1) over the stellar model for every energy ω and then computing the overall luminosity as an energy integral if the total luminosity is the desired quantity.
In the trapping limit, the contributing region (for every ω) is from a few optical depths to the surface and in this sense a volume integral of significant geometrical extent, not a thin shell near some hypothetical decoupling sphere. We have explicitly studied this question for the case of a "grey atmosphere," where the reduced absorption rate does not depend on energy. On the other hand, the emerging flux is surprisingly well accounted for by assuming it is emitted by a surface at optical depth τ = 2/3 with a flux given by the Stefan-Boltzmann law for a blackbody surface corresponding to the radius at τ = 2/3 with the local temperature of the background medium. The agreement is best if the temperature varies with a power law τ 1/4 as a function of optical depth. The Stefan-Boltzmann recipe has been often used and is surprisingly accurate.
When the reduced absorption rate depends on energy, possibly involving strong variations due to resonance effects, one could apply this approximation separately for every energy ω, but we have not studied how well it approximates the full integration. Of course, the Stefan-Boltzmann approximation is mostly useful as a quick estimate to avoid multidimensional numerical integrations that can become cumbersome. However, for the correct result one should simply perform the full volume integration based on the integral kernels that we have provided.
While our derivations and discussions are based entirely on standard radiative transfer theory, not all of our results can be found explicitly in the literature. In this sense we hope that our systematic exposition is useful to the astroparticle community and clarifies some issues that have emerged in the recent literature on FIB emission from stellar bodies.