Walker-like Domain Wall breakdown in layered Antiferromagnets driven by staggered spin-orbit fields

Within linear continuum theory, no magnetic texture can propagate faster than the maximum group velocity of its spin waves. Here we report a transient regime due to the appearance of additional antiferromagnetic textures that breaks the Lorentz translational invariance of the magnetic system by atomistic spin dynamics simulations. This dynamical regime is akin to domain wall Walker-breakdown in ferromagnets and involves the nucleation of an antiferromagnetic domain wall pair. Subsequently, one of the nucleated 180$^{\circ}$ domain wall creates with the original domain wall a 360$^{\circ}$ spin-rotation which remains static even under the action of the spin-orbit field. The other 180$^{\circ}$ domain wall becomes accelerated to super-magnonic speeds. Under large spin-orbit fields, multiple domain wall generation and recombination is obtained which may explain the recently experimentally observed current pulse induce shattering of large domain structures into small fragmented domains and the subsequent slow recreation of large-scale domain formation prior current pulse.

Within linear continuum theory, no magnetic texture can propagate faster than the maximum group velocity of its spin waves. Here we report a transient regime due to the appearance of addi- In a general physical context, a topological defect is characterised by a local/core region where order is highly altered whereas far from that core region the order parameter varies smoothly in the space 1 .
Topological defects can be found in a variety of fields such as in superfluid helium 2 named vortices, in periodic crystal structures, they are called dislocations [3][4][5] and in magnetism, magnetic domain walls 6,7 , vortices 8,9 and skyrmions 10 . In elementary dislocations, the underlying reason why these defects are called topological is because they cannot be made to disappear by a continuous deformation of the order parameter in infinite extended films. The measure of the distortion is associated to the so-called Burgers vector whose magnitude expresses the strength of the dislocation. The rule consists in taking a closed path within a chosen portion of the crystal lattice, moving in steps of the lattice parameter and keeping the number of steps in each direction equal as the circuit is completed. If the section of crystal is defect-free, the end point will then coincide with the starting point. If however there is a defect contained within the circuit, then the end-point will not be equal to the starting point. The vector from the end-point to the starting-point is then the Burgers vector and its magnitude is used as a measure of the strength of the dislocation. The Burgers vector is quantised to integer multiples of the atomic spacing of the lattice vectors 11 . Continuum-based textures such as magnetic textures like domain walls, vortices and skyrmions are often characterised by the so called winding number whose value stands in direct relation to the magnitude of the normalised Burgers vector. However, topology is not at odds with instability.
A relevant example of the previous is the plastic deformation and its role in dislocation dynamics 12 .
Under low applied stresses, there is a quasilinear dependence between the dislocation velocity and the applied stress described by the linear elasticity theory 13 . However, dislocation mobility is not restricted to subsonic regime under high strain-rates. Surpass the sonic barrier was first postulated theoretically [14][15][16] and then confirmed by atomistic simulations in bcc tungsten 17,18 and in iron screw dislocations 19 by means of introducing secondary kinks via mother-daughter nucleation process with opposite Burgers vector 20 .
Within the micro magnetic approximation, in a two-dimensional infinite magnetic system, the order parameter that breaks the symmetry is reflected by an average quantity: the magnetisation m(x, y) at each site in a given lattice. The connection between topology and magnetic media and its role in dynamics became apparent by the pioneering work of several authors [21][22][23] when studying the dynamics of magnetic topological defects. The analytical approach consists of assuming the centre of the texture as a dynamical variable where the equation of motion can be derived using Lagrangian formalism. Concretely, in ferromagnets, domain walls (DW) are defined by the transition region that separates two homogeneously and opposite magnetised regions. In a similar fashion as for dislocations, a quantised topological invariant is defined for the magnetisation field 24 . By framing the magnetic texture on a closed path, it allows to define a map of the magnetisation configuration. Distinct topological textures are characterised by different winding number, w, which counts the number of times the magnetisation is wrapped onto itself 25 . In the case of Bloch walls separating magnetic domains in ferromagnets topological defects arise in the form of Néel lines 26 . Since the nineties, the study of DW mobility in ferromagnets has been the flagship in Spintronics because of its potential in data storage 27,28 and logic devices 29 . However, a major difficulty to implement such technological proposals arise as a result of the intrinsic instability the DW structure presents when it overcomes a certain threshold velocity. In this Article we show that this interesting feature of domain walls is present irrespective of whether the magnetic material is a ferromagnet or antiferromagnet. In ferromagnets, this regime is known as the Walker limit 30 and it is characterised by a DW mobility that is the result of the combination of translational and oscillatory dynamics 31 leading to to periodic alteration of the DW chirality. Unfortunately, the Walker breakdown (WB) occurs at a low domain wall velocities and imposes a limit on a reliable control of the DW motion. An alternative to prevent the WB consists of using nanotubes or nano cylinders 32,33 . These particular geometries have unveiled the spin analog of the Cherenkov effect which involves the emission of electromagnetic waves by charge particles moving at speeds faster than the speed of light in the given medium 34 . In the case of antiferromagnets, several studies have provided significant insights on DW motion in antiferromagnets where high speeds have been predicted thanks to prevention of the WB 35, 36 yet there is an upper theoretical limit which can not be surpassed. Interestingly, the manner in which the spin-space is conformed in antiferromagnets obeys the relativistic kinematics of special relativity where the analogue role played by the photons in conventional space is played by the magnons in a magnetic system. Consequently, by comparison, the ultimate limiting velocity of any texture is given by the maximum group velocity of the magnons, v g . An interesting question involves whether, under high-strain, in the same fashion as it occurs for dislocations under plastic deformation, super-magnonic regime of motion could occur through a sequence of additional nucleation textures mimicking the so-called WB. This behaviour would be akin to supersonic motion of dislocations through nano-twin formation. To elucidate the existence of such type of dynamics, we combined one-dimensional atomistic spin dynamics simulations of a the layered antiferromagnet, Mn 2 Au 37 , with one-dimensional analytical theory. In order to induce magnetisation dynamics in Mn 2 Au (see Fig. 1 a), we make use of the predicted staggered field-like torque in such crystal structures 38 , where the effective magnetic field resulting from a staggered induced spin-density H so , possesses opposite signs at each sub-lattice and gives rise a spin-orbit torque (see Fig. 1 b). Some representative dynamics of the temporal evolution of DW position is obtained by extracting the spatiotemporal evolution of the winding number density shown in Fig. 1 c, d. These spatiotemporal evolutions were obtained by firstly ramping the Spin-Orbit field (SO-field) from 0 to its maximum in a time interval of 10 ps and then it is kept constant till the end of the simulation. The intrinsic inertial character of the translational motion of the DW is revealed by a non-linear trend that appears within the first 4-5 picoseconds. The ballistic nature of the DW motion is revealed by the narrowing of the DW as it approaches v g according to the relativistic kinematics framework. The DW width at a given velocity is given by: .8 nm is the DW width at rest. This contraction translates into an increase of the DW exchange energy. It is noteworthy to comment that in conventional relativity the-ory, the contraction observed from non-inertial reference frame does not result in any incremental (or detrimental) interaction among the constituents of the object. First, the DW is accelerated to a speed that approaches to v g . After this transient dynamics, two distinct behaviours are obtained depending on the applied SO-field. For low SO-fields (see Fig. 1 c), a saturation effect of the DW speed is observed, characterised by the slope of the maximum winding number density position (DW position), q, with respect to time, t with v = ∂ t q. A linear variation of the DW speed as a function of the SO-field is consistent with previous reports 36 (see Supplementary Information). In such a case, the energy dissipation induced by the viscous damping inhibits the DW to move faster as reported in previous work 35 . In contrast, at an onset of 65 mT, a different dynamical behaviour is observed which can be defined as the critical stress at which the DW suffers a threshold deformation which leads to the breeding of a DW-pair with trivial composite winding number, w, preserving the overall winding number (see Fig. 1

d, bottom panel) . This
strongly nonlinear regime which involves the generation of additional particles has been reported under several conditions in ferromagnets, for instance, domain wall WB mediated by the vortex-antivortex generation 33 , vortex-core reversal in spin-torque oscillators in nanodots 39 and nanocontacts 40 . However, its occurrence in antiferromagnets has never been reported nor proposed. In ferromagnets, the physical origin for the appearance of new magnetic textures, concretely in the case of vortex-core reversal has been associated to the emergence of a magnetic field whose origin is purely dynamical and can be derived from the kinetic part of the Lagrangian density, L = L kin − V . In our system, the DWs are located in the ferromagnetic basal atomic planes where the coupling between the adjacent planes is antiferromagnetic.
We can derive the temporal evolution of the kinetic field using the Lagrangian formalism.
where M s is the sub-lattice magnetisation. Although this emergent field has already been proposed to explain the reversal of the vortex core polarity in ferromagnets where the dynamics are gyrotropic, we use this formalism also to translational motion of DWs in layered antiferromagnets. In antiferromagnets due to the relativistic nature of the dynamics, the DW already suffers a substantial deformation due to its translational velocity. Interestingly, the extension and the kinetic field profile central location with respect to the moving DW depends upon the dynamical regime the DW is at a given instant (see Fig. 2  where A is the effective exchange stiffness per unit area and K is the uniaxial anisotropy (see Table 1).
The evolution of this reverse magnetic domain is govern by the competition between the exchange and anisotropy energies which results into breaking the reverse magnetic domain into a DW-pair with each of its constituents holding an opposite winding number (thus the sum of their windings is zero). We identify the time involved in the nucleation to be around 1 ps. Different type of magnetic distribution could in principle appear obeying the boundary conditions over the generated magnetic domain, however, we observed that the initially moving DW and the nucleated DW located just in front share the same winding number. The origin of this arrangement can be attributed to the minimisation of the Zeeman energy now comprised of the SO-field and the kinetic field (see Supplementary material). Therefore, while the torque provided by the kinetic field leads to a reverse magnetic domain in front of the moving DW with the appropriate extension to accommodate two DWs, the competition between the magnetic energies and the minimisation of the Zeeman energy gives rise to the generation of a DW-pair with a specific ordering.
The temporal evolution of the system after the DW-pair is present in the system is dramatically affected. Hereafter, we call the initially moving DW, DW 1 and the nucleated DWs, DW 2 and DW 3 being DW 2 always the closest to DW 1 . As can be seen in Fig. 3 a, immediately after the nucleation, DW 1 and DW 2 perform few oscillations and then get stuck evolving in time but not in space. As DW 1 and DW 2 have the same winding number, they repel each other due to the ferromagnetic exchange interaction (see Fig. Fig. 3 a). However, the region in between the pair (DW 1 and DW 2 ) has opposite orientation with respect to H so and therefore, in order to reduce the Zeeman energy provided by the SO-field (Note that the kinetic field is zero as the DWs are static), the distance between DW 1 and DW 2 has to be minimised.
Consequently, as a result of the competition between the two forces, a local minimum appears that corresponds to a certain distance between the DW-pair. We verified this interpretation by calculating the energy barrier related to the distance among two DW 1 and DW 2 (see Fig. 3 b). The energy minimum corresponds to a given separation between the DWs' center, which depends on the magnitude of the SOfield (see Supplementary material). From numerical simulations, we obtain a stable distance of 32 nm whereas from analytical calculation the distance is 29 nm for H so = 65 mT.
Another observation is that DW 3 propagates forward by undergoing a velocity boost that surpass the maximum group velocity of the magnons. This regime is transient and last for several tens of picoseconds which is sufficient for DW 3 to cover over three microns. This mobility regime can be explained as follows. While the DW-pair is nucleated, DW 1 and DW 2 repel each other thanks to the ferromagnetic exchange interaction between them leading to a potential motion in opposite directions for each DW. Despite this fact, they barely move (See inset Fig. 3 a) and therefore, all this latent momentum is transferred from DW 1 to DW 3 41 . Hence, the DW 3 mobility at the moment of the nucleation consist of the velocity provided by the DW 1 , in addition to the boost coming from the momentum transfer due to the repulsion between DW 1 and DW 2 . Using collective coordinates approach combined with the Rayleigh dissipation function and the Euler-Lagrange equation of motion, we obtain the DW 3 velocity as a function of its distance with respecto the de DW-pair to be where α is the Gilbert damping, ∆ DW 1 and ∆ DW 3 are the DWs' width for DW 1 and DW 3 at the nucleation extracted from the atomistic simulations respectively and x(t) is the distance between DW 1 and DW 3 .
The maximum speed according to Eq. 2 is circa 133 km/s which corresponds to a distance of 16 emission of radiation has also been observed 42,43 . We attribute the origin of the emitted spin-waves to a mixture of the so-called Bremsstrahlung effect 44 (or breaking radiation) and the spin-Cherenkov effect 45 .
The Bremsstrahlung effect arises due to the deceleration of a charge particle. Therefore, one could speculate that such staggering deceleration could lead to excitations of the spin medium. As shown in Fig.   3 d, DW 3 not only surpass the maximum group velocity of the magnons but it also enters in a mobility regime where it moves faster than the phase velocity of the magnons. It is difficult in our case to clearly separate the existence of Bremsstrahlung and spin Cherenkov radiation. We note in Figures 1 d and 3 a that the number of oscillations of DW 1 and DW 2 equals the number of spin-wave ripples travelling with DW 3 and further, that the decay of these ripples appear to stretch over a rather long time. Therefore, one possibility for co-existence of both types of radiation could be that the Spin-Cherenkov-originated spin-waves acts in an anti-damping fashion onto the Bremsstrahlung in conjunction with the observation that the radiative ripples to not propagate in opposite direction to DW 3 . It is pertinent to explore how the system reacts when the applied SO-field as the system will have at its disposal additional pumped energy.
We start by injecting an electric current with the time profile illustrated in Fig. 4 a; The electric current has a rising time of 5 ps up to a peak value of H so (t) = 100 mT. The value of H so is kept constant for the following 50 ps before reducing it to zero, with a falling time of 5 ps. The SO-field is set to zero for 50 ps before starting this pattern again three times in succession. The coloured circle in Fig. 4 a, represents the maximum π/2-windings, i/e/ the number of DWs in each pulse. In each cycle, we can observe an avalanche of DW-pairs. Once the primal DW nucleates the first DW-pair, the DW that propagates from the DW-pair (DW 1 and DW 2 ) at supermagnonic speeds becomes a new breeder and gives rise to new DW-pair. This phenomenon repeats for a 13 times leading to the appearance of 26 additional DWs in the system preserving at all times the overall topological charge. It is worth noticing that the DWs once nucleated do not reorganised into a more stable configuration while the SO-field is maximum giving rise to DW-lattice-like structure whose inter-DW distance depends upon the SO-field pulse pattern. The system under the current excitation circumstances could potentially exhibit chaotic behaviour. This will be the topic of future work however. We note at this stage, that recent work on current-induced resistance changes in the antiferromagnet CuMnAs was, in conjunction with imaging attributed to the fragmentation and recovery of the domain structure 46,47 . In particular they observed a gradual increase in the resistance with repeated pulsing and a slow relaxation of the resistance towards lower values than after the last pulse when turning off the excitation completely. However, provided that finite size-effects can be excluded and that such a system observes global conservation of the total winding number, it should be possible to return to the original state by applying a reverse field below the AFM Walker-breakdown field, in order to recover completely the initial state of resistance, in other words, a complete reset of the system. Further, it is noteworthy that the relaxed resistance value after a long waiting time was higher than the original starting resistance. This could mean that residual domain walls are present. If the interpretation that the observed resistance changes is due to domain fragmentation and recovery/recombination is correct we cannot help but to speculate about possible measurements of our system under study. Although the mechanism of domain fragmentation in our case is due to the AFM Walker-breakdown in a perfect crystal and the results in the reported experimental work likely has heat as the main cause of fragmentation this means that that such resistance variations in analogous systems where pinning and heating effects can be ruled out could be used as an indirect detection method for the occurrence of the AFM Walkerbreakdown process and its associated generation of supermagnonic textures. It must be pointed out that for this particular pulse ramp-time no signature of DW-pair nucleation is obtained well after the SO-field has reached H crit. = 65 mT, meaning that the critical field depends upon the ramping time conditions rather than the absolute value of the energy pumped into the system. Fig. 4 b shows the transition from a system with a DW-lattice with 13 DW-pairs towards a system with only one DW-pair due to the multiple recombination when the SO-field is zero. It is noted that there is a distribution of distances that separate consecutive DWs with opposite winding number. This gives rise to a manifold of recombination times as shown in Fig. 4 c . Analytically, it is possible to obtain the recombination time as a function of the distance by assuming that the decompression introduced by the absence of SO-field results in a transition phase from a DW-lattice to a DW-gas for a transient time. Therefore, the recombination time of a given DW-pair is only governed by intrinsic parameters such as, the exchange interaction, the damping and the distance between the DWs ruling out the effect of the presence of the rest of DW-pairs in the system. We can derive the recombination time from a general expression of the exchange interaction between two domain wall pairs (see Supplementary material) with opposite winding number as where ∆ is the DW-width of each DW, l represents the distance between the DWs that will recombine at t 0 , l 0 =1 nm and t r represents the recombination time. We note that the predicted recombination time reproduce quantitatively the recombination time values extracted from the atomistic spin dynamic simulations which validates the hypothesis of 1D DW-pairwise gas approximation. For distances in the range of few hundreds of nanometers the recombination time is in the range of few picoseconds as the exchange interaction is a short range interaction. However we can see from Fig. 4 b that there are DWs remaining in the system whose separation is in the range of few microns. For a given distance of 1.8 µm, the expected recombination time lies in the range of ∼ 6 days. This suggests that in an ideal scenario where pinning effects, thermal fluctuation effects are not present, the stability of such a configuration is granted thanks to the absence of long-range interaction. We can observed in Fig. 4 b that once the SO-field is set to zero, each of the recombination processes produces an excitation in the continuum spectrum. This highly nonlinear process yields to a deformation amplitude that oscillates in time with a very precise frequency. The emergence of this breather-like excitation is due to the excess of kinetic energy carried by each DW involved in the collision whose value is not large enough to escape the attractive potential provided by its anti-particle. By mapping the breather into a simple damped harmonic oscillator, we obtained a good quantitative agreement on the breather lifetime. The breather decay occurs in 20-30 ps and it is governed by Gilbert damping. Moreover, the characteristic frequency of the breather is determined to 544 GHz frequency, which lies within the linewidth of the calculated spin-wave band gap at H SO = 100 mT , i.e. the frequency at zero wave vector, k=0 at. We are thus compelled to assign the breather frequency to that of the band-gap.

Conclusion
In conclusion, we have studied by ASD and theory, the dynamical properties of staggered field-driven DWs. A Walker-breakdown-like process has been identified whereby the DW nucleates an additional DW pair with a total trivial winding. One of these nucleated domain walls forms a π DW together with the original one, whereas the other nucleated DW travels away from the breakdown site as speeds far exceeding the maximum group velocity of spin-waves in the medium, thus achieving supermagnonic speeds. At such high speeds we observe a radiative tail travelling together with the supermagnonic texture. According to the computed dispersion relation, for this system, it is only at supermagnonic speeds where the conditions for spin-Cherenkov radiation is met. Thus one contributing source to the observed radiation, we attribute to spin-Cherenkov radiation. At the same time, we would expect a breaking radiation to be present due to a rapid deceleration of DW 1 and DW 2 . We observe associated oscillations in position of the said DWs as they come to a halt, with the number of oscillations equaling that of the number of radiative ripples travelling with DW 3 , opening the possibility of co-existence of the two types of radiation with the spin-Cherenkov radiation acting in an antidamping manner onto the breaking radiation , which could explain the rather slow decay of the radiative ripples. At higher spinorbit field magnitudes, both multiple DW nucleations and DW recombinations occur, in order to keep conservation of the global winding number. Provided that the separation between nucleated textures is sufficiently large, residual domain walls can remain in the system after the pulse-excitation is turned off.
Provided that the origin of the measured resistance in CuMnAs is of magnetic origin, in particular due to magnetic domain fragmentation 46,47 , our results imply that no Joule heating or temperatures approaching the Néel temperature are required in order to observe such resistance characteristics.

Methods
We perform numerical simulation via atomistic spin dynamics simulations. For this, the full Mn 2 Au crystal structure is taken into consideration. Two formula units as that in Fig. (1a) where γ is the gyromagnetic ratio of a free electron (2.21×10 5 m/As), α is the Gilbert damping set here to 0.001 and H eff i is the effective field resulting from all of the interaction energies. The energies taken into account are the three exchange interactions (two antiferromagnetic and one ferromagnetic), magneto crystalline energy contributions and the spin-orbit field. The total energy, E considered is: The first term on the right-hand side is the exchange energy where J ij is the exchange coefficient along the considered bonds (see supplementary material). The second and third terms are the uniaxial hard and easy anisotropies of strengths K 2⊥ and K 2 , respectively, while the fourth and fifth terms collectively describes tetragonal anisotropy. For the in-plane part of the tetragonal anisotropy, u 1 =[110] and u 2 =[110].
Finally, µ 0 and µ s are the magnetic permeability in vacuum and the magnetic moment, respectively. We have used µ s = 4µ B , with µ B being the Bohr magneton. It is to be noted that the tetragonal anisotropy was included for sake of completeness as it is present in this material but its role in the high speed dynamics is negligible due to the weak magnitude of its anisotropy constants. The effective field is then subsequently evaluated at each point in time using Eq. (5) as H eff i = −1 µ 0 µ S δE δS i . The system of equations,