MERRILL: Micromagnetic Earth Related Robust Interpreted Language Laboratory

Abstract Complex magnetic domain structures and the energy barriers between them are responsible for pseudo‐single‐domain phenomena in rock magnetism and contribute significantly to the magnetic remanence of paleomagnetic samples. This article introduces MERRILL, an open source software package for three‐dimensional micromagnetics optimized and designed for the calculation of such complex structures. MERRILL has a simple scripting user interface that requires little computational knowledge to use but provides research strength algorithms to model complex, inhomogeneous domain structures in magnetic materials. It uses a finite element/boundary element numerical method, optimally suited for calculating magnetization structures of local energy minima (LEM) in irregular grain geometries that are of interest to the rock and paleomagnetic community. MERRILL is able to simulate the magnetic characteristics of LEM states in both single grains, and small assemblies of interacting grains, including saddle‐point paths between nearby LEMs. Here the numerical model is briefly described, and an overview of the scripting language and available commands is provided. The open source nature of the code encourages future development of the model by the scientific community.


Introduction
Paleomagnetic observations have contributed a wealth of information about the evolution of the Earth and other planetary bodies (Dunlop & € Ozdemir, 2001;Merrill, 1998), through the interpretation of remanent magnetization in naturally occurring magnetic minerals. The vast majority of the natural magnetic archives are recorded by nanosized particles whose magnetic properties initially could only be inferred from experimental observations, either on natural samples or on man-made analogues of bulk particle arrays which have relatively broad particle size distributions. Rock magnetic interpretation is typically based on the assumption that paleomagnetic samples are dominated by single-domain particles which are theoretically accessible by N eel's theoretical description (N eel, 1955).
With the advent of numerical micromagnetics and high-resolution imaging techniques, it has been demonstrated that many of the remanence carriers in natural samples are inhomogeneously magnetized pseudosingle-domain (PSD) grains, which in their simplest form take on a single vortex (SV) magnetic structure. The first three-dimensional micromagnetic models were introduced in the late 1980s (Schabes & Bertram, 1988;Williams & Dunlop, 1989), and over the last 30 years have significantly advanced our understanding of magnetic recording in interacting and noninteracting nonuniformly magnetized particles. In paleomagnetism, for example, such studies have demonstrated not only that small PSD grains, ubiquitous in rocks, primarily occupy SV states but that these states can provide reliable and stable recording of the ancient magnetic field that remain stable of over billions of years (Nagy et al., 2017).
Modern desktop computers and workstations are now powerful enough to perform micromagnetic simulations of a wide range of grain sizes, materials, and even clusters of grains that are of interest to rockmagnetists, paleomagnetists, and environmental magnetists. The currently available open source programs (e.g. OOMMF (Donahue & Porter, 1999), Magpar (Scholz et al., 2003), and NMag (Fischbacher et al., 2007)) are developed for applications in material science and physics and require substantial specialist knowledge to install, maintain, and apply in Earth related contexts. The objective of this article is to introduce an open source micromagnetic model that we believe can significantly increase the ability of the rock and paleomagnetic community to define the magnetic properties of their samples through forward modeling of the behavior of their possible domain states.
We describe here a finite element micromagnetic modeling package called MERRILL (Micromagnetic Earth Related Robust Interpreted Language Laboratory) in honor of the early work of Ron T. Merrill on micromagnetics in rock magnetism (Merrill, 1977;Moon & Merrill, 1984, 1985. MERRILL is a script based modeling program designed for the Earth science community that needs no specialist computing knowledge or proprietary additional software to run the models and visualize the magnetic domain structures. Yet it is a fully tested research-strength modeling platform that uniquely provides specific features relevant for natural samples and may even outperform some of the above mentioned software in certain applications. Both, precompiled binaries and the FORTRAN source code are freely available and run on LINUX, macOS, and Windows. They can be downloaded from the MERRILL homepage at http://www.rockmag.org.

Workflow
A typical workflow using MERRILL is outlined in Figure 1. A tetrahedral mesh of the geometry of interest must first be generated in an external program (although MERRILL contains routines for generating some commonly used geometries), representing the magnetic material. An ''MScript'' command file is then passed to MERRILL to drive the model, e.g., loading the mesh, setting material parameters, varying external fields, minimizing the micromagnetic energy, and outputting the magnetization to disk. When the magnetization has been output to disk, an external visualization program can be used to inspect the results. The visualizations of the geometries and magnetizations presented in this paper were generated with ParaView (Ahrens et al., 2005). This magnetization may also be used as a starting point in future models.
A productive workflow might involve running a micromagnetic model, inspecting the results visually, and then continuing the model using the previous results as the new starting point, until some desired result is achieved.

Micromagnetism
In micromagnetism, a real physical magnetic system is described in terms of a continuous vector valued function:M : R 3 ! R 3 ; x 7 !MðxÞ; whereMðxÞ represents a mathematical magnetization vector at the mathematical pointx 2 R 3 (Brown, 1963;Hubert & Sch€ afer, 1998). All physical energies and processes are then mathematically studied through this continuous function. To link the results to the physical situation, one can assume that the magnetization at a given mathematical point is constructed by averaging the discrete physical sources of magnetism, e.g., electron spins and orbitals, over a small volume centered at that point. This volume must be large enough that the behavior of individual atoms are averaged out, but small enough to resolve inhomogeneous magnetic structures such as domain walls, vortices, or flower states. It turns out that this continuum approximation represents real magnetization structures astonishingly well, even if the mathematical magnetization at a given point represents only the average over a few atoms. Quantum mechanical effects, which are essential in magnetism, are represented purely phenomenologically through material constants like the exchange constant or the magnetocrystalline anisotropy constant. In micromagnetic models, the magnetic structure in a particle gives rise to various contributions to the total free magnetic energy functional EðMÞ, or to their associated effective fields. The effective fieldH is given as a function ofM as The dynamics of a micromagnetic system is described by the Landau-Lifshitz-Gilbert (LLG) equation (Gilbert, 2004): where c is the electron gyromagnetic ratio and k is a material-dependent damping factor. This behaves like a dynamical system where a forceH eff acts on a system with a nonzero angular momentum vector pointing parallel toM. A sufficient condition for an equilibrium magnetizationM 0 is given bỹ This solution is local energy minimum (LEM), which we might associate with a remanent magnetization state, and so we will denote itM rem . From equation (3) we can see two types of solutions. EitherM 0 is parallel toH eff ðM 0 Þ orH eff ðM 0 Þ50. In MERRILL, we focus on the second solution, that is, we solve for the LEM solutions by optimizing the expression for the total free magnetic energy of the system. This is generally much more efficient at finding stable domain structures that the full solution to the LLG equation, and we are not normally interested in the details of the domain transition dynamics that the LLG describes.

Effective Fields
The total effective fieldH eff ðMÞ for a typical cubic ferromagnetic crystal involves four primary components: Zeeman, anisotropy, exchange and demagnetizing fields (Brown, 1963;Kittel, 1949).

Zeeman Field
The Zeeman field represents the interaction of external sources of magnetic field with the magnetic system under investigation. As such, it is often referred to as the ''external field.'' It is assumed that the system under investigation has no effect upon the external source.

Anisotropy Field
The anisotropy field couples the magnetization to the crystal lattice. It is the primary mechanism by which the symmetries of the lattice affect the magnetization. For a cubic ferromagnetic crystal, if the crystal axes are along the x, y, and z coordinate axes, it can be writteñ with K 1 the anisotropy constant. In magnetite and iron higher-order terms are much smaller and can be safely ignored for micromagnetic calculations. However, MERRILL does allow a value for K 2 to be set.
The vectorã here represents the directional cosines of the magnetization with respect to the crystal axes. For a crystal with cubic axesã;b, andc, the vectorã is defined: such thatã Áã51.
For other noncubic anisotropies, appropriate equations are used, e.g., for a crystal with a uniaxial symmetry, e.g., a tetragonal mineral like tetrataenite, a uniaxial anisotropyH anis ðMÞ5K 1 ðM ÁãÞ 2 is commonly used, with a the anisotropy axis.

Exchange Field
The exchange field serves to align nearest neighbor magnetizations. Although this incorporates quantum mechanical spin coupling, in the continuum approximation a micromagnetic expression can be given as Geochemistry, Geophysics, Geosystems

Demagnetizing Field
The demagnetizing field represents the magnetic field generated by the magnetic material itself, derived from the ''magnetic self-energy.'' It is called the demagnetizing field, because a higher magnetic self-energy represents a higher energy configuration, so it typically acts on the magnetization in such a way that it minimizes external flux and thereby itself.H dmag 5r/; /ð1Þ 5 0: Outside the magnetic material, this is the ''stray field.''

Minimum Energy Solutions Using Finite Elements
There are a variety of numerical micromagnetic approaches than can be used to solve for locally stable magnetic domain structures. Here we use a method that is both numerically efficient and robust while requiring the simplest input for representing the geometry of the grain. In most micromagnetic models, the primary consideration is the efficiency with which in the internal demagnetizing field can be computed. This field calculation scales as OðN 2 Þ, where N is the number nodes, but if the geometry is meshed using a regular grid this scaling can be reduced to OðNlog NÞ using FFT (Fabian et al., 1996;Wright et al., 1997). However, such regular grids do not easily account for arbitrary grain geometries, although some success has been achieved by relative scaling of surface elements (Witt et al., 2005), or irregular FFT techniques (Kritsikis et al., 2008).
In MERRILL, we employ the Finite Element Method (FEM), a standard technique for describing functions over a geometry and solving differential equations in terms of those functions (Davies, 2011). MERRILL uses arbitrarily shaped linear tetrahedral finite elements to describe the geometry of a particle and to solve for the LEM stable domain states. Some care is needed in the calculation of the demagnetizing field described by equations (8-10), which involves solving the Poisson equation for the magnetic scalar potential over an infinite space.
MERRILL makes use of a Boundary Element Method (BEM) technique (Fredkin & Koehler, 1990;Lindholm, 1984), which is a specialization of the FEM for homogeneous Poisson equations, particularly suited to problems defined over an infinite space. This method has significant advantage in that we need not create a mesh in the free-space region outside the geometry of the magnetic particle (even for multiparticle solutions) but the method does, however, increase the memory requirements of the programme. In our experience, this is acceptable for single-grain geometries requiring up to about a million elements. A more detailed account of different mciromagnetic methodologies can be found in Fidler and Schrefl (2000).
LEM states are found by solving for the minimum free magnetic energy of the system E52H eff ÁM where the total effective field is given byH eff ðMÞ 5H exch ðMÞ1H anis ðMÞ 1H zm ðMÞ1H dmag ðMÞ: To include further phenomena, like stress fields or surface anisotropy, in theory one need only derive the corresponding effective field, typically by taking the derivative of the energy with respect to the magnetization, and add it to equation (11).
In order to determine minimum energy solutions, by default MERRILL makes use of an accelerated adaptive step-size steepest descent algorithm across the energy landscape, optimized for micromagnetics, here called ''Hubert Minimizer'' (Berkov, 1998a;Ramst€ ock, 1997). A standard conjugate gradient optimizer is also available as an option, but in most cases is slightly less efficient in finding energy minima.

Mesh Generation
MERRILL requires that the geometry of the magnetic particle is described using a linear tetrahedral mesh (by default defined in micron units). The magnetization is thus specified only at the four vertices of each element and linearly interpolated at all other locations. MERRILL is able to generate suitable meshes for simple grain geometries such as cubes and spheres, however, more complex particle geometries require additional meshing software of which there are many free and commercial programmes. The only requirement is that the mesh is formatted according to ASCII text based PATRAN (.neu) standard. Most finite element meshing applications will support this and more detailed information about the PATRAN format can be obtained from the website of the MSC Software Corporation (https://simcompanion. mscsoftware.com).
Like all finite element models, the quality of the mesh will affect the convergence efficiency of the model. In most cases meshing software will take care to produce a qood quality mesh, but for highly irregular geometries problems may still arise. Discussion of mesh quality metrics can be found in many publications (e.g., Dai et al., 2014;Knupp, 2006) For micromagnetic applications, it is important to ensure that that mesh is fine enough to resolve the expected spatial variation of the magnetization within the model geometry. The maximum element size is usually described in terms of the ''exchange length,'' l exch (Rave et al., 1998), which is dependent on the magnetic material parameters: For iron and magnetite, for example, the exchange length is around 3 and 9 nm, respectively, at 308C and slightly larger near the Curie temperature.
Clearly, the bigger particle and finer the mesh the longer the model will take to converge to a LEM state. For large grains, the mesh size is usually set to the exchange length. However, small grain geometries will require a mesh that is finer than the exchange length in order adequately represent the grain shape.

Model Validation
Validation of micromagnetic models should ultimately be done against experimental observations. However, such direct validation has until recently been extremely difficult to achieve since the maximum grain size that could be modeled numerically was much smaller than that which could be directly observed experimentally. As a result, the earliest micromagnetic models could only be validated against bulk observations such as assemblies of sized particle fractions or magnetosome observations (Fabian et al., 1996;Williams & Dunlop, 1995;Witt et al., 2005) in the case of natural materials, or on thin film and particulate man-made recording media (Labrune & Miltat, 1990;Silva & Bertram, 1990). However, results from MERRILL have been directly compared with nanoscale experimental data via electron holography with good agreement (Almeida et al., 2016).
With an increasing number of micromagnetic models being published, a number of standard tests were developed in order to ensure that the models were at least self-consistent. Validation against these standard tests is now a prerequisite for any newly published micromagnetic code. One such test is lMAG Standard Problem 3 (see Appendix A), which tests for the critical edge length of a cube for transition between a flower state and a vortex state. Our solutions found the flower and vortex states had equal energies at an edge length of 8.47 l exch when extrapolated to an infinitely fine mesh, which is in good agreement with other submissions to the lMAG problem.
MERRILL also tests the effective field components against some analytic solutions. For example, the demagnetizing field for a uniform sphere can be writtenH dmag 5 1 3M and should be independent of the sphere size.

MScript: The MERRILL Scripting Language
MERRILL is run at the command line with an input text file containing a series of MScript commands for MERRILL to execute. These commands allow the user to interact with MERRILL in a number of ways, for Geochemistry, Geophysics, Geosystems 10.1002/2017GC007279 example, setting material constants, loading meshes or finding a local energy minimum. In this section, we will outline a few typical computer experiments that can be used to probe the behavior of the magnetization of magnetic materials using MERRILL. This will also serve as a brief tutorial on the MScript language and introduce some of features in MERRILL.
The default units for MERRILL are microns and degrees Celsius, although different units may be used as specified in the user documentation. A full list of MERRILL commands can be found in Appendix B and on MERRILL download page (http://www.rockmag.org) where it will be updated as MERRILL is developed further.
Some typical commands include (but are not limited to) Magnetite htemperaturei C Set material constants to magnetite at htemperaturei degrees Celsius.
Iron htemperaturei C Set material constants to iron at htemperaturei degrees Celsius.
GenerateCubeMesh hwidthi hedgeleni Generate a cubic geometry of width hwidthi and with the average length between mesh nodes of length hedgeleni.
Uniform Magnetization hx; y; zi Set the magnetization of the material tõ Run the minimizer to find the local energy minimum.

WriteMagnetization hfilenamei
Write the magnetization to disk. Two files are written, hfilenamei.dat and hfilenamei_mult.tec. The file hfilenamei.dat contains a list of vertex points and the magnetization at that point. This is suitable for use with the ReadMagnetization command (see below). The hfilenamei_mult.tec file is the mesh and the magnetization in a TecPlot file format. This can be read by a number of visualization tools, including the free, open source and cross-platform ParaView software (Ahrens et al., 2005).

ReadMagnetization hfilenamei
Read a previously written file hfilenamei from the WriteMagnetization command, and set the current magnetization based on that.
External Field Direction hx; y; zi Set the Zeeman field direction toĤ zm 5ðx; y; zÞ= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 1y 2 1z 2 p External Field Strength hmagnitudei mT Set the Zeeman field strength to jH zm j5hmagnitudei MScript also includes supports for variables and loops. For example, the script Loop myvalue 0 100 10 Print #myvalue EndLoop creates a loop with the local variable myvalue. The numbers in the loop command determine initial value, end value and step size, such that the above loop prints out the values 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, and 100. To distinguish between floating point, integer and string values, the variable myvalue can be referred to as %myvalue, #myvalue, and $myvalue$, depending on what is needed. For example, in Geochemistry, Geophysics, Geosystems

Minimization
A basic minimization from a random start, or a uniformly magnetized start can be a quick way to tell if a grain is in the SD or PSD size range. Moreover, it is the basic building block from which most other MERRILL modeling experiments come.
A simple minimization script to determine a LEM state from a random initial state for a meshed geometry can be given as follows: ! Setup material constants for Iron ! at 20 Degrees Celsius Iron 20 C ! Load the meshed geometry from a ! Patran Neutral file ReadMesh 1 mesh.neu ! Randomize the magnetization Randomize All Moments ! Run the minimizer Minimize ! Output the solution M to two files: ! soln.dat and soln_mult.tec. ! The soln.dat file can be used as the ! input for another run as, e.g., an ! initial guess. soln_mult.tec is a ! TecPlot format file which can be used ! to view the solution in a visualization ! program such as ParaView ! (which is free and open source!) WriteMagnetization soln The output of this script can be seen in Figure 2 using a mesh generated from experimental focused ion beam nanotomography data for an example iron inclusion in dusty olivine in a chondritic meteor sample (Einsle et al., 2016).
It is important to note that there may be several possible magnetization states the minimization might settle upon, representing multiple possible, and completely valid local energy minima states. However, from a given starting point, barring numerical noise, the minimization should always reach the same end point.
The LEM found can be highly dependent on the choice of initial state. A simple minimization script to determine a LEM state for a cubic grain of magnetite, starting with uniform magnetization, can be given as follows: ! Setup material constants for Magnetite ! at 20 Degrees Celsius Magnetite 20 C ! Generate a.08 micron cube using the Geochemistry, Geophysics, Geosystems

10.1002/2017GC007279
! built-in cube mesh with mesh size 5 $nm ! generator GenerateCubeMesh 0.080 0.005 ! A flower state will nucleate from a ! uniform magnetization along the ! easy axis Uniform Magnetization 1 1 1

Minimize WriteMagnetization flower_soln
The output of this script can be seen in Figure 3, representing a flower state.
By starting from a rough vortex state rather than a uniform state, we can instead nucleate a vortex LEM. This can be accomplished by Magnetite 20 C GenerateCubeMesh 0.080 0.005 ! A vortex state will nucleate from a ! rough vortex state with core aligned ! along the easy axis ! The 1 1 1 refers to the direction, ! 0.02 is the ''tightness,'' which should ! by manually tuned by the user, and ! LH will produce a left hand vortex. Vortex Magnetization 1 1 1 0.02 LH

Minimize WriteMagnetization vortex_soln
The output of this new script (seen in Figure 4) highlights the ''local'' aspect of a local energy minimum.

Hysteresis Loops
Hysteresis loops are a useful tool when a magnetic material has several remanent magnetization states for a given set of parameters. From a hysteresis loop, it is possible to deterministically move from one state to another, and find the tipping point where the variation of a given parameter will cause one state to spontaneously switch to the other. This also provides information on the range of values for the given parameter where both states can coexist. In other words, it provides some information about the stability of a given state when multiple valid states exist.
Since the minimization should always reach the same solution when given the same starting point, a hysteresis loop run with the same parameters and the same changes in parameters should always return the same results in MERRILL. 6.2.1. Magnetic Field Hysteresis A magnetic hysteresis loop has several uses. A single hysteresis loop can be useful to get a feel for the behavior of the magnetization of a system (i.e., its coercivity, if it is SD or PSD). The average of hysteresis loops in many directions can be used to compare simulations with experimental observations of magnetic characteristics such as coercivity.
A hysteresis loop is a quasi-static thermodynamic process. That is to say, at each point of the hysteresis loop, the system is assumed to be in equilibrium and the energy at a local energy minimum. A change from one point to the next in a hysteresis loop represents the system moving from what was a local energy minimum in the previous step to the nearest, new local energy minimum in the current step. The simulation of a hysteresis loop is accomplished by the following scheme: 1. Set Zeeman field to saturating value in a fixed direction: Z 5 Zmax. 2. Find LEM. 3. Update Zeeman field: Z 5 Z 2 a*Zmax. 4. If Z !5 2Zmax, go to 2.
For a saturating fieldZ max , a hysteresis loop will run fromZ max to 2Z max in small increments. Here ''small'' means sufficiently small enough that any phenomenon sensitive to the external field, e.g., magnetization switching or nucleation of vortices, are accounted for. In practice, ''small'' means just small enough to resolve changes in the magnetization you are looking for, but as large as possible to reduce the total number of steps needed for the loop. Otherwise a large amount of work may be done unnecessarily for little progress.

A MERRILL script which accomplishes this is
Magnetite 20 C GenerateCubeMesh 0.100 0.005 External Field Direction 1 0 0 External Field Strength 0 mT Loop field 2100 100 5 Randomize Magnetization 10 External Field Strength Minimize ! Write current value to disk, ! so we can inspect it later.

WriteHyst hyst_soln EndLoop
Since hysteresis is generally symmetrical to forward and reverse fields, it is only necessary to run either the upper or lower branch of the hysteresis loop, not both.
A plot of the upper branch is shown in Figure 5. This curve is different to the usual SD curve, where the valueM ÁH remains near the saturation value, until passing the coercive field. Here the grain enters a vortex state, and the valueM ÁH can minimized while keeping the anisotropy energy relatively low by varying the shape of the vortex. An example of the vortex state at jHj5230 mT is shown in Figure 6. 6.2.2. Size Hysteresis A hysteresis loop where the size of the grain is varied rather than an external field can be a useful tool for determining SD and PSD ranges, and the evolution of magnetic domain states. Since several remanent states can exist for a grain, particularly in the early PSD size range, it can be difficult to pinpoint exactly where that regime begins.
In a size hysteresis loop, the grain is started in an SD magnetization state. The size of the grain is increased until it spontaneously switches to a vortex PSD state. The size of the grain is then scaled down until it is in a SD state again. The branch of increasing size can tell us what the largest grain size is that supports an SD state, and the branch of decreasing size can tell us what the smallest grain size is that supports a PSD state.  Geochemistry, Geophysics, Geosystems

10.1002/2017GC007279
Note that the exact initial SD was not defined. Ideally, a size hysteresis should be run for each possible SD remanent state. For symmetric grains, symmetry should make many of these redundant. However, for asymmetric grains, the direction of the initial magnetization may greatly affect the stability of the SD solution.
When scaling sizes, it is important to consider the size of the elements of the mesh during scaling. Typically, the edges of a mesh should be around the exchange length. When the size of a mesh is increased, the average edge length should not exceed the exchange length. By starting with a larger mesh and scaling down, rather than the other way around, we can avoid this problem. However, finer meshes take longer to run. A compromise can be found using the ''Remesh'' command. A user might use a mesh that is suitable for scaling up to, e.g., 0.1 lm, and at that size, switch to a mesh suitable for up to, e.g., 0.2 lm.
So if we want a maximum node spacing of, e.g., 0.005 lm, we make two meshes of width 0.1 lm: one with node spacing 0.005 lm, and one with node spacing of 0.0025 lm. The first mesh will cover grains from 0 to 0.1 lm, while the second can cover grains from 0.1 to 0.2 lm. We make the meshes the same initial size so the Remesh command can interpolate the magnetization directly from one to the other.
In MERRILL, the ''Resize'' command can be used to scale the mesh. An example size hysteresis script incorporating all of this is ! Use magnetite material parameters Magnetite 20 C ! Ensure we can load at least 2 meshes ! at a time Set MaxMeshNumber 2      Figure 8, demonstrating the ESV and HSV states for the lower branch, and the corresponding FS states on the upper branch. 6.2.3. Temperature Hysteresis A temperature hysteresis can be used to observe the effects of heating, say iron, from room temperature to the Curie temperature and back. The approach MERRILL takes to this is to simply change the temperature dependant material parameters for each step. In this way, the model is of a ''cold'' material, since no other thermal effects are taken into account which might spontaneously and nondeterministically effect the magnetization. However, the changes in material properties can have significant impact on the shape of the energy landscape with respect to the depths, positions and even number of the local energy minima.
An example temperature hysteresis of iron is GenerateCubeMesh 0.18 0.005 Uniform Magnetization 1 1 1 Loop temperature 20 770 10 ! Set the material parameters for iron ! at the given temperature Iron %temperature C

! Find the LEM Minimize
! Write current value to disk, so we ! can inspect it later. WriteHyst hyst_soln_#temperature EndLoop 6.3. Energy Barriers MERRILL includes a new method for finding the minimum energy transition between two given LEM states (K. Fabian & V. P. Shcherbakov, Energy barriers in three-dimensional micromagnetic models and the physics of thermo-viscous magnetization in multidomain particles, arXiv 1702.00070v1, 2017). It uses a combination of the Nudged Elastic Band technique (Dittrich et al., 2002;Henkelman & J onsson, 2000) and an action minimization method (Berkov, 1998a(Berkov, , 1998bFabian & Shcherbakov, arXiv 1702.00070v1, 2017. It finally constructs a prescribed number of intermediary states between the given LEM states L 1 , L 2 , such that their interpolation represents the average physical switching path for a fully damped (k 5 0) switching process from L 1 to L 2 . The path with the lowest maximum energy barrier determines the thermal activation energy required to switch between these states. From that, the probability of switching and the relaxation time across the energy barrier can be obtained.
These quantities are of paramount importance in paleomagnetism, where small changes in size, shape, or material parameters can change the relaxation time from the order of seconds, to millions of years, and the corresponding paleomagnetic information is either completely unblocked, or can be regarded as stable or blocked. In contrast to, for example, the simulation of writing heads for magnetic hard drives, the nanosecond dynamics of the switching processes are of no interest in paleomagnetism, such that time-consuming full-fledged LLG models do not improve the result. On the contrary, the exact determination of the energy barrier through the saddle-point path provides a better estimate of the thermal relaxation time. Note in that respect that the LLG equation does not include temperature, and that thermal activation has to be Geochemistry, Geophysics, Geosystems 10.1002/2017GC007279 added in a phenomenological way, e.g., as a stochastic zero-mean Gaussian fluctuation field (Torres et al., 2001), if it is considered at all.
If morphology and composition of a natural magnetic mineral system is known MERRILL can provide quantitative insight to whether the signal recorded has remained stable since it was recorded (Nagy et al., 2017).

Practical Considerations
There are a number of practical considerations when using MERRILL. Presently, MERRILL is not parallelized. However, many instances of MERRILL may be run in task-farm parallel manner for reasonably small grains. This approach can be used to model, for example, large assemblies of noninteracting grains, and hysteresis loops about many axes in parallel. The maximum grain size that you will be able to model will be dependent on a number of factors, but primarily the number of elements in the model. MERRILL will comfortably cope with models of up to approximately 10 6 elements for single grain of simple geometry on a machine with access to 64 GB of RAM (see Figure 9). If the model consists of a cluster particles where the total grain surface to volume ratio is significantly greater than that of a single grain then the maximum number of elements in the model would need to be less than 10 6 . It is not possible to formulate any hard rules that cover the combination of different number of particles, grain sizes and mesh sizes, but the user should be aware that the memory requirement increases as the square of the number of surface nodes (which is typically proportional to the number of surface elements).
The time taken for a model to converge on a local energy minimum will again be largely dependent on the model size (see Figure 9), but also on how close your initial guess is to the final LEM solution. Thus solution times for a random initial guess and a uniform saturated state initial guess can have significantly different solution times. A similar situation occurs during a change in the magnetic phase, i.e., when a solution changes from SD to SV, or hard-aligned SV to easy-aligned SV.
In rare cases, the model might attempt a large number of iterations because of slow convergence to the final solution. In such cases it is better to place a limit in on the maximum number of allowed iterations (which has default value of 100,000). Often such slowly converging solutions can be avoided by trying a slightly different initial guess, slightly changing the model mesh, or adding a small random kick to the solution. Figure 9. A log-log plot of the time taken for MERRILL to find a local energy minimum of a cube, starting from a random state, versus the number of elements in the mesh.

10.1002/2017GC007279
As mentioned in section 4, the mesh shape and quality also impact on convergence times. Regular tetrahedra are the most numerically stable elements, so the more regular the tetrahedra in a mesh, the quicker the model will converge. The single most important concern, however, is to avoid slivers. Slivers are elements that are so flat, their volumes approach zero (Alliez et al., 2017). These lead to extreme numerical instability and very slow convergence times, if the model converges at all. Slivers might occur during Delaunay triangulation, as the algorithm is well known to be blind to a class of slivers where the vertices all lie on the same circumcircle. An additional mesh optimization step like Lloyd optimization, or ODT smoothing is often performed during mesh generation to eliminate slivers. Most finite element meshing packages will have algorithms that attempt to detect and avoid such elements, but in all cases when the model fails to converge it is important to check the mesh quality (e.g., using the ''Mesh Quality'' filter in ParaView).
Also mentioned in section 4 is that mesh edge lengths should be smaller than the exchange length. In most meshing programs, it is unclear what the defined ''mesh size'' is. In some cases, it is the upper bound of the radius of the circumcircle of an element allowed before the element is decomposed into smaller elements. In some cases, it is the diameter of this circumcircle. For some software, it is the target average edge length of the elements in the mesh. In any case, the ReportEnergy function of MERRILL reports the average edge length of the mesh in nanometres, which is then directly comparable to the exchange length reported by this function in nanometres. Users are encouraged to use these values as a useful and direct report of whether the mesh is fine enough for the given material.
The solvers used in MERRILL can be susceptible to finding unstable equilibria. Gimbal locking, for instance, can occur for field hysteresis loops directly along the hard-axis of the material. Two simple techniques can be used to avoid this. First, prefer not to perform hysteresis loops directly along the major axes. Add a slight perturbation, i.e., use the [1 1 1.0001] direction instead of [1 1 1]. Another technique is to add a small random ''kick'' to the magnetization. Typically, moving each magnetization randomly by around 5 or 10 degrees is not so large that it will move the state out of a deep, stable LEM, but it can move it off a local maximum, and also out of a shallow LEM, where the magnetization is around a SD/SV phase transition.

Discussion
The formulation of a micromagnetic model in an easy to use form can significantly enhance its application in paleomagnetic and rock magnetic investigations. MERRILL presents such a tool that is particularly focused on finding remanent states and studying their stability. The parallels and usefulness to paleomagnetic and rock magnetic studies should be clear. MERRILL has been used in a number of publications and talks using functionality not presented in this paper, e.g., behaviors of assemblages of interacting grains, large models, strongly anisotropic materials and detailed simulations of magnetostrictive effects (e.g., Almeida et al., 2016;Chang et al., 2012;Conbhu ı et al., 2016;Einsle et al., 2016;Li et al., 2013;Nagy et al., 2017;Williams et al., 2010). Not all of these, however, used the scripting interface presented here.
The simple scripting language makes it particularly friendly to nontechnical users. The fast and efficient minimization scheme means simple computer experiments can be run quite quickly. The scripts presented here represent the sort of scripts run by the authors day to day. Some effort has been made to make these copy/pasteable, but a curious reader is recommended to look in the ''demo'' directory of the MERRILL package for more information and examples.
If the scripting language is not up to a particular task, MERRILL can also be used as a library and called from a Fortran program. In addition, a plugin interface has been included so that users can compile libraries that can be loaded from the MScript interface and hook into the MScript parser to add commands and variables, and also add effective field calculators toH eff not originally shipped with MERRILL. This should allow MER-RILL to be adopted to a wide range of problems. Geochemistry, Geophysics, Geosystems 10.1002/2017GC007279

A1. Introduction
The lMAG Standard Problem 3 is a test for the critical edge length of a cube with uniaxial anisotropy for a change in the magnetic phase from a flower state to a single vortex state (http://www.ctcms.nist. gov/$rdm/mumag.org.html). It was introduced in 1998.
In Standard Problem 3, the material parameters of the cube: the saturation magnetization M s , and the uniaxial anisotropy constant K u , are related by the relation: ) and the exchange length l exch is given in terms of the uniaxial anisotropy constant, the exchange coupling A l exch 5 Existing results suggest the critical edge length, where a flower state and a vortex state should have the same energy is in the region of 8.453l exch to 8.55 3l exch .

A2. Method
The material parameters used were The exchange length was therefore l exch 59:587248310 23 lm: To nucleate a flower state, the mesh is scaled to 7.53l exch and a local energy minimum (LEM) is found, then rescaled to 8.453l exch , our expected lower bound, and the LEM found and saved to disk. By ''scaled to,'' we mean the mesh is resized until edge length of the cube is the given value. This two step approach to nucleation is needed, since a flower state is not guaranteed to nucleate at 8:453l exch , as this is around the critical edge length; reminimizing at 8.453l exch is done primarily to save time during minimization later.
For any particular domain state, its form will change slightly with grain size. A flower at a larger grain size, for example, will have a more divergent magnetization at the grain surface than a flower state at a smaller grain size. In that manner, states from grains closer in size tend to be more ''similar'' to each other than states from grains less close in size. Since the size 8.453l exch is closer in size to our expected critical region than the initial size of 7.53l exch , the flower state at 8.453l exch should be more ''similar'' to the states about the region of interest too. As a result, minimizing from the more ''similar'' state at 8.453l exch to states about the region of interest takes less time than minimizing from the state at 7.53l exch . An initial flower state at 8.453l exch is shown in Figure A1.  Geochemistry, Geophysics, Geosystems

10.1002/2017GC007279
To nucleate a vortex state, the mesh is first scaled to 113l exch to guarantee vortex nucleation and the LEM found, then scaled to 8.453l exch and the LEM found and saved to disk. Again, reminimization saves a lot of time as the vortex at 8.53l exch is more ''similar'' to the states around the region of interest than the state at 113l exch is. An initial vortex state at 8.453l exch is shown in Figure A2.
To evaluate, say, the flower state energy at a certain cube size, say L, the mesh is scaled to length L, the flower state at is loaded from disk, which for the flower was found for 8.453l exch , and the energy is minimized. The energy is then found and written to disk. This is done for the flower and vortex states for a range of cube sizes about the expected critical point of 8.53l exch .
The critical length for the cube is then found by finding the iteration just before the flower state energy passes the single vortex energy and the iteration just after, and using a linear interpolation between the energies of each iteration to find the edge length where they intercept. Denoting the first iteration's flower energy, vortex energy, and length scale as f 2 ; v 2 ; l 2 , and the second iteration's as f 1 ; v 1 ; l 1 , the intercept length scale, l 0 , is l 0 5l 2 1 ðl 1 2l 2 Þðf 2 2v 2 Þ ðf 2 2v 2 Þ2ðf 1 2v 1 Þ : This assumes the energies scale linearly with the edge length, which is incorrect, but for sufficiently small steps of the edge length, it should be accurate enough for this application.
This entire process is then repeated for an increasingly fine mesh. The critical edge length of the cube should scale as the square of the mesh spacing. This is due to higher resolution of the derivatives, which, for the demagnetization and exchange calculations are of the order two. By performing a linear fit on the (mesh spacing) 2 versus the intercept length, it is be possible to extrapolate the intercept length for an infinitely fine mesh, i.e., where the mesh spacing is zero. Figure A3. Cube mesh with a node spacing of 0.53l exch . Critical edge length versus node spacing. Figure A4. The critical length versus the node spacing of the mesh, and linear fit extrapolating the node spacing to 0. This shows the critical edge length for an infinitely fine mesh using MERRILL is ð8:46860:002Þ3l exch . Critical energy versus node spacing.

10.1002/2017GC007279
For each mesh, the energies were found for the flower and vortex states for 8.393l exch to 4.703l exch in steps of 0.013l exch . The mesh spacings used were 1.03l exch to 0.33l exch in steps of 0.13l exch . The meshes were generated using a Delaunay triangulation algorithm from the CGAL library (The CGAL Project, 2017). An example mesh at node spacing 0.53l exch is given in Figure A3.

A3. Results
From Figure A4, the critical edge length for an infinitely fine mesh was found to be ð8:46860:002Þ3l exch , and from Figure A5, the critical energy was ð0:3024260:00008Þ3K d V. A break down of the demag, exchange, and anisotropy energies at the critical length for the flower and vortex states is in Table A1. The average magnetization for these states at the critical length is in Table A2.
These measurements were also carried out with MERRILL's built in regular mesh generation for cubes, and the difference in results was within 1%.
A4. Discussion MERRILL's solutions to the lMAG Standard Problem 3 are in very good agreement with the other submissions. This is good evidence that the energy terms in MERRILL are correctly determined and converge to the true continuum value when mesh size approaches zero.
The error values presented here are for the extrapolation of the critical values MERRILL would find on an infinitely fine mesh. They are not errors in the actual values generated by MERRILL compared to the physical equivalent. There are many sources of error unaccounted for by this approach, including material parameters, mesh geometry, numerical errors, and approximation errors. For this reason, when comparing between MERRILL's results, and other submissions to the lMAG site, it may be more informative to directly compare the extrapolated values than to compare the extrapolated values within the presented errors.

Appendix B: MERRILL Commands
The basic functionality of MERRILL provides a boundary element micromagnetic energy calculation for arbitrary particle shapes. Figure A5. The critical energy versus the node spacing of the mesh, and linear extrapolating the node spacing to 0. This shows the critical energy for an infinitely fine mesh using MERRILL is ð0:3024260:00008Þ3K d V. To use MERRILL requires at least 1. A particle geometry (e.g., a cube or octahedron) that describes the particle to be modeled. This geometry is translated into a tetrahedral mesh with a large set of model nodes by a mesh generator. The distance between the nodes defines the resolution of the model. Several ready to use meshes are provided together with MERRILL, but arbitrary meshes can be loaded as long as they fulfill some minimal quality criteria. Besides commercial software, there exist several professional open source mesh generators that are more than sufficient for generating simple particle meshes. For example, http://engits.eu/en/engrid, http://sourceforge.net/projects/netgen-mesher/, and http://geuz.org/gmsh/, http://meshlab.sourceforge.net/. 2. Material properties of the magnetic material to be modeled and the scaling constant determining the particle size. MERRILL provides a database of material properties for several common natural magnetic minerals. 3. With the above prerequisites it is already possible to do very sophisticated model runs resulting, e.g., in micromagnetic energies and mean magnetization values. To visualize the resulting output files MERRILL is set up to easily interact with ParaView.

B1. Interpreter Language
MERRILL is operated through an interpreter language script that determines all aspects of the micromagnetic calculation.
The script is a simple ASCII file containing a sequence of lines. Empty lines, leading, trailing, and multiple spaces or tabs are ignored, as well as anything behind an exclamation mark (!). A number of keywords are used to call subroutines or perform simple assignments. All keywords are case insensitive, e.g., ''ReadMesh'' and ''readmesh'' are equivalent. The script file is parsed line by line. Each valid line is immediately interpreted and executed.

B1.1. Basic Commands
The following list contains the currently available commands.
Set AEvariableae AEvalueae is used to define global variables for the material, the geometry of the mesh, or program parameters. The following variables are supported: Ms saturation magnetization in A/m.
K1 anisotropy constant for uniaxial or cubic anisotropy in J/m 3 .
Aex exchange constant in J/m.
Ls inverse length scale 1/m. Internally Ls 2 is used.
mu related to permeability of free space via mu 5l 0 =4p.
CurvatureWeight Weight of curvature contribution for nudged-elastic band method (NEB).
MaxMeshNumber Maximal number of finite element meshes stored. Must be set once before loading meshes.
PathN Number of structures along the magnetization path. Warning: The mesh must have been defined previously! Use only after ReadMesh.
MaxRestarts Maximum number of restarts during energy minimization.
Zone Current Zone to be written into the TecPlot output file (double). Zone can be set before each output, or c an be used with automatic increment.
AllExchange Values greater zero imply that all exchange energy calculations are performed and logged. This slows the program down.
Magnetite AEtemperatureae C Defines material constants for magnetite at temperature htemperaturei to be used in subsequent calculations. The temperature is given in degrees Celsius and must be positive and below magnetite's Curie temperature at 5808C.
Resize AEold lengthae AEnew lengthae changes the length scale of the mesh such that after this rescaling the length hold lengthi will become hnew lengthi.
(Cubic j Uniaxial) Anisotropy Defines the symmetry of the anisotropy energy.
CubicRotation AE/;t;aae Rotates the cubic easy axes using the rotation matrix with the prescribed Euler angles.
CubicAxes AEa x ;a y ;a z ae AEb x ;b y ;b z ae AEc x ;c y ;c z ae Directly set the 3 cubic axes,ã;b, andc, for the cubic anisotropy. Users should ensure the axes are mutually orthogonal.
External field direction AEx;y;zae Determines the direction vector of an external homogenous magnetic field. Intrinsically sets the field to B 5 1 T.
External field strength AEBae AEunitae Determines the strength of the external homogenous magnetic field as hBi in units of huniti. Possible values for huniti are ''lT,'' ''mT,'' or ''T.'' Must be set after defining the direction. Subsequent calls reset the field to hBi without changing the direction. Can be used for hysteresis modeling.
ReadMesh AEindexae AEfilenameae Reads the Patran file hfilenamei, and stores the corresponding mesh and finite element arrays at location hindexi. The index must be less or equal to the previously set MaxMeshNumber.
LoadMesh AEindexae Loads a previously read mesh and its finite element arrays from location hindexi. This mesh will then be used in subsequent operations.
ReadMagnetization AEfilenameae Reads a magnetization file, hfilenamei, (.dat) into the current mesh magnetization array. Make sure that it was created for the currently active mesh! The magnetization read is used in subsequent operations.
Uniform magnetization AEx;y;zae; [AEbae] Creates a uniform magnetization for the current mesh pointing in the normalized direction hx; y; zi. The optional parameter hbi is the index of the block of nodes in the mesh that should be set. By definition block 1 contains the free nodes, while higher block numbers can be used to define fixed nodes. These blocks have to be defined in the Patran file. Any previous magnetization is lost. Invert magnetization AEs x ;s y ;s z ae Inverts the current magnetization structure by multiplying each component of ðm x ; m y ; m z Þ with the corresponding sign ðs x ; s y ; s z Þ. For example, Invert magnetization 1 21 1 on each node changes the current magnetization ðm x ; m y ; m z Þ7 !ðm x ; 2m y ; m z Þ.

Vortex magnetization
Randomize magnetization AEangleae Randomly changes each current magnetization vector by at most hanglei degrees. The previous magnetization is (partially) lost.
Randomize all moments Replaces the current magnetization by randomly distributed unit vectors. Any previous magnetization is lost!
ReMesh AEindexae Takes the current magnetization array and interpolates it at the nodes of the previously read mesh at location hindexi. This mesh from location hindexi is then loaded and will be used in subsequent operations.
ConjugateGradient Uses conjugate gradient steps during the accelerated descent.
SteepestDescent Uses normal gradient steps during the accelerated descent.
Geochemistry, Geophysics, Geosystems WriteMagnetization AEfilenameae Saves the magnetization and the mesh in two files: hfilenamei.dat contains vertex coordinates and magnetization vectors.
hfilenamei mult.dat TecPlot file for visualization using ParaView or TecPlot. Contains mesh geometry and one or more magnetization states.
WriteHyst AEfilenameae Saves hysteresis data in 5 columns ofH zm ÁM; jH zm j and the 3 components on the average unit magnetization vector, hm x i; hm y i; hm z i, where M is the magnetization and jH zm j is the magnitude of external field.M ÁH zm is normalized to the saturated magnetization in the direction of the applied field.
SystemCommand AEcommandae . . . Performs the system command in the remaining arguments as a line. No guarantee can be given for correct behavior. Uses FORTRAN's SYSTEM command.
KeyPause Pauses script evaluation and waits for Key 1 Enter for continuation.
GradientTest Test feature to check energy gradient calculation against several finite-difference estimates.
(Stop j End) Stops script evaluation.

B1.2. Path Related Commands
The following contains all path related commands available for NEB calculations.
MagnetizationToPath AEindexae Saves the current magnetization in the path at location hindexi. This allows to assemble a path from individual magnetization states that have to fit to the current mesh! After assembling a path it must be renewed before further operations can be performed.
PathToMagnetization AEindexae Moves the path magnetization state at location hindexi to the current magnetization. This allows to change individual magnetizations in the path. For example, Initial and final states of a path read from a file can be minimized for new material constants.
RenewPath Defines all path variables, like distances and tangent vectors, assuming that all magnetizations have been correctly filled.
RefinePathTo AEnewlengthae Refines the current path to a new number of states by linear interpolation in the magnetization angles. This also resets PathN to the new value and renews the path. Of course, the new number of states can also be less than the previous PathN.
WriteTecPlotPath AEfilenameae Exports the current path to a TecPlot file with name hfilenamei. All states along the path are individual zones in the TecPlotFile.
ReadTecPlotPath AEfilenameae Reads a new path from a TecPlot file with name hfilenamei. All states along the path are individual zones in the TecPlotFile. Because this also reads in the mesh, all mesh related quantities are recalculated. Sets PathN and allocates necessary space for mesh and path arrays.
Make sure that all material parameters are correctly assigned, since those are not read! MakeInitialPath Assumes that a path is defined by set PathN hnumberi and that the first and last magnetization patterns are defined. Then proceeds by stepwise minimization to construct an initial path for subsequent optimization by the NEB method.
PathMinimize Assumes that an initial path is defined and minimizes the action integral using a variant of the NEB method.
PathLogfile AEfilenameae Starts logging all subsequent path minimization calculations into three logfiles hfilenamei.enlog hfilenamei.grlog, and hfilenamei.dlog. They contain energies along the path, norms of the gradients along the path and cumulative distances along the path. Logging can be stopped by EndLog or CloseLogfile.

10.1002/2017GC007279
PathStructureEnergies [AEfilnameae] Calculates the energies for each structure along the currently loaded path (in Joule). If hfilenamei is omitted, the path is written to standard output.
Energy Calculates the energy (in Joule) of the currently loaded magnetization structure.

B1.3. Interpreter Language Control Constructs
The following commands are simple options of the interpreter language to implement loops or script variables. Note that they cannot be nested or used recursively.
Loop AEvariableae AEstartvalueae AEendvalueae [AEstepae] Takes all commands until the next EndLoop statement and performs a loop over the enclosed commands by replacing the variable hvariablei with values from hstartvaluei to hendvaluei in steps of hstepi. If step is not given step-size step 5 1.0 is assumed. Within the loop the string #hvariablei is replaced by the integer value of variable, the string %hvariablei is replaced by the double precision value of variable, and the string $hvariablei$ is replaced by a string of the value of variable. Nested loops are not yet supported. Warning: The Loop command itself must NOT contain other variables than the loop variable. This is so because currently the parsing for replacing variables is performed only after unravelling the loops.
EndLoop Delimits the set of commands inside the active loop. Figure C4. Exchange coupling versus temperature for iron. Figure C5. Crystalline anisotropy constant versus temperature for iron.

10.1002/2017GC007279
Define AEvariableae AEvalueae Defines a numeric variable that can be used like a loop variable.
AddTo AEvariableae AEvalueae Adds a number to a previously defined variable.
Undefine AEvariableae Forgets a previously defined variable.

Appendix C: Material Parameters
The micromagnetic models require values of the three temperature-dependent magnetic parameters of saturation magnetization (M s ), the crystalline anisotropy constant (K 1 ), and the exchange constant (A). Each parameter is represented by a polynomial, based on the best fit to experimental data.

C1. Magnetite
For magnetite the experimental data was derived from Heider and Williams (1988) for the exchange constant, Fletcher and O'Reilly (1974) for the crystalline anisotropy and Heider et al. (1987) for the saturation magnetization. The polynomial expressions used in MERRILL for each of these three magnetite parameters is where T is the temperature in Celsius and T C is the Curie temperature for magnetite taken as 580 C. These parameter polynomials are plotted in Figures 15-17.

C2. Iron
For iron the experimental data for M s is derived from Crangle and Goodman (1971) and Cullity and Graham (2008), where the sparsity of data near the Curie temperature in tabulated values of Crangle and Goodman (1971) are supplemented data digitized from the M s -T graphs of Cullity and Graham (2008). For K 1 , the experimental data are taken from Honda et al. (1928) reproduced in graphical form in Cullity and Graham (2008). Although this K 1 material parameters fall to zero as the temperature approached the Curie point, more recent studies (e.g., Muxworthy & Williams, 2015) use a slightly higher room temperature value. In order to reconcile this discrepancy, the data from Honda et al. (1928) was scaled by a factor (480/456). Finally the exchange constant, A, for iron was obtained following the method outlined in Heider and Williams (1988), using the stiffness constant values reported by Stringfellow (1968).