Texture fluctuations and emergent dynamics in coupled nanomagnets

We analyse the thermal fluctuations of magnetization textures in two stray field coupled elements, forming mesospins. To this end, the energy landscape associated with the thermal dynamics of the textures is mapped out and asymmetric energy barriers are identified. These barriers are modified by changing the gap that separates the mesospins. Moreover, the coupling between the edges leads to an anisotropy in the curvature of the energy surface at the metastable minima. This yields a dynamic mode splitting of the edge modes and affects the attempt switching frequencies. Thus, we elucidate the mechanism with which the magnons in the thermal bath generate the stochastic fluctuations of the magnetization at the edges.

Single domain magnetic nanoislands with binary magnetization states -mesospins -are used as building blocks for magnetic metamaterials.Contemporary examples are e.g.artificial spin ices (ASI) [1][2][3][4][5], Ising chains and lattices [6,7], and reconfigurable magnonic crystals [8][9][10].These magnetic metamaterials can be viewed as having thermal fluctuations, associated with switching of the magnetic states of the mesospins [11][12][13][14][15], in an activated process [12,13,16,17].However, the magnetic fluctuations are not solely restricted to the switching of rigid mesospins, but can also occur in the interior textures of the mesospin, which necessitates a more nuanced view of their thermal excitations.Furthermore, a finite temperature leads to the excitation of magnons, which take place on timescales of the magnetization dynamics, much smaller than those of the fluctuations of individual mesospins [18].
Resonances can occur in mesospins when the wavelength of magnons matches their extension.At frequencies below the uniform (Kittel) mode, there is a band gap, characterized by the absence of magnon modes in the interior of the mesospins.However, dynamic edge modes may occur at frequencies inside to the band gap, which can be used for spectral detection of topological defects in ASI [19].We have previously confirmed the existence of additional modes [17,20], associated with a binary stochastic fluctuation of the magnetization at the edges, from hereon called edge fluctuations.These edge fluctuations cause switching between S-or C-states [21], and their presence results in a residual entropy in magnetic metamaterials [17,22].They have recently also been confirmed experimentally, and were found to occur at shorter time-and length-scales than conventional mesospin switching [17].Thus so far, the intra-nanomagnet excitations are well understood, however an understanding of the coupling between intra-and inter-mesospin degrees of freedom has not been developed for magnetic metamaterials.
For Ising-like mesospins, the edge fluctuations of the texture occur in the transverse direction with respect to the switching of the net magnetization (e.g. in the example of Fig. 1a, m y is the relevant reaction coordinate in phase space).The height of the barrier associated with these fluctuations depends on the size of the nanomagnets, and the corresponding balance between exchange and demagnetization energies.In the case of a single mesospin, the relevant phase space for the edge fluctuations can be captured in a 1-dimensional energy landscape.When adding one more mesospin, the texture excitations are influences by inter-mesospin interactions, with the phase space extending over a 2-dimensional energy landscape.The two nanomagnet system is in the ground state when the two mesospins are magnetized in the same direction, in which case the flux exiting one mesospin will be absorbed by the other, which results in a reduction of the transverse magnetization at the edges.In terms of the energy landscape, this situation leads to a single, global energy minimum, which is associated with to the two edges being aligned.In this case, the coupling causes the splitting of dynamic edge modes into one in-phase and one out-of-phase mode.On the other hand, when the two mesospins are oppositely magnetized (as illustrated in Fig. 1a), there will be two degenerate metastable states for the edge magnetization, leading to non-trivial behaviour upon thermal excitation.
In this work, we discuss the fluctuations of two oppositely magnetized mesospins.To this end, we start by exploring and quantifying the energy landscape associated with the intermesospin excitations.We then perform micromagnetic simulations with a temperature implementation, in order to inspect the behaviour of the system under thermal excitation, for a range of different temperatures and spacings between the mesospins.The observed switching events are used to infer the energy barriers and attempt frequencies utilizing Arrhenius analysis [23].Finally, the Arrhenius prefactors are discussed, relating them to the dynamic modes of the system and the curvatures in the energy landscape.
For the numerical investigation, we performed micromagnetic simulations in MUMAX 3 [24], which is based on the use of a finite difference method to solve the Landau-Lifshitz-Gilbert (LLG) equation of motion.The size of the stadium shaped mesospin used in this work is l × w × h = 360 × 120 × 4 nm 3 , using rectangular cells of size 2.5 × 2.5 × 4 nm 3 .We use parameters that are relevant for permalloy, i.e. saturation magnetization M s = 1 • 10 6 A/m, an exchange stiffness of A ex = 1 • 10 −12 J/m, and a Gilbert damping of α = 0.001.More information on the details of the simulations can be found in Slöetjes et al. [20].
Two mesospins are placed in line, separated by a gap, g, that defines the edge-to-edge distance, as illustrated in Fig. 1a.We define the magnetization in the edges as the averaged transverse magnetization over the half circle that constitutes an edge, i.e. m y,1 = 1 V edge 1 m y (r,t)d 3 r for the edge of mesospin 1, and similarly for mesospin 2, where V is the volume of the mesospin edge.The lowest energy configuration, state A, has opposing transverse magnetization components, m y,1 = −m y,2 , and state B has aligned transverse components, m y,1 = m y,2 , see Fig. 1a.
We represent these states in an energy landscape, parameterized by m y,1 and m y,2 , see Fig. 1c.This is an approximate representation of the energy landscape, since the Zeeman energy, E Z , is included in the total energy, therefore it only serves to indicate the morphology.Aside from the states A and B (and the degenerate states A and B ) one can observe the saddle points at S and the top of the energy barrier at C. As the gap g between the mesospins increases, the points A and B shift towards the same energy, while the points S and C merge having the same height.Bringing the mesospins closer together causes the energy hill to emerge at point C, pushing the states A and B outwards, and breaks the two-fold symmetry.Thus, the gap between the mesospins provides a handle to tailor the energy landscape.
When switching from State A to the opposite State A , it is not energetically favorable to take the direct route A −C − A , corresponding to both moments switching simultaneously.Instead, it is more favorable to take the route A − B − A via the saddlepoints, which corresponds to consecutive switchings of the nanomagnet edges.Due to the symmetry of the energy landscape, it is only necessary to know the minimum energy path A − S − B for an understanding of the fluctuations.We map out the energy barrier along the path that is highlighted in Fig. 1c (see Appendix A for the method), the resulting morphologies are shown in Fig. 1b.As previously noted, the energy landscape becomes increasingly asymmetric as the gap between the mesospins is reduced.The states A and B are simultaneously pushed away from the center.Furthermore, the saddle point moves away from the center position, towards state B. From hereon, we refer to the energy difference between A and S as the high energy barrier (E H ), and the energy difference between B and S as the low energy barrier (E L ).
Next, we will discuss the thermally induced dynamics of the two oppositely magnetized mesospins.To this end, we use micromagnetic simulations at a finite temperature, as described in Leliaert et al. [25].The thermal bath is imitated by transforming the LLG equation to a Langevin equation using a stochastic thermal field, given by: where k B , T , γ, V , and ∆t are the Boltzmann constant, temperature, gyromagnetic ratio, volume of the cell, and the time step, respectively.In addition, η is a random vector which changes direction and size at every timestep.We run the simulations for temperatures from T = 25 to 450 K in steps of 50 K, and eight different gaps, ranging from g = 25 to 200 nm.We analyze the timetraces of the perpendicular component of the magnetization in the two edges m y, j (t), where j denotes the particular edge.Examples of these time traces are shown in Fig. 2a, for a large (g = 200 nm) and a small (g = 25 nm) spacing.Each edge can be seen to switch between the two different states, and the absolute value of the perpendicular magnetization component increases as the spacing becomes smaller.The thermal energy k B T required to switch the edge magnetization needs to be increased by more than an order of magnitude for g = 25 nm compared to g = 200 nm, in order to obtain a similar switching rate.Furthermore, the magnetization at the edges is weakly anti-correlated when the mesospins are placed far apart, and strongly anti-correlated when they are close together.This correlation can be quantified using the Pearson correlation coefficient, ρ(m y,1 , m y,2 ), which is shown in the Appendix C for different spacings and temperatures.
When the two mesospins are strongly coupled, the edges seemingly switch simultaneously, but zooming in on the switching process (upper panel of Fig. 2b) reveals that the two edges switch in a consecutive manner, as was predicted on the basis of the energy landscape.Another possibility for a fluctuation is that the edge magnetization switches back to its initial state, as seen in the lower panel of Fig. 2b.We have observed (see Appendix D) that the process A − B − A occurs more often that the process A − B − A, while in principle these two processes are the same in terms of energy barrier.This asymmetry becomes more pronounced as the gap between the nanomagnets is reduced.We will now set out to determine the thermal energetics of the textures in the mesospins.The average time in each state is determined, and we use Arrhenius law to extract the energy barriers and attempt frequency: where i indicates A or B for τ, and H and L for E respectively, with τ 0,i being the inverse attempt frequency.Generally, for inverse gaps larger than 1/g = 0.027 nm −1 , the errors for E L diverge, due to the high temperatures used for the relatively low barriers, resulting in erroneous estimations of the switching time.For the largest gap, 200 nm, there is no notable difference between E L and E H , and the energy barrier (E = 8 meV) is close to that of a single mesospin [20] (E = 5 meV).Bringing the mesospins closer together leads to an expected increase in E H and E L .Whereas E L scales with the inverse gap, E H shows two different scaling regimes, as shown in Fig. 2d.It can be seen that up to an inverse gap of 1/g < 0.02 nm −1 , E H shows a parabolic scaling with 1/g, whereas for values 1/g ≥ 0.02, there is a linear scaling.The transition between the two regimes occurs between 1/g = 1/65 nm −1 and 1/g = 1/50 nm −1 , which means that the radius of the half circle at the edge (60 nm) defines the scale of the transition length.Overall, good agreement can be found between the energies obtained via the stochastic simu-lations and the energy landscape method, confirming the Arrhenius law.
Next, we shall focus on the attempt frequency, f 0 , calculated as the inverse of τ 0,i .Similar to the energy barriers, the attempt frequencies for State A and B also show a diverging behaviour upon a decrease of the gap (see Fig. 3a).At the largest spacing, the attempt frequency at A and B become equal, although they are not equal to that of a single mesospin [20].Interestingly, it can be observed that f 0,A increases when the nanomagnets are being brought closer together, while f 0,B shows a slight decrease.In order to understand the scaling of the attempt frequencies, we shall inspect the mode frequencies at the edge of the mesospins, as well as the curvature in the energy landscape at A and B.
The spectrum of magnon modes at the A state is different from those at B, see Fig. 3b.The edge modes reside in the frequency band below the Kittel mode (∼ 6 GHz) due to the demagnetizing field counteracting the effective field at the edges, thus lowering the eigenfrequencies with respect to the modes inside the mesospin [20].We can observe, for both cases, that a mode splitting occurs as the mesospins are brought closer together.The mode splitting in B happens at a larger spacing than the splitting in A, which can be understood as being due to the difference in strength of the dipole coupling: in B, the sources of the stray field are closer together than in A. In both cases, the two resulting branches after the splitting have opposite phases in the in-and out-of-plane directions of the precession.By using phase maps of the mesospins with the full inner texture taken into account (see Appendix E), we determine the in-plane (∆ϕ ) and out-of-plane (∆ϕ ⊥ ) phase differences of the modes in the two edges.The out-of-plane components are very small compared to the in-plane ones, i.e. the precessional motion is highly elliptical, therefore we will only consider the in-plane component, where in-and out-ofphase oscillations are indicated by ++ and +−, respectively.From Fig. 3b, it can be seen that the phases of the high and low frequency branches in B are inverted with respect to the branches in A. Moreover, at State A, both modes increase in frequency as the gap is decreased, whereas in B, the −+ mode decreases in frequency, which suggests that the −+ mode is responsible for switching from B to A.
We turn to the energy landscape to investigate this scenario, where we map the trajectories of the modes at A and B, shown in Fig. 3c, for a gap of g = 25 nm.The −+ and ++ modes represent straight paths in the diagonal direction in the (m y,1 , m y,2 ) phase space.The mode splitting occurs because the curvatures of the (m y,1 , m y,2 )-energy landscape at the metastable minima become anisotropic, which is a result of the dipolar coupling between the edges.We plot the minimum energy path from A to B in the vicinity of the minima, and inspect the tangent with the directions of the dynamic modes.At point A, the path is mostly tangential to the ++ mode.In contrast, the path emerging from point B is mostly tangential to the −+ mode, which indeed shows that it is the −+ mode is the dominant contribution for switching from B to A. These findings suggests that the transitions from A to be B can be stimulated using an applied microwave field with a frequency in the ++ branch, and vice versa with a frequency in the −+ branch, although the latter will be problematic in a ferromagnetic resonance experiment, because of the symmetry of the mode.
Since the damping is low, we can assume that f 0,i ≈ f i , i.e. the attempt frequency is approximately equal to the eigenmode in the direction of the minimum energy path [26].The frequency is related to the curvatures in the energy landscape via f ∝ (K) [27], where K is the total curvature, given in terms of the principle curvatures by K = κ κ ⊥ , here κ j = d 2 E/dm 2 j , where j runs over the in-and out-of-plane components (see Appendix B for more details).Thus, upon decreasing g, the curvature in the −+ direction at point B should decrease, while at point A, both curvatures in the −+ and ++ directions increase.We found that the curvature associated with bringing the magnetization out of the plane, κ ⊥ , is more or less constant as a function of g, which is due to the strong out-of-plane demagnetizing field overwhelming any effects of proximity to the other edge.The in-plane component, κ varies strongly, as seen in Fig. 3d.At A, κ shows an increase with 1/g, in qualitative agreement with f 0,A , whereas the curvature at B remains constant, which is not in line with the decrease of f 0,B .A possible explanation for this discrepancy could be attributed to the path B − A being favored over B − A as the gap is reduced (see Appendix D), thus suppressing one of the channels for switching.The suppression of this channel effectively lowers the entropic contribution, by reducing the multiplicity.It has been seen before in works by Desplat and Kim [27], Desplat et al. [28], that entropic contributions can have large effects on the attempt frequency, which occur when the positive curvatures at the saddle point are small, allowing for an increased number of paths from A to B. The positive curvature at the saddle point increases with the height of the energy barrier, thus restricting the saddle point and lowering the amount of paths.However, this effect is compensated by the increase in the curvature at the minima.
In conclusion, we have mapped out the energetics associated with the coupling of intra-and inter-mesospin degrees of freedom in a system of two interacting mesospins.It was found that, upon thermal excitation, a stochastic switching occurs between two degenerate states, A and B, and that the reciprocal gap between the nanomagnets provides a handle on the symmetry and barrier height in the energy landscape, modifying in turn the correlations between the edges.These correlations feature a strong temperature dependency, that could potentially be harvested in nanomagnetic neuromorphic-like systems, having the opportunity to operate at sub-ps timescales controlled by light, as demonstrated by Pancaldi et al. [29], and Mishra et al. [30].Moreover, by associating the dynamic modes of the system with directions in the energy landscape, we were able to explain the trends in the attempt frequency upon a variation of the gap value.As such, we shed light on how the thermal bath, i.e. the magnonic excitations, couples to the stochastic fluctuations on the next length scale.The stochastic fluctuations of the edges may in turn affect fluctuations on the next length scale, i.e. the full switching of the mesospin.For example, a certain state A can induce a preferred helicity for switching the magnetization in the m x direction, similar to the symmetry breaking pointed out by Leo et al. [31].The double degeneracy of the metastable states of the edge magnetization can in this case cause an increase of a factor two in the switching rate of the mesospin island.Tunable magnetization texture states have also been recently shown to be useful building block of arrays, where their magnonic properties can be harnessed for hardware implementations of reservoir computing schemes [32,33].
The energy barrier morphology is mapped out by "pushing" the system to a point close to the saddle point using an external field, and then letting it relax according to the Landau Lifshitz Gilbert (LLG) equation, which is given by: where γ, α, m, and H are the gyromagnetic ratio, Gilbert damping constant, magnetization vector, and the effective field, respectively.The effective field acts as a proxy for including the different energy contributions in the equation of motion via where µ 0 is the permeability of free space, and M s is the saturation magnetization.The total energy is a sum of the separate energies: E = E demag + E exchange + E Zeeman .The first term on the right hand side of the LLG equation describes the precession of the moment around the effective field, and the second term describes the movement of the magnetization towards the effective field.For the LLG equation to correctly predict the minimum energy path, the precession needs to be suppressed, therefore we set the damping to the high value of α = 10.
The system is first configured in either the A or B state (see Fig. 4a).Then, the field is applied to region 1 to push one of the edges to a point close to the saddle point S, after which the field is taken away, and the system evolves for 60 ns to either local minimum, A or B. During this process, the energy is recorded every 1 × 10 −12 s, giving a highly resolved trajectory in phase space.The resulting curves can be seen in Fig. 4b,c.The minimum energy path indicates that the majority of the change occurs in m y,1 , with a small response of m y,2 .
A cusp is visible near the saddle point, which is a result of the system having relaxed in an external field, thus not being precisely positioned on the minimum energy path.Thus, the cusp signifies the system moving from the position in the field to the minimum energy path.This cusp is relatively small and can be removed, after which we fit a 10 th order polynomial to the curve projected on the (m y,1 , E)-plane.This fitting procedure mostly serves to bridge the gap around the saddle point (see Fig. 4d).Strictly, the fitting should take place in three dimensions (m y,1 m y,2 , and E), but since m y,2 varies very little, we are only interested in the height of the energy barrier and later on the curvature at S, fitting the curve only to the projection is justified.With this method we recover the minimum energy path from A to B, yielding the energy barriers, E L and E H .This method is relatively simple and can be entirely carried out using MUMAX3 in the LLG framework.The only drawback is that it demands prior knowledge of the minimum energy path through phase space, in contrast to more sophisticated methods such as the Nudged Elastic Band (NUB) method, and variations thereof [34].

APPENDIX B: CURVATURE CALCULATIONS
As stated in the main text, we evaluate the curvatures by taking the second derivative of the energy to the magnetization, κ = ∂ 2 E/∂ m 2 , at certain points in the energy landscape.

Curvature at S
From the minimum energy path, we obtain the in-plane curvature at the saddle point, κ S , for which ∂ E/∂ m y,1 = 0 at the top of the energy barrier.We evaluate the curvature only in the (m y,1 , E)-plane, as the change in the m y,2 -component is negligible at S.

Curvatures at A and B
Next we calculate the curvature at A and B via a different method, as the morphology that we previously mapped out only provides one half of the metastable minimum.Instead, we use a "guiding" field parallel to the edge magnetization to wiggle the edge around A and B. The guiding field method is similar to the one used for mapping out the energy landscape.However, to obtain the precise value of the curvature, we should eliminate the Zeeman energy, E Z , contribution.In order to eliminate this contribution, we extrapolate to the zero field value of the curvature after calculating the curvature for different sizes of guiding fields, from 15 mT to 1 mT in steps of 1 mT (see Fig. 5c).Mathematically, this can be expressed as: The extrapolation method is based on a spline method, performed in MATLAB.An example of such an extrapolation for a gap value of g = 25 nm, is shown in Fig. 5d.
Here, it is necessary to take into account both perpendicular components, m y,1 and m y,2 , since the angle of the path with the m y,1 axis can be quite substantial at A and B, and varies as a function of 1/g, especially at State B (see Fig. 5d).We can determine the curvature of the parabola projected on the m y,1axis, given by ∂ 2 E/∂ m 2 y,1 , which in turn provides the total curvature of the parabola through: where θ is the angle between the m y,1 and m y,2 components of the path.The square root of the in-plane curvatures for all the simulated gap values are shown in Fig. 5a.

Out of plane curvatures at A and B
The out-of-plane curvatures were determined by simply relaxing the magnetization to either A or B, by slightly pushing the magnetization out of the plane using an external field and measuring the energy versus the magnetization response, m z,1 .
Here, the contribution from the Zeeman interaction can be neglected, since the field is perpendicular to the magnetization, and E Zeeman = −µ 0 m • H.A parabola can the be fitted to find the out-of plane curvature, κ A,B ⊥ .The result is shown in Fig. 5b.
The results of this calculation for some values of the gap are shown in Fig. 6.It can be seen that the correlations decrease with increasing temperature, as well as an increasing gap.We fitted these data points using the following function [35]: This implies that the correlations are only dependent on the difference between the energy barriers, as well as the temperature.

APPENDIX D: RATIO OF A-B-A TO A-B-A' SWITCHING EVENTS
Figure 7 shows the ratio of switching processes in which the system starts at state A, propagates to state B, and then goes to either state A or A'.It can be seen that this ratio scales mostly with the inverse gap, the temperature not playing a large role.
T m y,1 m y,2 dt, and considering just two intervals, τ A , and τ B , in which case T = τ A + τ B .Now, τ A runs from t 0 to t 1 , and τ B from t 1 to t 2 , thus the integral can be split into the two contributions.Assuming that m y,1 m y,2 = −1 in State A, and m y,1 m y,2 = 1 in State B (the normalized case), we obtain ρ(τ A , τ B ) = τ B −τ A τ A +τ B , which can be seen to converge to -1 as τ B → 0, and 1 if τ A → 0. Substituting the Arrhenius law, and assuming equal attempt frequencies for A and B gives , which is equal to Eq. 8. [

FIG. 1 .
FIG. 1.(a) Definition of the the geometrical parameters and the edge magnetization states (A and B).(b) The energy barriers associated with switching the edge magnetization of one of the edges, with indication of the points A, S and B, for different gap values g.(c) Morphology of the energy landscape as a function of the perpendicular edge components for a g = 30 nm gap between the elements.

FIG. 2 .
FIG. 2. (a) Timetraces of the edge magnetization for two different gaps and temperatures.The shaded regions mark the system being in State B, the regions without shading mark state A. (b) Switching events for a gap of 25 nm and a temperature of 400 K. Upper panel: a true switching event, lower panel: a false switching event.(c) Arrhenius plots for the relaxation times associated with the high barrier, τ A .(d) The energy barriers E H and E L against the inverse gap values obtained via Arrhenius law, and direct mapping of the energy barrier.The dashed lines are linear fits in the case of E L and the low gap regime of E H , and a parabolic fit to the high gap regime of E H .

FIG. 3 .
FIG. 3. (a) Attempt frequency dependence on the gap for the states A and B. (b) Dynamic mode spectra for the states A and B, generated by a small thermal field.(c) The minimum energy path emerging from A and B (blue) with the corresponding direction of the in-phase and out-of-phase modes (black).The contours in the background are the same as in Fig. 1c (d) Curvatures in the direction of the minimum energy path.

FIG. 4 .
FIG. 4. (a) The two mesospins, prepared in state B. The field is applied to region 1, which is defined by the half circle that constitutes the edge on mesospin 1.This field pushes the state towards the saddle point S. (b) The paths from a point in the vicinity of S to A and to B, mapping out the majority of the minimum energy path from A to B. (c) The projections of the paths shown in (b) on the (m y,1 , E), (m y,2 , E), and the (m y,1 , m y,2 )-planes.(d) The simulated path with the cusps removed, after fitting to a 10 th order polynomial.

FIG. 5 .
FIG. 5. (a) In-plane curvatures of the minimum energy path at A, B, and S. (b) Out of plane curvatures at A and B. (c) The "wiggling" of the magnetization at A, for different applied field strengths, and (d) the curvatures for different guiding fields.The dashed line denotes the extrapolation to 0.

FIG. 6 .
FIG.6.Correlation values: The symbols are calculated from the time traces for the spacings indicated in the legend and in nm.The dashed lines are fits using equation 8.

FIG. 7 .
FIG. 7. The asymmetry between A-B-A' ('true switches') versusA-B-A ('back switches') against the inverse gap (upper panel) and temperature (lower panel).The two outliers at are likely the result of low statistics.