Analytical computation of the demagnetizing energy of thin film domain walls

TL;DR

This paper presents an analytical method to compute the demagnetizing energy of a straight domain wall in thin ferromagnetic films, improving understanding and modeling of domain wall properties.

Contribution

It introduces an analytical calculation of demagnetizing energy for a $ anh$-profile domain wall, enhancing the 1D model of domain wall dynamics.

Findings

Derived accurate expressions for domain wall width and Walker breakdown field.

Validated analytical results with micromagnetic simulations.

Improved modeling of domain wall behavior in thin films.

Abstract

Due to its non-local nature, calculating the demagnetizing field remains the biggest challenge in understanding domain structures in ferromagnetic materials. Analytical descriptions of demagnetizing effects typically approximate domain walls as uniformly magnetized ellipsoids, neglecting both the smooth rotation of magnetization from one domain to the other and the interaction between the two domains. Here, instead of the demagnetizing field, we compute analytically the demagnetizing energy of a straight domain wall described by the classical $\tanh$ magnetization profile in a thin film with perpendicular magnetic anisotropy. We then use our expression for the demagnetizing energy to derive an improved version of the 1D model of field-driven domain wall motion, resulting in accurate expressions for important properties of the domain wall such as the domain wall width and the Walker breakdown field. We verify the accuracy of our analytical results by micromagnetic simulations.

Equations

$$\frac{\partial \mathbf{m}}{\partial t}+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}=\gamma \mathbf{m}\times {\mathbf{B}}_{\text{eff}},$$

$${\mathbf{B}}_{\text{eff}}=-\frac{1}{M_s}\frac{\delta E}{\delta \mathbf{m}}.$$

$${\mathbf{H}}_d=-M_s\left(N_x{m}_x{\mathbf{e}}_x+N_y{m}_y{\mathbf{e}}_y+N_z{m}_z{\mathbf{e}}_z\right),$$

$$N_x^E\approx \frac{\delta }{\delta +w},N_y^E\approx \frac{\delta }{\delta +D}.$$

$$\mathbf{m}(\mathbf{r})=\tanh \left(\frac{y-Q}{D}\right){\mathbf{e}}_z+\frac{{\mathbf{e}}_x\cos \varphi +{\mathbf{e}}_y\sin \varphi }{\cosh \left(\frac{y-Q}{D}\right)},$$

$$E_d=-\frac{{\mu }_0{M}_s}{2}\int \mathbf{m}\cdot {\mathbf{H}}_d{d}^3\mathbf{r},$$

$${\mathbf{H}}_d(\mathbf{r})$$

$${\mathbf{H}}_d^V(\mathbf{r})$$

$${\mathbf{H}}_d^S(\mathbf{r})$$

$$E_d$$

$$E_d^V$$

$$E_d^S$$

$$g_{\delta }(r)$$

$$=2\delta \mathrm{asinh}\frac{\delta }{r}+2r-2\sqrt{r^2+{\delta }^2}.$$

$$f_{\delta }(r)$$

$$\nabla \cdot \mathbf{m}=\frac{\partial m_y}{\partial y}=-\sin \varphi \frac{\sinh \frac{y-Q}{D}}{D{\cosh }^2\frac{y-Q}{D}}.$$

$$E_d^V={\mu }_0{M}_s^2{\sin }^2\varphi \frac{D^3}{8\pi }{\int }_{-\frac{w}{2D}}^{\frac{w}{2D}}{\int }_{-\frac{w}{2D}}^{\frac{w}{2D}}{\int }_{-\frac{h}{2D}}^{\frac{h}{2D}}{\int }_{-\frac{h}{2D}}^{\frac{h}{2D}}\frac{\sinh (y-q)\sinh (y^{\prime }-q)}{{\cosh }^2(y-q){\cosh }^2(y^{\prime }-q)}{g}_{\sigma }(\sqrt{(x-x^{\prime }{)}^2+(y-y^{\prime }{)}^2})dxd{x}^{\prime }dyd{y}^{\prime },$$

$$u=x-x^{\prime },$$

$$v=y-y^{\prime },$$

$${\mathcal{E}}_V=$$

$$G_{\sigma }(v)$$

$${\mathcal{E}}_V=-\frac{D^2}{2\pi }{\int }_0^{\infty }\frac{v\cosh v-\sinh v}{{\sinh }^{2}v}\frac{\partial G_{\sigma }(v)}{\partial v}dv.$$

$$\frac{\partial G_{\sigma }}{\partial v}=4\sigma \mathrm{atan}\frac{v}{\sigma }-2\sigma \pi -4v\ln \frac{v}{\sqrt{v^2+{\sigma }^2}}.$$

$$\frac{\partial {}^2{G}_{\sigma }}{\partial v^2}=-4\ln \frac{v}{\sqrt{v^2+{\sigma }^2}}=2\ln \left(1+\frac{{\sigma }^2}{v^2}\right),$$

$${\mathcal{E}}_V=\frac{D^2}{\pi }\left[\sigma \pi -{\int }_0^{\infty }\frac{v}{\sinh v}\ln \left(1+\frac{{\sigma }^2}{v^2}\right)dv\right].$$

$${\mathcal{E}}_V$$

$$=\delta D+\frac{D^2}{\pi }\mathrm{Re}{I}_{\sigma },$$

$$I_{\sigma }=\sum _{n=1}^{\infty }\frac{(-i\sigma {)}^n}{n}{I}_n,I_n={\int }_C\frac{v}{\sinh v}{v}^{-n}dv.$$

$$I_1$$

$$I_2$$

$$I_n$$

$$=-2(i\pi {)}^{2-n}\left(1-2^{2-n}\right)\zeta (n-1),$$

$${\mathcal{E}}_V$$

$$=\left(\frac{{\delta }^2}{\pi }\ln 2-\frac{{\delta }^3}{18D}+\frac{3{\delta }^4}{8{\pi }^3{D}^2}\zeta (3)\right)+{\delta }^2\mathcal{O}\left({\sigma }^3\right),$$

$$H_y^d=-M_s{N}_y{m}_y=-M_s{N}_y\frac{\sin \varphi }{\cosh \left(\frac{y-q}{D}\right)}.$$

$$\frac{E_d}{w}$$

Full text

I Introduction

Domain walls (DWs) in low-dimensional ferromagnetic systems such as nanowires and thin films are an active field of study, with promising applications in spintronics such as memory [1] and logic [2] devices. These applications typically rely on magnetic fields [3, 4, 5] or spin-polarized electric currents [6, 7] to drive DW motion. Hence, accurate analytical and numerical descriptions of field and current-driven DW dynamics are essential for future device applications.

The basic description of such magnetic systems start with the Landau-Lifshitz-Gilbert (LLG) equation for the time-evolution of the magnetization vector ${\mathbf{m}}={\mathbf{M}}/M_{s}$, with $M_{s}$ the saturation magnetization, which in the case of field-driven magnetization dynamics reads [8]

[TABLE]

where $\alpha$ is the phenomenological Gilbert damping constant, and $\gamma$ is the gyromagnetic ratio. The effective field $\mathbf{B}_{\text{eff}}$ in Eq. (1) can be formulated in terms of the total energy $E$ of the system,

[TABLE]

The energy contains contributions from the exchange, anisotropy, and Zeeman energies, as well as the demagnetizing energy due to the long-range interaction between magnetization vectors. Numerical solutions of Eq. (1) using a given space discretization are referred to as micromagnetic simulations, and form an important part of studies of DWs and their dynamics.

From an analytical perspective, a class of widely used reduced models of DW dynamics is given by the so-called 1D models, describing the DW in terms of a smoothly varying one-dimensional magnetization profile parameterized by the DW position, width, and an angular variable describing the orientation of the magnetization inside the DW. Dynamical equations for these variables are derived from the LLG equation [9, 10, 11, 12, 13]. As a general feature, such models (as well as the corresponding micromagnetic simulations) exhibit a regime of steady DW dynamics for small applied fields $B_{a}$ with the DW velocity increasing with the field. For fields stronger than a specific driving field magnitude $B_{W}$, the internal magnetization of the DW starts precessing, resulting in an abrupt drop in the DW propagation velocity. This instability is related to the breakdown of the solution found by Schryer and Walker describing the steady field-driven propagation of an infinitely extended planar DW [3], and is referred to as Walker breakdown.

In ferromagnetic systems with reduced dimensions compared to the DW width such as nanowires and (ultra) thin films, demagnetizing effects due to the spatial confinement of the DW become important. The demagnetizing field $\mathbf{B}_{d}=\mu_{0}\mathbf{H}_{d}$, arising from the demagnetizing part of the energy in the expression (2), contains nonlocal contributions from magnetic volume charges $\nabla\cdot\mathbf{m}$ and surface charges $\mathbf{m}\cdot\mathbf{n}$ at the boundary of the system, and gives rise to effects such as shape anisotropy, which penalizes any magnetization normal to the boundary, and the restoring force, which pulls the DW towards the center of the sample to keep the net magnetization neutral. In general, the demagnetizing field poses the biggest challenge to understanding domain structures in magnetic systems due to its long-range nature. A direct computation of the $\mathbf{H}_{d}$ at any given point is often intractable except in very simple cases. One such case is that of a uniformly magnetized ellipsoid, where the demagnetizing field inside the sample can be given as [14]

[TABLE]

where the constants $N_{i}$, $i=x,y,z$, known as the demagnetizing factors, depend on the axes of the ellipsoid in question, and must satisfy $N_{x}+N_{y}+N_{z}=1$. The simplicity of the demagnetizing factors has motivated approximations where the demagnetizing field is assumed to follow the form (3) even when it is strictly speaking not applicable. For example, in order to study the effects of the demagnetizing field on DW motion in thin films, Mougin et. al. [12] modelled the DW as a uniformly magnetized ellipsoid with axes $(w,D,\delta)$, which (with $\delta\ll D<w$) results in demagnetizing factors $N_{i}^{E}$ given by

[TABLE]

While this allows a simple description of demagnetizing fields, it is a rather coarse approximation because it ignores both the rapid variation of magnetization inside the DW as well as the interaction between the two domains and the DW. This directly affects the accuracy of the resulting properties of the DW, such as the Walker breakdown field $B_{W}$ and the DW width.

In this paper, instead of working with the demagnetizing field itself, we compute the energy due to the demagnetizing field of a uniformly magnetized, straight DW in a thin film with perpendicular magnetic anisotropy (PMA). As we shall show, the demagnetizing energy is more analytically tractable than the field, and still lets us derive dynamical equations for the DW using a Lagrangian framework. We assume a straight, infinitely long DW with no variation in the direction $\mathbf{e}_{z}$ normal to the film, which is valid if the film has thickness $\delta\ll D$. The direction of the in-plane magnetization vector, measured by the angle $\phi$ that $\mathbf{m}$ makes with the $x$ axis inside the DW, is taken to be uniform. $\phi=0$ corresponds to a Bloch wall configuration (see Fig. 1a), which is energetically preferred due to the absence of magnetic volume charges $\nabla\cdot\mathbf{m}$. However, applying a magnetic field $B_{a}$ in the $z$ direction will cause the in-plane magnetization to rotate into the direction $\mathbf{e}_{y}$ normal to the wall, so that a moving DW is associated with $\phi\neq 0$ (see section VI). We therefore keep the value of $\phi$ general in the following. Uniform DW solutions of the LLG equation (1), located at $y=Q$, take the general form [14]

[TABLE]

(see Fig. 1b), where the DW width $D$ remains to be determined. The derivation of this solution ignores the nonlocal effect of demagnetizing fields, however we do not expect deviations from this form to be important. We will therefore use this form when computing the demagnetizing energy.

The rest of this paper is structured as follows: In Sec. II we derive a convenient form for the contributions to the demagnetizing energy due to the in-plane and out of plane parts of the magnetization vector, respectively. We then study each of these contributions separately in the following sections III and IV, before applying the results to determine the DW width $D$ in Sec. V, and to the motion of DWs in Sec. VI. Our results are verified by comparison with micromagnetic simulations in Sec. VII, before we conclude in Sec. VIII.

II Magnetostatic energy integrals

The demagnetizing energy due to a magnet with magnetization vector $\mathbf{m}(\mathbf{r})$ is given by

[TABLE]

where the demagnetizing field $\mathbf{H}_{d}$ is determined by Gauss’ law for magnetic fields, $\nabla\cdot\mathbf{B}=\mu_{0}(\nabla\cdot\mathbf{H}_{d}+M_{s}\nabla\cdot\mathbf{m})=0$, as well as Ampere’s law $\nabla\times\mathbf{H}_{d}=\mathbf{J}=0$. The solution of these equations can be given in terms of Green’s functions as

[TABLE]

where $\mathbf{m}^{\prime}=\mathbf{m}(\mathbf{r}^{\prime})$, $\nabla^{\prime}$ denotes differentiation with respect to $\mathbf{r}^{\prime}$, and $d\mathbf{S}^{\prime}$ is the surface normal element at the boundary of the system. In a thin film of thickness $\delta$ much smaller than the relevant magnetic length scales, we can assume that the magnetization is constant in the $z$ direction normal to the film surface. Taking the film to be large in the lateral directions, the only relevant boundaries are the two horizontal surfaces of the film at $z=\pm\frac{\delta}{2}$. Inserting the Green’s function integrals for $\mathbf{H}_{d}$ into the demagnetizing energy, integrating by parts, and using these assumptions, we can show that the energy takes the form

[TABLE]

where the integration now extends only over the two-dimensional area of the film. The in-plane interaction kernel $g_{\delta}$ is given by integrating out the $z$ direction,

[TABLE]

The out of plane interaction kernel, meanwhile, comes from a surface integral over the thin film boundary, so we must instead evaluate the two coordinate vectors at the upper and lower boundaries,

[TABLE]

We now consider the in-plane $E_{d}^{V}$ and the out of plane $E_{d}^{S}$ contributions to the demagnetizing energy separately.

III In-plane energy and the effective demagnetizing constant

The in-plane demagnetizing energy $E_{d}^{V}$ in Eq. (11) requires the divergence of the magnetization in Eq. (5), which is given by

[TABLE]

Inserting into Eq. (11) and scaling the coordinates by $\frac{1}{D}$, the interaction kernel $g_{\delta}$ will transform as $g_{\delta}(r)=Dg_{\frac{\delta}{D}}\left(\frac{r}{D}\right)$. Also defining the small aspect ratio $\sigma=\frac{\delta}{D}$, we find

[TABLE]

where $w,h$ are the linear sizes of the film in the $x,y$ directions, respectively, and $q=\frac{Q}{D}$. This integrand decays exponentially with $y$ and $y^{\prime}$, so we can take $h\to\infty$ and $q\to 0$ without changing the result significantly as long as the DW is located near the center of the film. On the other hand, by symmetry we expect the energy to increase linearly with the width $w$ of the system, so we absorb this and some other constants into the definition by setting $\mathcal{E}_{V}=\frac{E_{d}^{V}}{\mu_{0}M_{s}^{2}w\sin^{2}\phi}$ and working with the reduced energy $\mathcal{E}_{V}$ instead. We now substitute into relative coordinates given by

[TABLE]

which transforms the integration limits to $-\frac{w}{D}..\frac{w}{D}$ for the $u$ integral and $-\frac{w-D|u|}{2D}..\frac{w-D|u|}{2D}$ for the $U$ integral. Since the integrand is independent of $U$, the integration over this variable amounts to a factor $\frac{w}{D}-|u|$.

The hyperbolic functions are most easily transformed to these coordinates by writing them out using their exponential definititions. The resulting transformed integral is given by

[TABLE]

The integral over $V$ can be done by substituting $t=e^{2V}$ and performing a partial fraction decomposition. Using that the integrand is even in $v$ to restrict the limits to $0\ldots\infty$ and integrating by parts in $v$, we obtain

[TABLE]

After differentiation with respect to $v$ under the integral sign and taking $w\to\infty$, the integral over $u$ in the $G_{\sigma}$ function can be computed, giving

[TABLE]

This is further simplified by another differentiation with respect to $v$,

[TABLE]

so after another integration by parts the energy is simplified to

[TABLE]

To make further progress we will need to expand in the small parameter $\sigma$ and integrate term by term. However, a naive expansion of Eq. (23) leads to integrals which diverge at the origin. This is because the interchange of summation and integration is only valid if the integrand is an analytic function of $v$ on the entire contour of integration, but the integrand has a branch cut when $v$ goes from $-i\sigma$ to $i\sigma$, which includes $v=0$. In order to avoid this branch cut, we extend the integration limits back to $-\infty\ldots\infty$ and decompose the logarithm as

[TABLE]

where $I_{\sigma}$ is the integral keeping only the first term inside the square brackets. This isolates the branch cut to the negative imaginary half-plane, so we can deform the integration contour to the contour $C$ going from $-\infty$ to $-r$ with $r>0$ an arbitrarily small number, then around a semicircle of radius $r$ into the positive imaginary half-plane to avoid the origin, then from $r$ to $\infty$ (see Fig. 2). Expanding the logarithm in $\sigma$, we find

[TABLE]

These integrals can be solved by multiplying the integrand with $e^{i\epsilon v/\pi}$ for some $\epsilon>0$ to ensure convergence in the positive imaginary half-plane, then extending the integration contour with a counterclockwise semicircle of radius $R$, which gives a vanishing contribution when $R\to\infty$. Summing over the residues at $v=i\pi k,k\in\mathbb{N}$ and then taking $\epsilon\to 0$, we find the values

[TABLE]

where the last line, valid for $n>2$, uses a known relation between the Dirichlet eta function $\eta(s)=\sum_{k=1}^{\infty}\frac{(-1)^{s+1}}{k^{s}}$, and the Riemann zeta function [15]. Inserting back into $\mathcal{E}_{V}$, the first-order term will cancel with the $\delta D$ term, giving an in-plane reduced energy of

[TABLE]

recalling that the full demagnetizing energy is related to this quantity by $E_{d}^{V}=\mu_{0}M_{s}^{2}w\sin^{2}\phi\mathcal{E}_{V}$.

This energy can be interpreted in terms of an effective demagnetizing constant $N_{y}$ inside the DW. Such a demagnetizing constant would mean that the demagnetizing field is given by

[TABLE]

Inserting into Eq. (6), this leads to a demagnetizing energy given by

[TABLE]

Comparing with the energy we computed in Eq. (29), we see that $N_{y}$ must be chosen as

[TABLE]

This expression should be compared with the demagnetizing constant obtained by taking the DW as a uniformly magnetized ellipsoid, given in Eq. (4), which can be expanded in $\delta/D$ to give

[TABLE]

While this has the same qualitative behavior as our expression (32), quantitatively it is quite different. Our lowest order term is smaller by a factor $\frac{\ln 2}{\pi}\approx 0.22$, which has a direct effect on the motion of DWs (see Sec. VI). Indeed, in Ref. [5], the elliptic approximation was used to estimate the Walker field $B_{W}$ from experimentally measurable quantities, giving $B_{W}\approx 12\,\text{mT}$ for the $0.5\,\text{nm}$ thin film, while micromagnetic simulations of the same system instead gave $B_{W}\approx 2.7\,\text{mT}$ [16], which is reproduced by our analytical computation (see also Sec. VII). Other authors use $\pi D$ in place of $D$ in the expression for $N_{y}^{E}$ [17]. This gets closer to our result, but will still give a different second-order correction.

IV Out of plane energy and the restoring force

Inserting $m_{z}=\tanh(\frac{y-Q}{D})$ into Eq. 12 and scaling the coordinates by $1/D$, the interaction kernel transforms as $f_{\delta}(r)=\frac{1}{D}f_{\sigma}\left(\frac{r}{D}\right)$, giving

[TABLE]

where again $q=\frac{Q}{D}$. By contrast to the in-plane energy, this energy is not only due to the DW, but also contains significant contributions from the dipole charges of the domains themselves. We will therefore have to keep the system size $h$ and the position $q$ general, expecting in particular a quadratic dependence on the position $q$ for a linear restoring force. Changing variables to the relative coordinates of Eqs. (17), taking $w\to\infty$ and integrating over the $U,u$ variables, the reduced energy $\mathcal{E}_{S}=\frac{E_{d}^{S}}{\mu_{0}M_{s}^{2}w}$ takes the form

[TABLE]

where we also used that the integrand is even in $v$ to keep $v$ positive while integrating, simplifying sign issues. The integral over $V$ can now be done using similar techniques as for the in-plane energy. However, the more general limits of integration lead to a more complicated expression. Defining the small aspect ratio $\nu=\frac{D}{h}$, we find

[TABLE]

If the DW is close to the center of the sample, we will have $q\ll\nu^{-1}$, so the denominator of the logarithm can be simplified to $\cosh(\nu^{-1})$. We can then expand in $q$ to give $M(v)=M_{0}(v)+M_{2}(v)q^{2}$ with

[TABLE]

In this way, the energy is decomposed as $\mathcal{E}_{S}=\mathcal{E}_{M}+\mathcal{E}_{q}+\mathcal{E}_{a}$, with

[TABLE]

We first consider the position-dependent $\mathcal{E}_{q}$, which will determine the strength of the restoring force pulling the DW towards the center of the sample. We can note that $M_{2}$ is exponentially small unless $v$ is close to $\frac{\nu^{-1}}{2}$, so we can extend the limits of integration to infinity. Changing variables to $t=v-\frac{1}{2}\nu^{-1}$, we find

[TABLE]

where we expanded the logarithm in powers of $\nu$ and kept the lowest-order term. We will want to compute the other energy terms to the same order in $\nu$, so that we can compare energies accurately. Of these, $\mathcal{E}_{a}$ can be done simply using integration by parts, giving

[TABLE]

where we expanded the result to order $\nu^{2}$.

To compute $\mathcal{E}_{M}$, it helps to understand the $M_{0}(v)$ function. It is plotted in Fig. 3, where we see that it can be approximated well by two different linear functions of $v$ for the first and second half of the interval, respectively. Indeed, writing out the hyperbolic cosine and approximating $2\cosh\nu^{-1}\approx e^{\nu^{-1}}$, we find

[TABLE]

where deviations from the linear behavior, described by the logarithms, are significant only when $v$ is close to $\frac{\nu^{-1}}{2}$. We therefore decompose this energy further as $\mathcal{E}_{M}=\mathcal{E}_{L}+\mathcal{E}_{R}+\mathcal{E}_{d}$, where $\mathcal{E}_{L}$ and $\mathcal{E}_{R}$ capture the linear behavior at the left and right half, respectively, and $\mathcal{E}_{d}$ captures the deviations from linear behavior on both sides. Ignoring the $\coth v$ function in $\mathcal{E}_{d}$ and $\mathcal{E}_{R}$, this means that

[TABLE]

For $\mathcal{E}_{d}$ we use the same substitution and expansion as we did for $\mathcal{E}_{q}$ to find

[TABLE]

$\mathcal{E}_{R}$ can be treated similarly to $\mathcal{E}_{a}$, giving

[TABLE]

Finally, in $\mathcal{E}_{L}$ we can not ignore the $\coth v$ function. Instead we expand in $\sigma$ using the same techniques as we used for the in-plane energy, namely

[TABLE]

where $I_{\sigma}$ replaces the logarithm with $\ln\left(1+i\frac{\sigma}{v}\right)$. Deforming the contour to avoid the origin and expanding in $\sigma$, we obtain

[TABLE]

where the contour $C$ is as in the previous section except that the endpoints are at $\pm\nu^{-1}$ instead of $\infty$. This still gives significant contributions from large $v$ when $n$ is small, so we can not use the standard residue method here. Instead we treat the integral explicitly by dividing the contour into two parts along the real line and one small semicircle where we can substitute $v=re^{i\theta}$ and expand in $r$. This results in the values

[TABLE]

where $\lambda$ is a numerical integration constant given by

[TABLE]

For $n>2$, we can enlarge the integration contour to infinity and use the standard semicircle contour to find

[TABLE]

where $\kappa_{n}$ is the error due to enlarging the contour. This vanishes for $n$ odd due to the integrals on the real line cancelling each other, but for even values it can be approximated as

[TABLE]

where we again ignored the $\coth$ function in the integral. To order $\mathcal{O}(\nu^{2})$, only the $n=4$ value is important, with value $\kappa_{4}=-4\nu^{2}$.

Combining everything, the reduced out of plane energy is given to order $\nu^{2}$ and $\sigma^{4}$ by

[TABLE]

These energy terms can be divided into four different types:

Terms independent of $D$ and $Q$: These are dominated by the $\frac{\delta h}{2}$ term corresponding to the energy $\frac{1}{2}\mu_{0}M_{s}^{2}V$ of two isolated domains magnetized in the $z$ direction, with combined volume $V=\delta hw$. This of course ignores the excluded volume from the DW itself, which is accounted for by the $-\delta D$ term. Higher-order terms give corrections to this, resulting in an effective demagnetizing constant $N_{z}$ which is slightly below $1$. 2. 2.

The logarithmic term: This can be interpreted as an interaction between the two thin film domains. They will attract each other due to the oppositely directed magnetizations, with a force $-\dfrac{\partial{}E}{\partial{D}}\sim-\frac{\delta^{2}}{D}$. 3. 3.

Terms depending on $D$, but independent of $Q$: These correspond to the demagnetizing energy of the DW itself. 4. 4.

The term proportional to $Q^{2}$: This represents a harmonic restoring force pulling the DW back to the center at $Q=0$, so that zero net magnetization is preferred.

In refs. [18, 19], the energy due to the restoring force is assumed to take the form

[TABLE]

where $\left\langle-\right\rangle$ denotes averaging over space, and $\left\langle H_{d}\right\rangle$ is taken as $-M_{s}\mathcal{N}\left\langle m_{z}\right\rangle$ for some effective demagnetizing constant $\mathcal{N}$ describing the entire domain structure. Inserting for $m_{z}$, this gives

[TABLE]

Comparing with our result $E_{Q}=4\mu_{0}M_{s}^{2}\frac{w\delta^{2}}{\pi h^{2}}Q^{2}$, we see that the effective demagnetizing constant must be chosen as $\mathcal{N}=\frac{2\delta}{\pi h}$.

V Steady-state domain wall width

The DW width $D$ is a dynamical variable evolving with time. In equilibrium or steady-state motion, the steady-state DW width is the one that minimizes the energy at a given value of $\phi$. In addition to the demagnetizing energies we computed above, the Landau-Lifshitz energy includes contributions from the exchange energy and the anisotropy energy. Using the form of Eq. (5) for the DW, these energies are readily computed as

[TABLE]

where we added a constant energy density to $E_{a}$ to keep it finite when $h\to\infty$. To the lowest order in $\delta$, the only contributing term of the demagnetizing energy is $\delta D$ from the out of plane energy, giving a minimizing equation

[TABLE]

with solution $D_{0}$ given by

[TABLE]

Considering higher orders in $\delta$, it is convenient to define the exchange length $D_{\text{ex}}=\sqrt{\frac{2A_{\text{ex}}}{\mu_{0}M_{s}^{2}}}$. In terms of $D_{0}$ and $D_{\text{ex}}$, the minimizing equation is given to order $\delta^{3}$ by

[TABLE]

Instead of solving this cubic equation directly, we treat it perturbatively in orders of $\delta$, by expanding the solution $D_{\text{eq}}$ as

[TABLE]

and solving for each order in $\delta$ separately. To order $n=0$, this gives the value in Eq. (62), while for higher orders we find

[TABLE]

This can be compared with the result commonly obtained by using demagnetizing constants (setting $N_{x}^{E}=0$ and $N_{z}^{E}=1-N_{y}^{E}$ in this geometry) [12],

[TABLE]

which depends on the film thickness $\delta$ through $N_{y}^{E}=\frac{\delta}{\delta+D_{N}}$. Expanding to first order in $\delta$, we obtain

[TABLE]

Here the dependence on the angle $\phi$ is of order $\mathcal{O}(\delta)$, by contrast with our result which only depends on $\phi$ in the second order term $D_{2}$. The difference is that the derivation of $D_{N}$ ignores the fact that $N_{y}$ depends on $D$ when minimizing the energy. As we can see from Eq. (31), the in-plane demagnetizing energy due to $N_{y}$ is proportional to $DN_{y}$, which should be differentiated with respect to $D$ to find the minimum. However, the lowest order term of $N_{y}$ is generally proportional to $\frac{\delta}{D}$, so the lowest order term of $\dfrac{\partial{}(DN_{y})}{\partial{D}}$ cancels out, leaving only a term of order $\mathcal{O}(\delta^{2})$. Our expression for $D_{\text{eq}}$ takes proper account of the dependence of the energy on the DW width, giving the correct dependence on $\phi$ as well as more precise constants.

VI Domain wall dynamics

The demagnetizing energies we computed here can be used to derive an accurate 1D model for the motion of a uniform DW, as originally derived by Slonczewski [9]. Following Refs. [11, 20], we employ the Lagrangian formulation of the LLG equation. The conservative Landau-Lifshitz equation can be posed in a Lagrangian form for the angle $\theta$ between $\mathbf{m}$ and the $z$ axis, and the angle $\phi$ describing the in-plane component. The Lagrangian is given by $L=\frac{M_{s}}{\gamma}\int\dot{\phi}\cos\theta d^{3}\mathbf{r}-E$, where the energy $E$ also includes the Zeeman energy due to a constant applied field $B_{a}$ in the $z$ direction, given by $-M_{s}B_{a}\int m_{z}d^{3}\mathbf{r}$. Gilbert dissipation can be included using a Rayleigh dissipation functional given by $F=\frac{\alpha M_{s}}{2\gamma}\int|\dot{\mathbf{m}}|^{2}d^{3}\mathbf{r}$. Inserting the ansatz (5) into these functionals and integrating over space, we find the Lagrangian and dissipation functional governing the three variables $s_{i}=\mathinner{\left\{Q,\phi,D\right\}}$ describing the DW. These variables obey the dissipative Euler-Lagrange equation for each variable,

[TABLE]

which results in the equations of motion given by

[TABLE]

Here $B_{R}$ is the effective field corresponding to the restoring force, and $N_{y}$ and $\dfrac{\partial{}E}{\partial{D}}$ are given in Eqs. (32) and (63), respectively. Eq. (72) describes how the DW width $D$ relaxes towards the steady-state value $D_{\text{eq}}$. We can estimate how fast this relaxation is by linearizing around the steady state. To zeroth order in $\delta$, this gives an exponential approach with relaxation time

[TABLE]

with higher-order corrections derivable. This timescale is shorter than the timescale $\tau_{V}=\frac{D}{V_{W}}$ of fast DW motion (see below) by a factor $\frac{\tau_{D}}{\tau_{V}}\propto\alpha N_{y}\ll 1$. It is therefore common to ignore the dynamics of $D$ and set $D=D_{\text{eq}}(\phi)$ at each point in time [11].

Walker-like steady-state solutions are found by setting $\dot{\phi}=0$ in Eqs. (70–71). The resulting equations are solvable only if

[TABLE]

giving an expression for the Walker breakdown field $B_{W}$. Below this field, the steady-state solution gives a DW velocity of

[TABLE]

In particular, the velocity at Walker breakdown is given by

[TABLE]

Note that these quantities depend on the geometry through both the steady-state DW width $D_{\text{eq}}$ (64), and the effective demagnetizing constant $N_{y}$ (32).

VII Numerical verification

To verify our analytic computations, we performed micromagnetic simulations using the MuMax3 software package [21]. An initial Bloch-type DW configuration was generated by setting the magnetization to point along $+z$ for $0\leq y<h/2$ and $-z$ for $h/2<y\leq h$, with a small region pointing along $+x$ at the boundary between the domains (see Fig. 1a). After relaxing the initial configuration in zero field to an energy minimum, we applied external magnetic field and calculated the subsequent DW motion by numerically integrating the LLG equation (1). In all simulations we used micromagnetic parameters previously determined from experiments on Pt/Co/Pt films with perpendicular anisotropy [5], with exchange stiffness $A_{\text{ex}}=1.4\times 10^{-11}\,\text{J/m}$, saturation magnetization $M_{s}=9.1\times 10^{5}\,\text{A}/\text{m}$, uniaxial anisotropy constant $K=8.4\times 10^{5}\,\text{J}/\text{m}^{3}$, and dissipation constant $\alpha=0.27$.

The variation of the Walker breakdown field $B_{W}$ with the film thickness $\delta$ was computed on a rectangular grid of $64\times 128\times 1$ cells with in-plane cell size $\Delta x=\Delta y=1\,\text{nm}$ set well below $D_{0}\approx 6.6\,\text{nm}$, and the cell thickness $\Delta z$ ranging from $0.5-4.0\,\text{nm}$. Since there is always only one cell in the $z$ direction, this enforces the assumption that $\mathbf{m}$ is independent of $z$. We used periodic boundaries in the $x$ direction, and to remove the effect of the restoring force $B_{R}$ we continually updated the position of the simulation window along $y$ to keep the DW centered inside the window. For each film thickness $\delta=\Delta z$, simulations were carried out in a $1\,\text{mT}$ range centered on the Walker field given by Eq. (74). The resulting average DW velocities $V_{DW}$ for each simulation are shown in Fig. 4 (filled symbols). The DW velocity increases with the applied field $B_{a}$, until Walker breakdown $B_{a}=B_{W}$ when the average velocity drops abruptly due to precession. The peaks in $V_{DW}$ (open symbols in Fig. 4) are therefore numerical estimates of the Walker field for each film thickness.

The DW width is measured numerically by fitting a line $ay+b$ to values of $\operatorname{atanh}(m_{z})$ close to the DW position, which is defined as the location where $m_{z}$ crosses [math], and setting $D=|a|^{-1}$. The resulting values are compared with the analytical result (64) in Fig. 5. For increasing film thickness the DW width diminishes, due to the negative first-order correction $D_{1}$ [see Eq. (65)]. This is mainly due to the attractive logarithmic interaction between the antiparallel domains as evident in Eq. (56), so that the distance $D$ between the domains is reduced as the strength of interaction increases.

The numerical values for the Walker breakdown field $B_{W}$ and velocity $V_{W}$ are shown as a function of $\delta/D$ Fig. 6, and compared with analytical predictions given by Eqs. (74) and (76) to first, second, and third order. The numerical and analytical values show good agreement with each other, highlighting the impact of higher-order corrections to the Walker field, particularly for larger values of $\delta/D$.

Analytical predictions for the restoring force arising from demagnetizing effects were also validated by comparison with simulations on a grid of $128\times n_{y}\times 1$ cells with $64\leq n_{y}\leq 512$, and cell size $\Delta x=\Delta y=1\,\text{nm}$ and $\Delta z=4\,\text{nm}$. The Bloch wall initial condition was relaxed in an applied magnetic field of varying strength, displacing the DW from the center of the sample. We then computed the demagnetizing energy of the relaxed configuration, and subtracted the $Q=0$ value to isolate the quadratic dependence given by Eq. (43). The result is given in Fig. 7, and shows a good agreement with the analytical prediction, up to some deviation from the quadratic behavior at large $Q^{2}$.

VIII Conclusion

The magnetostatic energy of domain structures is a challenging mathematical problem in the general case. Here we have made progress on the idealized case of a uniform, infinitely long DW in a thin PMA film by deriving analytic expressions for the energy. This allows accurate predictions for important properties of the DW, such as the DW width (Sec. V), DW dynamics (Sec. VI), and restoring force (Sec. IV). Micromagnetic simulations were employed to verify our results.

The effect of the in-plane magnetization of the DW can be understood by using an effective demagnetizing constant $N_{y}$, whose precise value differs from the commonly used elliptic approximation. The out of plane energy, on the other hand, consists of separate contributions from the DW itself, the two domains, and the interaction between the domains, which would be difficult to disentangle from each other without the principled approach employed here. For example, the DW width is affected by the attractive interaction between the domains at either side, which the formalism of demagnetization constants fails to include. Using an explicit expression for the in-plane energy when finding the equilibrium DW width also avoids a subtle mistake due to the variation of the effective demagnetization constant $N_{y}$ with the DW width $D$.

For simplicity we considered the case of an infinitely extended thin film without vertical boundaries. This idealized case does not include the effect of disorder, which inevitably distort the shape of the DW in real thin films. The more realistic problem of a uniform DW in a nanostrip of width $w=\mathcal{O}(D)$ introduces further difficulties. Additional boundary integrals will need to be included in Eqs. (11–12). More severely, the technique of deforming the integration contour to integrate term by term in Eqs. (25, 52) relied on simplifying the interaction kernels $g_{\sigma}$ and $f_{\sigma}$ into a simple logarithmic form [see Eqs. (19–22) and (35)], which was obtained by taking the film width $w$ to infinity. It is unclear how to perform a similar expansion without this simplification. Some mathematical difficulties therefore stand in the way of extending these results to the case of nanostrips.

In order to consider thicker films with width $\delta\sim D$, it is necessary to allow for variations in the vertical direction $\mathbf{e}_{z}$ as well. In this case, the magnetization vector can deflect away from the vertical close to the surface, in order to reduce the energy penalty from the out of plane demagnetizing field. This can lead to the formation of complicated structures such as Bloch lines at the surface [14, 22]. Our analysis is therefore restricted to $\delta\ll D$, which avoids these complications.

Another avenue for generalization is to include effects such as the Dzyaloshinskii-Moriya interaction (DMI) and spin-transfer torque. The DMI is a local energy term which can be straightforwardly included in our analysis [23]. Spin-transfer torque, on the other hand, is a dynamical forcing mechanism and not an energy. It can however be included in the Lagrangian framework as a dynamical term in the Lagrangian [20].

In all, our analytical computations provide a much better understanding of the effect of long-range demagnetizing fields on the properties and motion of DWs in thin films. The methods we have used are quite general, and can maybe be employed to understand other similar systems, such as DWs in systems with in-plane magnetic anisotropy, and provide a solid foundation on which a principled understanding of more complicated DW behavior can be built.

Acknowledgements.

This work has been supported by the Academy of Finland through an Academy Research Fellowship (L.L.; project no. 268302).