Detecting axion dark matter with Rydberg atoms via induced electric dipole transitions

Long-standing efforts to detect axions are driven by two compelling prospects, naturally accounting for the absence of charge-conjugation and parity symmetry breaking in quantum chromodynamics, and for the elusive dark matter at ultralight mass scale. Many experiments use advanced cavity resonator setups to probe the magnetic-field-mediated conversion of axions to photons. Here, we show how to search for axion matter without relying on such a cavity setup, which opens a new path for the detection of ultralight axions, where cavity based setups are infeasible. When applied to Rydberg atoms, which feature particularly large transition dipole elements, this effect promises an outstanding sensitivity for detecting ultralight dark matter. Our estimates show that it can provide laboratory constraints in parameter space that so far had only been probed astrophysically, and cover new unprobed regions of parameter space. The Rydberg atomic gases offer a flexible and inexpensive experimental platform that can operate at room temperature. We project the sensitivity by quantizing the axion-modified Maxwell equations to accurately describe atoms and molecules as quantum sensors wherever axion dark matter is present.


I. Introduction
There is overwhelming astrophysical and cosmological evidence that approximately 85% of matter in the universe is in the form of non-luminous dark matter.Unfortunately, little is known about its nature beyond its gravitational influence on galactic and cosmological scales.While the search for historically popular models like Weakly Interacting Massive Particles (WIMPs) continues, it is important to expand experimental efforts to other well-motivated candidates.One such class of models, ultralight bosonic dark matter, is particularly favoured because of discrepancies that arise when simulations of structure formation with WIMP-like dark matter are compared to observations on galactic scales [1][2][3].These tensions are somewhat alleviated when the dark matter is modelled not as a WIMP-like particle but as an ultralight boson with de-Broglie wavelength comparable to small-scale galactic structures, which corresponds to a mass of order 10 −21 eV/c 2 .The axion is a famous candidate of bosonic dark matter originally predicted as the symmetry consequence of the Peccei-Quinn model as new physics beyond the standard model of elementary particle physics [8][9][10][11].It provides a "missing" elementary particle capable of naturally explaining why the charge-conjugation and parity symmetry (CP) are preserved in quantum chromodynamics (QCD) but are known to be violated in the electroweak interaction.This crucial conundrum has been coined the "strong CP" problem in elementary particle physics.While these QCD axions may be the dark matter [12][13][14], it is possible that the latter are pseudo-scalars which, however, do not solve the strong charge-parity problem.To distinguish them from the QCD axion, such pseudoscalars are called axion-like particles, and can usually be searched for with the same tools and techniques used to probe QCD axions.For the remainder of this article, we will use "axion" to refer to either a dark-matter candidate or a new QCD particle whenever the distinction is unimportant.It is interesting to note that, aside from their key role in high energy physics, axions have also been predicted to exist in condensed matter systems [15][16][17][18][19][20] where they manifest in the form of magneto-electric transport effects.Thus, a discovery of the elusive axion particle will not only resolve the strong CP problem in QCD, but could simultaneously track down the missing dark matter in the universe [21].We will focus on the latter aspect in this work.
The axion-modified Maxwell equations predict that a arXiv:2304.05863v2[hep-ph] 18 Jul 2023 FIG. 1. a) Magnetic field meditated axion-photon conversion.b) Typical resonator based setup for the detection of axion dark matter.c) Illustration of the axion-induced electric dipole transition in an atom mediated by a magnetic field.The dipole transitions can be detected spectroscopically.d) Overview of exclusion boundaries for ma and gaγγ.Experimentally excluded regions are marked by a solid line, while proposed experiments are marked by a dashed line.The projected exclusion bounds using axion-induced dipole transitions in Rydberg atoms achievable within one second, one hour, and one month measurement time are outlined.Also shown are constraints from: the CERN Axion Solar Telescope (CAST) [4], and the pulsar timing array [5].Other constraints shown on the plot for axions of a higher mass are from ADMX [6] and RBF-UF [7].
magnetic field converts axions of energy ω a m a c 2 + 1 2 m a v 2 (mass m a , velocity v) into photons of frequency ω a as sketched in Fig. 1(a) [11].A typical experimental setup exploiting this resonant axion-photon conversion for the detection of the galactic axion field consists of two inductively-coupled microwave resonators as depicted in Fig. 1(b) [22].A microwave photon, sourced by the axion field in the first resonator, is detected in the second resonator using a SQUID device, Rydberg atoms, or comparable single-photon detectors.The RBF-UF [7], ADMX [6] and HAYSTAC [23][24][25] experiments attempt to detect the axion field in the narrow 1-200 µeV mass range to which they are limited by construction.It should not be surprising then that most constraints on axions in the mass regime m a c 2 < 2µeV are astrophysical or cosmological, making alternative laboratory probes valuable.Current approaches, e.g.ABRACADBRA [26], ADMX SLIC [27] and SHAFT [28] using lumped-element circuits to compensate for the frequency mismatch lose sensitivity as the mass of the axion is lowered [29,30].Notably a resonant cavity for the detection of ultralight axions would be comparable in size to a small galaxy.Even if possible, the construction of a cavity is complex and cost intensive.
In this article, we analyze the possibility to search for the ultralight galactic axion field without relying on an advanced cavity resonator setup.Based on a rigorous quantization of the axion-Maxwell equations, we derive an effective Hamiltonian which describes dipole transistion induced by the axion-sourced electric field as sketched in Fig. 1(c).Thus, highly sensitive electric field detectors using quantum emitters and spectroscopic methods can be easily repurposed for the search of axion dark matter.
The axion-induced dipole transitions are particularly well suited to Rydberg states which feature long lifetimes, large electric dipole moments and polarizability [31].Being quantum sensors, Rydberg atoms can take advantage of quantized energy levels, quantum coherence, or manybody entanglement to enhance detection sensitivity compared to classical systems [32].Their aforementioned properties make them excellent candidates for electric field metrology in the radio-frequency regime, allowing sensitivities up to µVm −1 / √ Hz [33][34][35][36][37].
After introducing the axion-sourced electric field, we proceed to estimate the sensitivity of Rydberg atoms in a superheterodyne (superhet) detection configuration.For measurement campaigns over a period between seconds to months we project constraints in the ultralight mass regime that can outperform existing exclusion limits [see Fig. 1(d)].Details of the derivations can be found in the 'Methods' section and the Supplementary Materials (SM).

A. Axion-sourced electric field
Here, we give a heuristic derivation of the axionsourced electric field that induces dipole transitions in atoms and similar quantum emitters.A mathematically rigorous derivation is given in the 'Methods' section and the SM.The quantization of the axion-photon Lagrangian L = −g aγγ 0 µ0 aE • B yields, to linear order in the coupling constant g aγγ , the axion/light-matter interaction Hamiltonian where Ê and B, denote the electric and magnetic field operators, while â is the (pseudo-scalar) axion field, which interacts via a coupling constant g aγγ with Ê and B. 0 and µ 0 denote the vacuum permittivity and permeability.The axion field is measured in units of eV and the interaction constant g aγγ in units of (eV) −1 .Crucially, Ê is a physical observable and denotes the total electric field.In the presence of matter, one commonly introduces the displacement field that is sourced by the density of free charges ρ free in the electric Gauss equation, i.e., ∇ • D = ρ free .The electric field sourced by the bounded charges is described by the polarization density P .The displacement field can be understood as the electromagnetic field in the absence of matter and describes thus an external electric field or the quantized electric field in a cavity.The polarization density describes the electric field generated by the matter.Let us consider an ensemble of N quantum emitters and denote their eigenstates by |i, µ , where i labels the quantum emitter, and µ represents the collective electronic, vibrational and rotational quantum numbers.In the dipole approximation, the polarization density operator can be expressed in terms of the eigenstates as where δ(r) is the three-dimensional delta function, r i is the position of the i-th quantum emitter, and d (i) µ,ν is the transition dipole moment between the two eigenstates µ and ν.Assuming a stationary magnetic field, we can determine the dynamics of the displacement field sourced by the galactic axion field, which, when interacting with polarization, gives rise to the effective Hamiltonian and describes dipole transitions driven by the axionsourced electric field where the function S (x) describes the suppression of the axion-sourced electric field and depends on the axion mass m a and the length scale R B of the magnetic field.For small arguments, it scales with S (x) ∝ x 2 , i.e., the electric field is suppressed for small axion masses, constituting the main challenge to detect ultralight axions.In this work we assume S (x) = x 2 and a magnetic field length scale of R B = 0.1 m.Relation (5) implies that one can repurpose sensitive atomic and molecular detectors of electric fields to search for axion dark matter.The (classical) magnetic field can be considered as an experimental switch controlling the effective interaction strength between the axion field and the quantum emitters.This allows to distinguish between a possible axion signal and background electric fields.

B. Rydberg atoms as axion detectors
Ultralight axions in the galactic halo exhibit wave-like behavior and must be treated as a classical time-varying background field [38], a(t) = a 0 cos (ω a t + φ a ), where a 0 is the amplitude of the axion field, φ a its phase, and ω a m a c 2 + 1 2 m a v 2 its energy.Since the field has a small frequency dispersion described by the dark matter velocity distribution (average value 10 −3 c), a coherence time can be estimated to be τ C = 2 mac 2 10 6 , which defines the time scale over which the phase φ a can be considered to be constant [39][40][41].For axions in the ultralight mass regime, this time scale will always be much larger than the measurement time in question.The amplitude a 0 can be estimated via the galactic dark-matter energy density ρ = (m 2 a c 4 a 2 0 )/(2 3 c 3 ) = 0.3 • 10 15 eV m 3 .Resolving Eq. ( 5) for g aγγ , we find that the sensitivity for the interaction parameter is given by which assumes that axions constitute 100% of the dark matter in the universe.Equation 4shows that atoms with large transition dipole moments, in particular Rydberg states, couple strongly to the axion field.Using an advanced electromagnetically-induced transparence (EIT) based superhet detection protocol, an electric field of |E * | = 78• 10 −9 Vm −1 could be detected within 5000 s measurement time [33].Taking this setup as the basis for our sensitivity projection, we estimate that Rydberg atoms are able to detect minimal electric fields of Taking these values for reference, Fig. 1(d) shows the projected minimal detectable g aγγ, * evaluated via Eq.( 6) for a magnetic field of |B| = 5.6 T. It can be readily seen that Rydberg atoms can compete with the CAST helioscope bounds [42] in the mass regime m a c 2 = 5 • 10 −8 − 5 • 10 −6 eV.Also shown on the same plot are constraints inferred from the pulsar timing array [5].Rydberg atoms can thus set new leading experimental constrains while being operationally simple to realize.Due to their level structure, Rydberg atoms are particular sensitive in the small frequency regime.Their electric-field sensitivity is thereby relatively independent of the frequency for ω < 10 GHz.The exclusion bound established by the Rydberg atoms scales thus with g aγγ ∝ m −1 a , as the suppression factor in Eq. ( 5) scales with ∝ m 2 a and the axon field amplitude with ∝ m −1 a .For larger masses m a c 2 > 5 • 10 −6 eV, Rydberg atoms perform significantly inferior than the ADMX [6], and RBF-UF [7] experiments, whose frequency regimes allows to amplify the axion-sourced electric field in a resonator setup prior to detection.

C. Level structure of Rydberg atoms
Rydberg atoms are highly excited atoms featuring large dipole moments and polarizability.Each eigenstate can be characterized by the quantum numbers n, l, j, m.The energies µ depend mainly on the principle quantum number n ∈ {1, 2, 3, . . .} and the angular-momentum quantum number l = 0, . . ., n − 1 (also denoted by s, p, d, f . . .).The numbers j and m characterize the total angular momentum quantum number and its z projection, respectively.The spectrum of a Rubidium atom is depicted in Fig. 1(a), where each column represents a different l.The energy splitting between the highly excited Rydberg states n 1 rapidly decreases as n becomes large.In the presence of a static macroscopic magnetic field B e z , the eigenstates |µ are subject to a Zeemann shift, whose linear contribution reads as where µ B is the Bohr magneton, Lz is the z projection operator of the angular momentum, Ŝz is the projection of the spin, and g z is the corresponding gyromagnetic factor.The Zeeman effect can be used to tune specific energy states such that the monochromatic axion field fulfills a resonance condition.Because the electron is excited so far away from the nucleus, the transition dipole moments of two neighboring Rydberg states with principle quantum numbers n = n ± 1 are very large and scale as |d µ,ν | 2 ∝ n 2 , given that no optical selection rules are violated, c.f., Fig. 2(d).Thus, due to their energy spacing and the scaling of their dipole elements, Rydberg atoms are very sensitive to lowfrequency and quasistatic electric fields.

D. Sensitivity estimation
Recent experiments using EIT have demonstrated an outstanding sensitivity of Rydberg-atom superhet detectors for sensing electric fields [33].The superhet configuration consists of two low energy states |1 , |2 , and two Rydberg states |3 , |4 , whose energetic location are marked in Fig. 2 The coupled Rydberg states form the states |− and |+ , which exhibit a slow monochromatic energy splitting As ω a ≈ ω LO , we consider Ω(t) to be quasistatic.Both states are coupled to state |2 with Rabi frequency The resulting four-level system is depicted in Fig. 2(e).
The superhet detector senses the axion field via a change of the absorption of the probe laser traversing a cloud of Rydberg atoms, as sketched in Fig. 2(f).The linear absorption rate α(ω p ) for the resonance condition ω p = 2 − 1 is shown in Fig. 2(g) as a function of ∆ω C = ω C − ( 3 + 2 )/ .The absorption rate exhibits a dip around ∆ω C = 0 for Ω = 0 (not shown).This is the celebrated EIT affect, which is induced by a destructive interference between the transitions from |1 to |2 and |3 (|4 ) to |2 in the system's stationary state, preventing the absorption of the probe laser [43,44].
For a finite Ω, the atoms are not perfectly transparent closely around |∆ω C | ≈ 0. The energy splitting between the two states |− and |+ results into two new transparency frequencies at ∆ω C = ±Ω/2.Their interference shifts the ideal transparency, which would appear for Ω = 0 at ∆ω C = 0, and gives rise to high sensitivity of the absorption rate as highlighted in the left inset of Fig. 2(g), which compares the absorption rate for Ω = Ω LO and Ω = Ω LO + Ω a .A thorough analysis reveals that the signal-to-noise ratio (SNR) for a total measurement time t m is given by where N is the number of Rydberg atoms, τ = √ T c • t m is the effective measurement time, and T c is the effective coherence time ) is bounded by −0.5 < Γ < 0 and has a minor impact on T c .The Doppler effect thus enhances the dephasing rate γ 2 → γ 2 + ωpσ c Γ, which is discussed in detail in the SM.Equation ( 8) accounts for the projection noise in the atomic system, which sets the only fundamental limit for the SNR, while photon-shot and measurement noises have been neglected here [37].
As shown in the inset of Fig. 2(g), the SNR exhibits a turnover as a function of the local oscillator strength.For small Ω LO , it increases linearly with Ω LO until it reaches a maximum around Ω LO ≈ 5γ.For large Ω LO , the SNR vanishes with Ω −3 LO .The maximum value of the coherence time is T c ≈ (Ω p /Ω C ) 2 /γ which determines the optimal SNR.Using Eqs. ( 5) and ( 8), we find an explicit expression for the projected sensitivity of the axion field where K 3 , K 4 are defined in Eq. (7).The first fraction represents the inverse magnetic field B res required for establishing a resonance condition, which suggests using states featuring similar Zeemann parameters K 3 ≈ K 4 and a large detuning, as long as the external magnetic field B res < 10 T is experimentally feasible.In the proposed Rydberg-atom detector, the magnetic field fulfills thus two tasks: (i) it establishes the resonance condition between the Rydberg states, which enhances the sensitivity to the electric field; (ii) it generates the axion-sourced electric field according to Eq. ( 5).The second fraction represents the minimal electric field which can be detected with the superhet configuration.
To estimate the projected sensitivity achievable with Rydberg atoms, we consider the states Their energies and dephasing rates can be calculated using the ARC package [45].The dipole matrix element of the Rydberg states is |d 3,4 | = 6425ea 0 (Bohr radius a 0 ).In the ultralight axion regime the optimal magnetic field |B res | ≈ 5.6T determined by the first fraction in Eq. ( 10) is almost independent of the axion mass as m a c 2 d − c .The projected minimal electric fields and the related exclusion limits for g aγγ for the measurement times t m = 1 s, t m = 1 h and t m = 1 month have been discussed above.

III. Discussion
In this article, we have discussed the possibility to detect the galactic axion field by deploying quantum emitters such as atoms and molecules without using an advanced cavity resonator setup.Our proposal thus facilitates the search in the ultralight axion regime, where the required cavity length becomes unreasonable long.Based on a rigorous quantization of the axion-Maxwell equations, we have derived an effective Hamiltonian that microscopically describes dipole transitions in atoms, molecules, trapped ions driven by the axion-sourced electric field.This presents an exciting opportunity for repurposing existing electric field detectors based on atomic quantum sensors, which promise performance enhancement by means of quantum engineering [46], for axion detection.In this article, we proposed one such method using highly excited Rydberg states, whose large transition dipole elements make them excellent probes of axioninduced dipole transitions.We further conjecture that the dipole transitions can also enhance the sensitivity of other axion detectors, such as helioscopos [4] and phonon polaritons [47].We note that the direct detection scheme proposed here is different from previous protocols using Rydberg atoms, which either employ the atoms for an indirect detection of the axion-sourced photon in a cavity [48], or the coupling of the electron spin to the axion wind [49][50][51].
It is useful to consider further sensitivity improvement of our setup.The sensitivity may be improved for longer measurement campaigns over several year, similar to the one performed by the BACON collaboration [52] searching for the time variation of fundamental constants.In this case one would also have to carefully account for the stochastic fluctuations of the axion dark matter field.Another avenue might be to seek improvements in the measurement process itself to facilitate sensing of weaker electric fields, e.g., by using stronger probe fields.The sensitivity estimate of the superhet detector for weak probe fields shows that the SNR is proportional to Ω p /Ω C , i.e., it improves for stronger probe fields.However, this requires the development of non-perturbative theoretical methods, e.g., based on the photon-resolved Floquet theory [53], to accurately predict the spectroscopic signatures of the axion dark matter.One interesting possibility proposed in [54] is to use trapped ion crystals which can achieve electric field sensitivities of 100 nVm −1 .Moreover, trapping the Rydberg atoms in an optical lattice can further help to mitigate the Doppler effect, and would allow for complex quantum operations to mitigate measurement noise [32,55].

A. Quantization of the axion Maxwell equations
The axion field, being a pseudo-scalar, interacts with the electric and magnetic fields via the Lagrangian term L = −g aγγ 0 µ0 aE • B, where E, B, and a denote the classical electric, magnetic, and axion fields, respectively.After deriving the Euler-Lagrange equations, we obtain the axion-Maxwell equations [11,39] ∇ where J is the electric current density.The constants have been described in the 'Results' section.Clearly, These equations reduce to the common Maxwell equations for g aγγ = 0.As we show in details in the SM, the canonically quantized axion-light-matter Hamiltonian reads as where the canonical electric field operator has been distributed in the transversal and longitudinal fields Êc = Ê⊥ c + Ê c .The transversal contribution is defined by where φ is the electrostatic vector potential .The vector potential and the magnetic field operators are denoted by B and Â, respectively.The axion field is denoted by â, and π denotes its conjugated momentum.Particles with charge q η and mass m η at positions rη having momentum pη are labeled by η ∈ {1, . . ., N p }.The coupling between light and matter reflects the minimal coupling principle.The relation between the canonical and the physical (observable) electric field operators Ê is established via the relation The canonical operators for all other operators agree with the physical operators.The field operators are quantized as The photonic operators dk,λ quantizing the electromagnetic field are labeled by wavevector k and polarization λ ∈ { , ↔}.The frequencies of the photonic modes are given by ω k .The quantization volume is denoted by V .The unit vectors e k,λ describe the direction of polarization and are perpendicular to the wavevector, i.e., e k,λ • k = 0.This fact readily ensures that the electric and magnetic Gauss equations in Eqs.(11) and ( 13) are fulfilled.Using the Heisenberg equations of motion and the relation between canonical and physical electric fields in Eq. ( 17), we can derive the axion-Maxwell equations Eq. ( 11)-(15) in the classical limit by generalizing the common derivation of the Maxwell equations [56].

B. Multi-polar Hamiltonian
The minimal-coupling Hamiltonian in Eq. ( 16) is inconvenient to deal with because of the vector potential Â that appears quadratically.This causes more complicated expressions when calculating, e.g., the spectroscopic response of quantum systems.Moreover, the vector potential Â is not a physical observable.For this reason, we will bring the Hamiltonian into the so-called multi-polar form that contains only the electric and magnetic fields [56].This is done by means of the Power-Zienau transformation [57], which we review here shortly.More details can be found in the SM [39].Here, we restrict our investigation to systems with bounded charges, such as atoms, molecules, and similar quantum emitters.
We will specify to atoms in the following for clarity.Accordingly, the summation over η in Eq. ( 16) will be modified into a double summation over i ∈ {1, . . ., N }, labeling atoms, and j ∈ {1, . . ., N c }, labeling the charges associated with a particular atom.The center of mass of atom i will be denoted by R i , while the position of the charges belonging to this atom is denoted by r i,j .
The Power-Zienau transformation is defined by the unitary operator where denote the macroscopic polarization and the polarization generated by the charges η = (i, j), respectively.For simplicity, we assume that there are no free charges in the system, such that φ = 0.The Power-Zienau transformation leaves all operators in Eq. ( 16) invariant, except of the canonical electric field, which becomes and is called the displacement field D(r).Importantly, the displacement field is quantized in terms of the photonic operators dk,λ : and not the canonical electric field Êc (r).However, one should keep in mind that the electric field Ê is the actual physical observable.Away from the matter where P (r) = 0, the displacement field is equivalent to the canonical electric field D(r) = 0 Êc (r).
The transformed axion-light-matter coupling Hamiltonian (∝ a Êc • B) now reads as where we have used that c 2 = 1/( 0 µ 0 ).The first term is equivalent to H ALM in Eq. ( 1) after substituting Ê according to Eq. ( 2).As the second term is proportional to g 2 aγγ , it can be safely neglected.The other terms in the Hamiltonian in Eq. ( 16) transform according to the common Power-Zienau transformation and are given in the SM [39].

C. Rydberg atoms as superheterodyne detectors
At present, there is a variety of electric field detectors that deploy an EIT coupling scheme in Rydberg atoms.The common EIT three-level configuration is easy to realize experimentally and allows already for excellent sensitivity to weak electric fields [34][35][36][37].Yet, in a three-level configuration, the signal field does only perturbatively couple off-resonant Rydberg states to each other, which limits the optimal sensitivity.To further improve the sensitivity, the superhet detection scheme has been developed, which takes advantage of a resonance condition of two Rydberg states, that features a large transition dipole matrix element [33].
A superhet detector deploying Rydberg states can be described by an effective four-level system whose states are denoted by |1 , |2 , |3 , |4 .Their positions in the energy spectrum are shown in Fig. 2. The corresponding Hamiltonian reads as where x with x = 1, 2, 3, 4 denote the level energies.
A description of the parameters ω p , ω C , ω a , ω LO , Ω p , Ω C , Ω a , Ω LO is given in the 'Results' section.We consider the time-dependent Rabi frequency Ω(t) = Ω LO + Ω a e i(ωa−ω LO )t to be quasistatic, and use the difference of the measurement signal for Ω = Ω LO − Ω a and Ω = Ω LO + Ω a to estimate the SNR.We describe the dynamics of the four-level system in the presence of dissipation by the Bloch equation where the dissipator is defined by and γ 2 , γ 3 and γ 4 denote the inverse lifetimes of the states |2 , |3 , |4 , respectively.
Transforming the subsystem of the states |3 , |4 into a frame rotating with ω LO , the effective energies of the resulting mixed Rydberg states become − = 3 − Ω/2 and + = 3 + Ω/2, which are sketched in Fig. 2(e).Both mixed Rydberg levels couple to state |2 with Rabi frequency Ω C / √ 2 and have a dephasing rate of γ = 1 2 (γ 3 + γ 4 ).As common practice in EIT experiments, we aim to detect the presence of the axion field by measuring the absorption of the probe field.When certain resonance conditions are fulfilled by the probe and coupling fields, the atoms become transparent to the probe laser within a very narrow window around ω p = ( 2 − 1 )/ .As the resonance conditions are modified in the presence of the axion field, the absorption significantly depends on Ω a , heralding the presence of axion dark matter.
The absorption rate α(ω p ) of the probe beam in Fig. 2(f) is proportional to the susceptibility α(ω p ) ∝ Im χ(ω p ) [56].As explicitly shown in the SM, the susceptibility can be expressed in terms of the operator X = N i=1 Xi with Xi = tm 0 dt P(i) 2,1 (t)ie iωpt + h.c., where the subscript i labels the N atoms in the ensemble and t m denotes the measurement time.As Û (t) is the time-evolution operator associated with the Bloch equation in Eq. ( 30), X acts nonlocal in time.
Explicitly, the relation between the expectation value of X and the susceptibility in linear order of Ω p is given by where ρ N denotes the atom density.The expectation value • = tr [•ρ st ] is defined in terms of the stationary state of the system for Ω p = 0 that reads ρ st = |1 1|.
For this reason, X can be considered as a susceptibility operator, and be employed to estimate the SNR.A straightforward perturbative treatment in Ω p , which generalizes the standard treatment of the EIT [43] to the four-level system under investigation, and a velocity average to account for the Doppler effect, reveals that the susceptibility is given by where we have defined where erf (z) denotes the complex-valued error function, and we have introduced for a rotational reason.
The absorption rate (i.e., the imaginary part of the susceptibility) is depicted in Fig. 2(g) as a function of ∆ω C = ω C −( 3 + 1 )/ for ω p = ( 2 − 1 )/ , for a suitable local oscillator strength Ω LO , and two different Ω a = 0 (dashed) and Ω a > 0 (solid).In the inset we depict the susceptibility for frequencies close to |∆ω C | γ, where we observe a clear difference between the cases Ω a = 0 and Ω a > 0. For this reason, we choose the coupling frequency to be ∆ω C = 0, for which the signal for small Ω a is given by where we have defined ) .Likewise, we find that the variance of X for vanishing Ω a is given as where we have approximated Ω → Ω L0 .Putting everything together, we obtain the SN R ≡ in Eq. (8).

S1. Quantization of the axion-Maxwell equations
In the article, we propose to deploy Rydberg atoms as highly sensitive dark matter detectors.In order to describe these quantum sensor in a consistent framework, we here perform a rigorous quantization of the axion-Maxwell equations.The structure of the section is as follows: In Sec.S1 A we introduce the axion-Maxwell equation.In Sec.S1 B, we show how to accurately quantize these equations by microscopically deriving the Hamiltonian.In Sec.S1 C, we perform a Power-Zienau transformation, to bring it into a suitable form to describe spectroscopic experiments.

A. Classical axion-Maxwell equations
The relativistic Lagrangian density describing the dynamics of the pseudo-scalar axion field a is given by [58] where the covariant electromagnetic field tensor F αν = ∂ α A β − ∂ β A α can be expressed in terms of the covariant 4-vector potential A α = (φ, A).The labels α, β, γ, δ ∈ {ct, x, y, z} represent space-time coordinates.The covariant electric 4-current is given as J α = (cρ, J ).The dual tensor of the electromagnetic field is defined by F αβ = 1 2 αβγδ F γδ , where αβγδ is the dyadic tensor.The electric and magnetic fields can be obtained via Using the Euler-Lagrange equations and the Bianchi identity, we find the following relativistic equations of motion for the above-introduced fields which in terms of the electric, magnetic, and axion fields reads as These are the axion-Maxwell equations in SI units, which reduce to the common Maxwell equations for g aγγ = 0.The axion field is given in units of eV, while the axion-photon coupling is given in units of eV −1 .

B. Quantization
Quantization of the axion-Maxwell equations means finding the microscopic Hamiltonian and the canonical commutation relation of the operators.We first state the result and then prove its correctness by deriving the non-relativistic classical axion-Maxwell equations using the Heisenberg equations of motion.The derivation is a generalization of the derivation of the common Maxwell equations for g aγγ = 0, which can be found in, e.g., Ref. [56].
In the Coulomb gauge ∇ • A = 0, the canonically quantized axion-light-matter Hamiltonian reads as where the canonical electric field operator has been distributed in the transversal and longitudinal fields Êc = Ê⊥ c + Ê c .The transversal contribution is defined by ∇ • Ê⊥ c = 0, while the longitudinal contribution is given by Ê c = Êc − Ê⊥ c .The electric potential φ is defined via Ê c = −∇ φ.The vector potential and the magnetic field operators are denoted by B and Â, respectively.The axion field is denoted by â, and π denotes its conjugated momentum.Particles with charge q η and mass m η at positions rη with momentum pη are labeled by η ∈ {1, . . ., N p }.The coupling between light and matter reflects the minimal coupling principle.The commutation relations of the operators are given below.The relation between the canonical and the physical (observable) electric field operators Ê is established via the relation The canonical operators for all other operators agree with the physical operators.The field operators are quantized as The photonic operators dk,λ quantizing the electromagnetic field are labeled by wavevector k and polarization λ ∈ { , ↔}.The quantization volume is denoted by V .The unit vectors e k,λ describe the direction of polarization and are perpendicular to the wavevector, i.e., e k,λ • k = 0.This fact makes sure that the electric and magnetic Gauss equations in Eqs.(S4) and (S6) are fulfilled.The frequencies of the photonic modes are given by ω k .
We are now in position to derive the non-relativistic axion-Maxwell equations.The electric and magnetic Gaussian equations are fulfilled since where we have used that ∇ • Ê⊥ c = 0 and ∇ • E c = ρ/ 0 because of the parameterization in Eqs.(S14) and (S15).Using now the relation of the canonical and physical electric fields in Eq. (S10) we obtain the electric Gauss equation.For the same reason, we also find ∇ • B = 0, i.e., the magnetic Gauss equation in Eq. (S6).
The Faraday and Ampere equations can be constructed by means of the Heisenberg equations of motion The derivation follows essentially the same lines as the common Maxwell equations upon simply replacing Ê by Êc .
Using the commutation relations in Eq. (S18), we obtain In the magnetic equation, we can replace Ê⊥ c → Êc as ∇ × Ê c = 0 since Ê c = −∇ φ for the electrostatic potential φ defined in Eq. (S16).Moreover, the notion and derivation of J ⊥ and J , that fulfill J = J + J ⊥ , can be found in standard textbooks on quantum electrodynamics.We proceed to combine the longitudinal and transversal electric equations Resolving the magnetic equation in Eq. (S21) for cg aγγ ∇ × â B and inserting into the electric equation Eq. (S22), we obtain Using now the relation of the canonical and physical electric field operators in Eq. (S10), we readily find the Faraday and Ampere equations in Eqs.(S7), (S5) and (S8), respectively.We continue to construct the Heisenberg equations of motion for the axion field, which read Upon replacing the canonical electric field by the physical one according to Eq. (S10), we obtain Deriving Eq. (S24) with respect to time and inserting Eq.(S26), we readily obtain the axion equation of motion in Eq. (S8).

C. Multi-polar Hamiltonian
The minimal-coupling Hamiltonian in Eq. ( 16) is inconvenient to deal with because of the vector potential Â that appears quadratically.This causes more complicated expressions when calculating, e.g., the spectroscopic response of quantum systems.Moreover, the vector potential Â is not a physical observable.For this reason, we will bring the Hamiltonian into the so-called multi-polar form that contains only the electric and magnetic fields [56].This is done by means of the Power-Zienau transformation, which we review here shortly [57].To this end, we restrict our investigation to systems with only bounded charges, such as atoms, molecules, and similar quantum emitters.We will specify to atoms in the following for clarity.Accordingly, the summation over η in Eq. ( 16) will be modified into a double summation over i ∈ {1, . . ., N }, labeling atoms, and j ∈ {1, . . ., N c }, labeling the charges associated with a particular atom.The center of mass of atom i will be denoted by R i , while the position of the charges belonging to this atom is denoted by r i,j .
The Power-Zienau transformation is defined by the unitary operator where denote the macroscopic polarization and the polarization generated by the particles η = (i, j) carrying charge q i,j , respectively.The Power-Zienau transformation has the following effect on the system operators: Consequently, the electromagnetic field operators become Thus, only the transverse canonical electric field becomes modified and is shifted by the polarization.The transformed electric field is called the displacement field D(r).Importantly, the displacement field is quantized in terms of the new photonic operators d k,λ : and not the canonical electric field Êc (r).However, one should keep in mind that the electric field Ê is the actual physical observable.Away from the matter where P (r) = 0, the displacement field is equivalent to the canonical electromagnetic field D(r) = 0 Êc (r).In the absence of free charges, the displacement field is entirely transverse and reads as as the longitudinal polarization fulfills P c = 0 Ê c = 0 in the absence of free charges as considered here.The multi-polar Hamiltonian reads as where the transformed matter Hamiltonian is given by The last term represents the polarization self-energy.The electromagnetic field Hamiltonian is now given as where in the second equality, we have expressed it in terms of the new photonic operators.The transformed lightmatter coupling now reads as The first term describes the coupling of the electric dipole moment to the displacement field.The second term is the coupling of the magnetic dipole moment (thoroughly defined in Ref. [56]) to the magnetic field.The third term describes the coupling of the electric dipole density to the magnetic field.Yet, as it is divided by the rest energy of the charges m ij c 2 , this term can be safely neglected.
The free axion Hamiltonian remains unchanged, while the transformed axion-light-matter coupling Hamiltonian now reads as where we have used that c 2 = 1/( 0 µ 0 ).

A. Decoupling transformation
To simplify the following analysis, we consider a strong static magnetic field B, which is in agreement with our experimental proposal.The total magnetic field thus reads as where Bf denotes the fluctuations of the field.We now carry out the unitary transformation which has the following effect on the system operators Thus, the transformed Hamiltonian now reads as In the transformed Hamiltonian, the axion-polarization coupling is mediated only by the fluctuations of the magnetic field Bf , which can be thus neglected.As we see in the forth line, the polarization couples to the displacement field.For this reason, we have to determine its dynamics in order to estimate the axion field can induce some physically measurable effect on the polarization.The equations of motion of the displacement field reads as which is a linearized version of the axion-Maxwell equations in Eqs.(S4) -(S8).As c 2 is a large number, we have also constructed the equation of motion for the magnetic field fluctuations.For simplicity, we assume that the axion field is in a coherent state, such that its impulse has the following time evolution i.e., we consider it as a classical source term.

B. Solution of the axion-Maxwell equations
To solve the axion-Maxwell equations, we derive the first equation in Eq. (S42) with respect to time, and insert then the magnetic field expression of the second equation.Using identities of vector calculus, we readily arrive at the inhomogenous wave equation where the axion field acts a source term.This equation has the well-known solution which can be derived, e.g., using the Green's function formalism [59].In our case, the source term explicitly reads as Without loss of generality, we consider the position r = 0 in the following, as we can place the origin of the coordinate system at the position where we want to sense the displacement field.Expanding the magnetic field in terms of spherical harmonics Y n,l (θ , ϕ ), we find the following expansion for the displacement field Dl,m (0, t). (S48) The partial displacement fields can be expressed as where the angular dependence is described by the coefficients We emphasize that Eqs.(S48) and (S49) comprise an exact solution of the inhomogeneous wave equation in Eq. (S44).
For a practical reason, we want to bring Eq. (S48) into a more compact form.To this end, we assume that all B l,m (r ) are parallel to the magnetic field at the origin B(0).In this case, the displacement field can be expressed as D(0, t) = 0 cg a,γγ a 0 cos(ω a t + φ)B(0)S B (m a ) , (S51) where we have introduced which incorporates the form of the magnetic field.The phase φ will be neglected in the following.It is interesting to analyze the form factor F for small masses m a , such that e −imar /c ≈ 1. Rescaling the argument of the magnetic field B (r) = B(αr), we find that F → α 2 F .We thus conclude the following scaling behavior of the suppression factor where R B parameterizes the spatial extend of the magnetic field, and the suppression function scales as S(x) ∝ x 2 for small x.The scaling for large x depends on the particular form of the magnetic field.

C. Example
To illustrate the physical implications of the general solution of the wave equation in Eq. (S48) on a concrete example, we consider an exponentially decaying magnetic field, The maximum amplitude of this field has the following scaling properties: For small axion masses m a R B 1, this relation shows that the axion-sourced displacement field is strongly suppressed, in agreement with Eq. (S53).In contrast, for m a R B 1, the axion-sourced displacement field becomes independent of the axion mass and approaches the solution of the translationally invariant system.

D. Effective Hamiltonian
Based on the previous derivations, we introduce here an effective Hamiltonian, which we use to predict the sensitivity of Rydberg atoms.First, we neglect the fluctuations of the electromagnetic fields and the axion field in the Hamiltonian in Eq. (S41).In doing so, the Hamiltonian reduces to where the expectation value is defined by • = tr [•ρ st ] with the stationary state of the system for Ω p = 0 given by ρ st = |1 1|.Due to its close relation to the susceptibility (and thus to the absorption rate), we refer to X as the susceptibility operator in the following.In Sec.S3 H, we thus take the operator X as the basis to determine the SNR of the experimental setup.In doing so, we accurately account for the projection noise, which is the only fundamentally limiting noise source in the experiment [36].Other noises, like photon shot noise or detection noise are not considered here.

C. Susceptibility
The susceptibility can be expressed in terms of the linear response function S(r, t), which determines the polarization in response to an external probe electric field Êp (r, t) = E p,0 cos(kr − ω p t) [56], i.e., We assume here that the probe frequency is close to resonance with the transition between |1 and |2 in the Hamiltonian in Eq. (S63).For this reason, it is sufficient to parameterize the polarization operator as Applying standard first-order perturbation theory to the Hamiltonian in Eq. (S63), we find that the linear response function is given by [56] S(r, t) where ρ N denotes the atom density.Here, we have assumed that all atoms are equal, such that we can neglect the superscript i.The time-evolution is thereby determined for Ω p = 0.The linear susceptibility can be obtained from the linear response function via Fourier transformation For later purpose, we have introduced the correlation function Its Fourier transformation is defined via In the approximation in Eq. (S73) we have deployed the Kubo-Martin-Schwinger relation C(ω) = C(−ω)e β ω assuming room temperature 1/β ≈ 25meV and probe frequency ω p > 1eV.

D. Calculation of the signal
We calculate the expectation value of the operator Xi in Eq. (S68) in the presence of the probe field Êp (r, t) = E p,0 cos(kr − ω p t) in first-order perturbation theory in Ω p = |d 1,2 In agreement with the RWA in the Hamiltonian in Eq. (S63), we neglect the fast oscillating terms e ±2ωpt .In doing so, we find Setting ϕ = π/2 and comparing with Eq. (S73), we obtain which establishes the desired relation between the operator X and the susceptibility.

E. Calculation of the noise
In the same manner, we evaluate the variance of the observable X .As Ω p is assumed to be small, it has a minor influence on the result and we can set Ω p = 0. Explicitly, the variance can be evaluated to be Comparing this result with Eq. (S73), we infer that where we have again deployed the Kubo-Martin-Schwinger relation.

F. Electromagnetically-induced transparency in a four-level system
In Secs.S3 D and S3 E, we have derived formal expressions for the signal and the noise in terms of the susceptibility.It remains to express the susceptibility in terms of the microscopic system parameters.We achieve this by solving the Bloch equation in Eq. (S64) and finding an explicit expression for the matrix element ρ 12 (t).To connect this to the susceptibility, we use the relation in first-order perturbation theory where we have introduced the correlation function C(t) for a notation reason.Performing a (normalized) Fourier transformation, we find where normalization is necessary to avoid a divergence.We note that due to the Kubo-Martin-Schwinger relation and within the RWA approximation C(ω p ) ≈ C(ω p ).We thus obtain C(ω p ) and consequently the susceptibility in Eq. (S73) via To find an expression for ρ 21 , we generalize the standard treatment of the EIT in three-level systems (see, e.g., textbook of Scully [43]) to the four-level system in Eq. (S63).To start with, we transform the Hamiltonian into a frame rotating with frequencies ω p and ω C , such that the resulting time-independent Hamiltonian reads as We aim to find an expression for ρ12 in the leading order of Ω p .As inspection of the Bloch equations reveals that ρ11 ≈ 1 ρ22 , ρ2− , ρ−2 , ρ−− , ρ2+ , ρ−+ , ρ++ ∝ Ω 2 p , we neglect corresponding terms.The three remaining equations thus are which we have already transformed into Laplace space under the assumption that ρ 12 (0) = ρ 1− (0) = ρ 1+ (0) = 0. Operators in Laplace space are marked by an inverted hat, i.e., •.Moreover, we have defined To establish the second equality, we have transformed the density matrix into the lab frame via ρ 12 (t → ∞) = ρ12 (t → ∞)e iωpt , and carried out the scaled Fourier transformation in Eq. (S82).

G. Temperature dependence
Usually, the superhet detector is operated at a finite temperature (e.g., T = 300 K), which will deteriorate its sensitivity.Two effects have to be taken into account.Firstly, the dephasing rates of the excited states γ 2 , γ 3 , and γ 4 will increase.The temperature-dependent dephasing rates can be calculated using the Alkali-Rydberg-Calculator (ARC) package [45].Overall, this effect is relatively small.Secondly, the temperature-induced motion of the atoms will give rise to a Doppler shift concerning the probe and coupling beams.To mitigate this effect, the probe and coupling beams shall propagate in opposing directions, as sketched in Fig. 2(f) in the article.Consequently, both laser frequencies are Doppler shifted as where λ p = 2πc/ω p (λ C = 2πc/ω C ) is the wavelength of the probe (coupling) beam, and v is the velocity of an atom parallel to the beams.The Doppler-averaged density matrix elements can be obtained by evaluating the integral where m R denotes the mass of the Rydberg atoms, and T is the temperature of the thermal cloud.For a notation reason, we have introduced the parameters This integral can be solved analytically.Using that Im c 2 > 0, we can transform  Without local oscillator Ω LO = 0, the absorption rate exhibits a broad single dip at ∆ω C = 0, which is the celebrated EIT.Since Im χ(ω p ) ∝ γ, the EIT is not complete for T = 300 K, as γ increases for increasing temperature.The width of the dip scales with ∝ 1/γ 2 and thus decreases for larger temperatures.
For a finite local oscillator strength Ω = 0.1 MHz, the two levels |3 and |4 mix, which leads to constructive interference at ∆ω C = 0 and destroys the EIT.This generates a peak in the absorption rate, whose width is proportional to γ.The hight of the peak sensitively depends on Ω.

H. Signal-to-noise ratio
As the susceptibility sensitively depends on Ω at the resonance condition ∆ω C = 0, the superhet detector operators under these circumstances, in which (S98) Based on Eqs.(S78), (S80) and (S97), we can now evaluate the SNR for the operator X .Recalling that the total Rabi coupling between the Rydberg states consists of the local oscillator and the axion-induced dipole transitions Ω = Ω LO + Ω a , we find that the signal for small Ω a is given by Finally, we recall that we have considered Ω as quasistatic throughout the derivations.Taking the time dependence Ω(t) = Ω LO + Ω a e i(ωa−ω LO )t into account, we find that the measurement signal will slowly oscillate with frequency ω a − ω LO and amplitude in Eq. (S99).Thus, the SNR will remain unchanged in the case of heterodyning.
(a) and (b).The states |1 and |2 are coupled by the probe laser of frequency ω p and Rabi frequency Ω p , while the states |2 and |3 are coupled by the coupling laser of frequency ω C and Rabi frequency Ω C .The Rydberg states are coupled by the axion field of frequency ω a and Rabi frequency Ω a = g aγγ ca 0 |d 3,4 • B| / • S. The corresponding energies 3 and 4 are assumed to be close to resonance 4 − 3 ≈ ω a .The resonance condition can be adjusted using the Zeeman effect [see Fig. 2(c)].To enhance the signal, the superhet detector uses a local oscillator with frequency ω LO = 4 − 3 and Rabi frequency Ω LO Ω a which heterodynes the axion field.

FIG. 2
FIG. 2. a) Energy levels of a Rubidium atom for n ≥ 5 resolved for the angular momentum orbitals s, p, d, f .(b) depicts a magnification of highly-excited f orbitals, that are coupled to the Rydberg level |3 with quantum numbers n = 100 and l = 2 (d orbital).(c) is the same as (b), but for a finite external magnetic field B = Bres, which establishes a resonance condition between levels |3 with (n, l, j, m) = (100, 2, 3/2, 1/2) and |4 with (n, l, j, m) = (99, 3, 5/2, 1/2).Gray lines depict the allowed dipole transitions for a linearly-polarized driving field.(d) Transition dipole matrix elements between the d and the f orbitals.(e) Effective four-level system considered in the description of the Rydberg-atom superheterodyne detector.Levels |− and |+ appear upon coupling the states |3 and |4 .(f) Sketch of the spectroscopic setup, in which the probe and coupling fields propagate in the opposite direction to mitigate the Doppler effect as explained in the SM [39].(g) Absorption rate α(ωp) of the probe beam as a function of ∆ωC = ωC − ( 3 + 2)/ at resonance ωp = 2 − 1.The left inset depicts the absorption rate for small |∆ωC | γ.The right inset depicts the signal-to-noise ratio (SNR) as a function of the local oscillator ΩLO at ∆ωC = 0. Other parameters are Ωp = 1.7 MHz, ΩC = 23 MHz, Bres = 5.6 T, temperature T = 300K.

FIG
FIG. S1.(a) Absorption rate as a function of ∆ωC for three different temperatures.(b) Doppler-effect form function Γ(x) and its derivative Γ (x).(c) SNR as a function of local oscillator strength for three different temperatures.Parameters are Ωp = 1.7 MHz, ΩC = 23 MHz, and Bres = 5.6 T. If not specified differently, the parameters are the same as in Fig. (2)(g).