A comparison of numerical simulations and analytical theory of the dynamics of interacting magnetic vortices

Magnetic vortices, which involve the curling of the magnetic spins into a circular pattern with a central core that points out-of-plane, are often found in micron- and sub-micron-sized structures that are magnetically soft.^1,2 Vortices can interact dynamically with one another when confined in the same structure³ or in structures that are in close proximity, where coupling in the latter case occurs via magnetostatic interactions.^4–7 Modifications of the gyrotropic resonance frequency, first predicted to arise in a 2D square array of magnetic vortices,^8–10 have been observed experimentally for 2D square arrays^4,11,12 and in the simpler case of a pair of closely spaced magnetic squares.^5,13 In the case of the 2D arrays, the effect is primarily a line-broadening effect of ∼20% for an inter-disk spacing of 20% of the disk radius,⁴ whereas the gyrotropic resonance frequency splits into two distinct modes for a vortex pair, the frequencies of which depend on the vortex polarities.⁵ The frequency splitting tends to be more pronounced when the polarities are opposite to one another. Coupling effects are important for a variety of applications including vortex-based magnonics,⁸ for increasing the signal from vortex-based spin torque oscillators,^14,15 and they may also lead to new devices since it was recently shown that the dependence of the resonance frequencies on the polarities can be used to dynamically control and read the state of coupled vortices.¹⁶ Furthermore, recent studies have shown that signals can be transmitted from one disk to another via dipolar coupling in chains of two,^5,13,17 three,¹⁸ or more structures lined up in a row,¹⁹ and schemes to improve the signal transfer rate²⁰ and coupling efficiency in vortex chains have been investigated.²¹ Several theories of vortex coupling have been proposed^8–10,22,23 but have not yet been compared against full micromagnetic simulations.

Here, we conduct a detailed comparison of analytical theory of vortex dynamics including magnetostatic interaction effects with full micromagnetic simulations to determine the extent that the analytical expressions can be used to predict the behavior of interacting vortices. We use a geometry involving three disks arranged to form an equilateral triangle for this comparison, a configuration that is simple enough that it can be treated analytically and realized experimentally but complex enough to serve as a useful test. The magnetization signals are found to exhibit beating patterns in the time domain that vary with the polarities p and chiralities χ; however, the resonance frequencies only depend on whether all p are the same or whether one p differs from the others. The qualitative aspects of the simulations agree well with even the simplest analytical calculations of the dynamics based on purely dipolar interactions; however, more exact models of the dynamic vortex spin distributions will be required to accurately predict the magnitude of the splitting. Our recent time-resolved photo-emission electron microscopy results verify that the dipolar interactions between the vortices are strong enough to cause a splitting in the frequency spectrum that is consistent with the theory and simulation results.²⁴ 

The dynamic behavior of magnetic vortices in a system consisting of three Permalloy disks arranged in an equilateral triangle (inset of Fig. 1) was simulated using the Object Oriented Micromagnetic Framework (OOMMF).²⁵ Disks of diameter $2R=252$ nm were represented using a 2D mesh with cells of 4 × 4 nm² and thickness $t=20$ nm. Material parameters typical for bulk Permalloy were used: exchange stiffness constant of $A=1.3×10−11$ J/m, saturation magnetization $Ms=8.6×105$ A/m, gyromagnetic ratio $γ=1.76×1011 s–1 T–1$⁠, and the magnetocrystalline anisotropy was neglected. A variety of vortex χ and p combinations were considered for a fixed radius R to inter-disk center-to-center separation D ratio of $R/D=0.45$⁠. The dependence of the mode frequencies on D was simulated for selected χ and p combinations for $R/D$ from 0.20 to 0.49, which covers the range from where the disks are far enough apart that the magnetostatic interactions are small to where the disks are almost touching. The magnetization was first allowed to relax to an equilibrium configuration in the presence of a static field of $μoH=5$ mT applied along the direction connecting two disks centers (x-direction) using a large damping parameter $α=1$ for fast convergence, and then the static field was abruptly removed and the dynamics were modeled using a realistic $α=0.01$⁠. Additional simulations were conducted for selected configurations with $α=0.001$ to construct spatial profiles of the individual eigenmodes.

Fig. 1 shows the time-evolution of the magnetization of each disk and the corresponding Fourier spectra for a variety of p combinations and fixed, identical χ for $R/D=0.45$⁠. The details of the temporal responses vary with the choice of p; however, two main types of responses were observed: those associated with disks of all the same p and those with one p reversed. When χ and p are all the same (Figs. 1(a) and 1(e)), two frequencies are observed: a dominant resonance at 630 MHz and a weaker resonance at 690 MHz, which both differ from the single vortex gyrotropic frequency of 680 MHz (not shown). The same spectral response and overall temporal characteristics are present for other χ combinations as well (not shown). In contrast, when one p is opposite to the others, three distinct resonant frequencies are observed at 600, 660, and 750 MHz (Figs. 1(b)–1(d) and 1(f)–1(h)) and pronounced beating patterns are detected in the time domain. The phases and the relative mode amplitudes and hence the details of the temporal patterns depend on the exact p and χ configuration, but the frequencies depend only on the relative p of the three cores.

Figure 2 shows how resonance frequencies change as a function of D. When all vortices have the same p (Fig. 2(a)), the resonance splits into two different modes as the disks are brought closer together, whereas when one p is opposite to the others, three different frequencies are observed (Fig. 2(b)). In both cases, the magnitude of the splitting increases with D up to a maximum shift of 22% (R/D = 0.49). When the disks are further apart (R/D = 0.20), the vortices oscillate at the single disk gyrotropic frequency of 680 MHz.

From an analytical perspective, the Thiele equation, which describes the dynamics of a magnetic vortex in terms of the position of its core, serves as the basis for most models of vortex dynamics.²⁶ Theoretical models have been developed to describe the expected modes for vortices coupled via the magnetostatic fields that are generated by the moving vortices for a variety of geometries. Calculations have been conducted for vortices in arranged in a 2D square array,⁸ for the normal modes for two disks that are excited together,²³ and a more general treatment examines vortex pairs as well as 2D square and hexagonal arrays.²² The main focus of the previous work has been on the calculation of interaction integrals for specific models of the vortex spin distribution, e.g., the poles-free model.²² Here, we derive the theory of the expected modes for the triangular arrangement considered in the simulations using a simple dipole–dipole interaction model and compare the eigenfrequencies and eigenvalues to the simulation results.

For small displacements of a vortex $x→$ from its central equilibrium position, there is a restoring force $F→ =−k x→$ on the vortex that is mainly due to the increase in the demagnetization energy. Interplay between the restoring and gyrotropic forces gives rise to spiraling motion with a characteristic resonance frequency $ωo=k/G$⁠, where $G=2πtMs/γ$ is the gyrotropic constant that is independent of $p=±1$ and $χ=±1$⁠, where $χ=+1$ indicates a counterclockwise vortex. Reference 27 gives an expression for $ωo$ in terms of the magnetic and geometrical parameters of the disk. In the limit of small aspect ratio $β=t/R ≪1$⁠, $ωo=59πγμoMsβ$⁠, where $μo$ is the magnetic permeability of free space, which corresponds to $kapprox=10 μoMs2t29R$⁠. This estimate is reasonable for $β<0.1$⁠; for the geometry used here, $β=0.16$⁠, so we have calculated $k$ using the full expressions given in Ref. 27, which yields $k=0.81⋅kapprox$⁠.

Magnetic vortices in adjacent disks will be coupled via magnetostatic interactions, and the dominant term in a multi-pole expansion will be the dipole–dipole interaction term. The interactions are generally only significant when the space between the disks is considerably smaller than the disk diameter, and in this case, higher order terms in a multipole expansion are needed to fully describe the magnetostatic interactions.²² Nevertheless, as we will show below, the expressions obtained using just the dipole–dipole interactions are sufficient to describe the eigenfrequencies as a function of the inter-disk separation and the eigenvalues. For small displacements, the dipole moment $d→$ of a shifted vortex can be approximated as

where dimensionless factor $C$ is a function of $t/R$⁠. The dipole–dipole interaction energy between the disks is

where $D̂ij$ is the unit vector pointing from the center of disk i to the center of disk j, which can be written in terms of the vortex displacements as

where $k2=C2μoMs2t2R2/4πD3$ for circular disks. For two disks arranged along the x-axis, this reduces to $Eint=χ1χ2(k2x1x2−2k2y1y2)$⁠. This can be written more generally in the form $Eint=χ1χ2(nxx1x2+nyy1y2+O(| x→j |3))$⁠, as in Ref. 22, where $nx$ and $ny$ are interaction coefficients that can be calculated for not just dipolar but also higher order magnetostatic interaction effects. For purely dipolar interactions, $nx=k2$ and $ny=−2k2$⁠. Methods for calculating the coefficients $nx$ and $ny$ for two specific theoretical models of the vortex spin distribution, the rigid vortex model and the poles-free model, are given in Ref. 22 (see Eq. (13)) where the interaction energies are integrated over the full vortex magnetization distributions. Regardless of the form of the interaction energy, the Lagrangian for an ensemble of disks is given by

We can remove the chiralities by switching variables to $x⇀i→χix⇀i$⁠, which means that the frequencies are independent of $χi$⁠. They will, however, depend on p.

As mentioned above, there are only two possible spectra based on the vortex polarities. Without loss of generality, we shall set $p1=p2=1$⁠, and leave $p3=p$ to distinguish the two possibilities. Setting disks 1 and 2 along the x-axis, the magnetostatic energy of the system is

If we define new coordinates $x¯=(x1+x2)/2$⁠, $y¯=(y1+y2)/2$⁠, $d=(x1−x2)/2$⁠, and $Δ=(y1−y2)/2$⁠, then the equations of motion are given by

In the form where only dipolar coupling is considered, the equations of motion reduce to

The eigenfrequencies and eigenvalues can be obtained numerically from Eq. (6) or (7). These expressions are valid for any structure shape, not just circular disks, provided that the values for $k2$ (or $nx$ and $ny$⁠), which can be derived for just a pair of disks, and $k$⁠, which pertains to a single structure are obtained for the geometry of interest. These parameters do depend on the structure shape and dimensions and can be found either through the evaluation of appropriate interaction integrals or from simulations of just one or two disks. Furthermore, the method used to derive (6) and (7) can be easily adapted to any geometrical arrangement of an arbitrary number of structures.

The dependence of the eigenfrequencies of (6) on $k2/k$ are shown in Figs. 2(c) and 2(d), where $k2k=9C240π(RD)3$ (valid for circular disks). Note that the solutions are not all real for all $k2/k$⁠. For $k2/k>2/5$⁠, at least one frequency is imaginary, suggesting that the vortices are destabilized. In Figs. 2(c) and 2(d), the solid green line shows the dipole–dipole interaction results using a value $Cfit=2.82±0.06$ that was extracted from simultaneous fits to the results shown in both plots (where the full expression in Ref. 27 was used to find k). This $Cfit$ yields $k2/k≈0.070$ for the maximum $R/D=0.5$ (disks touching), which is below the instability threshold. Note that the resonance frequencies obtained for R/D = 0.49 were not included in the fits since the gap between the disks is as little as a single cell and edge effects may artificially increase the coupling effects; for $R/D=0.44$⁠, a < 2% difference was observed when the cell size used in the simulations was reduced from $4×4 nm2$ to $2×2 nm2$⁠.

The frequency spectra obtained from the simulations agree well with the theory (Eq. (7)), provided C is used as a fitting parameter (Figs. 2(c) and 2(d)). Equations (6) and (7) predict three distinct modes for both p combinations; however, when all p are the same, the two higher-frequency modes are quite close, which is why only two modes are identified in the simulations. The frequency splitting is more pronounced for $p=−1$⁠. The modes become distinguishable from one another at slightly lower R/D (larger D), and the magnitude of the splitting is larger, which is consistent with reports of stronger coupling for opposing polarities in experiments¹⁹ and simulations.²⁰ 

There are, however, discrepancies in the predicted magnitude of the splitting effects for all of the models considered. In Fig. 2, the parameter $C$ was used as a fitting parameter to set the scale of the dipolar interaction model; however, its value can be estimated from the simulations using Eq. (1) and fits to the core positions and magnetization values. This yields $Cdip=2.0±0.1$⁠, which is smaller than the value obtained from fitting, which was $Cfit=2.82±0.06$⁠, hence the dipole model underestimates the magnitude of the observed splitting. Figs. 2(c) and 2(d) also show the frequency splitting that is predicted from Eq. (6) using $nx$ and $ny$ for the two-vortex and poles-free models from Ref. 22, which include higher order interaction effects. The poles-free model provides closer agreement but still underestimates the splitting, whereas the rigid vortex model overestimates the effect. The functional form of $f/fo$ vs. $(R/D)3$ is similar for all of the models, which occurs because the ratio $nx/ny$ is close to the dipolar ratio of −1/2. The simple dipolar model therefore provides a useful functional form for the splitting; however, an intermediate model for the magnetization profile of a dynamic vortex is needed to correctly predict the magnitude of the effect.

Fig. 3 (Multimedia view) illustrates the normal modes for the system from Eq. (6). Similar results are obtained for the poles-free and rigid vortex models, again because $nx/ny∼−1/2$⁠. Figs. 3(a), 3(c), and 3(e) show the modes for $p=1$⁠. For the lowest frequency mode (Fig. 3(a)), the cores all have close to the same phase at any given time and attractive interactions occur whenever the dipoles are aligned along an axis that connects two disk centers. A spatially uniform field pulse applied to disks with the same χ provides a phase relationship between the cores that is similar to this mode, hence this is the dominant mode in Fig. 2(a). When the core in the top disk is reversed (Fig. 3(b)), there are still attractive interactions for lowest frequency mode but the phases are timed such that when the north or south pole of the top disk is at its lowest point it is attracted equally to both of the lower disks.

In contrast, the higher frequency modes involve mainly repulsive interactions or combinations of repulsion and attraction. Repulsive interactions occur for the two higher frequency modes for $p=1$ (Figs. 3(c) and 3(e)) but at different phases for each. For $p=−1$⁠, the middle mode (Fig. 3(d)) involves motion of mainly the two bottom vortices and is characterized by both attractive and repulsive interactions, whereas the highest frequency mode (Fig. 3(f)) mainly involves repulsive interactions that occur when like dipoles are aligned along an axis that connects two disk centers.

Spatial profiles of the normal modes (Fig. 4 shows modes for $p=−1$⁠) were obtained from the simulations by calculating the Fourier transform of the z-component of the magnetization for each cell, as described in Ref. 28. Plots of $A cos ϕ$⁠, where $A$ and $ϕ$ are the amplitude and phase, respectively, calculated from the simulations for R/D = 0.45 provide a snapshot of the core positions in each disk. The phase has been shifted so that white (black) coincides with the positive (negative) pole of the dipole in the theoretical modes shown in Figs. 3(b), 3(d), and 3(f) to show that the relative phases observed in the simulations agree well with the eigenmodes of Eq. (6). Similar agreement is observed for the low frequency $p=1$ mode (Fig. 3(a)), but it is difficult to obtain reliable mode maps for the two higher frequency $p=1$ modes because they are so close in frequency.

In conclusion, we have compared results from micromagnetic simulations and analytical theory for a system of interacting vortices in three magnetic disks arranged in an equilateral triangle. The theory and simulations agree on the main features of the observed dynamics. The vortex gyrotropic mode splits into three modes, the magnitude of the frequency splitting increases with decreasing D, and only two distinct frequency splitting behaviors are observed that depend on whether all p are the same or one p opposes the others, where the latter leads to stronger splitting. A simple model involving only dipolar interactions accurately captures the qualitative features of the coupled vortex dynamics including the frequencies and eigenmodes; however, a more detailed model of the dynamic vortex spin distribution will be required to improve predictions of the magnitude of the frequency splitting. Nevertheless, the simple dipole-dipole model, once calibrated for a particular disk geometry, provides an excellent description of the dynamics of interacting vortices that can be extended to more complex arrangements of magnetic vortices, including square or hexagonal arrays and chains of arbitrary geometries, which will be helpful for the development and optimization of new vortex-based devices.

This work was supported by the NSF under Grant Nos. NSF DMR-0907706 (CSU), NSF DMR-1053854 (Bryn Mawr), and NSF DMR-1104753 (Johns Hopkins); work at the Advanced Photon Source was supported by the U.S. DOE-BES, Contract No. DE-AC02-06CH11357.