Thermodynamic properties and switching dynamics of perpendicular shape anisotropy MRAM

The power consumption of modern random access memory (RAM) has been a motivation for the development of low-power non-volatile magnetic RAM (MRAM). Based on a CoFeB/MgO magnetic tunnel junction, MRAM must satisfy high thermal stability and a low writing current while being scaled down to a sub-20 nm size to compete with the densities of current RAM technology. A recent development has been to exploit perpendicular shape anisotropy along the easy axis by creating tower structures, with the free layers’ thickness (along the easy axis) being larger than its width. Here we use an atomistic model to explore the temperature dependent properties of thin cylindrical MRAM towers of 5 nm diameter while scaling down the free layer from 48 to 8 nm thick. We find thermal fluctuations are a significant driving force for the switching mechanism at operational temperatures by analysing the switching field distribution from hysteresis data. We find that a reduction of the free layer thickness below 18 nm rapidly loses shape anisotropy, and consequently stability, even at 0 K. Additionally, there is a change in the switching mechanism as the free layer is reduced to 8 nm. Coherent rotation is observed for the 8 nm free layer, while all taller towers demonstrate incoherent rotation via a propagated domain wall.


Introduction
In recent years, magnetic random access memory (MRAM) has been extensively studied and is considered one of the most promising emerging systems to replace current memory technologies.MRAM stores the binary digits as the magnetic states of a switchable layer (free layer) that can either align with a fixed layer (reference layer) or anti-align; one state represents a '0' and the other a '1'.Since the logic states are encoded by magnetic alignment, once the free layer has been saved into a state, no more power is drawn and thus this technology is non-volatile.Modern dynamic and static random access memory (RAM) (DRAM and SRAM) use electrical charge in capacitors to store the states and thus are volatile, leaving FLASH memory as the most significant non-volatile alternative [1].These volatile RAM devices draw a significant amount of power in HPC systems, and this is a large part of the appeal to introduce MRAM.This is coupled with MRAM's fast read/write speed, which makes potential performance comparable to SRAM.These features are often predicted to make MRAM the universal memory of the future, replacing both DRAM and SRAM on the market.However, for this to be a reality, MRAM needs to be scaled down to a sub-20 nm dimension to compete with the density of DRAM while maintaining a compromise of sufficiently high thermal stability for data retention and a low writing current [2][3][4][5].
MRAM consists of a CoFeB/MgO magnetic tunnel junction (MTJ), which contains a thin CoFeB reference layer with fixed magnetisation and a much thicker CoFeB free layer that can be switched, separated by a thin non-magnetic MgO layer.This choice of materials for the MTJ has been shown to have a high tunnel magnetoresistance ratio and a low damping constant to contribute towards a lower switching current.Additionally, the hybridisation of the Fe and O orbitals at the interface between the layers induces a perpendicular magnetic anisotropy (PMA) for additional stability [6][7][8].Initial studies of MRAM involved in-plane magnetisation and disc structures, whereby the width was significantly larger than the thickness of the free layer.However, it was shown that perpendicular magnetisation leads to an increased thermal stability factor, since the magnetisation aligns with the PMA direction.The thermal stability factor ∆ must be greater than 60 to meet the industry standard of a 10 year retention rate, while additionally being important for efficient magnetisation (i.e.minimal read/write error rate) [3,[9][10][11].
Various MRAM designs have been proposed to satisfy these criteria across past and current research.Examples include introducing dual MgO barriers and/or heavy metal capping layers for increased perpendicular anisotropy [12].However, dual MgO structures are not compatible with spin orbit torque (SOT) switching, which may be a beneficial technique for certain applications.It is therefore worth exploring additional design options [13].Another possible method to increase the anisotropy of the free layer is to create structures that introduce shape anisotropy in the same direction as the PMA.The disc-shaped structures mentioned previously had a shape anisotropy directed in-plane, and thus competed with the PMA, reducing the total anisotropy in the free layer.This is why those structures do not maintain a sufficient thermal stability at competing dimensions with DRAM.One solution is to create tower structures whereby the free layer thickness is larger than the tower's width.The resulting shape anisotropy now aligns with the PMA in the perpendicular direction, adding further stability to the system.These tower structures, known as perpendicular shape anisotropy MRAM (PSA-MRAM), are expected to scale to a sub-20 nm dimension while maintaining a sufficiently high thermal stability, without increasing the writing current prohibitively high [14][15][16][17].
In addition to the MTJ structure, we add two antiferromagnetically coupled CoPt layers below the reference layer to model a synthetic anti-ferromagnet (SAF).This is introduced for two reasons.Firstly, the CoPt layers provide a more stable reference layer, which ensures the reference layer does not switch throughout this study.Secondly, the CoPt layers can be engineered to minimise the stray field that may introduce a significant bias to the free layer [18,19].The details for these layers are discussed in the methodology section.
The magnetisation dynamics of nanoscale towers are very sensitive to precise dimensions and shape, and form the basis of this study.We aim to explore the effects of PSA on 5 nm diameter cylindrical tower structures as the free layer thickness is reduced to a sub-20 nm dimension.Taking a systematic approach, we also explore how temperature impacts the stability and switching mechanism as a function of free layer thickness.Since the inclusion of thermal fluctuations and the shape and interfacial anisotropy are essential to our goals, we use an atomistic spin dynamics approach implemented in the vampire software package.It is expected on this scale that finite size and edge effects will also be hugely influential on the material properties, and a micro magnetic discretisation will miss this essential physics [20][21][22].
The towers in this study consist of 7 nm of CoPt, with 5 nm magnetised downwards and 2 nm magnetised upwards.This is followed by a 1 nm CoFeB reference layer, which is pinned by the CoPt layers.A 1 nm MgO layer then separates the reference layer and the free layer, and then a CoFeB free layer has a changeable thickness.The thickness of the free layer is changed in steps of 10 nm, from 8 nm up to 48 nm.To capture the effect of the PMA, there is an enhanced CoFeB monolayer either side of the MgO.Shape anisotropy is captured via the demagnetisation tensor approach.The details of the whole tower is found in the methodology section.
We will start with susceptibility data obtained using the Monte Carlo method to explore how the dimensions of the device free layer changes the magnetisation with temperature.It is interesting to include these results, since they are in the absence of an external field.We will then focus on how the temperature and free layer dimension affect the stability and therefore the validity of these devices.To do this, we analyse hysteresis curves produced using the Landau-Lifshitz-Gilbert (LLG) equation.Starting with a free layer of 48 nm and decreasing in steps of 10 nm, we can see the trends that occur as we move towards the industry target of a sub-20 nm dimension free layer.This enables us to better understand what properties are dominating and emerging and how the shape anisotropy and edge effects begin to affect stability.Additionally, by doing these processes systematically at multiple temperatures, we can see the effects of thermal fluctuations on these properties.Finally, we look at the switching mechanism for different free layer dimensions.The hope is that these factors together will give a more complete understanding of how the thermal fluctuations, edge effects, and shape anisotropy, which have not been the focus of many previous studies, can be hugely important and even begin to dominate the behaviour of these nanoscale devices.The implications can then be considered in future design, development and creation of MRAM.This increases the chances of achieving an appropriately high thermal stability, with a low enough writing current, at cell volumes small enough to compete with the densities of DRAM and SRAM.This study switches the free layer with an external magnetic field, which will have a simplified domain wall propagation when compared to STT or SOT switching that may be used in commercial MRAM.This is due to the more advanced dynamics of a moving polarised current which can produce variation of rotational modes at the edges of the device compared to the centre [23].This study therefore explores the basic properties of PSA-MRAM, particularly how thermal fluctuations and shape anisotropy impact stability, to prepare for future work that introduces the more complicated physics of STT switching methods.

Methodology
The simulations have been carried out using the atomistic spin model from the vampire software package, which is based on the Heisenberg Hamiltonian [24,25].Thus, the spin Hamiltonian takes the Heisenberg form of exchange, where the terms on the RHS from left to right are the isotropic exchange energy, the magnetic anisotropy, the applied magnetic field from the zeeman interaction, and the demagnetising term [4].The two constants, k u and µ s , are the uniaxial anisotropy constant per atom and the atomic spin moment respectively.J ij is the exchange interaction between nearest neighbour spins and S i and S j are the normalised unit spin vectors on local sites i and j respectively [25,26].B ext is the applied external field vector and B demag is the long-range magnetostatic interaction [22,27].
We will start by modeling some equilibrium properties, and for this we use the Monte Carlo Metropolis algorithm.This will allow us to efficiently compute the temperature-dependent susceptibility of the system, compared to the LLG equation (below), which uses a small time step and has a slow convergence rate, and is therefore not as computationally timeefficient for this task [25].The first step of the Monte Carlo algorithm is a trial move, where a random spin i is selected and its spin direction S i is randomly changed to a new trial position S ′ i .The change in energy between to two states ) is then accepted with the probability by comparison with a uniform random number in the range 0-1.This process is repeated for N trial moves, where N is the total number of spins in the system, to form one Monte Carlo step [28].The trial move should be selected to yield as close to a 50% acceptance rate as possible, since a rate too high or too low will increase the number of Monte Carlo steps required.At low temperature, for example, the new spin direction is not permissible if its direction has jumped too far, since this would represent a large thermal fluctuation.In contrast, at increased temperature the allowable range of thermal fluctuation is much larger.To maintain the desired acceptance, we use an adaptive algorithm [29] with a tunable width, where each spin is moved according to where Γ is a randomly generated number and σ g is the acceptance-probability dependent width of the cone.
While the Monte Carlo Metropolis algorithm is best for equilibrium properties, it cannot model magnetisation dynamics since it does not provide an evolution in time.Instead, we use the atomistic LLG equation to model magnetisation dynamics, which is of the form where γ is the gyro-magnetic ratio, λ is the microscopic damping moment and B eff = − 1 µs ∂H ∂Si is the effective field acting on the local spin [4,30].Thus, the LLG equation describes how an atomistic spin moment interacts with an effective field.However, in its present form, the LLG equation is only applicable at zero K.The inclusion of thermal fluctuations is vital to this study to explore the validity of MRAM at operational temperatures, so we must add thermal fluctuations to this model.This is achieved by adding a thermal field to the effective field.We also add a demagnetising field for the inclusion of shape anisotropy since this is also not captured thus far.The equation for the effective field now becomes Thermal effects are included using Langevin Dynamics as developed by Brown, whereby thermal fluctuations on each atomic site are represented with a Gaussian white noise term [31][32][33].The magnitude of the thermal fluctuations is governed by the width of the Gaussian distribution Γ(t), which increases with temperature, where the distribution in three dimensions has a mean of zero [34].For each time step, ∆t, the thermal field for each spin i is given by where k B is the Boltzmann constant, λ is the Gilbert damping parameter, γ is the absolute value of the gyro-magnetic ratio, ∆t is the integration time step and µ s is the atomic magnetic moment.The effective field acting on each spin is thus temperature dependent, which naturally captures the loss of exchange and anisotropy at increased temperatures [25,35].The demagnetising field is a long range interaction and is therefore calculated separately using the inter-intra dipole tensor approach in vampire which discretizes the system into macrocells.These are treated as dipoles, calculating nearby macrocells to atomistic resolution, and follows from the work of Bowden [36].This magnetostatic energy in equation (1) takes the form 1

∑
p m p mc • B p demag where B p demag is the magnetostatic field in macrocell p and is given by where r is the distance between the dipole of p and q, r is the unit vector in the direction − → pq, V p mc is the volume of the macrocell p, and m p mc is the magnetic moment of macrocell p.The first term therefore describes the dipolar field acting on p from all other macrocells, while the second term represents the self demagnetisation of macrocell p on itself [22,25,36].Collating terms in equation (7), the pairwise interactions between the discretized cells can be put into the matrix M pq , for improved performance, The magnetostatic field in each macrocell p is then given by Since this is a long range interaction, all layers will contribute a small dipole field to all other layers in the system.For shapes with one dimension significantly larger than the others, the inclusion of the demagnetising field provides shape anisotropy.Since this study focuses on cylindrical shape PSA-MRAM, the demagnetisation field is a particularly essential addition to the effective field.The LLG equation is then integrated using the Heun predictor-corrector integration algorithm which is well suited to the stochastic nature with larger time-steps for computational efficiency [4,22,25].Now that we have discussed the mathematical models used in vampire, we now outline the atomistic details of the PSA-MRAM tower structures.An atomistic model such as vampire is an ideal tool for these towers.Its applicability sits between the short time and length scales of quantum mechanical abinitio models (DFT), and the large time and length scales best suited to continuum micromagnetic models.The PSA-MRAM tower stacks in this paper contain between 25 000 and 100 000 atoms, which is too many to be viable via computationally expensive DFT models.However, micromagnetic models that discretise the system into macrocells cannot capture edge effects, material boundary variations or temperature, which all have a significant effect on the magnetisation dynamics and observed behaviours.vampire satisfies these requirements, since shape anisotropy is included via the dipole interaction, thermal effects via Langevin Dynamics, and edge and material boundary effects are naturally captured due to its atomistic nature.
In an atomistic model, every atom is assumed to have a magnetic moment localised to the atoms lattice site, known as the spin moment.The schematic of our studied tower structure is shown in figure 1.It is a body centred cubic (bcc) structure with a lattice constant of a = 0.2866 nm.The towers are built with the easy axis along the positive z direction, with all layers' initial magnetisation pointing up or down along this axis.Each layer then has its own set of parameters shown in table 1 and discussed below.These include the exchange constant, uniaxial anisotropy and magnetic moment which we find from ab-initio and experimental data.The free layer is the only material to change volume in this study, starting at 48 nm thick, and being decreased to only 8 nm thick in A schematic of the PSA-MRAM tower structure being studied, consisting of two anti-ferromagnetically coupled CoPt layers followed by the MTJ with two high anisotropy CoFeB mono-layers on either side of the MgO separating layer.The free layer is studied with a thickness of 8 nm, 18 nm, 28 nm, 38 nm and 48 nm.We also note the easy axis for our structure is along the positive z-direction.
our smallest tower.The diameter of the cylinders are always 5 nm, which is particularly slim compared to other research on PSA-MRAM to date.This means there will always be a small amount of shape anisotropy in the free layer, since its thickness is always larger than its width.The tower is divided into seven sections corresponding to the labels in the figure.This includes two CoPt layers representing the SAF, the reference layer and free layer sandwiching a MgO layer, then two enhanced mono-layers.The details of these layer is now discussed.
The reference layer, free layer, monolayers and the MgO make up a standard MTJ.We consider an idealised case, where the uniaxial anisotropy is provided by just the monolayers either side of the MgO barrier.In vampire, this is achieved by simply creating a mono-layer of CoFeB sandwiching the MgO with enhanced values for the anisotropy and the exchange constants [37,38].The rest of the free layer and the reference layer are modelled as having zero anisotropy.This is a valid approach, as the uniaxial anisotropy and enhanced exchange arises from the hybridisation of the Fe and O orbitals at the interface, thus is a localised effect [27,39].The anisotropy of bulk CoFeB is negligibly small, so approximating the rest of the CoFeB to zero is valid [27,40].We also do not model the Co, Fe and B individually; instead we use an average magnetic moment for all the atoms in the CoFeB layers.The MgO layer separates the reference layer from the free layer but is not included explicitly in the calculations and interactions since it is non-magnetic.The reference layer is required to have greater stability than the free layer, since it must not fluctuate significantly from its magnetic state.At operational temperature, a thin 1 nm reference layer of CoFeB is not sufficiently stable to satisfy this requirement.We therefore include two additional CoPt layers below the standard reference layer to create a SAF structure.The CoPt in contact with the CoFeB layer is magnetically aligned, along the positive z-direction for this study.These two layers together now form the top half of the SAF.The purpose of this design is to provide stability to the CoFeB.The reference layer may then be considered as these two materials together, but throughout this paper the CoFeB alone is referred to as the reference layer.Since we do not care about the behaviours of these layers, so long as the reference layer is stable, this does not affect any results/discussion.In some studies, the CoFeB could be replaced entirely by CoPt as it has a greater stability.However, CoFeB has a higher spin polarisation than CoPt, so is a better choice material for STT or SOT MRAM.Future research using more complicated models including STT are currently in the works, so retaining this CoFeB structure as the reference layer is useful for future work.These two layers together, while stable, would produce a significant stray field in the free layer.This causes an undesirable bias in the switching field for the free layer.This is the reason for the bottom half of the SAF structure, where additional CoPt is anti-ferromagnetically coupled to the first CoPt.This coupling reduces the stray field via the demagnetising field outlined above in equation (9).The thickness of this anti-ferromagnetically coupled CoPt was adjusted until only a small field is present in the free layer.Since all layers contribute to the dipole field, it would be challenging to perfectly cancel the stray field in the free layer.Our dimensions shown in figure 1 minimises the switching field bias in the free layer, but will still be slightly present in the results.This is reasonable, since annihilating any stray field in experimental studies would be equally problematic, and a small bias is expected.To summarise these layers, the SAF structure has provided enough stability to the reference layer to avoid any switching, with minimised stray field on the free layer.The reference layer and the CoPt layers do not switch at any point in this study, due to this stability.
The free layer starts aligned with the reference layer in all the simulations in this paper.The material values are shown in table 1 [1,8,38,41].
The values of exchange and anisotropy for the CoFeB bulk and interface layers were taken from previous studies, where the non-magnetic Boron is neglected.This treats the system as crystalline CoFe, with a ratio of Co-Fe of 1:3 [27].The justification comes from studies of the annealing process, where the majority of the Boride atoms are diffused out of the free layer, normally towards a Ta layer, so we just model the resulting CoFe [35].We use a slightly higher magnetic moment, µ s = 2.5, corresponding to bulk CoFe, while some studies use a smaller value of 1.6 J/T [42].
The exchange constant for the CoPt layers had to be calculated using the mean-field expression where T c is the Curie temperature, z is the number of nearest neighbours and ϵ is the correction factor.The Curie Temperature for CoPt is found to be 723 K [43], z = 8 and ϵ = 0.766 for bcc structures [25] so using equation ( 10) the exchange constant for CoPt is J ij = 4.88 × 10 −21 J/link.The magneto-crystalline anisotropy was also calculated using where K u is the macroscopic anisotropy constant, a is the lattice constant and n atoms is the number of atoms per unit cell.For CoPt at zero K, K u = 2.83 × 10 6 Jm 3 [44] and n atoms = 2 so using equation ( 11) the magneto-crystalline anisotropy for the CoPt layers is k u = 3.33 × 10 −23 J/atom.This model does not capture any impurities or defects from the annealing process or the boundaries between materials.It also treats the MRAM tower as a perfect cylinder shape, with no intermixing or overlapping of atoms at the boundaries between layers.The effect of the MgO is included via the enhanced mono-layer as discussed, rather than any attempt to replicate the physical hybridisation of orbitals using DFT models.

Susceptibility
Here we explore the temperature-dependent properties of the MRAM towers with a focus on the effects of size and temperature on the switching of the free layer.Firstly, we explore the static equilibrium properties.In figure 2, we plot the susceptibility as a function of temperature for the free layer to compare five different free layer thicknesses.The mean magnetic susceptibility for the free layer is given by where we plot χ m , the longitudinal susceptibility, in figures 2 and 3.The tallest tower's peak in figure 2 was an order of magnitude larger than the smallest tower.We therefore decided to divide each susceptibility curve by the free layer thickness to Susceptibility against temperature for the free layer of each tower thickness.As the thickness decreases, the curves are smoother and the peak is lower.This suggests more correlated spins as the free layer volume decreases, which may predict a trend towards coherent rotation.
Figure 3. Susceptibility against temperature for the whole 48 nm FL (as in figure 2), compared to local 10 nm thick sections of the free layer.The localised curves are all the same order of magnitude and width (with the 8 nm section slightly reduced) when compared to the whole FL.This demonstrates that locally, the spins are correlated, so the noise in the whole FL is a consequence of variation in the spin direction or non-uniform excited modes.
normalise them.As a result, we can plot them all on the same graph for comparison, giving the units mT −1 nm −1 .Since the free layer magnetisation is saturated, the susceptibility slowly increases with temperature before reaching the Curie temperature.This is in contrast to experimental expectation, where the magnetisation decreases with temperature.Effectively, it is not possible to over saturate until the Curie temperature in this model.The focus of the susceptibility plots is therefore the peak at the Curie temperature.We see that as the free layer thickness increases, the peak gets larger, noisier, wider and occurs at a lower temperature.Together, these results suggest more correlated spins with reduced free layer volume.This is explained as a consequence of mixing the longitudinal and transverse susceptibility in the taller towers, since the data output is an averaged susceptibility of the entire free layer as shown in equation (12).The longitudinal susceptibility should describe fluctuations in exchange length, and will be present in all of the towers regardless of dimensions.This is expected to increase as temperature is increased, with a large peak at the material's Curie temperature.This is what we aimed to plot in figure 2, however, there is also transverse susceptibility which stems from fluctuations of the magnetisation perpendicular to the easy axis.The fluctuations in the magnetisation of each atomic site is due to thermal fluctuations, as outlined in the methodology.With greater thermal fluctuations, the spin direction may fluctuate by larger amounts.However, neighbouring atomic sites still experience a relatively strong dependence upon each other due to the exchange constant.This partially restricts the freedom of transverse movement, leading to predominantly longitudinal fluctuations in any localised enough section of the free layer.Even as the temperature increases towards the Curie temperature, only very small amounts of transverse fluctuation may occur over these short atomic ranges.However, the exchange interaction is not long range, so after several atomic spaces, the correlation in transverse direction begins to reduce.In other words, transverse fluctuations are minimal from one neighbouring spin to the next due to the exchange constant, but over a larger volume begins to exaggerate.We therefore expect the smallest tower is predominantly longitudinal susceptibility, since any fluctuations in magnetisation direction are negligible and equation (12) would predominantly give us longitudinal susceptibility.But an undesired consequence of increasing the thickness in this model is that it inadvertently adds transverse susceptibility too.Fluctuations in magnetisation direction at the bottom of the free layer (closest to the MgO layer) versus the top of the free layer may be much less correlated, presenting much larger peaks.The peak shifting to the left with increased free layer thickness is therefore a product of coherency.The peak of a small free layer is caused only by thermal fluctuations, while taller towers additionally suffer from a loss of coherency.This causes a lower Curie temperature, along with a taller peak.
To prove that this is the case, we took the tallest tower from figure 2 and divided its 48 nm thick free layer into 10 nm sections (with one 8 nm section to make 48 nm) and then ran the same program again.In figure 3, we plot the susceptibility against the temperature for these local sections and again we normalised by the free layer thickness to keep the units the same as figure 2. The region labeled 0 − 10 nm corresponds to the bottom of the free layer that is in contact with the MgO, while the 40 − 48 nm region corresponds to the top of the free layer.The interfacial anisotropy at the 0 − 10 nm region has no significant effect here, since variation in transverse susceptibility is dominated by the exchange constant.We confirm that locally, the susceptibility is behaving the same as the 8 nm case in figure 2, with the peak reaching the same magnitude.We conclude that the origin of the larger susceptibility curves in free layers of increased thickness is a consequence of increased transverse fluctuations, while longitudinal susceptibility dominates the smallest free layer.This data then suggests that the taller towers have non-uniform magnetisation modes throughout the stack, since we have shown increased transverse susceptibility in taller free layers.We can use this to predict the reversal mechanism for taller towers to be incoherent, whereas the smallest tower would be expected to be closer to coherent rotation as it is closer to uniform magnetisation.

Hysteresis
Next, we investigate the dynamic hysteretic properties of the free layer.This demonstrates how varying the thickness of the free layer and the temperature of the system affects the magnetic properties such as the coercivity.Identifying precisely how the shape (free layer volume) and the temperature impact the switching mechanism is crucial to understanding the limitations and validity of these devices [45].The coercivity is known to vary with the sweep rate in LLG models, so care must be taken choosing our parameters.For all hysteresis loops in this paper, the system is first equilibrated for 10 000 time steps (∆t = 1fs) in the absence of a field to reach a ground state.For the hysteresis data, the damping is changed to the critical damping, α = 1, to allow rapid relaxation towards a ground state without affecting the magnetic properties being explored.A magnetic field of 1T is then applied along the easy axis, and is swept to −1T and back to 1T in steps of 0.01T.At each field increment, 100 000 time-steps of 1fs occur, which given the large damping is sufficient for relaxation of the spins with the field.Increasing the times-steps or decreasing the field sweep rate by a factor of 10 results in the simulations becoming too computationally expensive.While our hysteresis loops will still be slightly larger than expected in reality, this setup is an appropriate compromise of accuracy and computation time and is similar to previous studies using vampire [13,22,25].
We first demonstrate the impact and importance of including thermal effects, which is a benefit of using an atomistic model, in figure 4. Thermal fluctuations are modelled with a Gaussian white noise term with a random number generator responsible for its distribution, which reflects the stochastic nature of thermal fluctuations.This means that each simulation is repeatable, since the sequence of numbers is the same, but this produces square loops representing a precise coercivity.In reality, thermal fluctuations should add some uncertainty to the precise coercivity.To capture the random nature of fluctuations, we run 40 independent loops with a different integration seed, then average those loops as shown in figure 4. We see that at 300 K, each individual run produces a square loop, representing a precise value of the coercivity for that random seed.However, given a different sequence of random numbers, this precise value can shift slightly.Across 40 independent loops, a normal distribution of switching fields appears.The range of switching fields is referred to as the switching field distribution (SFD) and represents an uncertainty in the precise switching field of the free layer.Averaging all of the individual Figure 4. Forty individual hysteresis loops with different integration seeds and then the averaged loop at 300 K. We see the randomness of thermal fluctuations changes the precise point the system switches, so the averaged loop has a curved shape representing a switching field distribution.The Gaussian overlays are calculated in section C and confirms that 40 loops are enough to produce a normal distribution for averaging purposes.
square loops produces a curved loop, which is also shown in figure 4. The curvature of the averaged hysteresis loop represents the SFD, corresponding to the small expected uncertainty in the systems switching field.In section 3.3 of this paper, the SFD for this curved averaged loop is extracted by means of an error equation.Details of this process are found in that section.The value for the SFD is then the difference between the leftmost and rightmost curve in the positive region of figure 4. We can then use this SFD and the extracted coercivity to plot a Gaussian distribution on the rightmost curve in figure 4.This demonstrates that the majority of the individual square loops lie within the main peak of the distribution curve.We therefore conclude that 40 loops is sufficient for this averaging process, which is in agreement with previous studies of MRAM using vampire [22].The use of the error equation is discussed later, but it is beneficial to the reader to have outlined the SFD at this stage of the paper.Since all other parameters included in the model were kept constant for the 40 independent loops, such that the only difference was the integrator seed for the Gaussian, we can conclude that this is entirely a thermal SFD, in agreement with [22].All subsequent hysteresis loops in this paper are averaged in this way.
We now compare the different free layer thicknesses at three different temperatures in figure 5 to systematically view the trends in magnetic properties for these systems.In the 0 K case, there are no thermal fluctuations, so no SFD is present.As a result, averaging over 40 independent loops was not necessary, all 40 loops would do the same thing.It is worth noting that there is still slight curvature in the loops at 0 K.This is due to the ballistic effect from a finite step size, and does not represent a switching distribution.It is still useful to include the 0 K case however, since it shows the relationship between different free layer thicknesses and additionally demonstrates the theoretical best for these tower structures.Hysteresis loops for the free layer of five towers with differing free layer thickness at (a) 0 K, (b) 150 K and (c) 300 K. We observe a reduction in the coercivity for all towers as temperature increases, but we also see that at all temperatures the 8 nm and 18 nm towers have reduced coercivity compared to the saturated 28 nm, 38 nm and 48 nm towers, due to a loss of shape anisotropy.This may result in a lack of thermal stability for MRAM volumes small enough to compete with the cell densities of modern RAM.
We see the tallest three towers have reached an asymptotic maximum value for the coercivity, with the 18 nm free layer slightly reduced followed by a comparatively large reduction for the 8 nm free layer.This is an interesting result, as the coercivity would initially be expected to increase linearly with volume.However, the coercivity is also related to the anisotropy of the system, and the shape anisotropy does not increase linearly.The asymptotic limit observed for the taller free layers is therefore a consequence of reaching maximum shape anisotropy for cylinders of this diameter.Below around 30 nm, the rapid reduction of stability as the free layer thickness is reduced could pose significant limitations on the usage for these 5 nm diameter devices.We also see a very small bias between the right hand and left hand branch, which is due to the long range demagnetisation field stemming from the CoPt SAF and reference layer [13].As discussed in the methodology section, the thickness of the bottom CoPt layer was adjusted to minimise the stray field present in the free layer.It is very difficult to completely remove any stray field, and not possible in the real systems we are trying to emulate.For this reason, a small stray field in the free layer is acceptable, which results in the very small bias.The free layer therefore has a very small preference for being aligned with with the reference layer.In other words, a very slightly stronger field is required to switch from aligned to anti-aligned compared to the reverse.Another possible contribution towards a bias could originate from distortions of the magnetisation at the MgO/CoFeB interfacial layer in anti-parallel alignment (e.g. a flower state) [46].However, in figure 5(a), we see no evidence of a significant flower state in our studied structures, so we conclude that the magnetisation can be considered essentially uniform across the x-y plane for our small lateral sizes.In reality, there is likely to be very small distortions, but due to the exchange constant being very large for so few atoms, the distortions are insignificant and this contribution to the bias is ignored for the rest of this study.It is worth noting that we would expect this effect to become more dominant in towers of greater diameter.Comparing the 150 K and 300 K cases in figures 5(b) and (c) respectively, we note that the trends in the 0 K case are still present.The tallest towers are still at an asymptotic limit and the greatest reduction in coercivity occurs as the free layer is decreased to the 8 nm thickness.The bias between the right hand and the left hand coercivity for all heights at all temperatures is 0.09 T to two decimal places.This is worth noticing, because the bias in the free layer stems from the stray field emanating from all other layers.Since this does not change to two decimal places at increased temperature, the reference layer/ SAF layers must not deviate significantly from their vertical magnetisation.Noticing this result is therefore a recognition of the stability of the reference layer, which was the purpose of the CoPt layers.Beyond two decimal places, there is a very slight decrease in the bias as temperature increases, of around 0.001 T between 0 K and 300 K.This is because, at finite T, the lower layers will always deviate a small amount from the idealised vertical magnetisation due to thermal fluctuations.As a result, a slightly smaller stray field is present in the free layer compared to 0 K.
We can now also compare the difference in temperature on the system and immediately notice the reduction of coercivity as temperature increases for all free layer thicknesses.The thermal fluctuations directly reduce the saturation Figure 6.Hysteresis curves with increasing temperatures for the 8 nm thick free layer.The curved shape of the loop becomes more exaggerated with temperature representing an increased switching field distribution.This suggests that thermal effects are driving the switching mechanism.magnetisation and anisotropy, which leads to the reduction of the coercivity for all of the towers.The reduced shape anisotropy of the smallest tower, possibly coupled with increased edge effects due to the loss of exchange at edges contributes to a significant reduction for the smallest tower.The curvature also increases as the temperature increases representing a larger thermal SFD due to thermally activated transitions over the energy barrier.The width of the SFD is largest for the smallest tower at both 150 K and 300 K due to the lower shape anisotropy.This leads to a larger relative impact of thermal fluctuations on the switching process.These results suggest that at operational temperatures, the 8 nm free layer has a greatly reduced stability compared to the taller towers that are still converged on a more stable asymptotic value.
To improve our discussion of thermal effects on the reduced volume, we plot hysteresis for the 8 nm free layer again but with more temperature iterations to better see the trend in figure 6.We can see that as the temperature increases, the shape of the hysteresis loop is converging towards the paramagnetic regime, whereby the system loses any remnant magnetisation and cannot be used for storage.However, at 300 K we still observe a distinct loop and have not reached the paramagnetic behaviour.It is also much easier to see the bias here, with the right-hand branch dropping below 0.5 T from 50 K, while the left-hand branch drops below −0.5 T at 150 K.This is an unavoidable consequence of providing the necessary stability to the reference layer and shows the importance of the inclusion of dipole-dipole interactions.The effect of thermal fluctuations on the hysteresis loops shows the importance of including temperature in an atomistic approach, since we have shown the SFD is not constant and in fact the reversal is partially driven by these effects.

Systematic trends in the coercivity
The interesting features of the hysteresis curves are both the coercivity and the SFD, since these tell the field range required to switch the free layer.A low coercivity and high SFD present challenges for stability and recording quality of MRAM devices [47].Here, the coercivity is the value of the magnetic field as the hysteresis curve crosses the x-axis.However, the width of the SFD demonstrates how thermal fluctuations can begin to drive the switching process at reduced field strength, before this coercivity is reached.Extracting the values for the coercivity and the SFD is therefore useful for analysing the hysteresis curves.We extract these values by fitting each individual switching branch of the hysteresis curves to an error function which takes the form where, because the curves are not centered on x = 0, x is given by x 0 is then the shift from x = 0, so provides the coercivity, while σ is the standard deviation.The error function is the integral of the normalised Gaussian function, thus provides the cumulative distribution of the Gaussian.From figure 4, the purpose of repeating for 40 independent loops was to ensure a normal distribution of the coercivity.The width of the error function, σ is therefore the width of the Gaussian distribution of coercivities, which represents the SFD.We demonstrate the application of equation ( 13) in figure 7.Here we fit the equation to the right-hand branch of the 8 nm thick free layer at 300 K.In this example, fit we extract a coercivity (x 0 ) of 0.28 T and a SFD (σ) of 0.07 T. The right side Gaussian that was added to the plot in figure 4 is given by 1 2π e − (x−x 0 ) 2 σ 2 , with this value for x 0 and σ.This fitting process is then repeated for every branch of all free layer hysteresis loops.We can then compare the coercivity and the SFD as a function of free layer thickness and temperature.
Repeating this process for all the temperatures and free layer thicknesses for the right hand branch of the hysteresis curves, we produce figure 8, which plots the coercivity against temperature.The crosses represent our data points and the error bars represent the width of the SFD curvature.The solid lines in this figure are the Sharrock fit, described in section 3.4.For this section, we focus on our output data (the crosses) and the error bars.Like in figure 5, we see the reduction of the coercivity as the temperature increases, but we further see that the trend is virtually identical regardless of the free layer thickness.This is an interesting observation, since it suggests the temperature dependence of the coercivity is not volume dependent.As a result, this is not something that may be engineered away with volume changes in the MRAM stack.It is challenging to observe this feature in figure 5, supporting our decision to express the coercivity in this way.The right-hand branch of the hysteresis curve for the 8 nm thick free layer at 300 K with its fit curve from equation ( 13).Here we observe the fit is reasonable for crossing the axis at the right point, and we extract a coercivity of 0.28 T and a standard deviation of 0.070 T to 2 s.f.We conclude that this method of extraction is valid to analyse the branches in figure 6.
The standard deviation extracted is shown using error bars and represents the SFD at each temperature.We notice two trends that are not so obvious in figure 5.The error bars that drop below the extracted coercivity represent a reduced field, below the coercivity, that could begin to switch the free layer due to the thermal fluctuations driving the switching process.There is therefore an uncertainty in the precise field required to switch the free layer, since it could theoretically happen at any field point within the error bar.The first observation, therefore, is that the three tallest free layers in figure 8 have overlapping SFDs for all temperatures 50 K and above, meaning they could all be switched at the same field.This is in contrast to the exact coercivity (the precise field strength in a square hysteresis loop) which demonstrates the importance of including thermal effects.The overlap for these tallest towers adds context to the asymptotic maximum observed in figure 5, since we now see that they do in fact overlap when including a SFD.It also demonstrates how significantly the coercivity drops as the free layer is further reduced to 18 nm and 8 nm thick.We see that the 18 nm thick free layer has a coercivity that is lower than the taller three towers even accounting for the SFD, confirming that the free layers are beginning to lose stability at these dimensions.The extreme drop in coercivity for the 8 nm free layer is then highlighted when we notice that the lowest field for the 18 nm free layer at 300 K is still higher than the 0 K best-case coercivity for the 8 nm.The second observation is that the SFD grows larger as the temperature increases, which is expected as larger thermal fluctuations cause larger fluctuations in the saturation magnetisation and begin to drive the switching process.For all five free layer thicknesses, the increase in the width of the SFD with each 50 K temperature increase is between 0.008 T and 0.016 T, with the largest increases at the sub-150 K data points.We see the coercivity decreases and the SFD increases with temperature.This also demonstrates the significant loss of shape anisotropy for the 8 nm free layer.
The large reduction of the coercivity for the 8 nm tower is due to the shape anisotropy.This is included in the model as a consequence of the demagnetisation factors, N xx , N yy , N zz , shown in table 2.Here we see the diagonal components of the matrix shown in equation ( 8) (since off-diagonal components are near zero) where these factors are determined from the ratios of the length of the axis [48].It is clear that the taller towers have significantly higher demagnetisation factors in the x and y direction, with comparatively small factors in the z direction.This produces a strong preference for the spins to align along the easy axis and results in the increased coercivity even at 300 K seen in figure 7. Conversely, the smallest tower free layer is only 8 nm tall, which is comparable to its diameter of 5 nm.This results in minimal shape anisotropy supporting the interfacial anisotropy from the mono-layer, which is comparable to previous studies of disc like structures MRAM [22].We also see that the reduction of the shape anisotropy is not linear, since the tallest three towers have a demag factor in x and y that is more than 4x greater than the demag factor in z.The 18 nm tower shows a small decline, but is still over 3.5x greater in the x and y direction compared to its z dimension, while the 8 nm free layer is less than 2x greater.A free layer of 5 nm thickness would be the same as the diameter and would have no significant shape anisotropy, so we see that the 8,nm tower is approaching this rapidly.(A 5 nm free layer would have a small shape anisotropy due to the towers being perfectly cylindrical rather than spherical, but it would be insignificant).
From these demagnetisation factors in table 2, the shape anisotropy energy for each height tower may be calculated.Since the structures in this paper are cylindrical, these ).Since the crystal anisotropy for our CoFeB free layers is negligibly small, these values for the shape anisotropy energy are also the anistropy energy K u to a good approximation.These values for K u for each height tower will be used in the following section.

A comparison with the Sharrock equation
It is interesting to compare the trend in our coercivity data to that predicted by the Sharrock law, which is an expression for the thermal dependence of the coercivity that comes from the Arrhenius-Neel law.The Sharrock equation takes the form where T is the temperature, k B is the Boltzman constant, K is the anisotropy constant, V is the volume of the switchable layer, τ is the relaxation time and f 0 is the frequency.Slightly less obvious, H a is described as the anisotropy field ) , corresponding to the precise field that reduces the energy barrier to zero, while H c is the actual coercivity [49,50].This is an important distinction, since at finite temperature H c will be lower than H a , as thermal fluctuations partially contribute to transitions over the energy barrier.At 0 K however H c = H a , which we exploit to compare our temperature dependence for the coercivity with that predicted of Sharrock.
The temperature dependence of the coercivity from our model was demonstrated in figure 8 for the different thickness free layers.The temperature dependence of the coercivity in the Sharrock model in equation ( 15) is of the form with A = H a and B = − √ kB KV ln(f 0 τ ).We then perform a fit of equation (16) to our data, obtaining optimised values for A and B. The obtained value for A is then the prediction for H a from the Sharrock model, which should be comparable to our value of H c (T = 0).To demonstrate, all of the fits are shown in figure 8. Our data points align reasonably well with the A comparison of the extracted coercivities from our 0 K hysteresis loops to the Sharrock approximation.We see the trend is similar for the two lines but the extracted coercivities are lower for all free layer thicknesses.The 8 nm free layer is closest to the Sharrock approximation, likely due to coherent rotation.Sharrock fit at finite temperature.All free layer thicknesses show small deviations from Sharrock fit, but the trend is the same.Small fluctuations are expected, as the Sharrock model is simplistic, assuming coherent rotation.We the believe that the uptick in the Hc values at low temperature is due to to Hc approaching Ha, leading to a small energy barrier (at the switching point) and breakdown of the high energy barrier approximation used by Sharrock.
The fitted H c (T = 0) as a function of free layer thickness is shown in figure 9 along with the value from the atomistic calculations.We see good agreement in the trend between our results and the Sharrock model fits.For all free layer thicknesses, the Sharrock model predicts slightly higher value for H a than we find in our model.This discrepancy is likely due to the simplified coherent rotation assumed in the Sharrock model.This suggestion is supported, as the 8 nm free layer displays coherent behaviour and has the closest agreement with the Sharrock equation (∆H 8nm a = 0.05T).All the thicker free layers display incoherent switching behaviour and have larger discrepancy, with the 18 nm having slightly larger (∆H 18nm a = 0.07T) and the three thicker free layer having a greater discrepancy (∆H 28 nm a = ∆H 38 nm a = ∆H 48 nm a = 0.08T).From the the Sharrock equation fits, we can further obtain an estimate for the thermal stability factor ∆. The thermal stability factor is given by ∆ = KV kBT [51].From fitting equation (16) to our data, a value for B is obtained, where B = − √ kB KV ln(τ f 0 ), so we can rearrange this to find kB KV = ln(τ f0 B 2 .We can therefore re-write the thermal stability factor as Where B is obtained from our fit, T is the temperature, and an estimate must be found for the value of ln(τ f 0 ).We obtain an estimate for ln(τ f 0 ) from our fitting of the Sharrock equation, as follows.
Since we are estimating a value for ln(τ f 0 ) using the Sharrock model, which assumes coherent rotation, our best estimate will come from the fit to the 8 nm free layer of this study, as it is also coherent.Firstly, by differentiating equation (15) with respect to √ T we find where we have introduced a term G = dHc d √ T for easier notation.Replacing the macroscopic anisotropy K = HaMs 2 and rearranging equation ( 18) we arrive at the expression The first fraction is known from our fit of equation ( 16), as it is straightforward to show that G = AB and H a = A. The second fraction contains known constants in our model, M s = µsn a 3 = 1.97 × 10 6 A m −1 , V is the volume of the free layer and k B is the Boltzmann constant.From the fit to the coherent 8 nm free layer, which is the lowest fit in figure 8, we obtain an estimate using equation (19) of ln(τ f 0 ) = 8.2251.With this value, we can further estimate the thermal stability for the 8 nm free layer at 300 K using equation (17), which yields ∆ = 24.7567.Using equation (17) for the other free layers, we obtain the thermal stability as a function of free layer thickness as seen in figure 10.
The prediction for the thermal stability stemming from the Sharrock equation as outlined above can be compared to the analytical model.Analytically, the thermal stability is given by ∆ = KuV kBT , where an estimate for K u for each free layer thickness is found in table 2. We therefore add the analytical value for the thermal stability to figure 10 to contrast with our Sharrock prediction.The analytical equation predicts a linear trend, with the thermal stability growing as the volume of the free layer grows, yielding a very large thermal stability (∆ ≈ 390) for the 48 nm free layer compared to the Sharrock prediction.However, this major discrepancy can be accounted for, as the analytical approach assumes coherent rotation, where the whole free layer is involved.In fact, reversal is via nucleation and propagation which means that only the top of the tower has to be switched.Comparing the thermal stability of the 48 nm free layer with the Sharrock estimate (∆ ≈ 85) gives an estimate of the height associated with the switching volume of 85 390 × 48 ≈ 10 nm.This is consistent with our findings that the switching becomes increasingly non-uniform for heights greater than about 8 nm.We can therefore add a third comparison to figure 10, using ∆ = KuV kBT , but with a constant volume, representing the volume associated with switching (10 nm thick).This shows strong agreement with the Sharrock model, predicting the same trend.
The Sharrock prediction and the analytical prediction with constant volume in figure 10 suggest that the 8 nm free layer would not satisfy the requirement for a thermal stability of ∆ = 60.The lack of shape anisotropy for such a small tower Figure 10.A of the calculations of the stability factor the different free layer thicknesses at 30 K. The analytical prediction, KuV kBT , show a linear trend as volume increases.By fixing the volume to be the volume affiliated with switching (10 nm thick), we see agreement between the Sharrock model and the analytical model.We observe the expected decrease in stability as the thickness is decreased towards 8 nm, with the stability plateauing beyond 28 nm.This suggests nothing is gained in terms of stability for increased towers and represents a limit for towers of 5 nm diameter.at these dimensions results in a loss of data retention due to thermal fluctuations.The 18 nm free layer is found to have a thermal stability ∆ = 60.6, just above the industry required threshold, while the taller towers are all similar at ∆ ≈ 80.Long term stability requires zero-field simulations of the energy barrier, which are beyond the scope of the current work.For this reason, we only focus on the trends predicted in figure 10 in this paper.
It can be seen that the thermal stability plateaus at around a 30 nm free layer thickness, suggesting that any growth above that in the z direction will not gain further stability and represents the maximum for towers of 5 nm diameter.We see that the reduction in the thermal stability factor is in agreement with our previous data, suggesting a significant and non-linear reduction of shape anisotropy as the tower free layer thickness is reduced.The loss of shape anisotropy begins as the free layer is reduced to thickness below 30 nm, with a greater reduction as the thickness approaches the width [5,22].This is expected, since a free layer of 8 nm is only marginally above the 5 nm width.Previous studies have already shown that below this threshold, when the width becomes greater than the thickness, the shape anisotropy acts perpendicular to the interfacial anisotropy direction and MRAM ceases to be viable [40,52].
Given an value of τ it is possible to make an estimate of the pre-exponential factor f 0 .This is complicated by the fact that the magnetisation is calculated in the atomistic model using a sweep rate process.The effective time for a swept field process is given by by Chantrell et al [53] as where R is the sweep rate, H a is the anisotropy field extracted from our fit, M s is the saturation magnetisation, and V is the volume of the free layer.With this we can use equation (19) to evaluate an estimate of f 0 .The calculated values for f 0 for the free layer thicknesses given in table 3, along with the corresponding values for V and H a .The values of f 0 are physically reasonable and also increase with tower height which is consistent with the expectation of the increase in f 0 with anisotropy predicted by the Brown relaxation time [31].We note that the estimate of the energy barrier from the Sharrock law can only give guidance.However, the plateau in the stability factor is most likely physically realistic on the basis that reversal begins with the nucleation of a small reversed region of magnetisation which is apparently at most weakly dependent on the tower height for values greater than 30 nm.A calculation of the exact energy barrier and its size and temperature dependence for the atomistic model is indicated, however this can only be done using the constrained Monte Carlo method [54] which involves very CPU intensive calculations.This is beyond the scope of the current paper.Several results in this paper thus far have suggested that the 8 nm free layer may be switching coherently, while the taller towers are incoherent.To demonstrate whether this is the case, our final results in figures 11 and 12 are snapshots of the tower structures at 0 K and 300 K respectively, produced using POV-Ray.We output the spin files that create the images at every complete field step, so each snapshot is separated by 0.0001 T, or 0.1 ns, and we are watching the free layer switch from its aligned state to anti-aligned.Firstly, we provide snapshots at 0 K in figure 11, where the switching mechanism is entirely a consequence of the magnetic properties of CoFeB, particularly the PMA and PSA.In the absence of thermal fluctuations, we observe a clear propagated domain wall for the larger free layers while the 8 nm free layer appears coherent.This is in agreement with the susceptibility results in figures 2 and 3, which suggest non-uniform magnetisation modes in free layers greater than 10 nm.Additionally, the 8 nm free layer displayed a closer agreement to the sharrock model in figure 9, where the model is used to treat coherent rotation and does not account for the incoherent rotation seen in the taller towers.We then compare this to the 300 K case in figure 12 where we have previously shown in figures 4-8 that thermal fluctuations lead to a SFD.While figure 12 shows the coherent rotation for the 8 nm free layer and the incoherent rotation for the taller free layers that we have already seen at 0 K, we now additionally see thermal noise as individual spins fluctuate from a unified direction.In figure 12(d) for instance, we see the 48 nm thick free layer switches in more steps than in figure 11(d), even though they are both in steps of 0.0001 T. This is because the thermal fluctuations begin the switching process at a lower field as there is a thermal SFD.Snapshots of the switching mechanism at 300 K for different free layer thickness.We observe coherent (8 nm) and incoherent (18 nm-48 nm) switching as in the 0 K case of figure 11, except now we have additional fluctuations of spin orientations from random thermal noise.

Conclusions
We have investigated the effect of free layer thickness and temperature on the magnetic reversal process for 5 nm wide cylindrical PSA-MRAM tower structures.We have included thermal fluctuations and magnetostatic effects using an atomistic spin model in vampire and focused primarily on the hysteresis and switching dynamics.From our susceptibility data, we found uniform magnetisation modes in thickness up to around 10 nm, where further extension causes non-uniform magnetisation modes as you traverse the stack, which may be responsible for non-coherent rotation.We then moved on to hysteresis plots, and found that thermal fluctuations are responsible for driving the switching process with significant SFDs, proving these effects cannot be ignored in future studies of these nanoscale systems.By extracting the coercivities and SFDs to observe the trend with temperature, we conclude that the reduction of shape anisotropy, as the thickness approaches the width is responsible for a significant loss of thermal stability in these devices.Given the engineering challenge of physically creating towers of 5 nm diameter, it is significant to show how rapidly these towers lose stability below 20 nm thickness.We compare the coercivity results of our model to the Sharrock model with surprisingly good agreement in the trends at 0 K.This trend suggests that the coercivity reaches an asymptotic limit when the thickness increases to around 30 nm, with any subsequent growth along the easy axis lacking any significant gain in stability, suggesting a fundamental limit on stability for towers as thin as 5 nm diameter.This is further supported by the trend found for the thermal stability factor, which levels off for a free layer around 30 nm thick.Our approximation for the thermal stability factor is not suspected to be an accurate value, so further study is required to determine whether these structures satisfy the industry requirement of ∆ > 60, but we do expect that the trend we predicted is reflective, such that the towers do plateau at a maximum stability.This is supported by a lack of increased shape anisotropy beyond 30 nm.Finally, we demonstrate the switching mechanism at 0 K is an incoherent propagated domain wall motion for all towers taller than 10 nm, but is coherent when further decreased below this threshold, which is in line with estimates of the domain wall length in previous studies.Further study should focus on STT switching of the free layer and improved estimates for the thermal stability factor based on an accurate prediction of the energy barrier.

Figure 1 .
Figure 1.A schematic of the PSA-MRAM tower structure being studied, consisting of two anti-ferromagnetically coupled CoPt layers followed by the MTJ with two high anisotropy CoFeB mono-layers on either side of the MgO separating layer.The free layer is studied with a thickness of 8 nm, 18 nm, 28 nm, 38 nm and 48 nm.We also note the easy axis for our structure is along the positive z-direction.

Figure 2 .
Figure2.Susceptibility against temperature for the free layer of each tower thickness.As the thickness decreases, the curves are smoother and the peak is lower.This suggests more correlated spins as the free layer volume decreases, which may predict a trend towards coherent rotation.

Figure 5 .
Figure 5. Hysteresis loops for the free layer of five towers with differing free layer thickness at (a) 0 K, (b) 150 K and (c) 300 K. We observe a reduction in the coercivity for all towers as temperature increases, but we also see that at all temperatures the 8 nm and 18 nm towers have reduced coercivity compared to the saturated 28 nm, 38 nm and 48 nm towers, due to a loss of shape anisotropy.This may result in a lack of thermal stability for MRAM volumes small enough to compete with the cell densities of modern RAM.

Figure 7 .
Figure7.The right-hand branch of the hysteresis curve for the 8 nm thick free layer at 300 K with its fit curve from equation (13).Here we observe the fit is reasonable for crossing the axis at the right point, and we extract a coercivity of 0.28 T and a standard deviation of 0.070 T to 2 s.f.We conclude that this method of extraction is valid to analyse the branches in figure6.

Figure 8 .
Figure 8.The coercivity as a function of temperature for each height tower.The data points are the extracted coercivity with the switching field distribution for each point shown as error bars.The solid lines are the fitted Sharrock equation, discussed in section 3.4.We see the coercivity decreases and the SFD increases with temperature.This also demonstrates the significant loss of shape anisotropy for the 8 nm free layer.

Figure 9 .
Figure 9.A comparison of the extracted coercivities from our 0 K hysteresis loops to the Sharrock approximation.We see the trend is similar for the two lines but the extracted coercivities are lower for all free layer thicknesses.The 8 nm free layer is closest to the Sharrock approximation, likely due to coherent rotation.

Figure 11 .
Figure11.Snapshots of the switching mechanism at 0 K for different free layer thickness.We observe the smallest tower switching coherently, while all other towers demonstrate an incoherent mechanism, with an observable propagated domain wall motion. τ

Figure 12 .
Figure12.Snapshots of the switching mechanism at 300 K for different free layer thickness.We observe coherent (8 nm) and incoherent (18 nm-48 nm) switching as in the 0 K case of figure11, except now we have additional fluctuations of spin orientations from random thermal noise.

Table 1 .
The exchange, anisotropy, and magnetic moment constants for the bulk CoFeB reference and free layer, the enhanced CoFeB monolayer and the CoPt layers (MgO is not directly involved in the magnetic properties and does not have these constants).

Table 2 .
Demagnetisation factors (Nxx, Nyy, Nzz) for the different free layer dimensions, showing the loss of a dominant direction with reduction in volume.This loss of shape anisotropy significantly reduces the coercivity and stability factor for the 8 nm tower and represents a limitation in scalability.The final column is the calculated anisotropy energy Ku for each height tower.

Table 3 .
(20)ntial components for the Sharrock model.The volume for each height is included for completness, as it is used in several equations.Ha is the value extracted by the Sharrock model in figure8, and f 0 is found using equations(19)and(20).