Tunable terahertz oscillation arising from Bloch-point dynamics in chiral magnets

Skyrmionic textures are being extensively investigated due to the occurrence of novel topological magnetic phenomena, and their promising applications in a new generation of spintronic devices take advantage of the robust topological stability of their spin structures. The development of practical devices relies on a detailed understanding of how skyrmionic structures can be formed, transferred, detected, and annihilated. In this work our considerations go beyond static skyrmions and theoretically show that the formation/annihilation of both skyrmions and antiskyrmions is enabled by the transient creation and propagation of topological singularities (magnetic monopolelike Bloch points). Critically, our results predict that during the winding/unwinding of skyrmionic textures, the Bloch-point propagation will give rise to an emergent electric ﬁeld with a substantial amplitude and in the terahertz frequency range. We also demonstrate ways for controlling Bloch-point dynamics, which directly enable the tunablility on generation of this signal, as well as its frequency and amplitude. Our studies provide a concept of directly exploiting topological singularities for terahertz skyrmion-based electronic devices. DOI: 10.1103


I. INTRODUCTION
Magnetic skyrmions are stable topological solitons with unique textures that cannot be continuously unwound into a collinear ferromagnetic state [1].Various methods of skyrmion formation have been proposed [2,3].Recent demonstrations of the existence of skyrmions at room temperature [4,5] and current-driven motion in racetrack-based nanostructures [6,7] have made skyrmion-based devices promising candidates for future nonvolatile data storage and spintronics-based nanocomputing.Importantly, the study of skyrmionic textures and their formation/annihilation provides a deeper understanding regarding the role of topology in nanomagnetic materials.
In magnetic materials with strong spin-orbit coupling, the competition between antisymmetric Dzyaloshinskii-Moriya interaction (DMI) [8,9] and Heisenberg exchange in conjunction with long-range dipolar interaction generates topologically nontrivial skyrmions.Different classes of DMI arise from different crystal symmetries and therefore a variety of skyrmionic textures can be stabilized, including Bloch [10][11][12] and Néel skyrmions [5,6].Moreover, the DMI in Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.materials with D 2d crystal symmetry supports antiskyrmions whose existence has recently been theoretically predicted [13,14] and experimentally confirmed [15], either as isolated solitons or in regular lattices.A single skyrmion/antiskyrmion has a quantized topological charge (winding number) [1], which is equal to with m the unit vector of magnetization.In particular, the sign depends on the polarity, and skyrmions and antiskyrmions of the same polarity can be classified by the chirality of in-plane magnetization [13].In a recent work [16] it has been demonstrated that the magnetic transition of an ideal hexagonal skyrmion lattice involves a transient antiskyrmion lattice state.During that process, the annihilation of skyrmions mediated by a transient state involving antiskyrmions, into a topologically trivial ferromagnetic state has been observed.These magnetic transitions involve complex dynamic variations of the local magnetization with associated topological winding/unwinding effects.These dynamic local magnetization variations induce changes in the electronic state, because electrons adiabatically moving through complex magnetic textures have their spins continuously adjusted to the local magnetization.Hence, the spatial and temporal variation of the magnetization will induce emergent electrodynamics of conduction electrons [17,18], which can be represented by an emergent magnetic field B e and an emergent electric field E e .In analogy to Faraday's law of induction, for example, a moving skyrmion with a finite velocity will induce an emergent electric field, which directly contributes to Hall measurements as has been experimentally demonstrated by Schulz et al. [19].During the formation and annihilation of skyrmions, the magnetization dynamics should also generate dynamic emergent fields; an idea which forms the core of this paper.These fields provide a signal that could allow for the direct electronic detection of skyrmionic switching processes, as well as technological applications of spintronic devices beyond data storage.The overall formation and annihilation processes involve phenomena that manifest topological singularities along with the associated emergent-field considerations, which can have implications from fundamental physics to novel applications in the terahertz range, and thus more detailed studies are essential.
In this work we start by analyzing annihilation phenomena of single skyrmions and antiskyrmions induced by a sweeping external magnetic field, and we show that within numerical models the transition processes and the topology changes are similar for both textures.During these annihilation processes, skyrmion or antiskyrmion tubes break via the creation and subsequent propagation of topological singularities (monopolelike Bloch points).This propagation exhibits a velocity modulation in the terahertz frequency regime, which in turn gives rise to strong high-frequency emergent electric fields.In addition, we propose a design of nanofabricated defects to enable deterministic controllability of the creation and propagation of Bloch points, and thus provide a route to the engineering of robust devices for ultrafast spintronic applications at terahertz frequencies.

A. Field-induced annihilation of single skyrmions and antiskyrmions
In order to investigate the intrinsic mechanisms of the transition process from skyrmionic to ferromagnetic states, we focus on the magnetic transitions of single skyrmions and antiskyrmions, thus not considering the interactions with other skyrmionic textures.Our numerical simulations are based on ideal thin films of FeGe with noncentrosymmetric crystal structure [10][11][12] (see Appendix D 1 for the material parameters and simulation methods).Bloch-type skyrmionic configurations can be stabilized numerically in both micromagnetic and atomistic simulations of FeGe by the bulk/intrinsic isotropic DMI.We also find metastable single antiskyrmion tubes in micromagnetic simulations of the same material.The magnetization in antiskyrmions has different windings along different radial directions away from the center [13], but bulk DMI only favors either clockwise or anticlockwise orientation of in-plane spins; in our case the stabilization of antiskyrmion tubes is expected to arise from contributions of magnetostatic interactions.
In both single skyrmions and antiskyrmions, when a magnetic field of opposite direction to the core polarity is increased in magnitude, the antiskyrmions annihilate at a smaller field than the skyrmions.This is because their more complex spin texture provides additional reversal routes that are not available to skyrmions, but even so, the two processes have similar features [Figs.1(a) and 1(b), and Video 1].If observed in a cross-sectional view [Fig.1(c), taking a single skyrmion tube as an example], the diameter of the skyrmion tube initially decreases, and then unwinds asymmetrically with the core detaching from the bottom surface, followed by the creation of a Bloch point paired with a chiral bobber [20] [Fig.1(d)]; as the Bloch point propagates through the whole layer towards the top surface, the tube gradually transforms into a trivial uniform ferromagnetic state [Fig.1(e)]; finally, the entire skyrmion/antiskyrmion tube annihilates when the Bloch point disappears at the top surface [Fig.1(f)].It is worth noting that because single skyrmions and antiskyrmions of the same polarity have different winding and opposite topological charge, their annihilation processes are mediated by two different types of Bloch-point singularities [21] [bottom right insets of Figs.1(a) and 1(b)].Moreover, the existence of a metastable chiral bobber has been observed in simulations of the annihilation of single skyrmions, either as isolated skyrmions [20] or in skyrmion lattices [12].Similar magnetic transition phenomena have also been found during the inversion and annihilation in ideal hexagonal FeGe skyrmion lattices (see Appendix A).However, this is the first observation of antiskyrmionic chiral bobbers.

B. Asymmetry in the creation of Bloch points
As our micromagnetic and atomistic calculations show, a single skyrmion/antiskyrmion tube starts to unwind from the bottom of the film during the annihilation process.Dominated by the isotropic DMI, the magnetization profile of a single skyrmion located at the center of the film adopts a Blochtype configuration.However, the spins near free surfaces experience interactions of neighbors from only one side and thus the effect of short-range interactions is weakened.This induces an inhomogeneous magnetization profile near the film surface [22].If the DMI constant is positive the magnetization will turn outwards near the top surface and inwards near the bottom surface, resulting in a pure Bloch-type skyrmionic structure in the middle of the film, and a mix of a Blochand Néel-type skyrmionic structures near the two surfaces [Figs.9(d)-9(f) in Appendix B].
As the applied magnetic field is increased in strength, the magnetization tends to align with the field direction in order to minimize the Zeeman energy, while the precession term makes the in-plane magnetization precess in the right-handed manner around the field direction, which increases the DMI energy at the bottom surface but decreases it at the top surface (see analytical explanations in Appendix B).At a critical field with a magnitude B c , the annihilation process is initiated, and the energy profile at the bottom contributes to a detachment of the skyrmion tube from the surface and creates a Bloch point.Similarly, if the DMI constant is negative then a Bloch point is created at the top surface and propagates downwards (Video 2).
In order to test the stability of the observed behavior against finite-temperature effects, we introduced a randomly fluctuating thermal field corresponding to T = 200 K (see Appendix D 2), which is below the ferromagnetic Curie temperature of 280 K for FeGe [23].We find that the annihilation follows a similar process (Fig. S1 in the Supplemental Material [24], and Video 3) as that at zero temperature, and as expected the critical field of the annihilation process is reduced to B c = 1.54 T, because thermal fluctuations assist in overcoming the energy barrier.

C. Detectable terahertz emergent electric field
As the annihilation process of the skyrmion tube is mediated by Bloch-point dynamics, it involves a rapid change of topological charge: the nontrivial skyrmionic configuration near the surface folds into a topologically equivalent Bloch-point singularity, and the topological charge density rapidly decreases when the Bloch point propagates through the film and leaves a trivial ferromagnetic state in its wake (see Appendix D 3 for the calculation of topological charge density).Importantly, the Bloch point is not localized in a single simulation cell because at the center of the Bloch-point texture the magnetization vanishes, and the assumption of a continuous magnetization field breaks down [1].Hence, for discrete atomistic simulations the center of a Bloch point has to lie between two atoms, since the magnetic moments of the atoms are assumed to be constant in magnitude.This is also valid in micromagnetic simulations where the Bloch-point center is among computational cells.
For single crystals with regular atomic arrangement, such as FeGe [10][11][12] that we simulated in this paper, the lattice structure induces potential wells among the neighboring atoms [25].When a Bloch point propagates through the film, the energy barriers among the potential wells have associated forces that cause the Bloch-point velocity to be modulated.We have performed atomistic simulations to model the dynamics of classical spins and their related magnetic moment in the FeGe system (see Appendix D 1 for methods of the atomistic simulations).The results show that when a propagating Bloch point experiences three-dimensional periodic potentials, its position/displacement also exhibits corresponding oscillations in the x, y, and z directions, rather than simply moving upwards with uniform velocity [Figs.2(a A real-space Berry phase can be regarded as the solid angle subtended by the spin of a conduction electron that is continuously adjusting to the direction of the local magnetization m [18].In a temporally inhomogeneous magnetic texture, the Berry phases result in emergent electromagnetic fields acting on the conduction electrons.The emergent electric field in direction i = x, y, z is given by For the specific case of a Bloch point, which is a topological point defect [1,26], we obtain a quantized topological charge upon integration over the unit sphere S 2 .When a conduction electron traverses a magnetic texture with a moving Bloch point, it experiences a net force F = qE e with magnitude F ≈ hν BP /2λ 2 , where v BP is the velocity of a Bloch point and λ is its approximated diameter (we assume λ = 1 nm in our simulations) [26].The effect of the emergent electric field on the conduction electrons is thus similar to that of a real electric field.Its strength is linearly dependent on the propagation velocity, and can be easily calculated using the temporal evolution of the Bloch-point position.We find that the velocity v BP is on the order of several kilometers per second and its associated electric field has a magnitude on the order of megavolts per meter [Fig.2(b)], similar to the values found in permalloy nanowires as reported by Charilaou et al. [26].The potential wells originating from the periodic lattice structure discussed above induce rapid oscillations of the Bloch-point velocity and thus the associated emergent field.The fundamental frequency of the emergent electric field extracted with the fast Fourier transform (FFT) method, lies in the terahertz range [Fig.2(c)].
During the final collapse of the structure the Bloch point moves at a velocity of several kilometers per second in a direction perpendicular to the plane of the magnetic layer, but the periodic lattice potential [inset of Fig. 2(a)] causes oscillations and the consequent emergent field of strength ∼2 MV/m is similar in all three dimensions [Fig.2(b)].This is different from previous related experimental work [19], which focused on a study of the electric field arising solely from the in-plane motion of entire skyrmions where the electric current is required to obtain a temporally varied magnetic texture.Here, in contrast, when the skyrmion annihilates via the Bloch-point-mediated mechanism, the radiating emergent electric field is solely excited by a uniform magnetic field.This suggests that the formation and annihilation processes of single skyrmions may be electronically detectable, and it also gives rise to a new mechanism for the generation of terahertz-frequency emergent electromagnetic fields.
In a micromagnetic model, one averages over a large number of atomic spins and thus the Bloch-point behavior and the resulting emergent field cannot be accurately predicted.However, the mesh-friction effect in micromagnetics [27], where potential wells are induced by finite-size meshes to discretize the magnetization in practical calculations may cause similar phenomena (see micromagnetic results in Appendix c), i.e., the Bloch point exhibits oscillatory motion that can generate an associated electric field [Figs.10(a) and 10(b) in Appendix C].As the mesh-friction effect that originates from the discrete lattice structure is relevant to the lattice constant, different cell size in micromagnetics compared to the lattice constant may cause the pinning potential to be significantly overestimated/underestimated.Therefore, we studied the effect of different cell sizes in micromagnetic simulations [Figs.10(d) and 10(e) in Appendix C].In these simulations the values of the fields are on the same order of magnitude as those predicted by the atomistic model, demonstrating that both models are able to produce comparable results.In particular, when the cell size is comparable to the lattice constant a 0 , the critical fields of the annihilation processes are very similar (1.86 T in the atomistic simulation and 1.82 T in the micromagnetic one), which suggests that the energy barriers in the two models produce comparable effects.Although it is generally difficult to assess the micromagnetic accuracy of the Bloch-point velocity modulation, such models can show trends upon the variation of the material parameters.Since micromagnetic simulations are computationally much more efficient than atomistic simulations, especially for the calculation of demagnetization fields in large systems, we mainly present micromagnetic models for the systematic investigations to reveal more Bloch-point dynamics in the subsequent sections of this paper, assuming qualitatively similar results compared to those obtained from atomistic simulations.

D. DMI and thickness-dependent behavior of Bloch points
As described above, we attributed the asymmetries in the creation of the Bloch point and the annihilation of the skyrmion tube to the distribution of the DMI energy profile through the film driven by the sweeping magnetic field.The surface at which a Bloch point is created is determined by the sign of the DMI in the material, i.e., if the DMI constant is positive, a Bloch point is created at the bottom surface and propagates upwards [Fig.3(a The magnetostatic energy plays a crucial role in the stabilization of skyrmionic textures in three-dimensional systems.In order to investigate the influence of the magnetostatic energy on the position where the Bloch point is created, we varied the film thickness and divided the system into 1-nmthick layers along the z axis.The corresponding magnetostatic energy of each layer is calculated at the applied field B ext with a magnitude 1.63 T, which is less than B c for any thickness (Fig. S2 in the Supplemental Material [24]).The demagnetization field provides greater stability for a skyrmion structure near the film surface than in the center of the film, whereas isotropic DMI provides greater stability for a skyrmion in the center of the film than at its surface [28,29].Thus, the combined effect of the two interactions results in an inhomogeneous energy distribution throughout the film, where the location of the maximum energy depends on the film thickness.
For a thin film (thickness below 45 nm in our case, around 0.6L D with FeGe helical period L D = 70 nm [30]), the magnetostatic energy inside the film is lower than at the surface, and the Bloch point is created at the bottom surface for positive DMI [Fig.3(a) and Video 4] or the top surface for negative DMI [Fig.3(b)], as we showed above.For a thick film (thickness above 55 nm in our case, around 0.8L D ), the energy is significantly higher in the middle than at the surface and a pair of Bloch points emerges close to the middle [Fig.3(d)].Then in the intermediate regime (thickness between 45 and 55 nm in our case, around 0.6 − 0.8L D ), the magnetostatic energy in the middle layer is comparable to the energy at the surface, and consequently, the skyrmion tube becomes thinner in the middle and a pair of Bloch points is created close to the middle at the same time as an individual Bloch point emerges at the surface [Fig.3(c)].It is worth noting that the Bloch-point pair is not created exactly in the middle of the film, but slightly closer to the bottom surface.This can be attributed to the influence of a positive DMI value, and conversely the Bloch-point pair emerges closer to the top surface for a negative DMI.On the other hand, in DMI-free systems, including centrosymmetric magnets with frustrated exchange interactions [31,32] which favor both skyrmions and antiskyrmions, the dynamic magnetization precession does not induce energy difference on the top and bottom surfaces.Then the creation of Bloch points will start in the middle of the film, dominated by magnetostatic interactions, and the phenomena will be similar to the Bloch-point-mediated switching processes in traditional magnetic bubbles, which have previously been reported [33].Here we demonstrate that the site of the Bloch-point creation is determined by the relative amplitudes of the thickness-dependent magnetostatic energy and the DMI energy term, elucidating the intimate relation between Bloch points and the destruction of skyrmion lines.

E. Tunable Bloch-point dynamics by Gilbert damping and magnetic field
Since Bloch-point creation and propagation are a series of dissipative dynamic magnetization phenomena, their behavior depends on the Gilbert damping parameter α, which is known to slow down the depinning of magnetic domain walls [34].In this paper, we use α = 0.3 for FeGe [35], but we have also analyzed these phenomena for a wider range of α to look at the universality of our results (Figs. 4(a) and 4(b), and Figs.S3(a) and S3(c) in the Supplemental Material [24]).Particularly noteworthy are the cases of ultralow damping, α = 6 × 10 −5 , where we find a large amplitude for the oscillatory phenomena on the skyrmion tube during the Bloch-point propagation [Fig.4(a)], and that of critical damping α = 1, where we observe a direct energy-minimization process [Fig.4(b)].The value α = 6 × 10 −5 is comparable to the damping parameter of sub-100-nm-thick yttrium iron garnet (YIG), which has the lowest known damping of any thin film [36].For a low damping (α < 0.1 in our case), the amplitude of the magnetization oscillation before the Bloch point approaches the top surface is sufficiently large to detach the skyrmion tube from the top surface.This leads to the creation of a new Bloch point (α = 0.01, from the second to third panel in Fig. S3(a) in the Supplemental Material [24]) below the top surface, or even a Bloch-point pair for an ultralow damping [α = 6 × 10 −5 , from the second to third panel in Fig. 4(a)].
We find that the Bloch-point propagation velocity increases with increasing α and saturates close to the critical damping [Fig.4(d)].This is to be expected since increasing the damping increases the rate at which the magnetization moves towards the minimum energy in the Landau-Lifshitz equation.This dependence of velocity on the damping parameter is opposite to the one observed in a current-driven skyrmion longitudinal motion along a racetrack [37], where the velocity is inversely proportional to α due to a different mechanism.It is worth noting that in the low damping region (α < 0.1), complex processes such as the creation of additional Bloch points or pairs do not significantly affect the average velocity.
The strength of the applied magnetic field affects the Bloch-point propagation velocity more straightforwardly and has a greater effect than damping (Figs.S3(d) and S3(e) in the Supplemental Material [24]).Whereas the velocity between the cases of ultralow damping and critical damping increases by only 10%. Figure 4(e) shows that for α = 0.3 the propagation velocity linearly increases with the magnetic field strength and can be doubled by applying a field that is around 15% higher than B c .It is worth noting that at around B ext = 2.25 T, there is a sudden velocity jump attributed to another Bloch-point pair breaking off from the chiral bobber as the skyrmion tube unwinds [Fig.4(c)].Interestingly, when a high magnetic field B ext = 1.80 T is applied in the ultralow damping case, two Bloch-point pairs are created, i.e., one Bloch-point pair during the propagation due to the application of the field, and another pair just below the top surface similar to the behavior of the ultralow damping case (Fig. S3(c) in the Supplemental Material [24], and see all the above-mentioned dynamic processes in the Video 5).
The Bloch-point velocity can be modulated by applying different strengths of the magnetic field, and the oscillation frequency of the emergent electric field is determined by the oscillatory velocity of the Bloch point when it propagates through a tilted "washboard" potential originating from the discrete nature of the crystal lattices [Figs.5(a) and 5(b)].Therefore, as the applied field increases so does the Bloch point velocity and the associated frequency.Their mutual relationship is illustrated in Fig. 5(c), which shows a linear dependence.Our findings, therefore, show that the Bloch-point propagation during the annihilation of skyrmions will result in a terahertz-level oscillating emergent electric field with a duration of several picoseconds, and the frequency of which can be controlled by an external magnetic field, prompting us to envisage possible applications in future spintronics devices.

F. The effect of an engineered notch
Racetrack systems have previously been designed with in-plane nonmagnetic notch structures to enable geometrically localized disordered spin textures and control the lateral position of domain walls, which may result in domain wall pinning [38] or skyrmion formation [39] when driven by electric currents.Here we propose using three-dimensional notch structures to realize the manipulation of Bloch-point dynamics, since the site of the Bloch-point creation can be effectively controlled by the notch structure.In this way we propose a new mechanism where the creation of Bloch-point singularities as well as the writing and deleting of individual skyrmions can be achieved by applying a magnetic field of appropriate strength.
We have simulated the effect of a conical notch near the top surface.The notch has a diameter of 10 nm and a depth of 5 nm, and the saturation magnetization M s was set to zero inside the notch region (Fig. 6).To determine the effect of the notch we studied the behavior of skyrmions in a uniform magnetic field, starting from the uniform ferromagnetic state, and then swept the external field from B ext = 1.50 to 0.50 T. The artificial defect changes the geometry of the system, effectively modifying the demagnetization field surrounding the notch region and consequently changing the distribution of magnetostatic energy density throughout the film.The local magnetostatic fields cause the spins to twist around the nonmagnetic notch, causing an accumulation of topological charge density around the notch which lowers the energy barrier for the creation of a Bloch point.In Fig. 6 we plot both the local magnetization (m z ) and the topological charge density (Q dens ) to illustrate the dynamic winding and unwinding of the spin textures around the notch structure.At B ext = 0.67 T, it detaches from the notch and creates a Bloch point, which propagates downwards until the whole skyrmion tube finally has switched back to the trivial ferromagnetic state.The global magnetization and topological charge are monitored during the whole notch-induced skyrmion formation and annihilation processes (Fig. S4 in the Supplemental Material [24], and Video 6).The results show that even in a material with a positive DMI, Bloch points will emerge at the tip of the notch near the top surface, in contrast to the creation at the bottom surface in a notch-free system.Hence, a notch structure allows us to tune the balance between DMI and magnetostatics.
Even though the notch depth (5 nm) is much smaller than the thickness of the films that we have studied (20-250 nm), the modifications to the magnetization profile and associated magnetostatic energy are substantial.The twisted magnetization state around the notch enables annihilation to occur at significantly lower fields than in films without the notch.Taking the 70-nm-thick sample as an example, the skyrmion tube starts to annihilate at B c = 1.66 T in the notchfree system [Fig.3(d)], in contrast to B c = 1.30T from the top side in the notch-based system [Fig.6(c)].Meanwhile, at this thickness the notch-free film would normally create a Blochpoint pair inside the film resulting in a complex emergent electric field with two superimposed oscillation frequencies, as the two Bloch points will propagate in opposite directions with slightly different velocities in the initial stage (Fig. S5(a) in the Supplemental Material [24]).In contrast, the notch generates only a single Bloch point that propagates through the film, giving rise to an emergent field with a more uniform frequency and amplitude (Fig. S5(b) in the Supplemental Material [24]).We also varied the film thickness to show that the critical field of annihilation processes in the notchbased system is almost independent of the thickness, i.e., the variation in the field strength [light blue region in Fig. 6(d)] is much smaller with a notch compared to the notch-free structure [light red region in Fig. 6(d)].
In order to show that these phenomena can be mainly attributed to the distortion of the spin texture around the notch, we have investigated how the depth of the notch structure affects the switching fields for the formation and annihilation of skyrmions.Figure 6(e) shows that at small notch depths, the surface charge and a resulting magnetostatic field have an influence on the magnetostatic energy distribution [40], which weakens the notch-induced localized distortion.However, the local distortion inside the system does not vary significantly if the notch depth is large enough (greater than 6 nm in our case), and the switching fields become almost independent of the notch depth.These findings show that relatively deep notches will provide sufficient stability for the control of Bloch-point creation as well as skyrmion writing/deleting in order to make the fabrication robust against engineering errors in the notch size.

III. CONCLUSION
In conclusion, we have identified the mechanisms of skyrmion annihilation in chiral magnetic thin films with isotropic Dzyaloshinskii-Moriya interaction.Below a certain film thickness, skyrmions are stable and continuous throughout the film.Upon application of an external field, a skyrmion tube breaks and terminates in the form of Bloch points that propagate below the film surface.In thicker films, the skyrmions break inside the film and pairs of Bloch points are created that move in opposite directions towards the surfaces of the film.In both cases, the Bloch points move at speeds of kilometers per second.Surprisingly, the speed is modulated at terahertz frequencies due to the rapid hopping of the Bloch point along the energy wells of the discretized lattice.The Bloch-point motion also generates emergent electric fields with a substantial and consequently measurable magnitude in the megavolt per meter range.
Our parametric studies also reveal the underlying mechanisms of Bloch-point dynamics, which are attributed to the competition among DMI, magnetostatic energy, and external magnetic field strength.We demonstrate that the frequency of the emergent electric field is proportional to the Bloch-point propagation velocity, which in turn depends on the external magnetic field.We propose that nanofabricated defects can deterministically localize topological charge, allowing us to tailor the creation and propagation of Bloch points, thus realizing controllable generation of the high-frequency electromagnetic signal.Since Bloch-point singularities are the smallest topological defects in spintronic systems, it is tantalizing to envisage their possible use in future electronic devices that exploit their emergent intriguing dynamics.system when the magnetic field is increased [Fig.7(b)].Despite the opposite winding of skyrmions and antiskyrmions, the annihilation of both the antiskyrmions [Fig.7(c)] during the first step inversion and the skyrmions [Fig.7(e)] in the second step occurs in a similar manner because the cores of both objects are polarized in the direction opposite to the field.Their diameters initially decrease and then, if observed in a cross-sectional view [Figs.7(d) and 7(f)], can be seen that both unwind asymmetrically, with their cores detaching from one surface and propagating through the whole layer towards the opposite surface, leaving a trivial ferromagnetic state behind.

APPENDIX B: AN ANALYTICAL MODEL FOR DMI-INDUCED ASYMMETRIC PHENOMENA DURING SINGLE-SKYRMION ANNIHILATION
In order to better understand the role played by DMI in asymmetric annihilation phenomena, an analytical model is developed with the energy density written as [1,41] where m is the magnetization unit vector, A is the exchange stiffness, μ 0 is the vacuum permeability constant, and H indicates the externally applied magnetic field.Here we use the function w DMI to indicate the DMI energy density, as the ratio of contributions between Rashba spin-orbit coupling (SOC) D ⊥ and Dresselhaus SOC D varies throughout the layer [42].Since a single skyrmion in an infinite thin film has a centrosymmetric profile, it is reasonable to introduce cylindrical polar coordinates r = (rcosψ, rsinψ, z) and spherical polar coordinates to describe the normalized magnetization m = (sinθ cosφ, sinθ sin φ, cos θ ).For each position, m is then described by θ (r, ψ ) and φ(r, ψ ) and Eq.(B1) can be reduced to where the subscripts indicate the partial derivatives.Because of the rotationally symmetric magnetization profile, the minimization of Eq. (B2) results in the solutions θ = θ (r) and φ = φ(ψ ).In order to explain the asymmetric phenomena through the layer from a viewpoint of the analytical model, we will in the following solve the DMI energy for a given z position by calculating the exact form of θ = θ (r) and φ = φ(ψ ).

Analytical solution of θ(r)
The out-of-plane component of the skyrmion magnetization has been shown to be similar to the profile of a 2π domain wall [1,43].Based on the results known for one-dimensional nanowires, the form of θ can be described by the ansatz [1,43,44] where R s is the variational parameter and δ s is the skyrmion radius.Fitting Eq. (B3) to the micromagnetic results (Fig. 8) shows that in 30-nm-thick films, at B ext = μ 0 H = 1.00 T, the out-of-plane profile does not vary significantly with the position in the layer.This suggests comparatively similar contributions to the DMI energy from the out-of-plane component at each value of z.

Analytical solution of φ = φ(ψ)
For a skyrmion of given chirality [Figs.9(a)-9(c)], the in-plane profile can be written as [41] with φ 1 expressing the chirality and favored by Lifshitz invariants that determine the energy via [45] and The in-plane profiles have been extracted from the micromagnetic results and are shown in Figs.10(d)-10(g).The middle layer of the skyrmion tube behaves like a pure Bloch skyrmion with φ 1 ≈ 90 • , and a Néel component is increasingly admixed when approaching the surface.

Calculation of DMI energy through the film
In order to analyze the variation of DMI energy through the whole film [Fig.9(h)], the volume integrals are calculated using Eq.(B7) over a layer of computational cells, with the thickness at a given z value being equal to the cell size of 1 nm.In the middle of the film (z = 15 nm) we have the lowest DMI energy because the in-plane magnetization has the profile of a bulk isotropic system.During the magnetic-field sweeping process, when the field is increased, dynamic magnetization precession attributes to the torque of the effective field and induces a clockwise in-plane rotation (when viewed from the top of the system).The precessional rotation simultaneously increases the in-plane angle φ 1 for all the layers, which effectively lowers the DMI energy of the top layer and increases that of the bottom layer [inset of Fig. 9(h)].In our simulations the increase in applied field is 0.001 T/step, i.e., around 0.05% of the critical field.The DMI-induced asymmetric energy change is balanced by various other energy terms (such as exchange and magnetostatic energy), and the annihilation starts from a layer where the energy increase is no longer compensated by other terms.

APPENDIX C: MICROMAGNETIC SIMULATIONS OF SINGLE-SKYRMION ANNIHILATION, AND ITS ASSOCIATED BLOCH-POINT DYNAMICS
The process of single-skyrmion annihilation is mediated by the creation and propagation of Bloch points.A Bloch point is defined as a monopolelike singularity in continuum theory, and based on discrete meshes with a fixed length of magnetic moments.The Bloch-point center is located between the mesh points.In atomistic simulations, during the Blochpoint propagation, the velocity periodically oscillates due to the discrete potential wells between neighboring spins caused by the crystal-lattice structure.In micromagnetics, similar to the atomistic model, the potential well can be induced by a mesh-friction effect [27], where the magnetization is discretized by the effects of the finite-size meshes.It causes apparently similar phenomena, i.e., the Bloch point exhibits oscillatory motion and an associated electric field is produced.Thus, in order to compare the micromagnetic and atomistic models, we performed micromagnetic simulations of singleskyrmion annihilation processes (without taking into account the demagnetization field in order to keep consistency with the atomistic simulations).When the cell size is equal to the lattice constant a 0 , the Bloch-point propagation exhibits periodic oscillations with a velocity on the order of a few From the FFT-amplitude spectrum, the peak frequency of the velocity oscillations is on the order of a few terahertz [Fig.10(c)].
When the cell size in micromagnetic simulations is much different from the lattice constant, the energy around the Bloch points may be significantly overestimated/underestimated, which further influences the energy barrier of magnetic switching and the mesh-friction effect.Therefore, we conducted simulations with different cell size to study these mesh-dependent phenomena [Fig.10(d)].All of the results are on the same order of magnitude as those from atomistic models (Fig. 2 in the main text), demonstrating that both models yield qualitatively comparable results.Especially when the cell size is close to the lattice constant a 0 , the critical fields of annihilation processes in the micromagnetic model reveal similar values (1.82 T) compared to those from the atomistic model (1.86 T).This suggests that both models produce comparable energy barriers for the magnetic switching.It is worth noting that the frequency observed in the micromagnetic model is linearly related to the critical field [Fig.10(e)], in agreement with Figs.4(e) and 5 in the main text, which show a similar relationship at different external fields but with the same cell size.This illustrates that the frequency intrinsically depends on the magnetic-field excitation, while the variation of cell size only changes the degree of the mesh friction (energy-barrier height), resulting in a different critical field.Our comparisons reveal that it is feasible to use micromagnetic models for parametric studies of Bloch-point dynamics.

Simulation methods and parameters
The dynamic micromagnetic simulations were carried out using Mumax 3 , a finite-difference GPU-accelerated package [46], whose underlying physics is based on the Landau-Lifshitz equation of motion [47]: where m is the magnetization unit vector, γ is the gyromagnetic ratio of the electron, and α is the dimensionless Gilbert damping parameter.Here we used α = 0.3 [35] if not specified otherwise.Furthermore, H eff is the effective magnetic field experienced by the magnetization M = M s m with constant magnitude M s , and is defined as H eff = −μ −1 0 δE /δM with the vacuum permeability constant μ 0 .The energy density E v is written as which consists of: (1) a symmetric exchange interaction with exchange stiffness A; (2) an isotropic Dzyaloshinskii-Moriya interaction (DMI) for crystallographic classes of D n symmetry as a sum of Lifshitz invariants [48] with a strength D; (3) a magnetocrystalline anisotropy with the first-order uniaxial anisotropy constant K u1 and its direction n; (4) a Zeeman coupling of the magnetization with the external magnetic field H; and (5) a magnetostatic interaction energy of the magnetization in the field H d generated by volume and surface magnetostatic charges.
In our micromagnetic simulations, Eq. (D1) is solved using the Dormand-Prince method (RK45 solver) with adaptive time steps.The magnetic parameters were adapted from the cubic FeGe chiral magnets [35], with exchange stiffness A = 8.78 pJ m −1 , DMI strength D = 1.58 mJ m −2 , and saturation magnetization M s = 384 kA m −1 .Bulk FeGe has an isotropic lattice structure and consequently no magnetocrystalline anisotropy.Although the existence of anisotropy has been experimentally validated in epitaxial thin films [12,49], our initial simulations indicate that the experimentally observed values do not produce qualitatively different changes on the simulated behavior of skyrmions.Furthermore, because the anisotropy of thin films is relevant to the substrate and deposition methods, we used the bulk values.The overall system was divided into 1 nm × 1 nm × 1 nm cubic cells.Infinite samples were simulated by setting the dimensions in the xy plane to 200 nm × 200 nm, with periodic boundary conditions in the plane, and in the z direction the default thickness is 30 nm unless specified otherwise.
The atomistic methods are based on the open-source spin simulation framework-Spirit [50].The ideal crystalline B20 FeGe is a simple cubic lattice structure with a lattice constant a 0 = 0.47 nm and four Fe atoms in the P2 1 3 space group, in Wyckoff position (x, x, x), (−x + 1/2, 1 − x, x + 1/2), (−x, x + 1/2, −x + 1/2), and (x + 1/2, −x + 1/2, −x), with x = 0.3865 [51].We calculated the exchange coupling between neighboring Fe atoms necessary to achieve consistency with the micromagnetic model as J = Aa 0 /n = 6.45 meV/pair, where n is the number of atoms per unit cell [52].We took n = 4 and only considered the first-nearest-neighbor exchange interactions.The DMI strength between first-nearest-neighbor spins is D a = 2π a 0 J/L D = 0.27 meV/pair, where the helical period of FeGe L D is 70 nm [30].The local spin moment of Fe is μ s = M s a 3 0 /n = 1.08 μ B /atom, similar to the theoretical value of 1.16μ B /atom [53], with μ B the Bohr magneton [52].The atomistic simulation system was 50 × 50 × 20 unit cells (thickness 9.4 nm) with periodic boundary conditions in the x and y (in-plane) directions.

Finite-temperature simulation
The stability of our results against nonzero temperatures is tested by introducing a randomly fluctuating thermal field to H eff , and is determined by Brown [54]: with a randomly oriented normal vector η, the Boltzmann constant k B , the temperature T (200 K in our simulations), the volume of single computational cell V, and the time step t.
The simulations are implemented using the sixth order Runge-Kutta-Fehlberg (RKF56) solver with adaptive time steps [55].

Calculation of topological charge density
The topological charge denotes the winding of spin texture in a mapping S 2 → S 2 , and the two-sphere S 2 can be either the compactified plane R 2 or a sphere surrounding a point [1].For a 3D skyrmion tube, it has a quantized charge Q = 1 in each slice of the xy plane perpendicular to the film.For an illustration of the topological-charge evolution during the annihilation process of a skyrmion/antiskyrmion tube, we plotted the topological charge density Q dens = 1 4π m × (∂ x m × ∂ y m).In our discrete micromagnetic and atomistic system, the equation is implemented with a MATLAB script, using a higher-order polynomial interpolation to approximate the partial derivatives ∂ x m and ∂ y m.The lattice version of topological charge introduced by Berg et al. [56] has also been implemented by calculating the solid angle from each neighboring triangle magnetization group.The results obtained from both methods exhibit a high level of consistency, confirming the accuracy of our calculations.

Calculation of Bloch-point position and velocity
For the calculation of the Bloch-point propagation velocity, the localization of the Bloch point is of vital importance.In our simulation we calculate the Bloch-point center using the method proposed by Thiaville et al. [27], in which the location of the Bloch point is determined by a polynomial interpolation of the magnetization from nearby mesh points.

FIG. 1 .
FIG. 1. Bloch-point dynamics during the field-induced annihilation of single skyrmions/antiskyrmions in FeGe materials.The magnetization is shown with arrows and the topological charge density is represented by the background color.A chiral bobber (dashed-line region) followed by a Bloch point (circled region) is created from a single (a) skyrmion and (b) antiskyrmion.The top view of magnetization is shown for a rectangular region in the top-right corner, and a schematic view of the Bloch point is illustrated in the bottom-right corner.(c)-(f) Snapshots of the dynamic annihilation process of a single skyrmion computed by atomistic simulations.
) and 2(b), see Appendix D 4 for the calculation of Bloch-point position and velocity].

FIG. 2 .
FIG. 2. Oscillatory Bloch-point motion during a single-skyrmion annihilation.(a) Temporal evolution of the Bloch-point displacement along x (orange line), y (blue line), and z (black line) direction.The displacements along x and y direction are also shown in the enlarged inset.(b) Velocity of the Bloch-point propagation with respect to (a), and associated emergent electric field.(c) FFT-amplitude spectrum of the velocity oscillation in (b).Note that the peak frequencies are in the terahertz range.

FIG. 3 .
FIG. 3. Thickness dependence of Bloch-point dynamic processes during skyrmion annihilation.(a) and (b) For a film thickness of 30 nm, the Bloch points emerge from (a) the bottom surface for positive DMI, and (b) the top surface for negative DMI. (c) For a thickness of 50 nm a pair of Bloch points is created inside the film together with an individual Bloch point at the bottom.(d) For a thickness of 70 nm, a pair of Bloch points is created inside the film.
)]. Conversely, if the DMI constant is negative, a Bloch point is created at the top surface and propagates downwards [Fig.3(b), and Video 2].

1 )FIG. 4 .
FIG. 4. Effect of Gilbert damping and external magnetic field on the Bloch-point behavior.(a)-(c) The number of Bloch points that are created and propagate depends on α and B ext .(d) and (e) Dependence of the Bloch-point velocity on (d) the damping parameter at B ext = 1.71 T and (e) the external magnetic field at α = 0.3.The number of Bloch points N BP is indicated by the symbols as shown in the legends.The solid lines connecting the data points are guides to the eye.

FIG. 5 .FIG. 6 .
FIG. 5. Dependence of the peak-oscillation frequency on the Bloch-point propagation velocity.(a) Schematics of the Bloch-point propagation with a tilted washboard potential originating from discrete crystal lattices.This results in (b) oscillation of the propagation velocity when a Bloch point is hopping through the film.(c) Peak frequency of the velocity as a function of the Bloch-point propagation velocity.The solid lines connecting the data points are guides to the eye.a Bloch point is created below the notch tip, and as the Bloch point propagates downwards, a chiral bobber is generated and merges with the propagating Bloch point until it reaches the bottom surface, where the Bloch point disappears and a com-

FIG. 7 .
FIG. 7. Topological phenomena in an FeGe skyrmion lattice.(a) Hexagonal skyrmion lattice at zero field.Inset: Enlarged panel showing that each skyrmion is surrounded with six skyrmionlike and six antiskyrmionlike configurations.(b) Field-induced magnetic transition of a skyrmion-lattice inversion involving the appearance of transient antiskyrmions at 0.73 T. The topological charge density distribution of the enlarged region is shown in the top-right corner.(c) Magnetization profile of an individual antiskyrmion and (d) its dynamic annihilation process at 0.73 T in cross section.(e) Magnetization profile of an individual skyrmion and (f) its dynamic annihilation process at 1.42 T in cross section.The simulation time during the annihilation process is shown in the top-right corner in (d) and (f).The magnetic field is swept from 0 to 1.42 T in the direction of the z axis.

FIG. 8 .
FIG.8.Fitting of the skyrmion out-of-plane profile with film thickness 30 nm.Skyrmion positioned (a) at the top with z = 30 nm, (b) in the middle with z = 15 nm, and (c) at the bottom with z = 1 nm.The ansatz is chosen from Eq. (B3), and applied field is B ext = 1.00 T.

)
Here the C n symmetry corresponds to the linear combination of Eqs.(B5) and (B6): w DMI = D 1 w 1 + D 2 w 2 with the angle φ 1 = arctan(−D 2 /D 1 ) and effective DMI strength D = D 2 1 + D 2 2 .As a result, we can simplify the general DMI energy density at a given value of z as w DMI (θ, r, φ 1 ) = −D θ r cos2φ 1 +

FIG. 10 .
FIG. 10.Bloch-point dynamics implemented by micromagnetic simulations.(a) Temporal evolution of the Bloch-point position.(b) Blochpoint velocity and associated emergent electric field as a function of time.(c) FFT-amplitude spectrum of the velocity in (b).(d) Critical fieldof annihilation processes as a function of the computational cell size (in the unit of the FeGe lattice constant, a 0 = 0.47 nm).The critical field shows a roughly logarithmic relationship to the cell size.When the micromagnetic cell size (black diamond) is equal to a 0 , the critical field is 1.82 T, while that from atomistic simulation (red diamond) is 1.86 T. (e) Velocity-oscillation frequency versus critical field for different computational cell sizes.The solid lines connecting the data points are guides to the eye.kilometers per second, with an associated emergent electric field on the order of megavolt per meter [Figs.10(a)and 10(b)].From the FFT-amplitude spectrum, the peak frequency of the velocity oscillations is on the order of a few terahertz [Fig.10(c)].When the cell size in micromagnetic simulations is much different from the lattice constant, the energy around the Bloch points may be significantly overestimated/underestimated, which further influences the energy barrier of magnetic switching and the mesh-friction effect.Therefore, we conducted simulations with different cell size to study these mesh-dependent phenomena [Fig.10(d)].All of the results are on the same order of magnitude as those from atomistic models (Fig.2in the main text), demonstrating that both models yield qualitatively comparable results.Especially when the cell size is close to the lattice constant a 0 , the critical fields of annihilation processes in the micromagnetic model reveal similar values (1.82 T) compared to those from the atomistic model (1.86 T).This suggests that both models produce comparable energy barriers for the magnetic switching.It is worth noting that the frequency observed in the micromagnetic model is linearly related to the critical field [Fig.10(e)], in agreement with Figs.4(e) and 5 in the main