Axion-matter coupling in multiferroics

Multiferroics (MFs) are materials with two or more ferroic orders, like spontaneous ferroelectric and ferromagnetic polarizations. Such materials can exhibit a magnetoelectric effect whereby magnetic and ferroelectric polarizations couple linearly, reminiscent of, but not identical to the electromagnetic $\boldsymbol{E}\cdot \boldsymbol{B}$ axion coupling. Here we point out a possible mechanism in which an external dark matter axion field couples linearly to ferroic orders in these materials without external applied fields. We find the magnetic response to be linear in the axion-electron coupling. At temperatures close to the ferromagnetic transition fluctuations can lead to an enhancement of the axion-induced magnetic response. Relevant material candidates such as the Lu-Sc hexaferrite family are discussed.


I. INTRODUCTION
Dark matter (DM) is believed to comprise over three quarters of the universe's mass-density, yet so far it eludes detection in spite of the intensive search using numerous detection schemes [1]. One possible explanation is that the DM particles are too light to strongly scatter off nuclei and electron-based detection schemes currently proposed, motivating the need for new ideas to explore the sub-GeV range of DM masses. In this regard, the O(meV) axion, a pseudoscalar boson introduced to solve the strong charge-parity problem in quantum chromodynamics [2][3][4][5], offers a particularly well-motivated DM candidate [6][7][8][9][10].
A new avenue to DM detection is offered by quantum materials [11] in which the energy scales of the excitations coincide with the requirements for lower-mass DM detection. Moreover, quantum materials possess entangled collective degrees of freedom that allow for an enhanced coupling to axions. Traditionally, the majority of axion detection proposals have focused on exploiting the axion-photon coupling in cavities or plasma haloscopes using strong external fields [12][13][14][15][16][17][18][19][20][21][22]. More recently, axion physics has been increasingly discussed in the condensed matter context [23][24][25][26] primarily in analogy to axion electrodynamics in topological insulators [27][28][29][30][31], but also in terms of particle axion detection via resonant coupling to quasiparticles [32,33]. In most of the existing proposals, the detection scheme uses DM particle absorption to induce single particle excitations in the sensor material, e.g., the particle-hole excitations in a small gap material.
FIG. 1. An incident axion a couples to multiferroic orders via the electron spin and momentum below the Curie temperature. Close to the onset of ferromagnetism, softness in the magnetization makes the system susceptible to the axion coupling. The magnetic material response is imagined sensed by a sensitive magnetometer such as a superconducting quantum interference device (SQUID). The inset shows a sketch of the A2 phase of multiferroic Lu1−xScxFeO3 (h-LSFO) [34].
An alternative approach is to detect DM by their coupling to collective modes of the quantum material. With the rapid progression of the field of axion electrodynamics, we anticipate more such proposals for axion sensing schemes will be explored.
Here we explore multiferroics (MFs) as a platform for DM detection using materials with parallel ferroic orders P M as an "impedance matched" medium to interact with axion DM. Here, P (M ) is the electric (magnetic) polarization vector. Multiferroics with P M act as the "condensate" axion field seen by an arriving external axion. The external axion field linearly couples to the polarization fields and both changes net magnetization and excites the collective modes of the MF. The mechanism for DM detection proposed here is qualitatively different from existing schemes in the sense that the axion DM affects the magnetization-a macroscopic observable-by modifying the parameters that control the macroscopic state. As such, the proposed method of detection is similar to the transition-edge sensors, except that our sensor (the multiferroic sample) remains in a macroscopically ordered state before and after interaction with DM.
To place the current discussion in a broader context: One finds two main developing threads for axion DM sensors. One approach is relying on the sensor device being placed in external electric (E), magnetic (B) or both fields simultaneously, whereby one breaks time reversal (T ) and parity (P) to induce coupling to the axion field. Another approach is to use systems that take advantage of spontaneous symmetry breaking (no external fields applied). In the phase with broken P and T symmetry coupling of the axion DM field to matter fields occurs naturally. The coupling induces collective excitations that are subsequently read out by appropriate sensing means. Here we propose a new realization of the latter platform where P and T breaking (though with P T being a symmetry operation), and hence axion electrodynamics, occurs naturally.
Our axion-matter sensing proposal is sketched in Fig. 1. Our central scheme is that in MFs with parallel ferroic orders in zero external fields one has a medium that linearly couples to an axion field a through an aP · M term in the Hamiltonian. This term allows for a novel axion DM detection scheme: Impinging axions will hybridize and directly couple with ferroic orders in MFs without requiring an intermediate photon coupling. This coupling between the axion field and matter fields (with P M ) is a bulk effect. The electrons in the MF produce a P and T breaking term that can be expressed in terms of P and M but which microscopically results in the bulk expectation value of the axion-electron coupling as described in Sec. III A. The macroscopically ordered state provides a set of sensors where each "atomic" unit interacts with the same axion field, and their contributions add coherently-the signal scales with a macroscopic number of atoms. Phrased differently, we obtain Avogadro scaling due to the P M coherent state of the material. The proposed scheme thus offers a quantum advantage at the macroscopic level. This "Avogadro advantage" is a unique feature of macroscopically coherent quantum matter where each unit cell, or each part of the coherent condensate, adds to increased detection sensitivity. We see this advantage operational for the numerous detection schemes relying on the macroscop-ically coherent states ranging from ferromagnets to superfluid/superconducting condensates and MFs. To discuss multiferroics as a proposed platform in this context we will use the acronym MERMAID (multiferroic matter axion detector).
We also present a known material possessing our required coupling to give materials-specific estimates of our MERMAID scheme. We focus on the recently discovered hexagonal LuFeO 3 thin films and the Lu 1−x Sc x FeO 3 crystal family, which have a ferromagnetic ground state for x ≈ 0.4 below T C M ≈ 160 K [34,35]. Sample candidates may be synthesized and characterized using THz magneto-optical spectroscopy to determine the nature of the collective excitations, including electromagnons, and the resulting P M coupling strength.
This paper is organized as follows. We provide background on axions in Sec. II. An effective theory of the axion-matter coupling in multiferroics is presented in Sec. III. In Sec. III A we discuss the effective axionelectron coupling in multiferroics, and we pursue an effective Ginzburg-Landau theory to describe the dynamical response in Sec. III B. Critical slowing down and the axion-induced magnetic response are discussed in Sec. III C and III D, respectively. In Sec. III E we discuss the advantages of Avogadro scaling in multiferroics. In Sec. IV we discuss material candidates and estimates from ab initio density functional theory calculations for a specific compound. We provide conclusions and an outlook in Sec. V.

II. BACKGROUND: AXIONS
The axion is an hypothesized pseudoscalar boson originally proposed to solve the strong charge-parity problem via the Peccei-Quinn mechanism in quantum chromodynamics (QCD) [2][3][4]36]. It was realized that the axion could be a viable DM candidate, which has significantly fuelled interest in it [6][7][8][9]. The axion couples to electromagnetism (in units where c = = 1) via where F µν is the electromagnetic tensor withF µν = 1 2 ε µναβ F αβ its dual (ε µναβ is the Levi-Civita symbol), and where a (m a ) is the axion field (mass), and g γγa is a dimensionful coupling strength.
Although axion-like particles with a Lagrangian like Eqs. (1) and (2) within a vast mass range could account for a partial DM abundance, QCD axions have the particular coupling strength g γγa = C a ma GeV 2 , where C a is a model dependent number that ranges between −0.39 and 0.15 for various theories [37]. Despite an in principle still vast (m a , g γγa ) parameter space, the axion mass is restricted by m a 0.1 eV from astrophysical bounds [14] and by m a 10 −12 eV from cosmological bounds [6][7][8]. Moreover, the coupling strength is roughly limited to |g γγa | < 10 −10 GeV −1 by helioscope experiments [14]. Interest in the mass range around meVs has been fuelled by extensive lattice QCD computations [38] and searches at the XENON experiment [39].
Moreover the axion couples to Dirac fermions via [40] where ψ is the Dirac spinor and g af is a dimensionless coupling strength. The coupling to electrons is limited by |g ae | < 3 × 10 −11 from solar neutrino experiments [41,42]. The relevant low-energy Hamiltonian resulting from Eq. (4) is derived in Appendix A. We note that this Hamiltonian contains an effective term of the form θP · M , which is present even in the absence of external fields (A = 0), see Eq. (A4).
Since the axion is expected to be light, we can treat a as a classical field with a long de Broglie wavelength compared to the typical experimental probe length scale to give an estimate of θ in Eq. (3) for QCD axions. We assume that the the axion field is oscillating with a frequency set by its mass, a = a 0 exp(im a t). Assuming that the axion makes up the local DM density [14,43,44], ρ DM ≈ 300 MeVcm −3 , the amplitude a 0 is fixed by ρ DM = 1 2 m 2 a |a 0 | 2 , yielding e.g., |θ| ≈ |a 0 g γγa | ∼ 10 −22 .

III. EFFECTIVE THEORY
In this section we develop a framework to study the magnetic response in MFs mediated by an axion-fermion coupling. In Sec. III A we discuss how an effective matter analogue of the axion E · B coupling emerges in MF systems, with details of the derivation listed in Appendices A, B, and C.

A. Multiferroics and effective matter coupling
With a matter analogue of Eq. (2) in mind, multiferroics are materials in which spontaneous electric and magnetic polarizations can develop [34,[45][46][47][48][49]. In solidstate crystals, ferroelectricity is normally driven by lone-pair, geometric, charge ordering, or spin-driven mechanisms [47], which result in a macroscopic breaking of inversion symmetry. As a result of spontaneous symmetry breaking in the MF bulk of both time reversal and inversion symmetry, there exist internal electric and magnetic fields (even in the absence of external fields) resulting from these symmetry-broken order parameters. Based on symmetry considerations one should expect an effective χθP · M coupling to be present in the MF phase, despite this term being commonly overlooked in the high-energy literature.
Microscopically, an effective P · M axion term can be derived from the coupling to electrons, as per Eq. (4), in the MF ground state. For details we refer to Appendices A and B. The key result is that the axion-electron coupling at low energy contains the term where |Ψ is the MF ground state. In the MF phase Eq. (5) is dictated by the distortions of electronic orbits (densities) such that one obtains a spontaneous magnetic and electric dipole moment throughout the whole sample. In a (rare) subclass of MFs these moments are parallel, as discussed in Sec. IV. Despite the m a /m e suppression, coming from the µ = 0 contribution to Eq. (4), we argue in Sec. III E that the energy of Eq. (5) can add coherently over a macroscopic volume limited by the MF domain size, which can result in an appreciable effect. To be clear, the interaction in Eq. (5), evaluated on electronic states, yields a bulk contribution that is evaluated by model calculations (Sec. III E) and ab initio means (Sec. IV). In Appendix A we discuss the relation to other couplings that are considered in existing detection schemes [26,50,51].

B. Ginzburg-Landau theory
Using an isotropic sample bulk model, we proceed with an effective description of the longitudinal dynamics of the internal MF degrees of freedom. This approach is valid close to the transition temperature where the free energy can be expanded in powers of the respective order parameters. In terms of the free energy density F , , in the absence of external electric and magnetic fields, we have (in SI units) where ε 0 (µ 0 ) is the vacuum permittivity (permeability), and where α's, β's, γ's, and λ's are phenomenological co-efficients, α ij is the material dependent ME tensor, and χ is a material dependent susceptibility. The transition temperature T C M (T C P ) denotes the magnetic (electric) Curie temperature above which the ferromagnetic (ferroelectric) order melts. In analogy to the P T -odd axion field, where P (T ) is the inversion (time reversal) symmetry operator, the ME tensor is P T -odd if the system (crystal) otherwise respects inversion and time reversal symmetry. The magnitude of the elements α ij are limited by thermodynamic stability of the theory [52]. We consider the case T ≈ T C M T C P (as applicable to the majority of MF materials). This makes P act as a "hard" and M as a "soft" field; for the latter fluctuations are greater. Hence we take the electric polarization to be close to its minimum energy value, and consider fluctuations, where δP are the ferroelectric fluctuations, P = P 0 +δP , and where m P is the effective mass. Below we relabel δP → P and bear in mind that P is a fluctuation field. When associating F in Eq. (6) with a Lagrangian density it yields the classical equations of motion for M and P : where c is the speed of light and α is the ME tensor with elements α ij . Turning off both α and θ reduces Eqs. (7) and (8) to two uncoupled Klein-Gordon equations when ignoring the cubic (Gross-Pitaevskii) term, which is justified when T > T C M . We emphasize that the above effective theory for the magnetic order is a classical and minimal description of the longitudinal magnon dynamics. In situations where the transverse modes are involved, dynamics from the Landau-Lifshitz-Gilbert equation is expected to play a role [53]. Finally, we note that in a crystal the Ginzburg-Landau theory away from low filling should respect symmetries of the crystal point group and partition the solution into irreducible representations. This layer of complexity is omitted in the isotropic bulk theory considered here and would constitute a natural extension of the theory.

C. Critical slowing down
On general grounds we expect to have two regimes of the coupling following the above equations. The axion field provides a periodic perturbation set by the time scale τ a ∼ 1/ω a , while the characteristic time of magnetization is set by the mass term τ ∼ |T − T C M | −1/2 , and more generally, the magnetization dynamics experiences a critical slowing down near T C M : for some dynamical critical exponent zν ∼ 1 determined by the nature of the transition [54][55][56]. The (i) adiabatic regime of the axion-magnetization dynamics occurs at τ a > τ . In this case the effective axion field seen by magnetization has a similar effect to a steady external field, and the response accumulates in time. In the opposite case of (ii) antiadiabatic evolution, τ a < τ , the magnetization "sees" a rapidly oscillating θ(t) term that effectively averages to zero. The most direct and strongest effects of the axion field will be seen in the temperature regime where the axion field dynamics is slow compared to the magnetization time, i.e., regime (i). This requirement places a constraint on how close one can approach criticality, as illustrated in Fig. 2.
To give a numerical estimate of this constraint we assume that τ = τ 0 (1 − T /T c ) −νz for T < T c , where τ 0 is material dependent. The crossover between the adiabatic and the antiadiabatic regime occurs at T /T c = 1 − (τ 0 ω a ) 1/(νz) , given that τ 0 ω a < 1. For m a = 1 meV this would require τ 0 0.5 ps. Experiments conducted on ultrathin Fe films and multiferroic MnWO 4 have found τ 0 in the range of τ 0 ∼ 1 ns [57,58], which if applied to our context would limit the crossover to occur for light axions m a 1 µeV. For instance, for m a = 0.1 µeV and νz = 1 we find T /T c ≈ 0.85. For these estimates, which are merely indicative, one would thus want to maintain the system close to the temperature T = 0.85 T c . In subsequent work we plan to undertake a more systematic analysis of relaxation times.

D. Magnetic response
The coupled dynamical equations (7) and (8) can be treated with linear response theory with θ = θ 0 exp(iω a t) acting as a weak time-dependent driving force. Here we consider the result of such a treatment in one spatial dimension with α ij = αδ ij assumed for simplicity, with the details listed in Appendix D. In short: We find that the axion acts as a periodic driving term to the dynamic ma-terial magnetization. The basic experimental signature of the axion would thus be an oscillatory magnetization of frequency ω a on top of the static (ferromagnetic) background.
When T < T C M T C P and the magnetization is soft compared to the electric polarization, the classical magnetic response resembles that of a driven harmonic oscillator. In frequency space the magnetization has a linear response function of the form with M = M 0 + δM , P ≈ P 0 , where M 0 and P 0 are the static polarizations, and δM are magnetic fluctuations. ω M is the magnetic (magnon) frequency, and η M is a damping coefficient of a term η M ∂ t M added to simulate energy dissipation in the system. Addressing how low η M can be pushed and whether it is a limiting factor experimentally is outside the scope of this paper. However, for spin waves to be resolved we require The expression Eq. (10) sheds light on the key factors in the response. First, there is the overall coefficient, i.e. the factor inside the absolute value. This shows that there are two contributions to the driving: the parallel magnetoelectric effect α, and the static electric polarization P 0 . Since the P 0 term can reasonably be expected to dominate in many cases [59], the ferroelectricity will be a controlling factor in the driving.
Secondly, Eq. (10) has the characteristic frequency dependence of a driven oscillator. Considering the response as a function of axion frequency, there are three regimes, each of which can be clearly seen in Fig. 3: 1. ω a ω M : The antiadiabatic regime discussed above. The response dies off quickly as the axion frequency becomes much larger than the magnon frequency, so this limit is not optimal.
2. ω a ω M : In this regime the response is controlled by the magnon resonance. The response will go as the size of the peak, which will be controlled by the damping coefficient η M .
3. ω a ω M : The adiabatic regime. The axion oscillates much more slowly than the magnetic response of the material, therefore it can be approximated as a static external field. Setting ω a = 0 in Eq. (10), the left hand side behaves as ω −2 M . The response will increase as the temperature approaches the critical magnetic temperature because ω M becomes small there.
As discussed in III C, avoiding the first regime introduces a limitation on what range of axion frequencies can be spanned for a given material and temperature (approximately ω a 0.1 eV). The axion is expected to be a long-wavelength field, and so the relevant magnon frequency is evaluated at the Brillouin zone centre, ω M = ω M (Γ) [60,61].
In the second regime, damping will be a critical factor. The third regime is the one most likely to be accessible in experiment, and the one on which the authors would like to concentrate. It has an attractive feature: The temperature can be used as a tuning dial to adjust both the range of accessible axion frequencies and the size of the response.
In Fig. 3 the frequency response Eq. (10) is shown for parameters motivated by h-LSFO. Both the static response enhancement near criticality and the resonance variation with temperature are clearly visible.
Focusing on the static regime, we would like to study the behaviour of the magnetic response near the critical temperature. As already mentioned, the static response is controlled by ω −2 M , which we must therefore write in terms of the Ginzburg-Landau parameters in Eqs. (7) and (8). To proceed we need the full non-linear solution for the magnetization as presented in Appendix E-here we summarize the key points. The static magnetization satisfies the cubic equation The solution can generically be expanded to linear order in θ as where f (T ), which is just the static limit of the linear magnetic response, depends on λ M , m P , α, P 0 , and γ M . There is one extra subtlety: owing to restored time reversal symmetry in the paramagnetic phase, α = 0 for T > T C M [62,63], and so α itself has a temperature dependence near criticality. We model this using This extra scaling makes the T → T C M limit somewhat subtle, but we find that f (T ) ∼ (T C M − T ) −2b/3 . Therefore the magnetic response grows in principal without bound as the critical temperature is approached. In practice, disorder will pose a challenge for temperature tuning, which motivates consideration of a broader range of magnetoelectric materials for dark matter detection.

E. Avogadro scaling
In "quantum sensing" devices one exploits quantum properties such as entanglement to enhance and measure some physical property [64][65][66][67]. In the case we consider the ground state of the sensor is a coherent MF state with a wave function given by |Ψ = |P x, M x in each unit cell (more precisely, each operational unit of the coherent state). The coherent state of the sensor occurs as a result of spontaneous symmetry breaking of both parity and time reversal symmetry without external electric or magnetic fields. This fact alone might provide an added advantage in designing robust DM sensing schemes.
Within a macroscopic coherent state scheme each unit cell is functioning as a sensor for DM detection, e.g., as illustrated above. Macroscopic coherence implies the synchronized response of the material. Each unit cell may have a very low sensitivity due to small coupling and size mismatch between the DM field and the unit-cell fields. In a macroscopically coherent quantum state, however, the unit-cell sensors add coherently and provide an Avogadro scaling enhancement of the signal-to-noise ratio, at least in principle.
In the case of MFs as a putative sensor the magnetic dipole moments, i.e., electron spins, are aligned with the electric dipole moment. The ground state of the MF is described by a coherent product state with a finite magnetic moment and polarization in each unit cell. The scalar contributions of Eq. (5) add constructively over a macroscopic volume ideally limited by the length scale min{λ a , L domain }, where L domain is the linear ferroic domain size and λ a is the axion de Broglie wavelength. Since this wavelength is on the order of meters or more [68], we anticipate that coupling to a macroscopic P M domain can still be ensured on appreciable sensor length scales in clean systems. The coherent nature of the response is sketched in Fig. 4. As an order-of-magnitude estimate (see Appendix B) we expect effective coupling to the axion to be χθ ae ∼ 7.7 × 10 −11 g ae , so that the total deposited energy becomes where we have assumed that the axion makes up the local DM density. Importantly, though the axion-induced magnetization will be on the order of δM ∼ O(1 aT) (in vacuum) due to the small coupling of Eq. (13), the generated magnetic flux scales with the area covered by the magnetometer, Φ ae ∼ L 2 domain δM , which can be enhanced experimentally. In practice, the size of the ferroic domains is likely dictated by Kibble-Zurek scaling as we elaborate on in the next Section.
In Appendix F we provide estimates of the signal-tonoise ratio for SQUID magnetometers, indicating that the Avogadro scaling makes a successful measurement of the response under ideal conditions conceivable. Yet, the tiny magnitude of the proposed effect makes it likely that background noise (e.g., thermal or quantum) will overwhelm the signal, and that material growth optimization and signal processing tools will be necessary ingredients if the axion-induced magnetization is to be successfully measured. Finally, we note that in reality decoherence and loss from various sources, such as acoustic phonons and (magnetic) impurities, may make a naive Avogadro scaling difficult to achieve.

IV. MATERIAL CANDIDATES
Since the discovery of multiferroicity in antiferromagnetic BiFeO 3 [69], the list of established MFs continues to grow and garner interest for their potential in next-generation microelectronic. However, ferromagnetic order in these materials, which is much soughtafter, is rare. Established ferromagnetic cases include CoCr 2 O 4 [70], DyFeO 3 [71], GdFeO 3 [72], Mn 2 GeO 4 [73][74][75][76], and more recently hexagonally grown Lu 1−x Sc x FeO 3 (h-LSFO) with x ≈ 0.4 [34,35]. Among these, the subset with P M is practically limited to the latter four, though the polarizations are not spontaneous in DyFeO 3 . In Sec. IV A we focus on h-LSFO by providing an overview of its properties and ab initio density functional theory calculations. Then, in Sec. IV B, we move on to discuss a broader range of other related material platforms, including nonspontaneously symmetry breaking cases.

A. Lu-Sc hexaferrite
The high-temperature paraelectric phase of the hexagonal ferrite h-LSFO adopts the P 6 3 /mmc space group, which has inversion symmetry. At the ferroelectric Curie temperature of T C P ≈ 1200 K, a spontaneous polarization emerges that breaks inversion symmetry. This symmetry-breaking structural distortion is caused by the condensing of a K 3 phonon mode manifesting in a trimerization of the Fe-O polyhedra, and a polar Γ − 2 phonon mode causing a rigid shift in the Lu ions. The resulting broken-inversion-symmetry crystal structure adopts the P 6 3 cm space group, as shown in Fig. 5 (a), in which the arrows indicate atomic displacement in the ferroelectric phase [34,77,78].
In contrast with conventional ferroelectrics, the hexagonal manganites and ferrites form topological defects (vortices) at their ferroelectric phase transition owing to the emergent U (1) symmetry of the structural trimerization at the onset of the ferroelectric phase transition [79,80]. This U (1) symmetry further breaks into a Z 6 ground state corresponding to the three trimerization directions of both the "up" and "down" polarization directions, the latter of which can be directly imaged using piezoresponse force microscopy [81]. Interestingly for our case, the domain sizes are determined by the Kibble-Zurek description of nonequilibrium dynamics during a driven phase transition: In this case the typical standard domain size of 10 µm can be controlled by the quench rate through the ferroelectric phase transition [80,82].
The spin structure of the trimerized Fe 3+ ions allows for a net magnetization in the c direction (in the so-called A 2 configuration) due to spin canting as illustrated with purple and blue arrows in Fig. 5(c). Ferromagnetic order develops below T C M ≈ 160 K, with much larger domains of typical size 100 µm [34].
To give materials-specific estimates of the expected axion coupling to the matter P and M we use ab initio calculations based on density functional theory (DFT) in the case of hexagonal LuFeO 3 . The calculation details are given in Appendix G. Here we summarize the key parameters used in the estimates of the strength of the effective matter coupling, namely the expectation value of the parallel polarization P and net ferromagnetism M .
We first confirmed that the A 2 magnetic ordering adopted by the P 6 3 cm polar structure is a frustrated 120°antiferromagnetic ordering with a slight out-of-plane canting. The calculated out-of-plane spontaneous magnetization was found to be 0.02 µ B /Fe, consistent with previous theoretical and experimental reports [78,83,84]. We find this value is relatively insensitive to choice of U eff in our DFT+U calculations for modest values of U eff . However, we note that an alternative magnetic order with A 1 symmetry is essentially degenerate with the A 2 magnetic order with a calculated energy of 0.06 meV/f.u. lower in energy than A 2 , which is also consistent with other experiments [60,85]. Such a small energy barrier between these two magnetic orders and the contradictory experimental measurements point to highly tunable magnetic orders that can be sensitive to synthesis conditions and strain, for example. We briefly note that this opens a possibility for increasing the observed out-of-plane net spontaneous magnetization through strain, chemical substitution, or doping.
Next we calculated the magnitude of the spontaneous polarization of the P 6 3 cm structure using the Berry phase approach. Debate over the definition of macroscopic polarization in bulk periodic systems by summation over local electric dipoles included whether it was truly a bulk effect, or a result of surface termination details, and how to resolve the multivaluedness of the polarization resulting from unit cell choices [86,87]. These issues were resolved with the introduction of the modern theory of polarization, which defined the bulk polarization from the Berry phase of the constituent wave functions, naturally prescribing the gauge freedom arising from the Berry phase [88,89]. To resolve this latter point in calculations, we generate a set of interpolated structures between the paraelectric P 6 3 /mmc and the ferroelectric P 6 3 cm to give a smooth interpolation of the Berry phase calculated electronic polarization. Our final calculated polarization was 10.5 µC/cm 2 in the out-ofplane direction.

B. Other candidates and nonspontaneous symmetry breaking
As a second promising material platform we mention the olivine compound Mn 2 GeO 4 [73,76]. In this material the ferroic orders originate from spiral spin-order, resulting from transverse conical chains that yield a net magnetization and electric polarization pointing along the c direction. Moreover, magnetic order here sets in at 47 K. Lowering the temperature leads to two additional first-order transitions, and below T C P ≈ 5.5 K both M P c are spontaneously generated. At low temperature the magnitude of the spontaneous electric and magnetic polarizations along c are about 6 µC/m 2 and 7 × 10 −3 µ B /Mn, respectively [73]. Ferroelectric domains on the scale of 500 µm have been observed in this compound [90].
Another possibility is to relax the assumption of spontaneous polarizations and consider material candidates where the magnetic A 2 state is induced by an external magnetic field, though we note that such cases lie beyond the scope of the theory in Eq. (6). The realization of a field-induced A 2 phase was proposed for the hexagonal manganites h-RMnO 3 with R = Ho, Er, Tm, Yb [91]. For instance, h-ErMnO 3 is believed to realize a fieldinduced ferromagnetic phase for a broad range of external magnetic fields. Furthermore, this may allow the reach of favourably large polarizations, of several µ B /Er and µC/cm 2 . However, the details of the magnetic phase diagram are complicated and still debated [92,93]. Antiferromagnetic materials naively appear unsuited for sensing the axion-matter coupling [of Eq. (5)] since the axion would couple with alternating sign to the antialigned spins and average out. However, ME field cooling techniques have demonstrated controlled switching of the antiferromagnetic state in Cr 2 O 3 heterostructures and made single ME domains possible in principle [94,95]. Hence, even antiferromagnetic magnetoelectrics may be surveyed for axion detection [96].
Finally, we mention the possibility of considering heterostructures with alternating layers of ferromagnetically and ferroelectrically ordered constituents, which also could enhance the space of relevant materials platforms.

V. CONCLUSION AND OUTLOOK
In this paper we have considered coupling between axion dark matter and macroscopic orders in multiferroics. The proposed scheme fits into the growing set of proposals to use quantum materials where the (broken) symmetry of the ground state enables the coupling between the axion and matter fields. Key advantages of the proposed MERMAID platform for DM detection include (i) reliance on the spontaneously broken symmetry state rather than external fields, and (ii) Avogadro scaling that allows long-range crystalline order in nature to build a coherent stack of sensors instead of experimenters.
Using materials-specific quantities, we estimated the axion-electron coupling in multiferroic materials with parallel magnetic and electric polarizations. The effective axion coupling in these materials is expected to be of the form P · M , similar to a linear magnetoelectric effect. We argued that the coupling is enhanced/resonant as the longitudinal magnon frequency matches the axion frequency. While we use an established MF material for these estimates, further studies should explore how to enhance the linear polarizations, and hence the coupling to axions, in real materials. Moreover, while here we focus on the case of spontaneously broken symmetries to induce electric-and magnetic-dipoles without using any external magnetic or electric fields, one can enhance the magnetization and polarization in the material with applied fields.
Temperature may in principle be used as a tuning knob to scan axion frequencies and to boost the response near critically. The sensing is based on extreme sensitivity of the multiferroic magnetic response near the ferromagnetic transition to the parameter changes induced by the axion field. In this respect the idea is similar to the transition-edge sensor.
We discussed Avogadro scaling as a result of the coherent response in materials with macroscopic quantum order. The quantum coherence of all unit cells means that the amplitude for scattering at each cell in the coherent state effectively enhances the signal-to-noise ratio. In the case of multiferroics this coherence emerges as a result of a macroscopic ferromagnetic and ferroelectric state.
High sensitivity magnetometers such as SQUIDs combined with discoveries of P M multiferroics, such as the ferrites Lu 1−x Sc x FeO 3 and the olivine compound Mn 2 GeO 4 , may provide a realization of the axion-matter coupling. In future work we would like to expand the sensitivity analysis and address the expected limitations in greater detail. As identification of suited condensed matter systems for axion and dark matter detection is in its infancy, exploitation of multiferroics for this purpose has the potential of opening a new avenue of research. There is also a possibility of considering dark matter detection schemes via the axion-matter coupling in a wider class of magnetoelectric materials than in the fairly limited class considered here. where the factor of m a appears when assuming an axion field a = a 0 exp(im a t), which is valid since the de Broglie wavelength is on the order of meters [68]. As explained in Sec. II the magnitude of a 0 is fixed by assuming that the axion make up the local DM density [14,43,44], ρ DM = 1 2 m 2 a |a 0 | 2 ≈ 300 MeVcm −3 . To quantify the size of the effect induced by Eq. (B1) one should ideally calculate

ACKNOWLEDGMENTS
where |Ψ is the many-body electronic ground state [i.e., Ψ = AΨ(r 1 , s 1 ; r 2 , s 2 ; . . . ; r N , s N ), A being the antisymmetrizing operator] with spin-orbit coupling included, which is a nontrivial problem. For the material estimates in Sec. IV and Appendix G |Ψ is determined with spin-orbit coupling included self-consistently using density functional theory. Here we take a simplified approach to give an order-of-magnitude estimate for the nonrelativistic outer electron orbitals. For the expectation value of Eq. (B1) one can generally insert a resolution of the identity such that In the subsequent steps we ignore the off-diagonal terms in Eq. (B3). For an insulating state the excited states are separated by the gap and their contributions to Eq. (B3) are expected to be small. A realistic estimate taking them into account will have to be done using ab initio methods and will be a subject of subsequent work. We thus approximate To have a ferroelectric polarization, i.e., to account for the misalignment of charge centres, one needs a linear combination of an even and odd orbital parity component, labelled by subscript ±: as illustrated in Fig. 6, resulting in FIG. 6. Illustration of the ferroelectric polarization P generated by a mixed orbital state, consisting of Ψ+ (s-wave) and Ψ− (pz-wave), with color representing charge density.
The electric polarization takes the form P = ρ e dr r|Ψ(r)| 2 = ρ e Ψ| r |Ψ , where ρ e is the charge density. This noticeably vanishes for any pure orbital state. However, for the mixed orbital state of Eq. (B5) we get if the (atomic) spatial variations take place over the scale of the Bohr radius, a B . Hence Eq. (B7) shows that essentially The spin expectation value defines a magnetic moment in the ferromagnetic phase, Ψ| σ |Ψ = − 2 gsµ B µ, where g s is the effective electron g-factor, µ B is the Bohr magneton, and µ is the total magnetic moment of the unit cell. The microscopic contributions are summed over individual spins in the unit cell (giving µ). In the ferromagnetic state we get δE ae ≈ a 0g 2 g s αe 2 P · µ = a 0g 2 g s αe 2 V uc P · M , where α is the fine structure constant and the subscript "uc" refers to the unit cell. This demonstrates how a macroscopic P · M coupling results from the axion-fermion coupling in a multiferroic. The above contributions add up coherently over a macroscopic number of unit cells, as dictated by the ferroic domain size. When comparing to the definition of χθ ae in Eq. (6), we arrive at the effective dimensionless coupling: χθ ae ≡ g ae a 0 m a m 2 e c 2 g s αe 2 ε 0 µ 0 ∼ 7.7 × 10 −11 g ae , (B10) Though the effective coupling is small, the total energy scales with volume and hence takes advantage of the coherent state of the material. For representative parameters g ae = 10 −13 , M 0 = 10 Oe, P 0 = 0.7 µC/cm 2 , and L domain = 1 mm the deposited energy from the axion is δE ae = 0.1 neV.
The full expressions for the first terms in the series expansion of M 0 (θ) as appearing in Eq. (12) are given by

(E1)
Appendix F: Signal-to-noise ratio To give an estimate of the signal-to-noise ratio for the schematic setup of Fig. 1, we here compare the magnitude of the axion induced magnetic field to presently achievable sensitivity levels for white noise in SQUIDs.
From the response of Eq. (10), we obtain a (vacuum) magnetic response of magnitude δM ≈ 0.35 aT (δM ≈ 0.31 aT) at ω a = ω M (ω a ω M ), i.e., on (off) resonance, at T − T C M = 1 µK and using the parameters of Fig. 3, which are motivated by the compound h-LSFO. We reiterate, however, that at resonance the response is controlled by the loss rate δM ∝ 1/η M , which experimentally may well be the limiting factor.
A key aspect of our proposal is that our signal will add coherently over the surface area of the sample, limited only by the maximum ferromagnetic domain size available. Given a linear domain size of 1 mm, we estimate the total magnetic flux measured over the whole area of the surface (1 mm 2 ). This gives on the order of 10 −8 Φ 0 , where Φ 0 ≡ (hc)/2e is the fundamental flux quantum.