Tailoring the anomalous Hall effect of SrRuO₃ thin films by strain: A first principles study

The anomalous Hall effect (AHE) plays a crucial role in modern condensed matter physics and material science research, as it often yields direct insight into fundamental physical properties.^1–3 In the absence of noncollinear magnetic order, the behavior of the Hall resistivity in an external magnetic field $H$ is expressed as $ρxy=R0H+RsM$⁠, where $R0$ and $Rs$ are ordinary and extraordinary Hall coefficients, respectively, and $M$ is the magnetization of the sample. The second term in the above equation is known as anomalous Hall contribution. The additional hump-like anomalies appearing on top of the expected behavior of Hall resistivity with increasing magnetic field are often interpreted as a signature of the formation of non-collinear magnetic structures, including topological chiral particles such as magnetic skyrmions, contributing to the Hall signal via the so-called topological Hall effect. The latter phenomenon, perceived as one of the manifestations of chirality-sensitive Hall effect of textures,^4–6 has gained considerable attention in the past few years as it allows for probing the existence of chiral spin textures by purely electrical means.^7,8 In this context, after the prediction of interface-stabilized skyrmion formation in SrRuO$3$/SrIrO$3$ heterostructures,^9,10 thin films of SrRuO$3$ (SRO)—a material that is historically of great general interest in the context of spintronics applications^11–13—have drawn an immense activity in the past few years aimed at the observation of skyrmion phase by the means of magneto-transport.

In particular, additional hump-like anomalies on top of the expected AHE signal were reported recently in ultrathin films of ferromagnetic SrRuO$3$⁠.^14–17 It was postulated that RuO$6$ octahedra-titling induced by the local orthorhombic-to-tetragonal structural phase transition at the SrRuO$3$/SrTiO$3$ interface breaks structural inversion symmetry and, mediated by symmetry lowering and intrinsic strong spin-orbit coupling of $4d$ SRO states, leads to a finite Dzyaloshinskii–Moriya interaction and stabilization of the skyrmion phase.¹⁵ However, this picture of skyrmion phase stabilization in the system was questioned by several works.^18–22 Very recently, the behavior of the Hall signal measured in SrRuO$3$ films deposited on Pr$0.7$Ca$0.3$MnO$3$ (PCMO)²¹—reminiscent of that driven the topological Hall effect—was interpreted as arising from a superposition of Hall effect contributions originated separately in tetragonal and orthorhombic SrRuO$3$ layers. Depending on strain^23–25 or thickness of PCMO,²¹ modifications in the RuO$6$ octahedral tilting were reported for ultrathin ferromagnetic SRO layers, while for larger strain or for thicker PCMO films, crystal symmetry of deposited ultrathin SrRuO$3$ was found to undergo an orthorhombic-to-tetragonal transition associated with a large change in the RuO$6$ octahedral tilting angle. Notably, an opposite sign of the AHE for tetragonal and orthorhombic SrRuO$3$ films was reported, and it was speculated that this was a probable reason for an observed additional peak in the field dependence of the AHE.^21,22

The existing controversy motivated us to study the AHE of ferromagnetic SRO films from first principles theory in order to gain a better understanding of the AHE in this system in relation to its structural properties. In this work, we investigate the electronic and Hall transport properties of ferromagnetic SRO trilayer by tailoring the oxygen-octahedra tilting angle which is driven by the substrate induced compressive strain. By performing first principles density functional theory (DFT) calculations, we find ferromagnetic SRO films grown on SrTiO$3$ (corresponding to a strain of $−0.47%$⁠) to have an orthorhombic (both tilting and rotation) distorted structure, and with an increasing amount of substrate induced compressive strain the octahedral tilting angle is found to be suppressed gradually, with ferromagnetic SRO films grown on NdGaO$3$ (NGO) (corresponding to a strain of $−$1.7$%$⁠) found to be stabilized in the tetragonal distorted structure (with zero tilting) with magnetic moments pointing out of the plane of the film in complete agreement with experimental observations.^24–28 

Our calculations of the AHE predict a positive value of the anomalous Hall conductivity (⁠$+76$ S/cm) at $−1.7%$ strain, whereas the anomalous Hall conductivity is found to be negative (⁠$−$156 S/cm) at $−0.47%$ strain, which is in good quantitative agreement with experimental data.^12,21,29,30 We attribute the found behavior of the AHE to the nodal point dynamics in the electronic structure arising in response to tailoring the oxygen octahedral distortion driven by the substrate induced strain. We also estimated the AHE as a function of reduced magnetization to mimic the experimentally observed temperature dependence of the AHE in the system, achieving good agreement with the experimentally observed trend.^21,29 We believe that our findings contribute significantly to understanding the physics of the anomalous Hall effect in this delicate material and point toward new routes of engineering the anomalous Hall effect properties in complex oxides.

DFT calculations were carried out with two different approaches: the full-potential linearized augmented plane wave (FLAPW) method as implemented in the Jülich DFT code FLEUR³¹ and the plane-wave projected augmented wave (PAW) method as implemented in Vienna ab initio Simulation Package (VASP).^32,33 The structural optimization of the bulk, as well as thin film structures, was carried out using the VASP code maintaining the symmetry of the crystal. The positions of the atoms were relaxed toward equilibrium until the Hellmann–Feynman forces became less than 0.001 eV/Å. The Monkhorst–Pack³⁴ $k$-point mesh of $8×8×6$ was used for structural optimization of bulk SRO. This choice of the $k$-mesh and a plane-wave cutoff of 500 eV were found to provide a good convergence of the total energy.

The Monkhorst–Pack $k$-point mesh of $10×10×1$ was used for the structural optimization of SrO-terminated SRO thin films, consisting of three unit cells of SRO along the $z$-axis, considering different strains. For the structural optimization of the thin film with different compressive strain in the plane-wave basis, we included 20 Å  of vacuum to minimize the interaction between periodically repeated images along the $z$-axis. Then, we carried out the structural optimization of the thin film structure by relaxing the internal positions allowing for tilting and rotation of RuO$6$ octahedra and keeping the in-plane lattice parameters fixed at different strain values.

Using relaxed atomic positions of SRO trilayer, total energy calculations of different structures, the electronic structure calculations including the effect of spin-orbit coupling (SOC) and the AHE calculations were carried out with the film version of the FLEUR code.³¹ For self-consistent calculations with the LAPW basis set, a plane-wave cutoff of $kmax=4.0$ a.u.$−1$ and a $k$-mesh of $12×12$ in the two-dimensional Brillouin zone were found to be sufficient for the convergence of the total energy. The plane-wave cutoff for the potential (g$max$⁠) and exchange-correlation potential (g$max,xc$⁠) was set to 12.6 and 10.5 a.u.$−1$⁠, respectively. The muffin-tin radii for Sr, Ru, and O were set to 2.80, 2.32, and 1.31 a.u., respectively. For calculations of the magnetic anisotropy energy, the effect of SOC was included self-consistently using $24×24$ $k$-points in the two-dimensional Brillouin zone. We used the Vosko–Wilk–Nusair (VWN)³⁵ exchange-correlation functional within the local density approximation (LDA) for self-consistent calculations.

Bulk SRO is stabilized in an orthorhombic crystal structure (Pbnm) below 850 K with GdFeO$3$-type distortion³⁶ characterized by the tilting of the RuO$6$ octahedra in alternate directions away from the $z$-axis and the rotation of the octahedra, as shown in Fig. 1. In the orthorhombic unit cell (Pbnm) of bulk SRO with experimental lattice parameters $a=5.5670$ Å, $b=5.5304$ Å, and $c=7.8446$ Å, the RuO$6$ octahedral rotation produces a distorted pseudo-cubic perovskite structure with the pseudo-cubic lattice parameter $apc=3.93$ Å.³⁷ In Glazer notation,³⁸ the orthorhombic SRO shows a$−$a$−$c$+$ rotation of octahedra. The orthorhombic unit cell contains four formula units (f.u.) of the ideal cubic perovskite. Sr and O(1) atoms occupy the 4c Wyckoff positions, while Ru and O(2) atoms occupy the 4a and 4d Wyckoff positions, respectively.

Bulk SRO is a ferromagnetic metal with a Curie temperature of 160 K and magnetic moment of 1.1–1.7 $μB$/f.u.^39,40 The optimized lattice parameters of the bulk SRO, obtained while keeping the symmetry of the structure fixed, were found to be in good agreement with the previous studies.^39–41 We find the distortion of RuO$6$ octahedra that manifests in unequal bond lengths and deviations of O-Ru-O bond angles away from 90$°$⁠: the optimized tilting angle, $(180−ϕ)/2$⁠, is found to be 10.56$°$ (corresponding to a Ru-O-Ru angle of 159$°$⁠) and the rotation angle, $(90−θ)/2$⁠, is found to be 7.56$°$⁠.

To simulate the thin film structure with different compressive strain, we consider a free standing tri-layer SrO-terminated SrRuO$3$ film grown on different substrate lattice parameters. We take the optimized lattice constant of bulk SRO as the reference value (i.e., having zero strain) and define the compressive strain value with respect to it. Starting from the SRO trilayer film, grown on the SrTiO$3$ (STO) substrate (corresponding to $−$0.47$%$ of compressive strain), strain is increased gradually up to the value of $−$1.7$%$⁠, which corresponds to SRO thin films grown on the NdGaO$3$ (NGO) substrate. For each value of strain, we first optimize the atomic positions, while keeping the symmetry of the structure fixed. When doing so, we observe a gradual suppression of RuO$6$ octahedra’s tilting angle compared to that of the $−$0.47$%$-strain situation, with the trilayer in which RuO$6$ tilting at $−$1.7$%$ of strain is completely suppressed found to be energetically most stable. In the remaining part of the discussion, we focus only on two values of compressive strain: $−$0.47$%$ and $−$1.7$%$⁠, which correspond to SRO thin films grown on SrTiO$3$ and NdGaO$3$ substrates, respectively.

Our calculations indicate that at the strain of $−$0.47$%$ the SRO trilayer with orthorhombic distortion of RuO$6$ octahedra (i.e., with simultaneous tilt and rotation) is by 48 meV/Ru lower than the SRO film with tetragonal distortion (i.e., only with octahedral rotation). By contrast, at $−1.7%$ strain, the tetragonal structure is by 25 meV/Ru lower in energy than the orthorhombically distorted structure, in agreement with the experimental observations.^23–25 For the calculations presented below, we thus assume that at $−$0.47$%$ strain the SRO trilayer has orthorhombic RuO$6$ octahedral distortion (i.e., having an orthorhombic structure), whereas it has a tetragonal distortion of RuO$6$ at $−1.7%$ strain (i.e., having a tetragonal structure). The changes in the octahedral distortion in terms of the compression of bond length along the $z$-axis and tilting angle for the strain values of $−0.47%$ and $−1.7%$ have a crucial impact on the energetic position of Ru-t$2g$ states, thus directly influencing the electronic structure for two cases, as discussed later.

For both cases of strain, i.e., at $−0.47%$ and $−1.7%$⁠, the ferromagnetic (FM) state was found to be more stable as compared to the anti-ferromagnetic or non-magnetic one. Assuming the FM spin structure, in Figs. 2(a) and 2(b) we show the orbitally resolved density of states (DOS) as calculated with LDA, for the SRO trilayer considering the strain of $−0.47%$ and $−1.7%$⁠, respectively. For both cases, states close to the Fermi level are dominated by the Ru-t$2g$ states, which exhibit a finite hybridization with the O-$p$ states predominantly in the lower region of the valence band (not shown in the plots).

We observe a clear difference in the DOS between the two strain cases. Due to the orthorhombic distortion of RuO$6$ at $−$0.47$%$ stain, the degeneracy of the Ru-$t2g$ is lifted, while for the case of $−$1.7$%$ strain, due to complete suppression of octahedral tilting (tetragonal distortion), the degeneracy is lifted only partially, with two degenerate Ru-$dxz/yz$ bands shifted in energy with respect to a single Ru-$dxy$ band. Due to the octahedral elongation along the $z$-axis with the application of compressive strain, $dxz/yz$ states are stabilized by being pushed down in energy leading in turn to the reordering of Ru-t$2g$ states with $dxz/yz$ states followed by the $dxy$ state.

To quantify the energy level splitting of Ru-t$2g$ bands, we compute the energy-level diagram of Ru-$d$ states, employing the technique of maximally localized Wannier functions (MLWFs),^42–45 considering only the Ru-$d$ Hamiltonian constructed out of non-spin-polarized LDA calculations. With the information on the energy level position of Ru-$d$ states, obtained from the real-space representation of the Hamiltonian in the MLWFs basis, we calculate the energy level difference between three Ru-t$2g$ states. In the inset of Figs. 2(a) and 2(b), we show the energy level splitting of Ru-t$2g$ states for the strain of $−0.47%$ and $−1.7%$⁠, respectively. At $−1.7%$ of strain, driven by the tetragonal distortion, Ru-$dxz/yz$ states become degenerate with an energy splitting of 80 meV between $dxz/yz$ and $dxy$ orbitals. On the other hand, at $−$0.47$%$ of strain, all three t$2g$ states are split off by an almost equal amount of 30 meV.

As also evidenced in Figs. 2(a) and 2(b), owing to the orthorhombic distortion of oxygen octahedra, the bandwidth is slightly narrower in the SRO film at $−0.47%$ strain as compared to that at $−1.7%$ strain. This suggests an increase in itinerancy of the SRO film at $−1.7%$ strain, which is consistent with the experimental observation of increased Curie temperature (155 K) compared to that of SRO films at $−0.47%$ strain (150 K). Within the mean-field theory approach, the Curie temperature, T$c$⁠, is proportional to $Δ$⁠, T$c∼Δ$⁠, where $Δ$ is the energy difference between the ferromagnetic and paramagnetic ground state. Using the mean-field theory approach, we estimate the ratio of the T$c$ of SRO films at $−0.47%$ and $−1.7%$ (1.06), finding it to be in close agreement with the experimentally measured value (1.03), indicating that the description of the magnetic ground state within LDA is quite accurate.

Next, we investigate the magnetic anisotropy energy (MAE) of the trilayer SRO films at $−0.47%$ and $−1.7%$ compressive strain. We obtain the MAE from the difference in the total energy of the ferromagnetic state with spins aligned along the crystal axes $x$ (⁠$y$⁠), with the total energy of the system when the magnetization is pointing along the $z$-axis. For the case of compressive strain studied here, irrespective of its magnitude, the out-of-plane direction is found to be the easy axis of the system in complete agreement with experimental observations.^24,25 Magnetic anisotropy energy of the trilayer SRO films at $−0.47%$ strain is found to be $−0.455$ meV/Ru, which is in qualitative agreement with the experimentally measured value.⁴⁶ At $−1.7%$ strain, the calculated value of the MAE constitutes about $−0.67$ meV/Ru. The calculated averaged spin moment at the Ru ion is found to be 1.20 (1.14) $μB$ at $−0.47%$ (⁠$−1.70%$⁠) strain, which is consistent with the low spin state of Ru$4+$ ion (⁠$d4$⁠:t$2g↑3$⁠,t$2g↓1$⁠). The orbital moments of Ru$2+$(⁠$d4$⁠) point in the same direction as the spin moments, which is expected due to more than half-filled Ru-t$2g$ sub-shell.

Having understood the electronic structure and the magneto-crystalline anisotropy energy of the SRO trilayer films at $−$0.47$%$ and $−$1.7$%$ strain, next we proceed to investigate the anomalous Hall effect of the system, motivated in part by the recently observed anomalies in the behavior of the AHE as a function of an applied magnetic field in this material^14,15,21 Here, we assess the intrinsic Berry curvature contribution to the AHE employing the Wannier interpolation technique.⁴⁷ To compute the Berry curvature, we first construct a tight-binding MLWFs Hamiltonian projected from the LDA+SOC Bloch wave functions. Atomic orbital-like MLWFs of Ru-t$2g$ and -e$g$ states were used to construct the minimal tight-binding Hamiltonian, which reproduces the spectrum of the system in a wide energy window around the Fermi energy. In Figs. 3(a) and 3(c), we show the comparison of the ab initio LDA+SOC band structure of SRO trilayer at $−0.47%$ and $−1.7%$ strain, respectively, with that obtained by diagonalization of Ru-$d$ projected Wannier Hamiltonian, finding excellent agreement in the region of $±3$ eV with respect to the Fermi energy for both strain cases. In the latter plots, one can see a clear difference in the band structure along the $Γ$X direction in the Brillioun zone for the SRO trilayer at $−0.47%$ and $−1.7%$ strain. Due to the suppression of the oxygen octahedral tilting at $−1.7%$ strain, the bands remain degenerate along $Γ$X, whereas for $−0.47%$ strain due to strong oxygen octahedral distortion the degeneracy is lifted. For the case of $−1.7%$ strain, we find that the number of majority and minority bands crossing the Fermi energy is larger than that of the SRO film at $−0.47%$ strain.

From the Wannier Hamiltonian, the Berry curvature was calculated according to the linear-response expression as given by

where $Ωn(k)$ is the Berry curvature of band $n$⁠, $ℏv^i=∂H^(k)/∂ki$ is the $i$th velocity operator, and $unk$ and $ϵnk$ are the eigenstates and eigenvalues of the Hamiltonian $H^(k)$⁠, respectively. Using the tight-binding Hamiltonian constructed from atomic orbital-like Wannier functions of Ru-$d$ states for the two strain cases, as shown in Figs. 3(a) and 3(c), we calculate the Berry curvature on a $50×50$$k$-mesh employing an adaptive $5×5$ refinement scheme⁴⁸ at points where the value of the Berry curvature exceeds 50 a.u. These numerical parameters provide well-converged values of the anomalous Hall conductivity (AHC) determined as

where $Ω(k)$ is the sum (for each $k$-point) of Berry curvature of the occupied states and the broadening of 45 meV was used. We have recently shown that the AHC values obtained in such a way for SRO films are stable with respect to the choice of the MLWFs reproducing the band structure in the whole energy window of occupied states.⁴⁹ 

Our calculations of the AHC of the FM SRO trilayer for two values of compressive strain are shown in Fig. 3(b). We observe that the sign of the AHC for the SRO films at $−0.47%$ of strain is negative, and it becomes positive at $−1.7%$ of strain, which is in complete agreement with the experimental observations.^12,21,29 To understand the impact of the octahedral tilting and rotation on the sign and strength of the AHE, we start with a tetragonal distorted (zero tilting) SRO film at $−1.7%$ strain with a positive AHC and gradually introduce the octahedral tilting. With the increasing value of the octahedral tilting that mimics the SRO films at $−0.47%$ strain, the computed AHC displays a gradual decrease from a positive value and changes sign for larger tilting. This indicates that symmetry lowering realized by octahedral distortion in terms of tilting and rotation followed by a corresponding redistribution of the bands plays a crucial role in shaping the strength and sign of the AHC in strained SRO films.

To gain a better insight into the microscopic origin of the observed sign changes in the Hall response for large compressive strain of SRO trilayer films, for which oxygen octahedral tilting is completely suppressed as compared to that of the SRO films at $−0.47%$ strain, we closely investigate the interplay of the Berry curvature with the electronic band structure for the two cases, with and without SOC. The electronic band structure of the SRO trilayer at $−0.47%$ and $−1.7%$ strain with and without SOC is shown in the upper and middle panel of Figs. 4(a) and 4(b), respectively. For the case of the trilayer at $−1.7%$ strain, we find that the number of majority and minority bands crossing the Fermi energy is larger than that of the SRO film at $−0.47%$ strain. In the absence of SOC, nodal points are formed due to the crossing of the majority (orange lines) and minority (gray lines) bands at the Fermi energy, in Figs. 4(a) and 4(b). When SOC is switched on, the degeneracy of the nodal points is lifted as shown with red dots and small gaps in the electronic structure, which give rise to large peaks in the Berry curvature, are formed. The sign and magnitude of the AHC depend on the energies of the gapped nodal points relative to the Fermi energy.²⁹ At the strain of $−0.47%$⁠, large negative contribution of the Berry curvature is induced by the gapped nodal points found below the Fermi energy close to the $Γ$-point, whereas at the strain of $−1.7%$ gapped nodal points lie above the Fermi energy (hole-like states) inducing large positive Berry curvature. And while the total value of the AHC emerges as a result of contributions coming from various parts of the Brillouin zone, as illustrated in the case of nodal points shown in Fig. 4, the sensitive redistribution of the states with respect to the Fermi energy as well as the changes in the details of their hybridization mediated by SOC and symmetry leads to a very large response of the AHC to applied strain in this system.

Finally, we address the behavior of the anomalous Hall conductivity of SRO trilayer for two compressive strain cases, with the magnitude of the magnetization, in order to make an attempt at understanding the experimentally observed temperature dependence of the anomalous Hall conductivity in SRO thin films. Within the physical picture that we assume, the part of the temperature-driven variation of the AHC that we try to account for is driven by the modifications in the magnitude of the exchange splitting of the states and corresponding ferromagnetic magnetization $M$⁠, rather than transversal fluctuations of the local Ru moments.¹² To do so, we optimize the atomic positions of the SRO trilayer at both values of strain with reduced magnetization $M$⁠, while keeping the symmetry of the structure fixed. Considering optimized structures with reduced magnetization for the two compressive strain values, we calculate the intrinsic contribution to the AHE in the way identical to that described above. Here also, we consider the atomic orbital-like Wannier functions, Ru-$t2g$ and e$g$ to construct the tight-binding Hamiltonian, projected from the LDA+SOC Bloch wave functions for each considered magnetization value.

In Fig. 4(c), we present the computed AHC as a function of reduced magnetization away from the “equilibrium” value, obtained by scaling down the magnitude of the exchange field inside Ru atoms, for the compressive strain of $−0.47%$ and $−1.7%$⁠, respectively. At $−$0.47$%$ strain, for large values of the magnetization, up to a magnitude of about 0.9 $μB$⁠, the AHC remains negative, whereas at the strain of $−1.7%$ it stays positive for magnetization values in a similar range. For both cases of compressive strain, the calculated AHC after the initial increase in the absolute value decreases sharply and changes its sign, with an overall behavior reflecting basic features of experimental findings.^12,21,29 The found correspondence between the experimentally observed data and our calculations suggests that in thin films of strained SRO longitudinal thermal fluctuations of Ru spin moments have a drastic impact on the AHE.

In this work, by performing first principles calculations we investigate the effect of compressive strain-induced oxygen octahedral distortion in the electronic structure and anomalous Hall response of the ferromagnetic SRO ultrathin films. We find a strong deformation of the oxygen octahedra (RuO$6$⁠) with an increasing amount of substrate induced compressive strain. The free standing SRO film grown on SrTiO$3$ (in-plane strain of $−0.47%$⁠) is found to have an orthorhombic (both tilting and rotation) distorted structure, and with an increasing amount of substrate induced compressive strain the octahedral tilting angle is found to be suppressed gradually, with SRO films grown on NdGaO$3$ (in-plane strain of $−1.7%$⁠) stabilized in the tetragonal distorted structure (with zero tilting) in complete agreement with the experimental observations.^24,25 This strong modification of the oxygen octahedral deformation with strain leads to the re-ordering of Ru-t$2g$ states, which has a crucial impact on the electronic structure of SRO films and hence on the anomalous Hall response.

From the Berry curvature calculation, we find a positive value of the anomalous Hall conductivity of $+76$ S/cm at $−1.7%$ strain, whereas it is found to be negative (⁠$−156$ S/cm) at $−0.47%$ strain, which is in good quantitative agreement with the experimental data.^21,29 We attribute the observed behavior of the anomalous Hall effect to the nodal point dynamics in the electronic structure arising in response to tailoring the oxygen octahedral distortion driven by the substrate induced strain. We also calculate strain-mediated anomalous Hall conductivity as a function of reduced magnetization obtained by scaling down the magnitude of the exchange field inside Ru atoms to understanding the experimentally observed temperature dependence of the anomalous Hall conductivity. The result is found to be in good agreement with the experimental observations,^12,21,29 indicating a strong impact of longitudinal thermal fluctuations of Ru spin moments on the anomalous Hall effect in this system.

Overall, our findings reveal the strong influence of the oxygen octahedra deformation, arising from the substrates induced strain, on the electronic structure and anomalous Hall response of SRO films. We believe that our results will help to gain a better understanding of the AHE dynamics in this fascinating system in response to its structural properties, thus further motivating the development of a spintronics paradigm whose functionality is based on engineering the topological band structure by controlling the oxygen octahedral deformations.

We acknowledge extensive discussions with Ionela Lindfors-Vrejoiu, Jairo Sinova, and Thomas Lorentz. We acknowledge funding from Deutsche Forschungsgemeinschaft (DFG) through SPP 2137 “Skyrmionics,” the Collaborative Research Center SFB 1238, and Project No MO 1731/10-1. The work was also funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) TRR 173-268565370 (Project No. A11) and TRR 288-422213477 (Project No. B06). Y.M and S.B. acknowledge the DARPA TEE program through Grant No. MIPR# HR0011831554 from DOI. Simulations were performed with computing resources granted by JARA-HPC from RWTH Aachen University and Forschungszentrum Jülich under projects jpgi11, jiff40, and jias1f.

The data that support the findings of this study are available from the corresponding author upon reasonable request.