Stray Magnetic Field Variations and Micromagnetic Simulations: Models for Ni_0.8Fe_0.2 Disks Used for Microparticle Trapping

Abstract

Patterned micro-scale thin-film magnetic structures, in conjunction with weak (~few tens of Oe) applied magnetic fields, can create energy landscapes capable of trapping and transporting fluid-borne magnetic microparticles. These energy landscapes arise from magnetic field magnitude variations that arise in the vicinity of the magnetic structures. In this study, we examine means of calculating magnetic fields in the local vicinity of permalloy (Ni_0.8Fe_0.2) microdisks in weak (~tens of Oe) external magnetic fields. To do this, we employ micromagnetic simulations and the resulting calculations of fields. Because field calculation from micromagnetic simulations is computationally time-intensive, we discuss a method for fitting simulated results to improve calculation speed. Resulting stray fields vary dramatically based on variations in micromagnetic simulations—vortex vs. non-vortex micromagnetic results—which can each appear despite identical simulation final conditions, resulting in field strengths that differ by about a factor of two.

1. Introduction

Micromagnetic particles yield a tool in biology for the isolation and quantification of biomolecules such as DNA [1,2], RNA [3], and proteins [4], as well as whole cells from a heterogeneous mixture [5]. Upon conjugation of a magnetic microparticle or particles to a target entity, applications of magnetic fields and field gradients are used to isolate a target from a mixture for the purpose of detecting disease, for example by isolating circulating tumor cells [6], RNA from H1N1 viruses [7], and RNA from SARS-CoV-2 [8], including in wastewater samples [9,10]. Furthermore, applying magnetic fields to particles in a fluid environment can be used to treat wastewater [11], to modify physical characteristics of particle populations such as inter-particle spacing [12], as well as to guide fluids in a lab-on-chip environment [13]. Also of interest is targeted drug delivery at specified locations in the body while minimizing unwanted exposure, done by employing magnetic implants to concentrate magnetic nanoparticles and deliver drugs at desired locations [14,15,16]. Additionally, fluid-borne magnetic disks have shown promise as a means of targeting cancer cells in the body via localized heating (i.e., magnetic hyperthermia), delivering antitumor drugs, or mechanical disruption of cancer cells (i.e., by applying forces and torques) [17]. The above uses of magnetic particles involve manipulation via magnetic interactions, necessitating a thorough understanding of magnetic fields and the field gradients that are the cause of these forces.

Traditionally, the manipulation necessary for isolating or applying forces to magnetic microparticles is performed using handheld and otherwise macroscopic magnets, but it has been shown that interactions between microscopic thin-film features in magnetic materials are sufficient to capture fluid-borne particles, using both discrete [18,19] and continuous (i.e., wire) [20,21,22] patterned elements. These trapping and transport mechanisms are advantageous because of their precision, as trap sizes are governed by the geometry of the patterned materials, as well as their scalability, as modern lithography techniques allow for many traps to be engineered on a surface.

By applying and varying magnetic fields, the strength of magnetic forces can be tuned carefully, allowing for the manipulation of magnetic particles along pre-programmed routes on a surface [23,24,25,26,27,28], the tuning of trap strengths by external fields [18], the tuning of particle velocities by applying external fields [29], and the varying of Brownian fluctuations by small changes in external field strengths [30]. As small changes in fields can lead to large ramifications in particle behavior, to understand and successfully predict particle motion, having robust mathematical models of magnetic fields from magnetic thin-film devices is essential.

Analytical models exist and are well-understood for disks uniformly magnetized in a direction perpendicular to the top/bottom face [31,32]; however, for the patterned disks discussed in this paper, the width-to-thickness aspect ratio is sufficiently large to restrict magnetization to the plane of the film. Additionally, for in-plane magnetized disks, semi-analytical techniques are employed to calculate magnetic fields when the magnetization is uniform [33,34]. However, at low externally applied fields, disk magnetization is far from uniform, so other methods are required for calculating disk fields.

In this study, we investigate the local magnetic fields produced by permalloy (Ni_0.8Fe_0.2) magnetic disks. In particular, as it has been shown that these disks are capable of microparticle transport at low external fields—less than 100 Oe—we investigate magnetic field profiles from disks in the presence of low applied fields and as such, we require the use of a realistic magnetization landscape in our models. We employ simulations based on the Object Orientated Micromagnetic Framework (OOMMF) [35] to calculate local stray fields from the magnetic disks. While micromagnetic simulations such as OOMMF are commonly used to describe particles trapped and transported by thin-film structures [27,28,29,36,37,38,39], we focus on the variations that arise from repeated OOMMF simulation runs and notice differences that arise in resulting disk fields. Additionally, as methods based on the OOMMF rely on 2D summation over a discretely defined magnetization landscape for an individual field calculation and are, hence, computationally slow, we discuss a fitting method for faster computation of fields.

2. Materials and Methods

2.1. Particle Trapping and Transport and Effects of Externally Applied Magnetic Fields

Upon magnetizing disks by the application of a magnetic field in an in-plane direction (i.e., to the right in Figure 1a), traps are formed at opposite ends of magnetic disks—magnetizable particles such as superparamagnetic microparticles, or “beads”, are capable of being trapped at one of two sides of the disk (i.e., the north or south pole), as seen in Figure 1c. Upon the application of an out-of-plane magnetic field $H_{ext,z}$, preference is given to one pole, as seen in Figure 1d, while the other is repulsive. To illustrate this, we consider a superparamagnetic bead that is linearly magnetizable by an external field such that its magnetic moment ${\overrightarrow{m}}_{bead}=\chi V\overrightarrow{H}$, where $\chi$ is the particle’s magnetic susceptibility, $V$ is the particle’s volume, and $\overrightarrow{H}$ is the net magnetic field the particle is in. The resulting force on the bead from a spatially varying magnetic field is [19]

$$\overrightarrow{F}=\frac{1}{2}{\mu }_0\chi V\overrightarrow{\nabla }{H}^2$$

(1)

and since generally force $\overrightarrow{F}=-\overrightarrow{\nabla }U$, the resulting potential energy $U$, then, is given by the equation below.

$$U=-\frac{1}{2}{\mu }_0\chi V{H}^2$$

(2)

We note that the potential energy is proportional to the negative square of the field $\overrightarrow{H}$ where $\overrightarrow{H}$ is defined as the net magnetic field from all sources—in this study, the sum of fields from disks (${\overrightarrow{H}}_{disk}$) and externally applied fields (${\overrightarrow{H}}_{ext}$). In Figure 1b we plot $-H^2$ vs. position along a line parallel to the disk diameter at a height of 1.4 µm, the radius of a commonly used magnetic bead, above the disk center. For a magnetic disk in the presence of an externally applied 35 Oe, rightward facing field, the plot shows two energy minima—one on each edge of the disk, shown in blue in Figure 1b. For a disk in the presence of a 35 Oe in-plane and 50 Oe out-of-plane field, one energy minimum is reduced, suggesting a stronger trap at the rightmost disk edge, and the other energy minimum is inverted, suggesting a repulsive disk edge. This effect is reversed if $H_{ext,z}$ is down instead of up.

For this study, we exclusively focus on permalloy disks of a particular size: 5 μm diameter, 40 nm thickness. As permalloy is a soft magnetic material, the application of a magnetic field is required to produce magnetic traps—when fields are turned off, trapped particles are released from traps on disk peripheries and are observed to be carried away from disks by Brownian fluctuations or even gentle fluid flow.

While the externally applied magnetic fields are easy to measure, it is difficult to measure small variations in fields near the surface due to the patterned entities. This study focuses on methods for calculating these local “stray” magnetic fields from these disks (${\overrightarrow{H}}_{disk}$). We also explore variations that arise due to differences in micromagnetic simulations that guide the field calculations.

2.2. Methods for Field Calculations

Magnetic fields from disks were calculated first by simulating micromagnetic details using the OOMMF, then by calculating fields a defined point away from the disk either by summing dipole fields or summing Coulomb fields from a bound magnetic charge density distribution. As an example of an OOMMF simulation used, Figure 2a shows a 5 μm diameter, 40 nm thick permalloy disk in an external field configuration ($H_{ext,x}=$ 35 Oe, $H_{ext,z}=$ 50 Oe) used for particle manipulation [40]. The figure shows a representation of the varying magnetization directions of portions, or “cells”, of the disk generated by the OOMMF simulation. We refer to this spatially varying disk magnetization as a “magnetization landscape”, and it is indicated by the arrows in Figure 2a. In this case, the magnetization landscape appears as a vortex pattern. In the figure, the x- and y-axes are defined. (The z-axis is normal to the disk plane). We place the disk such that the center of the bottom of the disk is at the origin of our coordinate system, and we define a vector $\overrightarrow{r}$ that describes the location of the field being calculated. Furthermore, we define a vector, ${\overrightarrow{r}}^{\prime }$, from a cell of the micromagnetic simulation to the point in space defined by $\overrightarrow{r}$. This vector ${\overrightarrow{r}}^{\prime }$ is varied as we sum over these cells when calculating fields.

We used the default parameters from the OOMMF for permalloy for our simulation: saturation magnetization $M_S$ = 8 × 10⁵ amps/meter, the exchange constant of $A$ = 13 × 10⁻¹² joules/meter, a damping coefficient of 0.5, and no effect from crystalline anisotropy (i.e., $K1$ = 0 joules/meter³). The simulation uses a cell size of 40 nm × 40 nm × 40 nm (arrows in the figure are averaged over several cells for display purposes) and each dipole has a vector magnetic moment of $\overrightarrow{m}$, where the moment’s magnitude is constant (and equals $M_{S}V$, where $M_S$ is permalloy’s saturation magnetization and $V$ is the volume of the cell, in this case, (40 nm) × (40 nm) × (40 nm)).

By summing the field contribution of each dipole, the stray field ${\overrightarrow{H}}_{disk}$ at position $\overrightarrow{r}$ can be calculated by the equation for a field from a dipole:

$${\overrightarrow{H}}_{disk}\left(\overrightarrow{r}\right)=\frac{1}{4\pi }\sum _{\begin{array}{c}all\\ cells\end{array}}\left[\frac{3\left({\overrightarrow{m}}_i\cdot {\overrightarrow{r}}^{\prime }\right){\overrightarrow{r}}^{\prime }}{{\left|{\overrightarrow{r}}^{\prime }\right|}^5}-\frac{{\overrightarrow{m}}_i}{{\left|{\overrightarrow{r}}^{\prime }\right|}^3}\right]$$

(3)

where ${\widehat{r}}^{\prime }$ is the vector from a cell to the point at which the field is calculated, as seen in Figure 2b, and ${\overrightarrow{m}}_i$ is the moment of the ith cell.

Similarly, the field can be calculated by considering a bound charge density, ${\rho }_m$, that arises when the magnetization ($\overrightarrow{M}$) has a nonzero divergence:

$${\rho }_m=-\overrightarrow{\nabla }\cdot \overrightarrow{M}=-\frac{1}{V}\left(\overrightarrow{\nabla }\cdot \overrightarrow{m}\right)$$

(4)

noting that magnetization is related to the magnetic dipoles by $\overrightarrow{M}=\overrightarrow{m}/V$. This divergence in the magnetic moment is approximated using nearest neighbor cells.

Figure 3a shows the whereabouts of the bound magnetic charge density used for field calculations. What results is a charge distribution that resides primarily near the disk periphery—white (positive bound charge density) on the right side of the disk, black (negative bound charge density) on the left, and gray (comparably little bound charge density) throughout most of the disk. While most bound charge density resides on the periphery, some bound charge density resides in the vicinity of the vortex center and stretches toward the periphery. These deviations caused by the vortex structure—that not all charge is on the periphery as it would be for a fully magnetized disk—are blurred when examining the resulting magnetic field 1.4 microns above the disk, as seen in Figure 3b.

Noting that the charge in each cell (indexed as $i$) is calculated as $q_i={\rho }_{m}V$, fields from these charges can be calculated using Coulomb’s law:

$${\overrightarrow{H}}_{disk}\left(\overrightarrow{r}\right)=\frac{1}{4\pi }\sum _{\begin{array}{c}all\\ cells\end{array}}\frac{q_i{\overrightarrow{r_i}}^{\prime }}{{\left|{\overrightarrow{r_i}}^{\prime }\right|}^3}$$

(5)

where $q_i$ is the calculated bound magnetic charge at the $i^{\mathrm{t}\mathrm{h}}$ cell.

Because these two methods for calculating fields involve summing over all cells each time the field at a particular point is to be calculated, each field calculation requires a large computational time, which is unwanted for the purpose of, for example, dynamic simulations meant to simulate particle motion [40]. To address this, we use a fit that can calculate fields more quickly. To do this, we calculate the magnetic scalar potential ${\psi }_m$ [41], valid as we have no free currents, from which magnetic fields can be obtained using $\overrightarrow{H}=-\overrightarrow{\nabla }{\psi }_m$. We calculate ${\psi }_m$ using the following formula.

$${\psi }_m=\frac{1}{4\pi }\sum _{\begin{array}{c}all\\ cells\end{array}}\frac{q_i}{\left|{\overrightarrow{r}}^{\prime }\right|}$$

(6)

We calculate the potential at 276 points (at a height $z$ above the disk) from a point directly above the disk center to a distance of 27.5 µm (or 11 disk radii away). The spacing between points is 100 nm. To these points, we fit the following function:

$${\psi }_{1Dfit}(x,z)=\frac{ax}{d{x}^3+c{x}^2+bx+1}$$

(7)

such that the sum of least squares of the difference was minimized and parameters $a$, $b$, $c$, and $d$ are calculated. This fit function was chosen to ensure that the potential was zero in the middle (at $x=0$) and to have the potential function tend like 1/(distance squared) at distances far from the disk, consistent with the behavior of potential from a dipole. As this function could not perfectly fit the calculated potential, a piecewise function was employed, where in all intervals the same function ${\psi }_{1Dfit}\left(x,z\right)$ above was used but with different parameters $a$, $b$, $c$, and $d$. In our case, for a 5 µm diameter disk, the intervals were 0 to 1.875 µm, 1.875 to 5 µm, 5 to 12.5 µm, 12.5 to 20 µm, and 20 to 27.5 µm.

From this fit potential function ${\psi }_{1Dfit}\left(x\right)$, representing a fit to the potential along the x-axis, we calculate the potential at any location in space ${\psi }_{fit}\left(\overrightarrow{r}\right)$ by assuming that the magnetic scalar potential varies with the cosine of the azimuthal angle (defined as $y/x$)—in other words,

$${\psi }_{fit}\left(\overrightarrow{r}\right)={\psi }_{fit}\left(x,y,z\right)=\left[{\psi }_{1Dfit}\left(x,z\right)\right]\cos \theta =\left[{\psi }_{1Dfit}\left(x,z\right)\right]\left(\frac{y}{x}\right)$$

(8)

It is worth noting that, unlike a 100% magnetized disk, the disk properties do not vary perfectly with the cosine of the azimuthal angle. However, in Figure 4 below, we show that both (a) the magnetic scalar potential and (b) the z-component of the magnetic field, both calculated using the charge density method without using a fit, vary with angle much like a cosine curve, suggesting that fitting calculated potentials to a cosine curve is a reasonable approach.

3. Results

3.1. Field Calculation Results

In Figure 5, we show the calculated field from the disk with magnetization shown in Figure 2a, where the fields shown represent the x- and z-components of the field a height of 1.4 um above the disk, along a line parallel to the diameter. The dipole and charge density methods yield fields that differ by <12% (compared to the maximum field in the graph). We suspect the field differences that occur are a result of pixelation effects—the OOMMF models a disk as broken into 40 nm cube cells. Note that the fields shown are the field contributions from the disks themselves and do not include the externally applied fields that magnetize the disk.

Interestingly, the maximum field magnitudes illustrated are only about 12 Oe, yet these fields (and the resulting field gradients) are strong enough to modify a magnetic field landscape such that magnetic microparticles in fluids can be trapped and transported [40].

It is worth noting that the charge distribution method is a bit computationally faster than the dipole method, and the fit method is much faster than the other two. To illustrate this, using a desktop computer with 16.0 GB of RAM and a 3.20 GHz processor, it took 4.53 s to produce the dipole field graph in Figure 5b—121 field calculations in all—for an average field calculation of 37.4 ms. For comparison, the charge distribution curve took 2.28 s (18.8 ms/field calculation) and the fit graph took 4.5 ms (37 µs/field calculation). Calculations necessary for Figure 5a required comparable computation times. All computations were performed using Python 3.8.10.

3.2. Variation in Stray Fields across Simulations

Field strengths vary dramatically for different disk magnetization landscapes, and OOMMF simulations can yield different results for different simulation runs. For comparison, we discuss three different magnetization configurations: (a) a disk initially randomly magnetized and allowed to relax to an applied 35 Oe in-plane field/50 Oe out-of-plane field, followed by a 360-degree counter-clockwise rotation of the field (in 90-degree increments), resulting in a “vortex” configuration, (b) a disk initially randomly magnetized and then allowed to relax when experiencing an applied 35 Oe in-plane field/50 Oe out-of-plane field, but without the subsequent rotation, denoted “wavy” because of the wavy resulting magnetization landscape, and (c) a fully magnetized disk—denoted “full mag”—with all cell magnetizations pointing horizontally, reminiscent of a disk magnetized by an arbitrarily large in-plane field. (Note that case (a), the vortex, was already used to produce Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5).

While the fully magnetized case is not physically realistic for the small ($H_{ext,x}$ = 35 Oe and $H_{ext,z}$ = 50 Oe) fields employed for the simulation, it was incorporated for the purpose of comparison: a fully magnetized disk would yield the highest possible stray fields.

Our in silico results show that the disk stray fields varied to a large extent even for configurations (a) and (b). This was the case even though the final external field employed by the OOMMF simulation was the same—different initial conditions but identical final conditions yielded maximum fields that differed by almost a factor of 2. (The maximum $H_{disk,z}$ in Figure 6e was 24.8 Oe for the wavy case and 13.2 Oe for the vortex case.) The result of the field rotation producing Figure 6a was the removal of energetically unfavorable domain walls and other locations where the magnetization had large divergences, resulting in a vortex void of all frustration points except the single vortex center. It is important to note that for magnetic particle transport on disks [40], fields are rotated as part of the transport process, so a configuration like 6a is a more likely configuration than what is shown in Figure 6b. Verifying this field calculation model would be possible by measuring fluid-borne magnetic particle responses—particle speed measurements, for example—upon application of trapping forces and manipulation forces from disks.

It is important to note that upon randomly initializing OOMMF simulations, different magnetization landscapes emerged even for identical final experimental magnetic field conditions. To characterize disk magnetization, we recorded “$M_x/M_S$”, a metric for how magnetized the resulting disk described by the OOMMF is. (In reality, “$M_x/M_S$” is the average of all $M_x$ values in a magnetization landscape divided by the saturation magnetization $M_S$). As an example, $M_x/M_S=$ 0.308 for Figure 6a, 0.533 for Figure 6b, and 1 for Figure 6c (by definition). Stray field maxima increase with $M_x/M_S$, as seen in Figure 7. It is worth noting the variation in simulation results even when the setup of the simulation was identical. All blue data points represent identical OOMMF simulations (except for the random magnetization initialization), and likewise for orange data points. Despite this, the standard deviation of $M_x/M_S$ values for “identical” simulations was rather large. For 10–15 simulations, $M_x/M_S$ = 0.286 ± 0.092 for the simulations with rotation and 0.536 ± 0.056 for simulations without rotation. The magnetization landscapes in plots 6(a) and (b) were chosen as representative cases near these averages.

4. Discussion

To summarize, we explore methods for calculating magnetic fields in the vicinity of permalloy disks by using micromagnetic simulations. Regardless of whether fields were calculated directly from dipoles within the material, from calculated bound charge distributions, or from the magnetic scalar potential fit method described, reasonably similar field calculations emerge. However, for the weak applied magnetic fields investigated here and used experimentally to trap and transport magnetic microparticles [40], small or seemingly insignificant differences in OOMMF simulations can yield drastic changes in resulting calculated disk magnetization landscapes and stray fields. The micromagnetic simulations studied suggest that more than just a final external field configuration is necessary for obtaining a reliable understanding of the micromagnetic properties of a disk—the steps taken to get to that result influence the result as well.