Observation of a charge-neutral muon-polaron complex in antiferromagnetic Cr$_2$O$_3$

We report a comprehensive muon spin rotation ($\mu$SR) study of the prototypical magnetoelectric antiferromagnet Cr$_2$O$_3$. We find the positively charged muon ($\mu^+$) occupies several distinct interstitial sites, and displays a rich dynamic behavior involving local hopping, thermally activated site transitions and the formation of a charge-neutral complex composed of a muon and an electron polaron. The discovery of such a complex has implications for the interpretation of $\mu$SR spectra in a wide range of magnetic oxides, and opens a route to study the dopant characteristics of interstitial hydrogen impurities in such materials. We address implications arising from implanting a $\mu^+$ into a linear magnetoelectric, and discuss the challenges of observing a local magnetoelectric effect generated by the charge of the muon.


I. INTRODUCTION
Positively charged muons implanted into semiconductors and insulators often form muonium (Mu ¼ ½μ þ e − ), a hydrogenlike charge-neutral bound state. It is conventionally referred to as a paramagnetic center as the bound electron is unpaired and its spin is decoupled from all the other electrons. Since the electronic structure of Mu in a solid is virtually identical to that of hydrogen, Mu has been studied extensively using muon spin spectroscopy (μSR) to learn about interstitial hydrogen, one of the most ubiquitous defects in semiconductors. In a μSR experiment, spinpolarized muons (μ þ ) are implanted into the sample of interest, and the subsequent time evolution of the spin polarization is observed, providing an accurate measurement of the magnetic coupling of the muon spin with its environment [1]. Mu centers have been studied in a wide range of materials, providing direct information about the electronic structure of hydrogen defects as shallow or deeplevel dopants [2][3][4][5]. So far, the study of such charge-neutral muon states has been limited to nonmagnetic materials, since paramagnetic Mu is widely assumed to be subject to strong depolarization in the presence of magnetic moments such as those on Cr 3þ [3], and, with the exception of MnF 2 [6], no Mu has been confirmed in magnetic materials.
Here, we present strong evidence for a charge-neutral muon state in the antiferromagnet Cr 2 O 3 . In particular, our data, in conjunction with detailed density-functional-theory (DFT) calculations, make a compelling case for the existence of a muon-polaron complex, where the positive muon is bound to an oxygen, and an excess electron localizes on a nearby Cr ion, changing its valence to Cr 2þ (3d 4 ). The degeneracy of the now occupied e g orbital is lifted via a lattice distortion, leading to a Jahn-Teller (JT) polaron [7][8][9] on the Cr ion. Crucially, the resulting JT-stabilized muonpolaron complex is not paramagnetic and therefore distinct from Mu, since the bound electron is strongly coupled to the 3d electrons of the Cr host ion. Therefore, no signatures conventionally associated with a charge-neutral state are displayed, concealing its existence. However, despite its inconspicuous signal, the presence of such a complex has a significant impact on the location and stability of muon stopping sites, and the local fields experienced there.
This discovery of a charge-neutral muon-polaron complex in Cr 2 O 3 suggests that neutral-charge states could form in other insulating magnetic materials as well, which has implications for the interpretation of a wide range of μSR data. Furthermore, analogous to Mu in semiconductors, the study of muon-polaron complexes in magnetic oxides may provide detailed information on the dopant characteristics of interstitial hydrogen, a good understanding of which is crucial for a precise control of charge carriers in such materials. Cr 2 O 3 is of additional interest due to its magnetoelectric properties. Being the first material predicted [10] and measured [11,12] to exhibit an induced linear polarization (magnetization) in response to a magnetic (electric) field, it is widely regarded as the prototypical linear magnetoelectric [13,14] and remains the subject of active research [15], directed primarily at exploiting its magnetoelectric properties for device applications [16][17][18]. In addition, there are unresolved fundamental questions raised by the recent prediction that an electric charge within a linear magnetoelectric is surrounded by a monopolar magnetic field distribution, and thus is subject to a magnetic force in an external magnetic field [19]. μSR is a unique way to investigate such predictions since the spin-polarized muon acts both as a test charge and a sensitive probe of the local magnetic field. However, as indicated by studies from the early days of μSR [20][21][22][23], the spectra in Cr 2 O 3 are complex, and their interpretation was inconclusive. Furthermore, given the weak magnetoelectric coupling in Cr 2 O 3 , only subtle changes to the local magnetic environment in response to the muon charge are expected, and a thorough understanding of the interaction between the implanted muon and its host material is required as a prerequisite for the search for any muoninduced magnetoelectric effects.
The paper proceeds as follows. In Sec. II, we briefly introduce the μSR technique and summarize the experimental conditions. Next, in Sec. III, we report the results of a comprehensive μSR study of Cr 2 O 3 under zero-field (ZF) conditions and in applied magnetic fields. The data are presented in three parts: (1) In ZF, up to three spinprecession frequencies are observed, indicating three distinct muon environments with different internal magnetic fields B int . (2) Weak external fields B ext ð≪ B int Þ split the observed frequencies into multiplets, providing detailed information on the orientation of the internal fields.
(3) Large applied fields (B ext > B int ) corroborate the weak-field results and reveal an additional frequency. Together, the data exhibit a rich variety of dynamic phenomena that we explain in terms of site metastability and muon dynamics (Sec. IV). Most importantly, above approximately 150 K, we observe both highly dynamic muons undergoing locally restricted hopping and muons that remain static in their site. In order to explain this surprising behavior, we turn to DFT to identify candidate muon sites for all three environments and conclude that the coexistence of site-stable and dynamic muons can be explained with the formation of a charge-neutral Jahn-Teller-stabilized muon-polaron complex (Sec. V). Finally, in Sec. VI, we discuss the implications of charge-neutral states in Cr 2 O 3 and its relevance for other magnetic oxides, as well as possible consequences arising from implanting positively charged muons into a linear magnetoelectric.

II. EXPERIMENTAL DETAILS
The μSR experiments we report here are carried out at the Centre for Molecular and Materials Science at TRIUMF (Vancouver, Canada), although initial spectra were taken at the GPS instrument at PSI (Villigen, Switzerland). The zero-and low-magnetic-field measurements are taken in the LAMPF spectrometer, and the high-magnetic-field data are acquired in the NuTime spectrometer. All data are acquired with the initial muon spin polarization P i perpendicular to the beam direction (ẑ).
In a μSR experiment, spin-polarized, positively charged muons are implanted into the sample, where they decay with a lifetime of τ μ ¼ 2.2 μs. The resulting decay positron is emitted preferentially along the muon spin direction and can be detected in a time-resolved manner with plastic scintillators placed in pairs around the sample. This anisotropic positron emission introduces a spin-dependent imbalance in the count rate, causing the asymmetry signal SðtÞ, the countnormalized difference of a counter pair, to be directly proportional to the spin polarization PðtÞ along the detector axis. A detailed description of the μSR technique can be found in Ref. [1]. In the presence of a magnetic field B, the muon spin precesses about the field direction with frequency f ¼ γ μ =2π · jBj, where γ μ ¼ 2π × 135.5 MHz=T is the muon gyromagnetic ratio. Thus the precession frequency is a direct measure of the local magnetic field experienced by the muon.
In a crystal lattice, the charged muon usually stops in one or more distinct sites that minimize the overall energy. At a given temperature, several crystallographically distinct sites may be populated, each of which causes a different time evolution of the spin polarization. For example, in magnetic materials, muons may experience different internal fields at inequivalent sites, causing spin precession at different frequencies. In this case, the observed signal SðtÞ is a sum of several components S i ðtÞ. In this paper, oscillatory signals are fit to exponentially damped cosines where the amplitude A i is a measure of the signal weight, f i is the frequency, ϕ i the initial phase, and λ i the relaxation rate. Nonoscillatory components are parametrized by simple exponentials of the form S i ðtÞ ¼ A i expð−λ i tÞ. All data are fit with the MUSRFIT analysis framework [24]. Several single-crystal specimens sourced from SurfaceNet (Rheine, Germany) are used: a 10 × 10 × 10 mm 3 single crystal (C1) with the c axis in plane and ½1120 out of plane, and 8 × 8 × 0.5 mm 3 (C2) and 5 × 5 × 0.5 mm 3 (C3) single crystals with the c axis out of plane. The zero-field data are taken on C1 with the c axis oriented alongx to coincide with the initial spin direction. Small external fields are applied to C1 (B ext k½1120) and C2 (B ext kc). High-field experiments are carried out on C3 (B ext kc).

III. RESULTS
The primitive unit cell of Cr 2 O 3 is rhombohedral and contains four Cr atoms and six O atoms; see inset 1 in Fig. 1(c). The Cr are arranged in two pairs along the rhombohedral 111 axis (c axis), with the oxygens forming two triangles, rotated 60°with respect to each other, between the Cr pairs. In the absence of magnetic order, the primitive unit cell is inversion symmetric, and there is threefold rotation symmetry around the c axis.
In oxides, muons are generally found to stop approximately 1 Å away from an oxygen, similar to the hydrogen in a hydroxyl OH bond [25,26]. Assuming this holds for Cr 2 O 3 , we can use symmetry arguments to make some general statements about potential muon stopping sites. All six oxygens are crystallographically equivalent; thus, any given muon site close to one oxygen can be projected by either inversion or 120°rotations about c into another equivalent site. Consequently, there are at least six (or integer multiples thereof) electrostatically equivalent stopping sites within the primitive unit cell, which, when projected onto the c plane through the inversion center, form a hexagon; see inset 2 in Fig. 1(c) and Fig. 8.
Due to localized electrons in the Cr 3d shell, there is a magnetic moment associated with each Cr. Below the Néel temperature T N ¼ 307 K, those moments align pairwise opposite to each other along the c axis [see inset 1 in Fig. 1(c)], causing an internal magnetic field B int at the muon stopping sites. The magnetic structure breaks the inversion symmetry, B int ðrÞ ¼ −B int ð−rÞ. As a consequence, the direction of the internal fields associated with the various electrostatically equivalent sites is different. However, as the Cr moments are parallel to the c axis, the magnitude jB int j at each of the sites is the same. Since only the magnitude determines the precession frequency, muons that stop in any one of the equivalent sites in zero external field precess with the same frequency and contribute to the same signal S i ðtÞ. From now on, we refer to an ensemble of electrostatically equivalent sites that have the same jB int j as a muon environment. Note that at a given temperature, muons may stop in different environments with distinct jB int j.
We start with a presentation of the ZF results. Then, the effects of external magnetic fields B ext are described, first for fields small compared to the internal field (B ext ≪ B int ), then for large fields (B ext > B int ).

A. Zero external field
A ZF μSR spectrum taken at T ¼ 2. E1 − E3 in order of appearance: The frequency observed up to T N is called f ZF E1 , whereas f ZF E2 can be seen only up to 190 K and f ZF E3 up to 60 K. All three precession frequencies increase with decreasing temperature, approximately tracking the sublattice magnetization.
The spectra contain both oscillating and nonoscillating components and are fit with up to three damped cosines, Eq. (1), a nonrelaxing component, and a relaxing component (nonzero only above 160 K).
The fit results for the oscillatory components associated with E1 − E3 are shown in Fig. 2. The amplitude A E1 is constant up to 200 K, above which it increases and approximately doubles at T N . Both A E2 and A E3 are approximately constant. The relaxation rates λ E2 and λ E3 increase sharply when approaching the temperature where their associated ZF frequency vanishes. While both phases ϕ E2 and ϕ E3 , can be considered constant, there is a pronounced peak in ϕ E1 between 200 K and T N ; see Fig. 2(c). As we discuss in Sec. IV B, such a change in phase is indicative of a transition from another (so far unspecified) environment into E1, a hypothesis supported by the increase of A E1 at the same temperature. Aside from the precession signals, there is a sizable nonoscillatory relaxing component that appears above approximately 160 K (not shown). This is attributed to the E3 Ã component that we discuss below in Sec. III B 2.

B. Weak external fields
There are two main effects caused by weak external magnetic fields (B ext ≪ B int Þ, (1) the degeneracy of jB int j for electrostatically equivalent stopping sites within one environment is lifted, and (2) a component precessing in B ext rather than B int appears.

Orientation of the internal magnetic field
The internal field direction at a stopping site can be described by two angles; θ is defined as the smallest angle enclosed by B int and the c plane, and φ as the azimuthal angle (enclosed by ½1120 and the B int projection onto the c plane); see insets in Fig. 1. From symmetry, stopping sites forming a given environment can be projected onto the c plane to form a hexagon. The six φ values of an environment are given by δ þ 0°; δ AE 60°; δ AE 120°and δ þ 180°, with δ being the smallest angle enclosed by ½1120 and a hexagon corner. As we note above, the internal field magnitude at the stopping sites forming a given environment is the same, but its direction is not. While this is inconsequential in ZF, the relative orientations of B int and B ext matter in the presence of external fields, where the precession frequency is determined by the magnitude of the vector sum jB int þ B ext j. Consequently, the application of B ext lifts the degeneracy of the precession frequencies within an environment, causing multiplet splittings, which are schematically illustrated in Figs. 3(a)-3(c) both for B ext jjc and B ext ⊥c. The FTs of the μSR spectra at low temperatures are shown in Fig. 3 for (d) B ext ¼ 30 mTkc [27] and (e) B ext ¼ 20 mTk½1120⊥c. Comparison with the ZF spectrum [ Fig. 1(b)] shows that for B ext kc, the E1 − E3 lines split into doublets, while for B ext ⊥c, more complex multiplets are observed. The spectra are fit with up to 12 oscillatory signals [Eq. (1)] and a small nonoscillating signal. The obtained frequencies (f exp ) are shown in Tables III and IV. Under the assumption that B ext does not induce changes of B int , all multiplet frequencies can be consistently described by the vector sum jB int þ B ext j, which allows extraction of the θ and δ values describing the orientation of B int in E1 − E3; see Table I   Appendix A for details. The obtained angles provide stringent criteria for comparison with the internal field of candidate muon sites calculated with DFT; see Sec. V.

Evidence for a signal component with zero internal field
Having discussed the effect of B ext on the precession frequencies, we now turn our attention to the nonprecessing component that appears in ZF above approximately 160 K. At coinciding temperatures and in both B ext ⊥c and B ext kc (not shown), there is a component of comparable amplitude that oscillates at the Larmor frequency of the external field (f ext ¼ γ μ =2πjB ext j), which is absent below 150 K; see Fig. 4. Spin precession about B ext rather than B int , despite ordered Cr moments, indicates that the muons giving rise to this signal are not subject to an internal field. The temperature dependence of the amplitude and relaxation rate of this signal, termed E3 Ã in anticipation of its interpretation in Sec. IV B, is shown in Fig. 2.

C. Large external fields
The ZF spectra indicate that the highest B int is about 0.75 T. Here, we apply external fields significantly higher than this. The temperature dependence of the FTs of μSR spectra taken in B ext ¼ 4Tjjc is shown in Fig. 5. Again, a doublet splitting is expected for each ZF frequency; however, since the positions of f AE with respect to f ZF depend on the relative strength of B ext to B int , and The two outer frequencies (colored in black) compose the E1 doublet f AE E1 , while the second (orange) and third (blue) highest frequencies correspond to f þ E2 and f þ E3 , respectively. The remaining line (uncolored) is a superposition of both f − E2 and f − E3 , explaining its large amplitude. The temperature evolution follows mostly what is expected from the ZF results. With increasing temperature, the E1 doublet splitting decreases as jB int j decreases. Above 200 K, its amplitude becomes larger, and the line broadens approaching T N . Likewise, f þ E2 follows the decreasing B int and disappears above approximately 185 K. For all temperatures where f þ E2 is observed, the uncolored line has a contribution from f − E2 . Above approximately 170 K, a large component close to f ext appears, and is, in accordance with Sec. III B 2, assigned to E3 Ã . The E3 doublet is observable only below 50 K.
Remarkably, there is an additional component between 50 and 170 K, which is termed E3 0 [Figs. 5(c)-5(e), green]. As we discuss in Sec. IV B, this additional line is strongly indicative of local hopping between adjacent E3 sites.
Details on the data analysis as well as fit results for the amplitudes can be found in Appendix B. From the multiplet splittings, θ values matching closely those obtained in low field are extracted; see Table V. This agreement indicates that even in large B ext , B int is not significantly affected (reasonable since the Cr Zeeman energy in 4 T is much smaller than the exchange coupling [28]), and the precession frequencies are well determined by vector addition. The fitted frequencies are shown in Fig. 6. The red lines represent calculated doublet frequencies f AE E1 ðTÞ and f AE E2 ðTÞ assuming constant θ; see Appendix B for details.
There is good overall agreement with the data, indicating that θ is largely temperature independent.

IV. EXPERIMENTAL EVIDENCE FOR SITE METASTABILITY AND DYNAMICS
Three ZF frequencies are observed [see Fig. 1(c)] and attributed to three distinct muon environments, E1 − E3. Each environment contains a number of electrostatically equivalent sites with the same magnitude but different directions of B int . Above approximately 160 K, a component E3 Ã precessing in the external rather than the internal field is observed both in low and high field, indicating an environment characterized by zero internal field.
The E2 and E3 signals disappear at different temperatures, while E1 is observed over the complete temperature range, indicating that each environment has a distinct potential energy. At low temperatures, E1 − E3 are all populated. Since site populations are determined by the epithermal implantation process rather than thermodynamic equilibrium, it is possible that a muon occupies metastable sites with higher energy than the ground state. If thermally activated transitions to a lower energy state are inaccessible within its short lifetime, the muon may remain in the metastable site and give rise to a distinct signal [29,30]. However, with increasing temperature, site changes either within one or into another environment may become possible.
Around 180 K, the E1, E2, and E3 Ã signals are observed and account for the full signal; see Fig. 2(a). Noting that the amplitudes of both E2 and E3 Ã are approximately temperature independent in the respective regions where they are observed, we conclude that (1) muons in E2 do not transition into E3 Ã ; (2) thus, the disappearance of E2 stems from a transition into E1 at sufficiently high temperatures, and (3) consequently, by conservation of total amplitude, muons that stop in E3 below 50 K must give rise to E3 Ã at higher temperatures.  In this section, we first present a model supporting the E2 → E1 transition. Then, we discuss the evolution of muons from E3 into E3 Ã in terms of local muon hopping between adjacent electrostatically equivalent E3 sites.

A. E2 → E1 transition
Here, we show that the disappearance of the E2 signal around 200 K and the subsequent increase in E1 amplitude is consistent with a metastable E2 environment that allows for transitions into E1. The following discussion is based on two assumptions: (1) The E2 → E1 transition can be described by a thermally activated, exponential rate of the form ΛðTÞ ¼ ν 0 expð−E a =k B TÞ, where E a and ν 0 are the activation energy and attempt frequency. (2) At the time of implantation, the probability for the muon to initially occupy a site in any of the three environments is temperature independent (i.e., at all temperatures, the same fraction of muons start out in E1, E2, and E3). We discuss this assumption in Sec. VI. While the initial fraction starting in E2 is independent of temperature, the actual time spent in this environment depends on the transition rate ΛðTÞ. If ΛðTÞ is much smaller than f E2 , the E2 muons precess for many periods with f E2 before transitioning, and oscillatory signals from both E1 and E2 with amplitudes A E1 and A E2 can be detected. In contrast, if ΛðTÞ is much larger than f E2 , the E2 muons change to E1 before the muon spin has a chance to precess with f E2 , and a single oscillatory signal at f E1 with a combined amplitude A ¼ A E1 þ A E2 can be observed. However, if ΛðTÞ is comparable to f E2 , the E2 muons may precess at f E2 prior to the transition and acquire a phase shift with respect to muons initially in E1. This phase shift causes a smaller apparent E1 amplitude A and an overall phase shift ΦðTÞ. In Appendix C, a model accounting for such a transition and expressions for AðTÞ and ΦðTÞ are described. In Fig. 7, the E1 phase and amplitude data of both the ZF signal and the high-field doublet are compared with the transition model, with activation energy E a ¼ 180 meV and attempt frequency ν 0 ¼ 8 × 10 11 Hz being shared parameters for the complete dataset.
There is excellent qualitative agreement between the model and the data; the peak in the phase, including the opposite direction for the high-field doublet, and the increase in amplitude are well described with the same model parameters. Thus, the proposed E2 → E1 transition with a barrier E a ¼ 180 AE 40 meV provides a consistent explanation for the disappearance of the E2 signal, its associated increase in relaxation rate, and the subsequent increase in E1 signal amplitude.

B. Local hopping
Next, we show that both the E3 frequency observed below 50 K and the E3 Ã signal precessing at the Larmor frequency of the external field arise from muons in the same environment. The appearance of E3 Ã above approximately 160 K [see Figs. 4 and 5] indicates that a fraction of muons experience no internal field. This is surprising since the simultaneous observation of the E1 signal clearly shows the presence of ordered Cr magnetic moments. Although there are high-symmetry sites along the c axis where the internal field precisely cancels, these sites are far (≫ 1 Å) from an oxygen and energetically unfavorable for an interstitial μ þ site, as we confirm by DFT in Sec. V. Instead, we consider in Fig. 8 an example configuration of an environment comprised of six electrostatically equivalent muon stopping sites (white spheres) based solely on symmetry considerations (see beginning of Sec. III) and a muon-oxygen distance of 1 Å. Given the close proximity of nearby electrostatically equivalent sites, thermally activated local hopping between adjacent sites seems plausible. For sufficiently fast intraenvironment hopping, the effective internal field experienced by the dynamic muon is the average over all sites. Noting B int ðrÞ ¼ −B int ð−rÞ and the symmetry of equivalent sites, it becomes clear that this average is zero. We hypothesize that muons stopping in E3 sites undergo such local hopping around a single hexagon at elevated temperatures, leading to the disappearance of the f ZF E3 signal and the subsequent observation of the E3 Ã signal in an external field.
This consistently explains the observed data: At low temperatures, muons stopping in E3 sites are quasistatic, i.e., no or only very slow hopping occurs, and each muon precesses predominately in the internal field of one site, giving rise to f ZF E3 . Above approximately 160 K, a relaxing nonoscillatory component appears in ZF, consistent with a fraction of muons that are not subject to any field. In the intermediate temperature region 60-160 K, no signal is observed, as the hopping is neither fast enough to efficiently average out the internal field nor slow enough to allow for the observation of coherent E3 oscillations. External fields, both small and large, cause multiplet splittings of f E3 at low temperature, consistent with muons being quasistatic, while above approximately 160 K, a signal precessing in B ext can be observed (E3 Ã ), since B int is averaged to zero and does not contribute to the field magnitude. In high field, an additional signal, E3 0 [see Fig. 5], is observed in the intermediate temperature region. A Monte Carlo simulation of the muon depolarization function in 4Tkc, assuming local hopping, identifies E3 0 as the average of the doublet frequencies f av and allows for an estimate of the energy barrier between E3 sites E b ¼ 42 AE 5 meV; see Appendix D for details on the simulation and an in-depth discussion. Overall, we clearly show that the E3=E3 0 =E3 Ã signals arise from different dynamic regimes of muons undergoing local hopping in a single environment.
Already in the first μSR paper on antiferromagnets, local diffusion between electrostatically equivalent sites was considered a possibility in Fe 2 O 3 [31]. Subsequently, local motion was speculated to occur in Cr 2 O 3 [23], and is suspected [25,32] and observed [33][34][35] in a range of materials. Here, we conclusively show that local hopping indeed occurs in Cr 2 O 3 by direct observation, identification, and consistent description of the distinct signals that arise as a result of restricted motion in a system with broken magnetic inversion symmetry, a hop rate changing several orders of magnitude over the observed temperature range, and various applied fields.
We further note that muons in both E1 and E2 are, apart from the E2 → E1 transition, site stable; i.e., no intraenvironment motion occurs. This is evident from the pronounced multiplet splitting that is observed at all temperatures where E1 and E2 signals are detected; see Figs. 4 and 5.

V. IDENTIFICATION OF MUON STOPPING SITES WITH DFT
Thus far, using simple models describing a thermally activated E2 → E1 transition and local hopping within the E3 environment, and without explicit knowledge of the stopping sites, we have explained the major features in the data. The coexistence of site-stable and highly mobile muons is intriguing, especially since there is no evidence for interexchange between dynamic E3 muons and static E1 or E2 muons, even in the presence of the E2 → E1 transition. To gain deeper insight into this surprising behavior, we turn to DFT to identify muon stopping sites. With the recent increase in availability and capability of computing resources countering the large computational demands of first-principles calculations, DFT has had great success in providing information about location and stability of muon stopping sites in a range of materials and is developing into an important new tool for μSR (see Ref. [36] for a review, and Refs. [37,38] for recent developments).
Since the inception of the μSR technique, knowledge of the location of the muon within the sample was of key importance. The main motivation for the early μSR studies of antiferromagnets, prompted by the first observation of ZF μSR signals in Fe 2 O 3 [31] and its isomorph Cr 2 O 3 [20], was to establish the muon as a sensitive and useful probe of the local magnetic properties of the host material by determining (1) where the muon stops and what its dynamic properties are with respect to site stability and diffusion and (2) if and under what conditions muonium is formed in insulating (anti)ferromagnets. Based on simple electrostatic considerations, two sets of possible stopping sites were found for the corundum structure, so-called Rodriguez (R) sites [31] located in the Cr gap close to the inversion center [see Fig. 8], and Bates (B) sites [22] in (B0), or slightly above and below (B1) the oxygen basal plane; see Fig. 2 in Ref. [23]. The internal magnetic field in these sites was estimated by summing over the dipolar contributions from surrounding Cr moments. Additionally, covalency effects were considered, and attempts to assign ZF frequencies to specific sites yielded some partial and approximate agreements, although the overall results remained inconclusive FIG. 8. Example configuration of electrostatically equivalent muon sites based solely on symmetry considerations and the constraint that muon stopping sites (white spheres) are located 1 Å away from an oxygen (red spheres). The close proximity of nearby sites suggests localized hopping. The electrostatic potential isosurface of the undistorted crystal structure is shown in blue, with the red patches indicating local electrostatic minima corresponding to the position of the so-called Rodriguez sites [31]. [22,23]. No evidence for Mu or a neutral-charge state was identified.
Here we calculate the muon stopping sites in Cr 2 O 3 using the Vienna ab initio simulation package [39][40][41]; see Appendix E for details. The positive muon is modeled as a hydrogen nucleus embedded within an 80-atom 2 × 2 × 2 rhombohedral supercell (SC) of Cr 2 O 3 . Two muon charge states are considered (1) the bare, positive muon, with a uniform charge background ensuring overall charge neutrality and (2) a neutral muon state allowing for the extra electron. As a first step, promising initial muon positions are identified. Given the muon's tendency to form a muon─O bond of length approximately 1 Å, a set of 99 initial configurations are generated with the muon positions equidistributed on a 1-Å sphere centered on one oxygen. Since all oxygen atoms are electrostatically equivalent, it is sufficient to perform this search on a single oxygen. Initial static calculations of the Hellmann-Feynman forces allow us to discard sites with forces larger than 10 eV=Å. For the remaining sites, the Cr 2 O 3 ions are relaxed while keeping the muon fixed. This initial step of relaxing only the lattice allows for a search of self-trapped metastable sites. Finally, both the muon and lattice are fully relaxed until the Hellmann-Feynman forces are below 5 meV=Å. Additionally, muons starting out in the R and B sites and in the unit-cell center, both the exact center (C E ) and slightly offset along the c axis (C O ), are considered. The structure files of all candidate stopping sites for both charge states are available in Ref. [42].
The total hyperfine field B tot at each site r μ has (1) a dipolar contribution B dip mainly from the Cr 3d electrons, and (2) a Fermi contact term B c from unpaired spin density ρ s ðr μ Þ at the muon stopping site. B dip is calculated by embedding the distorted 2 × 2 × 2 SC in a superstructure of undistorted SCs and summing over the dipolar contribution from the spin-density grid points. Despite a fine grid spacing of 0.055 Å, the finite grid causes artifacts for points in close proximity to r μ . This is mitigated by excluding grid points less than R ¼ 0.5 Å away from r μ . B c is calculated by [37] ) and contact (f c ) contributions, the vector sum of which is given as f tot , and angles θ and φ as defined in III B 1. ΔE is the energy relative to the ground state of each charge state, respectively, and d z 2 the spin state of the extra electron for D1 0 -D5 0 . Bold rows indicate candidates for E1 − E3.  where μ 0 is the vacuum permeability and μ B the Bohr magneton. ρ s ðr μ Þ is approximated by projecting the spin density within an R ¼ 0.5 Å sphere onto an s-wave state at r μ . Results for both charge states (positive and neutral) are shown in Table II. Calculated fields are given in units of frequency f i ¼ ðγ μ =2πÞjB i j. Note that f tot is obtained by vector addition jB dip þ B cĉ j. The energies ΔE are given with respect to the ground state and are comparable only within a given charge state. The stated uncertainties are estimated by varying the sphere radius R for both the contact and dipole term calculations in the range 0.3-0.7 Å.
For the positive-charge state (superscripted +), the B1 þ site cannot be stabilized, and muons placed in the primitive unit-cell center (C sites) can be discounted as viable muon stopping sites based on the magnitude and direction of B int and large ΔE. For the same reasons, B0 þ is unlikely to represent a muon stopping site. Unless purposefully placed in a B þ or C þ site, muons relax into a position close to but distinct from the R sites (the electrostatic minima of the undistorted lattice). We name this site D due to the doughnut-shaped potential energy surface formed by electrostatically equivalent sites. The difference between D and R arises predominantly from the muon-induced lattice distortion, which was not accounted for previously. The close proximity of adjacent D þ sites makes them excellent candidates for the E3 environment. We discuss possible explanations for the discrepancy of 22% between f ZF E3 ¼ 76.7 MHz and the calculated value below. While DFT calculations considering the positive muon can account for E3 muons undergoing local hopping, E1 and E2 are thus far unexplained, motivating a search for charge-neutral muon states. Results are shown in the bottom half of Table II. B0 0 , B1 0 , and C 0 E can be dismissed based on large ΔE, and C 0 O based on the large f tot . This leaves five variations of the D 0 site labeled D1 0 -D5 0 . Comparison with Table IV and consideration of the calculated energies suggest that D1 0 and D2 0 are candidates for E1 and E2, respectively: The measured and calculated frequencies of E1 (þ14%) and E2 (þ12%) and θ agree reasonably well, and ΔφðE2 − E1Þ is close to the expected 30°; compare Table I. Additionally, consistent with the proposed E2 → E1 transition, D2 0 has a larger energy than D1 0 . We discuss the differences between D1 0 -D5 0 and further aspects of the D2 0 → D1 0 transition below. Figure 9 shows the positions of the E1 − E3 candidate stopping sites D1 0 , D2 0 , and D þ .
A more detailed analysis of the charge-neutral D 0 states reveals that the extra electron localizes predominately on a nearby Cr, where it changes the valence from Cr 3þ to Cr 2þ and singly occupies the initially orbitally degenerate e g orbitals. This causes a Jahn-Teller distortion [7] by further elongating the Cr─O bond of the oxygen the muon is bound to from 2.18 Å for the bare muon to 2.43 Å in the chargeneutral case; see Fig. 10. The subsequent lowering of energy [see Fig. 10] stabilizes this charge-neutral complex of a negatively charged JT polaron [8,9] and positive muon. This proposed mechanism is supported by the extra electron occupying the lowered e g level d z 2 ; see yellow isosurface representing the charge density of the topmost occupied band in Fig. 10. We note the similarity to the paramagnetic Ti-O-Mu complex recently observed in (nonmagnetic) TiO 2 [43,44], where an unpaired electron sits on a nearby Ti atom, and the oxygen-bound muon forms a complex with the resulting small polaron. Crucially, however, the muonpolaron complex that we report on here is not paramagnetic and therefore distinct from Mu, since the bound electron is strongly coupled to the 3d electrons of the Cr host ion. We discuss this in detail in Sec. VI.
The position. A transition between D 0 states is mainly characterized by a change in position of the extra electron, rather than the muon. Note that the spin of the extra electron, being coupled to the 3d electrons of its Cr host, may be different for different D 0 states, as indicated in Table II. We propose that the higher energy states D3 0 -D5 0 are not occupied since they can easily transition into either D1 0 or D2 0 depending on their spin. However, going from D2 0 into the ground state D1 0 requires an electron spin flip (assuming the muon stays stationary); i.e., an additional energy barrier has to be overcome. While the precise process for the E2 → E1 transition is still under investigation, we tentatively attribute the observation of the metastable E2 environment to the existence of such a spin barrier.
In general, the energy barrier to move the joint muonpolaron complex is expected to be significantly larger than for the bare muon [45], providing a compelling explanation for the stability of E1 and E2, and the coexistence of sitestable muons and highly mobile muons in E3.
Note that DFT local field predictions depend on the approximation of the exchange-correlation functional and value of the U eff parameter. By comparing the LDA, PBEsol, and SCAN functionals with a reasonable range of U eff corrections, we find variations of 15% in the predicted frequency magnitudes, 60% in the θ angle and 10% in the φ angle with respect to the LDA þ U eff ¼ 4 eV values; see Ref. [42].
Finally, we note that an accurate consideration of the zero-point motion of the muon is necessary to calculate the formation energies of the different charge states [46], and the energy barriers for intra-E3 hopping and the E2 → E1 transition. While a full treatment of the quantum nature of the muon is beyond the scope of the present paper, work is currently ongoing and will be published separately.

VI. DISCUSSION
Paramagnetic Mu centers are expected to be subject to fast relaxation in magnetic materials [3]. In the previous sections, we presented strong evidence for the formation of a charge-neutral muon-polaron complex in Cr 2 O 3 that, while not exhibiting signatures conventionally expected from neutral-charge states, significantly influences the muon behavior and contributes a well-resolved signal. In particular, rather than giving rise to a well-defined spectrum of typically two or four precession frequencies that are determined by a spin Hamiltonian involving the muon and electron Zeeman energies and a muon-electron hyperfine interaction [1], the precession signal for the muon-polaron complex in Cr 2 O 3 consists of a single frequency much like the normal positive-charge state, i.e., the bare μ þ with no additional electron nearby. The reason for the different behavior compared to that of Mu is that the muon-polaron complex is not paramagnetic, since its bound electron is strongly coupled to the 3d electrons of the Cr host ion, which themselves are antiferromagnetically coupled to the ordered network of magnetic ions of the host. We note that the bound electron is not centered on the muon but localizes on a nearby Cr. This is not unique to the muon-polaron complex, since similar situations are found for the paramagnetic Mu complex in TiO 2 [4,43,44] and ZrO 2 [47], and bond-centered Mu in silicon [48].
The relevance of the discovery of a muon-polaron complex in Cr 2 O 3 may extend to other magnetic oxides, such as CuO [32,49], Fe 2 O 3 and FeTiO 3 [23], Fe 3 O 4 [50], CaMnO 3 [51] and the orthoferrites [26], which all show multiple zero-field precession frequencies that, while attributed to metastable sites, are not conclusively explained. In general, our study suggests that neutralcharge states and their potential impact on crystal-field levels have to be carefully considered in all insulating magnetic materials, in particular in transition-metal oxides where multiple oxidation states for the metal ion are possible. Detailed DFT calculations may be required to separate intrinsic magnetic properties from muon-induced effects; however, we note that since the muon may occupy metastable states, the formation energies alone are not sufficient to reliably predict the muon charge states in a given material, but transition barriers between charge states have to be taken into account as well.
We emphasize that the present muon-polaron complex is different from the controversial concept of a muon-induced magnetic polaron proposed by Storchak et al. [52,53] and disputed by others [54,55], which invokes a localized electron bound to a muon mediating a ferromagnetic coupling between neighboring magnetic ions, resulting in a "ferromagnetic droplet" characterized by a gigantic local spin.
We speculate that muon-polaron complex formation in Cr 2 O 3 occurs by a similar mechanism to Mu formation in semiconductors [5]. Upon implantation, the muon slows down by creating electron-hole pairs. Toward the end of its ionization track, it may capture an electron and subsequently form a charge-neutral complex. The implantation and electron capture processes are epithermal, and thus independent of sample temperature and thermodynamic equilibrium, providing justification for the previously stated hypothesis that E1 − E3 are populated with the same ratio at all temperatures. Muons may stop and selftrap in metastable (charge) states with energies larger than the ground state and deexcite only if thermally activated transitions, e.g., from E2 to E1, are accessible during the muon lifetime [29,30].
Paramagnetic Mu has been used extensively to investigate the dopant characteristics of hydrogen in a wide range of semiconductors including oxides [2][3][4][5]. This is because the electronic structure of Mu in a solid is virtually identical to that of hydrogen, aside from small differences caused by the larger zero-point motion due to the lighter muon mass. We propose that with the observation of a neutral-charge state in Cr 2 O 3 , μSR shows its ability to investigate the behavior of interstitial hydrogen in magnetic oxides as well [56]. A thorough understanding of unintentional hydrogen doping in such materials is crucial, since a wide range of technologically relevant fields such as dilute magnetic semiconductors for spintronics [57][58][59][60] and superconductivity depend on a precise control of charge carriers in magnetic oxides, or, in the case of multiferroics, require low-leakage thin films [61].
Interestingly, interstitial hydrogen is predicted to form a shallow donor state in Cr 2 O 3 , alongside a range of other materials including the aforementioned CuO, Fe 2 O 3 , and FeTiO 3 [62]. The observation of E1 up to T N suggests that the muon-polaron complex stays intact up to room temperature, indicating that hydrogen is instead a deep impurity. This is in contrast to an ongoing study of muon-polaron complexes in Fe 2 O 3 (to be published separately), which suggests complex ionization above 200 K, indicating shallow donor behavior. We note that Cr 2þ with 3d 4 high spin is strongly JT active, while Fe 2þ with 3d 6 high spin is weakly JT active and hypothesize that donor characteristics are determined, at least in part, by the strength of the JT effect. The role of magnetic interactions in the stabilization of the muon-polaron complex remains an open question.
The prediction of hydrogen-induced n-type conductivity in ZnO [46] and the subsequent experimental observation of the corresponding shallow Mu donor state [63] prompted a search for criteria to predict dopant behavior and charge state of interstitial hydrogen and muons, leading to generalized principles for elemental and binary semiconductors [60] and oxides [4,62,64,65]. We hope that our discovery of a JT-stabilized muon-polaron complex stimulates research to extend those principles to explicitly account for polaronic, and if required, magnetic contributions. Such adjusted criteria could provide valuable guidelines on whether or not neutral-charge states are expected in μSR experiments on insulating magnetic materials.
Lastly, having acquired a detailed understanding of how the muon interacts with Cr 2 O 3 , we address possible implications arising from implanting a point charge into a linear magnetoelectric (ME) material. Khomskii predicts that a point charge inside an isotropic linear ME is surrounded by a monopolelike magnetic texture and subject to a force in applied magnetic fields [19]. The monopolar magnetic field distribution is only expected sufficiently far from the charge where the bulk approximation holds, while no predictions are made for the immediate core region. The case of a muon inside Cr 2 O 3 is more complex since (1) the ME coupling is not isotropic [13], (2) some muons are bound in a muon-polaron complex, an effectively chargeneutral entity, and (3) the muon induces significant lattice distortions in its immediate vicinity and experiences a magnetic field dominated by the core region which cannot be treated in the bulk limit. Thus, the muon may sense a change in its magnetic environment in response to the electric field of its charge arising from changes in the position of magnetic ions and canting of magnetic moments in its immediate environment, i.e., a local ME effect; however, such an effect may not necessarily follow the bulk ME coupling. Results from preliminary noncollinear DFT calculations suggest only a small spin canting <0.3°i n the immediate neighborhood of the muon; however, the elevated ME coupling at higher temperatures is not yet taken into account. We note that while the ME response outside the core region is expected to have monopolar contributions only in linear ME materials, muon-induced local ME effects may occur in non-ME compounds as well. Finally, we comment on the prediction that a static charge inside a ME is subject to a force in a magnetic field [19]. Using a lower limit of 1 eV=Å 2 for the force constant of the potential experienced by the muon bound to an oxygen, the change in position in response to a field-induced force on the muon can be estimated to be smaller than 10 −5 Å in 4 T; i.e., it is negligible in the context of this experiment. Despite the challenges discussed above, the investigation of local ME effects induced by the muon in its duality as a test charge and sensitive probe for magnetism remains a fascinating area of research, and we hope that the present paper initiates studies both experimental and theoretical in this direction.

VII. CONCLUSIONS
In summary, we carried out a comprehensive μSR study of Cr 2 O 3 under zero-field conditions and in applied magnetic fields. In zero field, we observe three spinprecession frequencies attributed to three distinct muon environments E1 − E3 with different internal magnetic fields. Small applied magnetic fields along various symmetry directions split the observed frequencies into multiplets, providing detailed information on the orientation of the internal fields. The temperature dependence reveals a rich dynamic behavior that we explain in terms of a thermally activated transition between E2 and E1, and intra-E3 local muon hopping. Notably, we observe coexistence of highly dynamic E3 muons and site-stable muons in E1 and E2. Muon stopping sites and charge states for all three environments are determined using DFT, and the coexistence is explained by the formation of a chargeneutral, JT-stabilized muon-polaron complex. The identification of such a charge-neutral complex in the antiferromagnet Cr 2 O 3 has implications for other magnetic oxides, since the formation of muon-polaron complexes can significantly influence the stability and location of stopping sites, but its existence may be "hidden" since the behavior conventionally associated with neutral-charge states is not displayed. Furthermore, this discovery opens up a route to study interstitial hydrogen in magnetic oxides, where precise control of the carrier density may be critical for device functionalities. Given the technological importance of magnetic oxides, we hope this study stimulates a search for generalized principles describing the dopant behavior of hydrogen that account for polaronic, and if required, magnetic contributions. Crystal structures were drawn using VESTA [66].

APPENDIX A: DETAILS ON INTERNAL MAGNETIC FIELD ORIENTATION
Under the assumption that B ext does not induce changes of B int , the resulting frequency multiplets can be calculated by simple vector addition. For B ext jjc, only θ is relevant, and we expect the ZF frequency f ZF to be split into a doublet [see Fig. 3(a)] with the frequencies f AE given by where f ext ¼ γ μ =ð2πÞ · jB ext j. Given the equal number of sites with þθ and −θ, equal amplitudes for both doublet lines are expected (for B ext ≪ B int ). For B ext jj½1120⊥c, the multiplet splitting determined by both φ and θ is given by fðθ; φÞ ¼ ½½f ZF sinðθÞ 2 þ ½f ZF sinðφÞ cosðθÞ 2 þ ½f ext þ f ZF cosðφÞ cosðθÞ 2 1=2 : Only for δ ¼ 0 or AE30°, B ext causes a multiplet splitting with fewer than six lines but yields either a quadruplet for δ ¼ 0 or a triplet for δ ¼ AE30°, with an amplitude ratio of 1∶2∶2∶1 and 2∶2∶2, respectively; see Figs. 3(b) and 3(c). For B ext kc, Eq. (A1) is used to calculate the expected doublet frequencies (f calc ) with f ZF from Fig. 1(c). Minimizing jf calc ðθÞ − f exp ðθÞj þ jf calc ð−θÞ − f exp ð−θÞj yields θ values that produce excellent agreement between f calc and f exp for E1 − E3; see Table III.
Using those θ values, Eq. (A2) is used to investigate the B ext k½1120⊥c multiplet splittings. For E1 and E2, the measured frequencies are in very good agreement with the calculated values f calc for δ ¼ 0°and 30°, respectively; see Table IV. Furthermore, the amplitudes are close to the predicted 1∶2∶2∶1 and 2∶2∶2 ratios. For E3, δ ¼ AE17.5°y ields a reasonable agreement between f calc and f exp , assuming that the outer two lines on either side of the resulting sextet are not resolved but appear at their average frequency (shown in bold in Table IV) with twice the amplitude. From that, an amplitude ratio of 2∶1∶1∶2 is expected, which is indeed observed. Considering crystal alignment and the observed frequencies and amplitude overlap, six components are considered; all three doublets share amplitude, and the overlapping line is fit to two components which have the same phase and frequency, but share the amplitude and relaxation rate with f þ E2 and f þ E3 , respectively. At and above 50 K, the data are fit to up to five exponentially damped oscillatory functions, with shared amplitudes for the E1 and E2 doublets but separate relaxation rates.
Additionally, for the spectra at T ¼ 50; 110-155 K, the amplitude of the E3 0 component [see Fig. 5] is fixed to 0.0365, a value obtained by averaging over the amplitude of surrounding temperature points. This constraint is necessary in order to get a meaningful measure of the relaxation rate of this strongly damped component.
The fitted amplitudes of E1 and E2 are shown in Fig. 11(a). The temperature dependence tracks largely the ZF behavior: E2 is constant in amplitude, and E1 starts to increase above 200 K. The displayed values are the amplitudes of one of the doublets; the total amplitude of muons in either E1 or E2 is twice that. The remaining amplitudes associated with E3, E3 0 , and E3 Ã are displayed in Fig. 11(b). The points with a thick horizontal bar indicate temperature points where the amplitude is constrained. The E3 doublet amplitude accounts for only one of the doublets; the dashed line indicates the total (double) value for reference.
The frequencies obtained from fitting to five exponentially damped oscillatory components are shown in Table V.
Using Eq. (A1), f calc is calculated for the θ value that minimizes jf calc ðθÞ − f exp ðθÞj þ jf calc ð−θÞ − f exp ð−θÞj, yielding very good agreement with the data, and θ values matching closely those obtained in low field; see Table III.
The temperature dependence of f AE E1 and f AE E2 is modeled by Eq. (A1) with θ from Table V and an interpolated temperature dependence of f ZF E1 ðTÞ in Fig. 1(c). f ZF E2 ðTÞ is further approximated by scaling f ZF E1 ðTÞ by the lowtemperature frequency ratio f ZF E2 =f ZF E1 j 2.2 K ¼ 1.49. The red lines in Fig. 6 represent the calculated doublet frequencies f AE E1 ðTÞ and f AE E2 ðTÞ. There is good agreement with the data, indicating that θ is largely temperature independent.

APPENDIX C: DETAILS ON THE E2 → E1 TRANSITION
Here we describe a model for the E2 → E1 transition assuming a thermally activated, exponential rate of the form ΛðTÞ ¼ ν 0 expð−E a =k B TÞ, where E a and ν 0 are activation energy and attempt frequency. The following expression describes the observable signal precessing at f E1 (compare Refs. [68,69]): Muons starting out in E1 are described by the first term, whereas the second describes the E2 → E1 transition taking into account the phase acquired while evolving in E2. For t ≫ Λ −1 , the resultant combined amplitude A and phase Φ can be expressed as where The expressions above assume that B int in E1 and E2 are parallel and perpendicular to the initial spin polarization P i . In order to compare this model to data taken in both ZF and large B ext , small modifications, outlined below, are necessary.
In general, B is not perpendicular to the initial polarization P i , causing the component of BjjP i to act as a holding field. The polarization signal can be decomposed into oscillating and nonoscillating components [1] SðtÞ ∝ cosðθÞ 2 cosðγ μ jBjtÞ þ sinðθÞ 2 ; where 90°− θ is the angle enclosed by B and P i . Zero field.-Here, P i jjc, and for both E1 and E2, the internal field B int encloses an angle θ with the c plane; see Table III. Thus, B int · P i ≠ 0, resulting in a nonoscillatory signal component. The oscillatory f E1 signal amplitude represents only a fraction of cosðθ E1 ¼ 24°Þ 2 ¼ 0.83 of muons in E1, while for E2, cosðθ E2 ¼ 6°Þ 2 ¼ 0.99, and virtually the complete E2 component is oscillating. Thus, if the complete E2 component with observed amplitude A E2 coherently transfers to E1, the transferred amplitude observed at f E1 is only ½cosðθ E1 Þ= cosðθ E2 Þ 2 A E2 . Strictly, this is only valid for Λ −1 ≪ 1=f E2 ; furthermore, a change Δφ ¼ 30°of the internal field direction upon transition is neglected. A calculation addressing both issues was carried out and yielded slight improvements but no major deviations from the simple model and was not included for clarity. The temperature dependence of f ZF E1 ðTÞ and f ZF E2 ðTÞ is obtained by interpolating f ZF E1 shown in Fig. 1(c). High field.-As B ext ⊥P i > B int , the resultant internal field is in good approximation perpendicular to P i . B ext causes a frequency splitting; see Sec. III C. The doublet frequencies (Fig. 6, red lines) are given by Eq. (A1). Here, only transitions without sign change of θ, i.e., f þ E2 → f þ E1 and f − E2 → f − E1 , are considered. Note that the phase shift Φ has opposite direction for the two doublet frequencies due to a sign change of ζ; see Eqs. (C3) and (C4). The model predicts slightly different temperature dependences AðTÞ for the f þ E2 → f þ E1 and f − E2 → f − E1 transitions. Since the experimental frequency doublet is fit to a common amplitude, the model curve shown in Fig. 7(b) is the average of both contributions.
For a quantitative analysis, the following parameters (as obtained in the sections above) are used: in ZF, A E1 ¼ 0.082 and A E2 ¼ 0.090, in 4T, A E1 ¼ 0.0253 and A E2 ¼ 0.0275. The derivation of Eqs. (C2) and (C3) does not consider initial phases (ϕ E1 ¼ ϕ E2 ¼ 0). This is accounted for by shifting the model curves by the E1 initial phase obtained at low temperatures.

APPENDIX D: DETAILS ON LOCAL HOPPING SIMULATION
In order to investigate the E3 0 component and obtain a better understanding of the dynamic behavior over the full temperature range, a Monte Carlo simulation of the muon depolarization function in B ext ¼ 4Tkc assuming local hopping between adjacent sites is carried out. Stopping sites are arranged on a hexagon, with the in-plane component of B int pointing radially outward, and the signs of θ alternating. An exponential correlation time τ ¼ 1=ν, where ν is the average hop rate, is assumed. Polarization spectra are simulated in the range ν ¼ 10 −1 -10 5 MHz for 1 μs with a time step Δt ¼ 0.0001 μs and 5000 repeats, using f E3 ¼ 76.7 MHz for jB int j, θ ¼ AE4.69°, and f ext ¼ 541.98 MHz. Details on the general setup of the simulation can be found in Ref. [70]. The simulated spectra are fit to either a single oscillatory component, Eq. (1), for fast hopping, or to two oscillatory components with shared relaxation rate in the quasistatic regime. The resulting precession frequencies and relaxation rates are shown in Figs. 12(a) and 12(b).
For direct comparison, the simulation results are mapped onto the data assuming an Arrhenius-like activation T ¼ E b =ðk B ln½ν 0 =νÞ, where E b is an activation energy and ν 0 an attempt frequency. The experimental 4T E3=E3 0 =E3 Ã frequencies and relaxation rate are displayed in Figs. 12(c) and 12(d), alongside the simulation results for ν ¼ 5 × 10 5 MHz and E b ¼ 42 meV, showing excellent qualitative agreement. The E3 0 component is clearly identified as f av E3 , the average of the doublet frequencies; compare Fig. 6. Note that due to the approximation of a fixed internal field, and neglect of susceptibility contributions, the model loses validity with increasing temperature. Both at the low-and high-temperature end of the data, the relaxation rate is not well described, indicating that a simple Arrhenius activation model is insufficient to fully describe the data.
Overall, local hopping describes the E3=E3 0 =E3 Ã signals in the combined dataset taken in ZF, small and large B ext very well. There is convincing evidence that E3 Ã arises from muons undergoing thermally activated hopping between adjacent E3 sites, causing B int to average to zero. This is strongly supported by the observation of E3 0 in large B ext and its identification as f av E3 . Additionally, the energy barrier E b ¼ 42 AE 5 meV between sites is estimated.

APPENDIX E: DETAILS ON DFT CALCULATIONS
DFT calculations are carried out using the Vienna ab initio simulation package version 5.4.4 [39][40][41]. The positive muon is modeled as a hydrogen nucleus embedded within an 80-atom 2 × 2 × 2 rhombohedral supercell (SC) of Cr 2 O 3 . The local spin-density approximation as parametrized by Perdew and Zunger [71] with an additional Hubbard-like correction (LDA þ U) is used. The LDA þ U correction scheme of Dudarev et al. [72] is employed with a U eff of 4 eV applied to the Cr d states. This choice of U eff is found to provide a good description of the crystal and electronic structure of Cr 2 O 3 and is in line with values used in previous works [73]. Brillouin zone integrations are performed using the tetrahedron method with Blöchl corrections on a Γ-centered 4 × 4 × 4 Monkhorst-Pack grid [74] for the 2 × 2 × 2 SC. A plane-wave cutoff of 700 eV is used. With respect to a 900-eV cutoff and an 8 × 8 × 8 k-point mesh, energy differences are found to be converged to within 1 meV=formula unit and forces to within 3 meV=Å. The full convergence tests are available in Ref. [42]. Within the projector-augmented plane-wave (PAW) method [75,76], the 14 electrons for Cr (3s 2 3p 6 3d 4 4s 2 ) and six for O (2s 2 2p 4 ) are treated explicitly [77]. A collinear treatment of spins is adopted, and the well-established "þ − þ−" G-type antiferromagnetic order indicated in the inset of Fig. 1 is confirmed. Spin-orbit coupling has previously been found to be negligible in Cr 2 O 3 [73] and is therefore not included here. Supercells of up to 4 × 4 × 4 are tested in order to ensure that the imposed periodic boundary conditions do not introduce artifacts. In particular, we find negligible changes in the calculated contact and dipole contributions with respect to the results shown in Table II.