White dwarfs as a probe of exceptionally light QCD axions

We study the effects of exceptionally light QCD axions on the stellar configuration of white dwarfs. At finite baryon density, the non-derivative coupling of the axion to nucleons displaces the axion from its in-vacuum minimum which implies a reduction of the nucleon mass. This dramatically alters the composition of stellar remnants. In particular, the modifications of the mass-radius relationship of white dwarfs allow us to probe large regions of unexplored parameter space without requiring that axions are dark matter.


I. INTRODUCTION
Recent years have seen a resurgence of interest in the physics of the QCD axion, driven by a thriving experimental program in sync with a burst of novel theoretical ideas.One instance is the possibility of relaxing the standard relation between the axion potential and its defining couplings to gluons.Arguably the most exciting outcome of these models is the new set of signals they give rise to beyond the canonical QCD-axion phenomena.Here we present a novel implication associated with exceptionally light QCD axions: white dwarfs (WDs) of a certain size should not exist.This leads to novel bounds in the QCD axion parameter space, Fig. 1.
We consider models of the QCD axion, a pseudo-scalar field ϕ with a coupling to gluons where g s is the strong coupling, f is the axion decay constant and G µν is the gluon field strength.Below the QCD scale the axion obtains the potential [1] where m π ≃ 135 MeV is the pion mass and f π ≃ 93 MeV is the pion decay constant.A parameter ϵ ≤ 1 is introduced to tune the axion lighter than naively expected [2].For a symmetry based realization, see e.g.[3][4][5][6], where note that around the origin and at finite density, these models are well parameterized by Eq. (2).While in the main text we focus on the simple potential Eq. ( 2), a full analysis and discussion of the modifications due to the form of the ZN axion potential [4] is presented in App. C. In vacuum, the axion is stabilized at the origin, thereby solving the strong CP problem.The coupling to gluons Eq. (1) also induces an isospinsymmetric coupling to nucleons where σ N ≃ 50 MeV is the pion-nucleon sigma term.Together with the observed mass, m N ≃ 939 MeV, this coupling can be interpreted as an effective ϕ dependent nucleon mass, m * N (ϕ) ≤ m N .A finite nucleon number density, ρ N ≡ ⟨ N γ 0 N ⟩, implies a non-zero expectation value of ⟨ N N ⟩.In the nonrelativistic limit, these two quantities are approximately equal, ρ N ≃ ⟨ N N ⟩.
Interestingly, a finite ρ N can destabilize the axion from its in-vacuum minimum as soon as σ N ⟨ N N ⟩ > ϵ m 2 π f 2 π , with new minima appearing at ±πf [2,7].Furthermore, once the axion sits in its new minimum, the neutron mass is reduced by Reducing the constituent mass acts as additional binding energy in compact objects.The sourcing of ϕ happens in non-relativistic systems of size R when where λ ϕ ∼ m −1 ϕ (ρ N ) is the typical length scale of the axion at finite density [2,8,9], see also Eq. (13).Typical WD densities fall in the range ρ WD ≃ 10 −4 − 1 MeV 3 .
This implies that for systems with a small characteristic length scale, such as single nuclei, the axion remains stabilized at ϕ = 0 and therefore does not significantly affect forces between nucleons.However, for large and dense systems such as stellar remnants, the sourcing of ϕ implies dramatic changes in their FIG. 1. Constraints and future projections on the axion parameter space for the ϵ-model.Exclusions from modifications of the white dwarf M -R relation are shown in red.Note that the WD bound overlaps with bounds from the Sun in large parts of the region (slightly darker red).The observation of WDs close to the Chandrasekhar limit can further probe the parameter space until the orange dashed line.The solid black line shows the QCD axion with m ϕ f = mπfπ.For reference, we plot f = MP in gray.Further bounds originate from the sourcing in the Sun [2] (blue) and the gravitational wave signal of the NS binary GW170817 [10] (violet), which we both adapted at large f according to numerically inspired O(1) factors, the supernova 1987A [11] (green) , and black hole superradiance [12] (yellow).We would furthermore like to note that the pulsar bound of [2] goes away once all finite gradient effects are properly taken into account (besides lying within a region that is strongly dependent on the neutron star EOS) [9].Finally, we show which parameters lead to a new ground state accessible in neutron stars (dot-dashed purple); for more details see [9].
composition and hence their mass-radius (M -R) relation.A particularly clean and well-studied example of stellar remnants are WDs.The modifications of their M -R relation allow us to probe large regions in the light axion parameter space.For an exploratory study of axions and other particles influencing the stellar structure of neutron stars in the same manner, see [9,13].

II. WHITE DWARF MASS-RADIUS RELATION
WDs balance the force of gravity with the degeneracy pressure of electrons, while almost their entire mass comes from light but non-relativistic nuclei.Due to charge neutrality, the number density of electrons is related to the energy density of nucleons, ε ψ ≃ Y e m N ρ e , where Y e = A/Z is the ratio of nucleons per electron.Since WDs are composed of light nuclei, ranging from helium 4 He to magnesium 24 Mg, the ratio of nucleons per electron is wellapproximated by Y e ≃ 2, see also Fig. 2. The upper and lower bands correspond to the constitutions of light and more heavy nuclei, i.e. 4 He which corresponds to Ye = 2 and 56 Fe corresponding to Ye = 2.15 respectively, while the gray shaded area corresponds to intermediate values.In orange we show the two branches with an axion for ϵ = 10 −11 in the limit RWD ≫ λ ϕ for Ye = 2.The meta-stable branch follows the free Fermi gas line at large radii, while the new ground state phase has much smaller radii.Data points are taken from [14] (turquoise), [15] (blue), [16] (pink), [17,18] (red), [19] (green) and [20] (gray).One can clearly see the gap in the predicted M -R relation that is incompatible with data.

A. Non-interacting gas of electrons and nuclei
The equation of state (EOS) for a degenerate WD can be described by a Fermi gas of non-interacting electrons together with a gas of nuclei.For simplicity, we take positively charged non-relativistic nuclei, which we denote by ψ, with twice the nucleon mass m ψ = 2m N .The pressure is dominated by the electron contribution, p = p e + p ψ ≃ p e , while the nuclei constitute, to a good approximation, the entire energy density, ε = ε e + ε ψ ≃ ε ψ .Due to charge neutrality, ρ ψ = ρ e ≡ ρ, we relate the electron Fermi momentum to the energy density as k F = (3π 2 ε ψ /Y e m N ) 1/3 , and derive the EOS We work in the limit T → 0, which is justified since T /µ e ≪ 1 for typical WDs, where µ e is the electron chemical potential.Temperature shifts the M -R relation up to higher masses, an effect most relevant for the largest and most dilute WDs, see Refs.[19,[21][22][23][24][25][26] as well as Ref. [27] for a review.
The EOS completes the set of equations that describe the balance between the electron degeneracy pressure and gravity, the Tolman-Oppenheimer-Volkoff (TOV) equa-tions [28,29], where G = M −2 P is Newton's constant, M (r) the enclosed mass and all derivatives are taken with respect to the radial coordinate.Solutions with varying central pressures lead to a M -R relationship that agrees well with current data, as shown in Fig. 2.

B. The axion WD system: a new ground state
In the presence of the axion, the full system is described by a free Fermi gas of electrons, an ideal gas of nuclei with a ϕ-dependent mass [30], m * ψ (ϕ) = 2m * N (ϕ), the gravitational field g µν , and the axion ϕ.The gravitational field is sourced by an energy-momentum tensor The first term takes the form of an ideal fluid, T ψϕ µν = Diag (ε, −p, −p, −p) , with where we neglected the sub-leading contributions to the pressure p ψ (ϕ, ρ) ≪ p e (ρ).The second term in Eq. ( 8) contains the contribution of the axion gradient Using Einstein's and the scalar equations of motion, we find the following set of coupled differential equations Eq. (12a) is the static axion equation of motion coupled to gravity, while Eqs.(12b) and (12c) are the TOV equations in the presence of an axion.Note that, we recover the ordinary TOV equations, Eq. (7), in the limit ϕ = 0.While it is possible to numerically solve Eq. ( 12) using the shooting method, there exists a limit in which these equations simplify dramatically.The displacement of the axion at sufficiently high densities costs gradient energy and therefore it only occurs if balanced by the gain in potential energy.This leads to the typical scale on which the axion is displaced to be evaluated at typical WD densities.For R WD ≫ λ ϕ , the field essentially tracks the minimum of the effective in-density potential on stellar scales and is given by the solution to At the same time, the gradient terms in Eqs.(12b) and (12c) are confined to a small transition shell, where the field does not follow its minimum.However, this localized contribution is negligible as long as Therefore, for large systems we can neglect the axion gradient ϕ ′ ≃ 0. As a result, Eq. ( 12) decouples to give the regular TOV equations, Eq. ( 7), in addition to Eq. ( 14).Note that the latter is the same condition as the minimization of the energy density ε(ϕ, ρ) with respect to ϕ. Solutions ϕ(ρ) describe a thermodynamically stable EOS used to solve the regular TOV equations.
Interestingly, if the axion is destabilized in a WD, the energy per particle of the light nuclei ε(ρ)/ρ is not minimized when the nuclei are infinitely separated (ρ → 0), but rather at some finite density ρ * , which can be found numerically.This implies the existence of an energetically favored state of matter at ρ * , where the axion is at ⟨ϕ⟩ = ±πf .This new ground state is in fact reminiscent to strange quark matter [31].Note that the density of the new ground state is slightly larger than the density at which the destabilization occurs, ρ * > ρ c ≡ ϵ m 2 π f 2 π /2σ N .For low densities ρ < ρ c , matter is in a meta-stable state where the classical sourcing of the axion is not preferable.Once ⟨ϕ⟩ = ±πf , there is a range ρ c < ρ < ρ * , where the energy per particle decreases ∂ ρ (ε(ρ)/ρ) < 0, implying a negative pressure.At densities slightly above ρ c , the total pressure turns negative due to the onset of the the axion FIG. 3. The energy per particle ε/ρ as a function of number density (left) and the EOS (right).At low densities ρ < ρc, the system is in its meta-stable ϕ = 0 phase (dark blue).For ρc < ρ < ρ * , the system is unstable i.e. p < 0 (dashed line in left panel).At larger densities ρ * < ρ, the system is in its ϕ = πf phase (light blue), with a new ground state at ε * where p = 0. potential p = p e − V < 0. As ρ increases, V stays constant while p e increases, until finally the system stabilizes at p = p e (ρ * ) − V = 0, see Fig. 3.In this unstable phase the system contracts until it stabilizes in the new ground state.
This instability leads to a gap in the predicted M -R relationship as seen in Fig. 2. The position of this gap is ϵ dependent; the smaller ϵ is, the more the gap is shifted towards small masses and large radii.We use the position of this gap to probe the existence of light QCD axions.
Note that the simplified discussion above is only valid for R WD ≫ λ ϕ .For R WD ∼ λ ϕ we numerically solve the full coupled system in Eq. ( 12) and find that for large values of the axion decay constant, f , and small ϵ, the position of the gap is ϵ independent.This is understood as follows: on the stable branch, the gradient pressure, which is controlled by f , is relevant.If gravity is subdominant, this pressure fixes the central density of the star, ρ(r = 0) > ρ * .The maximal radius is then achieved when the gravitational pressure equals the gradient pressure.For the meta-stable branch, the minimum radius is set by R ∼ lim ϵ→0 λ ϕ .
Finally, in the limit R WD ≪ λ ϕ , the gradient energy is so large that the field cannot move away from its invacuum minimum and therefore has no influence on the structure of WDs.
The data of measured WD masses and radii is scattered broadly between radii of (5000−40000 km), which matches reasonably well with the free Fermi gas description.The notable deviation in mass found at large radii in Fig. 2 is due to finite temperature effects; µ e in these dilute stars is typically smaller, increasing the relevance of T /µ e corrections.Finite temperature effects lead to modifications of the EOS and to a slight modification of masses and radii, but still predict a continuous M -R curve.The same holds for other well-known corrections to the EOS, such as different compositions, electrostatic corrections, or nuclear reactions, see e.g.Refs.[27,39].While nuclear reactions change in the sourced phase e.g.due to a different mass of the pions, this has negligible effect on the static structure of white dwarfs [40] [41].
We perform a simplified statistical analysis to determine the compatibility of the observed WD radii with a gapped radius distribution hypothesis (marginalizing over mass and neglecting small theory systematics).We summarize here its main results, with the full details given in the App.B. For the purpose of the analysis, we calculate the position of the radius gap as a function of ϵ and f , relying both on numerical results, as well as on numericallyverified analytical estimates.In the region f ≪ 10 9 GeV, finite gradient effects are negligible w.r.t. the position of the gap, making it f -independent.For the simplified model described in the main text, we are able to exclude at the 2σ level the following interval in ϵ, see Fig. 1.On the other hand, for the Z N -model, our analysis (see App. C) leads to the exclusion of as shown in Fig. 12.The upper limit is set by the smallest, most massive WDs and our analysis effectively excludes all points in the axion parameter space that cannot predict a WD with a radius smaller than around ∼ 4000 km on the meta-stable branch.The lower limit is sensitive to the largest, extremely low-mass WDs [20].For even lower values of ϵ, the stable branch covers most of the range of observed radii.In this case, although all the observed WDs are sourcing the scalar field, the gap in radius is relegated to extremely large (and potentially unpopulated) WD radii.Note however that this region in parameter space is already ruled out by requiring that no sourcing occurs in our sun [2].
Conversely, in the region ϵ ≪ 10 −20 , or N ≫ 71, finite gradient effects are dominant w.r.t. the position of the gap, making it ϵ-independent.Using solutions of the coupled system, Eq. 12, we are able to exclude at the 2σ level the following interval in f , 5.5 × 10 9 < f /GeV < 1.1 × 10 16 (95% CL) , see Fig.
(1).The upper value represents the limit in which WDs are not large enough to source the scalar field, i.e. λ ϕ ≳ R WD .In the region f ϵ −1/3 ∼ M P , in which both gradient and finite ϵ effects are important, we verify numerically that the sourcing stops at lower values of f .Similarly to the lower bound on ϵ, the lower bound on f is sensitive to the largest, extremely low-mass WDs.We stress again that we do not expect our results to strongly depend on finite temperature effects.

III. CONCLUSIONS
The mass-radius relationship of white dwarfs is wellunderstood and has been observationally tested with increasing accuracy in recent years.We showed how light QCD axions change the structure of WDs, thus predicting the presence of a gap.We used existing data to place novel bounds on their parameter space.We stress that the bounds arising from the existence of a new ground state accessible in white dwarfs, and the corresponding gap in radii, are qualitatively very different than the strategy proposed in Ref. [2], which relies on the change of the properties of nuclei, and the corresponding change in Xray emission, when a (lighter) QCD axion is displaced to θ = π [40].
The QCD axion generically predicts a non-derivative coupling to nucleons.At finite baryon density this coupling can destabilize the axion from its in-vacuum minimum.If sourced, the non-zero axion expectation value reduces the mass of nucleons.For a large region of the parameter space, this leads to a new ground state of matter, which has less energy per particle than infinitely separated nucleons.If accessible in WDs, this drastically changes their M -R relation.Since the axion is sourced by the WD, this does not rely on the axion contributing to the dark matter relic abundance.
More precise tests of the WD M -R curve using the recent Gaia DR3 are expected in the near future and will further probe the parameter space of light QCD axions.
As a consequence of the new ground state of matter, we predict new small self-bound objects held together by the gradient pressure of the axion.These objects could give rise to novel signatures of exceptionally light QCD axions down to the QCD axion line.
where ρ 0 is the number density at the core (r = 0).The gradient pressure in the sourced phase (i.e.stars on the stable branch) is given by where ρ R is the number density at the edge (r = R − λ ϕ ≈ R) of the WD.The first line in Eq. (A3) represents the thin wall limit λ ϕ ≪ R, in which the gradient pressure is exerted at a small transition region at the edge of the star.
In the last step we use the definition of the in-medium wavelength of Eq. ( 13) assuming a negligible contribution from the scalar potential and neglecting O(1) factors.The second line in Eq. (A3) is the opposite regime λ ϕ ∼ R, where the gradient pressure is delocalized and is spread throughout the star.This is the typical edge case configuration in which the star is barely large enough to source the axion.In the unsourced phase, i.e. stars on the metastable branch, p grad = 0. See Ref. [9] for more details on the derivation of Eqs.(A2) and (A3).
On the other hand, we define the outwards pressure at the core balancing p 0 , see Eq. A1, as the contribution of the electron gas and the scalar potential, which can be written analytically in the non-relativistic (NR) and ultrarelativistic (UR) limits as ) where in the sourced phase we have by definition V (ϕ) = p e (ρ * ), while in the unsourced phase the scalar potential FIG. 4. Left panel: the ratio between the numerical results and the analytical estimates for the radii which correspond to the edges of the radius gap, as a function of ϵ.In (dashed) red we plot the ratio for R stable max divided by the NR (UR) estimate given in Eq. (A6).In green we plot the ratio for R meta min , divided by the NR estimate give in Eq. (A7).In all cases we match the O(1) prefactors to the numerical results.Right panel: the radius gap defined by R stable max (red) and R meta min (green) as a function of ϵ.Solid curves indicate region where numerical results are used while the dashed curve indicates where extrapolation (using the verified analytic estimates) is used.The gray region corresponds to the observed radii range of WDs.
Let us start by estimating the minimal radius on the meta-stable branch, which we denote by R meta min , and the maximal radius on the stable branch, which we denote by R stable max , in the negligible gradient regime, where ∆p grav ≫ p grad .In this limit, R stable max is the radius of the largest approximately constant energy density configurations.Therefore, we estimate it by setting ρ 0 ≈ c ρ * ≈ c ϵ 3/5 (m e m 2 π f 2 π ) 3/5 , where c ∼ O(1), and solving for R. The contribution to the pressure from V (ϕ) = p e (ρ * ) is neglected since it is at most of the same order as p e (ρ 0 ), and would therefore have at most an O(1) effect on the final result.We find where for brevity we denoted m 2 π f 2 π ≡ Λ 4 QCD .We omit the weak dependence on c, which amounts to an O(1) prefactor.In the left panel of Fig. (4) we compare the analytic estimates to the numeric results.We find that the NR estimation is in excellent agreement with the numerical results for ϵ ≲ 10 −13 (red curve).For larger values of ϵ, the smaller minimal radii on the stable branch correspond to denser configurations, where relativistic corrections become important.Thus, for ϵ ≳ 10 −11 the UR estimation agrees with the numerical results (red dashed curve).
The edge of the meta-stable branch R meta min is found by taking ρ 0 ≈ ρ c ≈ ϵΛ 4  QCD /(2σ N ) and solving for R in the NR approximation, resulting in A similar UR approximation is straightforward to derive.However, it is only valid for R ≪ M P /(m N m e ) ∼ 5000 km, which is outside our range of interest.In the left panel of Fig. ( 4), we compare the analytic numerical results of R meta min (ϵ) to the analytic estimate.We find good agreement in most of the calculated region, namely for ϵ ≲ 10 −9 .For larger values, relativistic corrections start becoming important and the UR approximation begins to degrade.
The radius gap in negligible gradient regime is plotted in the right panel of Fig. (4) as a function of ϵ.For the purpose of the analysis of Sec.(B), solid curves indicate the regions where we use our numerical results, while dashed lines indicates regions where extrapolation, based on the verified numerical estimate, is used.The gray region corresponds to the observed radii range of WDs.
In the region of parameter space where ϵ is negligibly small (to be determined below), the position of the radius gap is determined by finite gradient effects.On one side, the edge of the meta-stable branch indicates when a region of size λ ϕ with above-critical density is formed, which leads to an instability.On the other side, the largest configurations on the stable branch are those in which the gravitational pressure begins to dominate over the gradient pressure exerted at the edge of the star [9].
First we find R stable max for lower values of f where the thinwall approximation holds.by taking ρ 0 = ρ R ≡ ρ eq > ρ * , where ρ eq is found by solving ∆p grav (ρ eq ) = ∆p grav (ρ eq ) for ρ eq .We than plug ρ eq into Eq.(A1), and using the NR FIG. 5. Left the ratio between the numerical results and the analytical estimates for the radii which correspond to the edges of the radius gap, as a function of f .In red and dashed-red we plot the ratio for R stable max divided by the NR and UR estimates given in Eq. (A8) and Eq.(A9), respectively.green we plot the ratio for R meta min , divided by the NR estimate give in Eq. (A10).In all cases we match the O(1) prefactors to the numerical results.Right panel: the radius gap defined by R stable max (red) and R meta min (green) as a function of f .Solid curves indicate region where numerical results are used while the dashed curve indicates where extrapolations (using the verified analytic estimates) are used.The gray region corresponds to the observed radii range of WDs.approximation of Eq. (A5) solve for R and find where we again neglect the contribution from V (ϕ) as an O(1) correction at most.We compare this estimate with the numerical results in Fig. (5) (red curve).We find is is consistent with the numerical results in the region f ≪ 10 15 GeV.Above these values of f , the thin-wall approximation breaks down and R stable max (f ) can be estimated using the λ ϕ ∼ R expression for p grav and the UR expression for the electron pressure, which gives us We find this estimate consistent with the numerical results in the region f ≫ 10 15 GeV, see dashed curve in Fig. (5).
The edge of the meta-stable branch R meta min is found by first finding the critical density for which the while size of the star is of the order of the scalar in-medium wavelength, namely by solving λ ϕ (ρ) = f / √ ρδm N = R for ρ and plugging the result in Eq. (A1) using the NR approximation for the electron gas.We find We find this estimate consistent with the numerical results, see green curve in Fig. (5).The deviations from the analytical estimate at large value of f is expected, since for smaller and denser stars relativistic corrections to the EOS become increasingly larger.
To conclude we note that for the whole {ϵ, f } plane, we define the maximal radius of the gap as where the two radii coincide R meta min (ϵ) ∼ R meta min (f ) around the curve defined by M P ∼ M P , (A12) using the NR estimations for both expressions.This defines (a posteriori) the ranges of validity for the negligible gradient and negligible ϵ approximations for the determination of R meta min , e.g. the negligible gradient is a valid approximation in the region of parameter space where f ϵ −1/3 ≪ M P .Around f ϵ −1/3 ∼ M P , we computed a numerical solution using the appropriate {ϵ, f } values in order to obtain R meta min (ϵ, f ).In a similar fashion, we define for the whole {ϵ, f } plane the minimal radius of the gap as where the two radii coincide R stable max (ϵ) ∼ R stable max (f ) around the curve defined by using the NR estimations for both expressions.This defines (a posteriori) the ranges of validity for the negligible gradient and negligible ϵ approximations for the determination of R stable max .
FIG. 6. Left panel: χ(ϵ), as defined Eq. (B1) in the negligible gradient limit, normalized to the effective number of degrees of freedom N − 1 = 294, as a function of ϵ.Right panel: The p-value as a function of ϵ.For reference, we plot the 2σ threshold, equivalent to p = 0.05, as a gray dashed line.FIG. 7. Left panel: χ(f ), as defined Eq. (B1) in the negligible ϵ limit, normalized to the effective number of degrees of freedom N − 1 = 294, as a function of f .Right panel: The p-value as a function of f .For reference, we plot the 2σ threshold, equivalent to p = 0.05, as a gray dashed line.

Appendix B: Statistical analysis and bounds
The goal is to quantitatively determine how compatible are the observed WD radii with a gapped distribution.Our working assumption is that the variance in the observed mass can be explained by varying other important properties of WDs, such as temperature and composition, which for simplicity we kept fixed.Therefore our focus is strictly on the radius distribution, and in order to determine the bounds on ϵ and f , we perform a 1D goodness-of-fit test on the radius axis.Given the central values from a com-bined dataset of N = 295 observed WD radii {r i } and their corresponding uncertainties {σ i }[14-20], we calculate the sum of squares, which we denote by χ, for each point in the {ϵ, f } plane with the distance function We use the numerically-calculated values for {R stable max (ϵ, f ), R meta min (ϵ, f )} where available, and other- In Fig. (6) and Fig. (7) we plot the results of our statistical analysis in the negligible gradient and negligible ϵ limits, respectively.In the left panels we plot χ normalized to the effective number of degrees of freedom N − 1.
As a rough estimate, the range in which χ > 1 is considered incompatible with the gapped radii distribution hypothesis.A more refined statement can be made by calculating the corresponding p-values for each value of ϵ of f , shown in the right panel of Fig. (6) and Fig. (7), respectively.We are able to exclude at the 2σ level the following interval in ϵ 2 × 10 −20 < ϵ < 2 × 10 −7 (95% CL) . (B3) We are able to exclude at the 2σ level the following interval in f , 5.5 × 10 9 < f /GeV < 1.1 × 10 16 (95% CL) .

(B4)
Appendix C: ZN axion The potential of lighter than expected QCD axions due to a Z N discrete symmetry reads, in the large N limit [4], where N ≫ 1 is odd and z = m u /m d .The corresponding axion mass reads We can map the Z N axion to the potential in Eq. ( 2) by requiring equal axion masses.We find the relation While this relation fixes the axion mass to be the same in both models, the value of the potential at θ = π is different.We find where V (θ) is the potential from Eq. ( 2).Clearly, there is a large suppression for N ∼ O (10), which is the scenario we are interested in.At finite density ρ = k 3 F /3π 2 , the potential for the Z N axion is given by [4] V where β = 4z/(1 + z) 2 .The critical density, at which the minimum at the origin is destabilized, i.e. ∂ 2 θ V N (0, ρ c ) = 0 is which matches the critical density for the potential in Eq. ( 2), ρ c = ϵm 2 π f 2 π /2σ πN , with ϵ given by Eq. (C3).At densities below the critical density, V N (θ, ρ) has (N + 1)/2 minima in the interval θ ∈ [0, π], which we label by k = 1, . . ., (N + 1)/2.
At non-zero densities, minima with k > 1 have lower energy than the minimum at θ = 0 (or k = 1), but get destabilized at lower densities than the minimum at θ = 0.In particular, we would like to point out that θ = π is initially not the lowest energy configuration since it remains a maximum for these intermediate densities, hence it is not preferred over e.g.θ = (N − 1)π/N .We solve ∂ 2 θ V N (π, ρ π c ) = 0, and find that the density at which θ = π becomes the lowest energy configuration is FIG. 9. Left: Pressure as function of the energy density for N = 31.Rainbow color-coding marks the solutions to the (31 + 1)/2 = 17 minima from the one closes to π (most red) to the one at zero (most violet).The red curve never experiences negative pressure.
While the pressure decreases due to the displacement of θ from 30π/31, the Fermi pressure, due to electrons becoming relativistic, keeps the total pressure positive.Right: M -R curves for N = 31 in the same color-coding.At densities ρ π c ≤ ρ ≤ ρ c , two minima exist, at θ = 0 and θ = π, while at densities ρ > ρ c the only remaining minimum is at θ = π.
We show an N = 7 potential for zero and sub-critical densities in the left panel of Fig. 8.In the right panel we show solutions θ(ρ) that sit in the (N + 1)/2 minima between θ = 0 and θ = π for N = 7.
In the negligible gradient limit, we can derive the EOS as described in the main text, i.e. by solving Eq. ( 14).At subcritical densities we find (N + 1)/2 independent branches of the EOS, one for each meta-stable minimum.The field value of θ increases with k and so does the reduction of the mass.The lowest energy configuration is therefore the branch that starts at θ = (N − 1)π/N at low densities.
Solving the system for all branches of the EOS, we generally find a qualitatively distinct behaviour at small and large values of N .
1. Small N At small N , i.e.N ≤ 31, the electron Fermi pressure dominates over the negative contribution from the potential, such that the total pressure stays positive.This is because electrons become relativistic at the relevant densities, since the lower the value of N , the larger the critical density, see Eq. (C6).
Therefore, we do not find a gap in the EOS, which is qualitatively different from what we found for the potential in Eq. (2).We show the EOS, i.e. p(ε), in the left panel of Fig. 9 for N = 31.Note that the EOS for the branch that sits closest to π (red curve) does not experience a thermodynamic instability for any density, given that the pressure stays a monotonically increasing function of the energy density.Nevertheless, the onset of the non-zero potential once the field starts to be displaced reduces the pressure significantly, leading to a much softer EOS at large densities than the EOS for θ = 0.In the zero gradient limit, the M -R curve is readily found by solving the regular TOV equations Eq. ( 7) with the prescribed EOS.As expected from the EOS, we find continuous M -R curves.For very low values of N , sourcing happens only for the densest WDs.We show the family of M -R branches in the right panel of Fig. 9 for N = 31.
Since we are agnostic about the formation mechanism of the star, and thus of the branch it ends up after formation, a χ 2 in radii is not enough to exclude this prediction.In principle an exclusion might be possible for the highest N in this regime, i.e.N = 31.This requires a more sophisticated statistical analysis than done for the gapped curve and is beyond the scope of this work.

Large N
At large N , i.e.N ≥ 33, we find that the negative contribution of the potential dominates over the electron Fermi pressure.This implies a thermodynamic instability and a gap in the EOS, associated with a new ground state with density ρ * > ρ π c .Contrary to the case of Eq. ( 2) (see Fig. 3), this gap can be covered by the θ = 0 branch.This is because the θ = (N − 1)π/N minimum disappears (and the one at θ = π appears) before the minimum at θ = 0 is destabilized, as follows from Eq. (C7).In Fig. 10, for N = 33, it is shown that the region in which the instability occurs is covered for all densities by the meta-stable (θ = 0) branch.For such N , the new ground state density ρ * lies below the critical density, as can be seen in the left panel of Fig. 10.For N ⩾ 39, we find instead ρ * > ρ c , a scenario that parallels that discussed in the main text.
Analytic estimates for the density of the new ground As we have seen from the EOS, for larger values of N we find negative pressure phases in the (N − 1)π/N branch.We therefore expect constant density self-bound objects (SBOs) and a gap in radii.In the left panel of Fig. 10 we show the EOS for N = 33, which is the lowest N for which we find a new ground state.In the right panel of Fig. 10 we show the corresponding M − R curves.As can be seen, the new ground state leads to SBO solutions i.e. the red curve that connects to M WD = 0 and R WD = 0.This branch is disconnected by a gap in radii from the meta-stable branches at larger radii.
In the left panel of Fig. 11, show the gap in radii as a function of N for N ≥ 33 in we show the minimal radius of the meta-stable branch R meta , where for very large radii we used the estimate R meta min ≃ 10 4 2m N c 10 6 g cm −3 −1/6 km, (C11) see e.g.[39].In red we show the maximal radius of the stable branch R stable max for which we use our numerical results and a similar analytic estimate used.In the right panel of Fig. 11, we show the results of the χ 2 analysis, described in App.B. As can be seen, at the 95% confidence level, we are able to exclude 33 ≤ N ≤ 69.Last, we expect gradient effects to be of similar importance as in the detailed study above for only one minimum.
Thus we expect also for the Z N axion that the bound shuts down similarly around f ∼ 10 16 GeV.With an analogous statistical analysis as done above, we come to a slightly modified bound, as is shown in Fig. 12.

FIG. 2 .
FIG.2.White dwarf M -R relation with light QCD axions.Free Fermi gas of nuclei and electrons without an axion (black).The upper and lower bands correspond to the constitutions of light and more heavy nuclei, i.e.4He which corresponds to Ye = 2 and 56 Fe corresponding to Ye = 2.15 respectively, while the gray shaded area corresponds to intermediate values.In orange we show the two branches with an axion for ϵ = 10 −11 in the limit RWD ≫ λ ϕ for Ye = 2.The meta-stable branch follows the free Fermi gas line at large radii, while the new ground state phase has much smaller radii.Data points are taken from[14] (turquoise),[15] (blue),[16] (pink),[17,18] (red),[19] (green) and[20] (gray).One can clearly see the gap in the predicted M -R relation that is incompatible with data.

FIG. 10 .
FIG.10.Left: Pressure as function of the energy density for N = 33.The rainbow color-coding marks the different branches of the EOS for the (33 + 1)/2 = 17 minima of θ with the ones closer to π being more red.Importantly, here the branch closest to π experiences an instability.Right: M − R curves for N = 33 in blue (SBO) and orange (meta-stable) and for the corresponding ϵ (Eq.(C3)) in green (SBO) and red (meta-stable).

FIG. 11 .
FIG. 11.Left: Radius gap as a function of N .Gray shaded area marks radii populated by data.Right: χ 2 values over data points in the gap as a function of N .

FIG. 12 .
FIG. 12. Constraints and future projections on the axion parameter space for the ZN -model.Same as in Fig. 1.