We first consider the dynamics of a single vortex. As shown in Figure 1, a circular region with variable magnetic anisotropy (vortex region) is embedded in a slab with a fixed uniaxial OOP anisotropy (OOP region). When the uniaxial magnetic anisotropy constant of the vortex region (K _V) is initially set to be 0, the vortex region has an IP anisotropy, and the chiral vortex magnetic texture is stabilised by the interfacial DMI: H _DM = – D _ij · ( m _i × m _j ) where D _ij denotes the DMI vector of interfacial-type DMI, and the DMI energy is given by U D M I = D ∫ d V ( m z δ m x δ x − m x δ m z δ x + m z δ m y δ y − m y δ m z δ y ) . Unlike the symmetric exchange interaction that favours collinear alignment of magnetic moments, DMI prefers an orthogonal arrangement between the neighbouring spins [45], and here in the Pt/Co/AlO_x system it gives a fixed left-handed chirality determined by the negative sign of the Dzyaloshinskii vector D _ij. Hence, the chiral coupling generates skyrmion-like spin texture in the vortex region which induces a core pointing along the OOP direction opposite to that in the outside OOP regions [23, 50]. The equilibrium position of the core in the absence of external stimulus is located in the geometry centre of the vortex region for a relatively large DMI strength (| D _ij|≥0.75 mJ/m²) [50] (Figure 1A). When the DMI strength is too small, the equilibrium position of the core is at a finite distance away from the centre, whereas the vortex texture may turn into helical state for a large DMI strength [50, 51]. In order to stimulate the gyrotropic motion of the vortex core with STT, a spin-polarised current is perpendicularly injected in the vortex region. The self-sustained oscillation can only be obtained when the current is injected here in the +z direction, while the current injected in the -z direction breaks the balance between the forces (Figure 1B). Due to the rapid dissipation of the injected current and the much lower SOT current-to-spin conversion efficiency cf. the STT current-to-spin conversion efficiency [52, 53], the SOT spin current generated in Pt is negligible as compared to the spin current generated via STT and not considered in our structure. In order to quickly reach the steady precession state, a preset process is needed [54–56], which is done here by a applying a magnetic field of 5 mT parallel to the plane to preset the vortex core away from the geometry centre. The magnetic field is then removed. If there is no external stimulus, the vortex core relaxes back to the centre following a damped gyrotropic trajectory (black curve in Figure 2A and left inset of Figure 2B). In the presence of STT, the energy dissipation can be counter-balanced by injecting a continuous spin-polarised current in the plain normal direction (z direction) through an OOP magnetized top layer via a magnetic tunnelling junction in adjacent to the magnetic disk in experiments (Figure 1B). When the current density J is above a threshold, self-sustained oscillation of the core with a finite distance away from the geometry centre can be obtained (red curve in Figure 2A and middle inset in Figure 2B). As shown by the snapshots taken from different time for the oscillation at J = 5.5 × 10⁸ A/m² (insets of Figure 2A), all vortex cores keep the same distance away from the geometry centre. Figure 2B shows the spectra for oscillation frequency as a function of spin current J when K _V = 0. The oscillation frequency increases with the increase of J, while a larger current (above 1.35 × 10⁹ A/m²) destroys the self-sustained oscillation, leading to the collapse of the vortex structure (blue curve in Figure 2A and right inset of Figure 2B). When K _V increases to 4 × 10⁵ J/m³, the frequency of the whole spectra increases and the modulation of oscillation frequency occurs at a higher current (Figure 2C). The detailed mechanism for the modulation in oscillation frequency will be discussed later.

Note that the required STT current density (10⁸–10⁹ A/m²) is much lower than the spin-orbital-torque (SOT) current density (10¹¹ A/m²) in a similar structure [57] which is due to the following reasons: First, the charge-to-spin conversion efficiency in the STT (around 0.5 in Pt/Co system [52]) is usually higher as compared to SOT (around 0.1 in Pt [53]). Second, the spin current drives a gyrotropic motion of magnetic oscillator around the equilibrium spin texture, while the SOT switching needs to overcome the magnetic anisotropy and rotate the magnetisation for 180°. So the oscillation spin current density is lower than that for switching the magnetisation. Meanwhile, we have not considered the pinning effect which may increases the required spin current density in the real devices.

To study the influence of magnetic anisotropy on the dynamic behaviour of a single vortex, we fix the electric current J = 8.0 × 10⁸ A/m² and vary the magnetic anisotropy constant K _V from 0. When K _V is equal to 0, due to chiral coupling the vortex region forms a stable skyrmion-like spin texture with a small OOP core (left inset of Figure 2E) as a result of exchange singularity [23, 50]. The self-sustained oscillation, as shown in Figure 2D (curve in black for K _V = 0), has a constant amplitude with a fixed period of 32.6 ns, corresponding to a frequency of 306 MHz. With the increase of K _V, the oscillation frequency first increases slowly when K _V <3 × 10⁵ J/m³, and then increases drastically when K _V >3 × 10⁵ J/m³ (Figure 2E). Meanwhile, the size of the vortex core, as shown in the middle and right inset of Figure 2E (for K _V = 4 × 10⁵ J/m³), undergoes a significant expansion. The red oscillation curve in Figure 2D (representing the oscillation at K _V = 3 × 10⁵ J/m³), as compared to the black curve (K _V = 0) has a shorter oscillation period (25 ns) with a larger frequency (400 MHz). The oscillation amplitude, which is proportional to the average magnetisation of the vortex region, also decreases, due to the expansion of the core region and the reduction of the IP magnetised region. Interestingly, when K _V is beyond a critical point of ∼5.8 × 10⁵ J/m³, the spectra intensity gradually vanishes, indicating the breakdown of the self-sustained oscillation. As shown in the right inset of Figure 2E, the size of OOP core enlarges to be almost equal to the vortex region when K _V = 5.8 × 10⁵ J/m³, resulting in the annihilation of the vortex structure. A similar phenomenon is observed when J is increased to 13.5 × 10⁸ A/m² (Figure 2F), where the frequency increases less rapidly with the increasing of K _V, with a peak reaching at a similar K _V = 5.8 × 10⁵ J/m³ before the vortex structure collapses. The annihilation of the vortex structure due to the modification of magnetic anisotropy has been realized experimentally through VCMA effects [36, 37, 58, 59], where a well-defined IP anisotropy in the radial vortex region is crucial for maintaining the noncolinear spin texture of vortex structure through DMI. When K _V increases, the antisymmetric DMI is gradually taken over by the symmetric exchange interaction, while the latter favours a colinear and parallel spin alignment between adjacent regions, which eventually results in the collapse of the noncolinear vortex spin texture.

To further understand the tunable dynamic behaviour of the vortex oscillation under the influence of J and K _V, we refer to the centre-of-mass dynamics of the magnetic vortex oscillator, which has been well-modelled by the single collective variable Thiele equation [60–62]. G × d X d t + ∂ U ( X ) ∂ X + α D d X d t + F S T T = 0 (1)

The first term represents the gyrotropic force (also termed as Magnus force), where G = − 4 π Q μ 0 d M S γ e z is the gyrovector, with skyrmion number Q = ±1 for isolated skyrmion. M _S is the saturation magnetisation, d is the film thickness, γ is the gyromagnetic ratio, μ 0 is the vacuum permeability and e z is the unit vector. The rigid profile of the vortex core is described by its position X = X (x(t), y(t)) in the plane, and its velocity is given by d X d t . The gyrotropic force is radial outward which primarily acts to expel the core to the outer shell of the vortex [62]. The second term represents the repulsive (or restoring) force which originates from the edge effect: when the core moves close to the boundary of the vortex structure, the magnetic charges distributed at the cylindrical side wall will generate a centripetal confine potential U ( X ) to regulate the core inside the vortex region [63]. As can be seen from the form, this second term ∂ U ( X ) ∂ X represents a conservative force acting on the vortex core that is derived from the total free energy of the system [13]. The edge-induced repulsive force together with the remanent gyrotropic force provide the centripetal force to sustain the circular motion of the core [54], as illustrated in Figure 1B. The last two terms in Eq. (1) represent the damping force and spin transfer forces, respectively [54]. The spin polarised current perpendicular to the layers results in a torque force on the vortex core, which can be balanced by the dissipative force, and, in overall, they jointly support orbital motion of the vortex core [23]. It is worthwhile to mention that Eq. (1) describes an intrinsically non-Newtonian dynamics, which means that the response of the vortex core to a conservative force is gyrotropic, i.e. in a direction perpendicular to the force, instead of accelerating in a direction collinear with the force. This gyrotropic dynamics therefore underpins the persistent oscillations of the vortex core under the four different contributions [13], as depicted in Figure 1B.

We first consider the role of spin polarised current J. As discussed above, the confining potential provides an inward radial repulsive force (second term in Eq. 1) to confine the vortex inside the vortex region, while the spin transfer torque (STT force, the last term in Eq. (1)) arising from J tends to drive the core outward toward the edge, so does the gyrotropic force (first term in Eq. (1)). With the increase in J, the radius of the steady state orbit increases (as shown in the middle inset of Figure 2B in comparison with the left inset), as a larger dissipation rate is required to compensate the increase in spin torques [50]. The gyrotropc force increases as well to counterbalance the increase in the STT force, resulting in the increase in the velocity of the core d X d t [64]. This entails higher oscillation frequencies above the intrinsic resonance frequency (Figures 2B,C). On the other hand, the radial inward repulsive force has an upper limit because of the limited confine potential U ( X ) [62]. The increase in J therefore leads to a growing outward net force to push the core towards the boundary. When J reaches a critical value, the core cannot be confined inside the vortex region and eventually be pushed out (as shown in the right inset of Figure 2B), leading to the collapse of the vortex spin texture [62]. To more effectively confine the core and maintain the vortex, the upper critical current density should be increased. This can be done by increasing the edge potential using nanoring with high perpendicular magnetic anisotropy [62]. As discussed below, the increase in the K _V in the vortex region can serve as an alternative stabiliser to increase the confine potential.

We now turn to understand the role of K _V on the oscillation dynamics. The intrinsic oscillation frequency is determined by the curvature of the confining potential at the equilibrium core position in the absence of STT force [50]. We first calculate the relative potential energy of system (U-U ₀, where U ₀ is the potential energy when the core is in the centre) with respect to different core distances from the centre (| X |) at different K _V values in the absence of STT (Figure 2G). Since the equilibrium position of the core is in the centre, U-U ₀ increases when the core displaces away from the centre. Interestingly, when K _V increases, U-U ₀ increases more rapidly with displacement, especially at large K _V values, indicating that the potential energy becomes higher when the core has the same displacement. This, on one hand, directly increases the slope of the curve as observed in Figure 2G, corresponding to an enhanced ∂ U ( X ) ∂ X , and causes the increase of the repulsive force (second term in Eq. 1). The increase in the repulsive force can in turn be counterbalanced by increasing the gyrotropic force (first term in Eq. 1) when J is fixed, as long as self-sustained oscillation is maintained. The increase in gyrotropic force eventually leads to the increase in the velocity d X d t and hence the oscillation frequency of the vortex core increases (Figures 2E,F). On the other hand, the increase in confine potential caused by the increase in K _V requires a larger STT current to maintain the oscillation, which elevates the upper limit of the critical current density (J _C) for the self-sustained oscillation before the annihilation occurs. Figure 2H shows the J _C under a 150 ns running interval: when K _V increases from 0 to 5 × 10⁵ J/m³, J _C increases monotonically from 1.35 × 10⁹ A/m² to 4.25 × 10⁹ A/m².

The potential energy is closely correlated to the evolution of the vortex structure. As mentioned above, the vortex core expands rapidly when K _V increases (insets in Figure 2E), which leads to the shrinkage of the noncolinear spin region. The enlarged OOP core and the reduced vortex area make it more difficult for the core to move, and thus a higher exchange energy is needed, which results in the increase in the overall potential energy and the oscillation frequency. Hence, the tuning of the oscillation frequency by K _V is understood.

We next turn to investigate how the spin current and magnetic anisotropy affect the mutual interactions for a couple of oscillators. In the simulation, two vortex oscillators are placed with a given distance, while the J and K _V are varied (Figure 3A). As mentioned earlier, interacting non-identical oscillators with different intrinsic frequencies tend to undergo resonance at a common frequency when their oscillation frequencies are within an certain locking range, resulting in the synchronization phenomena [23, 50]. Figure 3B shows the synchronization spectra as a function of the current density for a pair of oscillators (named as one and 2) with a 100 nm distance when K _V = 0. When the spin polarised current of oscillator 1 (J ₁) is fixed at 8 × 10⁸ A/m² and that of oscillator 2 (J ₂) varies from 0 to 1.25 × 10⁹ A/m², two distinct frequency peaks are observed at small J ₂, and then they merge into one peak with an enhanced intensity, undergoing synchronization when J ₂ ranges between 6 × 10⁸ A/m² and 1 × 10⁹ A/m². The merged peak diverge again when J ₂ > 1 × 10⁹ A/m². In order to study the influence of magnetic anisotropy on the synchronization, the magnetic anisotropy constant in the vortex region of oscillator 1 (K _V1) is fixed at 1 × 10⁵ J/m³ while that of oscillator (K _V2) varies from 0 to 3 × 10⁵ J/m³. We first study the case where different currents are injected into the two oscillators: J ₁ (8 × 10⁸ A/m²) ≠ J ₂ (6.5 × 10⁸ A/m²). As shown in Figure 3C, the two oscillators synchronize with each other when K _V2 ranges from 6 × 10⁴ J/m³ to 1.5 × 10⁵ J/m³. The synchronized region ΔK _V2, defined as the length of the interval between the starting and ending point of K _V2 of the synchronization region, is estimated to be 0.9 × 10⁵ J/m³. It is noteworthy that the frequency spectra and synchronization region keep unchanged when the two different oscillators swap their relative positions, as shown in Figure 3D, indicating that the contributions from two oscillators are switchable with regard to position.

The synchronization of nano-oscillators in this system can be realized through propagating spin waves or dipolar interactions [14], where the distance between the oscillators plays a crucial role to determine the synchronization stability and coupling strength [50]. Here, to better understand the distance effect, the separation between the oscillators is varied from 60 to 150 nm. Figures 4A–C show the synchronization spectra with a distance of 60, 80 and 100 nm, respectively, for J ₁ (8 × 10⁸ A/m²) ≠ J ₂ (6.5 × 10⁸ A/m²) (termed as ‘nonequal case’) and Figures 4D–F show the respective spectra for J ₁ (8 × 10⁸ A/m²) = J ₂ (8 × 10⁸ A/m²) (termed as ‘equal case’). In both cases, the synchronization range ΔK _V2 decreases with the increase in separation distance (Figure 4G), which could be attributed to the decay in coupling strength [50]. Interestingly, the decrease in ΔK _V2 for the nonequal case is much more rapid than the equal case: at a small distance of 60 nm, ΔK _V2 of the nonequal case is similar to that of the equal case, both having a large synchronization interval of ∼3 × 10⁵ J/m³; in contrast, when the distance increases to 150 nm, ΔK _V2 for the equal case keeps at a sizable value 1.7 × 10⁵ J/m³, while for the nonequal case it decays to 0.1 × 10⁵ J/m³. This stark difference is a clear indication of the faster decay for the overall strength of oscillator interaction with the increase in separation distance when the two current densities are different, which could be due to the relative contribution from different source of interaction (dipolar interaction or spin wave propagation) changes when the current density and distance are different.

The synchronization of nano-oscillators is of great interest for neuromorphic computing. As a proof-of-concept demonstration, we design a neural network comprised of six vortex oscillators (A, B and 1, 2, 3, 4) (Figure 5A), and explore the ability of coupled vortex nano-oscillators to work as a model system for neuromorphic computing via anisotropy control. Figure 5B shows the schematic for the intercorrelation between oscillators, neurons, inputs and outputs. The oscillators are separated from each other with a certain distance, with 1–4 (output neurons) surrounded by A and B (input oscillators), so that all oscillators can be mutually interacted with each other. As shown below, the change in K _V of the two input oscillators can modify the synchronization configurations of the whole neural network, mimicking the conversion of encoded inputs to classified output signals via programmable weighting mechanism [30]. Hence, the synchronization configuration between the six oscillators is interpreted here as the output signal.

To begin with, we set up one input signal using the magnetic anisotropy of either oscillator A or B (K _VA or K _VB), and then increase to two input signals by adopting both K _VA and K _VB. The current density J for all the oscillators is set to be 8 × 10⁸ A/m², and the K _V for the output neurons 1, 2, 3, and four are set to be 1.1 × 10⁵ J/m³, 0.8 × 10⁵ J/m³, 1.2 × 10⁵ J/m³, 0.9 × 10⁵ J/m³, respectively. When K _VA and K _VB are varied from 0 to 4 × 10⁵ J/m³, the frequency patterns of the neural network are converged and allocated to different synchronization configurations appearing at different K _VA and K _VB values. For example, when K _VA = 2.4 × 10⁵ J/m³ and K _VB = 1.6 × 10⁵ J/m³, five oscillators (1, 2, 3, 4 and B) are synchronized into one common frequency (Figure 5C) during a running period of 50 ns, corresponding to an output synchronization code of (1234B). A total of 13 output patterns are obtained, as represented by 13 different colours organised into a two-dimensional colour map (Figure 5D), while only up to six configurations can be recognised with one input (either K _VA or K _VB). Interestingly, when J _B changes to 9 × 10⁸ A/m², the synchronization pattern is modified (Figure 5E): the synchronization codes (124B), (134, 2AB) and (124B, 3A) that appear in the case of J _B = 8 × 10⁸ A/m² now disappear, while another three codes (13A, 24), (134, 2A) and (124) emerge. This indicates that current can work jointly with the magnetic anisotropy to modulate the output synchronization patterns.

In the proof-of-concept neural network shown in Figure 5, we have six oscillators, each having independent current density and magnetic anisotropy. The neural network hence has twelve variables. Two variables (K _VA and K _VB) can act as inputs and tune the synchronization patterns as the outputs (Figure 5D). We also found that the synchronization pattern changes when the current density of one oscillator is changed (Figures 5E vs 5D), indicating that the output depends on the combination of all variables. When two variables (K _VA, K _VB) are defined as the inputs, the combination of the rest ten variables (K _V1 … K _V4, J ₁, … , J ₄ and J _A, J _B) can be defined as the weight and bias. When these ten parameters are changed, the weight and bias are changed, resulting in a different output pattern (as shown in Figures 5D,E). Therefore, we can train the neural network to get a target output pattern by adjusting the weight and bias parameters via a learning rule.