New bounds on light millicharged particles from the tip of the red-giant branch

Stellar energy loss is a sensitive probe of light, weakly coupled dark sectors, including ones containing millicharged particles (MCPs). The emission of MCPs can affect stellar evolution, and therefore can alter the observed properties of stellar populations. In this work, we improve upon the accuracy of existing stellar limits on MCPs by self-consistently modelling (1) the MCP emission rate, accounting for all relevant in-medium effects and production channels, and (2) the evolution of stellar interiors (including backreactions from MCP emission) using the MESA stellar evolution code. We find MCP emission leads to significant brightening of the tip of the red-giant branch. Based on photometric observations of 15 globular clusters whose bolometric magnitudes are inferred using parallaxes from Gaia astrometry, we obtain robust bounds on the existence of MCPs with masses below 100 keV.


I. INTRODUCTION
Stellar interiors are among the best places to probe weakly coupled extensions of the Standard Model (SM), including hidden sectors containing light degrees of freedom [1].Particles within hidden sectors can be produced from extremely rare processes in the thermal plasma of the stellar interior.Despite the weak couplings, the total integrated emission rates can be relatively large due to the high density, temperature, and volume of stars.Highly interactive particles, for instance photons from the SM, have a short mean free path in stellar interiors and diffuse out slowly over tens of thousands of years.In contrast, weakly coupled hidden-sector particles have a long mean free path in the SM plasma, traversing the star within a matter of seconds.Their unimpeded emission can therefore be an important source of stellar energy transport and loss, and the emission of these particles can affect stellar evolution in analogy to that of SM neutrinos.These considerations have provided some of the strongest constraints on low-mass axions, dark photons, scalars coupling to various fermions, millicharged particles (MCPs), and other particles beyond the SM [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].
Some of the existing stellar limits on hidden sectors rely on a heavily simplified treatment of stellar interiors, in some cases assuming a single plasma frequency ω p and temperature T averaged over an entire star or even over an entire population of stars.Some hidden sector emission mechanisms are extremely sensitive to these assumed properties; for instance, even in the SM the neutrino emissivity from plasma processes Q ν scales as Q ν ∼ ω 15/2 p T 3/2 [24].In this work, we identify MCPs as a particle whose emission from stellar interiors merits a more detailed analysis.Specifically, we consider interactions of the form for Dirac fermion MCPs χ.The fractional MCP charge FIG. 1. Limits on the fractional charge q as a function of mass m.The red shaded region shows the limits (at 95% confidence) obtained in this work while the dashed red line shows previous constraints from the TRGB [5,6].Other shaded regions correspond to existing constraints from Solar energy loss (which may be modified at large values of q due to magnetic trapping) [20] and Supernova 1987A [21].Above the purple line, MCPs would be overproduced as dark matter via the freeze-in mechanism [22].Dark blue lines are the projected sensitivities of a direct deflection setup searching for a solar basin of MCPs with perturbed (solid) and unperturbed (dotdashed) phase space distributions [23].
q could arise in various ways, for instance as the lowenergy limit of a theory with an additional light U (1) ′ gauge field that mixes with SM hypercharge; for the purposes of this work, we treat q and m as phenomenological parameters without specifying their origin.It is timely to more accurately quantify the stellar limit on the existence of MCPs at and above the ∼keV mass scale, since MCPs (which would be produced at an irreducible level by the freeze-in mechanism [22]) are a key target of the sub-MeV dark matter direct detection program, with proposed experimental methodologies coming to fruition in the next decades [25].Current stellar limits on MCPs are competitive with cosmological limits on MCPs as dark matter [26].However for m near or above the ∼keV scale, the dominant MCP emission channel is a Compton-like process with a rate that is exponentially sensitive to the assumed properties of the stellar interior.It is therefore worthwhile to assess the accuracy of the simplifying assumptions that were used in previous analyses.Most notably, the backreaction from the emission of MCPs (which was not considered in previous works) can alter the properties of the stellar interior, potentially amplifying or quenching the effects of MCP emission.A fully self-consistent treatment of the effects of energy loss requires the use of modern stellar evolution codes.
We focus on the effect of MCPs on the tip of the redgiant (RG) branch (TRGB), which is one of the most sensitive probes of physics in the cores of evolved stars.Our results are summarized in Fig. 1.Energy loss to neutrinos produced by plasmon decay is known to delay the onset of helium burning in evolved RG stars.Because the TRGB brightness is insensitive to uncertainties in stellar mass, metallicity and modeling details, it is a particularly clean diagnostic of physics beyond the SM.The robustness of the TRGB is what makes it suitable for use as a standard candle in the distance ladder [27].In this work, we account for the effects of MCP emission on the TRGB using a modified version of the Modules for Experiments in Stellar Astrophysics (MESA) code [28][29][30][31][32][33].Using this framework, we simulate stellar evolution varying over MCP parameters and stellar properties like the metallicity.In contrast to previous works that relied on a simple rescaling of the results presented in Ref. [1], we compare the predicted TRGB to a set of 15 recent TRGB measurements from globular clusters (GCs), with distances calibrated astrometrically using Gaia parallaxes [34].
The rest of this article is organized as follows.In Section II we review the computation of energy loss rate from plasmon decay to light particles.We focus on the special case where the MCP is heavier than half the plasma frequency in Section II B, requiring a proper treatment for the 2 → 3 MCP emission process, e − γ → e − χ χ.In Section III we describe the physical effects of MCPs in a RG star and provide details of our implementation of energy loss into the MESA stellar evolution code.Our results and constraints are presented in Section IV.Concluding remarks and discussion follow in Section V.

II. ENERGY LOSS RATES
For a given m → n process, the energy density loss rate from a thermal plasma into final state particle X can be expressed as where F is a phase space factor expressed in terms of the phase space density f of each particle, and where plus (minus) signs correspond to Bose (Fermi) statistics.Due to the rareness of the processes that produce MCPs, their occupation number in the final state is much smaller order unity, allowing us to approximate F ≈ m i=1 f i throughout this work (this criterion can be violated in a stellar basin of non-relativistic MCPs [23], but here we are concerned primarily with energy loss where the dominant contribution is from the emission of relativistic MCPs).Note that in general, Q X is a function of the ambient temperature and density; both of these can enter in thermal phase space factors and they additionally lead to non-trivial in-medium effects for non-relativistic plasmas.

A. Emission of low-mass MCPs
For MCPs that are lighter than the plasma frequency ω p , the leading production channel is plasmon decay, γ * → χ χ, which has no vacuum analog.Plasmons correspond to poles in the in-medium photon propagator, hence their properties depend on the properties of the background.Stellar interiors can be generally characterized as either classical non-relativistic plasmas or degenerate plasmas.We determine plasmon properties using the approximations developed in Ref. [24].We find that to very good accuracy (better than the few-percent level), for the temperatures and densities relevant to stellar interiors, the residues of the poles in the photon propagator are unity in the relevant classical and degenerate limits.This applies for both transverse modes and longitudinal modes up to a maximum momentum where the plasmon would cross the lightcone, k max ∼ O(ω p ).Furthermore, we approximate the plasma frequency as where p F = (3π 2 n e ) 1/3 is the Fermi momentum.In the nondegenerate limit where p F → 0, this reduces to the classical result, ω 2 p = 4παn e /m e .With the plasma frequency of Eq. ( 4), we find that to very good approximation the longitudinal and transverse plasmon dispersion relationships are ω L = ω p and ω 2 T = k 2 + ω 2 p respectively.
Matrix elements for plasmon decay to MCPs and energy loss rates to MCPs for light MCPs with masses below the plasma frequency (such that plasmon decay is the dominant production channel).The form of the plasmon phase space depends on the dispersion relation.For the transverse mode, ω = k 2 + ω 2 p and no further closed-form analytic progress is possible (note that the integrals for the transverse processes are closely related to the plasmon number density).Meanwhile, the phase space for the longitudinal mode is independent of k, since ω = ωp, and we integrate only over k < ω 2 p − 4m 2 .
Using the analytic approximations for the plasmon dispersion relations, we can compute the energy loss rate from plasmon decay to two MCPs from Eq. ( 2), Here, p 1 and p 2 are the four momenta of final state MCPs, K = (ω, ⃗ k) is the four momentum of the plasmon with ω and k being linked through the dispersion relations, and the factor of 2 accounts for the energy lost to the second MCP since the integral is identical if we relabel 1 ↔ 2. The matrix elements for the different photon polarizations and MCP spins is in Table I, along with the total energy loss rates per unit volume.

B. Emission of high-mass MCPs
At higher masses, when ω p < 2m in a given region of the stellar interior, plasmon decay is no longer kinematically allowed.The next leading order process for producing MCPs is a 2 → 3 process resembling Compton scattering but with the final-state photon replaced with an off-shell plasmon that branches into 2 MCPs, e − γ → e − χ χ.It has been previously shown that the energy loss through this process can be modelled in terms of the decays of a thermally distributed population of "offshell" plasmons [6].In fact, Eq. ( 5) can be modified to account for a population of plasmons with an arbitrary four momentum K as [6,20] where ρ L,T (ω, k) is the spectral function for the plasmon which can be understood as the probability that a plas-mon with the four momentum K exists in the plasma [35], Here, Π L,T is plasmon self-energy at finite temperature and density.The real part of Π L,T sets the dispersion relation of the associated mode, given by the approximations of the previous Subsection.On the other hand, the imaginary part of the self-energy specifies the damping rate.At the energies we consider, the transverse mode is primarily damped through Thomson scattering, In contrast, the longitudinal mode is either Landau damped (when k > k max ), or damped by Thomson scattering (when k < k max ) [24,36].Since we only integrate over momenta k < k max , the damping rate for longitudinal plasmons is defined by the corresponding Thomson scattering cross-section [36] Im Π L = 8πα 2 n e 9m e T ωk ω p .
Using these expressions for the real and imaginary part of the plasmon self-energy and the matrix element for plasmon decay, we can write a general energy loss rate for transverse and longitudinal plasmons into MCPs with an arbitrary mass.Integrating Eq. ( 6) over the MCP phase space, we obtain, where the momentum-averaged matrix element is TABLE II.General expression for the energy loss rate to MCPs using the matrix elements provided in Tab.I and the plasmon spectral density ρL,T (ω) as defined in the text.These expressions reduce to the ones in Tab.I in the limit ωp ≥ 2m.The spectral function, ρ L,T (ω, k) is highly peaked when the dispersion relation is satisfied, i.e., when the plasmons are on-shell, ω 2 − k 2 = ReΠ L,T , where Re Π T = ω 2 p and Re Π L = ω 2 p − k 2 .This implies that for ω p ≥ 2m, the integral over the spectral function is sensitive to this peak and therefore dωρ L,T (ω, k) is well approximated by a Dirac delta function.In this case, Eq. ( 10) reduces to the expressions given in Table I for both the longitudinal and transverse case.On the other hand, for ω p < 2m, the integral over ω always misses the peak of the spectral function and the energy loss rate is exponentially suppressed.Additionally, it becomes highly sensitive to the MCP mass and the temperature of the stellar interior.This has the potential to significantly impact the shape of the constraint depending on the assumptions about the stellar interior.We provide the complete expressions for Q MCP for the longitudinal and transverse case in Table II.
In Fig. 2, we show the energy loss rate through the decay longitudinal and transverse plasmon modes as a function of the MCP mass, for T = 100 eV (left) and T = 1 keV (right) and for ω p = 0.042 keV and ω p = 0.42 keV.The energy loss rate is clearly quite sensitive to the temperature, especially for off-shell plasmon decay to high-mass MCPs.Furthermore, for ω p ≥ 2m, the energy loss rate is independent of the MCP mass, while for ω p ≤ 2m, it is exponentially suppressed as discussed above.We note that although the transverse mode gives the dominant contribution to the energy loss in a wide range of parameter space, the contribution of the decaying longitudinal mode cannot be ignored and may in fact even exceed the transverse contribution for large masses and low temperatures (left panel of Fig. 2).

III. ENERGY LOSS IN RG STARS A. Effects of Beyond-SM Particle Emission
In standard stellar evolution, stars that have exhausted their core supply of hydrogen depart from the Main Sequence.As the inert helium core contracts under gravitational pressure, its internal temperature rises, driving the star upward in luminosity and causing the envelope to expand, cooling the surface of the star, and creating the RG branch in colour-luminosity space.At the same time, a shell surrounding the helium core continues to fuse hydrogen, enlarging the core over time.In low-mass stars, this trend terminates at the TRGB with the ignition of helium, which is almost solely dependent on the mass of the helium core and is largely independent of other stellar properties.
The required core mass for the onset of helium ignition is rather sensitive to energy loss.For instance, SM plasmon decays to neutrinos have the effect of enlarging the required core mass for helium ignition.Beyond-SM neutrino properties, such as a small charge or magnetic moment, can be constrained because this would yield an even brighter TRGB than observed [12][13][14][15][16].Many previous works use a simple criterion in setting limits on physics beyond the SM, demanding that the total energy loss to new states not exceed twice the standard neutrino luminosity of the core (as this would cause the core mass at helium ignition to change by more than 5% with respect to the standard theoretical expectation [1]).However, later work has included self-consistent stellar simulations with new physics.For instance, Ref. [16] simulated a population of stars, constraining the neutrino dipole moment based on the TRGB of the ω-Cen GC.More recently, Refs.[37,38] used MESA to simulate the evolution of GCs in the presence of axion-like particles and dark photons, respectively.They found that the more realistic treatment of energy loss can lead to bounds that are stronger than previously-published results for some masses, while opening up other areas of parameter space.There have also been recent developments on the observational side, for example Ref. [39] compared a set of TRGB calibrations from other galaxies to local GCs, updating I-band magnitudes using Gaia distances and combining observational and theoretical uncertainties in order to set limits on axion and neutrino properties.
Prior constraints on light MCPs from stellar energy loss have assumed an energy-loss mechanism similar to plasmon decay to neutrinos [2], i.e. leading to core cooling and a delayed (and therefore brighter) TRGB.Ref. [6] updated the computed emission rates by including offshell plasmon decay.However, constraints were set by a simple rescaling of the limits in Ref. [1], assuming a single plasma frequency across an entire star.Ref. [20] computed the effect of MCP emission in the Sun using the Garching Stellar Evolution Code (GARSTEC) [40] and set relatively robust bounds using helioseismology.In the calculations described below, we aim to apply many of these state-of-the-art developments in the field in order to modernize the bounds on MCPs from the TRGB.

B. Numerical modeling
We use MESA to account for the cumulative effect of MCP emission over the lifetime of the star.We use the run star extras module to include the extra energy loss rate from both the transverse and longitudinal modes as given in Table I at each time step, based on the stellar temperature and composition profiles.This is then fed into the stellar evolution code via the extra heat hook,  as a radial array of energies that contribute to the total energy budget of the star.We run our simulations starting with a pre-main-sequence model, and evolve them through to the horizontal branch phase.We take the maximum brightness as the L TRGB .The absolute magnitude M bol is then obtained via the standard definition where M and L are absolute magnitudes and luminosities, respectively.
The GC sample that we use contains clusters with ages ∼ 10-13 Gyr.In this range, the TRGB is represented by stars with masses in the range 0.8-0.9M ⊙ .Rather than simulate a full population for each cluster, we perform simulations using a set of fiducial values for the astrophysical and modelling parameters.We additionally vary these parameters to quantify the extent to which they affect our results.As a benchmark, we focus on a 0.8 M ⊙ star and set the initial helium fraction to Y = 0.241, the mixing length parameter α = 2.0 and we do not include mass loss.Throughout this work, all simulations shown assumed this same set of fiducial values unless otherwise specified.Because GCs cover a wide range of metallicities, we do vary [M/H] to cover the full range spanned by the data, and interpolate results to match the reported value of [M/H] for each cluster.
To verify the robustness of these assumptions, we have performed a set of simulations with and without the effects of plasmon decay to MCPs that span a range in mass-loss rate, mixing length parameter α, metallicity and initial helium fraction.The results are presented in Fig. 3 in terms of the variation of M TRGB with respect to the fiducial magnitude (when all the parameters are set to their default values).In Fig. 3, the fiducial value for [M/H] is set to −1.42, the value reported for ω-Cen [34].We also show the 68% containment on ω-Cen's TRGB luminosity as a grey band.
We vary the mass-loss rate by scaling the Reimers mass-loss factor η from zero to 1.6 × 10 −13 M ⊙ yr −1 , which is twice the typically-used value [41] of 8 × 10 −14 M ⊙ yr −1 , found to reproduce the observed envelope mass of RG stars.As shown in Fig. 3, even unphysically large mass-loss rates do not appreciably change the the TRGB magnitude.
Similarly, we vary the convective mixing length α and initial helium fraction over a wide range, finding no significant change in M TRGB , except for especially large values of Y ≳ 0.28 where dimming of the TRGB is observed.Because of the age of GC stars ∼ 13 Gyr, the initial helium content should be close to primordial Y P ≃ 0.24, so such a large Y would have to have resulted from an unrealistically large core dredge-up after H exhaustionthis is unlikely, as the effects of He dredge-up should be more-than-compensated by diffusion (see e.g.Ref. [42]).
Fig. 3 also shows the effect of varying the metallicity within the reported errors for ω-Cen.The small correlation that can be seen is consistent with the expected effect of metallicity on the TRGB, but is small enough within this range to justify using a fixed value of [M/H] for each individual cluster.
Finally, Fig. 3 shows the change in luminosity of the TRGB with stellar mass, becoming slightly dimmer for SM stars with increased mass; once energy loss from plasmon decay to MCPs is added, the TRGB actually becomes insensitive to the stellar mass.

IV. EFFECT OF MCPS ON THE TRGB
A. Simulated stellar interiors Fig. 4 shows the energy loss profile in a simulated 0.8 M ⊙ star due to SM processes and due to the effect of plasmon decay to MCPs, just before the TRGB.The key SM processes include neutrino production from nuclear fusion, plasmon decay, bremsstrahlung, photo-neutrino, pair annihilation, and recombination.The inclusion of MCP emission, which alters the free-electron density and temperature profile of the star shown in Fig. 5, also affects the SM energy-loss rates.For all the energy-loss channels, solid lines represent the SM while dashed lines represent the case where the emission of light MCPs with q 5 × 10 −14 is included self-consistently, lowering the energy loss from SM channels.The structure of the emis- sion region for plasmon decay to MCPs, which has prominent support in the envelope outside of the degenerate helium core, also differs substantially from SM plasmon decay to neutrinos which primarily occurs in the core.Thus, the previously used method of modelling MCP production as a scaling correction to the neutrino production rate leads to qualitatively incorrect modelling of the stellar structure.Fig. 5 shows the stellar temperature T , free-electron density n e , and MCP emission Q MCP profiles as a function of the enclosed mass for evolved stars on the RG branch, just before the TRGB.The lower panel shows the stellar radius as a function of enclosed mass.The stellar profiles are shown for different MCP masses at a fixed charge (top), and for different MCP charges at fixed mass (middle).The solid line represents the stellar structure in the absence of new physics.In all cases, the energy loss from MCPs in the core leads to a reduction in core temperature, a longer period of hydrogen shell burning, and the production of a more massive core before the onset triple-alpha ignition.
As expected, the effect on the core mass is larger for smaller m, where MCP emission in the entire core can occur.The high plasma frequency required for on-shell decay to more massive MCPs, discussed in Section II B, leads to a suppression of MCP emission, resulting in an additional step down of Q MCP in the inner core.The middle panel of Fig. 5 shows the effects of a larger value of q on a star for m = 7 keV.Here, q is increased to 7 × 10 −13 , around an order of magnitude larger than the MCP charge is fixed at q = 7 × 10 −13 assuming m = 7 keV (dashed), m = 10 eV (dotted).Middle: the MCP mass is fixed at m = 7 keV, with q = 7 × 10 −13 (dot-dashed), q = 5 × 10 −14 (dashed).In both panels, the SM prediction is shown as a solid line.Bottom: the stellar radius as a function of the enclosed stellar mass.FIG. 6. Stellar isochrones for a cluster of stars with stellar masses between 0.7 M⊙ and 3 M⊙ and a metallicity representative of ω-Cen, [M/H] = −1.42,evolved up to the TRGB (indicated with a star) without MCPs (solid) and with MCPs with m = 10 eV and q = 2 × 10 −14 (dashed).The isochrones were computed using the MESA Isochrones and Stellar Tracks package [43].
constraints that we obtain at this MCP mass.In this case, the core temperature drops by a factor of three, actually lowering MCP emission from the inner core in comparison to models with a lower charge.The subsequent core contraction ultimately changes the qualitative features of the stellar structure.A shallower density gradient develops, erasing the sharp boundary between the helium core and the envelope, producing a broad, hightemperature shell leading to enhanced MCP energy loss.The main features remain, notably a larger helium region, and a TRGB that is a full 0.8 Mag brighter than in the SM-more than three standard deviations away from the measured value for M92.
The net effects of MCP emission lead to a delay in helium ignition and thus a higher luminosity (lower bolometric magnitude) as a star reaches the TRGB.We illustrate the effect on a population of stars with isochrones representative of ω-Cen shown in Fig 6 .These colourluminosity diagrams show the distributions for a population of stars with [M/H] = −1.42 and masses between 0.7−3 M ⊙ , taken at age snapshots of 1, 5, 10 and 13 Gyr in the SM (solid lines) and in the presence of a q = 2 × 10 −14 MCP in the low-mass m ≪ 1 keV limit (dashed lines).For representative ages of Milky Way GCs (∼ 10 − 13 Gyr), the TRGB luminosity increase due to MCPs remains insensitive to the exact age.

B. Constraints
In order to set limits on the MCP charge q as a function of the MCP mass m, we compare the results given by our simulation to recently-reported TRGB magnitude measurements.Ref. [34] obtained bolometric mag- FIG. 7. Bolometric magnitude at the TRGB as a function of MCP charge q, for MCP masses at ranging from 100 eV to 100 keV, from left to right.Fiducial values, marked as red crosses in Fig. 3, are assumed in all simulations, with [M/H] = −2.05, the metallicity of M92.Also shown are the observed magnitude (dashed) and 2 σ observational uncertainties of the single cluster M92 on the TRGB luminosity from Ref. [34] (blue), and total error combined with the theory uncertainty (grey) recommended in Ref. [46].
nitudes for 22 GCs using photometric data from the Hubble Space Telescope and ground-based optical measurements.They report two distinct values of M bol based on different distance-measurement techniques: first via a zero-age horizontal branch (ZAHB) calibration, and second with distances measured using Gaia astrometric data [44] and line-of-sight velocities from groundbased instruments including Keck the Very Large Telescope [45].Even though parallax distances are only provided for 16 of the 22 clusters, they are reported with higher precision and with errors that are more reliably estimated for each cluster.Though we have not explicitly checked in the present work, it is also possible that the properties of the ZAHB luminosity will be modified by the additional MCP emission considered here, whereas astrometry would be unaffected.Therefore, while we compare constraints on MCPs obtained using both techniques (parallax and ZAHB) in Appendix A, for our main results we use the parallax distances.The astrometric data set does contain one outlier, NGC 6553, which has a measured magnitude of M bol = −4.76± 0.17.This is in significant tension with the reported ZAHB measurement of −3.93 ± 0.25 and is 0.7 dex away from any of the other reported magnitudes for other GCs.Here, the emission of MCPs would actually lead to a significantly better fit between theory and observation; however, producing a magnitude this low requires a charge that is well beyond values allowed by the combination of other clusters.In some cases, the changes to the stellar interior induced by such a high MCP charge (and therefore a high MCP emission rate) cause our models to fail to converge.these reasons, we omit NGC 6553 from our analysis.
We simulate stars on a grid of MCP masses, MCP charges, and stellar metallicities covering the range of clusters studied in Ref. [34].Fig 7 shows the TRGB magnitude predicted in a 0.8 M ⊙ star with the metallicity of M92 for a few MCP masses as a function of the fractional charge.This is overlaid with the 68% containment on the TRGB magnitude of M92 [34].Once the energy-loss rate becomes large enough, we observe a sharp increase in TRGB luminosity.The impact of increasing q is mitigated at higher MCP masses, as off-shell plasmon decay becomes thermally suppressed, requiring much larger values of q for any effect to be observed.
To evaluate the upper bound on MCP charge at a given MCP mass, we define a Gaussian likelihood where the index i runs over the list of GCs, M i (q, m) is the predicted bolometric magnitude of the TRGB as a function of MCP charge and mass at the metallicity of a given cluster, M i,obs is the bolometric TRGB magnitude inferred from measurement, and σ 2 i represents the uncertainties on the inferred bolometric magnitudes.Our uncertainties include the reported observational uncertainties given in Table 6 of Ref. [34], added in quadrature to the theoretical uncertainties on the bolometric magnitude predicted by stellar evolution codes.Following the recommendation of Ref. [46] we set this theoretical error to 0.12.Applying Wilks' theorem, we use a likelihood ratio test to find the limit on the MCP charge q lim for for a fixed MCP mass m.We set the limit by finding 2 log L({M i,obs }|q max , m) − log L({M i,obs }|q lim , m) = 2.71, where q max is the value of the MCP charge that maximizes the likelihood at MCP mass m and 2.71 corresponds to the one-sided 95% confidence level (CL) [47].We simulated MCP masses from 10 eV to 100 keV, and show the resulting limits on the charge q in Fig 1 .We additionally show prior limits on MCPs from Solar modeling [20], the absence of anomalous cooling of Supernova 1987A [21], and projected constraints from proposed direct deflection searches [23].As expected, our constraints plateau for low masses, as m ≪ ω p .For the lowest mass considered here, the resulting 95% CL limit is a little more than a factor of three stronger than limits based on those of Ref. [2].As m increases, the stellar region (notably, the He core) in which on-shell plasmon decay can occur shrinks, leading to a rapid loss of sensitivity above masses corresponding to half the lowest plasma frequency in the core.The "kink" seen around 10 keV corresponds to the transition between having kinematically accessible on-shell plasmon decay in some part of the star versus relying purely on off-shell decays throughout the entire star.At high mass, sensitivity is lost at an exponential rate because of the temperature scaling of the off-shell plasmon decay rate.Coincidentally, our limit at high masses is very similar to the estimate performed in Ref. [6], which was obtained by rescaling the low-mass results from Ref. [2] using a single plasma frequency for a single idealized TRGB star obtained from assuming SM-only properties of the stellar interior.The parts of our improved analysis that would enhance the strength of the constraint (including energy loss from the longitudinal mode, using modern simulations and observations, etc.) are compensated by the effect of MCP emission on the stellar interior, which reduces the energy loss from SM channels.

V. CONCLUSIONS
We have presented a self-consistent analysis of the effects of MCP emission (via on-shell and off-shell plasmon decay) on the TRGB by implementing this model into a modified version of MESA.As discussed in Section II, we have implemented the total energy loss through both the transverse and longitudinal plasmon decay channels, using the in-medium matrix element for the corresponding decay process, as well as the full spectral function for the longitudinal and transverse plasmons.The decay of the longitudinal mode, which we included for the first time, is the dominant contribution to energy loss for higher MCP masses and lower stellar temperatures.Using these rates, we simulate stars starting with a premain-sequence model and evolve until the TRGB.We find that excess energy loss from MCP emission can alter the density and temperature profiles of the stellar interior, lowering the energy loss rates from SM processes (e.g.plasmon decay to neutrinos), which partly compensates for the effect of MCP emission.We also find that the emission morphology of MCPs is quite distinct from SM emission channels, occurring outside of the degenerate helium core.The net effect is that if MCPs exist, the TRGB should look substantially brighter due to their emission.We have checked that the size of this effect is robust to varying assumptions about stellar ages and compositions, mass-loss, and convective mixing.To set a limit on the existence of MCPs from the TRGB, we use 15 GCs whose parallax distances are measured with Gaia.Our resulting limit is stronger and more robust against theoretical uncertainties compared to previ-ous analyses that made simplifying assumptions that are not supported by our simulations.
The analysis presented in this work paves the way for several avenues that may further improve the limits on MCPs.When setting limits on MCPs from the horizontal branch, previous analyses did not self-consistently account for deviations in the stellar density and temperature profiles compared to the SM-only stellar structure.Given that we see such deviations in our simulations, it will be worthwhile to evolve forward in time to see the progression and to determine whether the true MCP emission can be enhanced or quenched compared to the MCP emission assuming a SM-only stellar structure.Futhermore, we have not accounted for any trapping effects due to magnetic fields; such effects were previously estimated to be relevant in the Sun for large q, potentially leading to non-trivial energy transport [10].It will therefore be worthwhile to consider magnetic trapping for post-main-sequence stars, though little remains known when it comes to magnetic fields in stellar interiors.Finally, there may be an opportunity to set stellar bounds on MCPs from distinct stellar observables with different systematic uncertainties, including asteroseismology.Other stellar populations could provide additional information.For instance, younger clusters may provide more stringent information, as the "plateau" in the TRGB luminosity versus stellar mass does not dip as quickly as in the SM when MCPs are included.We leave consideration of all these directions to future work.

FIG. 3 .
FIG.3.Variation of the bolometric magnitude M TRGB of ω-Cen without (triangles) and with (dots) MCPs (m = 10 eV and q = 2 × 10 −14 ), compared to our fiducial model, when varying the convection mixing length parameter α, metallicity [M/H], initial helium fraction Y , and stellar mass.In each set of simulations, only one parameter varies at a time while others are held at their fiducial values (red × on the corresponding axes).The [M/H] axis covers the range of uncertainty in the metallicity measurement of ω-Cen.The 1σ uncertainty on M TRGB is indicated by the shaded band.