Quantum effects in muon spin spectroscopy within the stochastic self-consistent harmonic approximation

TL;DR

This paper introduces a quantum variational approach using the stochastic self-consistent harmonic approximation to accurately model the anharmonic vibrational effects of muons in materials, improving agreement with experimental data.

Contribution

It applies the stochastic self-consistent harmonic approximation to include anharmonic quantum effects in muon spin spectroscopy calculations, advancing beyond harmonic approximations.

Findings

Muon vibrational frequencies are strongly renormalized by anharmonicity.

Including quantum anharmonic fluctuations stabilizes the muon at the octahedral site in bcc Fe.

The method significantly improves agreement with experimental hyperfine field measurements.

Abstract

In muon spin rotation experiments the positive implanted muon vibrates with large zero point amplitude by virtue of its light mass. Quantum mechanical calculations of the host material usually treat the muon as a point impurity, ignoring this large vibrational amplitude. As a first order correction, the muon zero point motion is usually described within the harmonic approximation, despite the large anharmonicity of the crystal potential. Here we apply the stochastic self-consistent harmonic approximation, a quantum variational method devised to include strong anharmonic effects in total energy and vibrational frequency calculations, in order to overcome these limitations and provide an accurate ab initio description of the quantum nature of the muon. We applied this full quantum treatment to the calculation of the muon contact hyperfine field in textbook-case metallic systems, such as Fe, Ni, Co including MnSi and MnGe, significantly improving agreement with experiments. Our results show that muon vibrational frequencies are strongly renormalized by anharmonicity. Finally, in contrast to the harmonic approximation, we show that including quantum anharmonic fluctuations, the muon stabilizes at the octahedral site in bcc Fe.

Tables

Table 1. Table 1: Harmonic muon frequencies ω i h superscript subscript 𝜔 𝑖 ℎ \omega_{i}^{h} along the mode i 𝑖 i and harmonic zero-point energy E h = ∑ i = 1 3 ℏ ​ ω i h / 2 superscript 𝐸 ℎ superscript subscript 𝑖 1 3 Planck-constant-over-2-pi superscript subscript 𝜔 𝑖 ℎ 2 E^{h}=\sum_{i=1}^{3}\hbar\omega_{i}^{h}/2 , together with the SSCHA muon frequencies ω ~ i subscript ~ 𝜔 𝑖 \widetilde{\omega}_{i} and energy E ~ ~ 𝐸 \widetilde{E} at the minimum that includes the anharmonic contribution.

Host	${\omega }_x^h$ (cm-1)	${\omega }_y^h$ (cm-1)	${\omega }_z^h$ (cm-1)	$E^h$ (eV)	${\tilde{\omega }}_x$ (cm-1)	${\tilde{\omega }}_y$ (cm-1)	${\tilde{\omega }}_z$ (cm-1)	$\tilde{E}$ (eV)
Fe - bcc \tnotextnxx6	4364.01	2913.01	4364.62	0.72	4769.08	2572.58	5088.37	0.74
Fe - bcc \tnotextnxx62	1965.08$i$	1958.72$i$	6828.00	\tnotextnxx63	2005.24	2005.24	6364.81	0.53
Co - hcp	2930.41	2929.85	2752.25	0.53	3741.10	3741.10	3476.24	0.61
Co - fcc	2607.29	2607.02	2606.66	0.49	3424.16	3424.16	3424.16	0.56
Ni - fcc	2377.62	2377.60	2377.61	0.44	3317.78	3317.78	3317.78	0.53
MnGe	3123.70	3123.67	3123.66	0.58	3470.29	3470.29	3470.29	0.64
MnSi	3296.27	3296.32	3296.11	0.61	3685.25	3685.25	3685.25	0.67

Table 2. Table 2: Calculated contact hyperfine field for the point muon at the equilibrium position B c ​ ( 𝐫 μ e ​ q ) subscript 𝐵 𝑐 superscript subscript 𝐫 𝜇 𝑒 𝑞 B_{c}(\mathbf{r}_{\mu}^{eq}) , the calculated contact hyperfine field averaged over the spread of the muon wavefunction ⟨ B c ⟩ delimited-⟨⟩ subscript 𝐵 𝑐 \left<B_{c}\right> and experimentally observed values (Exp).

Host metals	$B_c\left({𝐫}_{\mu }^{eq}\right)$ [T] \tnotextnxx667	[T]	Exp
Fe-bcc \tnotextnxx666	-1.25	-1.07	-1.11 Nishida et al. (1977)
Fe-bcc \tnotextnxx668	-1.22	-1.13	-1.11 Nishida et al. (1977)
Co-hcp	-0.79	-0.64	-0.61 Graf et al. (1976)
Co-fcc	-0.73	-0.68	-0.58 Lindgren and Ellis (1982)
Ni-fcc	-0.15	-0.14	-0.071 Graf et al. (1979)
MnGe	-1.14	-1.07	-1.08 Martin et al. (2016)
MnSi	-0.22	-0.21	-0.207 Amato et al. (2014)

Equations

$$H_{tot}=T_e+T_{\mu }+T_N+V({\mathbf{r}}_e,{\mathbf{r}}_{\mu },{\mathbf{R}}_N),$$

$$H_{tot}|{\Psi }_{tot}⟩=E_{tot}|{\Psi }_{tot}⟩.$$

$$|{\Psi }_{tot}⟩=|{\psi }_e⟩|{\varphi }_{\mu }⟩|{\Phi }_N⟩.$$

$$H_e=T_e+V({\mathbf{r}}_e;{\mathbf{r}}_{\mu },{\mathbf{R}}_N).$$

$$H_{\mu }=T_{\mu }+V({\mathbf{r}}_{\mu };{\mathbf{R}}_N),$$

$$T_{\mu }=\sum _{\lambda =1}^3\frac{p_{\lambda }^2}{2m_{\mu }},$$

$$E=⟨\varphi \left|H\right|\varphi ⟩,$$

$$\tilde{E}=⟨\tilde{\varphi }\left|\tilde{H}\right|\tilde{\varphi }⟩.$$

$${\tilde{E}}_H[\tilde{H}]=⟨\tilde{\varphi }\left|H\right|\tilde{\varphi }⟩.$$

$$E\le {\tilde{E}}_H[\tilde{H}]$$

$${\tilde{E}}_H[\tilde{H}]=\tilde{E}+⟨\tilde{\varphi }\left|(V-\tilde{V})\right|\tilde{\varphi }⟩.$$

$$\tilde{\rho }({\mathbf{r}}_{\mu })=⟨{\mathbf{r}}_{\mu }|\tilde{\varphi }⟩⟨\tilde{\varphi }|{\mathbf{r}}_{\mu }⟩={\left|\tilde{\varphi }({\mathbf{r}}_{\mu })\right|}^2.$$

$${⟨A⟩}_{\tilde{H}}=\int d{\mathbf{r}}_{\mu }A({\mathbf{r}}_{\mu })\tilde{\rho }({\mathbf{r}}_{\mu })$$

$${\tilde{E}}_H[\tilde{H}]=\tilde{E}+\int d{\mathbf{r}}_{\mu }\tilde{\rho }({\mathbf{r}}_{\mu })(V({\mathbf{r}}_{\mu })-\tilde{V}({\mathbf{r}}_{\mu })).$$

$$\sum _{i=1}^3\frac{1}{2}\hslash {\tilde{\omega }}_i=\tilde{E}.$$

$$B_c({\mathbf{r}}_{\mu })=\frac{2}{3}{\mu }_0{\mu }_B\left[n_{\uparrow }({\mathbf{r}}_{\mu })-n_{\downarrow }({\mathbf{r}}_{\mu })\right],$$

$$⟨B_c⟩=\int d{\mathbf{r}}_{\mu }{B}_c({\mathbf{r}}_{\mu }){\left|\varphi ({\mathbf{r}}_{\mu })\right|}^2.$$

$$\int d{\mathbf{r}}_{\mu }{B}_c({\mathbf{r}}_{\mu })\tilde{\rho }({\mathbf{r}}_{\mu })\simeq \frac{1}{N_c}\sum _{n=1}^{N_c}{B}_c({\mathbf{r}}_{\mu }^n)\equiv {⟨B_c⟩}_{\tilde{H}}.$$

$$\tilde{H}=\sum _{\lambda =1}^3\frac{p_{\lambda }^2}{2m_{\mu }}+\frac{1}{2}\sum _{\lambda \nu }^3{K}^{\lambda \nu }(r_{\mu }-r_{\mu }^{eq}{)}^{\lambda }(r_{\mu }-r_{\mu }^{eq}{)}^{\nu },$$

$$\sum _{\nu =1}^3\frac{K^{\lambda \nu }}{m_{\mu }}{\widetilde{ϵ}}_i^{\nu }={\tilde{\omega }}_i^2{\widetilde{ϵ}}_i^{\lambda },$$

$$\tilde{\rho }({\mathbf{r}}_{\mu })=\frac{1}{\prod _{i=1}^3\sqrt{2\pi {\tilde{\sigma }}_i^2}}\exp \left(-\sum _{\lambda \nu i}^3\frac{{\widetilde{ϵ}}_i^{\lambda }{\widetilde{ϵ}}_i^{\nu }}{2{\tilde{\sigma }}_i^2}(r_{\mu }-r_{\mu }^{eq}{)}^{\lambda }(r_{\mu }-r_{\mu }^{eq}{)}^{\nu }\right),$$

$${\tilde{\sigma }}_i=\sqrt{\frac{\hslash }{2m_{\mu }{\tilde{\omega }}_i}}.$$

$$\tilde{E}=\sum _{i=1}^3\frac{1}{2}\hslash {\tilde{\omega }}_i.$$

$$(r_{\mu }^n{)}^{\lambda }=(r_{\mu }^{eq}{)}^{\lambda }+\sum _{i=1}^3{\widetilde{ϵ}}_i^{\lambda }{\tilde{\sigma }}_i{\xi }_{in}.$$

$$\begin{array}{rl}{\nabla }_K{\tilde{E}}_H[\tilde{H}]=& -\sum _{i\lambda \nu }\left({\widetilde{ϵ}}_i^{\lambda }{\nabla }_K\ln {\tilde{\sigma }}_i+{\nabla }_K{\widetilde{ϵ}}_i^{\lambda }\right){\widetilde{ϵ}}_i^{\nu }\times \\ & \int d{\mathbf{r}}_{\mu }[f^{\lambda }({\mathbf{r}}_{\mu })-{\tilde{f}}^{\lambda }({\mathbf{r}}_{\mu })](r_{\mu }-r_{\mu }^{eq}{)}^{\nu }\tilde{\rho }({\mathbf{r}}_{\mu }),\end{array}$$

Full text

I Introduction

In muon spin rotation ($\mu$SR) experiments, spin-polarized positive (anti)muons are used to probe the microscopic field distribution at the interstitial site(s) where the $\mu^{+}$ stop inside the sample under investigation. The extreme sensitivity of the muon to small magnetic fields as well as the absence of quadrupolar coupling makes this technique very effective in probing magnetic orders, offering a valuable alternative to neutron scattering. This approach, which shares many similarities with nuclear magnetic resonance, has the advantage of being applicable to virtually any material, but it has the drawback that the interstitial sites where the muon stops and the nature of muon interaction with the host are generally unknown. Here we discuss an improved method to tackle this problem based on computational chemistry methods.

An accurate, ab initio, description of the electron-muon interaction in periodic solids has been out of reach until a few years ago. The dramatic increase of both the computational power and the accuracy of first-principles calculations make this goal possible. Self-consistent electronic structure calculations, in particular those based on density functional theory (DFT), are already employed to study the muon implantation site, muon interaction parameters, and for understanding the muon-induced distortion in the lattice J.Rath et al. (1979); Möller et al. (2013, 2013); Bernardini et al. (2013); Bonfà and De Renzi (2016); Onuorah et al. (2018); Liborio et al. (2018). This turns out to be a very valuable tool for analyzing experimental data and interpreting the results Disseler (2014). The knowledge of the muon implantation site(s) and of the hyperfine field allows very important quantitative information, including the magnetic structure and the moment size, to be obtained from $\mu$SR experiments. Moreover, a reliable quantum calculation of the muon embedded in the system under investigation provides an estimate for its induced perturbation; the probe is an impurity and it may in principle alter the local electronic properties. Fortunately this is a very rare case, and yet assessing these rare cases Foronda et al. (2015); Dalmas de Réotier et al. (2018) is very important.

However, self-consistent DFT calculations often treat the muon as just another atom in the lattice, within the Born-Oppenheimer (BO) approximation Born and Oppenheimer (1927), without taking into consideration the quantum effect of the muon zero-point vibrations, which is sizable relative to those of heavier nuclei. The embedded muon, by virtue of its very light mass ($\sim 1/{9^{th}}$ the proton mass), is characterized by zero-point vibration with amplitude typically of the order of 1 Bohr radius J.Rath et al. (1979). The neglect of this effect may have two major consequences: inaccurate estimation of the contact hyperfine field, and/or incorrect identification of muon implantation sites. The former is due to neglect of the space extent of the muon wavefunction, whereas the latter happens when the quantum zero-point vibration energy is comparable with the energy difference between the various implantation sites Möller et al. (2013, 2013); Bonfà et al. (2015); Bonfà and De Renzi (2016).

Earlier approaches towards a quantum-mechanical description of the muon zero-point vibration include calculations within the harmonic approximation Boxwell et al. (1993); Möller et al. (2013). However, the muon potential has been discussed and shown to be anharmonic, for instance by total energy calculations with site exploration algorithms J.Rath et al. (1979); Porter et al. (1999); Herrero et al. (2006); Bonfà et al. (2015). Furthermore, a break down of the harmonic approximation takes place when within the range of the muon vibrations the potential is not dominated by the second-order term in its Taylor expansion.

Alternative methods do take into account the anharmonic nature of the crystalline potential. One of them consists in the potential exploration approach Bonfà et al. (2015). The non-BO methods represent another computationally demanding alternative, employing a linear combination of Gaussian basis functions to realize both the nuclear and the electronic degrees of freedom Tachikawa et al. (1998); Nakai (2001); Webb et al. (2002); Kerridge et al. (2004) and optimized local potentials to represent the nuclear-electron correlation N.I. and K. (2014). One of the most advanced approaches relies on ab initio path integral molecular dynamics, which allows for contextual quantization of both the muon and the electrons in the calculation of the electronic structure and of the interatomic forces  Valladares et al. (1995); Probert and Glover (2006); Herrero et al. (2006). However, computational resources required by this method grow exceedingly with the size of the cell.

In this paper, we describe a stochastic self-consistent harmonic approximation (SSCHA) that allows us to include the effects of anharmonicity in the muon vibrations Errea et al. (2013, 2014); Bianco et al. (2017); Monacelli et al. (2018). The SSCHA is a quantum variational method that efficiently calculates anharmonic free energies and phonon frequencies in a non-perturbative way. This approach has been very successful for calculating phonon frequencies and superconducting properties in hydrogen-rich materials, as well as in systems undergoing charge density wave (CDW) transitions, ferroelectrics, and thermoelectrics Errea et al. (2013, 2015); Leroux et al. (2015); Errea et al. (2016); Ribeiro et al. (2018); Aseginolaza et al. (2019); Bianco et al. (2019). For the muon, the SSCHA is variational in the muon (free) energy, with this energy evaluated stochastically from forces and energies calculated at a sufficient number of random muon configurations. The muon energy is minimized using trial harmonic wavefunctions that are Gaussian, while the minimization parameter is the width of the Gaussian. From the output of the minimization, muon frequencies including anharmonic contributions and the muon ground-state energy can be extracted.

With this approach, we demonstrate that there are anharmonic contributions to the harmonic muon vibrational modes, as expected for the muon due to its light mass. We further use the SSCHA muon wavefunction to refine the contact hyperfine field in a series of metals: Fe, Ni, Co, MnSi and MnGe, where the SSCHA improves the agreement of the calculated value with the experimental results, with respect to recent point impurity calculations Onuorah et al. (2018). Finally, the SSCHA together with energy curvature considerations Bianco et al. (2017) allows the stable occupation of the muon at the octahedral site in bcc Fe, which is unstable within the harmonic regime.

The paper has the following structure: Sec. II discusses the double Born-Oppenheimer approximation, which allows us to separate the muon degrees of freedom from those of the host nuclei and electrons. In Sec. III, we describe the working principles of the SSCHA, including the stochastic implementation. In Sec. IV, we discuss the muon zero-point energy calculation results using the SSCHA together with the stability of the muon at octahedral and tetrahedral sites in Fe(bcc). Finally, in Sec. VI we present the results of the quantum corrections in the calculation of the contact hyperfine field and then conclusions are given in Sec.VII.

II Double Born-Oppenheimer approximation

The BO approximation considers the nuclei frozen on the time scale of electron dynamics in view of their sufficiently large mass ratio Born and Oppenheimer (1927). Hydrogen is already sufficiently lighter than most other atoms to allow a further separation of time scales, and this holds a fortiori true for a positive muon. This allows for the quantum treatment of a single muon impurity in the crystal by employing the so-called double Born-Oppenheimer approximation (DBO) Soudackov and Hammes-Schiffer (1999); Porter et al. (1999); Bonfà et al. (2015). The muon dynamics ($m_{\mu}\sim 200m_{e}$) is much slower than that of electrons, thus justifying an electron structure obtained by DFT with frozen muon and nuclei. The same muon dynamics is still much faster than that of other nuclei, since transition metals are typically 400 times heavier than a muon (care must be taken when considering e.g., hydrogen, which is only nine times heavier than a muon). Therefore it is justified to use total DFT energy versus the muon configuration coordinates as a frozen potential energy landscape in which the muon dynamics takes place on its characteristic time scale. This allows us to consider the zero-point vibration of only the muon within the potential energy surface, drastically reducing the computational load requirements for the calculations.

The total Hamiltonian $H_{tot}$ describing the many-body interaction including explicitly the muon coordinates is written as

[TABLE]

with subscript $\mu$ describing the muon-related quantities while $e$ and $N$ describe those of the electrons and host nuclei respectively. $T$ and $V$ are the kinetic and potential energy, respectively. The Schrödinger equation is then written as

[TABLE]

This further allows us to write the DBO wavefunction as a product wavefunction of the electrons, the muon and the nuclei in the form:

[TABLE]

The Hamiltonian for the electronic problem can be re-written to specifically point out the presence of the muon position operator as

[TABLE]

Similar to the BO approximation, only the position operators of the muon and the nuclei enter in the eigenvalue problem of the electrons. The solution of the electronic problem gives the BO potential energy surface, $V(\mathbf{r}_{\mu},\mathbf{R}_{N})$, dependent on the muon and the nuclei position operators.

Hence, the ground-state Hamiltonian $H_{\mu}$ for the muon can be written as

[TABLE]

where the muon kinetic energy $T_{\mu}$ is defined as

[TABLE]

with $p$ the momentum operator along the Cartesian component indexes $\lambda$, while $m_{\mu}$ is the muon mass.

The acquisition of the DBO potential energy surface $V(\mathbf{r}_{\mu};\mathbf{R}_{N})$ for the solution of the Schrödinger equation Eq. 5, is still a long and difficult task. However, the DBO approximation is advantageous since it allows to consider separately only the degrees of freedom of the muon. For this reason, in the next section we revisit the SSCHA theory originally presented in Ref. Errea et al. (2013, 2014) specializing its application to the muon dynamics.

III Stochastic self-consistent harmonic approximation (SSCHA) for muons

To begin with the formal description of the SSCHA restricted only to the muon modes, let us write the muon Hamiltonian $H_{\mu}$, the muon wavefunction $\phi_{\mu}$, and the DBO potential energy surface $V(\mathbf{r}_{\mu};\mathbf{R}_{N})$, appearing in the previous section, simply as $H$, $\phi$ and $V(\mathbf{r}_{\mu})$ respectively.

The muon zero-point energy from the Hamiltonian $H$ is given as

[TABLE]

where $\left|\phi\right>$ is the muon ground-state wavefunction. Calculating $E$ is far from trivial since the form of the muon potential (Eq. 5) is not known. However, it is possible to establish a quantum variational principle for the muon ground state energy $E$, by replacing the exact muon wavefunction $\left|\phi\right>$ with the wavefunction $\left|{\widetilde{\phi}}\right>$ of a trial muon Hamiltonian $\widetilde{H}=T_{\mu}+\widetilde{V}(\mathbf{r}_{\mu})$ with energy

[TABLE]

This is such that one can define an energy functional of the trial Hamiltonian as

[TABLE]

The variational form of the muon ground state energy can be written as

[TABLE]

such that the equality holds when the true and trial potentials are the same.

By adding and subtracting Eq. 7 to Eq. 8, $\widetilde{E}_{H}[\widetilde{H}]$ is written in the form

[TABLE]

The above definitions allow to formulate a variational principle following the Gibbs-Bogoliubov inequality theorem Isihara (1968) at zero temperature, similar to the Rayleigh-Ritz inequality MacDonald (1933).

According to the trial wavefunction, the probability to find the muon in the position $\mathbf{r}_{\mu}$ is

[TABLE]

Thus, an observable A dependent only on $\mathbf{r}_{\mu}$ can be averaged statistically within the form of the corresponding Hamiltonian $\widetilde{H}$ as

[TABLE]

and the muon energy in Eq. 10 can be evaluated as

[TABLE]

With the above form of $\widetilde{E}_{H}[\widetilde{H}]$, the muon energy can be evaluated at each step during the variational minimization. One can directly see that the equality in the form of the variation in Eq. 9 holds if $V=\widetilde{V}$. Hence, with the variational principle, the ground state of the muon is determined if the potential $\widetilde{V}(\mathbf{r}_{\mu})$ that minimizes $\widetilde{E}_{H}[\widetilde{H}]$ is found.

To proceed with the minimization of $\widetilde{E}_{H}$, in the SSCHA implementation we restrict the muon wavefunctions only to the Gaussian form. The term harmonic in the technique refers to the fact that each Gaussian is the ground state of a trial harmonic Hamiltonian, with known analytic solutions (see Appendix A) where the trial potential is expressed in terms of a force constant matrix. Moreover, using Gaussian functions has the advantage of allowing to sample the wavefunction by extracting randomly distributed configurations without any Metropolis algorithm that requires long equilibration time and also provides an analytic expression for the kinetic energy.

Finally, the actual minimization is obtained using the conjugate gradient (CG) algorithm Hestenes and Stiefel (1952) which requires the evaluation of the energy gradient, whose analytic form is given in Ref. Errea et al. (2014) and in Appendix B for the muon case, and depends on the forces acting on the muon when displaced from the equilibrium position.

The evaluation of the quantities of interest at each minimization step, namely $\widetilde{E}_{H}$ and its gradient, is performed stochastically. One of the advantages of the stochastic sampling resides in the gradual optimization of the potential felt by the muon during the iterative process. This ensures that the entire BO landscape, beyond the harmonic component around the minimum, is sampled, hence capturing the anharmonic effects.

The stochastic sampling of the BO energy and of the forces acting on the muon and entering the energy gradient (see Appendix B) can be calculated with any ab initio method including DFT Kohn and Sham (1965) and Hartree-Fock Fock (1930); S (1947); Slater (1951) based approaches.

The evaluation of the forces and energies for the random muon configurations in the stochastic sampling represents the most computationally demanding task in the SSCHA minimization cycle. This effort can be partially alleviated with a re-weighting procedure based on importance sampling. The reader is referred to Ref Errea et al., 2014 for a detailed description of this additional detail.

When the energy gradient numerically vanishes, the $\widetilde{E}$ that minimizes $\widetilde{E}_{H}[\widetilde{H}]$ is the zero-point energy of the muon and the anharmonic vibrational frequencies $\widetilde{\omega}_{i}$ of the auxiliary Hamiltonian whose SSCHA wavefunction is the ground state are obtained, so that

[TABLE]

The formal description of the trial Hamiltonian and the trial wavefunction is given in Appendix A.

IV Muon zero-point energy

Let us first describe the zero-point energy of the muon obtained in the harmonic approximation, which is later used in comparisons with the anharmonic one.

The harmonic muon frequencies $\omega^{h}_{i}$ and the corresponding energies $E^{h}=\sum_{i=1}^{3}\hbar\omega_{i}^{h}/2$ were calculated by the finite difference method Kresse et al. (1995); Parlinski et al. (1997), which allows only the muon frequencies to be singled out, for all the materials under investigation, namely Fe, Co, Ni, MnSi and MnGe. These were also used to generate the starting wavefunctions for the SSCHA minimization except for the stability discussion in Sec. V with the muon at the octahedral and tetrahedral site in bcc Fe. Here, the density functional perturbation theory (DFPT) within the Quantum ESPRESSO suite of code Baroni et al. (2001); Giannozzi et al. (2009) was used to calculate the frequencies of the whole system, including those of the host Fe nuclei. The resulting harmonic muon frequencies from both methods in the two Fe systems are in good agreement.

For the SSCHA minimization and stochastic averaging (see Eq 17), hundreds (100 to 400) of random configurations were generated for the muon, while keeping the host atoms fixed, to ensure that the muon energy gradient vanishes. Their energies and Hellmann-Feynman forces Pulay (1969) were calculated by DFT as implemented in the Quantum ESPRESSO suite of code Giannozzi et al. (2009). The details of the muon site in these systems and DFT input parameters are contained in Ref. Onuorah et al. (2018). For all the systems, a 2$\times$2$\times$2 supercell constructed starting from the conventional unit cell was used for the harmonic frequency calculations, the SSCHA frequency minimization and the force calculation within DFT. Other DFT computational details are identical to those reported in Ref. Onuorah et al. (2018). To accommodate the muon impurity in the supercell, the forces introduced by the muon in the system were relaxed by DFT and the relaxed structures were used for the SSCHA calculations. Relaxations were converged with force and energy thresholds of 10*-3* a.u and 10*-4* Ry respectively.

Figure 1 shows the evolution of the SSCHA muon frequencies and energy during the minimization procedure. Significant anharmonic contributions to the resulting SSCHA frequencies can be deduced from the difference between the initial values, i.e. the starting harmonic guess, and the final converged results (the comparison with the anharmonic correction obtained for host atoms is presented in Appendix  C). The anharmonic correction to the harmonic frequencies is found to be in the range of 330 - 820 cm*-1* except for the muon at the octahedral site in Fe.

The stochastic implementation ensures that the effect of the muon vibrations, the effect of the chemical environment around the muon and anharmonic contributions to the forces acting on the muon are all incorporated in the muon ground state minimum.

Table 1 contains the harmonic frequencies $\omega^{h}_{i}$ and energies $E^{h}$, obtained with the finite difference method and used as the starting point of the SSCHA iterative process, and the SSCHA frequencies $\widetilde{\omega}_{i}$ and energies $\widetilde{E}$ at the end of the minimization. The error estimates of the reported muon energies are within the range of 0.1 meV. The results show the anharmonic effects in the muon vibrational frequencies. Notice that the muon at octahedral implantation site in Fe is unstable in the harmonic regime. For all other cases with positive harmonic frequencies for which $E^{h}$ can be defined, the difference between the SSCHA muon vibrational energies and the harmonic ones is in the range of 0.02 - 0.09 eV.

V Tetrahedral and octahedral muon site in Fe

Conflicting experimental and theoretical studies report the muon site in Fe to be either at the tetrahedral (T) or the octahedral (O) interstitial sites Seeger (1976); Graf et al. (1980); Lindgren and Ellis (1982); Estreicher and Meier (1982); Yagi et al. (1984); Kossler et al. (1985). From the point of view of the DFT total energy, the T site is 0.184 eV lower that the O site. This would indicate that the T site is the stable one. However, since the calculated muon zero-point energies (above 0.5 eV) are large relative to the DFT energy difference, the possible population of both sites cannot be excluded.

DFPT calculations of the muon frequencies provide further insight into the stability of the two candidate sites. Unphysical negative frequencies, generally a signal of instability, are obtained for the muon at the O site, as opposed to those of the T site, which are always positive. The harmonic approximation then appears to indicate an instability of the muon at the O site.

However, the anharmonic effects, fully captured by the SSCHA, yield positive frequencies also for the O site indicating that the instability is an artifact of the harmonic approximation. As the $\widetilde{\omega}_{i}$ frequencies are positive-definite by definition, this is not proof that the O site occupation is stable. Obtaining the frequencies from the energy curvature Bianco et al. (2017), which can correctly describe an instability, confirms, however, that the O site interstitial site is in fact stable. The SSCHA frequencies for the muon in the O site are larger than the frequencies resulting from those obtained from the curvature by only 0.53 % along the $x$,$y$ axis and 0.14 % along the $z$ axis.

The quantum correction with the SSCHA shows that both T and O are stable local minima. The vibrational contribution to the energy is 0.21 eV less for the O site than for the T site (see Table 1). Adding this to the static DFT contribution makes the O site energetically favored by approximately 0.03 eV over the T site, thus indicating that the two sites are basically degenerate, and possibly both occupied.

VI Quantum corrections on a muon contact hyperfine field

The contact hyperfine field $B_{c}(\mathbf{r}_{\mu})$ at the muon position $\mathbf{r}_{\mu}$ is computed ab inito by considering the imbalance in the spin density at the muon site Onuorah et al. (2018) given as

[TABLE]

where $\mu_{0}$ is the vacuum permeability, $\mu_{B}$ is the Bohr magneton and $n_{\uparrow}-n_{\downarrow}$ represents the spin polarization at the muon site $\mathbf{r}_{\mu}$ calculated here by DFT. $B_{c}(\mathbf{r}^{eq}_{\mu})$ has been calculated in this way for metals within a point impurity treatment of the muon Onuorah et al. (2018). We now calculate the effect of the muon quantum delocalization on its contact hyperfine field, using the muon SSCHA wavefunctions $\phi$ that already contain the anharmonic contributions.

The quantum expectation value, $\left<B_{c}\right>$ is given by

[TABLE]

where the probability density $\left|\phi(\mathbf{r}_{\mu})\right|^{2}$ has been defined in Eqs. 11 and is obtained from the SSCHA muon frequencies $\widetilde{\omega}_{i}$ according to Eq.  20.

The above integral can be evaluated in a post-DFT calculation by a statistical average performed stochastically, i.e. according to

[TABLE]

where the sum extends over a number of muon random configurations $N_{c}$ displaced from the equilibrium position $\mathbf{r}^{eq}_{\mu}$ and generated with the probability distribution of the muon wavefunction (see Eq. 23). The number of muon random configurations used is the same as in the SSCHA minimization of the muon wavefunction. However, the new muon random positions are generated considering the anharmonic corrected SSCHA muon wavefunction. Figure 2 shows the distribution of the 100 configurations used for fcc Co in the unit cell.

$B_{c}(\mathbf{r}_{\mu})$ was calculated by DFT for each of these random configurations within a 3$\times$3$\times$3 supercell for Fe, Co and Ni and 2$\times$2$\times$2 supercell for MnGe and MnSi, while other computational details are the same as reported in Ref. Onuorah et al. (2018).

Table 2 and Fig. 3 show the calculated contact field $B_{c}(\mathbf{r}^{eq}_{\mu})$ for a point-like muon Onuorah et al. (2018) and its stochastically averaged $\left<B_{c}\right>$ values together with the experimental values. For all the systems the statistical error for the stochastic sampling of $\left<B_{c}\right>$ is in the range of $\approx$ 1 mT. The contact hyperfine field including quantum correction within the SSCHA, $\left<B_{c}\right>$, improves the agreement with the experiments, thus underlining the importance of considering the finite muon wavefunction when computing muon hyperfine interactions. Admittedly, the correction to the contact hyperfine field appears to be less relevant than the outcome obtained on the stability of the muon at the octahedral site in Fe, still $|\left<B_{c}\right>|$ introduces a correction that ranges between 1 and 18%.

VII Conclusion

In conclusion, we have presented a general, effective and robust approach, based on the DBO approximation, to obtain the ground state wavefunction and zero-point energy of a positive muon embedded in a crystal from first principles. The adaptation of the SSCHA to the muon case allows us to evaluate the delocalized muon wavefunction including anharmonic contributions, that correct harmonic ones.

Moreover, the SSCHA circumvents the problem of directly reconstructing the potential energy surface by replacing this task with a variational problem, and more importantly, it provides a computationally tractable method to describe the zero-point energy of the muon. This leads to a number of important insights concerning the stability of the muon sites and its coupling with the surrounding electrons.

The first point has been discussed by considering the case of the muon site in Fe, where anharmonicity plays a crucial role in establishing the stability of the muon in the tetrahedral and octahedral sites.

We reformulated the calculation of the muon contact hyperfine field by including the effects of its anharmonic zero-point vibration, improving the agreement with experiments with respect to previous estimates based on the point impurity treatment of the muon. Even though the correction is small, in numerous cases the contact field is of the order of tenths of a Tesla, thus making the absolute value of the correction presented here quite relevant.

Finally, the clean iterative procedure of the SSCHA makes it rather straightforward to define standardized workflows to automate the computational procedure. This represents another step towards routinely supporting experimental data analysis with computational simulation results.

Acknowledgements.

RDR acknowledges grants from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 654000. RDR, PB and IJO also acknowledge computing resources provided by, the Swiss National Supercomputing Centre (CSCS) under Project ID sm16, CINECA under Project ID IsC58, the STFC Scientific Computing Department’s SCARF cluster and the HPC resources at the University of Parma, Italy. IE acknowledges funding from the Spanish Ministry of Economy and Competitiveness (FIS2016-76617-P). This work is part of the PhD thesis of IJO at the University of Parma, Italy.

Appendix A The trial muon harmonic Hamiltonian

The trial muon harmonic Hamiltonian is of the form

[TABLE]

where $\lambda$ and $\nu$ are Cartesian component indexes, $\mathbf{r}_{\mu}^{eq}$ is the muon equilibrium position, $m_{\mu}$ is the mass of the muon and $K^{\lambda\nu}$ is the muon $3\times 3$ force constant matrix. The force constant matrix $K^{\lambda\nu}/m_{\mu}$ can be constructed and diagonalized as

[TABLE]

where $i$ is the index of each of the orthogonal modes, $\widetilde{\epsilon}_{i}^{\nu}$ is the polarization vector and $\widetilde{\omega}_{i}$ is the muon frequency corresponding to the trial Hamiltonian $\widetilde{H}$ for each mode.

Assuming a trial harmonic potential, the probability to find the muon at $\mathbf{r}_{\mu}$ can be written simply as

[TABLE]

where $\widetilde{\sigma}_{i}$, the normal length for each of the modes $i$, is given as

[TABLE]

Using the quantum statistical averaging defined in Eq. 12, the energy of the trial harmonic Hamiltonian can be calculated as

[TABLE]

Appendix B Random configuration sampling and minimization details

The distribution for the generation of the random muon position configurations is realized using random numbers $\{\xi_{in}\}_{n=1,...,{N_{c}}}$ generated with the Gaussian distribution ${\widetilde{\rho}}(\mathbf{r}_{\mu})$ and re-scaled by the corresponding normal length modes $\widetilde{\sigma}_{i}$ and polarization vector $\widetilde{\epsilon}_{i}^{\lambda}$. The generated positions are thus obtained as

[TABLE]

This constitutes the set of points used in the stochastic evaluation of $\widetilde{E}$ and of the gradient of the energy functional, namely $\nabla_{K}\widetilde{E}_{H}[\widetilde{H}]$, with respect to the force constant $K$. The analytic form of this last term is written as (see also Ref. Errea et al. (2014))

[TABLE]

where $\mathit{f}^{\lambda}(\mathbf{r}_{\mu})$ is the muon force component in the $\lambda$ Cartesian direction for all muon positions $\mathbf{r}_{\mu}$ and ${\widetilde{\mathit{f}}}^{\lambda}(\mathbf{r}_{\mu})$ are the forces obtained with the ${\widetilde{V}}$ potential. The SSCHA minimization is performed respecting the symmetries of the crystal Errea et al. (2014).

We also add that with the SSCHA, it is possible to minimize the energy both with respect to the muon position $\mathbf{r}_{\mu}$ and also the force constant matrix $K^{\lambda\nu}$. However, for the materials considered in this paper, there is sufficient knowledge of the equilibrium muon position $\mathbf{r}_{\mu}^{eq}$. Hence, the muon energy is only minimized with respect to the force-constant matrix $K$. For the muon in a high symmetry position, the force-constant matrix is a 3$\times$3 matrix, with the diagonal elements of the matrix accounting for the dominant contribution.

Finally, it is important to note that, in order to obtain physical phonons from the ground-state minimized quantities provided by the SSCHA, the second derivative (curvature) of SSCHA energy at the minimum with respect to $\mathbf{r}_{\mu}$ has to be calculated Bianco et al. (2017), which includes a correction term to the force constants matrix $K^{\lambda\nu}/m_{\mu}$ in Eq. 19. We verified that for the cases under study here the muon frequencies are affected by less than a 1% by this extra correction. Thus, we can treat the $\widetilde{\omega}_{i}$ frequencies as the physical vibrational energies of the muons.

Appendix C Evolution of muon and host Fe SSCHA frequencies

The evolution of the frequencies in the SSCHA calculation including anharmonic effects both for the Fe host nuclei and the muon at the tetrahedral site in a 2$\times$2$\times$2 supercell is shown in Fig. 4. The figure indicates that there is a significant anharmonic contribution to the muon eigenfrequencies after several iterations (upper panel), whereas the lower frequency modes of the heavier Fe nuclei (lower panel) remain negligibly changed. This consideration together with the DBO approximation discussed in sec. II also supports separating and concentrating only on the muon degrees of freedom.