format: text

text:
# Langevin Simulation of Thermally Activated Magnetization Reversal in Nanoscale Pillars

## I Introduction

Several emerging technologies will incorporate fabricated magnets that are small enough to contain only a single magnetic domain. These nanoscale magnets will be essential for smaller components, lower power consumption, and completely new applications in fields such as information storage, integrated circuits, sensor technology, and microelectromechanical systems. Successful implementation of these technologies requires a fundamental understanding of the dynamics of the internal magnetic structure of the nanoscale magnets on time scales ranging from the nanoseconds associated with gigahertz applications to the years over which magnetic information storage must be stable.

An essential factor in many of these applications, especially in information storage, is the free-energy barrier that separates two antiparallel orientations of the magnetization. This free-energy barrier can be surmounted using thermal energy momentarily “borrowed” from the surroundings, and often device engineers strive for barriers of at least $`40k_\mathrm{B}T`$ to make this thermal bit-flipping a rare event. Here $`T`$ is the absolute temperature and $`k_\mathrm{B}`$ is Boltzman’s constant. The particle can be magnetizes in a specified direction by applying a field strong enough to remove the free-energy barrier. The smallest field sufficient to do this at zero temperature is called the coercive field $`H_\mathrm{c}.`$ Fields smaller than the coercive field cannot change the magnetization orientation deterministically, but they do cause lower free-energy barriers and make thermal crossing more probable. In fact, hybrid recording  uses lower-than-coercive fields and thermal barrier crossing to write data in high-coercivity magnetic media.

Two examples of particular thermally activated barrier crossings are shown in Fig. 1. Red is associated with the magnetization orientation antiparallel to the applied field, which is metastable. Blue is associated with the parallel magnetization, which is thermodynamically stable. Green and yellow are intermediate between the two orientations. For times before those shown in the figure, the magnet stays in the metastable state with end caps of non-uniform magnetization induced by pole avoidance. Eventually thermal fluctuations in the end caps carry the system over the free-energy barrier. After the barrier crossing, the entire magnet changes quickly to the stable equilibrium orientation.

Often, nanoscale magnets are assumed to be uniformly magnetized bodies. For instance, Garcìa-Palacios and Lázaro  considered the stochastic trajectories of isolated magnetic moments using the Landau-Lifshitz-Gilbert equation subject only to applied fields and a uniaxial anisotropy, which supplied the free-energy barrier. They measured the susceptibility of the magnetization in small sinusoidal probe fields. Multiple crossing of the free-energy barrier, a frequently noted consequence of the gyromagnetic motion in single-spin models, was also observed for their uniform-magnetization model.

When at least one dimension of the nanomagnet is greater than the exchange length of the material, nonuniform reversal modes become energetically possible. To consider such nonuniform magnetization reversal we use micromagnetic simulation with a large number of points inside the magnet. Specifically, nanomagnets with an aspect ratio of approximately $`17`$ are considered throughout this work, and are referred to as pillars. The numerical approach is discussed in Sec. II, with further details found in the appendices. Results for the estimation of the coercive field and the statistics of magnetization switching are discussed in Sec. III. A simple model for nonuniform switching is discussed in Sec. IV. Section V considers the effect of interactions between pillars in a one-dimensional array. A brief summary of the results is given in Sec. VI. Some preliminary results of this study were presented in Refs. .

## II Numerical Details

The basic approach of micromagnetic modeling is to consider a system of coarse-grained magnetization vectors $`\stackrel{~}{𝐌}(\stackrel{~}{𝐫}_i)`$, with $`\stackrel{~}{𝐫}_𝐢`$ indicating the location in space of the $`i`$-th spin. The vectors are assumed to have a fixed magnitude $`M_s`$, corresponding to the bulk saturation magnetization density of the material. This is a valid approximation for temperatures well below the Curie temperature.  The evolution of each magnetization vector is governed by the damped precessional motion given by the Landau-Lifshitz-Gilbert (LLG) equation

$$\frac{d\stackrel{~}{𝐌}(\stackrel{~}{𝐫}_i)}{d\stackrel{~}{t}}=\frac{\gamma _0}{1+\alpha ^2}\stackrel{~}{𝐌}(\stackrel{~}{𝐫}_i)\times \left(\stackrel{~}{𝐇}(\stackrel{~}{𝐫}_i)\frac{\alpha }{M_s}\stackrel{~}{𝐌}(\stackrel{~}{𝐫}_i)\times \stackrel{~}{𝐇}(\stackrel{~}{𝐫}_i)\right)$$
(1)
where the electron gyromagnetic ratio is $`\gamma _0=1.76\times 10^7\mathrm{Hz}/\mathrm{Oe}`$,  and $`\alpha `$ is a phenomenological damping parameter. For the sign of the undamped-precession term we follow the convention of Brown.  Here a tilde is used to distinguish dimensional quantities from their dimensionless counterparts, which will be introduced later.

The dynamics are controlled by $`\stackrel{~}{𝐇}(\stackrel{~}{𝐫}_i)`$, the local field at the $`i`$-th position, which, in general, is different at each lattice site. The local field mediates all of the interactions in the system with the contributions combined via linear superposition,

$$\stackrel{~}{𝐇}(\stackrel{~}{𝐫}_i)=\stackrel{~}{𝐇}_\mathrm{z}(\stackrel{~}{𝐫}_i)+\stackrel{~}{𝐇}_\mathrm{e}(\stackrel{~}{𝐫}_i)+\stackrel{~}{𝐇}_\mathrm{d}(\stackrel{~}{𝐫}_i)+\stackrel{~}{𝐇}_\mathrm{a}(\stackrel{~}{𝐫}_i)+\stackrel{~}{𝐇}_\mathrm{n}(\stackrel{~}{𝐫}_i),$$
(2)
where $`\stackrel{~}{𝐇}_\mathrm{z}(\stackrel{~}{𝐫}_i)`$ represents the externally applied field (Zeeman term), $`\stackrel{~}{𝐇}_\mathrm{e}(\stackrel{~}{𝐫}_i)`$ is due to exchange effects, $`\stackrel{~}{𝐇}_\mathrm{d}(\stackrel{~}{𝐫}_i)`$ is the dipole field, $`\stackrel{~}{𝐇}_\mathrm{a}(\stackrel{~}{𝐫}_i)`$ is due to crystalline anisotropy, and $`\stackrel{~}{𝐇}_\mathrm{n}(\stackrel{~}{𝐫}_i)`$ is the random field induced by thermal noise.  The estimates used for each of these fields are given below, except for $`\stackrel{~}{𝐇}_\mathrm{a}(\stackrel{~}{𝐫}_i)`$, which we have taken to be zero. In this article, a $``$ emphasizes that the symbol represents a quantity with units, the corresponding dimensionless quantity is written without the tilde. Material parameters, such as $`M_s,`$ always have units when appropriate.

We have chosen material parameters to match those of bulk iron; a summary of all parameters used here is given in Table 1. The saturation magnetization density is $`M_s=1700\mathrm{emu}/\mathrm{cm}^3,`$  while the exchange length, the length over which $`\stackrel{~}{𝐌}`$ can change appreciably, is $`l_\mathrm{e}=3.6\mathrm{nm}.`$  The damping parameter $`\alpha `$ has proven difficult to measure or estimate from ab initio considerations, and may even depend on numeric details.  Here we have chosen $`\alpha =0.1`$ to represent the under-damped behavior usually assumed to exist in nanoscale magnets. In the other extreme, the over-damped limit, the gyromagnetic motion intrinsic to the LLG equation is suppressed, and the system can be more efficiently simulated using Monte Carlo methods.

The magnetization was discretized on a cubic lattice with discretization length $`\mathrm{\Delta }\stackrel{~}{r}=1.5\mathrm{nm}`$ chosen to give sufficient resolution across the cross section of the pillar. The time discretization of $`\mathrm{\Delta }\stackrel{~}{t}=5\times 10^{14}\mathrm{s}`$ was chosen to avoid numerical instabilities in the integration. The details of the integration scheme are given below. The equation of motion was cast into dimensionless quantities by considering the rescaled length $`𝐫=\stackrel{~}{𝐫}/l_\mathrm{e},`$ time $`t=\gamma _0M_s\stackrel{~}{t},`$ field $`𝐇=\stackrel{~}{𝐇}/M_s,`$ and magnetization density $`𝐌=\stackrel{~}{𝐌}/M_s.`$ This rescaling is facilitated in cgs units,  where magnetic fields can be normalized by the magnetization density without the inclusion of any constants.

In a continuum model, exchange represents local differences in the alignment of the magnetization. The contribution to the local field from exchange interactions can be approximated by $`l_\mathrm{e}^2^2\stackrel{~}{𝐌}(\stackrel{~}{𝐫})`$ which has been implemented in the simulations by

$$\stackrel{~}{𝐇}_\mathrm{e}(\stackrel{~}{𝐫}_i)=\left(\frac{l_\mathrm{e}}{\mathrm{\Delta }\stackrel{~}{r}}\right)^2\left(6\stackrel{~}{𝐌}(\stackrel{~}{𝐫}_i)+\underset{|\stackrel{~}{𝐝}|=\mathrm{\Delta }\stackrel{~}{r}}{}\stackrel{~}{𝐌}(\stackrel{~}{𝐫}_i+\stackrel{~}{𝐝})\right),$$
(3)
where the summation is over the six nearest neighbors of $`\stackrel{~}{𝐫}_i.`$ Here the exchange length is defined in terms of the exchange energy  $`E_\mathrm{e}=(ł_e^2/2)𝑑𝐫\stackrel{~}{𝐌}^2\stackrel{~}{𝐌}.`$ The dimensionless form of the exchange contribution to the local field is

$$𝐇_\mathrm{e}(𝐫_i)=\left(\frac{1}{\mathrm{\Delta }r}\right)^2\left(6𝐌(𝐫_i)+\underset{|𝐝|=\mathrm{\Delta }r}{}𝐌(𝐫_i+𝐝)\right).$$
(4)
The computationally most intensive part of determining the local field is finding the contribution from dipole-dipole interactions, $`𝐇_\mathrm{d}(𝐫)`$. Efficient calculation, which uses the Fast Multipole Method (FMM),  is quite involved and is discussed in Appendix 1.

Thermal fluctuations are included in the LLG equation of motion by inclusion in the superposition of a random field $`\stackrel{~}{𝐇}_\mathrm{n}(\stackrel{~}{𝐫})`$ with Gaussian fluctuations whose first moments are zero and whose second moments obey the fluctuation-dissipation relation

$$\stackrel{~}{H}_{\mathrm{n}\mu }(\stackrel{~}{𝐫}_i,\stackrel{~}{t})\stackrel{~}{H}_{\mathrm{n}\mu ^{}}(\stackrel{~}{𝐫}_i^{},\stackrel{~}{t}^{})=\frac{2\alpha k_\mathrm{B}T}{\gamma _0M_sV}\delta \left(\stackrel{~}{t}\stackrel{~}{t}^{}\right)\delta _{\mu ,\mu ^{}}\delta _{i,i^{}},$$
(5)
where $`\stackrel{~}{H}_{\mathrm{n}\mu }`$ indicates one of the Cartesian components of $`\stackrel{~}{𝐇}_\mathrm{n}.`$ Here $`V=(\mathrm{\Delta }\stackrel{~}{r})^3`$ is the discretization volume of the numerical implementation and $`\delta _{\mu ,\mu ^{}}`$ is the Kronecker delta representing the orthogonality of Cartesian coordinates. This result was derived for isolated particles , and interactions between the discretization volumes can be important.  To make the simulations tractable, the effect of interactions on the thermal noise has been neglected. Here $`\delta (\stackrel{~}{t}\stackrel{~}{t}^{})`$ is the Dirac delta function, and its dimensions $`[\delta (tt^{})]=1/\mathrm{s}`$ are important. After the transformation to dimensionless quantities, the fluctuation-dissipation result is

$$H_{\mathrm{n}\mu }(𝐫_i,t)H_{\mathrm{n}\mu ^{}}(𝐫_i^{},t^{})=ϵ\delta \left(tt^{}\right)\delta _{\mu ,\mu ^{}}\delta _{i,i^{}},$$
(6)
with the dimensionless strength of the stochastic field given by

$$ϵ=\frac{2\alpha k_\mathrm{B}T}{M_\mathrm{S}^2V},$$
(7)
where the untransformed discretization volume and saturation magnetization are both used in cgs units to yield the dimensionless result. The dimensionless form of the LLG equation is

$$\frac{d𝐌(𝐫_i)}{dt}=\frac{1}{1+\alpha ^2}𝐌(𝐫_i)\times \left[𝐇(𝐫_i)\alpha 𝐌(𝐫_i)\times 𝐇(𝐫_i)\right].$$
(8)
This dimensionless stochastic differential equation is used for the numerical integration.

Since the LLG equation conserves the magnitude of the magnetization density, each integration step amounts to a rotation. The size of the tangential displacement at each integration step is given by

$$𝐑(𝐫_i,t)=\frac{1}{1+\alpha ^2}𝐌(𝐫_i,t)\times [(𝐈(𝐫_i,t)\alpha 𝐌(𝐫_i,t)\times 𝐈(𝐫_i,t)],$$
(9)
where $`𝐈`$ is the “impulse” over the integration step defined below. The magnetization after the integration step is given by

$$𝐌(𝐫_i,t+\mathrm{\Delta }t)=\frac{𝐌(𝐫_i,t)+𝐑(𝐫_i,t)}{\sqrt{1+𝐑^2(𝐫_i,t)}},$$
(10)
which ensures the conservation of unit magnitude for $`𝐌`$. If only the deterministic dynamics were important a high-order integration method could be used, but to correctly include the thermal noise a lower-order method must be used. Using first-order Euler integration, the impulse including the fluctuating field is

$$𝐈(𝐫_i,t)=\left[𝐇_\mathrm{a}(𝐫_i,t)+𝐇_\mathrm{d}(𝐫_i,t)+𝐇_\mathrm{e}(𝐫_i,t)\right]\mathrm{\Delta }t+\sqrt{ϵ\mathrm{\Delta }t}𝐠(𝐫_i,t),$$
(11)
where $`𝐠(𝐫,t)`$ is a random vector with each component chosen independently from a Gaussian distribution of zero mean and unit variance. This result, including the $`\sqrt{\mathrm{\Delta }t}`$ contribution characteristic of integration of stochastic processes, is explained in Appendix 2.

## III Isolated Nanomagnets

The numerical micromagnetic simulation methods described in the previous section have been applied to model nanoscale magnets inspired by those recently fabricated using a combination of chemical vapor deposition and scanning tunneling microscopy techniques by Wirth, et al.  The technique has been used to produce arrays of nanoscale magnetic particles with diameters down to about $`10`$ nm and lengths from $`50`$ to $`250`$ nm. The reversal of pillars with diameters more than twice the exchange length, such as these, have been found numerically to proceed by a mode where the magnetization is not constant across the diameter of the pillar.  The model particles considered here are rectangular prisms $`9`$ nm $`\times `$ $`9`$ nm $`\times `$ $`150`$ nm, which require $`N=4949`$ sites on the computational lattice. The long axis is taken as the $`z`$-axis, and the cross-sectional area is chosen comparable to that for a $`10`$ nm diameter pillar. The external field is always applied along the $`z`$-axis – the easy axis induced by shape anisotropy. As mentioned in Sec. II, the material parameters were chosen to match those of bulk iron. These are consistent with measurements on the experimentally produced iron nanoparticles.

### A Coercive Field

Images of the $`z`$-component of the magnetization, $`M_z`$, for such a particle as described above are shown in Fig. 2(a) for several fields during a hysteresis-loop simulation (in a one-quarter cut-away view). The color is a linear scale of $`M_z,`$ shown in the legend of Fig. 1 with the most negative values at the bottom and the most positive values at the top. For hysteresis-loop simulations, the field is always oriented along the $`z`$ axis. Starting from its most positive extreme, the value is varied sinusoidally. For the simulation shown in Fig. 2 a period of one nanosecond, a field amplitude of $`4000\mathrm{Oe}`$, and a temperature of $`0\mathrm{K}`$ were selected. The response of the initially upward-directed magnetization after the applied field is oriented downward can be clearly seen. First, end caps associated with pole avoidance form at both ends. These end caps are regions of high curl in the magnetization density, $`𝐂(𝐫)=\times 𝐌(𝐫)`$. The end caps are obvious in Fig. 2(b), where $`\mathrm{sign}\left(C_z(𝐫,t)\right)|𝐂(r,t)|`$ is shown on a linear color scale, so that right (left)-handed end caps are colored red (blue). At zero temperature, the two end caps then grow symmetrically until they meet near the midpoint of the particle. While they grow, the area of large curl is concentrated at the interface between the volume where the magnetization has already aligned parallel with the applied field and the volume where it remains antiparallel to the field. Each interface is composed of two regions with opposite $`z`$-component of the curl, with a curling type of singularity in the center of the pillar. These two regions are separated by a layer where the $`z`$-component of the curl is close to zero, and the magnetization vectors in the $`x`$$`y`$-plane have a negative divergence located at the center of the pillar. Because the end caps have opposite curls, in the region where the end caps come into contact the curl changes abruptly from large positive to large negative values. Some time is required for this defect to disappear. This reversal of pillars by nucleation of reversed volumes at the ends followed by growth of the reversed regions is consistent with minimization of the micromagnetic energy at zero temperature  and Monte Carlo simulations at finite temperature,  both under quasistatic field-sweep conditions.

The hysteresis loop associated with this simulation is shown in Fig. 3, with the field values shown in Fig. 2 indicated by large tick marks. Two periods are shown for the $`1\mathrm{ns}`$ loop. Aside from differences due to the initial alignment of the magnetization and the presence of the defect after a complete loop, the reproducibility at zero temperature is excellent. Hysteresis loops under the same conditions, but for periods of $`2`$ and $`4`$ ns, are also shown. These simulations show that the hysteresis loop becomes more square for longer periods, indicating that the magnetization is not following the applied field in a quasi-static fashion. This is quite reasonable given that the simulated loops correspond to microwave frequencies, but it makes the infinite-period coercive field difficult to measure.

The average energy density, along with its contributions, as functions of time appear in Fig. 4 for the $`4`$ ns hysteresis loop. The energy density for zero crystalline anisotropy is calculated as

$$E=\frac{1}{N}\underset{i=1}{\overset{N}{}}𝐌(𝐫_i)\left(\frac{1}{2}\left[𝐇_\mathrm{e}(𝐫_i)+𝐇_\mathrm{d}(𝐫_i)\right]+𝐇_\mathrm{z}(𝐫_i)\right),$$
(12)
where $`N`$ is the number of lattice sites. As the applied field begins to decrease and then become negative, it is the Zeeman energy that changes the most while the exchange and dipole-dipole energies remain nearly constant. Near the coercive field, the exchange energy rapidly increases as regions of reversed magnetization, and the interfaces associated with them, form at the ends of the pillar. Once these reversed regions become large enough, they grow spontaneously and the Zeeman energy rapidly decreases, while the exchange energy from the interface between reversed and unreversed regions remains roughly constant. The growth of the reversed regions does not take long, and the exchange energy quickly dissipates when the regions from the two ends merge at the middle and the reversal is complete.

The information in Fig. 4 can be used for extracting the coercive field of the nanomagnet. For infinitely slow variations in the applied field, $`dE/dt`$ is proportional to $`dE/dH_\mathrm{a}`$ using the chain rule. Thus portions of the figure with positive slope of the total energy correspond to fields for which the system is trapped behind an energy barrier, while portions with negative slope correspond to fields for which the system responds deterministically by finding a new energy-minimizing magnetization configuration. Fields for which the energy has a maximum should then correspond to the coercive field. Dynamic effects cause the field at which this maximum occurs to depend on the frequency.

To estimate the static coercive field we find $`H_c(R)`$ for fields that vary linearly with time, $`H(t)=Rt`$ for $`t>0`$, and then extrapolate to the $`H_c(R=0)`$ value. The energy density is presented as a function of the absolute value of the applied field in Fig. 5 for rates ranging from $`R=250\mathrm{Oe}/\mathrm{ns}`$ to $`R=10000\mathrm{Oe}/\mathrm{ns}.`$ The coercive fields estimated from the energy maximum are shown in the inset where the error bars are estimated from the second derivative near the energy maximum. The coercive field clearly decreases with $`R.`$ Extrapolation of a weighted least-squares fit yields a coercive field of $`H_c=1979\pm 14\mathrm{Oe}`$ at $`R=0.`$ This is a more accurate estimate of the static coercive field than those taken directly from hysteresis-loop simulations.

### B Switching in Constant Field

At finite temperature the magnetization can reverse, even when the applied field is weaker than the coercive field. In this case, thermal fluctuations carry the magnetization past the free-energy barrier and on to the new equilibrium configuration. Two examples of such thermally induced switching are shown in Fig. 1 for the conditions $`H_0=1800`$ Oe and $`T=20`$ K. For small positive times (not shown) the volumes associated with the end caps fluctuate due to the thermal noise. The reversal appears to proceed by a nucleation process, with the end caps serving as the seeds for heterogeneous nucleation. The initiation of the switching from the end caps is similar to results for reversal induced by applied fields greater than $`H_c`$ in experiments,  and simulations with $`T`$$`=`$$`0\mathrm{K}`$ where reversal begins at the end caps  or corners.

In classical droplet theory, in which boundaries such as the sides of the magnet are ignored, nucleation of the equilibrium magnetization is governed by the competition between the favorable Zeeman energy due to alignment with the applied field and the unfavorable exchange energy due to the interface with the majority volume of the misaligned orientation. The exchange energy dominates for small droplets, which tend to shrink. Large droplets tend to grow because the Zeeman energy dominates. These two regimes are separated by a critical droplet size, the saddle point associated with the free-energy barrier, where the tendencies toward growth and shrinkage are balanced. In nanoscale magnets the situation is more complicated because the free energy cannot be easily separated into surface and volume contributions. Nevertheless, the important aspect remains that a free-energy barrier exists that must be crossed for reversal to begin. The end caps fluctuate until a succession of highly unlikely fluctuations carries the end-cap configuration past the free-energy barrier. After that, growth of the reversed region is energetically favorable and occurs rapidly, on the order of $`0.2`$ to $`0.3`$ ns for the pillars.

For long pillars, nucleation events at the two ends occur approximately independently of each other. In addition, the growth of a supercritical region takes a significant amount of time: enough time for the other end cap to have a reasonable probability to nucleate and begin to grow. One example of both ends nucleating at different times is shown in Fig. 1(a), while Fig. 1(b) shows nucleation of one end that grows to switch the magnetization of the entire pillar. For the present applied field the latter situation is somewhat rare, from the theory presented below it can be inferred that it occurs in about 10% of simulated switches. From Fig. 1(b) it can also be seen that the nucleation of the two end caps is not completely independent. As one reversed region grows along the pillar, the free-energy barrier at the other end cap decreases. Eventually the free-energy barrier for the second end becomes zero and it will also begin to grow. If a less negative field is applied, the free-energy wells of the metastable states will be deeper and the first reversed region will have to grow further along the pillar before this occurs.

The mean switching time, $`\overline{t}_{\mathrm{sw}}`$, defined here as the first time when $`M_z=0,`$ is shown in Fig. 6(a) as a function of applied field for $`T=20\mathrm{K}.`$ The difference in the error bars reflects the amount of statistical sampling at the different conditions. The dependence of $`\overline{t}_{\mathrm{sw}}`$ on temperature is shown in Fig. 6(b) for applied fields of $`H_0=1850\mathrm{Oe};`$ fields are always applied along the pillar. While there is a clear dependence on field and temperature, we find no clear functional form for the dependence. Notably, the exponential dependence of $`\overline{t}_{\mathrm{sw}}`$ on inverse temperature, expected when the barrier is high, is not seen. This implies that the free-energy barrier is low. In fact, this picture is consistent with the switching statistics described next, where the dynamics are consistent with a biased walk in a shallow well.

The statistics of magnetization switching are measured by $`P_{\mathrm{not}}(t)`$, the probability of not switching before time $`t`$. The simulation results for $`P_{\mathrm{not}}(t)`$ for $`200`$ switches at $`H_0=1800\mathrm{Oe}`$ and $`T=20\mathrm{K}`$ are presented in Fig. 7(a). The result is clearly not exponential, the expected functional form when the free-energy barrier is large. Experimental results for single-domain magnets, in which $`P_{\mathrm{not}}(t)`$ is not exponential, have been reported recently.

A simple theory that assumes completely independent nucleation at the two ends of the pillar produces a reasonable description of the simulation result. The theory for $`P_{\mathrm{not}}(t)`$ requires two parameters. The first is the constant nucleation rate, $``$, for the formation of a supercritical reversed region at one isolated end of a pillar. The second is the rate, $`v`$, at which a single growing region changes the normalized magnetization of the pillar. Labeling the two independent nucleation times at the top and bottom $`t_\mathrm{t}`$ and $`t_\mathrm{b}`$, respectively, the positive values of the magnetization of the pillar can be described by

$$M_z(t)=\{\begin{array}{cc}1\hfill & \hfill t<t_1\\ 1v(tt_1)\hfill & \hfill t_1t<t_2\\ 1v(tt_1)v(tt_2)\hfill & \hfill t_2t\end{array},$$
(13)
where $`t_1=\mathrm{min}(t_\mathrm{t},t_\mathrm{b})`$ and $`t_2=\mathrm{max}(t_\mathrm{t},t_\mathrm{b})`$. We define switching to occur at $`M_z(t_{\mathrm{sw}})=0`$, when roughly one-half the volume of the nanomagnet is oriented in the equilibrium direction. The switching can occur either through one or both ends nucleating, and the actual switching time for a particular $`(t_\mathrm{t},t_\mathrm{b})`$ is given by

$$t_{\mathrm{sw}}=\mathrm{min}(2t_0+t_1,t_0+\frac{t_1+t_2}{2}),$$
(14)
where $`t_0=1/(2v)`$ is the earliest time at which the pillar can switch. For times $`t<t_0`$ the probability of switching is zero, even if both ends nucleate immediately after the reversal. The situation in the $`t_\mathrm{b}`$-$`t_\mathrm{t}`$ plane is shown in Fig. 8, where the dashed curves correspond to two different conditions of constant $`t_{\mathrm{sw}},`$ $`t_{\mathrm{sw}}/t_0`$$`=`$$`3/2`$ and $`3.`$ The solid lines divide the plane into regions with different switching histories. In the triangle near the origin (I), switching occurs for times $`t_0<t_{\mathrm{sw}}<2t_0`$, which can happen only if both ends nucleate. Between the parallel solid lines (II) switching occurs by double nucleation, while outside them (III & IV) a single nucleation causes the switching. Since each nucleation process has a constant rate, the corresponding probability density function for the nucleation time is exponential, $`\mathrm{exp}(t)`$. Integrating along curves of constant $`t_{\mathrm{sw}}`$, the probability of not switching before time $`t`$ is found to be

$$P_{\mathrm{not}}\left(t\right)=\{\begin{array}{cc}1\hfill & \hfill t<t_0\\ e^{2(tt_0)}\left[1+2(tt_0)\right]\hfill & \hfill t_0t<2t_0\\ e^{2(tt_0)}\left[1+2t_0\right]\hfill & \hfill 2t_0t\end{array}.$$
(15)
We will refer to this as the “two-exponential” decay model. In the limit of infinitely fast growth, $`t_0=0,`$ only regions III and IV remain, and this becomes a simple exponential decay associated with thermally activated crossing of a free-energy barrier. For $`t_0>0`$, $`P_{\mathrm{not}}(t)`$ is quadratic in $`tt_0`$ for $`tt_0^+`$, has a discontinuity in the slope at $`t=2t_0`$, and is exponential for $`t>2t_0`$.

The fit of the two-exponential decay model of Eq. (15) to the simulation results appears in Fig. 7(a) as the dashed curve. The dotted curve is a fit to the error function

$$P_{\mathrm{erf}}(t)=\frac{1}{2}\frac{1}{2}\mathrm{Erf}\left(\stackrel{~}{}(t\stackrel{~}{t}_0)\right),$$
(16)
which was chosen as an empirical form because it also has two adjustable parameters, chosen here to have a similar interpretation as the theory given above. The parameters for both theoretical forms were determined by matching the first and second moments of the theoretical form to the moments of the simulation data, and the parameters are $`4.049`$, $`t_00.527`$ and $`\stackrel{~}{}4.176`$, $`\stackrel{~}{t}_00.790`$. For these conditions, both forms, Eqs. (15) and (16), seem to fit the simulation data about equally well. The disagreement with Eq. (15) probably occurs because the condition that the free-energy barrier is much greater than $`k_\mathrm{B}T`$ is not met. When the free-energy barrier is much less than $`k_\mathrm{B}T,`$ the magnetization should be more like a biased random walk. In that case the most noticeable effect of thermal fluctuations would be a symmetric smearing of the switching times around $`\overline{t}_{\mathrm{sw}}.`$

The question of statistical sampling is also important. Specifically, one expects to have undersampling of the population because the finite-sampling $`P_{\mathrm{not}}(t)`$ goes to zero faster than the infinite-sampling $`P_{\mathrm{not}}(t)`$ with increasing $`t`$. In fact, for times less than the maximum observed switching time in the long-time tail of the distribution, the error function underestimates the population, while the two-exponential theory overestimates it. To this extent, it is possible that data with better sampling may shift towards better agreement with Eq. (15). Indeed, experimental switching probabilities  in single-domain magnets have been observed to have exponential tails, which is inconsistent with Eq. (16).

Better statistics and a higher free-energy barrier are needed to validate the theoretical $`P_{\mathrm{not}}(t)`$ given in Eq. (15). For the present model, better statistics would be prohibitive; the $`200`$ switches considered here took approximately $`20`$ weeks to generate using $`4`$ processors on an Origin 2000. A less realistic but numerically faster simulation, in which the pillars are modeled as one-dimensional stacks of cubes, is presented in the next section. The advantage of this simplified model is that it allows many more switching events to be observed in a reasonable amount of computer time, even with much higher free-energy barriers. This improves the statistics for the very rare events with $`t_{\mathrm{sw}}2t_0`$ and $`t_{\mathrm{sw}}t_0`$, which are important to distinguish between alternative theoretical forms for $`P_{\mathrm{not}}(t)`$.

## IV A Simple Model

In order to generate a sufficiently large sample of switching events in nanoscale magnets, we consider a simpler model in which the pillar consists of a one-dimensional stack of cubic cells, which is the model introduced by Boerner and Bertram.  The cubes have a linear size $`\mathrm{\Delta }\stackrel{~}{r}=2\mathrm{}_\mathrm{e},`$ here $`5.2\mathrm{nm},`$ since the material parameters for iron are used (Ref.  used parameters for nickel). To keep the aspect ratio of the nanopillars as close as possible to those considered above, a pillar composed of $`17`$ cubes is used. For this simple model, the local field due to the dipole interactions is calculated as

$$\stackrel{~}{𝐇}_\mathrm{d}^{}(\stackrel{~}{𝐫}_i)=(\mathrm{\Delta }\stackrel{~}{r})^3\underset{ji}{}\frac{3\widehat{𝐫}_{ij}\left(\widehat{𝐫}_{ij}\stackrel{~}{𝐌}(\stackrel{~}{𝐫}_j)\right)\stackrel{~}{𝐌}(\stackrel{~}{𝐫}_j)}{\stackrel{~}{𝐫}_{ij}^3},$$
(17)
where $`\stackrel{~}{𝐫}_{ij}`$ is the displacement vector from the center of cube $`i`$ to the center of cube $`j,`$ and $`\widehat{𝐫}_{ij}`$ is the corresponding unit vector. The factor of volume $`(\mathrm{\Delta }\stackrel{~}{r})^3`$ results from integrating over the constant magnetization density in each cell. This factor can be combined with $`\stackrel{~}{r}_{ij}^3`$, so that the denominator depends only on $`𝐱_{ij}=\stackrel{~}{𝐫}_{ij}/(\mathrm{\Delta }\stackrel{~}{r})`$, which is a vector of integers denoting the difference in cell indices. It can readily be seen that for a uniformly magnetized pillar the summation over index integers leads to different values of $`𝐇_\mathrm{d}^{}(𝐫)`$ for different choices of discretization length along the long axis. This happens because self-contributions have been ignored and the far-field result has been used for neighboring cells. However, since this model is being used to investigate only the statistics of switching, and not for estimating physical values, the results should be qualitatively correct. Using the scaling to dimensionless units, the approximate dipole field is

$$𝐇_\mathrm{d}^{}(𝐫_i)=\underset{ji}{}\frac{3\widehat{𝐫}_{ij}\left(\widehat{𝐫}_{ij}𝐌(𝐫_j)\right)𝐌(𝐫_j)}{𝐱_{ij}^3}.$$
(18)
We have verified that our implementation of the model agrees with that of Boerner and Bertram by reproducing the observed coercive field for the nickel pillars discussed in Ref. . Because of the different approximations, the coercive field for the present iron pillars is about $`1500\mathrm{Oe}.`$

The probability of not switching for $`2000`$ switches for this model of iron pillars is shown in Fig. 7(b), along with the two theoretical forms, Eqs. (15) and (16), which have been fit using the same procedure as in Sec. III. The kink in the two-exponential theoretical form at $`2t_0`$ is quite noticeable, but cannot be seen in the simulation data. The kink has its origin in the complete suppression of nucleation for negatives times, and the absence of the kink in the simulation data may stem from the finite chance of nucleation at negative times due to the way the field is reoriented, as described in Sec. III A. Another possibility is that the interactions between the ends smooth out the difference between the one- and two-nucleation decay modes. Despite the kink, the theoretical form Eq. (15) does a good job of describing both the exponential tail and the rounding at early times of $`P_{\mathrm{not}}(t)`$. It is clearly superior to the error-function form, as well as to a displaced, single exponential (not shown).

## V Arrays of nanomagnets

The long-ranged dipole-dipole interactions that contribute to the shape anisotropy of nanoscale magnets also cause interactions between nanomagnets. This is especially true in most potential applications, where miniaturization will drive devices to high densities. In addition, the interactions between nanomagnets in arrays could be the basis of device applications, prototypes of which have already been investigated.  Regular arrays of nanomagnets have already been used experimentally to provide magnetic signals strong enough to be measured  in the experiments our nanomagnets are modeled after. Our simulations show that even for very wide spacings the magnetic interactions between nanomagnets have significant effects on the switching properties.

To allow for analytic treatment, we consider only the two leading contributions to the local dipole field $`𝐇_\mathrm{d}(𝐫)`$ from the dipole and quadrupole moments of the source magnet. The former contributes uniformly throughout the observation volume, while the latter changes linearly in each Cartesian direction. Specifying the distance between pillars as $`d`$ and considering the dipolar and quadrupolar moments of the source magnet, $`M_1^0`$ and $`M_2^0`$, respectively, the leading contributions to the observed demagnetizing field are $`𝐇_\mathrm{d}^{(1)}(𝐫)=\widehat{𝐞}_\mathrm{z}M_1^0/d^3`$ and $`𝐇_\mathrm{d}^{(2)}(𝐫)=9(2z\widehat{𝐞}_\mathrm{z}x\widehat{𝐞}_\mathrm{x}y\widehat{𝐞}_\mathrm{y})M_2^0/(4d^5),`$ respectively. Simple expressions for the multipole moments can be calculated in the following way. Assume a pillar with regions of uniform magnetization oriented in the $`+\widehat{𝐞}`$ direction in the middle and in the $`\widehat{𝐞}`$ direction at the ends, as shown in Fig. 9. With the total pillar length $`L`$, the top-region length $`\mathrm{}_\mathrm{t}`$, and the bottom-region length $`\mathrm{}_\mathrm{b}`$ all measured along the long axis of the pillar, the dipole moment for a reversing pillar is

$$M_1^0=M_\mathrm{S}A\left(L2\mathrm{}_\mathrm{t}2\mathrm{}_\mathrm{b}\right)$$
(19)
and the quadrupole moment is

$$M_2^0=2M_\mathrm{S}A\left(\mathrm{}_\mathrm{t}^2\mathrm{}_\mathrm{b}^2L\left(\mathrm{}_\mathrm{t}\mathrm{}_\mathrm{b}\right)\right),$$
(20)
where $`A`$ is the cross-sectional area of the pillar. The dimensionless form is achieved by setting $`M_\mathrm{S}`$ to unity and using dimensionless lengths.

Using the simplified expression for the dipole moment in Eq. (19), the contribution from one upward-magnetized pillar at the nearest-neighbor position in the array is $`𝐇_\mathrm{d}^{(1)}\widehat{𝐞}_\mathrm{z}1`$ Oe. The quadrupole contribution from Eq. (20) at the end of one nanomagnet that is the nearest neighbor to a switching pillar with $`\mathrm{}_\mathrm{t}=L/2`$ is $`𝐇_\mathrm{d}^{(2)}\widehat{𝐞}_\mathrm{z}0.4`$ Oe. This is a rough estimate of the maximum interaction through the quadrupole moment, and the sign indicates a tendency for neighboring pillars in the array to switch at opposite ends.

As a simple initial investigation, we consider a linear array of four of the rectangular nanomagnets described in Sec. III, with the nanomagnets oriented perpendicular to the substrate. Since dipole-dipole interactions within each nanomagnet are calculated using the Fast Multipole Method described in Appendix 1, the multipole moments for each nanomagnet are readily available. These moments can be used to quickly calculate the interactions between nanomagnets in the array, under the constraint that the far-field description is appropriate. To ensure this we have considered only a spacing between pillars of two pillar lengths, $`300\mathrm{nm}.`$ Note that this situation is quite different from Fast Fourier Transform approaches, in which the calculation must be carried out on a lattice that also fills the space between the nanomagnets.  To do that practically, the space between the magnets must be kept small.

To study arrays of nanomagnets, systems consisting of four $`9`$ nm $`\times `$ $`9`$ nm $`\times `$ $`150`$ nm parallel pillars arranged in a line perpendicular to their long axes and spaced $`300`$ nm apart were simulated using the Fast Multipole Method truncated at $`p=3.`$ Hysteresis loops with periods on the order of a few nanoseconds for the individual pillars in the array look similar to those for isolated pillars in Fig. 3, with no observable difference between pillars on the outside and those on the inside of the array. The symmetry equivalence for the two pillars on the inside, and for the two pillars on the outside, will be used throughout to double the statistical sampling.

The probability of not switching for $`H=1800`$ Oe and $`T=20`$ K is shown for $`40`$ array switches in Fig. 10(a) for pillars on the inside and outside of the array, as well as isolated pillars. No significant difference can be seen between the three curves. However, the coupling between the pillars can be seen in the difference between the $`P_{\mathrm{not}}(t)`$ for inside pillars with one or both nearest neighbor pillars switched, shown in Fig. 10(b). Here, $`t`$ is the time difference between $`t_{\mathrm{sw}}`$ and the last time a neighboring pillar switched. From this data it can be seen that of pillars with two neighbors, those with only one neighbor switched tend to switch earlier than those with both pillars switched. The effect is even more pronounced in simulations of the simple model at $`H`$$`=`$$`1000\mathrm{Oe}.`$

## VI Summary

Numerical simulation of the Landau-Lifshitz-Gilbert micromagnetic model has been used to investigate spatially non-uniform magnetization switching in nanoscale magnets at the nanosecond time scale. We have focused on iron pillars $`9\mathrm{nm}\times 9\mathrm{nm}\times 150\mathrm{nm}`$ because such pillars and arrays have been constructed and measured experimentally.  The zero-temperature static coercive field has been estimated numerically to be $`H_c=1979\pm 14\mathrm{Oe}`$ by finding the field of maximum energy for fields swept at constant rate, and then extrapolating to find the zero-rate estimate. Simulations of thermally activated magnetization switching are possible at fields well below the coercive value. For the pillars studied here, reversal occurs through nucleation at the ends of the pillars. The probability of not switching, $`P_{\mathrm{not}}(t),`$ is not well described by a delayed exponential, and a theoretical form, Eq. (15), based on independent nucleation at the ends of the magnet is developed here. The agreement with results for intensive, fully three-dimensional simulations with an applied field near the coercive value are reasonable, but an ad hoc error function gives similar agreement. The agreement with Eq. (15) is much better when the field is well below the coercive value and the statistics are better, which currently we have only studied with less-intensive simulations that discretize the nanomagnet only along its long axis.

The fully three-dimensional micromagnetics program has been developed for massively parallel computers and implemented using the Fast Multipole Method. These two features make it feasible to simulate nanomagnets in widely spaced arrays. Even though the interactions between the magnets are quite weak for the specific linear array considered here, there are significant effects on the statistics of the magnetization switching. Specifically, there is a dependence of $`P_{\mathrm{not}}(t)`$ for a given pillar on the orientation of the magnetization of its nearest neighbors in the array. The nature of the cooperative reversal mode observed here is a topic for future research.

## Acknowledgments

This work was supported by NSF Grant No. DMR-9871455, Florida State University (FSU) Center for Materials Research and Technology, the FSU Supercomputer Computations Research Institute (U.S. DOE contract number DE-FC05-85ER25000), and the FSU School of Computational Science and Information Technology. Extensive supercomputer resources were provided by FSU Academic Computing and Network Services and by the U.S. Department of Energy through the U.S. National Energy Research Scientific Computing Center.

## Appendix 1

Calculating the dipole-dipole interactions is the most intensive part of the numerical calculation. The magnetic potential approach used here involves defining a magnetic charge density $`\rho _\mathrm{M}(𝐫)=𝐌(𝐫).`$ (We present this only in terms of the dimensionless quantities.) This charge was evaluated on a cubic lattice dual to that of $`𝐌(𝐫_i),`$ see Fig. 11, using equally weighted two-point differences, specifically

$$\rho _\mathrm{M}(𝐫_\mathrm{d})=\underset{\mu =1}{\overset{3}{}}\left(\frac{1}{4\mathrm{\Delta }r}\underset{|𝐡|=\sqrt{3}\mathrm{\Delta }r/2}{}\mathrm{sign}\left(𝐡\widehat{𝐞}_\mu \right)𝐌\left(𝐫_d+𝐡\right)\widehat{𝐞}_\mu \right),$$
(21)
where $`𝐫_\mathrm{d}`$ is a position on the dual lattice, $`_𝐡`$ is the sum over the corresponding corners of the cube from the direct lattice, and $`\widehat{𝐞}_\mu `$ are the Cartesian unit vectors. If we were to apply Eq. (21) with $`𝐫_\mathrm{d}`$ just outside the magnetic material, it would give a nonzero charge there. To avoid this unphysical result, we have moved this charge to the surface by defining a surface magnetic charge density $`\sigma _\mathrm{M}(𝐫)=\widehat{𝐬}(𝐫)𝐌(𝐫),`$ where $`\widehat{𝐬}`$ is the unit vector directed out of the surface. We have considered surface charges only on the surface of the model magnet, and we evaluate them at the centers of the squares defined by adjacent points on the surface of the direct lattice. The four corners are equally weighted so that

$$\sigma _\mathrm{M}(𝐫_\mathrm{d}^{})=\frac{1}{4}\underset{|𝐡|=\sqrt{2}\mathrm{\Delta }r/2}{}\widehat{𝐬}𝐌(𝐫_\mathrm{d}^{}+𝐡),$$
(22)
where $`_𝐡`$ runs over the four corners. The numerical approach of Eqs. (21) and (22) ensures that there is no net magnetic charge on the system as a whole.

The magnetic potential $`\varphi _\mathrm{M}(𝐫)`$ is found by integrating over both the volume and surface charges,

$$\varphi _\mathrm{M}(𝐫)=_\mathrm{V}𝑑𝐫^{}\frac{\rho _\mathrm{M}(𝐫^{})}{|𝐫𝐫^{}|}+_\mathrm{S}𝑑\widehat{𝐬}^{}\frac{\sigma _\mathrm{M}(𝐫^{})}{|𝐫𝐫^{}|}.$$
(23)
Numerically, such an operation can be quite expensive since unsophisticated algorithms will require $`O(N^2)`$ operations, where $`N`$ is the number of lattice sites. An algorithm that remains reasonable for large systems is the Fast Multipole Method.

We have chosen to calculate the magnetic potential using the Fast Multipole Method (FMM) because it has several advantages over the more traditional Fast Fourier Transform (FFT) approach. The biggest difference between the two methods is that the FMM makes no assumptions about the underlying lattice, while the FFT method assumes a cubic lattice with periodic boundary conditions. One consequence of this assumption is that numerical models of systems without periodic boundary conditions require empty space around the magnet so that the boundary conditions do not affect the calculation. The FFT also requires the lattice to continue into regions of empty space that lie between elements of an array of magnets. By contrast, a FMM implementation only needs to consider volumes occupied by magnetic material. It does not need any padding. In addition to these advantages, the FMM is more efficient for large numbers of lattice sites.

The popularity of the FFT approach stems from the fact that it takes $`O(N\mathrm{ln}N)`$ calculations to evaluate Eq. (23). The FMM has a larger overhead, but requires only $`O(N)`$ operations to calculate the same potential.  This means an FFT approach to Eq. (23) makes sense for small, cubic lattices, but that the FFM approach will be more efficient for large, irregular, or incomplete lattices.

The FMM algorithm exploits the fact that $`\varphi _\mathrm{M}`$ at each lattice point can be expanded in terms of spherical harmonics,

$$\varphi _\mathrm{M}(r,\theta ,\phi )=\underset{i=0}{\overset{\mathrm{}}{}}\underset{j=i}{\overset{i}{}}\left(L_i^jr^i+\frac{M_i^j}{r^{i+1}}\right)Y_i^j(\theta ,\phi ),$$
(24)
where the $`L_i^j`$ terms can be used to represent the potential close to the lattice point, and the $`M_i^j`$ terms can be used to represent the potential far from it, but not both simultaneously. In this context, near $`𝐫_i`$ means being closer to $`𝐫_i`$ than any other lattice point, and far means distances more than twice the largest distance from the center of cell $`j`$ to any of its boundary points. Following Greengard,  we define the spherical harmonics

$$Y_i^j(\theta ,\phi )=\sqrt{\frac{(i|j|)!}{(i+|j|)!}}P_i^{|j|}(\mathrm{cos}\theta )e^{Ij\phi },$$
(25)
with $`P_i^j(x)`$ the associated Legendre polynomial, and $`I=\sqrt{1}`$. Actually implementing this approach requires a truncation of the expansion in $`i`$ at order $`p`$. We have found that the demagnetizing field for $`p=3`$ is within $`1\%`$ of the exact value for our simulations.

Our implementation of the FMM algorithm starts by partitioning the model space into a system of cubes, an example is shown in Fig. 12. The length of the side of the smallest cubes, which are separated by dotted lines, is the same as the discretization length. For our lattice system we have found that the most efficient choice is to use cubes centered on the direct lattice. For each cube, the multipole expansion coefficients of the far-field potential are calculated for the specific configuration of magnetic charges $`\rho _\mathrm{M}`$ and $`\sigma _\mathrm{M}`$

$$M_i^j=_V𝑑𝐫\left[\rho _\mathrm{M}+\sigma _\mathrm{M}\right]\varrho ^iY_i^j(\theta ,\phi ),$$
(26)
with the coordinates centered on the lattice site. Note here that $`\varrho ^i`$ is the distance from the center of the cell raised to the $`i`$-th power. For our cubic lattice, each quadrant of the cube is contributed by a different region of constant $`\rho _\mathrm{M}`$ from the dual lattice (similarly for $`\sigma _\mathrm{M}`$). With the geometry of the lattice fixed, the multipole expansion coefficients are easily calculated and summed to yield the total expansion coefficients for the lattice site.

Each level of the hierarchy involves grouping cells into successively larger cubes that completely contain cubes at the lower level. The obvious hierarchy with each larger cube containing eight of the smaller cubes apparently works best. In Fig. 12, cubes of the second level are separated by dashed lines and those of the third level by solid lines. Since the number of nodes in any direction along the simulation lattice is not restricted to a power of two, cubes do not always contain eight smaller cubes. The $`M_i^j`$ for the larger cube can be rapidly evaluated from those of the smaller cubes using the rule for translation of a multipole expansion

$$M_i^j=\underset{k=0}{\overset{i}{}}\underset{l=k}{\overset{k}{}}O_{ik}^{jl}\frac{J_l^{jl}A_k^lA_{ik}^{jl}}{A_i^j}\varrho ^kY_k^l(\theta ,\phi ),$$
(27)
where $`O_i^j`$ are the expansion coefficients for the smaller cube and the spherical coordinates $`(\varrho ,\theta ,\phi )`$ here are for the vector from the origin of the large cube to that of the small cube. The relations for the new factors are

$$A_i^j=\frac{(1)^j}{\sqrt{(ij)!(i+j)!}}$$
(28)
and

$$J_i^j=\{\begin{array}{cc}(1)^{\mathrm{min}(|i|,|j|)}\hfill & \hfill \mathrm{if}ij<0\\ 1\hfill & \hfill \mathrm{otherwise}.\end{array}$$
(29)
The construction of the hierarchy terminates when the entire system is enclosed by a single cube. The partitioning can be generalized to noncubic rectangular prisms, but the restriction that the multipole expansion is only valid at distances larger than twice the diagonal of the rectangle will complicate the algorithm.

A downward pass through the hierarchy involving three types of operations is required to construct $`\varphi _\mathrm{M}`$ as represented by the local expansion coefficients $`L_i^j`$ for each cube. The first operation is the translation of the local expansion for the encompassing bigger cube (if one exists) that includes all the contributions from its far field, i.e., those cubes on its level that are not its neighbors. This translation is accomplished by the rule

$$L_i^j=\underset{k=i}{\overset{p}{}}\underset{l=k}{\overset{k}{}}O_{k,l}\frac{\overline{J}_{ki,lj}^lA_{ki}^{lj}A_i^j}{A_k^l}\varrho ^{ki}Y_{k1}^{lj}(\theta ,\phi ),$$
(30)
where $`O_i^j`$ are the local expansion coefficients of the larger cube, and the spherical coordinates $`(\varrho ,\theta ,\phi )`$ here are for the vector from the origin of the smaller cube to that of the larger cube. $`A_i^j`$ is defined in Eq. (28), while

$$\overline{J}_{i,j}^k=\{\begin{array}{cc}(1)^i(1)^j\hfill & \hfill \mathrm{if}jk<0\\ (1)^i(1)^{kj}\hfill & \hfill \mathrm{if}jk>0\mathrm{and}|k|<|j|\\ (1)^i\hfill & \hfill \mathrm{otherwise}.\end{array}$$
(31)
The next operation incorporates the contributions from areas in the near region of the larger cube, but in the far region of the smaller cube. This is accomplished by transforming the multipole expansion of the source into a local expansion of the smaller cube using

$$L_i^j=\underset{k=0}{\overset{p}{}}\underset{l=k}{\overset{k}{}}O_k^l\frac{\stackrel{~}{J}_k^{j,l}A_k^lA_i^j}{A_{i+k}^{lj}}\frac{Y_{i+k}^{lj}(\theta ,\phi )}{\varrho ^{i+k+1}},$$
(32)
where the spherical coordinates are for the vector from the origin of the smaller cube to that of the larger cube and

$$\stackrel{~}{J}_i^{j,k}=\{\begin{array}{cc}(1)^i(1)^{\mathrm{min}(|j|,|k|)}\hfill & \hfill \mathrm{if}jk>0\\ (1)^i\hfill & \hfill \mathrm{otherwise}.\end{array}$$
(33)
The third type of operation is used for termination of the algorithm at the lowest level, where the near-region contributions for the smallest cubes must be evaluated exactly. Since a regular lattice is used, the point $`𝐫_i`$ where $`\varphi _\mathrm{M}`$ is being calculated always lies at the corner of the neighboring cube of constant $`\rho _\mathrm{M}`$, and the contribution is simply

$$\frac{\left(\mathrm{\Delta }r\right)^2}{8}\left(6\mathrm{sinh}^1\left(\frac{1}{\sqrt{2}}\right)\frac{\pi }{2}\right)\rho _\mathrm{M}\left(𝐫_d\right),$$
(34)
because only one quadrant of the region of constant $`\rho _\mathrm{M}`$ needs to be considered in this exact manner. The contribution for squares of surface charge that touch $`𝐫_i`$ is

$$\mathrm{\Delta }r\mathrm{sinh}^1\left(1\right)\sigma _\mathrm{M}\left(𝐫_d\right).$$
(35)
Similar expressions can be calculated for general rectangular prisms.

The FMM is an efficient way of determining the magnetic potential $`\varphi _\mathrm{M}(𝐫_i)`$ associated with a particular configuration of a dipole field. The local observed magnetic field due to the other dipoles is obtained using

$$𝐇_\mathrm{d}(𝐫_i)=\varphi _\mathrm{M}(𝐫_i)$$
(36)
from potential theory. Numerically, this gradient was estimated with a centered difference of the nearest-neighbor sites, except on the boundaries, where forward or backward differencing was used. Note that the discretizations of all operators have been chosen for consistency between surface and volume charges such that the model magnet can be padded by volumes of lattice sites with $`𝐌=0`$ without changing the results.

The FMM was implemented using C++ and named Hierarchy.h. Fundamental to the implementation is the class sh\_expansion whose instances are a $`(p+1)^2`$ array of complex numbers representing either local or multipole expansion coefficients. The class includes methods for indexing expansion coefficients within an expansion, evaluating the expansion at a specified point, and transforming an expansion using Eqs. (27), (30), and (32). In addition, the class encapsulates the precomputation of quantities that do not depend on the displacement vector, and it can return a pointer to a table of precomputations that do depend on a fixed displacement. For instance, for efficient evaluation of the translation of multipole expansions, Eq. (27), the class contains the $`(p+1)^4`$ values

$$C_{i,k}^{j,l}=\frac{J_l^{jl}A_k^lA_{ik}^{jl}}{A_i^j}$$
(37)
and an indirect-indexing array that specifies the sequence of $`O_{k,l}`$ after the loops have been unrolled. In addition, a pointer to the $`(p+1)^2`$ precomputed values of $`\varrho ^kY_k^l(\theta ,\phi )`$ can also be returned. Thus the transformation can be evaluated efficiently and with minimal overhead from loop-control variables.

The computer memory requirements to store the precomputations that depend on spatial relationships can be greatly reduced for implementations that assume a regular lattice for spatial decomposition. For each level of the hierarchy, there will be a fixed number of displacements between the cells in that level and with cells on the parent level and the child level. The class sh\_expansion encapsulates this efficiency at each request for a precomputation by searching a linked list of previous results and creating a new result only when no previous result exists. A similar scheme is also necessary for the compilation of multipole expansion coefficients for the lowest level of the hierarchy, Eq. (26). Since the precomputations represent a significant amount of the memory requirements of the overall simulation, these memory savings can greatly increase the number of lattice points used to represent the model magnet.

The partitioning of space is accomplished through the class hbox. Instances of this class contain geometric information about the decomposition cell, expansions for both the potential and charges within the cell, as well as links to other cells. These links point to the encompassing cell, the encompassed cells, a list of same-level cells in the far field (the “interacting cells”), and a list of nearest-neighbor cells. The hboxs for each level of the decomposition hierarchy are held in a container class hlevel, and the entire hierarchy is maintained as a linked list. Our implementation of hlevel is particular to the cubic-lattice decomposition of space, and irregular geometries have to be padded with empty space. (This is also true for our implementation of the classes VectorField and ScalarField from VectorField.h used throughout our numerical integration of the LLG equation.)

## Appendix 2

A useful way to represent the thermal Landau-Lifshitz-Gilbert equation is in a form with the deterministic and stochastic parts separated

$$d𝐌(𝐫_i,t)=𝐁\left(𝐌(𝐫_i,t)\right)𝐇_{\mathrm{det}}(𝐫_i,t)dt+\sqrt{ϵ}𝐁\left(𝐌(𝐫_i,t)\right)d𝐖(𝐫_i,t),$$
(38)
where $`𝐇_{\mathrm{det}}(𝐫_i,t)`$ is the deterministic part of the local field at $`𝐫_i`$ and $`𝐇_n\left(𝐫_i\right)=\sqrt{ϵ}d𝐖(𝐫_i,t)`$ is the stochastic part. The matrix $`𝐁`$ is given by

$$𝐁\left(𝐌(𝐫_i,t)\right)=\frac{1}{1+\alpha ^2}\left(\begin{array}{ccc}\alpha \left(M_y^2+M_z^2\right)& M_z\alpha M_xM_y& M_y\alpha M_xM_z\\ M_z\alpha M_xM_y& \alpha \left(M_x^2+M_z^2\right)& M_x\alpha M_yM_z\\ M_y\alpha M_xM_z& M_x\alpha M_yM_z& \alpha \left(M_x^2+M_y^2\right)\end{array}\right),$$
(39)
where $`x`$, $`y`$, and $`z`$ represent the Cartesian coordinates and the space and the time dependence of the $`M_\mu `$ have been omitted for clarity. The stochastic nature of the field results from the Wiener  process $`𝐖(𝐫_i,t)`$ which has the properties

$$W_\mu (𝐫_i,t)=0W_\mu (𝐫_i,t)W_\mu ^{}(𝐫_i^{},t^{})=(tt^{})\delta _{\mu ,\mu ^{}}\delta _{i,i^{}}.$$
(40)
The stochastic differential equation (38) can be treated numerically using first-order Euler integration. For small $`\mathrm{\Delta }t`$, the deterministic integral is

$$𝐈_{\mathrm{det}}(𝐫_i,t)=𝐁\left(𝐌(𝐫_i,t)\right)𝐇_{\mathrm{det}}(𝐫_i,t)\mathrm{\Delta }t.$$
(41)
The integral of the stochastic part,

$$𝐈_{\mathrm{sto}}(𝐫_i,t)=\sqrt{ϵ}_t^{t+\mathrm{\Delta }t}𝐁\left(𝐌(𝐫_i,t^{})\right)𝑑𝐖(𝐫_i,t^{}),$$
(42)
takes more consideration since it involves the product of the magnetization with the Wiener process. In such cases of multiplicative noise, different methods for evaluating Eq. (42) correspond to different Fokker-Planck equations.  There are an infinite number of ways to interpret Eq. (42), but usually only the two extreme cases, the Itô and Stratonovich interpretations, are considered. The Fokker-Planck equation considered by Brown  only has the proper equilibrium properties when interpreted in the Stratonovich sense.  This is complicated by the fact that numerical implementation of the stochastic integral is particularly convenient in the Itô interpretation. Then the discretized integral is

$$𝐈_{\mathrm{sto}}(𝐫_i,t)=\sqrt{ϵ\mathrm{\Delta }t}𝐁\left(𝐌(𝐫_i,t)\right)𝐠(𝐫_i,t),$$
(43)
where each component of $`𝐠`$ is a random number from a Gaussian distribution with zero mean and unit variance. Fortunately, changing interpretations can be accomplished through a correction discussed below. Then, Eq. (43) can be combined with Eq. (41), and the result rearranged, to give the “impulse” during the integration step, Eq. (11).

Changing from the Stratonovich interpretation of the stochastic integral to that of Itô requires the addition of a deterministic term. Specifically, a Stratonovich interpretation of a multivariate-Langevin equation of the form

$$d𝐱=𝐚(𝐱,t)dt+𝐛(𝐱,t)d𝐖,$$
(44)
is equivalent to the Langevin equation

$$d𝐱=\left(𝐚(𝐱,t)+𝐜(𝐱,t)\right)dt+𝐛(𝐱,t)d𝐖$$
(45)
in the Itô interpretation  where the components of the new drift term are

$$c_\mu =\frac{1}{2}\underset{\nu ,\nu ^{}}{}b_{\nu ^{},\nu }\frac{b_{\mu ,\nu }}{\nu ^{}}.$$
(46)
For the present system, the additional drift term is readily found to be

$$𝐜(𝐫_i,t)=\frac{ϵ}{(1+\alpha ^2)}𝐌(𝐫_i,t),$$
(47)
which is equivalent to the result given in Ref. . This drift term is always directed along the local magnetization density. The process of normalizing the magnitude of the spins during each integration step, Eq. (10), is also directed along the magnetization, and essentially takes the correction into account.

| $`\mathrm{}_e`$ | $`3.6\mathrm{nm}`$ |
| --- | --- |
| $`M_s`$ | $`1700\mathrm{emu}/\mathrm{cm}^3`$ |
| $`\mathrm{\Delta }\stackrel{~}{𝐫}`$ | $`1.5\mathrm{nm}`$ |
| $`\mathrm{\Delta }\stackrel{~}{t}`$ | $`5.0\times 10^5\mathrm{ns}`$ |
| $`\gamma _0`$ | $`1.76\times 10^7\mathrm{Hz}/\mathrm{Oe}`$ |
| $`\alpha `$ | $`0.1`$ |

Table 1. Parameters used in the micromagnetic simulations described here.


text:
# Topological Quantum Computation

## Implementation of a quantum computer

It is not possible to realize a unitary gate precisely, and even when we do not intend to apply a gate, the environment will interact with the qubits causing decoherence. Both imprecision and decoherence can be considered in a unified way as $`\mathrm{`}\mathrm{`}`$errors” which can be quantified by a fidelity distance $`(\mathit{17})`$ or a super-operator norm $`(\mathit{19})`$. A crucial step in the theory of quantum computing has been the discovery of error-correcting quantum codes $`(\mathit{28})`$ and fault-tolerant quantum computation $`(\mathit{25},\mathit{29})`$. These techniques cope with sufficiently small errors. However, the error magnitude must be smaller than some constant (called an accuracy threshold) for these methods to work. According to rather optimistic estimates, this constant lies between $`10^5`$ and $`10^3`$, beyond the reach of current technologies.

The presumptive path toward a quantum computer includes these steps: (1) build physical qubits and physical gates; (2) minimize error level down below the threshold; (3) implement decoherence-protected logical qubits and fault-tolerant logical gates (one logical qubit is realized by several qubits using an error-correcting code.) As a counterpoint, the theme of this article is that implementing physical qubits might be redundant. Indeed, one can $`\mathrm{`}\mathrm{`}`$encode” a logical qubit into a physical system which is not quite a set of qubits or even well separated into subsystems. Such a system must have macroscopic quantum degrees of freedom which would be decoherence-protected. A super-conducting phase or anything related to a local order parameter will not work for if a degree of freedom is accessible by local measurement, it is vulnerable to local error. However, there is another type of macroscopic quantum degree of freedom. It is related to topology and arises in collective electronic systems, e.g. the fractional quantum Hall effect $`(\mathit{13})`$ and most recently in $`2D`$ cuprate superconductors above $`T_c`$ $`(\mathit{12},\mathit{23})`$.

Though much studied since the mid-1980’s the connection between fractional quantum Hall effect and quantum computation has only recently been realized $`(\mathit{11},\mathit{20})`$. It was shown by $`(\mathit{4})`$ that the ground state of the $`\nu =\frac{1}{3}`$ electron liquid on the torus is 3-fold degenerate. This follows from the fact that excitations in this system are abelian anyons: moving one excitation around another multiplies the state vector by a phase factor $`e^{i\varphi }`$ (in this case $`\varphi =\frac{2\pi }{3}`$). The process of creating a particle-antiparticle pair, moving one of the particles around the torus, and annihilating it specifies a unitary operator on the ground state. By moving the particle in two different directions, one obtains two different unitary operators $`A_1`$ and $`A_2`$ with the commutation relation $`A_1A_2A_1^1A_2^1=e^{i\varphi }`$, implying a ground state degeneracy. This argument is very robust and only requires the existence of an energy gap or, equivalently, finite correlation length $`l_0`$. Indeed, the degeneracy is lifted only by spontaneous tunneling of virtual excitations around the torus. The resulting energy splitting scales as $`e^{\frac{l}{l_0}}`$, where $`l`$ is the size of the system. Interaction with the environment does not change this conclusion, although thermal noise can create actual excitation pairs rather than virtual ones.

Both the ground state degeneracy on the torus and the existence of anyons are manifestations of somewhat mysterious topological properties of the $`\nu =\frac{1}{3}`$ electron liquid itself. Anyons can be regarded as $`\mathrm{`}\mathrm{`}`$topological defects” similar to Abrikosov vortices but without any local order parameter. The presence of a particle enclosed by a loop on the plane can be detected by holonomy—moving another particle around the loop.

At $`\nu =\frac{1}{3}`$, the electron liquid on the torus could be used as a logical $`\mathrm{`}\mathrm{`}`$qutrit” (generalized qubit with 3 states). Unfortunately, this will hardly work in practice. Besides the obvious problem with implementing the torus topology, there is no known way to measure this logical qutrit or prepare it in a pure state.

A more flexible and controllable way of storing quantum information is based on nonabelian anyons. These are believed to exist in the $`\nu =\frac{5}{2}`$ fractional quantum Hall state. According to the theory $`(\mathit{22},\mathit{24})`$, there should be charge $`\frac{1}{4}`$ anyonic particles and some other excitations. The quantum state of the system with $`2n`$ charge $`\frac{1}{4}`$ particles on the plane is $`2^{n1}`$-degenerate. The degeneracy is gradually lifted as two particles come close to each other. More precisely, the $`2^{n1}`$-dimensional Hilbert space $`_n`$ splits into two $`2^{n2}`$-dimensional subspaces. They correspond to two different types of charge $`\frac{1}{2}`$ particles which can result from fusion of charge $`\frac{1}{4}`$ particles. Thus observing the fusion outcome effects a measurement on the Hilbert space $`_n`$. This model supports adaptive quantum computation when surfaces of high genus are included in the theory and admits a combinatorial description $`(\mathit{1})`$ apparently in the same universal class as the fractional quantum Hall fluid.

Beyond this, a discrete family of quantum Hall models exists $`(\mathit{26})`$ based on $`k+1`$-fold hard-core interaction between electrons in a fixed Laudau level which appears to represent the same universality class as Witten-Chern-Simons theory for $`SU(2)`$ at level=$`k`$ $`(\mathit{33})`$. Anyons in these models behave as topological defects of a geometric construction $`(\mathit{10})`$ and their braiding matrices have been shown to be universal $`(\mathit{8},\mathit{9})`$ for $`k3,k4`$.

## Code spaces and quantum media

Even after the particle types and positions of anyons are specified, there is an exponentially large (but finite dimensional) Hilbert space describing topological degrees of freedom. In combinatorial models, this Hilbert space becomes a $`\mathrm{`}\mathrm{`}`$code subspace” $`W`$ of a larger $`\mathrm{`}\mathrm{`}`$quantum media” $`Y`$. Thus a fundamental concept of cryptography is transplanted into physics. Let $`V`$ be a finite dimensional complex vector space, perhaps $`^2`$, and $`Y=V\mathrm{}V`$ an $`n`$-fold tensor product (where $`n`$ is typically quite large.) Let $`WY`$ be a linear subspace. We call $`WY`$ $`k`$-code if and only if:

$$W\stackrel{\mathrm{\Pi }_W𝒪}{}W$$

is multiplication by a scalar whenever $`𝒪`$ is a $`k`$-local operator (an arbitrary linear map on any $`k`$-tensor factors of $`Y`$ and the identity on the remaining $`nk`$ factors ) and $`\mathrm{\Pi }_W`$ is the orthogonal projection onto $`W`$. We think of such spaces as resisting local alteration and in the usual interpretation of $`Y`$ as the Hilbert space of $`n`$ particles, quantum information stored in $`W`$ will be relatively secure. It is a theorem $`(\mathit{14})`$ that the quantum information in $`W`$ cannot be degraded by errors operating on fewer than $`\frac{k}{2}`$ of the $`n`$ particles.

Let us define a (discrete) quantum medium to be a tensor product $`Y=_iV_i`$ as above, where now the set of indices $`\{i\}`$ consists of points distributed on a geometric surface $`T`$, together with a local Hamiltonian $`H=H_i`$ (each $`H_i`$ is a Hermitian operator defined only on those tensor factors whose index is within $`ϵ0`$ of the $`i`$-th point in the geometry of the surface.) Local Hamiltonians $`H`$ have been found $`(\mathit{10},\mathit{20})`$ with highly $`d`$-degenerate ground states corresponding to modular functors $`(\mathit{31},\mathit{32})`$, (and thus braid group representation $`(\mathit{16})`$ and link polynomials $`(\mathit{18})`$.) In these cases, the ground state $`G`$ of $`H`$ will be $`k`$-code for $`k\text{injectivity radius of T}\sqrt{\text{area T}}`$. The topological degrees of freedom referred to above reside in $`G`$. But we do not attempt a precise mathematical definition of topological degrees of freedom since we would like it to extend beyond discrete system, e.g., to fractional quantum Hall ground states.

In the case when $`T`$ is a disk $`D`$ with punctures—physically anyonic excitations—a sequence of local modifications to $`H`$ (see $`(\mathit{10})`$) effects a discrete 1-parameter family $`H_t`$ of Hamiltonians, where the ground states $`G_{t_i}`$ and $`G_{t_{i+1}}`$ at consecutive time steps differ by a $`\frac{k}{2}`$-local operator $`𝒪_i`$, $`𝒪_i(G_{t_i})=G_{t_{i+1}}`$. Note that if $`G_{t_i}`$ and $`G_{t_{i+1}}`$ are both $`k`$-code that for $`𝒪_i`$ as above, $`𝒪_i|_{G_{t_i}}`$ must be unique up to a scalar (for a distinct $`𝒪_i^{}`$, consider the restriction of the $`k`$-local operator $`𝒪_{i}^{}{}_{}{}^{}𝒪_i^{^{}}`$ projected to $`G_{t_i}`$, $`\mathrm{\Pi }_{G_{t_i}}𝒪_{i}^{}{}_{}{}^{}𝒪_i^{^{}}`$). This uniqueness property forces this discrete-cryptographic transport of code spaces to coincide (up to phase) with the differential geometric notion of adiabatic transport—integration of the canonical connection in the $`\mathrm{`}\mathrm{`}`$tautological” bundle of $`d`$-planes in $`Y=^{2^n}`$. If the 1-parameter family is a closed loop, a projective representation of the braid group on $`𝕌(\text{code})`$ is obtained. By choosing $`H`$, these can be engineered to be precisely the Hecke algebra representations $`\{\rho _\lambda \}`$ associated to the Jones polynomial. From $`H`$ one can build a concrete model of quantum computation; the model and its connection to the Jones polynomial are described below. Although the $`H`$ found in $`(\mathit{10})`$ is enormously cumbersome, it appears to share the universality class $`\mathrm{`}\mathrm{`}`$Witten-Chern-Simons theory of $`SU(2)`$ at level $`3`$” with a simple $`4`$-body Hamiltonian $`(\mathit{26})`$, which has been proposed as a model for the fractional quantum Hall plateaus at $`\nu =\frac{13}{5}\text{and}\frac{12}{5}`$.

## An anyonic model for quantum computation

A family of unitary representations of all mapping class and braid groups with certain compatibility properties under fusion is known as a unitary topological modular functor. To define our model, we take only the planar surface portion of the simplest universal modular functor, Witten-Chern-Simons $`SU(2)`$ modular functor at level 3, and this reduces to the Jones representations of braids $`\{\rho _\lambda \}`$. For an appropriate family of local Hamiltonians $`H_t`$, these representations describe the adiabatic transport of the lowest energy states with $`n`$ anyonic excitation pairs, these states form a subspace $`W`$ of dimension Fibonacci(2n)— as the $`2n`$ anyonic excitations are $`\mathrm{`}\mathrm{`}`$braided” around each other in 2+1 dimensional space time. In this theory, (we denote it $`CS5`$ because of its link to the Jones polynomial at a fifth root of unity), there are 4 label types 0,1,2,3 corresponding to the complete list of irreducible representations of the quantum group $`SU(2)_5`$ of dimensions 1,2,3 and 4 (or equivalently the irreducible positive energy loop group representations of $`SU(2)`$ at level $`3`$). We initialize our system on the disk $`D`$ in a known state by pulling anyonic pairs out of the vacuum. This theory is self-dual so the two partners have identical types. Pairs are kept or returned to the vacuum according to the results of local holonomic measurement. Finally we have a known initial state in the disk with $`2n`$ punctures where each puncture has label=1 and $`D`$ has label=0. We assume $`n`$ is even and group the punctures into $`n/2`$ batches of 4 punctures each. Similar to $`(\mathit{8})`$, each batch $`B`$ is used to encode one qubit $`^2`$: the basis $`\{|0,|1\}`$ is mapped into the type ($`0`$ or $`2`$) of the 2-fold composite particle (round circle), which would result from fusing a (fixed) pair of the type 1 particles within $`B`$. Initially, both the double and 4-fold composites (ovals) have type 0. By maintaining, at least approximately, this condition on the ovals after the braiding is complete, we define a $`\mathrm{`}\mathrm{`}`$computational summand” of the modular functor: it is spanned by $`n`$-bit strings of $`0`$’s and $`2`$’s residing on the 2-fold composites (round circles).

Figure 1

Now as in the quantum circuit model, a classical poly-time algorithm looks at the problem instance $`(F,y)`$ and builds a sequence of $`\mathrm{`}\mathrm{`}`$gates”, but this time the gates are braid generators (right half twist between adjacent anyons) $`\sigma _i,\mathrm{\hspace{0.33em}1}i2n1`$, and a powerful approximation theorem $`(\mathit{19},\mathit{30})`$ is used to select the braid sequence which approximates the more traditional quantum circuit solving $`(F,y)`$. So the topological model may be described as:

1. Initialization of a known state in the modular functor.
2. Classical computation of braid $`b`$ effecting a desired unitary transformation $`X`$ of the computational subspace of the modular fucntor.
3. Adiabatic implementation of the braid by (somehow) moving the anyons in $`D`$ to draw $`b`$ in space-time. (Here we keep the anyons separated by a scale $`l`$.)
4. Application of a projection operator $`\mathrm{\Pi }`$ to measure the type ($`0`$ or $`2`$) of the $`\mathrm{`}\mathrm{`}`$left most” composite particle (as seen in Figure 1.).

The last step is the direct analogue of measuring the first qubit in the quantum circuit model and the same formula (1) applies: the probability of observing type 0 (the null particle) is:

$$\text{prob}(0)=0|X^{}\mathrm{\Pi }X|0,$$

where the $`0`$’s on the right hand side represent our carefully prepared initial state with $`2n`$ type 1 excitations.

To close the topological discussion, we note that the previous formula can be translated (using the $`S`$ matrix) into a plat closure (See Figure 2) of a braid $`L=\text{plat}(b^1\gamma b)`$, where $`\gamma `$ is a small loop inserted to measure the left-most qubit, and now the outcome of the quantum circuit calculation, $`\text{prob}(0)`$, becomes a Jones evaluation

Figure 2

$$\text{prob}(0)=\frac{1}{1+[2]_5^2}\left(1+\frac{(1)^{c(L)+w(L)}(a)^{3w(L)}V_L(e^{2\pi i/5})}{[2]_5^{m(L)2}}\right);$$

where $`[2]_5=\frac{1+\sqrt{5}}{2}`$; $`c,m`$, and $`w`$ are the number of components, number of local minima, and writhe of $`L`$ respectively; and $`a=e^{\pi i/10}`$. Details of this calculation will be posted at url: www.tqc.iu.edu.

The braiding disturbs and reforms composite particle types with sufficient subtlety to effect universal computation. To reduce this model to engineering, very significant obstacles must be overcome: stable quantum media must be maintained in a suitable phase, e.g., $`CS5`$; excitations must be readily manipulated, and electrical neutral particle types 0 and 2 must be distinguished, presumably by holonomy experiments. Although these challenges are daunting, they are, perhaps, less difficult than a head-on assault on the accuracy threshold in the quantum circuit model.