Temperature increase in STT-MRAM at writing: A fully three-dimensional finite element approach

The writing process in spin transfer torque magnetoresistive random access memories (STT-MRAM) is facilitated by elevated temperatures. In this work we investigate the temperature distribution and development in the free layer (FL) of an STT-MRAM during switching. With our fully three-dimensional (3D) finite element method simulation approach, we numerically solve the heat transport equation coupled to the electron, spin, and magnetization dynamics and demonstrate that the FL temperature is highly inhomogeneous due to the non- uniform magnetization of the FL during switching for an STT-MRAM with large magnetic tunnel junction (MTJ) diameter. While the average temperature in the FL can be obtained based on an average current density and an averaged potential drop across the tunnel barrier, a fully 3D model is required to evaluate the large local temperature variations. The temperature variations in the FL further increase, when realistic device structures with an MTJ coated by a passivation layer and contacted to broad electrodes are considered. Lowering the FL thermal conductivity increases the temperature variations even further. However, the variations are strongly reduced with the decrease of the MTJ


Introduction
The ongoing miniaturization of semiconductor components has pushed the chip technology to its limits due to rapidly increasing leakage currents. Moreover, in the commonly employed von Neumann architecture, where the processing unit is separated from the memory, significant energy losses due to the data transfer back and forth, the so called von Neumann bottleneck, exist. The magnetoresistive random access memory (MRAM) is a promising emerging candidate to overcome these issues. MRAM is complementary metal-oxide semiconductor (CMOS) compatible [1][2][3] and has zero stand-by power consumption as it is intrinsically nonvolatile. It can also be integrated into logic [4], offers a wide temperature operation range [5], and has very good endurance.
A magnetic tunnel junction (MTJ), the basic building block of MRAM cells, consists of three layers: a pinned ferromagnetic layer, an oxide barrier, and a free ferromagnetic layer (FL). In spin-transfer torque MRAM (STT-MRAM) [6], relatively high current densities through the structure are required to switch the magnetization of the FL. This results in an increased temperature of the STT-MRAM cell, which mediates the switching of the FL magnetization [7,8]. On the other hand, the increased FL temperature caused by self-heating can result in an information loss, as it compromises the thermal stability [8]. To preserve the data, the temperature must be rapidly relaxed after writing. In [9], a heating asymmetry in the MTJ was observed for reversed current direction. This asymmetry was further numerically studied in [10], showing a non-linear increase of the saturation temperature with increasing power.
In this work we investigate the inhomogeneity of the temperature in the FL during switching. This inhomogeneity can be expected due to the fast magnetization dynamics of the STT-MRAM cell, which results in a significant inhomogeneity of the current density across the FL plane caused by a non-uniform magnetization [11][12][13][14]. Moreover, during switching, the current densities change significantly due to the rapid dynamics of the magnetization direction. Therefore, to model the temperature in STT-MRAM, the heat transport must be coupled to current, magnetization, and spin dynamics in a fully three-dimensional (3D) model. To the authors' knowledge, full 3D simulations coupling temperature, spin, and magnetization dynamics of STT-MRAM during switching have not been performed. The previous studies utilize either a simple one-dimensional (1D) model [15] for the heat or employ power (current) averaging over the barrier [7]. This manuscript is a follow-up article of a previously published conference contribution [16], which extends the previously published work to more realistic structures with broadened contacts and an oxide passivation layer, and investigates the effects of lower FL thermal conductivity and longer STT-MRAM contacts.
In Section 2, the general model used to simulate temperature, magnetization and spin dynamics is introduced together with the simulated structures and used parameters. Section 3 presents and discusses results from switching simulations under a constant voltage applied on the STT-MRAM cell. Temperature profiles of the FL for different structures are shown and analyzed. Section 4 summarizes the findings of this work.

Temperature modelling
To describe the dynamics of the temperature T in the structure at time t and position r, the heat flow equation is used.
c v , ρ, and κ stand for the heat capacity, mass density, and heat conductivity of the material, respectively. q is the heat source term. In MTJs, two main heat sources can be identified. The first heat source is attributed to the Joule heating in the ferromagnetic layers and in the metal contacts. It can be expressed as: q r (r, t) = j 2 (r,t)ρ E , where ρ E is the material resistivity and j is the current density.
The second heat source is associated with electrons tunneling through the insulating barrier [7,8]. When a potential difference is applied across the MTJ, an electron starts tunneling from the low potential side (source) and arrives at the high potential side (receiver) as a hot electron. Its energy is above the receiver Fermi energy and must be dissipated through various scattering processes. Similarly, the hole left behind the electron at the source must be filled by an electron from a higher energy level, the excess energy of which is therefore released.
Setting the x-axis to coincide with the main axis along the structure, the hot electron/hole heat source is described by [10].
where j x and ΔU stand for the x-component of the current density and the potential drop at the position (y, z) along the FL plane. The x-coordinate of the position of the free/pinned-barrier interface is denoted by x F/P and λ is a characteristic length at which hot electrons/holes lose their energy. The heat production imbalance between the receiver and the source side is accounted for by the asymmetry coefficient α(ΔU) [10]. In this work α is set to zerono asymmetry between the receiver and source sides is consideredin order to fully investigate the details of inhomogeneous temperature development due to the non-uniform current density. The system relevant parameters are listed in Tables 1 and 2.

Magnetization dynamics modelling
To describe the magnetization dynamics in the micromagnetic model the Landau-Lifshitz-Gilbert (LLG) equation is used.
Here m is the normalized magnetization, γμ 0 is the scaled gyromagnetic ratio, α is the Gilbert damping, and M S stands for the saturation magnetization. The effective field H eff consists of several components, with the demagnetization fieldH demag , the anisotropy field H aniso and the exchange field H ex contributing the most. In the following, H demag is determined by an optimized hybrid FEM-BEM approach [18], and H aniso and H ex are calculated with the parameters listed in Table 2. T S stands for the STT exerted by the spin current on the magnetization. To determineT S , the spin accumulation in the structure is computed [11,13]. The current densities and potentials within the structure, further used in the heat simulations, are determined by solving the coupled spin, charge, and magnetization dynamics in an MTJ [12]. The two main approaches to include temperature in the LLG equation introduce scaling of M S with temperature [19] or a random thermal field contribution to H eff [20]. The former seems not to be suited for STT-MRAM modelling as the switching remains deterministic, therefore, it is not fully usable for statistical switching simulations of STT-MRAM as it does not represent the correct dynamics of the system. The latter requires a careful temperature rescaling to account for the mesh size [20]. This method mainly reduces the incubation phase and introduces thermal fluctuations to the system, but does not have significant effects during the main switching phase. Therefore, the coupling of temperature to the LLG was omitted as the focus of the work is the FL temperature profile which is not significantly affected by the random thermal field inclusion.
The described approach was implemented in a fully 3D finite element method solver based on an open-source library MFEM [21]. For the time integration, an implicit Euler method was used.

Simulated structures
For the simulations, three different structures were considered. In Fig. 1, a simplified mesh of a cylindrical five-layer STT-MRAM structure is shown. The three-layer MTJ stack consisting of CoFeB(1 nm)/MgO(1 nm)/CoFeB(1.2 nm) is connected to non-magnetic metal (NM) contacts Taken from [15]. B. Taken from [17].
Parameter Value (30 nm). Both contact ends are kept at a constant temperature. No heat transfer is assumed across the shell of the cylinder. The structure diameter is 40 nm. In a real device, the MTJ stack is connected to wider contacts and is surrounded by an oxide passivation layer [15,22]. Fig. 2 shows such a structure with broadened contacts and SiO 2 serving as the passivation layer (yellow). The contacts consist of three different regions: 20 nm tall cylinder with a diameter of 200 nm, 10 nm tall conical tapering, and a 20 nm tall cylinder with a diameter of 40 nm. The MTJ diameter is again 40 nm (same layer color coding as in Fig. 1). The hollow space between the broadened contacts is filled with SiO 2 . There is no heat transfer assumed across the shell of the outer cylinder and both ends are kept at constant temperature (top and bottom of the structure in Fig. 2). The third structure is identical to one in Fig. 2, but the MTJ and the narrow contact ends have a smaller diameter of 20 nm.

Results
When a voltage is applied to the ends of the contacts, an electric current starts to flow through the structure, exerting an STT on the FL magnetization. For a high enough voltage, the FL magnetization is eventually flipped. In the simulations, 2 V and − 2 V between the contacts were applied to achieve a reliable parallel to anti-parallel (P-AP) and anti-parallel to parallel (AP-P) switching, respectively. We first consider the structure shown in Fig. 1. The averages of the magnetization components of the FL for P-AP (AP-P) at switching are displayed in Fig. 3 in the left (right) panel. In both plots, three different switching regions can be identified: Initial oscillations of the magnetization in the y-z plane, fast x magnetization component (m x ) change, and final magnetization oscillation in the y-z plane. While the fast m x changes are comparable for both P-AP and AP-P switching, the initial and final oscillatory phases vary significantly. In the P-AP switching, the initial (final) oscillations are much longer (shorter) lasting than for the AP-P switching. In both simulations, the FL magnetization is tilted by 5 • from the easy axis to eliminate the slow incubation phase. Fig. 4 shows the FL temperature increase development for P-AP (left panel) and AP-P (right panel). The contact end temperature is set to 300 K. In the following, we refer to the temperature increase of the FL above the end contact temperature as the temperature of the FL, but indicate it with δT to prevent any confusion. The average temperature δT avg (solid green), the minimum temperature δT min (dot-dashed gray), and the maximum temperature increase δT max (dashed black) of the FL obtained from the fully 3D model are shown. The average temperature increase δT avg− 1D calculated with an average current density and an average potential drop across the barrier is also shown (dotted orange) and coincides with δT avg from the 3D simulations. The fast initial heating phase (~150 ps), after which the FL temperature saturates for the given magnetization orientation, is followed by a slow temperature change due to dynamic magnetization behavior during switching. In the beginning, δT min and δT max nearly coincide. This indicates a homogeneous FL temperature profile. Later, δT min and δT max vary significantly and coincide again towards the end of the switching. Fig. 5 shows the temperature profile (left) next to m x (right) of the FL at 1 ns for the AP-P switching. The low (high) temperature region clearly correlates with the positive (negative) value region of m x . The heating in the STT-MRAM during switching is dominated by the tunneling electrons/holes (ca. 93-96 % of the total heat produced).
Another correlation is found when comparing Fig. 3 to Fig. 4. The maximum temperature difference at the FL,ΔT max ≡ δT max − δT min , is significantly increased when oscillations in the y-z plane decrease. This is further visualized in Fig. 6, where the ratio of ΔT max to δT avg is shown. In both switching scenarios, P-AP and AP-P, this ratio exceeds 30% in the middle switching region. During the switching, δT avg follows the increase of the average magnetization m x with approximately a 20-60 ps delay.
The more realistic structures described in Section 2.3 with broadened contacts and the SiO 2 passivation layer are now considered. For AP-P switching, − 2 V is again applied between the contacts and a switching comparable to Fig. 3 (right) is observed (not shown). The FL temperature in a cut through the center of the structure with the 40 nm MTJ at 1.75 ns is shown in Fig. 7. The warmest area is confined around the MTJ stack. Fig. 8 shows the temperature development of the FL for the MTJ structure with a diameter of 40 nm and 20 nm, respectively. For the MTJ with the diameter of 40 nm (left plot), the general temperature behavior is similar to that observed for the simpler structure (Fig. 4). The main difference can be seen in the last switching phase, where δT max and δT min remain different, in contrast to the results for the simpler structure. This is due to the presence of the passivation layer. The edges of the FL are cooled and have lower temperatures than the FL center. For the MTJ with the diameter of 20 nm (right plot) the situation results are dramatically different. BothδT max ,δT min , and ΔT avg nearly coincide. In this case, the MTJ is so narrow that even the cooled FL edges have almost identical temperature to the FL center. This fact is clearly reflected in Fig. 9, where the ratio of ΔT max to δT avg is shown for both structures. For the broad MTJ (left plot), ΔT max exceeds 45% ofδT avg . When the narrow MTJ is considered, the ratio ΔT max to δT avg stays below 5% for all  switching phases. Due the increased heat flux from the FL for both structures with the broadened contacts and the passivation layer, the initial saturation is reached faster (~100-120 ps) than for the simple scenario. Now we consider lower thermal conductivity in the FL -10 % of the original onefor the structure in Fig. 2. The rest of the parameters is kept unchanged. When the thermal conductivity is reduced, the temperature inhomogeneity further increases as show in the left plot in Fig. 10. The ratio of ΔT max to δT avg exceeds 50 % for the fast switching phase (not shown). As the last scenario we investigate the FL temperature dependence on the contact length. The left plot in Fig. 10 shows δT avg of the FL at 4 ns (end of switching) for different lengths of the narrow part of the contacts of the structure in Fig. 2. δT avg increases linearly with the increasing length of the contact. The linear fit crosses l = 0 nm at 32.4 K due to the remaining wide and tapered contact parts. ΔT max did not show any clear trend and remained within ±1 K of the Fig. 3. Averages of the normalized FL magnetization in x-, y-and z-direction during switching for the structure in Fig. 1. The left (right) panel shows switching from the initial parallel to anti-parallel (antiparallel to parallel) magnetic orientation. The parallel to anti-parallel switching shows a longer initial oscillating behaviour, whereas in the parallel to antiparallel switching, strong oscillating behavior is visible at the final switching phase. The initial magnetization was tilted by 5 • in the z-direction from the x-direction to eliminate the slow incubation phase. Fig. 4. Temperature in the FL during switching for the structure in Fig. 1. On the left (right) switching from the initial parallel to anti-parallel (anti-parallel to parallel) magnetic orientation is shown. An average temperature increase δT avg (solid green), maximum temperature increase (dashed black), minimum temperature increase (dotdashed grey) and δT avg− 1D (dotted orange, coincides with green) -calculated using averages of current densities and potential dropare shown. All temperatures are displayed with respect to the ambient temperature (300 K). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)    original value and therefore is not shown.

Conclusion
We have employed a fully 3D finite element model, coupling heat, charge, and spin transport to magnetization dynamics in order to study the FL temperature development in STT-MRAM during switching. At the beginning of the P-AP switching, the FL temperature reaches higher values than those for the AP-P switching, which will accelerate the initial switching phase. We also confirmed that a strong temperature inhomogeneity is developed across the FL, originating from the inhomogeneous current density due to non-uniform magnetization in a cylindrical STT-MRAM with a diameter of 40 nm. The ratio of the maximum temperature difference to the average temperature increase at the FL exceeds 30% for the simple five-layer structure in the rapid switching phase. While the maximum and minimum FL temperature cannot be reproduced using an average potential and current density across the FL, the average FL temperature coincides with the results of the fully 3D model.
When broadened contacts and a passivation layer are considered, the temperature variation in the FL with a diameter of 40 nm further increases due to some additional cooling at the FL edges caused by the passivation layer. The edge cooling will be further increased when a passivation material with a higher thermal conductivity, such as Si 3 N 4 , is used. The FL thermal conductivity reduction increases the temperature inhomogeneity in the FL even further, whereas a prolongation of the contacts only results in a linear increase of the absolute FL temperature. In the MTJ with a diameter of 20 nm, however, the temperature variation in the FL is strongly reduced. This suggests that a simpler 1D current model may be used with sufficient accuracy to calculate the FL temperature for MTJs with smaller diameters.
As the heat flux from the FL increases, when broad contacts and a passivation layer are used, the introduction of an additional insolation oxide layer between the MTJ and the contact may be desirable to trap the heat produced in the FL to facilitate the switching. However, this has to be done with caution as it will slow down the cooling of the FL, which may compromise the retention due to a possible data loss after writing.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Fig. 9. Ratio of the maximum temperature difference ΔT max of the FL to the average FL temperature increase δT avg for structures with broadened contacts and SiO 2 passivation layer. The structure with MTJ diameter of 40 nm (20 nm) is shown on the left (right). Due to the passivation, ΔT max is increased and exceeds 45 % of the average FL temperature increase for the wide MTJ, whereas it stays below 5% for the narrower MTJ. In addition to his expertise in condensed-matter physics, Tomáš has previous experience working in the fields of numerical acoustics and robotics. Tomáš joined the Christian Doppler Laboratory on Nonvolatile Magnetoresistive Memory and Logic in July 2020, where has been working towards his doctoral degree. His main focus is on numerical simulations of heat production, transport and thermal effects in the magnetoresistive memories.