Superradiance in rotating stars and pulsar-timing constraints on dark photons

In the presence of massive bosonic degrees of freedom, rotational superradiance can trigger an instability that spins down black holes. This leads to peculiar gravitational-wave signatures and distribution in the spin-mass plane, which in turn can impose stringent constraints on ultralight fields. Here, we demonstrate that there is an analogous spindown effect for conducting stars. We show that rotating stars amplify low frequency electromagnetic waves, and that this effect is largest when the time scale for conduction within the star is of the order of a light crossing time. This has interesting consequences for dark photons, as massive dark photons would cause stars to spin down due to superradiant instabilities. The time scale of the spindown depends on the mass of the dark photon, and on the rotation rate, compactness, and conductivity of the star. Existing measurements of the spindown rate of pulsars place direct constraints on models of dark sectors. Our analysis suggests that dark photons of mass $m_V \sim 10^{-12}$ eV are excluded by pulsar-timing observations. These constraints also exclude superradiant instabilities triggered by dark photons as an explanation for the spin limit of observed pulsars.


I. INTRODUCTION
The nature of dark matter is one of the biggest open questions in physics. Broadly speaking, there are two approaches to explain the gravitational anomalies that indicate the existence of dark matter. The first is to change the way that gravity works on large scales while preserving the short-distance behavior, e.g. MOND [1]. However, theories of modified gravity still require the addition of a dark matter particle to explain large scale structure [2,3]. The second approach postulates that the gravitational anomalies are due to dark matter. The most popular of these explanations advocates the existence of new degrees of freedom beyond the Standard Model (SM) that form a dark sector.

A. Ultralight bosonic fields
Some popular candidates for dark sector matter are ultralight bosonic fields. Indeed, bosonic fields are a generic feature of many theories [4,5]. A well-motivated scalar candidate is the QCD axion, a light bosonic degree of freedom introduced in physics to explain the smallness of the neutron electric dipole moment, years before the dark matter problem was fully appreciated [6][7][8]. In addition, a plethora of new light scalars was predicted to arise in the String Axiverse [4], making them important potential dark matter candidates.
Vector candidates are equally well-motivated. Additional U (1) gauge sectors arise in many string-motivated extensions to the SM [5,9,10]. In these scenarios, there can be extra degrees of freedom which are charged under both the U (1) hypercharge of the SM and a "hidden" U (1) , known as dark photons. This has motivated the study of kinetic mixing of the hidden sector with the SM.
Many of these searches have been focused on eV-GeV scales using direct detection and low-energy accelerator experiments (see e.g. [11] for a summary of current efforts).

B. Superradiant instabilities and ultralight fields
However, ultralight (i.e., sub-eV) fields which are weakly coupled to SM particles are difficult to probe with traditional colliding beam, fixed-target, and direct detection experiments. Instead, one can search for their imprints through their gravitational effects. A promising mechanism to probe bosonic fields is rotational superradiance [12][13][14][15]. Superradiance affects all known free, bosonic fields and has been well-studied for black holes. In this context, low-frequency wavepackets of bosonic fields are amplified upon scattering off rotating black holes, when the frequency of the field wave satisfies ω < mΩ, where m is the azimuthal number and Ω is the angular velocity at the event horizon [15]. When the bosonic field is massive, the effects of superradiance turn the entire system unstable [16][17][18][19][20][21][22][23], and the instability gives rise to a slowly-spinning black hole surrounded by a cloud of bosonic field [cf. Ref. [15] for an overview]. This cloud has a time-dependent quadrupole moment, and slowly dissipates through gravitational waves producing a monochromatic signal, which is a promising channel and smoking gun for new physics [24][25][26][27]. Furthermore, because superradiance drives the spin down, observations (either in the electromagnetic or gravitational-wave spectrum) of the spin-mass diagram of black holes may also bring convincing evidence for new physics [28]. Finally, it is also possible that superradiant effects are directly observable through enhanced scattering of electromagnetic or gravitational waves [29,30], or even through in-stabilities triggered in interstellar plasma environments surrounding black holes [31,32].

C. Superradiance in stars
However, superradiant effects are not limited to rotating black holes and in fact can appear in any classical system that is able to absorb radiation [12,14,15,22,33,34]. In this work, we show that superradiance also occurs in the presence of rotating and conducting spheres, and most notably in (rotating) stars with nonzero conductivity. This seemingly classical problem in electromagnetism has never -to the best of our knowledge -been worked out. We find that rotating stars amplify lowfrequency photons, whenever their frequency satisfies the usual superradiant condition, ω < mΩ, where now Ω is the rotational velocity of the fluid.
These superradiant effects may have interesting implications for theories of dark photons as well as more complicated hidden sector theories. We find that massive dark photons trigger an instability of rotating and conducting stars, analogous to the black hole case. Further-more, the superradiant effects may be entirely contained within the dark sector, but have observable consequences that are worthy of further investigation. The most direct signature of this scenario is the spindown of pulsars due to the superradiant instability. As we discuss, existing pulsar-timing measurements of the spindown rate of pulsars already constrain these models. Because these pulsar-timing constraints are rather stringent, they also exclude superradiant instabilities triggered by dark photons as an alternative explanation for the spin limit of observed pulsars (cf., e.g., Refs. [35][36][37][38][39] for a discussion on proposed limiting mechanisms on the spin of pulsars). Throughout this work, we use G = c = 1 units and unrationalized Gaussian units for the charge.

A. Maxwell and Proca theory in curved spacetime
To understand the effects of superradiance in stars, we work with the theory involving one vector field A µ with mass m V = µ V , where F µν ≡ ∇ µ A ν − ∇ ν A µ is the field strength. The vector A µ can describe either Maxwell theory with the standard massless photon, in which case µ V = 0, or a Proca theory in which the vector field is massive. We will show below that in both cases there are nontrivial superradiant effects around rotating stars. The theory above is a toy model designed to capture the main features of a general relativistic theory where a (possibly hidden) vector field is minimally coupled to the geometry.
The resulting field equations are G µν = 8πT µν matter + 16π where G µν is the Einstein tensor and T µν matter is the standard stress-energy tensor of matter fields, the latter being collectively described by S m in action (1).

B. Background: slowly-rotating, conducting star
Because the star is assumed to be uncharged, fluctuations in the vector A µ affect the geometry only at the quadratic order. Thus, to linear order, we can consider a standard general relativity background as a fixed geometry, around which the vector field evolves. We will always neglect backreaction of the vector field on the geometry. This is a reasonable approximation for all known astrophysical setups.
We consider a slowly-spinning star and neglect quadratic-or higher-order corrections in the spin. To linear order in the spin, the background line element is described by and the star's four velocity reads where Ω is the rotational velocity of the fluid. The slowrotation approximation requires Ω Ω K , with being the mass-shedding frequency, whereas M and R are the star's mass and radius, respectively. In the exterior, F = B = 1 − 2M/r and ζ = 2J/r 3 , where J is the angular momentum of the star. The interior depends on the type of matter and it is described by the classical Tolman-Oppenheimer-Volkoff equations for a perfect fluid with T µν matter = (P + ρ)u µ u ν + P g µν , namely where we defined F = e 2Φ , B = 1 − 2M(r)/r and := Ω − ζ(r). Assuming a barotropic equation of state in the form P = P (ρ), these equations can be integrated numerically with standard methods. For simplicity, we will focus on backgrounds describing a constant density, perfect-fluid star. In this case, the static part of the metric has an exact solution, where ρ = 3M/(4πR 3 ). The equation for (and there-fore for ζ) cannot be solved analytically for generic values of the compactness, whereas in the Newtonian limit yields ζ(r) = 2J/R 3 = const in the interior, which smoothly connects to ζ(r) = 2J/r 3 in the exterior. Finally, the vector A µ is evolving in the vicinities of an uncharged, rotating star made of material with conductivity σ and proper charge density ρ EM . We assume that the coupling between the vector and the material is given by the constitutive Ohm's law, which in covariant form reads [14], where all quantities are computed in the frame of the material whose 4-velocity is u α . This relation should be accurate for weak fields and represents the lowest order term in the family of possible couplings between the material and the vector field.
C. Perturbations of a spinning, conducting star in Maxwell and Proca theory An uncharged star in electrovacuum (A µ = 0) is a trivial solution to the previous equations. We now wish to understand linearized fluctuations around this background. We start by expanding the vector field A µ in 4-dimensional vector spherical harmonics, The first term on the right-hand side has parity (−1) l+1 and the second term has parity (−1) l , m is an azimuthal number and l is the angular number. Likewise, we expand the charge density in scalar spherical harmonics, ρ EM (t, r, ϑ, φ) = e −iωtρ EM (r)Y lm . Because the background is not spherically symmetric, the above decomposition introduces couplings between polar and axial modes and between perturbations with different harmonic indices [40]. To linear order in the spin, the coupling between polar and axial modes can be consistently neglected and one is left with an "axial-led" and a "polar-led" system of ordinary differential equations [19,20,41]. The decoupling procedure is given in Appendix B. Here, we report only the final result for the axial-led system to linear order in the spin, where dr/dr * = √ BF . Note that, within our slow-rotation approximation, σ can be a generic radial function. For simplicity, we take σ = const.
The polar sector is more involved and we leave it for future work. Here, we briefly mention that in the massless case (µ V = 0) the polar sector can be reduced to a single second-order differential equation by using some gauge freedom, whereas in the Proca case the polar sector describes the propagation of two physical degrees of freedom, and one is left with a system of two, coupled, second-order, differential equations. In both cases, the charge densityρ EM (r) is fixed in terms of σ and of the perturbations of the electromagnetic field by the field equations, similarly to the fluid density ρ which is fixed by the Tolman-Oppenheimer-Volkoff equations in terms of the pressure once an equation of state is given. This can be also understood by the fact that an applied electric field will modify the charge distribution, even when the object is globally neutral.

III. SUPERRADIANT SCATTERING FROM SPINNING AND CONDUCTING STARS
We now consider a scattering experiment. We focus on the axial sector, but the computation for the polar sector, although more technically involved, follows similarly. In the axial sector, the solutions to Eq. (13) behave asymptotically as We have selected the regular solution at the center of the star. From our conventions for the time-dependence of the fields, it follows that this state is composed of a piece, A in e −iωr , which is an ingoing wave and is scattered by the star, giving rise to an outgoing component A out e +iωr . It is also easy to verify that the incoming and outgoing fluxes at infinity are proportional to |A in | 2 and |A out | 2 , respectively [13]. We thus define the superradiant factor We have computed the superradiant factor Z numerically, by integrating Eq. (13) from the center of the star, outwards to some finite but large value of the radial coordinate r, where the numerical solution is matched against a higher-order version of expansion (16). The numerical results are shown in Fig. 1. As expected, Z > 0 when the superradiant condition is satisfied, ω < mΩ. The amplification factor grows with σ, until it saturates in the large-σ limit displaying a sharp maximum at ω mΩ. Although not shown, the amplification grows with the compactness and with the spin of the object.
We can also gain some analytical insight on the superradiant amplification. In the Newtonian limit, the external solutions are linear combinations of Bessel functions √ rJ l+1/2 (ωr) and √ rY l+1/2 (ωr). In the interior, and for small conductivities, the only regular solution admissible is Matching the functions and their derivatives at the surface of the star and expanding for small frequencies, we find (18) The above expression agrees remarkably well with the exact numerical result up to M/R ∼ 0.2 and for σM 1. This relation is also interesting, as it extends an observation made in Ref. [34]: one can try to naively compute the superradiant amplification factors of Kerr black holes by letting R = 2M and 1/σ = M , as this is now the only possible time scale in the problem. With this substitution, the above relation predicts that slowly rotating black holes in general relativity amplify l = 1 scalar fields On the other hand, a matched-asymptotic expansion calculation in full general relativity yields the same result to within an order of magnitude (the coefficient turns out to be 8/9 instead of 64π/45) [15]. As we show in Appendix A, one can improve on this relation by using the membrane paradigm for describing horizons [42]. In this framework, horizons are endowed with a surface conductivity of 1/4π, and a simple Newtonian analogue recovers exactly the general relativistic prediction.
For large conductivities, we have been unable to find concise analytical expressions, but in the Newtonian limit our results are well approximated by in the superradiant regime, with k 1 ∼ 0.78, k 2 ∼ 0.09 and c 1 ∼ 2, c 2 ∼ 25. The amplification factor is peaked at ω − mΩ ∼ 1/(σR), and bounded. The analytical expression above is not very accurate close to the peak of the amplification factor, but we find numerically that, for σM 1, the l = m = 1 peak is described by where, interestingly, the prefactor decreases at large compactness.

IV. SUPERRADIANT INSTABILITIES OF SPINNING AND CONDUCTING STARS
In analogy with the black hole case, we expect that the mass term for the Proca field can lead to superradiant instabilities in conducting stars. We show this explicitly by solving the perturbation equations numerically as an eigenvalue problem, and computing the quasinormal modes of the system, ω lmn = ω R + iω I , where n is the overtone number. In our notation, an instability corresponds to ω I > 0, and τ ≡ 1/ω I is the instability time scale. The parameter space of the spectrum is large and complicated, since -even for fixed "quantum" numbers (l, m, n) -it still depends on four dimensionless parameters, namely (µ V M, M/R, Ω/Ω K , σM ).
In the axial case, our results for the l = m = 1 fundamental unstable mode are well approximated in the small µ V M limit and to linear order in Ω/Ω K by (22) where α i are dimensionless constants that depend on the compactness and also on Ω since the combination Ω/σ is not necessarily small. Besides the prefactor in square brackets in Eq. (22), the functional form of the superradiantly unstable modes is the same as that found for a black hole [19][20][21]43]. The dependence of the prefactor in Eq. (22) on σ and M/R are presented in Fig. 2, which confirms the linear behavior in σ at small conductivities and the ∼ σ −1/2 behavior at large conductivities. Furthermore, the dependence on the compactness is monotonic at small conductivities, but it is more complicated at large conductivities, in line with our findings for the amplification factor of massless fields [see discussion around Eq. (19)]. Note that, because ω R ∼ µ V , the small-rotation approximation together with the superradiant condition requires µ V mΩ mΩ K , which implies µ V M 1. To avoid the factor (µ V M ) 8 in Eq. (22) to be exceedingly small, we consider in Fig. 2 a large rotation rate, Ω ∼ 0.9Ω K , although we stress that our results are also valid for smaller values of Ω.
In Appendix A, we discuss a simple model that shares many features with our numerical results.

V. PHENOMENOLOGICAL IMPLICATIONS
We now discuss some potential phenomenological implications of the superradiant instability of stars. We begin with a discussion of the standard (i.e., electromagnetic) conductivity of a neutron star and then generalize the discussion to the conductivity of a hidden sector. Finally, we discuss the implications of the superradiant instability of pulsars for models of dark photons.

A. Conductivity in Maxwell theory
The conductivity of a material can be estimated by a simple Drude model, where n e , e, and m e are the number density, the charge, and the mass of the charge carriers, and τ ft is the mean free time between ionic collisions. The standard charge carriers are electrons and the ionic collisions are between the electrons and protons through electromagnetic interactions. The interaction between electrons and neutrons is small as it proceeds solely through the neutron magnetic moment. More generally, the expression for τ ft will depend on all possible interactions of the electron with protons within the conducting material, where is the proton Fermi temperature, m p is the proton mass, q is the momentum transfer of the collision, and k B is the Boltzmann constant. |M| 2 is the proton-electron scattering matrix element, and in the limit where the electron energies are much smaller than the proton mass, it is given by the Mott formula where k FT is the Fermi-Thomas screening wavenumber for the system. In a neutron star, the protons are much more polarizable than electrons and so k FT corresponds to the contribution of protons alone, i.e.
We assume the star to be electrically neutral, n e = n p . To first order in k FT /k F 1, Together with Eq. (23), this yields [44] σ EM ∼ 2 3 π 3/2 4 (m p n e ) 3/2 where we included the label "EM" to distinguish the above electromagnetic conductivity from the hidden conductivity discussed below. For a typical neutron star with mass density m p n e 10 13 g/cm 3 and T 10 8 K, the above formula yields σ EM 6 × 10 22 s −1 , which in our units translates to σ EM M 10 17 for a typical neutron star mass. In this scenario, where σ EM M 1, we obtain from Eq. (22) a typical instability time scale Therefore, even when σ EM M ∼ 10 17 , the instability timescale can be smaller than a typical accretion time scale, τ Salpeter 4.5 × 10 7 yr. Note that the above estimate was in the regime where the stellar angular velocity is close to the mass-shedding limit, Ω = 0.9Ω K , and the compactness corresponds to the strongest instability [cf. Fig. 2]. In this case, the superradiant instability timescale of neutron stars is actually shorter than that of nearly-extremal BHs [45]. As discussed below, the measured spin of neutron stars is at least a factor of two smaller than the mass shedding limit. Since the timescale will increase for lower angular velocities and for other values of the compactness, Eq. (29) can be taken as a lower limit.

B. Conductivity in Hidden Sectors
We now extend the above discussion to include models of a secluded U (1) [9,46,47] with a massive vector boson X. For this scenario, we will consider the low-energy effective Lagrangian, where F µν is the field strength of the Maxwell vector A µ , X µν is the field strength of the new U (1) gauge boson X µ , m X is the mass of X, and is the kinetic mixing between the two sectors. One can rotate away the kinetic mixing term by working in the mass basis with A µ → A µ + X µ , but this induces a new term j µ X µ in the Lagrangian. The physical consequence is that particles with electric charge also carry a hidden charge e. Therefore, Eq. (11) is modified with σF µν u ν → σF µν u ν +σ X µν u ν , where σ = σ EM to leading order in . For sub-eV m X , the primary constraints on are from stellar production of the vector [48,49], precision tests of electromagnetism [50][51][52], and distortion of the cosmic microwave background (CMB) due to conversion of γ → X [53], which sets < O(10 −7 − 10 −5 ), depending on m X . One can further limit < O(10 −12 −10 −8 ) by constraining the cosmic abundance of X through CMB distortions due to the conversion of X → γ [54,55], while proposed electromagnetic resonator technologies can potentially probe even smaller values of [56]. Thus, the effective conductivity in these models can be much smaller than in Maxwell theory, σ σ EM . In this context it is also relevant to estimate plasma effects, since neutron stars will be surrounded by plasma in various forms. In Maxwell theory, ordinary photons propagating in a plasma acquire an effective mass given by [57] where n plasma is the electron number density in the plasma. In the millicharged cases, we should replace e → e in the above equation. In the context of superradiance [31,32], plasma effects can be neglected as long as ω p µ V . As discussed below, the relevant range of dark-photon masses is µ V ∼ 10 −12 eV. Therefore, if < O(10 −12 ) plasma effects are negligible whenever n plasma < 10 21 cm −3 .
We can also consider a case of a more complicated hidden sector in which the conductivity σ is set by the interactions between particles of opposite U (1) charge, which we denote as hidden electrons and hidden protons, with the hidden electrons serving as the charge carriers (cf., e.g., Refs. [58,59]). Here, j µ → j µ = σ X µν u ν + ρ u µ , which is entirely contained within the hidden sector. The calculation of σ requires the replacement 1 in Eq.
Taking for instance e = 0.01e [59,60], m p = 100 TeV, and assuming that the mass density of hidden protons inside the star is 1% (0.1%) of the mass density of ordinary protons, we estimate a conductivity for hidden electrons σ 10 −16 σ EM (3 × 10 −18 σ EM ), i.e. σ M 10 (0.3). In other words, models of hidden U (1) sectors 1 In the context of superradiant mechanisms, the relevant Compton wavelength of dark photons is much larger than the mean free path of the hidden electrons in the stars. Thus, the mass of the mediator has a negligible effect on the conductivity calculation.
above the TeV scale can have dramatically smaller values of neutron-star conductivity for the hidden electron than that of ordinary electrons, and values σ M ∼ O(1) are allowed. Thus, in our estimates we will consider σ as a free parameter.

C. Instability Time Scale
As discussed in Refs. [19,20], the minimum instability time scale τ ≡ 1/ω I can be estimated by computing the value of µ V which corresponds to the maximum value of ω I . From Eq. (22), dω I /dµ = 0 yields µ min V = 8Ω/9, which corresponds to τ min = 387420489 16777216 The minimum instability time scale is shown in Fig. 3 as a function of the ratio σ/σ EM where σ EM is the typical conductivity of ordinary electrons in a neutron star. As expected, τ min diverges both when σ → 0 and when σ → ∞, and it displays a minimum at σ 10 −17 σ EM , which corresponds to σM 1. Note also that τ min depends strongly on Ω. In Fig. 3, we considered the extreme case Ω = 0.9Ω K , but τ min roughly scales with (Ω K /Ω) 9 . Thus, the time scale for Ω = 0.3Ω K will be roughly 5×10 4 times longer than that shown in Fig. 3.

D. Pulsar-Timing Constraints on Dark Photons
Various arguments [15] suggest that the superradiant instability extracts angular momentum from the central object, spinning it down until the superradiant condition is saturated, µ V ∼ ω R = mΩ (this was recently confirmed by the first numerical simulations 2 of massive vector fields around a spinning black hole [63]). The superradiant instability develops by extracting energy away from the spinning object and depositing it on a bosonic condensate (or a "cloud") outside the object. This cloud has, in general, a time-varying quadrupole moment and will slowly dissipate through emission of gravitational waves. On very long timescales, the end product is an object spinning so slowly that the instability is no longer active.
Because angular-momentum extraction occurs on a time scale τ = 1/ω I , the observation of an isolated compact object with spindown time scale τ spindown excludes superradiant instabilities for that system, at least on time scales τ < τ spindown . Therefore, compact objects for which a (possibly small) spindown rate can be measured accurately are ideal candidates to constrain the mechanism and, in turn, the dark-sector models discussed here.
In Fig. 4, we show the excluded regions in the conductivity vs dark-photon mass plane obtained by imposing τ < τ spindown for three known sources, namely pulsars J1938+2012 [68] and J1748-2446ad [69], and pulsar binary B1957+20 [70]. The first one is representative of a pulsar with an exceptionally long spindown time scale (τ spindown 1.1 × 10 11 yr), but with a moderately large spin (f spin 380 Hz, which corresponds to Ω/Ω K ≈ 0.28 assuming M = 1.4M and M/R = 0.15). The second one is the fastest pulsar known to date (f spin 716 Hz, corresponding to Ω/Ω K ≈ 0.53 for M = 1.4M and M/R = 0.15), but only an upper bound on its spin derivative is available, from which we infer τ spindown > 7.6 × 10 7 yr. The last one is representative of a pulsar with very large spin (f spin 622 Hz, which corresponds to Ω/Ω K ≈ 0.46 again assuming M = 1.4M and M/R = 0.15), but moderately long spindown time scale (τ spindown 3 × 10 9 yr). Furthermore, because our fits for α 1 and α 2 appearing in Eq. (22) are independent Exclusion plots in the σ/σEM vs. mV plane obtained from the measurements of spin and spindown rate of pulsars J1938+2012 (orange) [68] and J1748-2446ad (green) [69], and of the pulsar binary B1957+20 (blue) [70]. In all cases we assumed M = 1.4M and two values of the compactness, namely M/R = 0.15 (solid) and M/R = 0.2 (dotted). The shaded areas correspond to regions excluded by the superradiant instability because τ < τ spindown for a given pulsar (i.e., the pulsar is observed to spin down at much longer rate than that predicted by the superradiant instability in that region of the parameter space). The horizontal dashed line corresponds to where σ = σEM. We only display the region where σ Ω. In the opposite limit, the instability time scale grows as τ ∼ 1/σ [cf. Eq. (22)] and eventually τ > τ spindown for sufficiently small σ, cf. discussion in the main text. The shaded gray region is excluded for σ from distortions of the CMB blackbody from γ → X photon depletion [53].
of Ω only for Ω/σ 1, in Fig. 4 we show only values of the conductivity which satisfy σ Ω. The exclusion plot shown in Fig. 4 is obtained as follows. For a given measurement of the spin frequency of a pulsar, f spin , we can estimate Ω and compute the instability time scale as a function of σ and µ V through Eq. (22). Furthermore, the measurement of a spindown timescale for a pulsar, τ spindown , implies that a faster spindown rate caused by the superradiant instability would be incompatible with observations. Thus, imposing τ < τ spindown yields an excluded region in the σ-m V plane. Fastly spinning pulsars constrain the rightmost part of the σ-m V diagram because the instability requires µ V ∼ ω R < mΩ. On the other hand, pulsars with longer spindown time scale correspond to higher threshold lines in the leftmost part of the σ-m V diagram.
whereṀ is the mass-accretion rate. Since the above time scale is much less than the age of a typical LMXB, many accreting neutron stars in weakly magnetic LMXBs should be observed rotating near the mass-shedding frequency, Ω K /(2π) 1 kHz. The lack of observed systems with f spin 700 Hz has motivated various limiting mechanisms for the maximum spin of a pulsar, many of them involving gravitational-wave dissipation -either through an accretion-induced mass quadrupole on the crust [71], a large toroidal magnetic field [72], or through the excitation of the unstable r-modes [73,74] -and more recently advocating the disk/magnetosphere interaction as leading spindown mechanism [37]. One might wonder whether -besides placing direct constraints on models of dark photons -the superradiant instability of neutron stars could also provide an alternative (albeit exotic) explanation for the spin limit of observed pulsars. For fixed values of σ and µ V , our model predicts that an accreting pulsar in a LMXB (for which the superradiant instability is initially effective) would reach a critical angular velocity such that τ (Ω) = τ spinup in a small fraction of its age. However, because τ spinup τ spindown for the observed pulsars discussed in the previous sections, the threshold line τ = τ spinup is already excluded by pulsar timing. In other words, only models that are already excluded by Fig. 4 would produce a superradiant instability strong enough to overcome accretion at a time scale given by Eq. (34).

VI. DISCUSSION AND FUTURE WORK
The scattering of light by rotating, conducting spheres is a classical problem in electromagnetism, and can lead to superradiant effects. Yet, to the best of our knowledge, a thorough understanding of this problem has not been framed within the context of superradiance. Superradiance in stars may have interesting and important applications in astrophysics and particle physics: stars are made of materials with small but nonvanishing resistivity in the standard Maxwell sector, leading to the amplification of low-energy pulses. In the context of ultralight dark photon models, any nonzero conductivity of stars in the dark sector will lead to superradiant instabilities that drive the star to lower rotation rates. In other words, the superradiant mechanism leads to potentially observable consequences, which can be used to constrain the dark sector.
We have shown that a direct signature of superradiant instabilities in stars is the spindown of pulsars in the presence of ultralight dark photons. As we discussed, existing measurements of the spindown rate of pulsars already place some stringent constraint on models of dark photons and of the hidden U (1) sectors. Although superradiance is typically weaker for stars than for black holes, the spindown rate of pulsars is measured with great precision and it is typically very low (i.e., τ spindown is very long), leading to direct constraints which are much more robust than those coming from mass-spin distributions in the so-called black-hole Regge plane. For example, our preliminary analysis suggests that ordinary models (σ σ EM ) of dark photons with mass m V ∼ 10 −12 eV are excluded by pulsar-timing observations.
There are many interesting follow-up questions to the effect of superradiance in stars. One of them concerns the polar sector of vector perturbations. Previous studies of black hole superradiance show that the vector sector triggers instabilities with much shorter time scales [19][20][21]. If such a result generalizes to conducting stars, the constraints on dark photons will certainly improve. We hope that the promising results of our exploratory study shown in Fig. 4 will stimulate further investigation on this problem, including a complete analysis of the constraints that can be placed on dark-photon models with pulsar timing. From a theoretical perspective, another interesting open issue concerns the functional dependence of the amplification factor on the frequency. Previously, effective field theory approaches have investigated the frequency dependence in the context of black holes [22]. It would be interesting to extend such an approach to stars. Furthermore, in this work we modelled the conductivity with a simple Drude model, in which electrons only scatter with protons. This gives us an order of magnitude of the constraints that one can impose via superradiance, and motivates a more complete calculation (e.g. [75,76]).
In the scenario in which the dark photon couples to Maxwell vectors, superradiance could work in more intricate ways: on the one hand both vectors are superradiantly amplified by the star's material, potentially leading to a stronger effect. On the other hand, Maxwell fields are massless and could easily escape, not being subject to the confinement necessary to create the instability. How exactly the mechanism proceeds depends on this interplay and depends on more detailed calculations. Furthermore, it would be interesting to explore the coupling to plasma. Equation (31) shows that even ordinary photons would acquire an effective mass ω p ∼ 10 −12 eV when propagating in a plasma with electron number density n plasma ∼ 10 −2 cm −3 . This might give interesting superradiant effects for ordinary photons [31,32] or also alter the instability for dark photons if the latter are coupled to plasma sufficiently strongly.
Our analysis also shows that it is, in principle, possible to generalize a number of results in the literature concerning black hole superradiance [15]. For example, for complex, massive vector fields there should exist new stationary solutions describing a star surrounded by a Proca condensate. This would be a natural generalization of the hairy black hole solutions found recently [77][78][79]. Likewise, imprints of superradiance in the luminosity of pulsars or black hole binaries [29,30] should also be present when the companion is a star, instead of a black hole. Finally, the development of the instability will certainly lead to nontrivial gravitational-wave emission. In the black hole case, the emitted signal can be used to impose interesting constraints on the models [15,[24][25][26][27]. On the other hand, stars have typically lower masses than black holes, and it remains to be understood if gravitationalwave emission is relevant in this case.

ACKNOWLEDGMENTS
We are indebted to Leonardo Gualtieri for suggesting a possible connection to the spin limit of pulsars, and to Masha Baryakhtar, Sam Dolan, Mauricio In the membrane paradigm [42], a black hole can be interpreted as a one-way membrane endowed with various properties. In particular, the surface resistivity reads R H = 4π 377 Ohm, so that the surface conductivity iŝ σ = 1/(4π).
Within our framework, a similar model can be investigated by considering a conducting thin shell in vacuum, so that the (volume) conductivity reads σ(r) =σδ(R−r). For simplicity, we consider the Newtonian limit, in which F = B = 1, and restrict ourselves to small frequencies, so that ωζ in Eq. (13) is negligible. In these approximations, axial perturbations reduce to Bessel's equation where m = 1, 0, −1. In the nonrotating case, this result is valid also beyond the small-frequency regime and, interestingly, it agrees exactly with that obtained in black hole perturbation theory [cf. Ref. [15], Eq. (3.103)] upon identification ofσ = 1/(4π) and R = 2M . Thus, a by-product of our analysis is the proof that the black hole membrane paradigm works also for linear electromagnetic perturbations.
The shell toy-model is also useful to understand the results for the instability. Instead of a massive field, we consider a spinning shell of radius R surrounded by a nonspinning perfect conductor of radius R 2 . The characteristic modes of the system can be found by imposing the above junction condition and a(r = R 2 ) = 0. For large values of R 2 /R, the l = 1 fundamental mode reads where γ 0 satisfies tan γ 0 = γ 0 and Υ = γ 0 9 − 16π 2σ2 Ω 2 R 2 − ΩR 2 9 + 16π 2 R 2σ2 Ω 2 (A5) Note that Υ → (γ 0 − ΩR 2 )/9 whenσΩR 1 (thus recovering the superradiant condition, ω R < Ω), whereas Υ ∼ −1/(σ 2 Ω 2 R 2 ) whenσΩR 1. At fixed rotation rate, the peak of the instability occurs atσ ∼ 3/(4πΩR). Therefore -at least qualitatively -this simple model shows the same features that we observe numerically, in particular the fact that the instability decreases asσ → 0 and at very largeσ, and it also informs us on the Ω dependence. Finally, if we substitute R 2 → 1/µ 2 V as discussed in Refs. [15,80], we recover the mass dependence presented in the main text for µ V Ω.
as well as of the identities 4l 2 −1 . By using the above relation, we obtain the following radial equations: Note that Eqs. (B10)-(B12) can be written in the schematic form 0 = A l + a mĀ l + a (Q lPl−1 + Q l+1Pl+1 ) , (B13) 0 = P l + a mP l + a (Q lÃl−1 + Q l+1Ãl+1 ) , (B14) where a is a bookkeeping parameter for the expansion in the angular momentum, A l ,Ā l are linear combinations of the axial perturbations with multipolar index l; similarly, P l ,P l are linear combinations of the polar perturbations with index l.

Axial-led and polar-led perturbations
We expand the axial and polar perturbation functions (schematically denoted as a l and p l , respectively) that appear in Eqs. (B13) and (B14) as a l = a where the first equation is solved to first order in the spin, whereas the second and the third equations do not contain zeroth-order quantities in the spin, i.e. p L±1 = O( a ). The truncation above is consistent because in the axial equations for l = L the polar source terms with l = L ± 1 appear multiplied by a factor a , so they would enter at second order in the rotation.
Similarly, another consistent set of solutions of the same system has a (0) L±1 ≡ 0. The corresponding "polarled" system reads    P L + a mP L = 0 A L+1 + a Q L+1PL = 0 A L−1 + a Q LPL = 0 . (B17) Interestingly, within this perturbative scheme a notion of "conserved quantum number" L is still meaningful: even though, for any given L, rotation couples terms with opposite parity and different multipolar index, the subsystems (B16) and (B17) are closed, i.e. they contain a finite number of equations which describe the dynamics to first order in the spin. Finally, note that the first set of equations in the axialled system (B16) and in the polar-led system (B17) do not involve couplings between axial and polar modes. Once the first set of equations in the system (B16) [or in the system (B17)] is solved, the remaining two equations can be solved separately. Therefore, if one is interested in the linear spin corrections to axial or polar perturbations with a given harmonic index L, one can solve only the first set of equations in the system (B16) or (B17), respectively.
a. Final equations for the axial-led system By using the coefficients given in the Supplemental Material, it is easy to show that the first equation of the system (B16) reduces to Eq. (13) in the main text.

b. Final equations for the polar-led system
The polar-led system is more involved. In general, one of the polar equations fixes the proper charge densityρ EM in terms of the other perturbation functions, even when σ = const. In the Proca case, by using the coefficients in the Supplemental Material, the system can be reduced to three differential equations that can be schematically written as where u 1 ≡ f lm , u 2 ≡ h lm and u 3 = k lm . Note that the first equation above does not contain u 3 . Therefore, it is possible to write a system of two second-order, radial equations for u 1 and u 2 simply by solving the second equation above for u 3 , differentiate it with respect to r, and then using the third equation above to eliminate u 3 . The final result is not shown explicitly and a detailed investigation is left for future work. Note that in the Maxwell case (µ V = 0) the usual gauge freedom can be used to eliminate one spurious degree of freedom.
Consequently, the Maxwell polar sector propagates only one degree of freedom, described by a second-order field equation.