Axion sourcing in dense stellar matter via CP-violating couplings

We study axions in the high density regime, assuming that they interact only with neutrons, protons, and electrons. The CP-conserving axion Lagrangian density reads [28, 34, 35]

In Eq. (1), $a$ is the axion field, $\mathcal{V}(a)$ denotes the axion potential, $F_{\mu\nu}$ is the electromagnetic strength tensor and $\tilde{F}_{\mu\nu}$ its dual. The index $j$ runs over the fermionic degrees of freedom interacting with axions (neutrons, protons, and electrons in our model, denoted by the fields $\psi_{j}$).

The quantity $i\bar{\psi}_{j}\gamma^{5}\psi_{j}$ has contributions from the divergence of the spin density in the non-relativistic limit, which in turn depends on the magnetic field strength, and to a lesser extent from the anomalous magnetic moment of the fermion (in case there is a non-vanishing electric field). Since the CP-conserving term couples the axion to the spin density of particles, one would need a high degree of polarization for the CP-conserving term to be comparable to the CP-violating term (discussed below). In the magnetar scenario, we expect $B\approx 10^{14}$ G and $T=10^{8}-10^{9}$ K = $10-100$ keV, so that $e\hbar B/m_{n}c\approx 1$ keV $\ll~{}kT$. Thermal fluctuations wash out any significant spin-polarization of nucleons, while for electrons one has $e\hbar B/m_{e}c\approx$ few MeV. However, the CP-conserving term is several orders of magnitude smaller than the CP-violating term. Thus, in the following, we focus on axion sourcing only via the CP-violating axion-fermion interactions.

If there are sources of CP-violation in the Standard Model, axions acquire CP-violating couplings leading to the following Yukawa terms [28, 34, 35]

where $\bar{g}_{aj}$ are the CP-violating coupling constants of axions with nucleons. Such couplings have been recently considered in the context of compact stars in [29], which studied the effect on the mass-radius relation of NSs in the case of a light scalar linearly coupled to standard matter. Note that the QCD axion interacts only derivatively with electrons (i.e. only with CP-conserving interactions). We assume that ALPs participate in the same interactions as QCD axions to compare the two cases directly.

We can now derive the Klein-Gordon (KG) equation for the axion field, reading [28]

where $m_{a}$ denotes the axion mass, and $\mathbf{E}$, $\mathbf{B}$ denote respectively the electric and magnetic fields (we neglect the terms proportional to $g_{aj}$). Since axions are expected to be light ($m_{a}\lesssim 0.1$ eV), their typical Compton wavelength is macroscopic, and they interact with the mean-field values in dense matter of the field bilinears (i.e. $\langle\bar{\psi}_{j}\psi_{j}\rangle$).

In the following, we use the notation $n_{S}=\sum_{j}\langle\bar{\psi}_{j}\psi_{j}\rangle$ for the scalar density in the KG equation. This can be obtained from [36]

where $M^{*}_{N}$ denotes the nucleon effective mass in dense matter, and $E_{N}^{*}(k)=\sqrt{k^{2}+M^{*2}_{N}}$. The prefactor $\gamma$ is the spin-isospin degeneracy factor [36], and $k_{F}$ denotes the Fermi momentum.

For NSs, we adopt the GM1A equation of state (EoS) [37], matched with the SLy4 equation of state in the crust [38]. We focus on a NS model with mass $M_{NS}=1.49\ M_{\odot}$ and with radius $R_{NS}=13.8$ km, containing only nucleons and leptons. For sufficiently low densities however (typically in the NS crust), the ground state of dense matter is given by nuclei, and we approximate the scalar density with its non-relativistic limit, i.e. the baryon density, similarly to the case of axion fields sourced by test masses in laboratory experiments (e.g. [39]). For WDs, we use a relativistic polytropic equation of state. We adopt a WD model with mass $M_{WD}=0.7\ M_{\odot}$ and $R_{WD}=7.98\times 10^{3}$ km. In this case, $n_{S}$ in Eq. (4) is approximated with the baryon density.

For ALPs, the axion mass $m_{a}$ and the couplings $\bar{g}_{aj}$ can be treated as free parameters. Due to the high electrical conductivity of both NSs and WDs, the electric field is weak, with $E\approx(v/c)B$ (where the drift velocity of the electrons is $1\ \textrm{km}\ \textrm{Myr}^{-1}\lesssim v\lesssim 10^{3}\ \textrm{km}\ \textrm{Myr}^{-1}$ in the crust). In such conditions, the electromagnetic source term ($\propto\mathbf{E}\cdot\mathbf{B}$) in the KG equation is subdominant, and we neglect it. Note however that the electromagnetic term can give relevant contributions in the magnetosphere of both NSs and WDs [40]. A full solution of the KG equation including the electromagnetic term requires accurate modeling of magnetospheric fields and is left for future work. We compare however our results obtained without the electromagnetic source terms with the ones found in [40] in the following sections.

To discuss in detail how the solutions of the axion field scale with the relevant parameters, it is convenient to introduce the energy scale¹¹1The following results are valid for both NSs and WDs, and we use the generic notation $M,R$ to denote the mass and radius of a compact star. When considering specific NS and WD models, we write $M_{NS},R_{NS}$ and $M_{WD},R_{WD}$ respectively.

where $R$ is the star’s radius and the scalar density is evaluated at the center of the star ($r=0$). We seek solutions of the stationary KG equation in terms of the dimensionless variable $\varphi=a/e_{a}$. Expressing the radial coordinate in units of $R$ (that is, hereafter $r\equiv r/R$), the KG equation can be written as

where $f(r)$ is the source term divided by $e_{a}$, which by definition decreases monotonically from $f(r=0)=1$.

We solve Eq. (6) with the appropriate boundary conditions both, numerically and exactly. The exact formal solution is ²²2Another useful way to write the solution is in terms of the Green’s function: $$\varphi(r)=-\frac{1}{2m_{a}Rr}\int_{0}^{1}dr^{\prime}r^{\prime}f(r^{\prime})\left(e^{-m_{a}R|r^{\prime}-r|}-e^{-m_{a}R|r^{\prime}+r|}\right)$$

with

Since the source function $f(r)$ vanishes for $r\geq 1$, the solution outside the star reduces to

The constant $A$ can be evaluated for any given density profile and sets the surface value of the axion field. We can also evaluate the limit to obtain the value at the origin $\varphi_{0}$: ³³3 Since $f(r^{\prime})$ is a decreasing, positive defined function with $f(0)=1$, we have $|\varphi_{0}|\leq\int_{0}^{\infty}xe^{-x}dx=1$

One can recover the original dimensions for the axion field by multiplying $\varphi_{0}$ by the energy scale $e_{a}$.

It is interesting to discuss how the amplitude of the field behaves in the different axion-mass limits. In the interior of the star, we find

where $\beta(r)\lesssim\frac{1}{2}$ is a function of order unity; for the particular NS model employed in this work we have $\beta(r=0)=0.32$. Note that for $m_{a}R\gg 1$ the maximum amplitude of the field scales as $(m_{a}R)^{-2}$.

Here we take $\bar{g}_{aN}=10^{-23}$, which is one order of magnitude smaller than the lower bounds reported in [28, 34] (see also references therein and the limits reported in [41]) for light ALPs. In NSs, this results in $e_{a}\approx 7.9\times 10^{13}~{}\mathrm{GeV}$.

In Figure 1 we consider the NS case and show results of the stationary KG equation varying the parameter $m_{a}R_{NS}$ in the range $1\leq m_{a}R_{NS}\leq 100$, corresponding to $10^{-11}\lesssim m_{a}/\textrm{eV}\lesssim 10^{-9}$. In the left panel of Figure 1, we compare the radial profiles of the dimensionless axion field $(m_{a}R_{NS})^{2}|\varphi(r)|$. The dotted curve represents the normalized scalar density profile of the NS model adopted in this work, which is almost indistinguishable from the axion field profile for $m_{a}R_{NS}\gtrsim 20$. The amplitude of the field in physical units (GeV) is plotted in the central panel, which shows that the field decays exponentially outside the star. We also note that further decreasing the dimensionless axion mass $m_{a}R_{NS}$ does not lead to an arbitrarily large axion field. Instead, it converges to a solution that tries to minimize gradients (which become the dominant contribution to the total energy in the regime $m_{a}R_{NS}\ll 1$) inside the star. In the intermediate regime, $m_{a}R_{NS}\lesssim 1$, the solution of Eq.(6) finds an optimal balance between gradient and potential terms. In the opposite limit ($m_{a}R_{NS}\gg 1$), the gradient terms can be neglected and the axion field follows the relation reported in Eq. (13). In the right panel, we show the scaling of the central value ($\varphi_{0}$) with $m_{a}R_{NS}$ and find that $\varphi_{0}\approx 0.32$ for $m_{a}R_{NS}\lesssim 1$ and $\varphi_{0}\propto(m_{a}R_{NS})^{-2}$ for large values. We show a wide range of $m_{a}R_{NS}$ for completeness, although values of $m_{a}$ in the range $10^{-13}\ \textrm{eV}\lesssim m_{a}\lesssim 10^{-11}$ eV (that is, $10^{-2}\lesssim m_{a}R_{NS}\lesssim 1$) are, in principle, excluded by black hole superradiance arguments [42, 43].

We now determine the energy stored in the ALP field sourced by NSs. The energy density reads $\mathcal{H}=(m_{a}a)^{2}/2+(\nabla a)^{2}/2$. From the behavior of the solutions described in Eq. (13), we find

at $r=0$. Note that the maximum value of the energy density is smaller than $\left(e_{a}/R_{NS}\right)^{2}$, which is also of the same order as the interaction energy $\mathcal{H}_{int}=\bar{g}_{aN}n_{S}(r=0)a\approx e_{a}a/R_{NS}^{2}$. Using Eq. (5), we have

which is comparable to the magnetic energy density of NSs (given by $B^{2}/8\pi$) with inferred fields in the range $B\in[10^{12},10^{13}]$ G for $m_{a}R_{NS}\lesssim 1$. For $m_{a}R_{NS}\gg 1$, the axion energy density scales as $(m_{a}R_{NS})^{-2}$.

The left panel of Figure 2 compares different radial profiles of the energy density for different choices of the parameter $m_{a}R_{NS}$. The typical energy density stored in the axion field is of order $10^{-2}\left(e_{a}/R_{NS}\right)^{2}\approx 10^{24}~{}\mathrm{erg~{}cm}^{-3}$ for $m_{a}R_{NS}$ of order unity, and scales as $(m_{a}R_{NS})^{-2}$ for heavier masses, as discussed above. The right panel of Figure 2 shows the total energy integrated over the star volume

which is a quantity of the order of $10^{-1}e_{a}^{2}R_{NS}\approx 7\times 10^{43}$ erg for the NS model considered here and $m_{a}R_{NS}=1$, which is comparable to the total magnetic energy in NSs with poloidal-dipolar fields and inferred fields of the order $B\approx 10^{13}$ G [44]. For completeness, we also include in the figure the total energy of the axion field outside the star

We note that the maximum energy available outside the star is of the order of $10^{44}$ erg, which can be relevant for magnetospheric physics but still irrelevant for long-range gravitational effects (compared to $M_{\odot}c^{2}$).

One can obtain a rough estimate of the total axion energy for the axion field sourced by WDs as follows. Since in WDs one has $m_{a}R_{WD}\gg 1$, we obtain

where we used $m_{a}=10^{-11}$ eV and $n_{s}=10^{30}$ cm^-3. Therefore, the energy of the axion field in WDs is $E_{in,WD}\approx 10^{37}\mathrm{~{}erg}$. While the axion energy is considerably lower than the NS case, it is still comparable to the magnetic energy expected in some white dwarfs. The magnetic fields of WDs are in the range $10^{3}\lesssim B/\textrm{G}\lesssim 10^{9}$, with approximately $90\%$ of the WD population having fields below $10^{6}$ G [45, 46]. This leads to an estimated magnetic energy of $10^{38}\left(B/10^{6}{\mathrm{G}}\right)^{2}$ erg assuming for simplicity a poloidal-dipolar magnetic field.

We now turn to the QCD axion. We calculate finite density corrections to the axion mass following [22, 26, 25]. The QCD axion potential in vacuo reads [47]

(we subtract a constant so that the potential vanishes when $a=0$), where $m_{\pi}=135$ MeV and $f_{\pi}=93$ MeV denote the neutral pion mass and pion decay constant respectively, $m_{u}$ and $m_{d}$ are the up and down quark masses in vacuo respectively and $f_{a}$ is the axion decay constant. The QCD axion potential in vacuo $\mathcal{V}_{0}$ in Eq. (21) is proportional to the vacuum quark condensate $\langle\bar{q}q\rangle_{0}$ (where $q$ denotes the quark field) via the Gell-Mann-Oakes-Renner relation

In dense matter, the expectation value of the quark condensate $\langle\bar{q}q\rangle_{n}$ changes with respect to the in-vacuo condensate. Using the Hellmann-Feynman theorem [48], one gets

where $m_{q}$ denotes the average, bare mass of up and down quarks and $\mathcal{E}$ is the energy density of the system, which depends on the EoS. For the GM1A EoS and the NS model studied in this work, one obtains

In Eq. (24), we take $\sigma_{N}=59$ MeV [48, 26], and use $M_{N}$ to denote the bare mass of nucleons. $g_{\sigma N}$ and $m_{\sigma}$ are the coupling constant of the $\sigma$ meson with nucleons and the mass of the $\sigma$ meson respectively, and $g_{3}$ and $g_{4}$ denote the coupling constants for cubic and quartic self-interactions of the $\sigma$ meson respectively. The last two lines in Eq. (24) include the contributions of the vector mesons $\omega$ and $\rho$, with masses $m_{\omega}$ and $m_{\rho}$ and couplings $g_{\omega N}$ and $g_{\rho N}$ to nucleons respectively. For the quantities $\chi_{l}$ (with $l=\sigma,\omega,\rho$), we follow [48] and set $\chi_{l}\approx m_{l}/M_{N}$. $n_{p}$, $n_{n}$, and $n_{N}$ represent the proton, neutron, and nucleon number densities respectively (in particular, $n_{N}$ coincides with the baryonic density in the absence of hyperons). For simplicity, in the following we use the notation $\langle\bar{q}q\rangle_{n}/\langle\bar{q}q\rangle_{0}=1-F_{n}(r)$.

Eq. (24) can be used to determine how finite density corrections change the QCD axion potential $\mathcal{V}(a)$ in dense matter, i.e.

The effective axion mass is then

with $z=m_{u}/m_{d}=0.48$. Given $m^{*}_{a}$, the radial profile of the QCD axion field in dense matter is determined using Eq. (25) to solve the KG equation.

We now calculate the magnitude of the QCD axion field in NSs. The CP-violating couplings for the QCD axion are expected to be in the range [33]

The resulting coupling strength $\bar{g}_{aN}$ is weaker than in the ALP case studied above, the sourced axion field is in general smaller, and the corresponding phenomenology is harder to test (as discussed in Section IV).

The natural energy scale is now $f_{a}$. We define $\phi=a/f_{a}$, and the KG equation becomes:

with $\phi=a/f_{a}$ and $e_{s}=\sqrt{e_{a}f_{a}}$.

In the following, we take $\bar{g}_{aN}=10^{-12}~{}\mathrm{GeV}/f_{a}$ for the QCD axion, i.e. the upper bound in Eq. (27). With our choice of the CP-violating couplings, the energy scale $e_{s}$ is fixed to $e_{s}\approx 2.8\times 10^{12}~{}\mathrm{GeV}$.

We also note that the coefficient of the last term on the right-hand side of Eq. (28) is independent of the scales $R$ and $f_{a}$, i.e.

For the EoS adopted in this work, $(m_{a}^{*})^{2}$ remains positive ($1-F_{n}>0$) and the solution can be well approximated expanding the QCD axion potential in Eq. (21) and using only the $(a/f_{a})^{2}$ contribution. Then, as with ALPs, we can consider the limiting case $m_{a}^{*}R\gg 1$. For $m_{a}^{*}R\gg 1$, since the solution to the KG equation satisfies $\phi\ll 1$ ($a\ll f_{a}$), it can be accurately approximated by

(note that now $m_{a}^{*}$ is not constant but a function that varies with $r$).

The case with $1-F_{n}=0$ at some critical radius $r_{c}$ in the star interior has been studied in the literature [26] (see [22, 23, 24] for a non-standard QCD axion model). In such conditions, the axion condenses, and the solution to the KG equation must also be obtained numerically. However, in the limit $m_{a}^{*}R\gg 1$, the solution approaches a step function, where a condensate of amplitude $a=\pi f_{a}$ in the core ($r<r_{c}$) drops exponentially for $r>r_{c}$.

We now report the field profiles and total energy stored in the axion field for the QCD axion model. In Figure 3 we study the QCD axion sourced by NSs for representative values of $f_{a}$. In the left panel we vary $f_{a}$ in the range $f_{a}\in[2\times 10^{16},4\times 10^{17}]$ GeV. We also report the dimensionless quantity $m^{out}_{a}R_{NS}$, where $m^{out}_{a}$ denotes the mass of the QCD axion in vacuo. Contrarily to the ALP case, varying $f_{a}$ affects not only the axion mass but also the couplings (and hence the magnitude of the source term in the KG equation). As a consequence, we find that the QCD axion field at the center of the star is largest for $f_{a}=8\times 10^{16}$ GeV. For higher values of $f_{a}$, the axion-fermion couplings ($\propto f^{-1}_{a}$) decrease, reducing the source term in the KG equation and the magnitude of the axion field. Compared to the ALP field calculated in Section II.2, the QCD axion field is approximately six orders of magnitude smaller. In the right panel of Figure 3 we calculate the total energy stored in the axion field for several values of $f_{a}$. Our results show that the total energy is orders of magnitude smaller than the ALP case, attaining approximately $10^{31}$ erg for $f_{a}=8\times 10^{16}$ GeV.

In general, the smaller the axion field and its gradients, the lower the energy stored in the axion field. This limits the possible effects of the axion field on typical observables of compact stars, as discussed below.