Energy barriers in three-dimensional micromagnetic models and the physics of thermo-viscous magnetization in multidomain particles

TL;DR

This paper presents a micromagnetic model and a new computational method to analyze viscous remanent magnetization in particle ensembles, explaining stability and acquisition dynamics in paleomagnetic studies.

Contribution

It introduces a fast relaxation algorithm combining nudged elastic band and action minimization to compute energy barriers in micromagnetic particles.

Findings

Identified 60 local energy minima and calculated all energy barriers between them.

Estimated energy barriers in weak external fields for VRM transition matrices.

Demonstrated that VRM acquisition can overshoot and occur faster than decay in PSD particles.

Abstract

A first principle micromagnetic and statistical calculation of viscous remanent magnetization (VRM) in an ensemble of cubic magnetite pseudo single-domain particles is presented. This is achieved by developing a fast relaxation algorithm for finding optimal transition paths between micromagnetic local energy minima. It combines a nudged elastic band technique with action minimization. Initial paths are obtained by repetitive minimizations of modified energy functions. For a cubic pseudo-single domain particle, 60 different local energy minima are identified and all optimal energy barriers between them are numerically calculated for zero external field. The results allow to estimate also the energy barriers in in weak external fields which are necessary to construct the time dependent transition matrices which describe the continuous homogeneous Markov processes of VRM acquisition and decay. By spherical averaging the remanence acquisition in an isotropic PSD ensemble was calculated over all time scales. The modelled particle ensemble shows a physically meaningful overshooting during VRM acquisition. The results also explain why VRM acquisition in PSD particles can occur much faster than VRM decay and therefore can explain for findings of extremely stable VRM in some paleomagnetic studies.

Equations

$$S(x(t)):=\underset{0}{\overset{t_{end}}{\int }}⟨\dot{x}+\nabla E,\dot{x}+\nabla E,⟩dt.$$

$$S_{\min }(p)=S(x_p(t))=\underset{0}{\overset{L(p)}{\int }}⟨\frac{dx}{ds}v+\nabla E,\frac{dx}{ds}v+\nabla E,⟩\frac{ds}{v},$$

$$\frac{d}{dv}\left(\frac{1}{v}⟨\frac{dx}{ds}v+\nabla E,\frac{dx}{ds}v+\nabla E,⟩\right)=0.$$

$$v=\Vert \dot{x}\Vert =\Vert \nabla E\Vert {\Vert \frac{dx}{ds}\Vert }^{-1}=\Vert \nabla E\Vert .$$

$$S_{\min }(p)=2\underset{0}{\overset{L(p)}{\int }}\Vert \nabla E\Vert \Vert \frac{dx}{ds}\Vert +⟨\frac{dx}{ds},\nabla E⟩ds=2\underset{0}{\overset{L(p)}{\int }}\Vert \nabla E\Vert +⟨\frac{dx}{ds},\nabla E⟩ds.$$

$$\underset{0}{\overset{L(p)}{\int }}⟨\frac{dx}{ds},\nabla E⟩ds=\underset{E(x_0)}{\overset{E(x_1)}{\int }}dE=E(x_1)-E(x_0)=:\Delta E.$$

$$S_{\min }(p)=2\Delta E+2\underset{0}{\overset{L(p)}{\int }}\Vert \nabla E\Vert ds.$$

$$S_{\min }(p)=2\Delta E+2\underset{0}{\overset{1}{\int }}\Vert \nabla E\Vert \Vert \frac{dx}{dw}\Vert dw.$$

$$\frac{d^{2}x}{d{s}^2}=\nabla \log \Vert \nabla E\Vert .$$

$$\dot{x}=\pm \nabla E$$

$$S_{\min }(p)\ge 2\Delta E+2\underset{0}{\overset{1}{\int }}\left|⟨\nabla E,\frac{dx}{dw}⟩\right|dw.$$

$$S_{\min }(p)\ge 2\Delta E+2\sum _{k=0}^K\left|E(x(w_{k+1}))-E(x(w_k))\right|.$$

$$N_{\min }+N_{\max }-N_{\mathrm{saddle}}={\chi }_{\mathrm{Euler}}.$$

$$\sum _{k=0}^N(-1{)}^k{N}_k={\chi }_{\mathrm{Euler}}.$$

$$N_{\min }+N_{\max }+N_{\mathrm{even}}-N_{\mathrm{odd}}=2^N,$$

$$\delta S_{\min }=2\underset{0}{\overset{L}{\int }}\left(\nabla \Vert \nabla E\Vert -\Vert \nabla E\Vert \frac{d^{2}x}{d{s}^2}\right)\delta xds.$$

$${\delta }^{\ast }x=\frac{\nabla \Vert \nabla E\Vert -\Vert \nabla E\Vert \frac{d^{2}x}{d{s}^2}}{{\Vert \nabla \Vert \nabla E\Vert -\Vert \nabla E\Vert \frac{d^{2}x}{d{s}^2}\Vert }_2}.$$

$$x_n:=x_{n-1}-\alpha {\delta }^{\ast }{x}_{n-1},$$

$$\frac{\partial x}{\partial u}(u,s)=-\nabla E+\frac{⟨\frac{\partial x}{\partial s},\nabla E⟩}{\Vert \frac{\partial x}{\partial s}{\Vert }^2}\frac{\partial x}{\partial s},$$

$$x(0,s)=x_0(s),x(u,0)=x_0(0),x(u,L)=x_0(L).$$

$$d(m_1,m_2):={\left[\frac{1}{V}\underset{V}{\int }{\arccos }^2(m_1(r)\cdot m_2(r))dV\right]}^{1/2}.$$

$$v(m_1,m_2):=-\nabla d(.,m_2)(m_1).$$

$$m^{k+1}(s):=m^k(s)-{\alpha }_k\left[\nabla E(m^k(s))-\left(\nabla E(m^k(s))\cdot t^k(s)\right)t^k(s)\right],$$

$$E_{\Delta }^{\ast }(m)=E(m)+\mu [d(m,m_A)-\Delta {]}^2+\frac{\beta }{d(m,m_B)+\epsilon }.$$

$$M(\Delta t)\equiv \mathbb{P}\left[S(t)=S_j\wedge S(t+\Delta t)=S_i\right]=\exp \left[\mu \Delta t\right].$$

$${\mu }_{ij}$$

$${\mu }_{ii}$$

$$\rho (t)=\exp \left[\mu t\right]{\rho }_i^0$$

$$m(t)=\sum _i{m}_i{\rho }_i(t).$$

$$E_b^{ij}(H)=E_b^{ij}+(m_j-m_{ij}^{\max })H,$$

$${\mu }_{ij}(H)$$

Full text

1 Introduction

1.1 Aim and outline of the article

The geomagnetic field has been perpetually recorded by magnetic remanence carriers in newly formed rocks throughout the Earth’s history. Therefore, crustal rocks form a paleomagnetic archive which is accessible through rock magnetic measurements. Yet, their interpretation requires a thorough understanding of the physical processes occuring during remanence acquisition. Thermo-viscous magnetization of natural pseudo-single or multidomain particles is the most abundant remanence in paleomagnetism, although more reliable single domain (SD) remanence carriers are preferred, and theoretical interpretation is based on paradigms developed from SD theory. The main aim of this article is to propose new computational and conceptional methods to obtain a physical understanding of the remanence acquisition in multidomain particles. The second section will introduce a new technique to determine energy barriers in micromagnetic models of inhomogenously magnetized particles. Although this is technically challenging, it is an essential prerequisite for a quantitative study of thermo-viscous magnetization. In the third section we use the new computational method to calculate all energy barriers for a three-dimensional model of metastable flower and vortex-states in a cubic magnetite particle. The fourth section introduces the theoretical background for a statistical analysis of viscous remanence acquisition and decay in weak fields, based on the computed energy barriers between the metastable states. In the fifth section this theory is applied to the energy barriers calculated in section 3, and the physical meaning of the model results obtained is discussed.

1.2 Micromagnetic modeling

Micromagnetic modeling is now a standard technique to determine stable and metastable magnetization states in small ferro- and ferrimagnetic particles. It is used to calculate and analyze magnetization structures in natural and synthetic magnetic nanoparticles. This size range is of special interest in rock magnetism where the reconstruction of the Earth magnetic field depends critically on the reliability of remanence carriers in natural rocks. The grain size distribution of these remanence carriers rarely is confined to the relatively small SD-size range. Accordingly, detailed knowledge of the physical mechanism of magnetization change in larger nanoparticles is needed to assess and evaluate the magnetic measurement results from natural materials. Because of its importance for understanding remanence acquisition in natural rocks, rock magnetic studies were among the first to apply numerical micromagnetic models. The first approach to estimate barriers between single-domain and two-domain states used a one-dimensional model of magnetization change [5]. When three-dimensional models were developed to understand inhomogenous remanence states [18], it was immediately a main interest to obtain energy barriers to model the acquisition of thermoremanence [6, 17, 19, 13]. Knowing the energy barriers between different magnetization states within a single particle also leads to a quantitative prediction of magnetic viscosity and magnetic stability of remanence information, even over geological time scales. An important result of early micromagnetic calculations was that beyond the regime where exchange forces dominate, i.e. beyond length scales of several exchange lengths $\sqrt{A/K_{d}}$, there exist a multitude of local energy minima corresponding to meta-stable magnetization structures [18, 8, 16]. In the context of thermo-viscous remanence, the most important property of meta-stable magnetization structures $m$ is their residence time $\tau(m)$. It denotes the expectation value of the time during which the system remains in state $m$, if it initially is in this state at time $t=0$. The residence time $\tau(m)$ directly depends on the transition probabilities $p(m,m^{\prime})$ between $m$ and all other LEM $m^{\prime}$, which in turn depend upon the possible transition pathways. To determine the transition probability $p(m,m^{\prime})$ in very good approximation, it is sufficient to find the most likely transition path between $m$ and $m^{\prime}$, which is the path with the lowest energy barrier. This path runs across the saddle-point with lowest energy of all which connect $m$ and $m^{\prime}$. Therefore, the problem of finding the transition probabilities $p(m,m^{\prime})$ is closely related to finding saddle-points in high-dimensional micromagnetic energy-landscapes.

1.3 Statistical theory

Given these transition probabilities, the geologically important mechanism of thermoremanence acquisition can be described as a stochastic process of magnetization change in a temperature dependent energy landscape. Its transition matrix is related to $p(m,m^{\prime})$, determined by the energy barriers between the possible states [7].

2 Action and path integrals

2.1 Micromagnetic modeling

Berkov [2, 3, 1] developed a numerical method to evaluate the distribution of energy barriers between metastable states in many-particle systems which determines the optimal path between the two given metastable states by minimizing the action in the Onsager Machlup functional [15] for the transition probability. This method essentially performs a local saddle-point search in an high-dimensional energy landscape. Mathematically similar problems exist in several disciplines of physics and chemistry. In the last years, several new methods to locate saddle-points have been developed in these fields [10, 9, 14]. Based on such algorithms, an improved elastic band technique for micromagnetics was presented by [4].

The main problem in energy-barrier computation is that micromagnetic structures $m$ are described by many variables and accordingly energy $E$ is a function of $m$. Minimizing $E(m)$ requires sophisticated algorithms, but for energy-barrier calculations it is even necessary to determine saddle points in this high-dimensional energy landscape.

Several approaches are available for this task, but because saddle-point calculation is equivalent to minimizing $(\nabla E(m))^{2}$ all rapidly converging methods require second derivatives of $E$. This however is rather to be avoided if the calculation should be performed effectively.

The present study develops a combination of several of the above cited techniques to efficiently calculate energy barriers in micromagnetic models.

2.2 Action minimization

Berkov [2] introduced action minimization as a tool for finding optimal transition paths in thermally driven micromagnetic systems. He discretized the time dependent action of the the magnetic particle system and used a numerical quadrature representation for direct numerical minimization. This rigorous approach is complicated by its explicit dependence upon transition time. However, transition paths turn out to be geodesics of the energy surface in the limit of infinite transition time, where energy barriers are lowest.

Dittrich et al. [4] make use of this fact by directly searching for geodesic paths using a modification of the nudged elastic band (NEB) algorithm of Henkelman et al.[10, 9].

A problem of this algorithm is that it involves the numerical solution of a large system of ordinary differential equations. Moreover, there is a tendency of the NEB algorithm to produce spurious up-down-up movements along the gradient (kinks) which cannot be completely removed in all cases.

Here, we combine both approaches by designing a path relaxation algorithm similar to NEB, but constraint to decrease the action at each step. The algorithm performs a fast gradient-like relaxation from an initial path towards the optimal transition path. It detects and avoids the development of kinks, and does not involve numerical solutions of differential equations.

The important problem of finding an initial path which is likely to lie in the basin of attraction of the optimal transition path under the proposed relaxation scheme is also investigated.

2.3 Geometric action

Here we define the geometric action for a path $p$ as the minimal action for any transition along this path.

In a general mechanic system the action of a transition process $x(t)$ from state $x(0)=x_{0}$ into $x(t_{end})=x_{1}$ is defined by

[TABLE]

The probability that this transition process occurs depends monotonously on $\exp(-S(x(t)))$. In the next section we will be looking for the optimal transition path in the energy landscape determined by $E$. The quality of any given path $p$ is defined as its geometric action: the action of the transition process $x_{p}(t)$ along $p$ which minimizes (1).

We start with a canonical parametrisation of $p$ by arc length $s$ and try to find a reparametrisation $s(t)$ which minimizes $S$ along $p$. For this optimal transition process we then have

[TABLE]

where $L(p)$ is the arc length and $v(s)~{}=~{}\frac{ds}{dt}(s)$ is the local velocity of the optimal transition at arc length $s$. Finding $s(t)$ thus is reduced to the variational problem of finding the function $v(s)$ which minimizes (2). The corresponding Euler-Lagrange equation is

[TABLE]

A short calculation confirms that it’s solution is

[TABLE]

The last equality uses the fact that for the arc length parametrisation $\|{dx}/{ds}\|~{}=~{}1$. Inserting this result into (2) yields

[TABLE]

This integral can be simplified further by noting that

[TABLE]

Accordingly, one obtains the geometric action of $p$ as

[TABLE]

2.4 Finding the optimal path by variation of the geometric action

It is possible to find the Euler-Lagrange equations for the optimal path by variation of the geometric action $S_{\rm min}(p)$ with respect to $x$.

To this end we reparametrize (7) by $w(s)=s/L(p)$ and obtain

[TABLE]

The Euler-Lagrange equation of the variational problem $\delta S_{\rm min}(p)~{}=~{}0$ after some simplification has the form

[TABLE]

The details of the calculation are given in the appendix.

2.5 The optimal path is a geodesic

A path along an energy surface which fulfills

[TABLE]

is a geodesic. In the one-dimensional case (4) directly implies that the optimal transition path is a geodesic. In the multidimensional case this not simply follows from (4) which is valid for any geometric transition path. Yet, by applying the Cauchy inequality to (8) one obtains for the optimal path $p$

[TABLE]

The integration interval $[0,1]$ can be divided into finitely many parts $[w_{k},w_{k+1}]$ with alternating constant sign of $\langle\nabla E,\,dx/dw\rangle$. Accordingly, $\nabla E(x(w_{k}))~{}=~{}0$ and

[TABLE]

Here the right hand side is a lower limit of $S_{\rm min}$. Accordingly, a geodesic which fulfills (10) achieves equality in (11). It therefore coincides with the least action path between the prescribed endpoints. Equality in (12) means that the least action depends only upon the energies at the traversed critical points.

2.6 Morse theory

Topologically different critical points in multi-dimensions are distinguished by their Morse index, which is defined as the dimension $n_{-}$ of the sub-manifold on which the Hessian is negative definite. Intuitively, the Morse index of the highest saddle point along the optimal transition path should not be too large. This is because , apart from singular cases, the action is minimized along an only one-dimensional manifold, i.e. an isolated path parallel to the gradient which connects initial and final minima. Therefore, if at the highest saddle point the Morse index is $n_{-}>1$, the other $n_{-}-1$ descending directions should lead into different LEM states. The choice for these LEM states should not be too large whenever the initial and final minima are close to the global energy minimum.

On the other hand, Morse theory [12] implies that the total number of saddle-points will be huge in realistic micromagnetic calculations. For a two-dimensional surface with Euler characteristic $\chi_{\rm Euler}$, the numbers $N_{\rm min}$ of minima, $N_{\rm max}$ of maxima, and $N_{\rm saddle}$ of saddle points are connected by the relation

[TABLE]

For a sphere is $\chi_{\rm Euler}=2$. The generalization of (13) to a finite-dimensional compact manifold is the Morse relation

[TABLE]

Here $N_{k}$ is the number of critical points with Morse index $k$, i.e. where the negative definite sub-manifold has dimension $n_{-}=k$. In case of the micromagnetic energy depending on $N$ magnetization directions the manifold is a direct product of $N$ two-dimensional spheres, therefore $\chi_{\rm Euler}=2^{N}$. On this $2N$-dimensional manifold we get from (14)

[TABLE]

where $N_{{\rm even}}$ and $N_{{\rm odd}}$ are the number of true saddle points with even or odd Morse index, respectively.

If this manifold describes an ensemble of interacting SD grains each grain has at least two critical points, minima or maxima, in zero external field. If interaction is weak, the total energy of the system inherits almost all these minima and maxima as saddle-points. Thus, it is understandable that the total number of critical points exceeds even a huge figure like $2^{N}$.

The situation is different if the $2N$-dimensional manifold describes an exchange coupled grain with inhomogeneous magnetization structure, like in most micromagnetic applications. In this case, it turns out that only a limited number of minima and maxima, like flower or vortex states, exist. According to (15), there must appear an enormous number of true saddle points with even Morse index – i.e. not one-dimensional lines.

Thus, inevitably there exists a large number of paths, with different action $S$, which connect the local minima and complicate the search of the lowest action path. Another circumstance causing difficulties to find path lines by numerical means, is that the field lines passing through the saddle points, form sub-manifolds of zero measure. In other words, the saddle points are unstable in the sense, that almost all field lines in their vicinity are deflected aside (except for those which go directly through the saddle point).

2.7 Euler-Lagrange path relaxation

To numerically calculate the optimal geometric transition path one would like to start from an arbitrarily chosen initial path $p_{0}$ which then is iteratively improved by some updating scheme. An heuristic procedure is known as ’nudged elastic band’ technique. Here we derive an optimal updating scheme which corresponds to a gradient minimization of minimal path action. To this end we follow the standard derivation of the Euler-Lagrange equation to prove that the change $\delta S_{\rm min}$ in action due to a variation $\delta x$ in the path is given by

[TABLE]

The direction of maximal increase in $S_{\rm min}$ now is a function $\delta^{*}x\in L_{2}([0,L])$ with $\|\delta^{*}x\|_{2}~{}=~{}1$ which maximizes (16). By Cauchy’s inequality

[TABLE]

The optimal Euler-Lagrange relaxation towards a minimal action path accordingly is given by

[TABLE]

for sufficiently small positive step size $\alpha$. Two points are noteworthy. First, the updating depends only on $\|\nabla E\|$ and not explicitly on $E$, and second, the term $d^{2}x/ds^{2}$ leads to reduction of local curvature in regions where $\|\nabla E\|>0$. This prevents the updating scheme from producing kinks. Only at points with $\|\nabla E\|=0$ the final solution may not be differentiable.

2.8 Relation between action and path relaxation

There exists a heuristic scheme to find a path between two minima across a saddle-point in multidimensional energy landscapes, which makes use of the fact that there is such a path which always runs along the energy gradient.

Starting from an arbitrary initial path $x_{0}(s)$ it is thus attempted to adjust the path in the direction along the negative energy gradient, but perpendicular to the local tangent. This relaxation scheme can be analytically described as the solution of a boundary value problem for a system of nonlinear partial differential equations

[TABLE]

with the boundary conditions

[TABLE]

3 A modified relaxation method to determine transition paths from micromagnetic models

3.1 Definitions

A micromagnetic structure $m$ is determined by $K$ magnetization vectors on a spacial grid over the particle. Each magnetization vector is a unit vector determined by two polar angles $\theta,\phi$. The distance $d(m_{1},m_{2})$ between two magnetization structures is defined by

[TABLE]

The local direction vector from $m_{1}$ to $m_{2}$ is the gradient

[TABLE]

If two magnetization structures $m_{0}$ and $m_{1}$ contain no opposite directions, i.e. for all $r\in V$ we have $m_{0}(r)\neq-m_{1}(r)$, then it is possible to linearly interpolate between $m_{0}$ and $m_{1}$ by defining $m_{t}(m_{0},m_{1})(r)$ as the intermediate vector on the smaller great circle segment connecting $m_{0}(r)$ and $m_{1}(r)$ which has angular distance $t\arccos(m_{0}(r)\cdot m_{1}(r))$ from $m_{0}(r)$.

By minimization of $E(m)$ using gradient information $\nabla E(m)$ an initial minimum $m_{A}$ and a final minimum $m_{B}$ are found.

3.2 Outline of the relaxation procedure

In summary, the results of the previous section show that the transition probability between $m_{A}$ and $m_{B}$ is in very good approximation determined by the minimal energy barrier between them. This barrier is achieved along some optimal geometrical least action path $m(s)$ with $m(0)~{}=~{}m_{A}$ and $m(1)~{}=~{}m_{B}$ which also represents the most likely transition path. The state of maximal energy along this path is a saddle point of the total energy function $E(m)$, and the least action path is everywhere parallel to $\nabla E(m)$. To find this optimal transition path, we propose a relaxation method which combines the advantages of the NEB technique of [4] with the additional action minimization of [3].

Two techniques are required for finding the optimal path by means of iterative relaxation:

An updating scheme which determines an improved transition path $m^{k+1}(s)$ from a previous path $m^{k}(s)$. 2. 2.

A method of finding an initial transition path $m^{0}(s)$ within the basin of attraction of the optimal path.

The here proposed updating scheme starts from an initial path $m^{0}(s)$ which is determined by interpolating $N$ intermediate states $m^{0}(s_{j})$ which correspond to the magnetization structures at the $s$-coordinates $0=s_{1}<s_{2}<\ldots<s_{N}=1$ for $j=1,\ldots,N$. For the interpolation to be well-defined, $N$ is required to be large enough to ensure that neighboring structures $m^{k}(s_{j})$ and $m^{k}(s_{j+1})$ never contain opposite magnetization vectors.

Similar to the NEB method, in step $k$ the path $m^{k}(s)$ is changed according to

[TABLE]

where $t^{k}(s)$ is the tangent vector to the path $m^{k}(s)$ at $s$ and $\alpha_{k}>0$ is a real number. This updating scheme moves the path downward along the part of the energy gradient which is perpendicular to path itself. This algorithm converges to a path which is almost everywhere parallel to $\nabla E$. However, (23) describes not a true gradient descent for the action, and the final path may not achieve minimal action due to the formation of kinks during the minimization [10, 9, 4].

The here proposed method differs from previous NEB techniques in two details: First, $\alpha_{k}$ is chosen such that $S(m^{k+1})<S(m^{k})$. This ensures that the action decreases in each step. Second, $\alpha_{k}$ is dynamically adapted to achieve rapid convergence. The following procedure to choose $\alpha_{k}$ fulfills both aims.

Adaption of $\alpha_{k}$

Starting with the initial value $\alpha_{0}=1$, it is evaluated after each step whether the action of the updated path is decreased, i.e. $S(m^{k+1}(s))<S(m^{k}(s))$. While this is true, the new value $\alpha_{k+1}=\alpha_{k}$ is kept constant, but only for at most five steps. In this phase (23) performs a quasi-gradient descent ’creeping’ towards the optimal path.

Afterwards, if $S$ still decreases, we set $\alpha_{k+1}=2\,\alpha_{k}$ in each following step until some $S(m^{k+1}(s))>S(m^{k}(s))$. This phase can be interpreted as an ’accelerated steepest descent’.

If at any time $S(m^{k+1}(s))>S(m^{k}(s))$, the path $m^{k+1}(s)$ is rejected and (23) is evaluated for a new $\alpha_{k+1}$-value of $\alpha_{k+1}=1/4\,\alpha_{k}$. This behavior avoids ’overshooting’ of the gradient descent and

All our tests show that this iterative adaption of $\alpha_{k}$ leads to a much faster convergence than choosing any fix value $\alpha_{k}=\alpha$.

By comparing the achieved action $S(m^{k}(s))$ to $S_{\rm min}$ from (12) during the relaxation, it is possible to detect the formation of kinks and to decide when the minimization succeeds. This action criterion is better than testing whether the final path is parallel to $\nabla E$, since the latter is also true for paths with kinks.

3.3 Determination of the initial path

Since the relaxation scheme works similar to a gradient minimization algorithm, it adjusts to the next local optimum of the action function. Therefore the choice of the initial path is crucial for obtaining the globally optimal path. Here we propose a method to find a good initial path $m^{0}(s)$ by a sequence of minimizations of modified energy functions.

For parameters $\mu,\beta,\varepsilon$ and $\Delta$ we consider the modified energy function

[TABLE]

The parameters $\mu,\beta,\varepsilon$ are chosen such that for $\Delta=d(m_{B},m_{A})$ the final state $m_{B}$ is a unique optimum of $E^{\ast}_{\Delta}(m)$, while for $m$ with $d(m,m_{A})\ll d(m_{B},m_{A})$ the last term of $E^{\ast}_{\Delta}(m)$ is small in comparison to $E(m)$. The value of $\mu$ should be large enough to ensure that the distance between $m_{A}$ and the minimum of $E^{\ast}_{\Delta}$ is indeed close to $\Delta$.

Now a sequence $m^{0}_{j},~{}j=1,\ldots,J$ of magnetization structures is iteratively obtained by setting $m^{0}_{0}~{}=~{}m_{A}$ and $m^{0}_{j}$ to the result of minimizing $E^{\ast}_{\Delta_{j}}(m)$, where $\Delta_{j}~{}=~{}j/J\,d(m_{B},m_{A})$.

Interpolating between the structures $m^{0}_{j}$ determines an initial path which (1) starts in $m_{A}$ and ends in $m_{B}$, (2) has relatively equally spaced intermediate states $m^{0}_{j}$, and (3) prefers intermediate states $m^{0}_{j}$ at distance $\Delta_{j}$ with low energy $E(m)$.

4 Energy barriers in a cubic pseudo-single domain particle

4.1 Material constants

The numerical calculations have been performed in terms of the reduced material parameters $Q$ and $\lambda$ [11]. The magnetic hardness $Q$, in the case of cubic magneto-crystalline anisotropy, is the numerical ratio of $Q=K_{1}/K_{d}$. Here $K_{d}$ is the characteristic magnetostatic self-energy, which in terms of the saturation magnetization $M_{s}$ is defined as $K_{d}=1/2\mu_{0}M_{s}^{2}$. The exchange length $\lambda=\sqrt{A/K_{d}}$ determines the characteristic length scale above which magnetostatic self-energy is able to overcome exchange coupling, represented by the exchange constant $A$. For magnetite these material constants are $M_{s}=480$ kA/m, $A=1.32\cdot 10^{-11}$ J/m, and $K_{1}=-1.25\cdot 10^{4}$ J/m3. Accordingly, $K_{d}=145$ kJ/m3, $\lambda_{ex}=9.55$ nm and magnetic hardness is $Q=-0.0863$. Unless otherwise stated, all results in the following subsections are obtained for the magnetite-like value $Q=-0.1$. Lengths $\lambda$ are given in units of $\lambda_{ex}$.

4.2 Potential barriers in cubic particles

$\lambda$ Initial state Initial energy $E$ Intermediate state Barrier $\varepsilon$

4.0 F-111 0.286404 F-110 0.007825

4.5 F-111 0.28308 V-011 0.005417

4.5 F-111 0.28308 F-011 0.008172

4.5 F-111 0.28308 V-001 0.011041

Table  1 : Energy barriers for numerically found optimal transitions between listed initial and final states for $Q=-0.1$ and different values of $\lambda$.

If to compare the barriers with the SD case, where $\varepsilon=1/120=0.008333$, the actual barriers are less, especially for $\lambda=4.5$. But it must be kept in mind that for $\lambda=4.0$ the transition is quasi-SD as the saddle point is F-110 type and almost homogeneous. For $\lambda=4.5$, however, the minimum energy saddle point is V-110 type, and development of vortex evidently considerably reduces the barrier.

Figure  1 : Energy variation across the optimal transition from a F-111 flower state to a F-1̄11 vortex state at $\lambda=4.0$. Each circle represents an intermediate magnetization state used for the calculation. The dashed line corresponds to the absolute value of the energy gradient along the transition path. The maximum energy state along the path is the F-011 flower state. The three-dimensional micromagnetic calculation ($\varepsilon_{FF}$)leads only to a minor decrease of the energy barrier with respect to coherent rotation ($\varepsilon_{coherent}$). This improvement results from the small spin deflections close to the particle surface.

For $\lambda=5.0$ the situation is most complex as there are a number of competitive LEM-states with similar energy. V001 with E = 0.27079, V110 with E = 0. 276931, V111 with E = 0.277701, F111 with E = 0.279527. As the result, there are no direct jumps between topologically identical states like V-001 and V-010. Such transitions were the only ones observed for lambda $\lambda=4.0$ and $\lambda=4.5$. For $\lambda=5.0$ the state V-001 can change to V-010 only indirectly, either through V-111, or V-011, or F-111. Thus, the diagonal elements of the matrix below are empty.

V-001 V-110 V-111 F-111

V-001 0 0.000123 0.001389 0.00349

V-110 0.006263 0 0.000598 0.002925

V-111 0.008297 0.001368 0 0.003836

F-111 0.012117 0.005521 0.005665 0

**Table  2 : Energy barriers for optimal transitions between a reduced set of initial and final states for $Q=-0.1$ and $\lambda=5.0$. Transitions with the same initial state are listed in the same column. Transitions with the same final state are listed in the same line. The corresponding complete set of 60 LEM states is obtained by taking into account cubic symmetry leading to 8 states of class F-111 or V-111, 6 states of class V-001, and 12 states of class V-110. In addition, all vortex states exist in two varieties of different helicity. **

Figure  2 : Energy variation across the optimal transition from a F-111 flower state to a V-111 vortex state at $\lambda=5.0$. Each circle represents an intermediate magnetization state used for the calculation. The dashed line corresponds to the absolute value of the energy gradient along the transition path. Because the energies of F-111 and V-111 are different, also the energy barriers for a transition from F-111 to V-111 ($\varepsilon_{FV}$) and from V-111 to F-111 ($\varepsilon_{VF}$) differ.

For $\lambda\geq 6.0$ the situation becomes simple again as now the only stable states are of type V-001.

$\lambda$ Initial state Intermediate state Barrier $\varepsilon$

6.0 V-100 V-111 0.009486

7.0 V-100 V-111 0.006372

8.0 V-100 V-111 0.005876

**Table  3 : Energy barriers for optimal transitions between vortex states for $Q=-0.1$ and $\lambda=6.0-8.0$. While V-111 is a marginally stable LEM at $\lambda=6.0$, it is unstable for $\lambda=7.0,8.0$. Types of V-100 vortex states are the global energy minima at these grain sizes and the energy barrier refers to a symmetric transition between two adjacent states of this type, e.g. V-100 to V-010. **

**Figure  3 : Energy variation across the transition from a V-111 flower state to a V-100 vortex state at $\lambda=6.0$. The top figures show the magnetization structure of the V-111 state as seen along the $\langle 111\rangle$-direction and the V-100 vortex as seen along the $\langle 100\rangle$-direction. In the bottom diagram again each circle represents an intermediate magnetization state used for the calculation. The dashed line corresponds to the absolute value of the energy gradient along the transition path. The tiny energy barrier $\varepsilon_{V111-V100}$ in relation to $\varepsilon_{V100-V111}$ indicates that the V-111 state is a very unstable LEM as compared to V-100. However, it is important for the transition between the more stable vortex states. E.g. the optimal transition from V-001 to V-100 is a combination of the symmetric vortex rotations V-001 to V-111, and the shown transition V111 to V100. **

Consider now the thermo-activation barrier $E_{b}=\varepsilon\,\lambda^{3}$ for the most stable LEM states: F-111 for $\lambda\leq 4.5$, and V-001 for $\lambda\geq 5.0$.

$\lambda$ $\varepsilon$ $E_{b}$ $m$

4.0 0.007825 0.5008 0.991297

4.5 0.005417 0.49362 0.986527

5.0 0.006263 0.782875 0.705928

6.0 0.009486 2.048976 0.509672

7.0 0.006372 2.185596 0.310451

8.0 0.005876 3.008517

Table  4 : Minimal energy density $\varepsilon$ or absolute energy $E_{b}$ necessary to leave the global energy minima for different values of $\lambda$. $m$ is the reduced magnetization at the global minimum.

**Figure  4 : Schematic representation of transition paths in magnetization space for a cubic particle with $\lambda=5.0$. Each sphere corresponds to the magnetization of an LEM state. The cubic structure reflects the cubic symmetry of the particle. Grey arrows indicate two of the many possible transition paths (or decay modes): (1) A direct decay form F-$\bar{1}\bar{1}\bar{1}$ to V-$\bar{1}\bar{1}\bar{1}$; (2) Indirect decay from F-$\bar{1}\bar{1}\bar{1}$ over an instable intermediate F-$\bar{1}0\bar{1}$ state into V-$\bar{1}0\bar{1}$ and further to either V-$\bar{1}00orV-00\bar{1}$. Note, that each sphere in principle represents two vortex states of inverse helicity (R and L). However, a transitions between vortex states of different helicity have large energy barriers and can be neglected. **

4.3 Statistical theory of MD VRM

The first application of the above presented calculation of energy barriers is a complete description of low-field viscous magnetization processes in a micromagnetically modelled cubic particle. Knowing the optimal transition paths between all LEM structures $S_{i}$ of the investigated particle allows for calculating the zero field temporal isothermal transition matrix $M(\Delta t)$, which describes the continuous homogeneous Markov process of random thermally activated transitions between all possible states:

[TABLE]

Here the matrix elements $\mu_{ij}$ of the infinitesimal generator of the semigroup $M(t)$ are given by the relative outflow from $S_{j}$ to $S_{i}$ for $i\neq j$. The relative inflow from all other states determines the diagonal element $\mu_{ii}$.

[TABLE]

Once, the $\mu_{ij}$ have been calculated, it is easily possible to determine the viscous decay of any initial probability distribution $\rho_{i}^{0}~{}\equiv~{}\mathbb{P}\left[{S_{0}=S_{i}}\right]$ by multiplication with the time evolution matrix exponential

[TABLE]

Multiplication by the corresponding magnetizations $m_{i}$ of states $S_{i}$ yields the viscous evolution of remanence

[TABLE]

When a small field $H$ is applied, the energy barrier $E_{b}^{ij}$ in first order changes according to

[TABLE]

where $m_{ij}^{\max}$ denotes the magnetization at the maximum energy state along the optimal transition path from $S_{i}$ to $S_{j}$. The approximation used to obtain (30) assumes that $H$ is so small that it does not change the magnetization structures of the LEM and saddle-point states noticeably. Only the field induced energy is taken into account. It is easily seen that all other energy changes are of second order in $H$.

Using the in-field energy barriers it is straightforward to determine the matrix exponential which governs VRM acquisition. By defining

[TABLE]

the above zero-field theory automatically extends to the weak field case.

In case of our cubic PSD particle, all matrices are of size $60\times 60$ and the calculations have been performed by a Mathematica (©Wolfram Research) program.

4.4 Viscous remanence acquisition and decay in an ensemble of cubic PSD magnetite

Using the above mathematical methods it is possible to calculate the statistics of viscous remanence acquisition and decay for our single PSD particle with respect to any field vector of sufficiently small length $H$. In order to model an isotropic ensemble, it is necessary to average the VRM properties over all possible field directions. This has been approximated by drawing 20 random directions from an equi-distribution over the unit sphere and averaging the modelled VRM acquisition and decay curves. For room temperature this yielded the ensemble curve as shown in Fig. 4.4.

Figure  5 : Modelled acquisition of viscous magnetization in the cubic particle with $\lambda=5.0$. The initial state is an equi-distribution over all possible LEM states with zero net magnetization. In a small external field the first acquisition process is the immediate decay from V-$110$ type states into V-$100$ type states, which occurs within about $10^{-9}$s. Due to the field induced asymmetry of the energy barriers, a remanence is acquired during this process. The second process is a decay of F-$111$ type states into V-$100$ type states. This occurs between about $10^{2}$s and $10^{3}$s and shows an intermediate overshooting of remanence.

The left hand side of Fig. 4.4 shows the remanence acquisition in a modelled field of $H=60\mu$T when starting from an initial state $\rho_{0}$ at $t=0$ which assigns equal probability to all existing LEM states. Already within $10^{-9}$s the remanence increases rapidly due to the immediate depletion of the nearly unstable V-$110$ vortex states which decay into the stable V-$100$ states (see Table 4.2). The remanence forms because in zero field there are two equally probable transitions, e.g. V-$110~{}\rightarrow$ V-$100$ and V-$110~{}\rightarrow$ V-$010$. Within the external field one of these decay paths becomes more probable which leads to a relative overpopulation of the field aligned V-$100$ type states. Nearly synchronously there occurs a two step process V-$111~{}\rightarrow$ V-$110~{}\rightarrow$ V-$100$.It is controlled by the somewhat slower transition V-$111~{}\rightarrow$ V-$110$, but still both take place within the first few $10^{-9}$s. The last VRM acquisition processes occurs only after a much longer waiting time of $10-10^{3}$s. First the initial F-$111$ type states transform via F-$110$ type states into V-$110$ type states which then immediately decay into V-$100$ (Fig. 4.2). This last process produces an astonishing remanence overshooting as displayed in Fig. 4.4: The remanence during the VRM acquisition process is for a certain time higher than the finally obtained equilibrium VRM. In the next section we will show that is is not an artifact of the modelling, but can be explained by a real physical process.

The right hand side of Fig. 4.4 shows that when the field is switched off after VRM acquisition, the obtained remanence is carried only by extremely stable V-$100$-type states which require a theoretical waiting time of $10^{15}$s to equilibrate into a zero remanence state.

5 Viscous magnetization

5.0.1 Viscous magnetization anomalies

One of the most astonishing results of this study with respect to viscous magnetization is the predicted transient increase of remanence during the VRM acquisition. To understand this phenomenon more closely, we give a physical explanation of this effect in terms of a simplified model.

The in field potential barrier is

[TABLE]

Here $E_{s}$ is the energy at the saddle point and $E_{F}$ is the zero field energy of the flower state. $W+$ and $W-$ are the number of grains (probability) in the vortex $V-100$ parallel and antiparallel to the external field $h$ in $x$-direction.

5.0.2 High stability of PSD VRM

Extremely stable VRM has been often observed in paleomagnetic studies. The above mechanisms give a first theoretical explanation why high stability of VRM should occur in PSD samples. The basic process is the relaxation of naturally produced metastable states into stable ones.

This typically occurs for TRM acquisition in PSD ensembles where the cooling rate is fast enough to stabilize metastable flower states. Then a long term VRM, acquired in the field after cooling and carried by newly formed vortex states is extremely stable and can significantly bias any paleomagnetic measurement, especially paleointensity determinations.

Laboratory AF demagnetization rather leads to a more stable LEM state (perhaps even the GEM) because the magnetization structure is provided with a lot of energy which is stepwise reduced. Thus AF-demagnetization would rather end up in a vortex state for a PSD particle.

During natural chemo-viscous (VCRM) magnetization by crystal growth a grain changes sequential from the SP state into a stable SD and later a PSD state. The first stable SD state is almost homogeneously magnetized along a $\langle 111\rangle$-axis. It then transforms in a more developed flower state which then becomes metastable as soon as the vortex has lower energy. At this point, the process of VRM acquisition will starts to produce extremely stable remanences.

6 Appendix

Details of calculation

This uses tensor calculus and Einstein’s sum convention and for comparison with classical physics interprets the variable $w$ as a time $t$ and $x$ as a generalized variable $q$. We have to apply the Euler-Lagrange operator $\frac{\partial}{\partial q}-\frac{d}{dt}\frac{\partial}{\partial\dot{q}}$ to the Lagrange function $L(q,\dot{q})~{}=~{}\sqrt{g_{i}\,g_{i}}\,\sqrt{\dot{q}_{i}\dot{q}_{i}}$, where $g_{i}~{}=~{}\partial_{i}E(q)$. This needs the following expressions

[TABLE]

[TABLE]

The last equation uses the fact that in the chosen parametrisation the tangent vector $\dot{q}$ has constant length along the path and thus is perpendicular to $\ddot{q}$. Further,

[TABLE]

Putting this together results in

[TABLE]

Using the arc length parametrisation where $\|\dot{q}\|=1$, this finally leads to the Euler-Lagrange equation

[TABLE]