Vortex lines in a cubic magnetic nanodot: structure and dynamics

. Langevin simulations of cubic magnetic nanodots were performed using the Landau-Lifshitz equation with exchange and dipolar interactions. Vortices tend to organize as lines: we establish the structure and dynamics thereof for a large range of the dipolar versus exchange ratio d . These lines tend to be bent and twisted. For large values of the dipolar interaction, a complex network of vortex lines arises. Dynamics evidences low frequency collective gyrotropic motions of vortex lines which maintain their distance during motion.


Introduction
A tremendous amount of both experimental an theoretical work is being done on magnetic nanodots for obvious practical purposes such as data storage and manipulation, but also from a fundamental viewpoint for their often challenging behavior.
Since the first observations of Weiss domains [1], many investigations on the nature of these complex magnetic structures [2] and of their dynamical properties [3] were carried out. Magnetic nanodots, through increased performance and miniaturization, follow the same trend as electronic devices with Moore's law. After magnetic tapes, ultrathin quasi-2D films were produced [4], then magnetic vortices were predicted theoretically to occur under the competition between exchange and dipolar interactions [5,6], and observed experimentally in nanodisks [7,8]. New kinds of magnetic memories linked with the polarity of these nanodisks were introduced [9]. Studies with different nearly 2D nanodots showed that under the influence of anisotropic Dzyaloshinsky-Moriya interaction, magnetic structures rather similar to asymmetric vortices are predicted [10], observed and called skyrmions [11], the advantage of skyrmions with respect to standard vortices being their reduced size [12].
The competition between long-ranged interactions such as dipolar ones, and short-ranged interactions such as exchange leads to the domain organization of magnetic samples [2]. Shape and size effects result from this competition. The shape effect, also called magnetostatic effect, has been long well-known in compasses. More recently, miniaturized needles are moved by means of distant magnetic fields, with mechanical or thermal effects for instance applied to dealing with cancer tumors [19]. The size and shape of elongated nanofibers can be controlled by means of chemistry in polyols [20][21][22] and by means of electrochemistry with further thermal treatment [23], with many applications. A new application of such electro-chemical * e-mail: depondt@insp.jussieu.fr treatments [23] enlarges these fibers in some places by thermal heating to form blobs. Such a connected set is similar to a set of neurons and axons : a useful analogy in the now very active field of neuromorphic spintronics applications [24].
The dynamics of these many new magnetic structures is obviously of interest. The study of magnetic excitations, spin waves in thin films, which were used in magnetic tapes, soon revealed a specific behaviour among other excitations: a wavelength dependent velocity [25]. This result introduced further questions on the possibility of 2D magnetism [26]. Now, the introduction of complex magnetic material sets open large possibilities for finding new excitation spectra and thus metamaterials. An original spectrum is observed for fractal sets when only exchange and a saturating magnetic field occur: in these connected samples, a very complex spectral structure with many gaps is observed [27][28][29]. These gaps and their spectral location strongly depend on the specifics of this potentially rich lacunary structure. The introduction of dipolar interactions in such magnetic structures leads even in standard samples to new specific modes which occur at rather low frequency, in the terahertz domain as early observed from Brillouin scattering [30,31]. With such possibilities of applications as metamaterials in this very interesting spectral domain, magnonics, the observation of spectral features in such new materials, received a lot of attention [32][33][34][35]. A great advantage of extended structures effectively submitted to dipolar interactions is the coupling of physically disconnected parts. It opens the path many studies in this very rich field of lacunary systems. This choice of structures considerably enlarges the field and so the possibility of searching for new metamaterials where gaps are crucial [36]: for instance, special orderings such as fractal ones [33] or orders induced by Fibonacci series [37] were recently considered with evidence for gaps in the spectrum.
of Langevin dynamics at low temperature, i.e. far from Curie and Kosterlitz-Thouless transitions, in presence of exchange and dipolar interactions. In this basic study a large variation of the ratio dipolar interaction versus exchange is considered. This trick enables us to consider different materials, and also to simulate different effective sizes from our study of a (64 × 64 × 64) cube. The use of this cube as a mesh opens the path for a full study from an atomic nanocube up to a very large cubic sample, analysing in this way the transition from a single domain without vortex, up to a full set of vortices and domains, both on the viewpoints of structure and dynamics, when taking into account the complete screening effect of a complex magnetic structure. Of course, this arrangement of domains and vortices is expected to occur [38][39][40], but the perfect lack of impurities and defects would be impossible to observe experimentally. So, this numerical experiment which avoids the nuisance of impurities brings a new paradigmatic brick to magnetism. A later comparison with experimental results, even with defects, sounds attainable since recently a high level of resolution has been obtained from electron spin holography [41,42] and even with time resolution [43] in different ways.
The long ranged dipolar interaction has several analogues in physics. Similar situations are issued from solutions of a Coulomb like equation with a linear Dirac like singularity at the origin. In mechanics the interaction between dislocations lies among these analogues. In the case of dislocations, atoms are not rigid and the mobility along dislocations, an essential feature, plays a main part. A strong mobility forbids a comparison with the case of magnetic vortices where matter remains perfectly static. However, in the very special case of hard materials which include numerous useful impurities which stack the possible displacements along dislocations in a quite rigid manner, such a comparison makes sense. And in this case a complex network of dislocations has been recently observed to occur with different crystallites between dislocations [44]. This result is quite comparable to our present observation of differently oriented magnetic domains parting magnetic vortices as an effective topological repulsion. Moreover, the atomic mobility which occurs in elasticity could introduce, in the case of spins, the principle of a synaptic like memory process, an Hebbian link, for neuromorphic spin applications.

Model
The model system consists in a rigid (64 × 64 × 64) simple cubic lattice ( i.e. 2 18 = 262144 spins), on which classical spins are free to rotate following the Landau-Lifshitz equation of motion at finite temperature T [71]: in reduced units and where s is the spin on site , α a damping coefficient and H th is a Langevin-style random thermal field with zero mean and normal distribution such that: is the local field felt by spin . The damping coefficient α is necessary as it makes the system relax toward thermodynamic equilibrium but it should be chosen small enough to have little influence on the dynamics of the system as it will tend to broaden vibrational peaks. The local field is generated by the interactions with all other spins in the system, no external field being introduced in the problem in this work: it is time-dependent through the motion of the other spins. The local field consists here in two terms representing the exchange and dipole-dipole interactions: where J is a positive constant for a ferromagnetic interaction and { , } is a sum restricted to nearest neighbours. The dipolar field writes: where d is the dipole-dipole interaction strength and r , the vector connecting sites and . We choose J = 1 so that d is also the dipole/exchange ratio. Other interaction types could easily be introduced but we concentrate here on the effects of the dipole-dipole interaction.

Integration of the spin precession motion
The instantaneous motion produced by equation (1) is a precession around the effective field H eff which conserves |s|, so that it is adequate to use the Rodrigues rotation matrix [72] to integrate the equation of motion: where δ t is the integration time-step, with Coordinates h x,y,z are those of the unit vector h = H eff / H eff . The so-called "improved Euler" or Heun method [66] is used to integrate equation (1) numerically and the damping coefficient α is chosen small enough not to perturb the dynamics of the system in the relevant frequency range.

Dipole-dipole interaction
Since the dipole-dipole interaction decreases with distance as 1/r 3 , each spin is under the influence of all others and cutoff radii and other such approximations cannot be used. Writing the interaction as a convolution and then using the convolution theorem and fast Fourier transforms significantly accelerates computations [66]. The introduction of dipole-dipole interactions will quite often generate vortices and anti-vortices, the cores thereof must be detected. If we consider one plane (e.g. (x, y) of the cubic lattice, the vortex-cores reside within the squares of the lattice, and arise when the projections of the four neighbouring spins on the the square's plane "wrap around" the square (figure 1). One can thus compute the four angles ϕ i , i ∈ [1, 4] that these projections make with a given axis. In the sum

The search for vortex cores
if Ω equals 2π then a vortex-core is present on that site, while if the sum equals −2π it is an anti-vortex; otherwise there is no such singularity on that site. This search must be done in the three relevant planes (x, y), (y, z) and (z, x) of the cubic lattice: the polarization of the vortex will then be along z, x or y respectively.

Vortex core line detection
The result of such a simulation with vortex-core detection will consist in a list of vortex-core sites which changes with time: in principle, to analyse such complex data, graph-theory based methods [73] could be used, however we adopt here a simpler procedure based on the specific physics of our problem. In three dimensions, vortex-cores are apt to gather up in lines which can have more or less complicated shapes: they can be more or less sinuous or "thick". In practise, these lines are obtained through a neighbour-of-neighbour-of-neighbour type of procedure: one chooses a vortex-core, looks for neighbours, then neighbours of these neighbours, etc., in a recursive manner. The simplest approximation of one vortex line is a straight line (figure 2). The orientation of that straight line will provide the overall orientation of the sinuous vortex core line, while the differences will yield the internal structure and dynamics of the line. A vortex line is defined by n vortex cores at locations: The center of the vortex line is thus The coordinates of a point on the line are therefore: r = r 0 + ξe where ξ is a scalar. We want that line to be as close as possible to the vortex cores. Let ξ i be the value of ξ for which the segment that joins the line to the i th core is normal to the line direction e: where: The squared distance between the set of vortex cores and this line is which should be minimized with the constraint that e is a unit vector: where λ is a Lagrange multiplier. Thus: which is a symmetric eigenproblem which can easily be solved by standard methods. Since this is a (3 × 3) problem, we obtain 3 eigenvalues and as many eigenvectors. The eigenvalues λ, divided by n, are the variances of the distance from x 0 to the actual points projected on the directions defined by the eigenvectors: this defines an ellipsoid, the longest axis of which is the overall orientation e of the line, and the two others the dispersion in the normal directions. This procedure is clearly very similar to finding the inertial tensor of a set of masses [74]. Figure 2 shows a set of points and the three eigenvectors multiplied by the associated standard deviation. Of course the largest eigenvalue gives the overall direction e of the line.
The time evolution of r 0 (t) should yield the translational motion, if any, of the vortex line, that of e(t)), the rotational motion of the line as a whole. The separations d i between the straight line and the actual vortex core positions r i should provide the internal motions of the vortex line.
The separations d i write: These vectors can be expressed in a reference frame (e 1 , e 2 , e 3 ), where e 3 is equal to e et e 1 and e 2 are the two other eigenvectors. This yields: If we want to include the position of the vortex upon the line, we can just as simply write δ i in the (e 1 , e 2 , e 3 ) reference frame: We therefore now have the means to attempt to understand vortex line dynamics.

Simulation parameters
All  (1) was set to α = 0.01 in reduced units which is a compromise between reasonably rapid thermalization and not too much interference with the system's dynamics.
The system was left to relax to whatever configuration it chose and then averages obtained. Figures 3 and 4 show two typical configurations for d = 0.005 and 0.1. Both were obtained after relaxation of identical initial conditions chosen as four straight vortex-core lines parallel to the z axis. The colors indicate the outof-plane component of the spins while the arrows show the in-plane components thereof: vortex-cores are visible through the in-plane components but also through the out-of-plane component: this is a well-known effect of the exchange-dipole competition. The first observation to be made from these figures is that the vortex-core location for each "slice" is approximately the same, meaning that vortex-cores indeed line up (in this case, parallel to z due to initial conditions, but the three axes are of course equivalent). The effect of sample boundaries can also be observed in these figures: spins tend to remain parallel to surfaces and edges resulting in competitions which clearly are not easily solved by the system.  For d = 0.005, the system spontaneously evolves towards a single vortex line configuration with quite a lot of disorder, while for d = 0.1 the four vortex line configuration remains. If the "width" of the line is measured by the extent of the area in which the out-of-plane component is non-zero, one observes, firstly that this width tends to decreases with increasing dipolar interaction, which is no surprise, while it changes quite a lot within the bulk of the sample: narrow on the surface (layers 0 and 63 in both figures), broad in the middle (layers 27 and 36). Flux closure induced by the dipole-dipole interaction implies that the spins on the surface layers tend to lie parallel to that surface: the out-of-plane component therefore is energetically not favourable on surface layers so that the vortex core extent as defined by the out-of-plane component will be reduced near surfaces, and not so within the bulk.

Configurations
We need to address the issue of whether these are equilibrium configurations. Simulations with different initial  Table 1. Effect of the initial configuration choice. Initial configurations with uniform magnetization, one vortex-core line or four vortex-core lines were tested and results in either uniform magnetization, one line or four lines. The last line of the table shows which appears to be the lowest energy configuration and therefore a good candidate for being the ground state.
configurations but otherwise identical were done: figure 5 shows the time dependence of the energy of these simulation runs for a sampling of values of the dipolar to exchange ratio d. In some situations (e.g. d = 0.005) all initial conditions end up producing the same final energy which turns out to correspond to equivalent configurations (e.g. for d = 0.05 all three initial configurations yield a final state with one vortex-core line). Table 1 summarizes the effect of the choice of the initial condition. By comparing figure 5 and  [75]. These results also show that metastable states are easily reached, at low temperature at least, and the notion of ground-state should be taken with caution, because of, among others, surfaces on the one hand, and vortex-core pinning which was observed in two-dimensional simulations [67], on the other.

Vortex cores
Direct observation of the configurations, while instructive, sometimes does not easily yield relevant information, so that we now focus on vortex-cores and forget about other details. Vortex-cores and anti-vortex-cores are detected via equation (3). Figures 6-8 show situations of increasing complexity with increasing dipole-dipole interaction. We can first confirm the previous remarks, that is the presence of one vortex line for small d and four vortex lines for larger d. The single anti-vortex line between the four vortex lines is also visible when applicable. We also confirm that the polarization of the vortex core is parallel to the overall line orientation. However, it now appears clearly that these lines are twisted and that additional features are present: for d = 0.1 and more so for d = 0.4, other shorter vortex lines spontaneously arise in a normal direction: one may surmise that for larger samples three dimensional arrays of such lines may be relevant. For the higher dipoledipole interaction case, it also seems that an effective re-   pulsion between vortex lines is present inducing additional bending of the lines.

Vortex core line dynamics
3.3.1 Line precession Figure 9. For d = 0.005 the above mentionned vortex-cores for which the color indicates time. Figure 9 represents the time evolution of a vortex-core line and shows that a dynamics of these are to be observed: this can be done for all d values. Again, we have a tremendous amount of data, out of which we would like to extract a few useful items. The next step is thus to carry out the analysis explicitly in terms of vortex lines as in section 2.4 and equation (4).
The first issue that must be addressed is the definition of when a vortex core is to be considered a neighbour of another. If we limit the search to nearest neighbours, it means in practise that only vortex-cores situated strictly on top of each other will be considered part of the same vortex line, so that only perfectly strait lines will be accepted. Some limit must then be defined, next nearest neighbours, or beyond: if one decides that the limit should be chosen at a distance less that √ 3, the two vortex-cores should be part of the same cube. This still seems a bit drastic since a vortex-core is a rather large object as can be seen in figures 3 and 4. At this stage, the threshold is set at 4 after gradually increasing it until the vortex-core lines thus detected correspond to those that the human eye will spontaneously establish: more rigorous procedures, e.g. Voronoi tesselation-based, are certainly possible. Figure 10 shows the line dynamics for d = 0.005: a case for which we expect a one-line situation, in which the line is roughly parallel to z. The straight blue (z component of r 0 ) line at z = 32 in the top panel indicates that the z coordinate of the middle of the line lays midway in the sample, but the presence of other points for z 16 and 48 also shows that sometimes the line gets cut somewhere in its middle, which lead to the suspicion that the vortex core density might be lower in the middle: this is consistent with the fact that that the vortex-core sizes tend to be larger in the middle of the sample: the number of detected vortex-core lines thus detected varies in this case between one and two. The x and y components oscillate which confirms lateral motion of the line taken as a whole. The bottom panel of figure 10 shows the time dependence of the coordinates of the line orientation e. Unsurprisingly, since we expect the line orientation to be parallel to z, the z component (in blue) remains close to one. The other two components however oscillate showing a precession motion of the line.   approximately as the logarithm of d as shown in figure 12. This allows to define a characteristic dipole/exchange ratio d 0 0.005 for which the precession vanishes; a characteristic frequency f 0 0.0013 can also thus be defined.

Intra-line dynamics
Using equation (6), one potentially has access to the internal structure and motion of a vortex line. This however reaches beyond the scope of the present paper.

Discussion
Vortex-core motion, usually the gyrotropic motion, is well known, it has been observed and calculated in various situations, and many times interpreted through the generic Thiele equations [64,67,69,76].
Miniaturization has another advantage, when dealing with useful active impurities scattered through a sample, in the possibility of reducing the averaging process over a large number of active "impurities". Finally, the possibility of dealing with a highly sensitive single impurity without any statistics can be reached. For instance, in colour centres, an impurity acts as a perfectly localized extended atom or molecule, where the spatial extension is basically driven by the host dielectric constant [14], and similarly the energy difference between electronic levels is reduced. Doping centres in semiconductors also act as localized atoms with their different electronic levels and ionization [15]. So, the purpose of miniaturization up to a single impurity sounds to be quite fruitful in various fields. It was reached for colour centres in the case of nitrogen vacancies in diamond, with very low doping and nanoparticles [16]. Then the local magnetic field can be detected by means of magnetic resonance on the Zeeman split levels [17]. With such a tool the resolution of a single spin near the surface can be reached [17,18]! Such studies open the path for an ultimate resolution of magnetic surface structures down to the atomic level. And quite numerous colour centres can be tried for such high-resolution applications.