Detectability of Axion Dark Matter with Phonon Polaritons and Magnons

Collective excitations in condensed matter systems, such as phonons and magnons, have recently been proposed as novel detection channels for light dark matter. We show that excitation of i) optical phonon polaritons in polar materials in an ${\mathcal O}$(1 T) magnetic field (via the axion-photon coupling), and ii) gapped magnons in magnetically ordered materials (via the axion wind coupling to the electron spin), can cover the difficult-to-reach ${\mathcal O}$(1-100) meV mass window of QCD axion dark matter with less than a kilogram-year exposure. Finding materials with a large number of optical phonon or magnon modes that can couple to the axion field is crucial, suggesting a program to search for a range of materials with different resonant energies and excitation selection rules; we outline the rules and discuss a few candidate targets, leaving a more exhaustive search for future work. Ongoing development of single photon, phonon and magnon detectors will provide the key for experimentally realizing the ideas presented here.

level axion interactions of interest are: where the three terms are the axion's electromagnetic, wind and electric dipole moment (EDM) couplings, respectively. In the nonrelativistic limit, the effective interaction Hamiltonian is 1 These couplings can be further matched onto axion couplings to low energy degrees of freedom in a crystal. In particular, phonon excitation results from couplings to atomic displacements u lj = x lj −x 0 lj , where l labels the primitive cell, j labels the atoms within each cell, and x 0 lj are the equilibrium positions, while magnons can be excited via couplings to the (effective) spins of magnetic ions S lj .
An axion field oscillating with frequency ω = m a and wavenumber p = m a v a is represented by a(x, t) = a 0 cos (p · x − ωt) , where the field amplitude is related to the energy density via ρ a = m 2 a a 2 0 /2. The resulting effective 1 The coupling to the axial current also generates a term proportional to mas f · v f , we neglect this term since its coupling to collective spin excitations is suppressed compared to the one generated by the ∇a · s f term.
Hamiltonian relevant for phonon and magnon production takes the general form δĤ = δĤ 0 e −iωt + c.c. , with the effective couplings, f j , proportional to a 0 and the relevant axion coupling. While our focus here is axion DM, the same equations hold for general field-like DM candidates.
In Sec. II, we derive rate formulae for single phonon and magnon excitations starting from the general form of couplings in Eq. (4). In the case of phonon excitation, the true energy eigenmodes in a polar crystal, at the low momentum transfers relevant for dark matter absorption, are phonon polaritons due to the mixing between the photon and phonons. We take this mixing into account while still often referring to the gapped polaritons as phonons since their phonon components are much larger. The final results for phonon and magnon excitation rates are Eqs. (18) and (27). Depending on the couplings f j and symmetries of the target system, it often happens that excitation of some of the phonon or magnon modes is suppressed, reducing the sensitivity to DM. We discuss this problem and possible ways to alleviate it in Sec. III.
Then it remains to determine the effective couplings f j in terms of particle physics parametersg aγγ , g af f , etc. -in the case of axion DM. The effective couplings can receive multiple contributions, some of which rely on the presence of an external field. We discuss the various possibilities for axioninduced single phonon or magnon production in Sec. IV. Among them, two are particularly promising: the coupling of the gradient of the axion field to the electron spin, g aee , allows for magnon excitation, while the axion-induced electric field in the presence of an external magnetic field, due to the axionphoton coupling g aγγ , can excite phonon polaritons. These processes are summarized in Table I. We present our numerical results for the projected reach via these processes in Sec. V. We find that, when the axion mass is well-matched to phonon polariton or magnon resonances in the target material, the QCD axion can be easily within reach. The sensitivity is inherently narrow-band for any specific target material, with the axion masses covered limited by the resonance widths. However, combining the reach of a set of judiciously chosen materials with different phonon and magnon frequencies can offer a broader coverage. Finally, we conclude in Sec. VI and discuss future interdisciplinary work needed to better understand and realize the potential of the ideas presented in this work. Eq. (18) Axion → magnon ∇a · s e f j = − i √ 2 g aee (g j − 1) √ ρa me v a Eq. (27) TABLE I. Summary of the potentially detectable channels identified in section IV. The axion field a is given by Eq. (3), ρ a is its energy density, and v a is its velocity. The axion couplings g aγγ and g aee are defined in Eqs. (1) and (2), and given by Eqs. (31) and (32) for the QCD axion. ε ∞ is the high-frequency dielectric constant due to electronic screening, Z * j is the Born effective charge tensor of the ion, and g j is the Landé g-factor. ε here.

II. GENERAL FORMALISM FOR ABSORPTION RATE CALCULATIONS
In this section, we adapt the DM scattering calculations in Refs. [35,37]  should be treated as a classical field. Within the coherence time τ a = (m a v 2 a ) −1 ∼ 10 −7 s (10 meV/m a ), its effect can be modeled as a harmonic perturbation on the target system as in Eq. (4). In this work, we focus on configurations with no external AC electromagnetic fields, so that ω = m a . An AC external field with frequency ω e would generate perturbations with ω = |m a ± ω e |, for which the calculations in this section also apply.
Phonons and magnons arise from quantizing crystal lattice degrees of freedom, displacements u lj and effective spins S lj respectively, which DM can couple to, as mentioned in the Introduction -see Eq. (4). The effective couplings f j depend on the atom/ion types, hence the subscript j. We will keep f j general in this section, and derive their expressions for the case of axion DM in Sec. IV.
We assume the target system is prepared in its ground state |0 at zero temperature. The transition rate from standard time-dependent perturbation theory reads Strictly speaking, since phonons and magnons are unstable particles, the sum over final states f should include multi-particle states resulting from their decays. In practice, however, when ω is close to a phonon/magnon resonance, we can simply smear the delta function to the Breit-Wigner function and sum over single phonon/magnon states: 2 where |ν, k is the single phonon/magnon state on branch ν with momentum k, and γ ν,k is its decay width. Away from a resonance, the lineshape deviates from Breit-Wigner, and depends on the details of phonon/magnon interactions. Since we are interested in sub-eV DM candidates, the momentum deposited is limited to m a v a < ∼ meV -the DM field drives phonon/magnon modes close to the center of the first Brillouin zone (1BZ). Finally, averaging over the DM velocity distribution f (v), we obtain the expected total rate: We take f (v) to be a boosted Maxwell-Boltzmann distribution, with parameters v 0 = 230 km/s, v e = 240 km/s, v esc = 600 km/s. The local axion DM density ρ a , which enters the effective couplings f j (see Table I), is assumed to be 0.3 GeV/cm 3 .
In the following subsections, we derive the rate formulae for single phonon and magnon excitations, respectively. For easy comparison, we present both derivations in as similar ways as possible.

A. Phonon excitations
We begin by calculating the absorption rate from couplings to phonons. The target Hamiltonian results from expanding the potential energy of the crystal around equilibrium positions of atoms: where p lj = m julj , and the force constant matrices, V ll jj , can be calculated from ab initio density functional theory (DFT) methods [40][41][42][43].
2 In deriving Eq. (6) we have assumed the observation time t > ∼ γ −1 ν,k , for which the transition rate Γ is time-independent. We also assume that the line width of the axion, ∆ω ∼ mav 2 a ∼ 10 −8 eV (ma/10 meV) is smaller than excitation linewidth, γ ν,k , which is true as long as γ ν,k is greater than ∼ 10 −6 times the resonance frequency.
To diagonalize the Hamiltonian, we expand the atomic displacements and their conjugate momenta in terms of canonical phonon modes: where ν labels the phonon branch (of which there are 3n for a three-dimensional crystal with n atoms per primitive cell), k labels the phonon momentum within the 1BZ, N is the total number of primitive cells, m j is the mass of the jth atom in the primitive cell, andâ † ν,k ,â ν,k are the phonon creation and annihilation operators. The phonon energies ω ν,k and eigenvectors ν,k,j = * ν,−k,j are obtained by solving the eigensystem of V ll jj , for which we use the open-source code phonopy [44]. The target For a polar crystal, since the ions are electrically charged, some of the phonon modes mix with the photon. This mixing has a negligible impact in most of the 1BZ where k ω, and in particular does not affect the DM scattering calculations in Refs. [32,33,37,38]. However, near the center of the To account for the photon-phonon mixing, we write the total Hamiltonian of electromagnetic fields coupling to the ions in the target crystal, and diagonalize its quadratic part via a Bogoliubov transformation. We explain this procedure in detail in Appendix A. The resulting diagonal Hamiltonian At each k, there are (3n+2) modes, created (annihilated) byâ † ν,k (â ν,k ), which are linear combinations of 3n phonon modes and 2 photon polarizations. Among them, 5 are gapless at k = 0, including 3 acoustic phonons and 2 photon-like polaritons. The number of gapped modes, 3n − 3, is the same as in the phonon-only theory, but their energy spectrum is shifted, {ω ν=(6,...,3n+2),k } = {ω ν=(4,...,3n),k }.
The original phonon modes are linear combinations of the phonon polariton eigenmodes: For DM coupling to the atomic displacements u lj , the perturbing potential is given by Eq. (4) and thereforeδ where |ν , p =â † ν ,p |0 . To arrive at the second equation, we have used the identity l e i(p−k)·x 0 lj = N δ k,p (for k, p ∈ 1BZ). It follows that where we have swapped the dummy indices ν and ν . The DM absorption rate per unit target mass, where ω = m a , p = m a v a , and m cell is the total mass of the atoms in a primitive cell. For our numerical calculations, we use the phonopy code [44] to process DFT output [33,37,38] to obtain the unmixed phonon energies and eigenvectors ω ν ,p , ν ,p,j as mentioned above, and then compute the polariton-corrected energy eigenvalues ω ν,p and mixing matrices U, V via the algorithm of Refs. [45,46].
We relegate the technical details to Appendix A, and review the diagonalization algorithm [45,46] in Appendix C.

B. Magnon excitations
We now move to the case of magnons. The target Hamiltonian is a spin lattice model, with the following general form:Ĥ = ll jj where l, l label the magnetic unit cells and j, j the magnetic ions inside the unit cell, µ B = e 2me is the Bohr magneton, B is an external uniform magnetic field, and g j are the magnetic ions' Landé g-factors.
In the simplest case of the Heisenberg model, J ll jj ∝ 1 for pairs of lj and l j on (nearest, next-tonearest, etc.) neighboring sites. One material which is well described by this simple model is yttrium ion garnet (YIG) [47,48], which has already been considered for DM detection [13-15, 19, 35, 49].
However, as we will see below, materials with spin-spin interactions beyond the simplest Heisenberg type can be useful for enhancing DM-magnon couplings.
The spin-spin interactions in Eq. (19) can result in a ground state with magnetic order. Here we focus on the simplest case of commensurate magnetic dipole orders, for which a rotation on each sublattice can take S lj to a local coordinate system where each spin points in the +ẑ direction: YIG and many other magnetic insulators have commensurate magnetic order. The calculation can be easily generalized to single-Q incommensurate orders, as we discuss in Appendix B.
For a magnetically ordered system, the lowest energy excitations are magnons. To obtain the canonical magnon modes, we first apply the Holstein-Primakoff transformation to write the Hamiltonian in terms of bosonic creation and annihilation operators: where S ± lj = S x lj ± iS y lj . The Holstein-Primakoff transformation ensures that the spin commutation relations [S α lj , S β l j ] = δ ll δ jj i αβγ S γ lj are preserved when the usual canonical commutation relations [â lj ,â † l j ] = δ ll δ jj are imposed. As in the phonon case, translation symmetry instructs us to go to momentum space:â where k ∈ 1BZ. The quadratic Hamiltonian, whose detailed form can be found in Appendix B, only couples modes with the same momentum, i.e.â j,k andâ j ,k ,â † j ,−k . A Bogoliubov transformation takes the quadratic Hamiltonian to the desired diagonal form: At each k, there are n magnon modes with n the number of spins per magnetic unit cell. These energy eigenmodes are created (annihilated) byâ † ν,k (â ν,k ), which are related to the unprimed creation and annihilation operators byâ For DM coupling to the effective spins S lj , the interaction is given by Eq. (4), and we find, in complete analogy with Eq. (16), where We can now obtain the DM absorption rate per unit target mass: where n 0 is the number of gapless modes, which depends on the material (in particular, on the symmetry breaking pattern). Similarity to the phonon formula Eq. (18) is apparent. We again use the algorithm of Refs. [45,46], reviewed in Appendix C, to solve the diagonalization problem to obtain the magnon energies ω ν,p and mixing matrices U, V.

III. SELECTION RULES AND WAYS AROUND THEM
Depending on the DM couplings f j , excitation rates for some of the phonon or magnon modes can be suppressed. In the context of DM scattering, it has been known that acoustic and optical phonons are sensitive to different types of DM couplings [37,38]: if DM couples to the inequivalent atoms/ions with the same sign (different signs), the single phonon excitation rate is dominated by acoustic (optical) phonons, corresponding to in-phase (out-of-phase) oscillations of the atoms/ions.
The same considerations apply to absorption of DM, though here the gapless acoustic phonons are kinematically inaccessible, and therefore only gapped optical phonons can be excited. Thus, the rate is suppressed if all f j point in the same direction. As an extreme example, consider f j = m j f , with f a constant vector. Up to the photon-phonon mixing (which mostly shifts the energy eigenvalues while leaving the factor (U * ν ν,p + V ν ν,−p ) close to δ ν ν ), the rate in Eq. (18) is proportional to j √ m j f · * ν,p,j 2 . However, one can show from translation symmetry that √ m j f can be written as a linear combination of the polarization vectors ν,p→0,j with ν ∈ acoustic -see for example the explicit discussion of [50], consistent with the earlier results in Refs. [32,33]. This is not surprising since gapless acoustic phonons are Goldstone modes of the broken translation symmetries. Thus, by the orthogonality of the phonon polarization vectors, optical phonons do not contribute to the rate in the p → 0 limit, and only higher order terms in the DM velocity can give a nonzero contribution. In the next section, in the context of axion DM, we will encounter both cases where the suppression due to the f j 's being aligned is present and absent, and will identify a process free of the suppression as a viable detection channel.
Additional selection rules may be present among the optical phonons. For example, sapphire has 27 optical phonon branches, but we find that near the 1BZ center, only 10 of them couple to the axion-induced electric field (in the presence of an external magnetic field). Furthermore, 8 of the 10 modes are degenerate in pairs, reducing the total number of distinct resonances to 6, as seen in Fig. 1. This is consistent with the well-known fact that, due to crystalline symmetries, sapphire has 6 infrared-active phonon modes [51][52][53]. Thus, despite the existence of many optical phonon modes, sapphire does not really offer broadband coverage of the axion mass. The same observation, that only a subset of gapped phonon modes couple to axion DM, also holds for the other targets considered: SiO 2 and CaWO 4 -see Fig. 1 (GaAs has only 3 optical phonon modes which are all degenerate and can couple to axion DM). To broaden the mass coverage, it is therefore necessary to run experiments with several target materials with distinct phonon frequencies.
There are also selection rules in the case of magnon excitations. It has been pointed out that, assuming the absence of an external magnetic field, for a target system described by the Heisenberg model with quenched orbital angular momentum, such as YIG, only gapless magnons can be excited in the zero momentum transfer limit [35,54]. To understand why, let us review and quantify the semiclassical argument given in Ref. [35]. Within a coherence length, the DM field couples to the spins as a uniform magnetic field, causing all the spins to precess. As a result, the rate of change in the Heisenberg interaction energy between any pair of spins is proportional to which vanishes for f j = f j . Therefore, if the target system is described by the Heisenberg Hamiltonian, and all f j are equal (which is quite generic since they all originate from DM-electron spin coupling), the total energy cannot change in response to the DM field, and no gapped magnons can be excited. In other words, the DM field only couples to gapless magnons. In the case of scattering, the rate is not severely suppressed by this fact since the scattering kinematics allows access to finite momentum magnons, and hence sufficient energy deposition to be detected.
The situation is much worse in the case of absorption, because the momentum transfer is small in comparison to the DM mass due to its small velocity v ∼ 10 −3 . Therefore, gapless modes cannot be excited due to kinematics. As a result, the detection rate is severely suppressed by powers of DM  Heisenberg interactions, which we include up to third nearest neighbors [48]. The ground state has ferrimagnetic order, where the 12 magnetic ions on the tetrahedral sites and the 8 magnetic ions on the octahedral sites have spins pointing in opposite directions [47], taken to be ±ẑ. The symmetry breaking pattern is SO(3) → SO (2), and hence there are two broken generators S x , S y . There is however just one Goldstone mode (with quadratic dispersion) due to the nonvanishing expectation value of the commutator between the broken generators, [S x , S y ] = (i/2) S z = 0 [55][56][57][58]. Thus, among the 20 magnon branches, only one is gapless. We find that at zero momentum, the j sum in the rate formula Eq. (27), j S j · U * jν,0 r j + V jν,0 r * j , indeed vanishes for all but the gapless mode (ν = 1), confirming the argument above that the DM field only couples to gapless magnons. b. Ba 3 NbFe 3 Si 2 O 14 . This is an example of materials with incommensurate magnetic order. We discuss the generalization needed in the rate calculation in Appendix B, with the final result given in Eq. (B11). The magnetic unit cell contains three magnetic ions Fe 3+ with spin 5/2, which form a triangle in the x-y plane. The crystal consists of layers stacked in the z direction. The antiferromagnetic Heisenberg interactions result in a frustrated order with 120 • between the three spins that are nearest neighbors. Further, the chiral structure of inter-layer Heisenberg exchange couplings results in a rotation of the order in the z direction, with a wavevector that is irrational, Q 0.1429 (2π/c)ẑ where c 5.32Å is the inter-layer lattice spacing [45]. This is known as a single-Q incommensurate order.
All 3 generators of SO(3) are broken while the ground state has zero total magnetization, so there are 3 Goldstone modes. These appear at k = 0, ±Q, which are also the momenta near which the axion coupling is nonzero due to (generalized) momentum conservation. We find that at all three momenta, the j sum in the generalized rate formula Eq. (B11) is nonzero only for ν = 1, i.e. the gapless modes, again confirming the argument above that the DM field only couples to gapless magnons.
Nevertheless, there are several possibilities to alleviate the problem. First, one can consider targets involving additional, non-Heisenberg interactions. These additional terms can explicitly break the rotational symmetries, causing the otherwise gapless Goldstone modes to become gapped, and match the DM absorption kinematics. Concretely, we can identify two ways of implementing this idea: • An external magnetic field B = 0 can generate a gap for the lowest magnon branch equal to the Larmor frequency, assuming g j = 2 for all j. The QUAX experiment [13][14][15] makes use of this to search for axion DM with m a ∼ O(0.1 meV), currently in the regime where the magnon number is large and a classical description can be used. Recently, a calculation to exploit this effect in the quantum regime, similar to the derivation in Sec. II B, was carried out in Ref. [49]. However, sensitivity of such a setup is limited to sub-meV DM by the achievable magnetic field strengths, and only the lowest magnon mode(s) can be excited.
• There are materials with anisotropic interactions where the number of gapless Goldstone modes is reduced. In Sec. V B, we consider a concrete example, NiPS 3 , where the J ll jj matrices in the spin Hamiltonian Eq. (19) have unequal diagonal entries [59]. In this case, two gapped magnon modes at 12 meV and 44 meV can couple to axions. There are also materials with nonzero off-diagonal entries in J ll jj , arising from e.g. Dzyaloshinskii-Moriya interactions, that could be used to achieve the same result.
An orthogonal route to solve the problem is to use targets where the DM-spin couplings f j are nondegenerate.
• For materials with nondegenerate Landé g-factors, even a uniform magnetic field can drive the magnetic ions differently and excite gapped magnon modes; the same is true for a uniform DM field. The basic reason for this is the presence of spin-orbit couplings that break the degeneracy between the DM couplings to the magnetic ions' total effective spins. Concretely, the Landé g-factors are given by where s j and j are respectively the spin and orbital angular momentum components of the total effective spin S lj . In the simplest and most common case of magnetic ions with quenched orbital angular momenta (i.e. j 0), we recover the usual result g j = 2. Breaking the degeneracy requires the magnetic ions to have different spin and orbital angular momentum compositions.
We demonstrate how this allows for axion couplings to gapped magnons in Sec. V B.
As in the case of phonons, there are usually additional selection rules due to crystalline symmetries.
As a result, the strategies discussed above usually generate axion couplings to only a subset of gapped magnon modes -see Fig. 1. Therefore, multiple target materials which cover complementary ranges of magnon frequencies are desirable.

IV. AXION COUPLINGS AND DETECTION CHANNELS
The derivation and discussion in the previous two sections apply to general field-like DM candidates.
We now specialize to the case of axion DM. Our goal in the present section is to identify the most promising detection channels involving phonon or magnon excitation via order-of-magnitude estimates.
We then examine these processes quantitatively in the next section.
The axion couplings of interest are already given in Eqs. (1) and (2). For the QCD axion, we have [60] g aγγ = C γ α 2πf a = 2.03 × 10 −12 C γ m a 10 meV g anγ = −g apγ = (3.7 ± 1.5) × 10 −3 1 f a 1 GeV = (6.5 ± 2.6) × 10 −12 m a 10 meV where we have denoted C γ ≡ E/N − 1.92 (4)  In the following subsections we consider axion couplings independent of external fields and in the presence of a magnetic field. 3 In each case, we discuss the phonon and magnon excitation processes that are allowed, and identify those with potentially detectable rates. The results of this exercise are summarized in Table I.

A. Axion couplings independent of external fields
The axion wind coupling to electron spin leads to a coupling to the spin component of S lj . From s lj + lj = S lj and 2s lj + lj = g j S lj , we see that the axion wind couples to s lj = (g j − 1) S lj . Thus, In the notation of Eq (4), we thus have For an order of magnitude estimate of the rate, let us note that the mixing matrices U, V in Eq. (27) generically scale as 1/ √ n with n the number of magnetic ions in a primitive cell. The maximum rate is obtained on resonance, which is parametrically given by where n s and ρ T are the spin and mass densities of the target, taken to be (5Å) −3 and 5 g/cm 3 , respectively (close to the values for YIG), in the estimate. We see that, with single magnon sensitivity, interesting values of g aee may be reached with less than a kilogram-year exposure.
The axion wind also couples to nucleon spins. However, these couplings do not excite magnons, since magnetic order originates from electron-electron interactions, meaning that the effective spins of However, the rate suffers from multiple suppressions. First, coupling to atomic displacements relies on the spatial variation of ∇a · S N , which brings in an additional factor of v a on top of the gradient: Second, there is a further suppression for exciting optical phonons since f j are approximately aligned with acoustic phonon polarizations if all S N,j point in the same direction (see discussion in Sec. III). Even without taking into account the second suppression, the estimated on-resonance rate using Eq. (18), can be sizable only for uninterestingly low f a . Therefore, we conclude that axion wind couplings to nucleon spins do not offer a viable detection channel.
Beyond the axion wind couplings, the axion field also turns magnetic dipole moments from s f into oscillating EDMs, and one may consider phonon and magnon excitation by the resulting electromagnetic radiation fields. In the case of the electron, since the EDM coupling is perturbatively generated by the aFF coupling, the process mentioned above is essentially converting the crystal magnetic field into electromagnetic radiation, and should be less efficient than applying an external magnetic field (discussed below in Sec. IV B). In the case of nucleons, the EDM coupling in Eq. (1) is not much larger than the perturbative contribution from aFF , so the same conclusion applies.
In sum, for axion couplings independent of external fields, we have identified magnon excitation via the axion wind coupling to electrons as the only viable detection channel.

B. Axion couplings in a magnetic field
In the presence of a DC magnetic field B, the axion field induces oscillating electromagnetic fields via the aFF coupling. Solving the modified Maxwell equations (see e.g. Ref. [61,62]), we find the induced electric field is E a = −g aγγ a ε −1 ∞ · B. Note that the high frequency dielectric constant, ε ∞ , which takes into account screening effects from the fast-responding electrons while excluding ionic contributions, should be used when solving the macroscopic Maxwell equations -essentially, since we are concerned with how the axion induced electric field acts on the ions, the ion charges should appear in the source term rather than being coarse grained into a macroscopic electric polarization. In the long-wavelength limit, the axion-induced electric field E a couples to charged ions via an effective dipole coupling: where Z * j is the Born effective charge tensor of the jth ion in the primitive cell -it captures the change in macroscopic polarization due to a lattice displacement, δP = e Z * j · δu lj /Ω, and is numerically close to the ionic charge, Z * j Q j 1. It follows that Noting that the phonon polarization vectors scale as 1/ √ n with n the number of ions in the primitive cell, and assuming photon-phonon mixing gives just a small correction, we can estimate the onresonance rate from Eq. (18) as follows: where we have taken Z * /ε ∞ ∼ 1 and m ion ∼ 20 GeV. We see that, with single phonon sensitivity, there is excellent potential for reaching the QCD axion coupling if the axion mass is close to a phonon resonance.
The axion-induced magnetic field is much smaller, B a ∼ E a v a . An order of magnitude estimate tells us that the magnon excitation rate by B a is much smaller than the phonon excitation rate by In sum, we have identified phonon excitation via the axion-photon coupling in a magnetic field as the only novel viable detection channel when considering an external DC magnetic field.

V. PROJECTED SENSITIVITY
We now compute the projected sensitivity for the two detection channels identified in the previous section (see Table I). In both phonon and magnon calculations, an important but elusive parameter is the resonance width, γ ν,p , of each mode. While all other material parameters entering the rate calculation (equilibrium positions, phonon energies and eigenvectors, magnon energies and mixing matrices) can be computed within the quadratic Hamiltonian, the widths involve anharmonic interactions and are not always readily available. In what follows, we present the reach with reasonable assumptions for γ ν,p . Our goal here is to demonstrate the viability of phonons and magnons for detecting axion DM, and motivate further study in the condensed matter and materials science community on phonon and magnon interactions. This will be crucial both for having more accurate inputs to the DM detection rate calculation, and for designing detectors to read out these excitations.

A. Phonon excitation via the axion-photon coupling
When axion absorption excites a phonon polariton, it cascade decays to a collection of lower energy phonons and photons on a timescale of γ −1 ∼ ps (meV/γ). Since the polariton is a phonon-like state which decays via anharmonic phonon couplings, theoretically it is most efficient to read out phonons (heat) in the final state. However, phonon readout (e.g. through a transistor edge sensor) is complicated by the strong external magnetic field needed for the axion absorption process. One possibility is to detect the phonon via evaporation of helium atoms deposited on the surface of the crystal. The helium atoms are then detected well away from the crystal, such that the magnetic field is isolated from the sensor. 4 Alternatively, if a target material can be found in which the photon yield from phonon polariton decays is substantial, photon readout becomes a viable option. In this case, the photons may be focused by a mirror and lens and directed by a waveguide onto a single photon detector (e.g. a superconducting nanowire) placed in a region of zero magnetic field. Both the phonon and photon readout possibilities sketched above will be studied in future work. In what follows, we simply assume 3 single phonon polariton excitation events per kilogram-year exposure (corresponding to 90% C.L. for a background-free counting experiment) when presenting the reach.
To compute the projected sensitivity to the axion-photon coupling, we use the rate formula Eq. for a sapphire target, whenb is parallel (perpendicular) to the crystal c-axis, chosen to coincide with the z-axis here, only 2 (4) out of the 6 resonances appear. This observation provides a useful handle to confirm a discovery by running the same experiment with the magnetic field applied in different directions.

B. Magnon excitation via the axion wind coupling
To compute the magnon excitation rate, we substitute the coupling f j in Eq. (35), into the rate formula Eq. (27). In Sec. III, we discussed three strategies to alleviate the suppression of axion-magnon couplings due to selection rules: external magnetic fields, anisotropic interactions, and nondegenerate g-factors. In this subsection, we show the projected reach for each of these strategies. The results are summarized in Fig. 5, assuming 3 single magnon events per kilogram-year exposure. Absent a detailed study of anharmonic magnon interactions, we take the resonance widths to be a free parameter, and show results for γ/ω = 10 −2 and 10 −5 , consistent with measured phonon width values on the high end and YIG's Kittel magnon width on the low end. We see that, on resonance, all methods could reach axion-electron couplings predicted by QCD axion models. In the following, we expand on the calculation for each strategy.
a. External magnetic field. The idea of using an external magnetic field to lift the gapless mode is the one adopted in the QUAX experiment [13][14][15]. In Ref. [14], a classical calculation was used to estimate the axion absorption rate. Our formalism allows to compute the same rate in the quantum regime, and agrees with the recent computation carried out in Ref. [49]. The projected reach we obtain for a YIG sample in a 1 T field is shown in Fig. 5, where the resonance is at 0.12 meV (see Eq. (29)).
For comparison, we also overlay the projection in Ref. [49] based on scanning the resonance frequency by changing the magnetic field, with a total integration time of 10 years.
b. Anisotropic interactions. As another way to lift the gapless magnon modes, so that the absorption kinematics are satisfied, we may use materials with anisotropic exchange couplings. As an example, we consider NiPS 3 , which has a layered crystal structure [65]. The magnetic ions are spin-1 Ni 2+ . Following Ref. [59], we  lifting gapless magnon modes by an external magnetic field (YIG target in a 1 T magnetic field, compared to the scanning scheme of Ref. [49]), anisotropic interactions (NiPS 3 target), and using targets with nondegenerate g-factors (hypothetical toy models based on YIG, referred to as YIG o and YIG t ). For all the cases considered we assume 3 events per kilogram-year exposure, and take the magnon width to frequency ratio γ/ω to be 10 −2 (solid) or 10 −5 (dashed).
c. Nondegenerate g-factors. Finally, we consider coupling the axion to gapped magnon modes in the presence of nondegenerate g-factors. We are not aware of a well-characterized material with nondegenerate g-factors so, as a proof of principle, we entertain a few toy models, where a nondegenerate component is added to the effective spins S in YIG. In reality, all the magnetic ions Fe 3+ in YIG have ( , s, S) = (0, 5/2, 5/2); the orbital angular momenta of 3d electrons are quenched. In Fig. 5, we show the reach for two toy models, with either the octahedral sites or the tetrahedral sites modified to have ( , s, S) = (1, 5/2, 7/2). In each case, only one of the 19 gapped magnon modes, at 7 meV and 76 meV respectively, is found to contribute to axion absorption. This is because, to preserve the lattice symmetries, we have modified all the effective spin compositions on tetrahedral or octahedral sites in the same way. The sole gapped mode that couples to axion DM corresponds to out-of-phase precessions of the tetrahedral and octahedral spins.

VI. CONCLUSIONS
In this paper we showed multiple ways collective excitations in crystal targets are, in principle, sensitive to QCD axion DM. Specifically, we identified two novel detection possibilities: axion-induced phonon polariton excitation in an external magnetic field, and axion-induced magnon excitation in the absence of external fields. These paths are complementary as each probes a different axion coupling, g aγγ and g aee in the phonon case and magnon case respectively.
In the phonon polariton case, we considered several example targets -Al 2 O 3 , CaWO 4 , GaAs and SiO 2 -and showed that on resonance, per kilogram-year exposure, they can reach g aγγ ∼ 10 −12 GeV −1 , as shown in Figs. 3 and 4. This outperforms the leading constraint in this mass window from stellar cooling, and reaches below the QCD axion band. Carefully choosing a set of target materials with different phonon frequencies is key to covering a broad range of axion masses.
Previous proposals for probing g aee via absorption onto magnons, which underlies the QUAX experiment, considered targets in an external magnetic field, which would lift the lowest magnon mode and kinematically allow sub-meV axion absorption [13-15, 19, 49]. In contrast, we focused on the O(1-100) meV axion mass window, and showed that without an external magnetic field, materials with anisotropic exchange interactions, e.g. NiPS 3 , and materials with nondegenerate g-factors can host gapped magnons coupling to the axion. On resonance and with kilogram-year exposure, they can have sensitivity to the DFSZ model and down to g aee ∼ 10 −15 , shown in Fig. 5.
Realizing the exciting potential of discovering O(1-100) meV axion DM via phonons and magnons hinges upon the ongoing effort to achieve low threshold single quanta detection. One possibility is to evaporate helium atoms from phonon interactions at the surface of the crystal and then detect the evaporated helium atoms in a region separated from the magnetic field region; R&D is underway for this direction. Another route is to read out photons produced from the decay of a phonon polariton.
Single photon detectors (e.g. superconducting nanowires) may operate in a field-free region, away from the crystal target and connected to it via a waveguide. Finally, magnons are read out in a resonant cavity in the QUAX setup in the classical regime [13][14][15]. Work is underway to detect single magnons in YIG by coupling cavity modes to a superconducting qubit [66], though as in other resonant cavity searches, axion masses are best produced near the inverse cavity size; for larger axion masses the virtual cavity modes will be off-shell and readout efficiency is suppressed. On the materials side, we would like to make more accurate predictions for the detection rates via an improved understanding of phonon and magnon resonance lineshapes, and explore the possibility of scanning the resonance frequencies by engineering material properties, in order to fully exploit the discovery potential of an axion DM search experiment based on phonon and magnon excitations. As outlined in Sec. II A, mixing between long-wavelength photon and phonon states needs to be taken into account when computing DM absorption rates in a polar crystal. Our starting point is the Lagrangian for an ionic lattice coupling to electromagnetism: In the long wavelength limit, the leading electromagnetic coupling is via the electric dipole, as shown in the third term. Plugging in E = −∇A 0 −Ȧ and integrate by parts, we can write it in the familiar when expanded to linear order in u. The last term in Eq. (A1) results from integrating out electron response. As explained in detail in Ref. [37], the photon self-energy Π µν can be related to the dielectric tensor of the medium, ε (taken to be the electronic contribution, usually denoted by ε ∞ , in the present case), and the photon Lagrangian can be written as It is convenient to choose a generalized Coulomb gauge, ∇ · ε · A = 0, and since the A 0 field is non-dynamical, it can be immediately integrated out. We thus obtain To derive the Hamiltonian, we note the canonical momenta are: Therefore, where Note that while we have written A and P as 3-vectors, they are implicitly assumed to satisfy the gauge condition. This means that, in Eq. (A9), ε should be projected onto the subspace satisfying the gauge condition before the inverse is taken.
We now consider the four terms in turn. First, as mentioned in Sec. II A, we use the phonopy code with DFT calculations of the force constants V (2) lj,l j to diagonalize the lattice (phonon) Hamiltonian H ph , giving Second, the Coulomb term, H Coulomb , becomes more transparent when written as a momentum space integral: with the charge densityρ Expanding u lj in (â ν,k +â † ν,−k ) as in Eq. (11) and summing over l picks out the k = k modes. With the momentum integral discretized, where As a technical note, while there is an option in phonopy (non-analytic correction) to also include H Coulomb in the diagonalization calculation, it seems to work only at k > ∼ ω. Therefore, in our calculation, we use phonopy to diagonalize only H ph , and include H Coulomb separately.
Next, we also write the photon Hamiltonian H EM in momentum space: where We decompose the photon fieldÃ into two orthogonal linear polarizations e 1,k ⊥ e 2,k which satisfy the gauge condition, k · ε · e λ,k = 0 (λ = 1, 2). We choose the basis in which the projection of ε onto the two-dimensional subspace, e λ · ε · e λ , is diagonal, with eigenvalues ε 1 , ε 2 . Denote the projection of K 2 in this basis by e * λ ,k · K 2 · e λ,k ≡ K 2 λ λ .
Finally, the photon-phonon mixing term H mix can be written in terms of the creation and annihilation operators according to Eqs. (12) and (A19): The matrix h k can be written as: where A k is a 3n × 3n matrix given by: √ ω ν,k ω ν ,k (k · ξ * ν,k )(k · ξ ν ,k ) k · ε · k , while B k is a 3n × 2 with the following structure and C k,νν is a 2 × 2 matrix given by We review a general algorithm for diagonalizing such Hamiltonians in Appendix C.

(B4)
The matrix h k can be written in terms of n × n sub-matrices as: where, by defining J lj,l j = R (x lj ) T · J lj,l j · R (x l j ), we have The procedure to diagonalize Hamiltonian of this kind is reviewed in Appendix C.
By using the commutation relations in Eq. (C2), it can be easily shown that, up to constant terms, Eq. (C3) is equivalent to Eq. (14) and Eq. (23) with the U and V matrices implicitly defined as Such diagonalization procedure (usually called para-unitary diagonalization) can be achieved if h(k) is positive definite. If the spectrum contains zero energy modes, the h(k) matrix will be positive semidefinite and such para-unitary diagonalization may not exist. This problem can be cured by adding a small value to the diagonal components of h(k) [46]. This introduces a small and negligible shift in the spectrum but makes h(k) positive definite and allows for a para-unitary diagonalization.
Ref. [46] provides a simple three-step algorithm to find the linear transformation T and the associated eigenvalues: • A Cholesky decomposition is applied to find a complex matrix K k such that h k = K † k K k .
• The eigenvalue problem for the Hermitian matrix K k gK † k is solved, and the resulting eigenvalues used to form the columns of the matrix U k . The order of the columns is chosen such that the first N elements of the diagonalized matrix L = U † K k gK † k U are positive and the last N negative.
• Finally, the matrix E k in Eq. (C4) is simply related to L k by E k = gL k , and the T k matrix is