Superconductivity assisted change of the perpendicular magnetic anisotropy in V/MgO/Fe junctions

Controlling the perpendicular magnetic anisotropy (PMA) in thin films has received considerable attention in recent years due to its technological importance. PMA based devices usually involve heavy-metal (oxide)/ferromagnetic-metal bilayers, where, thanks to interfacial spin-orbit coupling (SOC), the in-plane (IP) stability of the magnetisation is broken. Here we show that in V/MgO/Fe(001) epitaxial junctions with competing in-plane and out-of-plane (OOP) magnetic anisotropies, the SOC mediated interaction between a ferromagnet (FM) and a superconductor (SC) enhances the effective PMA below the superconducting transition. This produces a partial magnetisation reorientation without any applied field for all but the largest junctions, where the IP anisotropy is more robust; for the smallest junctions there is a reduction of the field required to induce a complete OOP transition (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_\text {OOP}$$\end{document}HOOP) due to the stronger competition between the IP and OOP anisotropies. Our results suggest that the degree of effective PMA could be controlled by the junction lateral size in the presence of superconductivity and an applied electric field. We also discuss how the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_\text {OOP}$$\end{document}HOOP field could be affected by the interaction between magnetic stray fields and superconducting vortices. Our experimental findings, supported by numerical modelling of the ferromagnet-superconductor interaction, open pathways to active control of magnetic anisotropy in the emerging dissipation-free superconducting spin electronics.


Correction of the OOP field misalignment
In order to ensure that the magnetic field is perfectly aligned with respect to the FM layers, we calibrate the angle between the V/MgO/Fe plane and the magnetic field created by the vector magnet superconducting coils by performing sub-gap conductance measurements at T = 0.3 K for different field directions around each axis. Figure S1 describes the calibration process in detail. The results are robust throughout the studied samples, with a misalignment of 8 ±1 degrees with respect to the X-axis superconducting coil, which is accounted for in the OOP experiments. There is no observed in-plane misalignment (Y and Z axis coils).
For the highest OOP-AP transition fields measured (about 2000 Oe as shown in Figure 1c in the main text), the uncertainty of 1 degree in the misalignment calibration for the X coil could result in an IP component of the applied field of 35 Oe, which is about 20 times smaller than the IP coercive field of the hard FM layer 1 . Therefore, undesired IP field due to misalignment can be ruled out as an explanation for the observed AP transitions.

Numerical simulations
Micromagnetic simulations were carried out using the MuMax3 2 software in order to estimate (i) the possible influence of Meissner effect on the OOP transition dynamics, (ii) the role of magnetic inhomogeneities caused by defects on the volatility of the OOP magnetic state (characterized by the OOP to IP transition field, H OOP-IP ) and (iii) the OOP reorientation dynamics. Additional simulations of the OOP transition were performed varying the lateral size of the simulated samples, in order to contrast the results with the observed enhancement of the OOP transition field H OOP with the junctions' lateral size. The following typical Fe magnetic parameters were implemented on the system: saturation magnetization M S = 2.13 T (1.7 × 10 6 A/m), exchange stiffness A exch = 2.1 × 10 −11 J/m, damping α = 0.02 and first order cubic anisotropy K C1 = 4.8 × 10 4 J/m 3 . All simulations were made with T = 0 K. The simulation cells of the system were set at 128 for the IP components and at 16 for the OOP component. The lateral size of the simulated samples was 0.4 × 0.4 µm 2 , with a thickness of 10 nm. The system was additionally surrounded by a vacuum box of an added 100 nm (50 nm on each side) and 3 nm on the top and bottom of the Calibration measurements for the X axis SC coil in a 20 × 20 µm 2 junction at T = 0.3 K. (a) shows the conductance measured at V = 0.3 mV, inside the SC gap, for different values of the ϕ angle (defined in the inset, which is a sketch of the sample holder situation with respect to the three SC coils). (b) shows two conductance curves at H = 0 and H OOP = 300 Oe. The dashed line indicates the applied voltage during the calibration process. When an OOP magnetic field is present, the SC gap diminishes and the conductance increases. The misalignment angle is therefore obtained as the one which minimizes the conductance (ϕ = 82 deg in (a)).
sample, leaving the discretization size of the simulation at 3.9 × 3.9 × 1 nm. Higher discretization in the OOP direction has been chosen on purpose to observe the OOP effects with high accuracy in the simulations. Perpendicular magnetic anisotropy (PMA) was introduced on the top and bottom layers of the Fe as surface anisotropy, with a value of K S1 = 8.32 × 10 −3 J/m 2 . Due to computational limitations on the size and detail of the simulations, the results discussed in this section should be taken as qualitative support for the experimental results, rather than quantitative estimations. The influence of superconductivity on the OOP transition has been studied by performing micromagnetic simulations with MuMax3 for 10 nm thick, 0.4 × 0.4 µm 2 Fe(001) films under the influence of a superconducting vanadium layer. We simulated OOP hysteresis cycles where a correction to the applied field was added based on a typical Meissner effect (ME) hysteresis cycle (obtained from 3 , shown in Figure S2b inset), scaled for different values of field contribution from Meissner effect and adapted to the first and second critical fields of vanadium (correspondingly H c1 and H c2 ). The contribution from superconducting vortices was taken into account by using an in-group developed program that numerically solves the time dependent Ginzburg-Landau equations in order to simulate the behaviour of type II superconductors under magnetic fields 4 . The initial stray fields from an in-plane saturated FM simulation were used to generate a distribution of vortices, and then the fields generated by those vortices were calculated 5 and added into the corrected hysteresis cycle. Our simplified numerical model, although limited, provides qualitative support for the mutual magnetostatic interaction between the FM and SC as the possible origin of the behaviour of the H OOP field in the superconducting state. We also note that simulations with a contribution exceeding 7% of Meissner effect resulted in the OOP state being volatile in the hysteresis cycle (i.e. the magnetization returns to an IP configuration before returning to zero field), in contradiction with the experimental observations. Consequently, we have not considered larger contributions of Meissner effect as a possible explanation for the observed behaviour of H OOP . We also underline that a complete numerical solution is a great challenge which is outside our current capabilities, as the problem should be solved self-consistently, so the results should be understood in a qualitative way.  . An asterisk has been manually added with the transition field for the simulation with no defects. For simulations with more than 1% of defects, the OOP-IP transition becomes volatile (i.e. it happens before the field changes from positive to negative), in contradiction with the experimental results. The solid line is a guide for the eye, while the inset shows an example image of the surface defects introduced.

Influence of Meissner effect on the OOP magnetization reorientation
Another concern about the OOP reorientation is the possibility of it being a trivial effect produced at the edges of the FM layers. As mentioned in the main text, there is experimental evidence pointing against this possibility. However, a more thorough study has been performed in order to fully discard this scenario. First, we simulated the field-induced OOP reorientation in 0.4 × 0.4 micron Fe films. As shown in figure S3a, the reorientation seems to be triggered by OOP oriented domains in the 3/7 interior of the film.
On the other hand, we wanted to see the influence of any possible dominating edge effects on the underlying SC layer. Simulations of 0.4 × 0.4 micron V films were made using the code for type II superconductors described above at T = 2 K, with OOP magnetic fields applied at the edges of the simulated samples. As shown in figure S3b, the stationary state is reached when vortices fill the interior of the film. This would produce high OOP stray fields affecting the interior of the Fe layer, triggering a reorientation in the whole film rather than limiting it to the edges.

Influence of defects on the OOP-IP magnetization reorientation
In order to better understand the experimentally observed non-volatility of the OOP state in the junctions, we have simulated numerically the influence of randomly distributed surface magnetization defects within the bottom and top layers of the 10 nm thick Fe layer. The defects are introduced in the simulations as spots of enhanced surface saturation magnetization (M S (defects) = 1.25 × M S (Fe)). We have found that the introduction of a small number (about 10 −3 %) of magnetic defects per layer does not affect the non-volatity, but only varies the characteristic field H OOP-IP of the transition from the OOP state to the IP alignment that takes place after the initial magnetization saturation ( Figure S4). Above some critical defect number of about 2 × 10 −3 % defects per layer the OOP-IP transition becomes volatile. As long as we always observe experimentally non-volatility of this transition in our junctions, we can conclude that there is a relatively small number of magnetic defects present in the epitaxial MTJs under study. shows the transition field for micromagnetic simulations of different lateral sizes. Note that, due to computational limitations, the simulated sizes are smaller than the actual samples. However, the trend is in qualitative agreement with the experimental results. The lines are linear fittings of the experimental/simulation points, which should serve as guides for the eyes. This trend is accomplished both with or without considering periodic boundary conditions (PBC).

Dependence of H OOP with the junctions lateral size
The behaviour of the OOP transition field H OOP has been studied as a function of the samples lateral size, in a total of 14 different samples varying from 10 × 10 to 40 × 40 µm 2 . We found an increasing trend of H OOP with the lateral size, as shown in Figure S5a. Micromagnetic simulations of the same transition have been performed in Fe films with lateral sizes of 0.1 × 0.1 to 0.4 × 0.4 µm 2 (as the real dimensions were computationally prohibitive to simulate) with a qualitative agreement to the experimental results, as shown in Figure S5b. These simulations have been made in the absence of previously discussed phenomena such as defects or Meissner effect.

Bogoliubov-de Gennes theoretical model
In order to demonstrate how the superconducting proximity effect can cause a change in the magnetic anisotropy, we consider a tight-binding Bogoliubov-de Gennes model for the V/MgO/Fe heterostructure. This model enables us to calculate the free energy of the system, so that we can study how the energy cost for reorienting the magnetization changes with temperature. We here only focus on the superconductivity-induced decrease in H OOP observed for 10 × 10 µm junctions. Vortex formation and the size of the lattice is not taken into account here, and instead discussed in section S2.1 and in the microscopic model section The first term describes the nearest-neighbor hopping, where t is the hopping integral. The second term describes the chemical potential µ i at each lattice site i, and the potential barrier V i > 0 present in the insulating MgO layers. The third term gives rise to an attractive on-site interaction in the superconducting V layer described by the onsite potental U i > 0. The fourth term introduces a local magnetic exchange field h i giving rise to ferromagnetism in the Fe layer. The Pauli matrices are contained in the vector σ . The last term describes the Rashba spin-orbit coupling boosted by the MgO layers, where the spin-orbit field has a magnitude λ i and is directed along the interface normal n. The vector d i, j connects site i and j. In the above Hamiltonian, c † i,σ and c i,σ are the second-quantization electron creation and annihilation operators at site i with spin σ , and n i,σ ≡ c † i,σ c i,σ is the number operator. The superconducting term is treated by a mean-field approach assuming that c i,↑ c i,↓ ≈ c i,↑ c i,↓ + δ , where terms to the second order in the fluctuations δ are negligible. The superconducting gap is defined as ∆ i ≡ U i c i,↑ c i,↓ and must be treated self-consistently. The above model is valid in the ballistic limit, but since the effects considered here depend on the formation of s-wave odd-frequency triplets that are robust to impurity scattering we would obtain qualitatively the same results in the diffusive limit.
We consider a cubic lattice of size N x × N y × N z with an interface normal along the x axis. We assume periodic boundary conditions in the y and z directions, and apply the Fourier transform along these axes. To simplify notation we have defined i ≡ i x , j ≡ j x , i || = (i y , i z ), and k ≡ (k y , k z ). We also use that The Hamiltonian can be written on the form where the basis is given by and where the Hamiltonian matrix H k consists of N x × N x blocks with row and column indices (i, j). Above,τ iσ j ≡τ i ⊗σ j is the Kronecker product of the Pauli matrices spanning Nambu and spin space,τ ± ≡ (τ 1 ± iτ 2 )/2, and The constant term is given by By diagonalizing H k , we obtain eigenvalues E n,k and eigenvectors The diagonalized Hamiltonian can be written as

5/7
where the marked sum goes over {n, k y , k z > 0}, {n, k y > 0, k z = 0, −π}, and {n corresponding to E n,k y ,k z > 0, k y = 0, −π, k z = 0, −π}. Expectation values of the new operators can now be evaluated according to where f (E n,k ) is the Fermi-Dirac distribution. The new quasi-particle operators are related to the old operators by The eigenenergies E n,k and eigenvectors Φ n,k obtained in this diagonalization, can be used to calculate physical observables for the system. The superconducting gap is given by and is treated self-consistently. We can calculate the critical field H OOP for reorienting the magnetization from an IP to an OOP orientation. The Zeeman energy of an external magnetic field H is given by where µ 0 is the vacuum permeability, µ tot is the total magnetic moment, and H is the applied field. If we consider a system where the free energy is minimal for an IP magnetization and maximal for an OOP magnetization, and we want to find the external magnetic field needed to reorient the magnetization to the OOP direction, we must require that |F Zeeman | ≥ F OOP − F IP . We can then calculate the critical field from To take into account other anisotropy contributions not covered by this model, we let F OOP → F OOP + K anis . Above, the free energy is given by where β = (k B T ) −1 . The total magnetic moment of the system for an OOP magnetization is given by Re(u * i,n,k v i,n,k ) f (E n,k ) + Re(x * i,n,k w i,n,k ) 1 − f (E n,k ) , when the interface normal is directed along the x axis. Since the lattice must be scaled down in order to make the system computationally manageable, we choose the magnitude of the on-site coupling potential U i so that the superconducting coherence length is comparable to the thickness of the V layer. The superconducting coherence length is given by ξ =hv F /π∆ 0 , where v F = (1/h)dE k /dk k=k F is the Fermi velocity calculated for the normal-state eigenenergy E k = −2t[cos(k x ) + cos(k y ) + cos(k z )] − µ, and ∆ 0 is the zero-temperature superconducting gap.
We determine the superconducting critical temperature by a binomial search, where we decide if a given temperature is above or below T C . This is decided by finding whether the superconducting gap measured in the middle of the superconducting region increases towards a superconducting solution or decreases towards a normal state solution from the initial guess ∆ ∆(T = 0) under iterative recalculations.
In Fig. 4 in the main text, we have used the parameters t = 1, µ i∈S = 0.9, µ i∈SOC = µ i∈F = 0.8, V i∈SOC = 0.79, U = 1.4, λ = 0.4, h = 0.8, N S x = 28, N SOC x = 3, N F x = 8, and N y = N z = 50. This gives a coherence length of 21 lattice sites. All length scales are scaled by the lattice constant a, all energy scales are scaled by the hopping parameter t, and the magnitude of the spin-orbit coupling is scaled by ta. 6/7