Phase-field simulation of magnetic double-hole nanoring and its application in random storage

ABSTRACT

As an ideal high-density storage unit, magnetic nanorings have become a research hotspot in recent years. We can both study the evolution of microscopic state of magnetization and acquire macroscopic magnetic properties by micromagnetic simulation, which has thus been widely used. However, traditional micromagnetism cannot simulate complex stress state. Due to the introduction of microelasticity theory, the phase field method for magnetic materials can be used to calculate the coupling effect of stress and magnetic field. However, the computing model usually needs to satisfy periodic boundary condition. In this paper, the phase field simulation combined with the finite element method is employed. By using user defined element, the evolution of magnetic domain structures of the double-hole nanorings has been studied. In different diameter of the holes and external magnetic field direction, we have found seven kinds of magnetic domain evolution mechanism. Among them, the twin-vortex evolution mechanism with high stability and low demagnetization interference characteristics of advantages, has good application prospect in magnetic random-access memory (MRAM) unit.

Graphical Abstract

KEYWORDS:

• Micromagnetic simulation
• phase field method
• finite element method
• magnetic random-access memory
• magnetic nanoring

Introduction

In recent years, magnetic random memory (MRAM) based on GMR (or TMR) effect has attracted great interest from researchers in various countries because of its nonvolatile, high density and high-speed reading. As an ideal high-density storage unit, magnetic nanorings have become the focus of research in recent years. In 2000, J G, Zheng Y, Prinz G A [1] proposed a random-access memory VMRAM design based on micromagnetic simulation analysis, which utilizes the macromagnetic resistance effect of vertical magnetic multi-layer, has good robustness by making the storage unit into a circular magnetic multi-layer stack, and estimates the limit storage surface density of 400Gbits/in². In 2001, Schneider M, Hoffmann H, Zweck J [2] prepared a sub-micron Pomo alloy disc using electron beam lithography and liftoff techniques, observed the magnetization flip mechanism and remaining magnetic configuration using the Lorenz transmission electron microscope, and revealed the possibility of adjusting the magnetization spin by introducing slight geometric asymmetry. In 2001, Rothman J, Kläui M, Lopez-Diaz L, et al [3]. studied the evolution of Co(100) magnetic ring magnetism using magnetic measurements and micromagnetic simulations, and discovered a new stable double domain state, the onion state. In 2002 Lai, Mei-Feng et al [4]. simulated the soft magnetic material Pomo alloy nanoring with micromagnetic, discovered a new stable twin-vortex state, and compared it with magnetic force microscope (MFM) images. In 2005, Moore T A, Hayward T J, Tse D H Y, et al [5]. used magneto-optical Kerr effect (MOKE) magnetic forces to determine the hysteresis return of the polycrystalline Permalloy ring with an outer diameter 2 μm and 5 μm, and used micromagnetic simulations to reveal the onion-to-vortex flipping mechanism, pointing out that the positive and reverse flipping of the vortex state is random, depending on historical configuration and thermal fluctuations. In 2005 M. Rahm et al [6]. found that by introducing hole defects, they could increase the number of hysteresis loops. They then investigated multiple holes and applied magnetic field control to make the magnetic domain state switch between multiple steady states, which they wanted to achieve multi-process storage. Furthermore, they designed a logical memory because their magnetic field flips do not need to be applied to saturation, and their flip fields are much smaller. In addition, they studied magnetic field flip fields and logical states. In 2006, Zhu F Q, Chern G W, Tchernyshyov O et al [7]. prepared asymmetric nanorings to regulate the vortex formation process by changing the direction of the external magnetic field. In 2006 Mani A S, Geerpuram D, Baskaran V S et al [8]. studied the effects of shape on vortex formation and vortex chirality, and found that the introduction of asymmetrically provided an accurate control of the flip mechanism. In 2007, Yang T, Hara M, Hirohata A et al [9]. prepared Co nanorings of horizontal size of deep sub-microns, and experiments revealed that the giant magneto-resistance effect of vertical ring planes can effectively describe a very small magnetic nanoring, and found that the transition between the nanoring onion state and the vortex state can be achieved not only by applying an external magnetic field but also by vertically through direct current. In 2012, Huang C H, Cheng N J, Wu F S et al [10]. used micromagnetic software to simulate asymmetric Pomo alloy rings, demonstrating that asymmetric rings can regulate not only the vortex chirality but also the vortex formation and vanish field.

With the improvement of computer storage and computation speed, micromagnetic simulation method has been developed rapidly and widely used. The main calculation methods can be divided into two broad categories, namely, the finite differential method(FDM) proposed by Hughes [11] in 1983, and the finite element method(FEM) proposed by Fredkin and Koehler [12] in 1987. Since the demagation field is along-range interaction field and its calculation are time-consuming, Bertram and Zhu [13] propose to introduce a fast Fourier transformation (FFT) in FDM calculations to calculate the demagation field. Magnetic scalcation potentials ϕ are usually introduced in finite element simulations to calculate static magnetic fields. Because the boundary of the magnetic field is an open Dirichilet boundary condition, the magnetic extraterritorial infinite area can be mapped to a finite area by angle-protected transformations [14], but these transformations typically apply only to the boundaries of a particular shape. Since the regression field decays rapidly in the air, the demagation magnetic field can also be approximated by taking the outer region large enough to approximate the de-magnetic field [15], but this will increase the amount of storage and calculation. Fredkin [16] proposed in 1990 to combine FEM with boundary integrals to calculate the demagation magnetic field, which greatly improved the computational efficiency and gradually developed the FEM-BEM (Boundary Element Method) hybrid calculation method.

Phase-field simulation methods were developed in the 1970s, initially to bypass the difficulty of tracking solid–liquid interfaces when simulating solidified tissue, and were later widely used in various tissue evolution problems. The phase field model describes different phases by introducing continuously changing field variables – order parameters in the entire calculation area. The order parameters change sharply at the interface, and there is no need to track the geometry of the interface. Therefore, the phase field is a diffuse interface model. The phase field theory is based on statistical physics. According to the Ginzburg-Laudau phase transition theory, it reflects the combined effects of diffusion, ordering potential and thermodynamic driving force through differential equations. As a powerful computational simulation method, phase field method is widely used in the simulation and prediction of material-based morphology and microstructive evolution [17]. In recent years, this method has been used to study the micro-evolution and macro-performance of ferromagnetic, ferroelectric and magnetic composite materials. In 2003, Koyama [18] used the phase field method to study the microstructure evolution of ferromagnetic shape memory alloys (MSMAs) Ni2MnGa under applied stress and magnetic field. Its free energy includes chemical-free energy, gradient energy, elastic energy, and magnetic energy. Among them, the magnetic energy only considers anisotropic properties. In 2012, Miehe [15] used node displacement, magnetic potential, and magnetization vector as field variables, established corresponding finite element equations based on elastic mechanics, magnetostatic theory and TDGL evolution equations, compiled 2D magnetic elements, and analyzed the regulation of the magnetic field and stress on the micro/nano structure of permalloy. When solving the demagnetization problem, Miehe takes a large enough external free space, and assume that the magnetic potential at the boundary decays to zero. Without considering the air element, Wang [19] in 2013 established a 3D model and analyzed the regulation of the nano-magnetic structure by strain, magnetic field and their coupling effects. Using this method can not only calculate the boundary of complex shapes, but also calculate the force-magnetic coupling response of materials under complex stress. However, this method is not yet mature enough and needs to be further developed and improved.

In this paper, the method of combining finite element and phase field [20] is used to derive the corresponding incremental variational principle from the equations of mechanics, magnetostatics and magnetization vector motion. The finite element equation is established, and the two-dimensional is compiled through the ABAQUS software subroutine development module. The ferromagnetic element can simulate the microstructure of any boundary under the action of magnetic field and complex stress field. Through the written user element, we study the evolution mechanism of the magnetic domain of the double-hole nanoring, and focus on the regulation mechanism of the geometric aperture and the direction of the external magnetic field to the double-hole ferromagnetic nanoring magnetic domain. Under the parameters of different geometric aperture and external magnetic field direction, we found seven kinds of magnetic domain evolution paths: twin-vortex mechanism, left and right vortex mechanism, ‘S’ evolution mechanism, ‘S’ + ‘C’ mechanism, ‘S’ and ‘twin-vortex mechanism’, ‘S’ + left vortex mechanism, left vortex mechanism. Among them, the twin-vortex mechanism has a good application prospect in random storage unit because of its advantages of high stability and low demagnetic interference.

Phase field model and finite element formulation

Basic equations and boundary conditions

Mechanical equation

Under the linear elasticity and small deformation hypothesizes, the basic equation of elasticity can be used to describe the behavior of magnetic materials. The mechanical equilibrium equation can be expressed as:

(1) $\mathrm{\nabla }\cdot \text{sigmab}+{\mathbf{F}}_b=\rho \ddot{\mathbf{u}}(V)$(1)

$\sigma$ represents stress and $V$represents elastomer. $\ddot{u}$ is the second derivative of $\mathbf{u}$ with respect of time $t$. Since the body force ${\mathbf{F}}_b$ is not considered in this problem, and the effect of inertial force $\rho \ddot{u}$ is very small, which can be ignored, formula (1) can be simplified as [15]

(2) $\mathrm{\nabla }\cdot \text{sigmab}=\mathbf{0}(V)$(2)

Boundary conditions can be represented as [15]

(3) $\sigma \cdot \mathbf{n}=\mathbf{t}(A)$(3)

$A$ represents the magnetic elastomer boundary, $\mathbf{t}$ represents the surface force on the boundary, and $\mathbf{n}$ represents the direction vector of the common line outside the boundary. Geometric equations based on the assumption of small deformation can be represented [15]

(4) $\epsilon =\frac{1}{2}\left(\mathbf{f}+{\mathbf{f}}^T\right)$(4)

where $\epsilon$ represents strain, $\mathbf{u}$ represents displacement, $\mathbf{f}\mathbf{=}\mathrm{\nabla }\mathbf{u}$. By (2)-(4) one can get the displacement variational equation of elasticity:

(5) ${\int }_V\sigma :\delta \epsilon dV={\int }_A\mathbf{t}\cdot \delta \mathbf{u}dA$(5)

Magnetostatics equation

The total magnetic field can be decoupled into two parts, the de-magnetic field and the external magnetic field, recorded as [20]

(6) $\mathbf{h}\mathbf{=}{\mathbf{h}}_d\mathbf{+}{\mathbf{h}}_e$(6)

$\mathbf{h}$ represents the total magnetic field, ${\mathbf{h}}_d$ represents the de-magnetic field, is determined by the geometry of the structure and the magnetization state, and ${\mathbf{h}}_e$ represents the external magnetic field. Degenerate Maxwell equations can get magnetostatic equations [19(19) $E_{exch}=A\left(m_{1,1}^2+m_{1,2}^2+m_{1,3}^2+m_{2,1}^2+m_{2,2}^2+m_{2,3}^2+m_{3,1}^2+m_{3,2}^2+m_{3,3}^2\right)$(19) ]

(7) $\mathrm{\nabla }\cdot \mathbf{b}=0$(7)

Among them $\mathbf{b}$ is magnetic induction strength, which can be expressed as the relationship between it and the magnetic field [20]

(8) $\mathbf{b}={\mu }_0\left({\mathbf{h}}_d+{\mathbf{h}}_e+m_s\mathbf{m}\right)$(8)

The saturation magnetization strength $m_s$ represents the magnitude of the magnetization vector $\mathbf{m}$, which is the material constant and represents the direction cosine of the magnetization vector. When the external magnetic field ${\mathbf{h}}_e$ is uniform,

(9) $\mathrm{\nabla }\cdot {\mathbf{h}}_d=-m_s\mathrm{\nabla }\cdot \mathbf{m}$(9)

The relationship between the magnum potential $\varphi$ and the demagnetic field can be expressed [20]

(10) ${\mathbf{h}}_d=-\mathrm{\nabla }\varphi$(10)

Boundary conditions at the magnetic material interface can be represented

(11) $\begin{array}{c}{\varphi }^{+}={\varphi }^{-}\\ \frac{\mathrm{\partial }{\varphi }^{+}}{\mathrm{\partial }n}-\frac{\mathrm{\partial }{\varphi }^{-}}{\mathrm{\partial }n}=-m_s\mathbf{m}\cdot \mathbf{n}\end{array}$(11)

where ${\varphi }^{+}$,${\varphi }^{-}$ represents the magnetic potential of the magnetic body and the outer boundary, respectively. The first formulation in (11) indicates that the magnetic potential is continuous at the magnetic body boundary, and the condition is automatically satisfied because it is a field function in finite element simulation. The second means that the flux flowing in and out of the boundary is equal, i.e., the static flux is 0. Magnetic potential needs to be satisfied in infinity

(12) $r\to \mathrm{\infty },\varphi \to 0$(12)

Because the de-magnetic field decays rapidly in the air, a large enough area can be taken outside the magnet to make the magnetic potential at the boundary of the region 0, thus replacing infinity boundary conditions (12) with an approximate boundary condition [15]

(13) $\varphi =0(B)$(13)

$B$is the boundary of the external air area. According to (9)-(11) and (13) the magnetostatic variational equation can be derived

(14) ${\int }_{\mathrm{\Omega }}{\mathbf{h}}_d\cdot \delta {\mathbf{h}}_{d}dV=-{\int }_V{m}_s\mathbf{m}\cdot \delta {\mathbf{h}}_{d}dV$(14)

where the air magnetizes the vector $\mathbf{m}=\mathbf{0}$.

Motion equations of magnetization vector

The motion equation of the magnetization vector usually adopts the LLG equation, but for the study of static magnetic properties, the TDGL equation can also be used [19]

(15) $\frac{\mathrm{\partial }\mathbf{M}(\mathbf{x},t)}{\mathrm{\partial }t}=-L\frac{\delta E}{\delta \mathbf{M}(\mathbf{x},t)}$(15)

where$L$is the motion coefficient, $E$is the density of free energy, and [15]

(16) $\mathbf{M}(\mathbf{x},t)=m_s\mathbf{m}(\mathbf{x},t)$(16)

The density of free energy in the magnetic body can be expressed as follows:

(17) $E=E_{anis}+E_{exch}+E_{elas}+E_{mag}+E_{cons}(V)$(17)

where $E_{anis}$ is the Anisotropy energy, $E_{exch}$ is the exchange energy, $E_{elas}$ is the elastic energy, $E_{mag}$ is the sum of the magnetostatic energy and the external magnetic energy, and $E_{cons}$ is the constraint term, which is called the penalty function in the finite element method. The specific forms of each item will be given below. For cubic anisotropic crystal, its magnetocrystalline anisotropy can be expressed as [19]

(18) $E_{anis}=K_1\left(m_1^2{m}_2^2+m_2^2{m}_3^2+m_3^2{m}_1^2\right)+K_2{m}_1^2{m}_2^2{m}_3^2$(18)

where $K_1$, $K_2$ are material constant, respectively. The exchange energy is related to the change gradient of the magnetization intensity and can be represented as follows [19]:

(19) $E_{exch}=A\left(m_{1,1}^2+m_{1,2}^2+m_{1,3}^2+m_{2,1}^2+m_{2,2}^2+m_{2,3}^2+m_{3,1}^2+m_{3,2}^2+m_{3,3}^2\right)$(19)

where $A$ is the exchange constant. Elastic energy expressions can be divided into two parts:

(20) $E_{elas}=E_{elas1}+E_{elas2}$(20)

$E_{elas1}$ represents energy items related only to elastic deformation, and $E_{elas2}$ energy items related to magnetic elastic coupling, as follows [19]:

(21) $E_{elas1}=\frac{c_{11}}{2}\left({\epsilon }_{11}^2+{\epsilon }_{22}^2+{\epsilon }_{33}^2\right)+c_{12}\left({\epsilon }_{11}{\epsilon }_{22}+{\epsilon }_{22}{\epsilon }_{33}+{\epsilon }_{11}{\epsilon }_{33}\right)+2c_{44}\left({\epsilon }_{12}^2+{\epsilon }_{23}^2+{\epsilon }_{13}^2\right)$(21)

(22) $\begin{array}{rl}& E_{elas2}=-\frac{3}{2}{\lambda }_{100}\left(c_{11}-c_{12}\right)\left({\epsilon }_{11}{m}_1^2+{\epsilon }_{22}{m}_2^2+{\epsilon }_{33}{m}_3^2\right)\\ & -6{\lambda }_{111}{c}_{44}\left({\epsilon }_{12}{m}_1{m}_2+{\epsilon }_{23}{m}_2{m}_3+{\epsilon }_{13}{m}_1{m}_3\right)\end{array}$(22)

where $c_{11}$, $c_{12}$, $c_{44}$ and are the anisotropic elastic constants of cubic crystals, ${\lambda }_{100}$ and ${\lambda }_{111}$ are the magnetostriction coefficients. For an ideal crystal of cubic crystal system, its elastic properties are completely described by a matrix of elastic constants (stiffness coefficient) $c_{ij}$, which contains three independent elastic constants $c_{11}$, $c_{12}$, and $c_{44}$. Among them, only the constant $c_{44}$ has a direct physical meaning as a measure of the resistance of the crystal to the shear force in the <100> cube plane along any direction lying in this plane. The matrix coefficients $c_{11}$ and $c_{12}$ have no such simple explanation. Only by using the corresponding Hooke’s law symbols can we understand their meaning [21]. For cubic materials, there are two independent magnetostriction constants ${\lambda }_{100}$ and ${\lambda }_{111}$. And ${\lambda }_{100}$, ${\lambda }_{111}$ are the saturation magnetostrictions along the <100> and <111> axes of the crystal [22]. And

(23) $\begin{array}{rl}& E_{mag}=-{\mu }_0{m}_s\left(m_1{h}_{d1}+m_2{h}_{d2}+m_3{h}_{d3}\right)-{\mu }_0{m}_s\left(m_1{h}_{e1}+m_2{h}_{e2}+m_3{h}_{e3}\right)\\ & -\frac{{\mu }_0}{2}\left(h_{d1}^2+h_{d2}^2+h_{d3}^2\right)\end{array}$(23)

where ${\mu }_0$ is the vacuum magnetic conductivity. Magnetization vectors should meet constraints $\left|\mathbf{m}\right|=1$, which are constrained by the introduction of energy terms $E_{cons}$ whose expressions can be written as [19]

(24) $E_{cons}=A_s(\left|\mathbf{m}\right|-1{)}^2=A_s{\left(\sqrt{m_1^2+m_2^2+m_3^2}-1\right)}^2$(24)

This term is equivalent to the penalty function term in the finite element method. The value of $A_s$ is not the larger the better, too large will easily cause singularity and cause the calculation to fail to converge. It should be determined according to the specific situation. Without considering the surface anisotropy, the magnetization vector usually satisfies the following equation at the boundary [20]:

(25) $\frac{\mathrm{\partial }\mathbf{m}}{\mathrm{\partial }n}=\mathbf{0}or\mathrm{\nabla }\mathbf{m}\cdot \mathbf{n}=\mathbf{0}$(25)

Regarding the variation of the free energy density $E$ to the magnetization vector $\mathbf{m}$, there is the following relationship [15]:

(26) $\frac{\delta E}{\delta \mathbf{m}}=\left\{\frac{\mathrm{\partial }E}{\mathrm{\partial }\mathbf{m}}-\mathrm{\nabla }\cdot \left(\frac{\mathrm{\partial }E}{\mathrm{\partial }\mathrm{\nabla }\mathbf{m}}\right)\right\}$(26)

According to (15), (25), (26), we can get

(27) ${\int }_V\left(m_s^2\frac{\mathrm{\partial }\mathbf{m}}{\mathrm{\partial }t}\cdot \delta \mathbf{m}+L\frac{\mathrm{\partial }E}{\mathrm{\partial }\mathbf{m}}\cdot \delta \mathbf{m}+L\frac{\mathrm{\partial }E}{\mathrm{\partial }\mathrm{\nabla }\mathbf{m}}:\delta \mathrm{\nabla }\mathbf{m}\right)dV=0$(27)

The basic equations and boundary conditions related to mechanics, magnetostatics, and magnetization vector motion are given above, which are the basis for deriving the finite element form.

Phase field finite element form

The weak form of equivalent integrals

Since there is a demagnetizing field in the air element, its free energy density can be expressed as

(28) $E=-\frac{{\mu }_0}{2}\left(h_{d1}^2+h_{d2}^2+h_{d3}^2\right)(\mathrm{\Omega }/V)$(28)

Then according to the free energy density (17) in the ferromagnetic body, the total-free energy within the model can be obtained [20]:

(29) $G={\int }_{\mathrm{\Omega }}EdV$(29)

External forces can be expressed as:

(30) $W={\int }_A\mathbf{u}\cdot \mathbf{t}dA$(30)

The dissipation potential of the magnetization vector motion is introduced [15]:

(31) $\mathrm{\Phi }\left(\dot{\mathbf{m}}\right)={\int }_V\frac{\eta }{2}{\dot{\mathbf{m}}}^{2}dV$(31)

where $\eta$ is the inverted number of the motion coefficient. In a time step$\left[t_n,t_{n+1}\right]$, set the interval to $\tau$, and the total potential energy of the system can be written as [15]

(32) $\prod =G-W-(G_n-W_n)+\tau \cdot \mathrm{\Phi }\left(\dot{\mathbf{m}}\right)$(32)

where$G_n$, $W_n$ represent the free energy of the moment $t_n$ and the external force work, which are constants, so that let the constant$C=G-W$. There are also the following relationships

(33) $\dot{\mathbf{m}}=\frac{\mathbf{m}-{\mathbf{m}}_n}{\tau }$(33)

${\mathbf{m}}_n$ is the cosine of the magnetization vector direction at time $t_n$. (32) can be rewritten as [20]

(34) $\begin{array}{rl}& \prod =G-W+{\int }_V\frac{\eta }{2\tau }{\left(\mathbf{m}-{\mathbf{m}}_n\right)}^{2}dV-C\\ & ={\int }_{\mathrm{\Omega }}EdV+{\int }_V\frac{\eta }{2\tau }{\left(\mathbf{m}-{\mathbf{m}}_n\right)}^{2}dV-{\int }_A\mathbf{u}\cdot \mathbf{t}dA-C\end{array}$(34)

From the potential energy functional form in (34), the differential equations (2)(2) $\mathrm{\nabla }\cdot \text{sigmab}=\mathbf{0}(V)$(2) , (9(9) $\mathrm{\nabla }\cdot {\mathbf{h}}_d=-m_s\mathrm{\nabla }\cdot \mathbf{m}$(9) ), (15(15) $\frac{\mathrm{\partial }\mathbf{M}(\mathbf{x},t)}{\mathrm{\partial }t}=-L\frac{\delta E}{\delta \mathbf{M}(\mathbf{x},t)}$(15) ) and boundary conditions (3), (11), (25) can be derived. Among them, there is a relationship

(35) $\eta =\frac{m_s^2}{L}$(35)

Therefore, the total potential energy given in (34) corresponds to these basic equations and boundary conditions. The condition for the extreme value of the functional is

(36) $\delta \mathrm{\Pi }=0$(36)

This results in a weak form of equivalent integral for the problem:

(37) $\begin{array}{rl}& {\int }_V\left[\frac{\mathrm{\partial }E}{\mathrm{\partial }\text{varepsilonr}}:\delta \text{varepsilonr}+\left(\frac{\mathrm{\partial }E}{\mathrm{\partial }\mathbf{m}}+\eta \frac{\mathrm{\partial }\mathbf{m}}{\mathrm{\partial }t}\right)\cdot \delta \mathbf{m}+\frac{\mathrm{\partial }E}{\mathrm{\partial }\mathrm{\nabla }\mathbf{m}}:\delta \mathrm{\nabla }\mathbf{m}\right]dV+{\int }_{\mathrm{\Omega }}\frac{\mathrm{\partial }E}{\mathrm{\partial }{\mathbf{h}}_d}\cdot \delta {\mathbf{h}}_{d}dV\\ & ={\int }_A\mathbf{t}\cdot \delta \mathbf{u}dA\end{array}$(37)

In addition, it is easy to obtain the equivalent integral weak form in (37) by properly processing and adding the variational forms of (5), (14), (27).

Finite element equation

In this work, a multi-field coupled finite element was developed via the ABAQUS user subroutine module. In order to establish the finite element equation, the field variable of the node is first determined, that is, the general displacement variable is as follows [20]:

(38) $\mathbf{U}=\left\{\begin{array}{c}{\left[\mathbf{u},\varphi ,\mathbf{m}\right]}^T(V)\\ \left[\varphi \right](\mathrm{\Omega }/V)\end{array}\right.$(38)

The field variables of the magnetic element include displacement, the magnetic scale potential, and the magnetization direction vector, and the field variables of the air element contain only the magnetic scale potential. The general strain corresponding to the general displacement can be represented [20]

(39) $\mathbf{C}=\left\{\begin{array}{c}{\left[\text{varepsilonr},{\mathbf{h}}_d,\mathbf{m},\mathrm{\nabla }\mathbf{m}\right]}^T(V)\\ \left[{\mathbf{h}}_d\right](\mathrm{\Omega }/V)\end{array}\right.$(39)

general stress can be expressed as

(40) $\mathbf{S}=\left\{\begin{array}{c}{\left[\text{sigmab},-\mathbf{b},\frac{\mathrm{\partial }E}{\mathrm{\partial }\mathbf{m}},\frac{\mathrm{\partial }E}{\mathrm{\partial }\left(\mathrm{\nabla }\mathbf{m}\right)}\right]}^T(V)\\ \left[-\mathbf{b}\right](\mathrm{\Omega }/\mathrm{V})\end{array}\right.$(40)

The element node displacement matrix can be represented [20]

(41) ${\mathbf{d}}_e=\left\{\begin{array}{c}{\left[{\mathbf{u}}^1,{\varphi }^1,{\mathbf{m}}^1,{\mathbf{u}}^2,{\varphi }^2,{\mathbf{m}}^2,\dots ,{\mathbf{u}}^I,{\varphi }^I,{\mathbf{m}}^I\right]}^T(V)\\ {\left[{\varphi }^1,{\varphi }^2,\dots ,{\varphi }^I\right]}^T(\mathrm{\Omega }/V)\end{array}\right.$(41)

where $I$ is the number of element nodes. The generalized displacement field variables can be obtained by interpolating the nodes [20]

(42) $\mathbf{U}=\mathbf{N}{\mathbf{d}}_e$(42)

The expression that can be generally strained is obtained accordingly [20]

(43) $\mathbf{C}=\mathbf{B}{\mathbf{d}}_{\mathbf{e}}$(43)

Based on the above definition, the equivalent integral weak form (37) can be rewrite into the following matrix form [20]:

(44) $\sum _e{\int }_{{\mathrm{\Omega }}_e}\delta {\mathbf{C}}^T\mathbf{S}dV=\sum _e\delta {\mathbf{d}}_e^T{\int }_{{\mathrm{\Omega }}_e}{\mathbf{B}}^T\mathbf{S}dV=\sum _e\delta {\mathbf{d}}_e^T{\int }_{\mathrm{\partial }{\mathrm{\Omega }}_e}\mathbf{T}dA$(44)

Suppose ${\mathbf{d}}_{\mathbf{e}}=\mathbf{Q}\mathbf{d}$, where $\mathbf{d}$ is the overall nodal displacement matrix, then (44) can be rewritten as [20],

(45) $\mathbf{R}=\sum _e{\mathbf{Q}}^T\left({\int }_{{\mathrm{\Omega }}_e}{\mathbf{B}}^T\mathbf{S}dV-{\int }_{\mathrm{\partial }{\mathrm{\Omega }}_e}\mathbf{T}dA\right)=\mathbf{0}$(45)

which is the residual equation of the coupled problem and can be solved by the Newton iteration method [20]:

(46) ${\mathbf{R}}^{\left(k+1\right)}={\mathbf{R}}^k+{\left.\frac{\mathrm{\partial }\mathbf{R}}{\mathrm{\partial }\mathbf{d}}\right|}^k\mathrm{\Delta }\mathbf{d}=\mathbf{0}$(46)

The label $k$ represents the $k$th iteration step. The above formula can be written as [20]

(47) ${\mathbf{K}}^k\mathrm{\Delta }\mathbf{d}=-{\mathbf{R}}^k$(47)

where [20]

(48) $\mathbf{K}=\frac{\mathrm{\partial }\mathbf{R}}{\mathrm{\partial }\mathbf{d}}=\sum _e{\mathbf{Q}}^T\left({\int }_{{\mathrm{\Omega }}_e}{\mathbf{B}}^T\mathbf{D}\mathbf{B}dV\right)\mathbf{Q}$(48)

(49) $\sigma =\mathbf{D}\cdot \epsilon$(49)

$\mathbf{D}$ is the tangent stiffness matrix. Thus, the form of the element residual matrix and the stiffness matrix can be obtained [20]:

(50) ${\mathbf{R}}_e={\int }_{{\mathrm{\Omega }}_e}{\mathbf{B}}^T\mathbf{S}dV-{\int }_{\mathrm{\partial }{\mathrm{\Omega }}_e}\mathbf{T}dA$(50)

(51) ${\mathbf{K}}_e={\int }_{{\mathrm{\Omega }}_e}{\mathbf{B}}^T\mathbf{D}\mathbf{B}dV$(51)

Through the residual and stiffness matrix of each element, the overall stiffness matrix and nodal load array are assembled, and the nodal value at the next time step can be obtained by iterative solution in one time step. Through time integration, iterative solution for the next time step is performed until the system reaches a stable equilibrium state.

Specific form of 2D element

In this paper, the 2D triangle element is used, and the magnetic element adopts the Hammer integral of four integral points, whose general displacement and general strain can be written as follows [20]:

(52) $\mathbf{U}={\left[u_1,u_2,\varphi ,m_1,m_2\right]}^T$(52)

(53) $\mathbf{C}={\left[{\epsilon }_{11},{\epsilon }_{22},{\gamma }_{12},h_{d1},h_{d2},m_1,m_2,m_{1,1},m_{1,2},m_{2,1},m_{2,2}\right]}^T$(53)

The element node displacement matrix can be represented

(54) ${\mathbf{d}}_{ei}=[u_1^i,u_2^i,{\varphi }^i,m_1^i,m_2^i{]}^T(i=1,\dots ,3)$(54)

The relationship between general displacement, general strain and node variables is [20]

(55) $\mathbf{U}=\sum N_i{\mathbf{d}}_{ei}$(55)

(56) $\mathbf{C}=\sum {\mathbf{B}}_i{\mathbf{d}}_{ei}$(56)

Among them,

(57) ${\mathbf{B}}_i={\left[\begin{array}{c}{N}_{i,x}0N_{i,y}00000000\\ 0N_{i,y}{N}_{i,x}00000000\\ 000-N_{i,x}-N_{i,y}000000\\ 00000N_{i}0N_{i,x}{N}_{i,y}00\\ 000000N_{i}00N_{i,x}{N}_{i,y}\end{array}\right]}^{\mathrm{T}}$(57)

The general stress matrix corresponding to the general strain can be represented

(58) $\mathbf{S}={\left[{\sigma }_{11},{\sigma }_{22},{\sigma }_{12},-b_1,-b_2,\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_1},\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_2},\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_{1,1}},\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_{1,2}},\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_{2,1}},\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_{2,2}}\right]}^T$(58)

The specific expression is shown in (Table 1). For a cubic magnetostrictive material, the deformation associated with the local magnetization is described by the stress–free strain [23]. Thus, for the 2D condition we assume ${\epsilon }_{33}={\epsilon }_{33}^0=\frac{3}{2}{\lambda }_{100}\left(m_3^2-\frac{1}{3}\right)=-\frac{{\lambda }_{100}}{2}$ for the description of stress-free strain. The generalized shear modulus matrix $\mathbf{D}$ is defined by $\sigma =\mathbf{D}\cdot \epsilon$, and the specific form is shown in (Table 2). Since $\mathbf{D}$ is a symmetric matrix, only the upper triangular part of the expression is written here, and the unspecified terms are all zero.

Table 1. Concrete expression of 2D magnetic element S matrix

${\sigma }_{11}=\frac{\mathrm{\partial }E}{\mathrm{\partial }{\epsilon }_{11}}=c_{11}{\epsilon }_{11}+c_{12}{\epsilon }_{22}+c_{12}{\epsilon }_{33}^0-\frac{3}{2}{\lambda }_{100}\left(c_{11}-c_{12}\right)m_1^2$

${\sigma }_{22}=\frac{\mathrm{\partial }E}{\mathrm{\partial }{\epsilon }_{22}}=c_{11}{\epsilon }_{22}+c_{12}{\epsilon }_{11}+c_{12}{\epsilon }_{33}^0-\frac{3}{2}{\lambda }_{100}\left(c_{11}-c_{12}\right)m_2^2$

${\sigma }_{12}=\frac{\mathrm{\partial }E}{\mathrm{\partial }{\epsilon }_{12}}=c_{44}{\gamma }_{12}-3{\lambda }_{111}{c}_{44}{m}_1{m}_2$

$\frac{\mathrm{\partial }E}{\mathrm{\partial }{h}_{d1}}=-{\mu }_0{m}_s{m}_1-{\mu }_0{h}_{d1}$

$\frac{\mathrm{\partial }E}{\mathrm{\partial }{h}_{d2}}=-{\mu }_0{m}_s{m}_2-{\mu }_0{h}_{d2}$

$\frac{\partial E}{\partial m_1}=2K_1{m}_1{m}_2^2-3{\lambda }_{100}\left(c_{11}-c_{12}\right){\epsilon }_{11}{m}_1-6{\lambda }_{111}{c}_{44}{\epsilon }_{12}{m}_2-{\mu }_0{m}_s{h}_{d1}-{\mu }_0{m}_s{h}_{e1}+\eta {\dot{m}}_1+2A_s{m}_1-\frac{2A_s{m}_1}{\sqrt{m_1^2+m_2^2}}$

$\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_2}=2K_1{m}_1^2{m}_2-3{\lambda }_{100}\left(c_{11}-c_{12}\right){\epsilon }_{22}{m}_2-6{\lambda }_{111}{c}_{44}{\epsilon }_{12}{m}_1-{\mu }_0{m}_s{h}_{d2}-{\mu }_0{m}_s{h}_{e2}+\eta {\dot{m}}_2+2A_s{m}_2-\frac{2A_s{m}_2}{\sqrt{m_1^2+m_2^2}}$

$\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_{1,1}}=2A{m}_{1,1}$

$\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_{1,2}}=2A{m}_{1,2}$

$\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_{2,1}}=2A{m}_{2,1}$

$\frac{\mathrm{\partial }E}{\mathrm{\partial }{m}_{2,2}}=2A{m}_{2,2}$

Table 2. 2D magnetic element D matrix upper triangular expression

$D_{11}=c_{11},D_{12}=c_{12},D_{16}=-\frac{3{\lambda }_{100}}{2{M_s}^2}(c_{11}-c_{12})M_1$

$D_{21}=c_{12},D_{22}=c_{11},D_{27}=-\frac{3{\lambda }_{100}}{2{M_s}^2}(c_{11}-c_{12})M_2$

$D_{33}=c_{44},D_{36}=-\frac{3{\lambda }_{111}}{2{M_s}^2}{c}_{44}{M}_2,D_{37}=-\frac{3{\lambda }_{111}}{2{M_s}^2}{c}_{44}{M}_1$

$D_{44}=-{\mu }_0,D_{46}=-{\mu }_0$

$D_{55}=-{\mu }_0,D_{57}=-{\mu }_0$

$D_{66}=\frac{K_1}{{M_s}^4}2{M_2}^2-\frac{3{\lambda }_{100}}{2{M_s}^2}(c_{11}-c_{12}){\epsilon }_{11}+2A_s-\frac{2A_s{M}_s}{\sqrt{{M_1}^2+{M_2}^2}}+\frac{\eta }{\tau }$

$D_{67}=-\frac{3{\lambda }_{111}}{2{M_s}^2}{c}_{44}{\gamma }_{12}$

$D_{77}=\frac{K_1}{{M_s}^4}2{M_1}^2-\frac{3{\lambda }_{100}}{2{M_s}^2}(c_{11}-c_{12}){\epsilon }_{22}+2A_s-\frac{2A_s{M}_s}{\sqrt{{M_1}^2+{M_2}^2}}+\frac{\eta }{\tau }$

$D_{88}=D_{99}=D_{\left(10\right)\left(10\right)}=D_{\left(11\right)\left(11\right)}=2A$

For the air element outside the ferromagnetic body, the corresponding general displacement, generalized strain, and generalized stress matrix can be represented

(59) $\mathbf{U}=\left[\varphi \right],\mathbf{C}={\left[h_{d1},h_{d2}\right]}^T,\mathbf{S}={\left[\frac{\mathrm{\partial }E}{\mathrm{\partial }{h}_{d1}},\frac{\mathrm{\partial }E}{\mathrm{\partial }{h}_{d2}}\right]}^T$(59)

Its node displacement matrix $\mathbf{d}$ and matrix $\mathbf{B}$ are

(60) ${\mathbf{d}}_{ei}=\left[{\varphi }^i\right](i=1,3),{\mathbf{B}}_i={\left[-N_{i,x}-N_{i,y}\right]}^T$(60)

The tangent stiffness matrix can be represented

(61) $\mathbf{D}=\left[\begin{array}{c}-{\mu }_{0}0\\ 0-{\mu }_0\end{array}\right]$(61)

The magnetic potential $\varphi$ remains continuous on the interface between the magnetic element and the air element, on the outer boundary of the air element $\varphi =0$. Only the magnetic potential is used as a variable for the magnetic field of the air element. By coupling the magnetic element to the air node on the boundary, the potential is continuous. In the simulation process, the backward Euler method is used for time integration.

Magnetic domain evolution mechanism of double-hole ferromagnetic nanorings

Through the written user 2D element in ABAQUS software, the magnetic domain evolution mechanism of the double-hole nanoring under a magnetic field is studied, the evolution mechanism of the magnetic domain of the double-hole ferromagnetic nanoring with the geometric aperture and the direction of the external magnetic field is analyzed, and its application prospect on random storage units is explored.

Simulation model and parameter settings

The two-hole ferromagnetic nanoring calculation model is shown in Figure 1. The magnetic easy axes are along x-axis directions. Where the outer diameter of the nanoring is D = 0.1 μm, the two round holes are symmetrical to the left and right, both in diameter d = 0.025 μm, and the center distance is 0.05 μm. The direction of uniform magnetic field application and the angle of y-axis direction are recorded, the hysteresis return line in different hole aperture and direction of different external magnetic field is calculated, and the evolution of magnetization state is examined. The calculation parameters of the nanorings are based on the FeGa material parameters [23] with the Zener anisotropy factor $\tilde{A}=\frac{2C_{44}}{C_{11}-C_{12}}=6.15$ [24,25], as (Table 3) shows. The thickness of the simulated magnet has strong influence on the micromagnetic structure especially in the thin film. However, we consider it that since the magnetic easy axes of FeGa material are along <100> directions [23], we performed the simulations in essentially 2D systems. There are several researches obey the rules in computing in 2D simulations. In 2D simulations, the field parameters in the z-axis consider to be uniform, which contains no information of the thickness. The further effects of thickness could be discussed in our future work.

Figure 1. Schematic diagram of calculation model of symmetric double-hole nanoring

Table 3. Material parameters (FeGa) [23] for the simulation

$C_{11}(N/\mu m^2)$	$C_{12}(N/\mu m^2)$	$C_{44}(N/\mu m^2)$	$K_1(N/\mu m^2)$	${\lambda }_{100}$	${\lambda }_{111}$
0.196	0.156	0.123	2×10^−8	2.64×10^−4	0
$m_s(A/\mu m)$	${\mu }_0(N/A^2)$	$\eta (N\cdot s/\mu m^2)$	$A_s(N/\mu m^2)$
1.432	$4\pi \times {10}^{-7}$	2.0×10^−13	5×10^−6

The revolution mechanism of the geometric aperture and the direction of the external magnetic field on the magnetic domain of the double-hole ferromagnetic nanoring

By defining the aperture ratio:

(62) $\alpha =\frac{d}{D},$(62)

changing the geometric aperture ratio of the ring$\alpha$ and the angle of the external magnetic field $\theta$, through static cyclic loading by magnetic field, we calculate the hysteresis return line under different regulatory parameters, and observe the evolution of the magnetic field. The magnetization direction is randomly oriented in the initial calculation. (Figure 2) shows a summary diagram of the evolution mechanism of symmetrical double-hole ferromagnetic nanorings. The horizontal coordinates in the figure are the angle between the direction of the applied magnetic field and the direction of the $y$ axis, and the ordinate is the aperture ratio $\alpha$. In the calculation parameter grid of aperture ratio $\alpha$ from 0.10 to 0.40 with range interval 0.05, applied external magnetic field angle $\theta$ from 0° to 90° with range interval 10°, we found a total of 7 evolutionary mechanisms, namely: twin-vortex mechanism, left and right vortex mechanism, ‘S’ evolution mechanism, ‘S’ + ‘C’ mechanism, ‘S’ + twin-vortex mechanism, ‘S’ + left vortex mechanism, and single left vortex mechanism.

Figure 2. Summary of the evolution mechanism of the symmetric double-hole ferromagnetic nanoring (TV: twin-vortex mechanism, L+ RV: left and right vortex mechanism, S: ‘S’ evolution mechanism, S + C: ‘S’ + ‘C’ mechanism, S+ TV: ‘S’ + twin-vortex mechanism, S+ LV: ‘S’ and left vortex mechanism, and LV: left vortex mechanism)

Twin-vortex mechanism

We found a new stable vortex state in the regulatory parameters: aperture ratio $\alpha =0.20$, applied external magnetic field angle $\theta$ from 40° to 90° range, twin-vortex state (TV), as shown in (Figure 3). The twin-vortex magnetic domain structure is similar to but different from the vortex closed domain structure of a single-hole ring: the same vortex-closed domain structure as the single-hole ring with a small demagnetizing field; the vortex domain structure of the single-hole ring surrounds the hole to form a stable vortex domain structure, but the double-hole ring twin-vortex domain structure surrounds the two holes to form two vortex domain structures with opposite rotation directions.

Figure 3. Twin-vortex magnetization state diagram

For twin-vortex magnetization state, Nipun A [26] in 2007 demonstrated twin-vortex in 30-nm-thick Co nanopatterned double pac-man (DPM) elements (with Outer Diameter (OD): 400 nm, notch radius: 120 nm, and Slot Angle 60°), fabricated by electron-beam lithography (EBL) and liftoff processing, when the field is applied in an orthogonal direction, using Lorentz microscopy and off-axis electron holography.

(Figure 4) shows the hysteresis loop of the twin-vortex evolution mechanism. The horizontal coordinates in the figure are the size of the external magnetic field $\mathbf{H}$, and the ordinates are the ratio of the size of the magnetization vector $\mathbf{m}$ to the intensity of the saturation magnetization $m_s$. The external magnetic field of the double-hole ring is greater than the strength of the positive saturated magnetic field and less than the strength of the reverse saturation magnetic field is in the onion state (ON), when the magnetic field decreases the magnetization state of the ring evolves from the onion state to the twin-vortex state (TV), i.e., there is an ON-TV mechanism, similar to the onion state (ON) – vortex state (VT) mechanism of the single-hole ring. Due to the closure of the twin-vortex magnetic field, when the external magnetic field is reduced to zero, the demagnetization field is very small, which retains the advantages of the closed single-hole ring vortex magnetic domain, so that the adjacent magnetic ring does not interfere with each other directly, which has great advantages for improving the density of the storage unit.

Figure 4. Twin-vortex mechanism hysteresis loop

In this paper, the critical flip field of the onion state-vortex state is called the vortex-vortex nuclear field $H_n$, and the critical flip field of the vortex-onion state is called the vortex vanishing field $H_{an}$, then

(63) $\mathrm{\Delta }H=H_n-H_{an}$(63)

is the magnetic field range of the twin-vortex state, the size of which shows the stability of the twin vortex state under the change of the magnetic field: the larger $\mathrm{\Delta }H$, the more stable the twin-vortex state; vice versa.

Figure 5. Variation curve of $H_n$, $H_{an}$, $\mathrm{\Delta }H$ under different applied magnetic field angle $\theta$

(Figure 5) shows the vortex nucleation field $H_n$, vortex vanishing field$H_{an}$, and magnetic field range $\mathrm{\Delta }H$ in the twin-vortex state as a function of the applied magnetic field angle $\theta$. The vortex nuclear field $H_n$ decreases with the increase of the magnetic field application angle $\theta$, the vortex vanishing field $H_{an}$ increases with the increase of the magnetic field application angle $\theta$, and the magnetic field range $\mathrm{\Delta }H$ decreases with the increase of the applied magnetic field angle $\theta$. That is, the stability of the twin-vortex state decreases as the magnetic field application angle $\theta$ increases.

Figure 6. Left and right vortex magnetization state diagram

Left and right vortex mechanism

The control parameters of the left and right vortex have a relatively wide range of values: the aperture ratio $\alpha$ range from 0.10 to 0.40, the angle of the applied external magnetic field θ from 0° to 60°, which is a relatively common evolutionary mechanism in the double-hole ring, as (Figure 2) shown.

The left and right vortex mechanism is: when the magnetic field is large, it is in a positive saturation state; when the magnetic field is reduced to a certain value, the magnetization state of (Figure 6(a)) will appear, that is, the vortex structure surrounding the small ring on the left in a counterclockwise direction. At this time, there is a positive net magnetization; when the reverse magnetic field increases, a vortex structure surrounding the small ring on the right will appear as a whole counterclockwise, as shown in (Figure 6(b)), at this time there is a reverse net magnetization; When the magnetic field increases to a certain value, the reverse saturation magnetization state will appear; when the reverse magnetic field decreases, the vortex structure surrounding the left small ring and the right small ring will appear in sequence, but the rotation direction is opposite, changing from counterclockwise It becomes clockwise, as shown in (Figure 6(c,d)). The hysteresis loop diagram of the left and right vortex evolution mechanism is shown in (Figure 7). We also mark the vortex nucleation field $H_n$ and the vortex vanishing field $H_{an}$, and observe the changes of the vortex nucleation field $H_n$, vortex vanishing field $H_{an}$, and magnetic field range $\mathrm{\Delta }H$ with the applied magnetic field angle $\theta$ and aperture ratio $\alpha$.

For Left and right vortex states, Rahm M [27]. in 2004 demonstrated the existence of two metastable vortex states with the vortex center either pinned at the left or right antidot in disks consist of Permalloy, 30 nm thick, with a 10 nm protective Ti layer, with diameters of 500 nm, containing two antidots with center-to-center distances d of 150 and 200 nm, measured with micro-Hall magnetometry.

(Figure 8) shows the vortex nucleation field $H_n$, vortex vanishing field $H_{an}$, and magnetic field range $\mathrm{\Delta }H$ in the left vortex state as a function of the applied magnetic field angle $\theta$. The vortex nucleation field $H_n$ decreases with the increase of the magnetic field angle $\theta$, and the vortex vanishing field $H_{an}$ also decreases with the increase of the magnetic field angle $\theta$. The magnetic field range $\mathrm{\Delta }H$ is not greatly affected by the change of the applied magnetic field angle $\theta$, that is, the stability of the left vortex state is not greatly affected by the change of the applied magnetic field angle $\theta$.

(Figure 9) shows the vortex nucleation field $H_n$, vortex vanishing field $H_{an}$, and magnetic field range $\mathrm{\Delta }H$ in the left vortex state as a function of the geometric aperture ratio $\alpha$. The vortex nucleation field $H_n$ first decreases and then increases as the aperture ratio $\alpha$ increases. The vortex vanishing field $H_{an}$ also decreases first and then increases with the aperture ratio $\alpha$. The magnetic field range $\mathrm{\Delta }H$ is not greatly affected by the change in the geometric aperture ratio $\alpha$, namely the stability of the left vortex state is not greatly affected by the change of the geometric aperture ratio $\alpha$.

In short, the stability of the left vortex state is not affected by the change in the angle of the applied magnetic field $\theta$ and the ratio of geometric aperture $\alpha$.

Figure 7. Hysteresis loop of left and right vortex evolution mechanism

Figure 8. Variation curve of $H_n$, $H_{an}$, $\mathrm{\Delta }H$ under different applied magnetic field angle $\theta$ ($\alpha =0.25$)

Figure 9. Variation curve of $H_n$, $H_{an}$, $\mathrm{\Delta }H$ under different aperture ratio $\alpha$ ($\theta =0^{\circ }$)

‘S’ evolutionary mechanism

The ‘S’ evolutionary mechanism regulates the distribution of parameters in small areas, only in cases where the aperture $\alpha$ is quite large and relatively small ($\alpha =0.10\sim 0.20,0.35\sim 0.40,\theta ={20}^{\circ }\sim {30}^{\circ }$), and the angle of the magnetic field $\theta ={40}^{\circ }\sim {30}^{\circ }$ is applied.

Figure 10. ‘S’ state magnetization state diagram

Figure 11. Hysteresis loop of ‘S’ evolution mechanism

The ‘S’ state evolution mechanism includes a new state of the double-hole ring: the ‘S’ state. The magnetization state of the ‘S’ state is that the two-hole ring contains two 180° domain walls, forming a magnetization vector from the right side of the right hole around the right hole to the middle and then around the left hole, and finally flowing out on the left side of the left hole (flowing from right to left, and vice versa from left to right), as shown in (Figure 10).

(Figure 11) shows the hysteresis loop of ‘S’ evolution. When the forward magnetic field is unloaded, the magnetization state of the double-hole ring changes from the positive saturation state to the ‘S’ state flowing from right to left. When the magnetic field continues to unload and enters the reverse loading, the magnetization state of the double-hole ring turns from the ‘S’ state magnetization flowing from right to left to the ‘S’ state flowing from left to right, and the magnetization state of the reverse magnetic field continues to increase and it will enter the reverse saturation state.

The ‘S’ + ‘C’ mechanism

The ‘S’ + ‘C’ evolution mechanism has a relatively small regulatory parameter area, which appears when the aperture ratio $\alpha$ is relatively large ($\alpha =0.30,0.40$) and the applied magnetic field angle $\theta ={40}^{\circ }\sim {60}^{\circ }$ . The evolution path of the ‘S’ + ‘C’ mechanism consists of the ‘S’ state and the ‘C’ state. The magnetization state of the ‘C’ state is shown in (Figure 12).

Figure 12. ‘C’ state magnetization state diagram

For ‘S’ and ‘C’ states, Jonathan K H [28] in 2003 simulated permalloy disks and observed in-plane magnetization buckle states, including first-order buckle state, C state, observed in a diameter of 200 nm and a diameter-to-thickness ratio d/t of 20 at 0-mT field, as well as second-order buckle state, S state, observed in a diameter of 500 nm and a diameter-to-thickness ratio d/t of 3 at 3-mT field.

(Figure 13) shows the ‘S’ + ‘C’ mechanism evolutionary hysteresis loop. When a large positive magnetic field is applied, the forward magnetization of the double-hole magnetic ring is saturated, and as the magnetic field decreases, it enters the ‘S’ flowing from right to left. When the magnetic field is reduced to zero and a reverse magnetic field is applied, the ‘S’ state does not flip and jumps to the counterclockwise ‘C’ state. The reverse magnetic field continues to increase, and the magnetic ring changes from the counterclockwise ‘C’ state to the reverse saturation state.

‘S’ + twin-vortex mechanism

The ‘S’ + twin-vortex mechanism only appears when the aperture ratio $\alpha =0.10$ and the applied magnetic field angle $\theta ={60}^{\circ }$, 80°. (Figure 14) shows the hysteresis loop diagram of ‘S’ + twin-vortex mechanism. Like the ‘S’ + ‘C’ mechanism, when the positive magnetic field is large, the positive magnetization of the double-hole magnetic ring is saturated. As the magnetic field decreases, it enters the ‘S’ state flowing from right to left. When the reverse magnetic field is applied, the ‘S’ state does not turn over and jumps to the twin-vortex magnetization state of the upward state. The reverse magnetic field continues to increase. The magnetic ring enters the reverse saturation state from the twin-vortex magnetization state of the upward state.

Figure 13. ‘S’ + ‘C’ mechanism evolution hysteresis loop diagram ($\alpha =0.40,\theta ={40}^{\circ }$)

Figure 14. Hysteresis loop diagram of ‘S’ + twin-vortex mechanism (α = 0.10)

‘S’ + left vortex mechanism

The ‘S’ + left vortex mechanism appears when the aperture ratio α is larger or smaller, $\alpha =0.15,0.30\sim 0.40$, and the applied magnetic field angle $\theta ={60}^{\circ }\sim {90}^{\circ }$. The hysteresis loop diagram of the ‘S’ + left vortex mechanism is shown in (Figure 15). Like the ‘S’ + ‘C’ mechanism, when the positive magnetic field is large, the positive magnetization of the double-hole magnetic ring is saturated. As the magnetic field decreases, it enters the ‘S’ state that flows from right to left. When the magnetic field is reduced to zero and a reverse magnetic field is applied, the ‘S’ state does not flip and jumps to a counterclockwise left vortex state, and the reverse magnetic field continues to increase. The magnetic ring changes from a counterclockwise left vortex state to reverse saturation state.

Figure 15. ”S” + Left vortex mechanism hysteresis loop diagram (α = 0.40)

Figure 16. Hysteresis loop diagram of left vortex mechanism

Left vortex mechanism

The left vortex mechanism only appears when the aperture ratio $\alpha =0.15$ and the angle $\theta ={80}^o$ at which the magnetic field is applied. The hysteresis loop diagram of the left vortex mechanism is shown in (Figure 16). When the applied forward magnetic field is large, the two-hole magnetic ring is saturated with positive magnetization, as the magnetic field decreases into the counterclockwise left vortex state, when the magnetic field decreases to zero and the reverse magnetic field is applied, the counterclockwise left vortex state does not flip to the counterclockwise right vortex state, but continues to increase with the reverse magnetic field into the reverse saturation state.

The magnetoelastic coupling has a significant effect on magnetization especially in multiferroic [29,30] and giant magnetostrictive materials [31,32]. Using the formula given by Zhang J X [23] in 2005, with the magnetic elastic coupling energy items taken into consideration, the homogeneous strain can be roughly evaluated as:

(64) $\begin{array}{rl}& {\bar{\epsilon }}_{11}=\frac{3}{2}{\lambda }_{100}\left(\overline{m_1^2}-\frac{1}{3}\right),\\ & {\bar{\epsilon }}_{22}=\frac{3}{2}{\lambda }_{100}\left(\overline{m_2^2}-\frac{1}{3}\right),\\ & {\bar{\epsilon }}_{12}=\frac{3}{2}{\lambda }_{111}\overline{m_1{m}_2},\end{array}$(64)

where $\overline{m_i^2}$ and $\overline{m_i{m}_j}$ represent the volume average of $m_i^2$ and $m_i{m}_j$ respectively. However, we consider it that in this paper the revolution of the magnetization is of key focus. The further effects of magnetoelastic coupling could be discussed in our future work.

The energy explanation of the evolutionary process of the magnetic domain

By extracting the energy density of various energies from the evolutionary process of the twin vortex mechanism, the left and right vortex mechanism and the ‘S’ mechanism, we analyze the demagation energy and exchange energy, investigate the change process of various energies in the evolutionary process, and explain the evolutionary process. We compare the twin-vortex reversal demagnetization energy with single hole ring (the magnetic field is gradually reduced from the forward maximum), as shown in (Figure 17). Blue line is a two-hole ring twin-vortex flip and green is a single-hole ring flip. The single-hole ring and the double-hole ring form a closed vortex structure with a small magnetic field, which reduces the demagation energy, and the magnetic field of the twin-vortex has a wider range and is less affected by the external magnetic field, but a larger magnetic field needs to be applied when flipping. We compare twin-vortex flip exchange energy with single-hole ring, as shown in (Figure 18). When the magnetic domain flips, the exchange energy will jump, indicating that the magnetic domain structure has changed.

Figure 17. Comparison of demagnetization energy of double-hole ring twin-vortex reversal mechanism and single-hole vortex reversal mechanism demagnetization energy (Blue line is double-hole ring and twin-vortex flip, green is single-hole ring vortex flip, the same below)

Figure 18. Comparison of the exchange energy of the double-hole ring twin-vortex flip mechanism and the exchange energy of the single-hole vortex flip mechanism

We also compare left and right vortex reversal demagnetization energy with single hole ring, shown in (Figure 19). The energy of the left vortex and the right vortex of the double-hole magnetic ring is not much different, but adding holes makes the vortex state stay longer on the magnetic ring and the magnetic field range of the vortex state becomes larger. We compare exchange energy between left and right vortex flipping and single hole ring in (Figure 20). The exchange energy of the double-hole magnetic ring undergoes three jumps during the unloading process of the magnetic field, and the magnetic domain structure undergoes three sudden changes, corresponding to the three jumps experienced by the four stable states of the hysteresis loop. Since the two 180° domain walls of the ‘S’ state rotate synchronously with the magnetic field loading on the edge of the ring, the magnetic domain area changes, so the magnetization of the ‘S’ evolution hysteresis loop changes with the magnetic field unloading. The magnetic field decreases, because the energy of the domain wall increases with the increase of the ring width. When rotating to a certain angle, the energy of the domain wall causes a sudden magnetization reversal, which turns from right to left to flow from left to right.

Figure 19. Comparison of the demagnetization energy of the double-hole ring left and right vortex reversal mechanism and the demagnetization energy of the single-hole vortex reversal mechanism (Blue is the double-hole ring with left and right vortex inversion, green is single-hole ring with vortex inversion, the same below)

Figure 20. Comparison of the exchange energy between the left and right vortex flip mechanism of the double-hole ring and the exchange energy of the single-hole vortex flip mechanism

Application prospect of double-hole nanorings on random storage

Zhu J G [33] in 2008 reviewed the different types of MRAM designs. Among them, there is a ring-based MRAM design, avoiding affection of geometric variation of the ends in actual fabrication on variation of the demagnetization field within the element. In the ring-based MRAM design, the magnetostatic flux forms a perfect closure, and the entire element has no demagnetizing field. In order to read the storage state, that is, the magnetization chirality in the storage layer, a bidirectional current pulse is injected into the ring element. The circular current field generated by the current will reverse the direction of the magnetization chirality in the integration according to the current layer. The magnetization chirality of the storage layer can then be determined by the corresponding magnetoresistance changes of the elements.

In this paper, we compare the magnetic field range $\mathrm{\Delta }H$ of the double-hole ring twin-vortex state (TV) with the magnetic field range $\mathrm{\Delta }H$ of the single-hole ring vortex state (VT) with the same outer diameter and different ring widths, as shown in (Figure 21), and find that the magnetic field range $\mathrm{\Delta }H$ of the double-hole ring twin-vortex state (TV) is greater than the magnetic field range $\mathrm{\Delta }H$ of the single-hole ring vortex state (VT), that is, the double-hole ring twin-vortex state (TV) is more stable than the single-hole ring vortex state (VT). The higher stability of the twin-vortex state (TV) of the double-hole ring makes the double-hole ring have greater advantages in high-density storage applications, such as: stored data is less likely to be lost and storage retention time is longer. Thus, twin-vortex evolution mechanism has the advantages of high stability and low demagnetization interference characteristics, and has good application prospects in magnetic random-access memory (MRAM).

Figure 21. Comparison of vortex stability of double-hole ring with twin-vortex and single-hole ring

Conclusion

In this paper, through the user subroutine element combining finite element and phase field, the magnetic domain evolution mechanism of the double-hole nanoring is studied, focusing on the mechanism of the magnetic domain affected by the geometric aperture and external magnetic field direction. Under different geometric apertures and external magnetic field direction parameters, we have found seven ways of magnetic domain evolution: twin-vortex mechanism, left and right vortex mechanism, ‘S’ evolution mechanism, ‘S’ + ‘C’ mechanism, ‘S’ + twin-vortex mechanism, ‘S’ + left vortex mechanism, left vortex mechanism. Due to the advantages of high stability and low demagnetization mutual interference, the twin-vortex mechanism has a good application prospect in random memory cells.

Acknowledgments

The authors acknowledge the support by the National Natural Science Foundation of China (Grant Nos. 12025201, 11521202, 11890681, 11522214). Calculations are supported by High-Performance Computing Platform of Peking University, China.

Disclosure statement

The authors declare that they have no conflict of interest.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [12025201, 11521202, 11890681, 11522214].