Local electronic structure of dilute hydrogen in $\beta$-Ga₂O₃ probed by muons

Gallium trioxide ($\beta$-Ga₂O₃) is attracting attention as a material for high-voltage power devices and other applications because of its large band gap ($E_{\mathrm{g}}\sim$4.9 eV) and associated high critical electric field Galazka (2018), where control of electrical activity by $p$/$n$-type carrier doping is one of the critical issues. While Sn or Si doping is known to induce $n$-type conductivity Ueda et al. (1997); Orita et al. (2000), it has been pointed out that impurity hydrogen (H) is the possible cause of unintentional carrier doping in as-grown samples Galazka et al. (2020).

From infrared spectroscopic studies of $\beta$-Ga₂O₃ containing macroscopic amounts of H, a variety of H states including O-H defects have been observed Weiser et al. (2018); Qin et al. (2019). However, their relationship to carrier doping is not always clear due to the lack of information on their local electronic structures and H valence. Furthermore, since isolated H exists only in trace amounts in solids, experimental means to study its microscopic electronic state are limited. In this regard, theoretical studies based on density functional theory (DFT) calculations have played an important role, and for H in $\beta$-Ga₂O₃  DFT calculations have shown that interstitial H may be indeed the origin of $n$-type doping Varley et al. (2010, 2011); Li and Robertson (2014). This prompted us to introduce muons as pseudo-H into $\beta$-Ga₂O₃ for investigating the corresponding interstitial H states by analyzing their electronic and dynamical properties in detail by muon spin rotation ($\mu$SR).

The positive muon ($\mu^{+}$) is an elementary particle with a mass of $0.1126m_{\mathrm{p}}$ (where $m_{\mathrm{p}}$ is the proton mass). The local electronic structure of muon in matter is determined by the muon-electron interaction which is practically equivalent to H, as can be seen from the fact that the difference in reduced electron mass between a muon-electron bound state (muonium) and a neutral H atom is only 0.43%. Therefore, Mu in matter can be regarded as a light isotope of H ($=^{0.1126}$H; we introduce the elemental symbol Mu below to denote muon as pseudo-H).

On the other hand, when interpreting the results of the $\mu$SR experiment based on the results of the DFT calculations, there is an important caveat that should be pointed out at this stage. Generally, muons are implanted as an ion beam of relatively high-energy (typically $\sim$4 MeV), and the associated electronic excitations in the insulator crystals produce free carriers and excitons ($\sim$10³ per muon) Thompson (1974); Alig and Bloom (1975); Itoh (1997). These often propagate rapidly in the crystal, and there is experimental evidence that Mu acts as a capture center for them Prokscha et al. (2007). Consequently, the electronic state of Mu is a relaxed-excited state generated by the interaction with free carriers and excitons, which can be different from the electronic state under thermal equilibrium as expected from the thermodynamic double charge conversion level ($E^{+/-}$) obtained by DFT calculations.

Another well-known evidence that Mu is in a relaxed-excited state is that two or more different electronic states of Mu are often observed simultaneously in the same material Lichti et al. (2008). In order to understand the origin of this phenomenon, we have recently conducted an extensive investigation to explore the regularities between the Mu valence states observed in various oxides and the results of previous DFT calculations. As a result, assuming that these relaxed-excited states are associated with acceptor ($E^{0/-}$) and donor ($E^{+/0}$) levels that are predicted by DFT calculations as metastable states, we found that the observed Mu states can be coherently explained by the relationship between acceptor/donor levels and band structure Hiraishi et al. (2022). We call it the “ambipolarity model” since such behavior of Mu is a manifestation of the ambipolarity of H through the relaxed-excited states.

In this study, we show that Mu in single-crystalline $\beta$-Ga₂O₃ exhibits two different electronic structures corresponding to the relaxed-excited states respectively associated with the donor and acceptor levels in the ambipolarity model. One is in the OMu-bonded state (Mu${}_{1}^{+}$) corresponding to H serving as a donor, and another is the hydride-like state in rapid motion (Mu${}_{2}^{-}$). In particular, the Mu₂ state is a component overlooked in the previous $\mu$SR experiments King et al. (2010); Celebi et al. (2012), and it is speculated to be a transient state diffusing rapidly along the $\langle 010\rangle$ axis while undergoing the charge exchange reaction, Mu${}_{2}^{-}\rightleftarrows{\rm Mu}_{2}^{0}+e^{-}$. The occurrence of charge exchange is supported by the observation that the temperature dependence of the fractional yield of Mu₂ exhibits a strong correlation with that of the bulk carrier electron mobility and density. This also implies that the interstitial H can take hydride state under electronic excitation and may exhibit fast diffusion motion, depending on the bulk electronic properties of the host.

The sample used in the $\mu$SR experiment was a slab of single crystal (sc-Ga₂O₃, 10$\times$15$\times$0.6 mm³) with $\langle 001\rangle$ plane synthesized by the edge-defined film-fed (EFG) method (provided by Novel Crystal Technology, Inc.) Kuramata et al. (2016). It is reported to have no twin boundary and the lowest carrier density ($N_{\mathrm{e}}\sim 2\times$10¹⁷cm^-3) commercially available Kuramata et al. (2016). The electrical conductivity and Hall effect measurements were performed using the PPMS (Quantum Design Co.). Ohmic contacts were formed from vacuum-deposited Ti at room temperature Víllora et al. (2008); Irmscher et al. (2011) and gold paste (Seishin Trading Co. LTD., No. 8556). The impurity H content in the sc-Ga₂O₃ sample was estimated to be 3.5$\times 10^{18}\mbox{cm}^{-3}$ by thermal desorption spectrometry, which is sufficient to explain the above described $N_{\mathrm{e}}$. Details are described in Supplemental Material (SM) sm .

In this study, $\mu$SR measurements and data analysis were also performed on a powder sample (99.99%, provided by Rare Metallics Co.) for comparison, and the results were found to be significantly different from those for single crystal. However, these results are excluded from the discussion in this paper, because it is difficult to measure the bulk electronic properties of powder samples which is necessary to consider the cause of the difference. Instead, they are presented with a brief interpretation in SM sm .

Conventional $\mu$SR measurements were performed using the S1 instrument (ARTEMIS) at the Materials and Life-science Experiment Facility, J-PARC Kojima et al. (2014), where high-precision measurements over a long time range of 20–25 $\mu$s can be routinely performed using a high-flux pulsed muon beam ($\sim$$3\times 10^{4}$ $\mu^{+}$/cm²/s for the single-pulse mode at a proton beam power of 0.8 MW). The $\mu$SR spectra [the time-dependent decay-positron asymmetry, $A(t)$] which reflects the magnetic field distribution at the Mu site, was measured from room temperature to 4 K under zero field (ZF), weak longitudinal field (LF, parallel to the initial Mu polarization ${\bm{P}}_{\mu}$), and weak transverse field (TF, perpendicular to ${\bm{P}}_{\mu}$), and were analyzed by least-squares curve fitting Suter and Wojek (2012). The background contribution from muons which missed the sample was estimated from $\mu$SR measurements on a holmium plate of the same geometry and subtracted from the asymmetry.

Since the $\mu$SR spectra were found to be dominated by signals from the diamagnetic Mu (i.e., Mu⁺ or Mu^-), the data under ZF/LF conditions were analyzed using the dynamical Kubo-Toyabe (KT) function, $G_{z}^{\mathrm{KT}}(t;\Delta,\nu,B_{\mathrm{LF}})$, where $\Delta$ is the linewidth determined by random local fields from nuclear magnetic moments, $B_{\mathrm{LF}}$ is the magnitude of LF, and $\nu$ is the fluctuation rate of $\Delta$ Hayano et al. (1979). The KT function is expressed analytically in the case of static ($\nu=0$) and ZF ($B_{\mathrm{LF}}=0$) conditions,

The magnitude of $\Delta$ for a given Mu site is evaluated by calculating the second moments of dipolar fields from nuclear magnetic moments by the following equation,

where $\gamma_{\mu}/2\pi=135.539$ [MHz/T] is the muon gyromagnetic ratio, ${\bm{r}}_{mj}=(x_{mj},y_{mj},z_{mj})$ is the position vector from the muon site to the $j$-th nucleus, ${\bm{\mu}}_{m}=\gamma_{m}{\bm{I}}_{m}$ is the nuclear magnetic moment of the atom with the natural abundance of $f_{m}$, $\Theta_{j}$ is the polar angle of ${\bm{r}}_{mj}$ Hayano et al. (1979). Since ⁶⁹Ga with $f_{1}=0.604$ and ⁷¹Ga with $f_{2}=0.396$ have nuclear spins of $I_{m}=3/2$, respectively, they have electric quadrupole moments. In this case, ${\bm{\mu}}_{m}$ is subject to electric quadrupolar interactions with the electric field gradient generated by the point charge of the diamagnetic Mu, leading to the reduction of the effective ${\bm{\mu}}_{m}$ (by a factor $\sin^{2}\Theta_{j}$) and to the modification of $\Delta$ that depends on the initial muon spin direction relative to the crystal axis. (In the powder sample, the spatial averaging leads to $\langle\sin^{2}\Theta_{j}\rangle=2/3$.) The $\Delta$ of the candidate Mu sites inferred from DFT calculations were evaluated using the Dipelec code Kojima with Eq. (2) implemented.

DFT calculations were performed to investigate in detail the local structure of H-related defects using OpenMX Ozaki (2003) and VASP codes Kresse and Hafner (1993). Structural relaxation calculations with H (to mimic Mu) using the GGA-PBE exchange correlation function were performed on a $1\times 3\times 2$ superlattice with cutoff energies of 200 Ry for OpenMX (520 eV for VASP), and $K$ points were set to $3\times 4\times 3$. Structures were relaxed until the the maximum force on each atom was less than $3\times 10^{-4}$ (Hartree/Bohr) for OpenMX (0.01 eV/Å for VASP). Additional structural relaxation calculations using the Heyd-Scuseria-Ernzerhof (HSE06) hybrid functional implemented in VASP were performed with $K$ points of $2\times 3\times 2$ to confirm the structure of the minimum formation energy Krukau et al. (2006). The Hartree-Fock mixing parameter was set to 0.35, which reproduces the experimental value of the band gap Varley et al. (2010, 2011).

First, we investigated the initial asymmetry [$A(0)$] of the $\mu$SR spectrum under TF = 2 mT and found that it is independent of temperature ($T$) and nearly constant within experimental error over the entire measured temperature range $4\leq T\leq 300$ K [see Fig. S5(a) in SM sm ]. Since the absence of a paramagnetic Mu at 300 K is confirmed from the TF data, the $T$-independent $A(0)$ values indicates that the muon is mostly in diamagnetic state(s) at all temperatures. (For the possibility of Mu⁰ with extremely small hyperfine parameters to exist, see Section IV.)

The normalized $\mu$SR time spectra [$G_{z}(t)\equiv A(t)/A(0)$] at 300 K under ZF and LF ($\parallel\langle 001\rangle$) are shown in Fig. 1(a). The ZF-$\mu$SR spectrum shows that it consists of two components, one we call Mu₁, which partially recovers at $t\gtrsim 12$ $\mu\mathrm{s}$ following slow Gaussian relaxation, and Mu₂, which exhibits almost no relaxation. Note that the recovery of $G_{z}(t)$ characteristic to the 1/3-tail of the Kubo-Toyabe function in Eq. (1) is clearly visible only for $t\gtrsim$15 $\mu$s. This could not have been detectable by previous $\mu$SR measurements under TF alone King et al. (2010) or those using a continuous muon beam Celebi et al. (2012), and has been first revealed by precise ZF-$\mu$SR measurements using the high-flux pulsed muon beam at J-PARC. The suppression of the Gaussian relaxation by the weak LF indicates that the relaxation for the Mu₁ component is induced by the quasi-static random local fields from Ga nuclear magnetic moments.

It is clear from Fig. 1(b) that the relative yield of Mu₁ vs Mu₂ apparently depends on $T$ at lower temperatures. Furthermore, the lineshape exhibits change from a Gaussian to exponential-like behavior at low temperatures [see $G_{z}(t)$ for $t\lesssim 6$ $\mu$s at 30 K]. Considering these features, we analyzed the ZF and LF spectra by global curve-fits using the following function,

where $f_{i}$ ($i=1,2$) are the relative yields of the Mu_i components. The exponential damping with the relaxation rate $\lambda$ is introduced to describe the possible influence of the fluctuating magnetic fields from unpaired electrons (including excited carriers and excitons localized nearby Mu). As shown by the solid lines in Figs. 1(a) (b), curve fit using Eq. (3) provides reasonable agreement with the data, and $\Delta_{1}=0.140(1)$ MHz is obtained from the fit at 300 K. The $T$ dependencies of $\Delta_{1}$ and $\nu_{1}$ obtained from the fit are shown in Figs. 1(c) and (d), and those of $f_{2}$ and $\lambda$ in Figs. 2(a) and (c), respectively.

Above $T_{\gamma}\simeq 80$ K where the $f_{1}$ ($=1-f_{2}$) is nealy independent of $T$, $\Delta_{1}$, $\nu_{1}$ and $\lambda$ also show similar trends; $\Delta_{1}$ increases slightly to $\sim$0.15 MHz as $T$ decreases, while $\nu_{1}$ and $\lambda$ shows a constant value close to zero ($\ll\Delta_{1}$). This indicates that the Mu₁ component is quasistatic for $T\gtrsim T_{\gamma}$. Meanwhile, the rapid motion of Mu₂ is inferred from the fact that there is no interstitial sites free of local magnetic fields from Ga nuclear magnetic moments (with 100% natural abundance); $\Delta$ calculated at any site is always larger than $\sim$0.15 MHz for the unrelaxed lattice, which is only compatible with the situation where the relaxation due to the local field (with a linewidth $\Delta_{2}\gtrsim 0.15$ MHz) is suppressed by the motional averaging, namely,

with $G^{\rm KT}_{z}(t;\Delta_{2},\nu_{2})\simeq 1$ for $\nu_{2}\gg\Delta_{2}$.

Below $T_{\gamma}$, $\nu_{1}$ increases with decreasing $T$ in correlation with $f_{2}$, and tends to approach a constant value or decrease for $T\lesssim T_{\alpha}\simeq 30$ K. A similar correlation is observed between $f_{2}$ and $\lambda$; $\lambda$ gradually increases with decreasing $T$ to exhibit a small peak around $T_{\gamma}$, then decreases towards $T_{\beta}$, followed by an increase to reach another maximum near $T_{\alpha}$. Here, we point out that, as shown in Fig. 2(b), the complicated $T$ dependence of $f_{2}$ and the increase of $\lambda$ occur in correlation with the decrease of $N_{\mathrm{e}}\mu_{\rm e}$ (with $\mu_{\rm e}$ being the mobility) of the sample used for $\mu$SR measurements. More specifically, $N_{\mathrm{e}}$ exhibits gradual decrease around $T_{\gamma}$ with decreasing $T$, followed by a sharp decrease below $T_{\beta}$ and an increase below $T_{\alpha}$ (see Fig. S2 in SM sm ). As shown in Fig. 2(d), $\mu_{\rm e}$ decreases with decreasing $T$ below $\sim$120 K, reaching a minimum around $T_{\beta}$, and then increases to a maximum at $T_{\alpha}$. The behavior of $\mu_{\rm e}$ for $T>T_{\gamma}$ is qualitatively in line with previous reports Parisini and Fornari (2016); Kabilova et al. (2019). This implies that the electronic state of Mu is highly sensitive to the quality of samples, which is also supported by the fact that the $T$ dependence of $f_{2}$ (and $f_{1}$) in powder sample differs significantly from that in single crystal (see Fig. S5(b) in SM sm ).

In the ambipolarity model, the electronic state of the interstitial Mu is determined by where the donor and acceptor levels associated with H, predicted from the Fermi energy ($E_{F}$) dependence of defect formation energy ($\Xi^{q}$) for H^q ($q=\pm,0$), is located in the energy band structure of the host material. Figure 3(a) shows $\Xi^{q}(E_{F})$ for each valence state of H obtained by previous DFT calculations Li and Robertson (2014), where H is substituted with Mu. The relationship between the donor/acceptor levels and band structure predicts that Mu can take two different diamagnetic states: a donor-like state in which it releases an electron into the conduction band to become Mu⁺, and a hydride (Mu^-) state, which is qualitatively consistent with the present observations. Interestingly, the $E^{+/-}$ level, which determines the charge state of H at thermal equilibrium, is in the conduction band as is the donor level, so the behavior of H is expected to be the same as that of the donor-like Mu.

To determine the local defect structure of the Mu₁ state, $\Delta_{1}$ was compared with those calculated for the candidate sites, and the results are summarized in Table  1. The value for the H_I site bonded to the three-coordinated O_I [see Fig. 4(a)] is in good agreement with $\Delta_{1}$ at 300 K, indicating that H${}_{\rm I}^{+}\approx{\rm Mu}_{1}^{+}$ in the dilute limit ($\sim$10⁵ cm^-3). This confirms the earlier point that, aside from not knowing the exact location of $E^{+/-}$, H_I can serve as donor Varley et al. (2011). The existence of H in the H_I structure hss been also inferred from the infrared spectroscopy of hydrogenated $\beta$-Ga₂O₃ Qin et al. (2019). The possibility that Mu₁ corresponds to H trapped in Ga and O vacancies is excluded by the comparison of $\Delta_{1}$ with those for Mu in these vacancies (shown in Table S1 in SM sm ).

As seen in Fig. 1(c), $\Delta_{1}$ exhibits a slight increase with decreasing $T$, approaching the value at the H_II/H_III sites [see Fig. 4(b), (c)]. This can be interpreted as reflecting the fact that the initial population of relaxed-excited states immediately after muon implantation is a random sampling of available metastable sites. Considering that the difference in the formation energy estimated by DFT calculations is small among these sites (e.g., $\sim$30 meV between H_I and H_II), the decrease in $\Delta_{1}$ with increasing $T$ above $T_{\gamma}$ suggests that the Mu site distribution approaches the lowest energy state by the annealing process. A similar $T$ dependence of $\Delta$ exhibited by donor-like Mu has been reported for InGaZnO₄ in which multiple H sites are available Kojima et al. (2019).

With Mu₁ identified as being in the donor-like state, the paired Mu₂ component is presumed to be in the acceptor-like diamagnetic state (Mu${}_{2}^{-}$) from the ambipolarity model. According to DFT calculations, the corresponding H^- state is bonded to the two nn Ga ions, as shown in Fig. 4(d) Li and Robertson (2014). The $E^{0/-}$ level is near the conduction band, and it can still accommodate two electrons as it is located below the charge neutral level ($E_{F}^{\rm int}$, see Fig. 3). This means that it is possible for Mu₂ to exchange electrons with the conduction band. Moreover, the conduction band has a relatively large dispersion (effective mass $\sim$0.28$m_{e}$) around the $\Gamma$ point Peelaers and Van de Walle (2015), allowing fast migration of electrons.

While Mu${}_{2}^{-}$ is unlikely to exhibit fast diffusion by itself due to the bonding to cations, there are many reported examples of fast diffusive motion of acceptor-like Mu⁰ in semiconductors and alkali halides Patterson (1988); Kadono et al. (1994); Gxawu et al. (2005); Chow et al. (2000). Since $N_{\rm e}\mu_{\rm e}$ and $f_{2}$ show almost the same $T$ dependence for $T>T_{\beta}$, it is inferred that the Mu₂ state diffuses rapidly for $T>T_{\gamma}$ through the neutral state temporally attained by the charge exchange reaction Mu${}_{2}^{-}\rightleftarrows{\rm Mu}_{2}^{0}+e^{-}$. The reaction rate is given by

where $r_{-0}$ ($=1/\tau_{-0}$) is the ionization probability ($\propto e^{-\varDelta E/k_{B}T}$, with $\varDelta E\sim E_{c}-E^{0/-}$ and $E_{c}$ the CBM energy) and $r_{0-}$ ($=1/\tau_{0-}$) is the capture rate of the free carriers ($\propto\sigma_{\rm c}N_{\mathrm{e}}\mu_{\rm e}V_{r}$, with $\sigma_{\rm c}$, $\mu_{\rm e}$, and $V_{r}$ being the cross section, carrier mobility, and electric field exerted from the Mu-Ga complex state, respectively). Note that $\lambda$ due to the charge exchange is also quenched by the motional averaging when $r_{\mathrm{ex}}$ is much greater than the hyperfine parameter. Such motional effects have also been observed in SiO₂ and are thought to be responsible for the disappearance of anisotropy in the hyperfine interaction of Mu⁰ with increasing $T$ Chow et al. (2000). The existence of the non-relaxing Mu component has also been reported for variety of materials including NaAlH₄, LaScSiH_x, and FeS₂ Kadono et al. (2008); Okabe et al. (2018); Hiraishi et al. (2021), implying the ubiquitous nature of acceptor-like Mu states.

It is natural to assume that Mu₂ diffuses along the $\langle 010\rangle$ direction because $\beta$-Ga₂O₃ has an open channel structure along the $\langle 010\rangle$ direction, In fact, secondary ion mass spectrometry analysis of ²H implanted $\beta$-Ga₂O₃ Reinertsen et al. (2020) revealed that ²H diffuses more easily in the $\langle 010\rangle$ direction than perpendicular to the $\langle-201\rangle$ surface. It has also been reported for rutile TiO₂ that hydrogen diffusion is more likely to occur in open channels along the $c$-axis than in the $a$-$b$ plane Bates et al. (1979).

There are three different open channels in the $\langle 010\rangle$ direction in $\beta$-Ga₂O₃. The potentials associated with diffusion in each channel were investigated by total energy calculations (see Sect. III in SM sm for details), and it was found that hydrogen tends to diffuse in the same channel as the Mu₁ site. Since the two charged states Mu${}_{1}^{+}$ and Mu${}_{2}^{-}$ accompany mutually different lattice distortions, it is assumed that the potential barriers caused by these distortions allow Mu₂ to diffuse while maintaining its state. However, because the adopted calculations do not take structural relaxation into account, it is a future task to evaluate using the nudged-elastic band (NEB) method.

The details of what causes the decrease in $f_{2}$ (increase in $f_{1}$) and recovery below $T_{\gamma}$, which is apparently correlated with the onset of slow dynamics for Mu₁ (as inferred from the $T$ dependence of $\Delta_{1}$, $\nu_{1}$, and $\lambda$), are currently unknown. Nevertheless, when analyzing the time spectra in this $T$ range using Eq. (4) with $\Delta_{1}$ and $\Delta_{2}$ as free parameters, $\Delta_{i}$ converge to mutually close values ($\sim$0.13–0.14 MHz). A curve fit performed with $\Delta_{1}$ fixed at 0.136 MHz and $\Delta_{2}$ at 0.345 MHz (corresponding to the stationary hydride state) and the other parameters free fails to reproduce the data, as indicated by a large reduced chi-square (see Fig. S6 in SM). Thus, the decrease in $f_{2}$ suggests that a transition from Mu₂ to Mu₁ occurs in the relevant $T$ range. One possible scenario is that Mu${}_{2}^{-}$ can be ionized by thermal excitation and transitions to Mu${}_{1}^{+}$ for $T<T_{\beta}$, but the ionization is prevented by increasing $N_{\rm e}\mu_{\rm e}$ for $T>T_{\beta}$ [i.e., $r_{0-}\gg r_{-0}$ in Eq. (5)], allowing it to continue to exist as a metastable state again above $T_{\gamma}$.

Finally, we comment on the relationship between our experimental results and the prior $\mu$SR studies on $\beta$-Ga₂O₃. In Ref. [King et al., 2010], TF-$\mu$SR measurements are reported for a powder sample from a commercial vendor, and it is observed that about 10% of the signal exhibits slight increase of relaxation rate from 0.08 MHz to 0.12 MHz below 50–100 K. Ref. [Celebi et al., 2012] reports results for a single-crystal sample (orientation unknown) and observes a 4–6% decrease in the initial asymmetry below 50–100 K in both ZF and TF-$\mu$SR measurements and an increase in $\Delta$ of the component described by $G_{z}^{\rm KT}(t)$ from 0.10–0.11 MHz to 0.15 MHz in the ZF-$\mu$SR spectra. Although detailed comparisons with our results are not possible because no time spectra are available for either case, these changes appear to correspond to the increase in $f_{1}$ at $T<T_{\gamma}$ that we observed in both the powder and single-crystalline samples [see Fig. S5(b) in SM sm ]. Meanwhile, these reports are silent about the component corresponding to Mu₂. We stress that the existence and origin of Mu₂ are clearly identified for the first time in this study with the help of the ambipolarity model.

The relaxed structures obtained by DFT calculations are shown in Fig. S3, and the simulated muon spin relaxation rates $\Delta$ in the Kubo-Toyabe function for each structure are shown in Table S1.

To investigate the Mu₂ diffusion path, $\langle 010\rangle$ axis position $y$ dependence of the total energy $E_{\mathrm{tot}}$ was calculated in a $1\times 2\times 1$ superlattice. Figure S4 shows the variation of the total energy,

where $\bm{r}$ is the position vector of channel $c_{i}$ $(i=1,2,3)$ and $E_{\mathrm{min}}$ is the global minimum of $E_{\mathrm{tot}}(\bm{r})$. Here, ${\bm{r}}$ for $c_{1}=(0.25,y,0)$, $c_{2}=(0.5,y,0.5)$ and $c_{3}=(0.5,y,0)$, respectively. The $\Delta E_{\mathrm{tot}}$ of $c_{2}$ and $c_{3}$ show large changes of more than 0.5 eV, whereas the $\Delta E_{\mathrm{tot}}$ of $c_{1}$ is almost independent of $y$ and the difference is only 0.04 eV, suggesting that $c_{1}$ is a the possible diffusion path for Mu₂.