Bloch points in nanostrips

Complex magnetic materials hosting topologically non-trivial particle-like objects such as skyrmions are under intensive research and could fundamentally change the way we store and process data. One important class of materials are helimagnetic materials with Dzyaloshinskii-Moriya interaction. Recently, it was demonstrated that thin nanodisks consisting of two layers with opposite chirality can host a single stable Bloch point of two different types at the interface between the layers. Using micromagnetic simulations we show that FeGe nanostrips consisting of two layers with opposite chirality can host multiple coexisting Bloch points in an arbitrary combination of the two different types. We show that the number of Bloch points that can simultaneously coexist depends on the strip geometry and the type of the individual Bloch points. Our simulation results allow us to predict strip geometries suitable for an arbitrary number of Bloch points. We show an example of an 80-Bloch-point configuration verifying the prediction.

Magnetic quasiparticles with non-trivial topology [1], such as vortices and skyrmions, are of great fundamental interest and could play an important role in novel technological applications [2].One of the quasiparticles is the Bloch point: a single point of vanishing magnetisation [3,4].A stable and manipulable Bloch point of two different types was predicted in helimagnetic two-layer nanodisks [5] where the two layers have opposite chirality, and the Bloch point nucleates at the interface between the layers.
In this work, we demonstrate that two-layer nanostrips can host multiple Bloch points in an arbitrary combination of different Bloch-point types.Using finite-difference micromagnetic simulations, we explore the parameter space of strip geometry to understand for which geometry constraints such magnetisation field configurations can be realised.The Bloch point as a (meta)-stable topological excitation and quasiparticle opens up new avenues of fundamental research.In particular, we can now start to investigate individual and collective behaviour of Bloch points towards discovering Bloch-point-based meta-materials.We conclude our study with a demonstration of encoding a 10-byte string using 80 Bloch points: we identify one Bloch-point type with the binary "1" and the other type with the binary "0" to encode and store the equivalent of an 80-bit long sequence.
The concept of a Bloch point in a two-layer system [5] is explained in Fig. 1, where we start from vortex configurations in single-layer materials, and then stack them on top of each other to obtain the Bloch-point configuration of the magnetisation field.Figure 1a shows schematically the four possible vortex configurations we can encounter in a thin layer of ferromagnetic material due to the competition between ferromagnetic exchange and demagnetisation energy.The vortex core (polarisation P ), pointing in the out-of-plane (z) direction, can either be pointing to +z (P = +1) or −z (P = −1) direction, and -independently -the circularity c of the magnetisation around the vortex core can either be clock-wise (c = −1) or counter-clock-wise (c = +1) [6].By adding the Dzyaloshinskii-Moriya (DM) interaction to the system, the relation between polarisation and circularity is fixed through the chirality, i.e. the sign of the DM interaction strength D: for a given D only two of the four vortex realisations are energetically favourable.
Figure 1b shows how a Bloch-point configuration can be realised by stacking vortex configurations with the same circularity and opposite polarisation on top of each other.The Bloch point emerges at the interface between the two vortex cores of opposite polarisation, and the top-and bottom-layer materials must have opposite signs of D to stabilise the Bloch-point configuration.Despite the high exchange energy density at the Bloch point, the configuration is stabilised by the exchange coupling across the comparatively large interface area between the layers that ensures the same circularity in both layers.Two different types of Bloch points can be realised with the magnetisation vectors of the vortex cores pointing either inwards (head-to-head Bloch point, Fig. 1b HH) or outwards (tail-to-tail Bloch point, Fig. 1b TT).
Figure 1c shows the magnetisation vector field for a single head-to-head Bloch point from our finite-difference micromagnetic simulations, and the colour represents the z component of the magnetisation.Figure 1d shows the corresponding plot for a tail-to-tail Bloch point.In the simulation, the Bloch point forms at the centre of eight discretisation cells as shown in the insets of Fig. 1c and d where the magnetisation of the discretisation cells is visualised by the cones.The Bloch point is located at the intersection of the three isosurfaces m x,y,z = 0 where the magnetisation vanishes.
The paper is organised as follows.First, we discuss two different configurations containing two Bloch points each in Sec.I A. We distinguish between configurations of two same-type Bloch points and two opposite-type Bloch points.These two pairs of Bloch points exhibit all features of a configuration containing multiple Bloch points.In Sec.I B, we demonstrate that multiple Bloch points can coexist in rectangular two-layer nanostrips using micromagnetic simulations.We find that all possible sequences of head-to-head and tail-to-tail Bloch points can be realised.Different combinations have different energy densities depending on the number of neighbouring Bloch points of the same type.We focus on systems containing up to eight Bloch points initially.Based on our results, we can predict a suitable strip geometry for an arbitrary number of Bloch points.To verify this prediction we present one example configuration containing 80 Bloch points in Sec.I C.

A. Two Bloch points
A pair of neighbouring Bloch points in multi-Bloch-point configurations can occur in two fundamentally different combinations.The Bloch points can either be of the same type, for example: a head-to-head (HH) Bloch point next to another HH Bloch point (HH-HH) as shown in the right column of Fig. 2. Alternatively, they can be of opposite type, for example a HH Bloch point next to a tail-to-tail (TT) Bloch point (HH-TT) as shown in the left column of Fig. 2.
The topmost row (Fig. 2a, b) shows a schematic drawing of the nanostrip highlighting the two-layer structure and the geometry.Additionally, position and type of the Bloch points visible in the simulations are indicated with arrows, where the colour of the arrows encodes the z component of the magnetisation (red: +z, blue: −z).
Figure 2c and d show 3D renderings of the simulation results.The isosurfaces show m z = ±0.9(m is the normalised magnetisation), colour indicates the z component.The isosurfaces above/below the Bloch point have a paraboloidallike shape pointing towards the Bloch point (similar to the single-Bloch-point simulation results in Fig. 1c, d).The Bloch point itself is not directly visible in this visualisation.The configuration in Fig. 2d additionally contains one antivortex between the two Bloch points.The m z = 0.9 isosurface of the antivortex extends throughout the whole thickness (z direction) of the two-layer system.The antivortex core shrinks towards the top sample boundary.
To reveal the full three-dimensional structure of the magnetisation field surrounding the Bloch points the magnetisation of each configuration is plotted in five different cut planes for each column (as indicated in the schematic drawings Fig. 2a, b).Four different cut planes show the magnetisation in the xy plane (Fig. 2e -l), at the top sample boundary (z = 10 nm) in Fig. 2e and f, above the interface (z = 1 nm) in Fig. 2g and h, below the interface (z = −1 nm) in Fig. 2i and j, and at the bottom sample boundary (z = −20 nm) in Fig. 2k and l.Colour encodes the z component of the magnetisation vector field, arrows the in-plane component.
Figure 2m and n show the z component of the magnetisation in an xz cut plane going through the Bloch points at y = 50 nm.Magnified subplots show the full magnetisation around the Bloch-point positions and in the centre region between the two Bloch points.The colour of the cones in the magnified areas encodes the y component of the magnetisation vector field.
The results shown in Fig. 2 show that Bloch points form at x ≈ 100 nm and x ≈ 300 nm in both cases, i.e. in both columns.Bloch-point pairs of the opposite type (Fig. 2, left column) show opposite circularity of the magnetisation within the xy plane around the Bloch-point cores.In this case, the in-plane magnetisation (xy component) between the two Bloch points (from x ≈ 100 nm to x ≈ 300 nm) shows a smooth transition from one Bloch point to the other.In contrast, an additional antivortex forms between two same-type Bloch points (Fig. 2, right column) at x ≈ 200 nm to mediate between the incompatible magnetisation configurations that originate from the Bloch points.Differing from the magnetisation of the Bloch-point cores, the magnetisation of the antivortex core (at x ≈ 200 nm) does not change significantly along the z direction (inset in Fig. 2n).

B. Parameter-space diagram and energy density
The spatially averaged energy density of a Bloch-point configuration depends on the number of Bloch points, their individual types, and the strip geometry.Furthermore, different spatial arrangements can be realised, e.g.four Bloch points on a line, or in the corners of a rectangle or diamond shape.Here, we only consider magnetisation configurations containing between one and eight Bloch points in a row, distributed in x direction (strip length) and centred in y direction (strip width).
In Fig. 2 we have seen the two fundamentally different configurations containing two Bloch points, head-to-head and tail-to-tail (HH-TT), and head-to-head and head-to-head (HH-HH).Now we investigate a system containing three Bloch points.In total, eight configurations can be realised.Three configurations are fundamentally different, namely HH-HH-HH, HH-HH-TT, and HH-TT-HH, because they contain distinct numbers of additional antivortices.The other five configurations are equivalent either because HH and TT swap roles (e.g.TT-TT-TT), or because of the symmetry of the system geometry along the x axis (e.g.TT-HH-HH).
The fundamentally different configurations (HH-HH-HH, HH-HH-TT, HH-TT-HH) are shown in Fig. 3d in a system with strip length l = 400 nm.We find one and two additional antivortices (AVs) for configurations HH-HH-TT and HH-HH-HH, respectively.The table in Fig. 3 lists all eight configurations and the respective number of antivortices.
For each of the three fundamentally different configurations, we compute the spatially averaged energy density (E/V ) and plot the three values in Fig. 3a (at l = 400 nm).We find that the HH-TT-HH configuration has the lowest energy density and the HH-HH-HH configuration the highest energy density.We will later discuss in more detail that the presence of antivortices between Bloch points generally increases the energy density of the system.The alternating configuration HH-TT-HH does not contain any antivortices, as these are only needed to mediate between neighbouring Bloch points of the same type.
The three yellow lines (filled and open diamonds) in Fig. 3a show the spatially averaged energy density for the three different configurations as a function of strip length.Not all configurations are stable for all strip lengths: for example the HH-HH-TT configuration is only stable for l ≥ 300 nm.If we try to create the three-Bloch-point configuration HH-HH-TT in a shorter nanostrip, e.g. at l = 275 nm, then the configuration is not stable and will change into a lower-energy configuration, in this case the HH-TT configuration containing only two Bloch points.We can see that the energy generally increases with increasing number of antivortices as mentioned in the previous paragraph.However, there is a deviation visible for l ≤ 225 nm where the HH-HH-HH configuration has a lower energy density than the HH-TT-HH configuration.This deviation is a result of the short strip length near the stability limit.We exclude these regions near the minimal stability strip length in the rest of our discussion.
The blue filled and open squares in Fig. 3a show the energy density for a system containing only two Bloch points.The corresponding magnetisation field for l = 400 nm is visualised in Fig. 3c, and in more detail in Fig. 2. The green circles in Fig. 3a show the energy density for a single Bloch point, and its magnetisation configuration for l = 400 nm is shown in Fig. 3b.For a given strip length l we describe the configuration with the lowest energy density as the energetically most favourable configuration: below l = 250 nm a single Bloch point (green circles) has the lowest energy density.Two opposite-type Bloch points (blue squares) have the lowest energy density for 250 nm < l <= 400 nm and three Bloch points of alternating opposite type (yellow diamonds) have the lower energy density above l = 400 nm.
Figure 3e shows a representation of the energy densities for all possible three-Bloch-point configurations at l = 600 nm.As already discussed, there are three fundamentally different configurations characterised by the number of additional antivortices contained in the configuration (as show in the table in Fig. 3).Different realisations of the same configuration type (swapping HH and TT or using the strip symmetry) exhibit the same energy density.
In Fig. 3a we have seen that the number of Bloch points in the energetically most favourable configuration changes depending on the strip length l. Figure 4 contains a parameter-space diagram showing the energetically most favourable configuration as a function of the strip length and strip width, using the Bloch point number as a label.To create Fig. 4, we ask for each strip length l and a given strip width w, which configuration has the lowest energy density.For example: all data points in Fig. 3a are for a width of w = 100 nm.Close to l ≈ 400 nm, we see that for l ≤ 400 nm the two-Bloch-point configuration HH-TT (blue squares) has the lowest energy density but that for l > 400 nm the three-Bloch-point configuration HH-TT-HH (yellow diamonds) has the lowest value.Figure 4 shows (for w = 100 nm on the y axis) that the two-Bloch-point configuration has the lowest energy density up to l ≈ 400 nm, and the three-Bloch-point configuration for larger l (up to l ≈ 600 nm).All configurations with lowest energy density are of the alternating Bloch-point type, i.e. left and right neighbours of a HH Bloch point are always of type TT, and vice versa (see discussion below).
Figure 4 shows that with increasing strip length, the number of Bloch points that are present in the lowest-energydensity configuration (as shown Fig. 3a for l ≤ 600 nm) increases: for nanostrips with lengths above l ≈ 1300 nm and width w = 100 nm, we find eight Bloch points.Furthermore, Fig. 4 shows that increasing the width of the nanostrip leads to a reduced number of Bloch points in the energetically most favourable configuration.
Figure 4 also shows magnetisation profiles for selected configurations revealing the similarities in the magnetisation profile in different strip geometries.The isosurfaces show m z = ±0.9,colour indicates the z component.All configurations shown in Fig. 4 contain Bloch points of alternating opposite type, i.e. all the lowest-energy configurations do not contain antivortices.
In the discussion of the fundamentally different configurations containing three Bloch points, we have noted that the different configurations can be characterised by the number of additional antivortices contained in the structure.Figure 3f summarises similar findings for eight Bloch points where configurations can contain between zero antivortices (Bloch points of alternating opposite type) and seven antivortices (all Bloch points of the same type).In total, 256 configurations can be realised.The plot in Fig. 3f shows the data for all configurations and a linear fit to the data.Different realisations with the same number of antivortices cannot be distinguished in this plot as their energies are nearly identical but cause broadening of the marker symbols for intermediate numbers of antivortices.The energy density increases linearly with the number of antivortices.
There is an important difference between three Bloch points (Fig. 3e) and eight Bloch points (Fig. 3f).For a fixed number of antivortices, all different three-Bloch-point configurations are equivalent because of the system's symmetry and therefore must have the same energy density.However, for eight Bloch points we additionally find that different configurations that are not related via symmetry also have the same energy density if they contain the same number of antivortices.For example, the configurations HH-HH-HH-HH-HH-HH-HH-TT and HH-HH-HH-HH-TT-TT-TT-TT both contain 6 antivortices (located between neighbouring same-type Bloch points) but cannot be transformed into each other using the symmetries we discussed.Yet, they exhibit the same spatially averaged energy density.Our findings suggest that the energy density of the Bloch points is independent of the configuration of other Bloch points in the system.The energy density can be obtained from a configuration containing Bloch points of alternating opposite type with additional contributions originating from the additional antivortices between neighbouring same-type Bloch points.This is the reason for the lowest-energy-density configurations shown in Fig. 4 consisting of pairs of Bloch points of alternating type: for same-type neighbours an antivortex is required to mediate the magnetisation between the Bloch points of the same type, and the presence of such antivortices would increase the spatially averaged energy density.
We can make one additional observation in Fig. 3a.The energy density of a configuration changes as a function of strip length l.All configurations containing two or three Bloch points have one energy minimum at a certain length that we call the optimal length l o .For example, the optimal length for the HH-TT configuration (blue filled squares in Fig. 3a) is l o ≈ 275 nm.

C. Predicting strip geometries for larger systems
So far, we have focused on small systems containing at most eight Bloch points.Based on this information we can predict strip geometries suitable for an arbitrary number of Bloch points.
Figure 3a shows that meta-stable configurations containing multiple Bloch points can be realised over a broad range of strip lengths but need a certain minimal strip length.This minimal length depends on the number of Bloch points and additional antivortices in the configuration.Furthermore, Fig. 3a shows that all configurations have a minimum in the energy density at a certain optimal length l o .
To predict strip geometries suitable for an arbitrary number of Bloch points we focus on configurations containing up to eight Bloch points of alternating opposite type.We find that the optimal length l o increases linearly with the number of Bloch points (Fig. 5c) with the slope defining the optimal Bloch point spacing s o .Furthermore, we find that s o increases linearly with increasing strip width (Fig. 5d).We obtain s o,w=100 nm ≈ 165 nm and s o,w=200 nm ≈ 272 nm with an estimated accuracy of δs o ≈ 3 nm.These observations lead to our working hypothesis that the ideal Blochpoint spacing s o is independent of the number of Bloch points in the system and suitable to predict geometries for more than eight Bloch points.This prediction can be used for arbitrary configurations not only alternating opposite-type Bloch points.
As an illustrative example, we simulate one specific configuration containing 80 Bloch points, encoding the 10character word Blochpoint in ASCII code (eight bits per letter).We simulate a strip with the predicted length l = 80s o = 13.2 µm at a width of w = 100 nm and with s o = 165 nm.
We minimise the energy of a suitable initial configuration resulting in the 80-bit configuration as shown in Fig. 5a, b.The cross sections show the xy plane at z = 1 nm (Fig. 5a) and the xz plane at y = 50 nm (Fig. 5b).Note that the aspect ratio is not correct in order to improve visibility.Figure 5e shows contour lines for m z for a part of the nanostrip (correct aspect ratio) as indicated in Fig. 5a.Bloch points in Fig. 5e are located at the small red and blue dots.The larger red circles (m z > 0.5) show antivortices between same-type Bloch points.
To test the stability of the 80-Bloch-point configuration we apply a short magnetic field pulse in the +y direction (H = 25 mT/µ 0 , applied for t = 0.5 ns).The modified magnetisation field configuration at the end of the 0.5 ns period is shown in Fig. 5f.Then, we set the applied field back to zero, and let the system evolve freely by carrying out a time-integration.We find that the magnetisation converges back to the initial state: Fig. 5g shows the configuration after t = 5 ns of free relaxation.
To understand the robustness of the predicted geometry, we vary the strip length l and find that the desired 80-Bloch-point configuration can be stabilised over a range of strip geometries.The minimal strip length is around 0.66l o the maximal strip length around 4l o .
Within the range of stability of the 80-bit configuration (0.66l o ≤ l ≤ 4l o ), we find that the length l o is closer to the lower stability boundary (≈ 0.66l o ) than to the upper limit (≈ 4l o ).This is consistent with the energy density curve for the HH-TT configuration in Fig. 3a (blue filled squares) where we see that the energy density as a function of the strip length is asymmetric, and that its energy minimum, located at strip length l o , is located at a comparatively small strip length within the range of possible strip lengths over which the configuration is meta-stable (stability limits are not visible in Fig. 3).

II. DISCUSSION
The Bloch-point configuration originates from vortices with identical circularity, but opposite polarisation, which are stabilised through the DMI of the material, which fixes the core orientation relative to circularity through the leftor right-handed chirality.The Bloch points form an interesting topological excitation in a helimagnetic system, which  extends the set of well-known magnetic structures such as domain walls, vortices, and skyrmions.In the geometry described here, the Bloch points are in equilibrium and can be manipulated.We have found remarkable features of multiple interacting Bloch points in two-layer nanostrips.The two different types -head-to-head (HH) and tail-to-tail (TT) -can be geometrically arranged in any arbitrary order, and this magnetisation configuration resembles a meta-stable configuration (within certain constraints on the strip width and length).The spatially averaged energy density for a system with n Bloch points increases in fixed steps.The number of steps scales -within our accuracy -exactly linearly with the number of antivortices in the configuration (or equivalently: the number of neighbouring same-type Bloch points).We can determine an optimal Bloch-point spacing s o between Bloch points within a line of Bloch points (corresponding to a distance over which a Bloch point extends).
In the following, we speculate about possible future applications of Bloch points.One key-feature distinguishing Bloch points from many other particle-like magnetic configurations is the demonstrated coexistence of Bloch points of two different types in a single sample making Bloch points an interesting candidate for binary data representation.In the racetrack-like designs [7,8], when realised with magnetic excitations of which only one type exists -such as skyrmions -we need to ensure that skyrmions keep their relative positions to be able to interpret the presence of a skyrmion as 1 and the absence of a skyrmion as 0. The two different types of Bloch points presented here could be used to encode binary data without the need to rely on fixed spacing of magnetic objects: a HH configuration could represent "1" and a TT configuration could represent "0".In the context of skyrmion-based realisation of the racetrack approach, other ideas to overcome the fixed-spacing requirement include the use of a combination of skyrmion tubes and chiral bobbers [9] and the two-lane racetrack memory [10].
On the way towards possible applications of Bloch points many more questions need to be addressed.These are related to Bloch point manipulation (movement, switching, creation/annihilation), sensing of Bloch points, and to the thermal stability of Bloch points in general and energy barriers between different configurations containing multiple Bloch points.
In summary, we have demonstrated that two-layer FeGe nanostrips can host multiple Bloch points in any combination of head-to-head and tail-to-tail.Based on our simulations containing up to eight Bloch points, we can predict strip geometries suitable for an arbitrary number of Bloch points.We have verified this prediction by studying a system containing 80 Bloch points.
All results obtained in this work can be reproduced from the repository in Ref. [11] which contains Jupyter notebooks [12] to rerun the micromagnetic simulations and recreate all plots.In the repository pre-computed datasets are also available.

A. System
We simulate rectangular two-layer nanostrips with opposite chirality (opposite sign of D) in the two layers.We vary strip length and width, the thickness of both layers is fixed (bottom layer: 20 nm, top layer: 10 nm).We focus on up to eight Bloch points and accordingly choose nanostrips with lengths between 100 nm and 1400 nm, and widths between 100 nm and 200 nm.The energy equation contains exchange energy density w ex , bulk Dzyaloshinskii-Moriya energy density w dmi , and demagnetisation energy density w d .The magnetisation dynamics is simulated using the Landau-Lifshitz-Gilbert equation [13,14]: where γ * = γ(1 + α 2 ), with γ being the gyromagnetic ratio and α Gilbert damping.Material parameters are based on FeGe [15]: A = 8.87 pJ m −1 , D = 1.58 mJ m −2 , M s = 384 kA m −1 , α = 0.28.We use finite-difference micromagnetic simulations to minimise the energy.All simulations are done using Ubermag [16][17][18] with OOMMF [19] as computational backend and an extension for DMI of crystalclass T [20,21].

B. Simulation procedure
All simulations in this study follow a three-step initialisation and minimisation scheme: (i) initialisation, (ii) fixed minimisation, (iii) free minimisation.In the micromagnetic framework the system is studied at zero temperature, i.e. without thermal fluctuations.Therefore, it is only possible to find local minima that are accessible from the initial configuration.Starting from experimentally feasible initial configurations, such as full saturation, we are able to find magnetisation configurations containing a single or multiple Bloch points depending on the strip geometry.
To facilitate the process of studying arbitrary Bloch-point configurations, independent of the strip geometry in a systematic way, we have developed a simulation scheme that guarantees a magnetisation configuration containing a predictable number of Bloch points.We note that this scheme can probably not be applied directly to an experimental set-up.
For the initialisation, step (i), we start by dividing the nanostrip into equally sized regions (in x direction), one region per Bloch point.To enforce the formation of a Bloch point, the magnetisation in each region is initialised as follows: for a head-to-head Bloch point we initialise the centre region of the topmost layer of cells with m = (0, 0, −1) and all other cells with m = (0, 0, 1).A region hosting a tail-to-tail Bloch point is initialised with reversed z component of the magnetisation (see supplementary Fig. 1 for a schematic plot of the different subregions).We then minimise the energy in two steps (supplementary Fig. 2).During the first energy minimisation, step (ii), we keep the magnetisation of the topmost cells -initialised with reversed magnetisation -and a similarly-sized layer of cells at the bottom sample boundary fixed.This ensures the formation of a Bloch point at the interface between the two layers.The second energy minimisation, step (iii), is done without any fixed cells, i.e. magnetisation in all cells can freely change, and Bloch points could move in any direction to further minimise the energy of the configuration.In this step, the system can find the local energy minimum.

C. Classification
In the micromagnetic framework, it is not possible to directly observe Bloch points because of the fixed norm of the magnetisation vector.A single Bloch point is characterised by the integral value of the topological charge density over a closed surface A surrounding the Bloch point [22]: where F is the emergent magnetic field [23,24].The components of F are defined as: where (i, j, k) is an even permutation of (x, y, z).To detect a single Bloch point in a sample the integral can be computed over the whole sample surface and the exact position of the Bloch point does not need to be known.This method is not directly applicable to multiple Bloch points when their positions are unknown: the sign of the topological charge of a Bloch point depends on its type (HH: S = −1, TT: S = +1).Therefore, contributions to the surface integral from Bloch points of opposite type cancel out. Figure 6 shows the divergence of the emergent field ∇ • F for a HH and a TT Bloch point (a) and two HH Bloch points (b), the two configurations discussed in Fig. 2.
To classify nanostrips that potentially contain multiple Bloch points we compute the convolution of the divergence of the emergent magnetic field with a Heaviside step function Θ: Due to numerical inaccuracies the result of the integral deviates from integer values.By translating the surface integral into a volume integral over the divergence of the emergent magnetic field using the divergence theorem the accuracy can be improved by roughly one order of magnitude.
In our set-up Bloch points are expected to be distributed along x following the strip geometry which justifies computing S as a function of x.This convolution can be interpreted as computing a series of integrals over increasing subvolumes V of the nanostrip starting at the left boundary (x = 0 nm).We round S(x) to integer values and count steps ∆S in this function.

VI. COMPETING INTERESTS
The authors declare no competing interests.

I. INITIALISATION AND ENERGY MINIMISATION SCHEME
In the main text we have explained that we divide the nanostrip into multiple subregions to initialise the system and enforce the formation of a specific number of Bloch points during the energy minimisation.Here, we show two different examples for configurations containing one and two Bloch points, respectively.
Figure 1 shows the eight subregions that are used to obtain a configuration containing two Bloch points.Four subregions are located below z = 0 and have a negative DM energy constant D, and four above z = 0 with positive D. The four small subregions, shown with solid lines, are located at the top and bottom sample boundary.Each of the small subregions is contained within one surrounding subregions shown with dashed lines.The strip geometry in Fig. 1 is stretched in the z direction for better visibility.
FIG. 1. Subregions used for the initialisation and fixed energy minimisation, and to define the two-layer geometry with opposite chirality (opposite sign of D) in the two layers for a nanostrip that shall contain two Bloch points.The magnetisation inside the subregions shown with solid lines is kept fixed during the first energy minimisation.Magnetisation in the top fixed subregions is initialised with reversed z component of the magnetisation (see Fig. 2).In the plot the z axis is stretched for better visibility.The fixed subregions are located at the strip boundary.
As outlined in the main text, the system initialisation and energy minimisation is done in three steps.The magnetisation after each step is show in Fig. 2 for a configuration containing a single HH Bloch point.The system geometry is 100 nm × 100 nm × (10 + 20) nm.This system only needs four subregions, i.e. one "column" of subregions in Fig. 1 (e.g.only the subregions for x < 125 nm).
During the initialisation step the magnetisation is initialised uniformly throughout the sample with m = (0, 0, 1).The only exception is the small subregion at the top sample boundary where the z component of the magnetisation is reversed, m = (0, 0, −1) (Fig. 2a).During the first energy minimisation the magnetisation in the small subregions at the sample boundary, shown with solid lines in Fig. 1 and highlighted in Fig. 2b, is fixed.This enforces the formation of the Bloch point in a controlled manner.The first energy minimisation leads to the magnetisation configuration shown in Fig. 2b.During the second energy minimisation, the magnetisation in all cells can freely change to allow the system find the local energy minimum, leading to the magnetisation shown in Fig. 2c.

II. CLASSIFICATION ACCURACY
In the methods section of the main text we have introduced our classification method.In short, to classify nanostrips that contain multiple Bloch points we compute the convolution of the divergence of the emergent magnetic field F with a Heaviside step function Θ: The components of F are defined as: where (i, j, k) is an even permutation of (x, y, z).The integral generally has non-integer values due to numerical inaccuracies resulting mainly from the finite discretisation cell size.To simplify the classification and counting we round S(x) to integer values whereby we obtain sharp steps (∆S = ±1) at Bloch-point positions.
To justify the rounding to integer values we compare the deviations of S(x) from integer values for different cell sizes.Figure 3a shows a configuration containing eight Bloch points in the pattern TT-HH-TT-TT-TT-TT-HH-TT.We compute S(x) for three different cell sizes with cubic cells with edge lengths l c = 5 nm (Fig. 3b), l c = 2.5 nm (Fig. 3c), and l c = 1 nm (Fig. 3d).For each cell size we show the integral result (solid lines) and rounded values (dashed lines).We can see that the difference (i.e. the inaccuracy) significantly decreases with decreasing cell size.For l c = 1 nm integral and rounded values cannot be distinguished visually.All simulations in the paper are done with a cell edge length l c = 2.5 nm.With decreasing cellsize -cell edge lengths: 5 nm (b), 2.5 nm (c), and 1 nm (d) -the difference between the integral S(x) and the corresponding integer values decreases.

FIG. 1 .
FIG. 1.(a) In a single layer of magnetic material four different vortices, with polarisation P = ±1 and circularity c = ±1, can form as a consequence of the competition between exchange energy and demagnetisation.Adding DMI couples circularity and polarisation.(b) By stacking two layers with opposite sign of the DM energy constant D, a Bloch point can be stabilised.The Bloch point can be of type head-to-head (left) or tail-to-tail (right).In the figure, the two layers are, for better clarity, separated in z direction as indicated by the grey dashed lines.(c, d) Simulation result for a single head-to-head (c) and tail-to-tail (d) Bloch point.The isosurfaces (of paraboloidal-like shape) near the centre show mz = ±0.9,colour indicates the z component.They are convenient to locate the Bloch point that is situated between them.The insets show three isosurfaces for mx = 0, my = 0, and mz = 0, respectively.The Bloch point is located at the intersection of the three isosurfaces where the magnetisation vanishes.The cones in the insets indicate the magnetisation directions in the eight discretisation cells surrounding the Bloch point.

FIG. 2 .
FIG. 2. Magnetisation profile of the two fundamentally different configurations containing two Bloch points: opposite-type Bloch points (head-to-head and tail-to-tail, left column) and same-type Bloch points (head-to-head and head-to-head, right column).The 3D renderings in (c) and (d) show isosurfaces for mz = ±0.9,colour indicates the z component.Several different cut planes in xy and xz are shown to reveal the three-dimensional structure of the Bloch points forming at the interface (z = 0 nm).For the xz plane (subfigures m and n) we also show enlarged plots around the Bloch point and antivortex position.The cones in (m) and (n) are coloured according to their my component, as indicated by the small colour bar in the right bottom corner of the figure.

FIG. 3 .
FIG. 3. (a) Energy densities for energetically different Bloch-point configurations containing at most three Bloch points for different strip lengths l at a strip width w = 100 nm.Simulations have been performed in steps of δl = 25 nm, the solid lines are shown to guide the eye.(b -d) Magnetisation profiles for the six different configurations shown in (a) at l = 400 nm.Isosurfaces show mz = ±0.9,colour indicates the z-component.(e) Energy density for all possible configurations containing three Bloch points.The first three configurations are show in (d) as indicated with the distinct marker symbols (at a different strip length).The table in (a) lists all different configuration containing three Bloch points, highlighting the number and position of the additional antivortices in the different configurations.(f) Energy densities for all configurations containing eight Bloch points.Energy densities for a fixed number of antivortices are nearly identical but cause some "smearing" of the marker symbols.

FIG. 4 .
FIG. 4. Parameter-space diagram showing the energetically most favourable Bloch point number as a function of length l and width w.All configurations contain Bloch points of alternating opposite type.Magnetisation profiles for selected configurations reveal the similarity of the different configurations, isosurfaces show mz = ±0.9.

FIG. 5 .
FIG.5.(a, b) ASCII encoding of the string Blochpoint using 80 Bloch points.Cross sections show (a) the xy plane at z = 1 nm and (b) the xz plane at y = 50 nm.The strip length is chosen according to the predicted value for a nanostrip with width w = 100 nm.Labels on the x axis mark blocks of eight Bloch points, i.e. individual bytes.Note that the aspect ratio is not correct in order to improve visibility.(c) Optimal lengths lo for two to eight Bloch points.The fit is used to predict lengths for more than eight Bloch points.(d) Optimal Bloch point distance so as a function of strip width w. (e -g) show an enlarged part of the nanostrip (correct aspect ratio) as highlighted in (a) to demonstrate the stability of the configuration: (e) initial configuration after energy minimisation; (f) an external magnetic field H = 0.25 T/µ0 is applied in the +y direction for 0.5 ns; (g) after removing the external field the system evolves freely and converges back to the initial state (snapshot after 5 ns).(e -g) show contour lines of the mz component to improve visibility of the disturbance introduced by the external magnetic field.Bloch points are located at the small red and blue dots, the larger red circles show the additional antivortices.

FIG. 6 .
FIG.6.Classification of the two configurations containing two Bloch points shown in Fig.2.The divergence of the emergent magnetic field for the two opposite-and same-type Bloch points is shown in panels (a) and (b), respectively.The xy plane visualised here is located at z = 1 nm, just above the interface.(c, d) The result of the convolution(5) which is used to identify the occurrence of Bloch points and their type due to the steps ∆S = ±1.
Figure 6c and d show S(x) for the two example configurations.A head-to-head Bloch point is identified by ∆S = −1, a tail-to-tail Bloch point by ∆S = +1 corresponding to the topological charge of a Bloch point being S = ±1.Rounding to integer values is justified because deviations from integer values in the integral are a direct V. AUTHOR CONTRIBUTIONS M.L., M.B. and H.F. conceived the study, and M.L. performed finite-difference micromagnetic simulations and analysed the data.M.L., M.B., and H.F. developed the Ubermag software package to drive the finite-difference solver OOMMF.M.L., with the assistance of H.F., M.B., and O.H. prepared the manuscript.

FIG. 2 .
FIG. 2. Initialisation and energy minimsation for a single head-to-head Bloch point (cross-section at y = 50 nm).(a) The magnetisation is initialised with m = (0, 0, 1) with the exeption of a top subregion where m = (0, 0, −1).(b) During the first energy minimisation the magnetisation is kept fixed inside the yellow-highlighted subregions (see Fig. 1 for a 3D plot).The formation of a Bloch point is enforced by the opposite mz values in the two fixed subregions.(c) During the second energy minimisation magnetisation in all cells can freely change.

FIG. 3 .
FIG. 3. Accuracy of the classification for a configuration containing eight Bloch points.(a) shows a cross-section of the magnetisation of configuration TT-HH-TT-TT-TT-TT-HH-TT in the xz plane located at the Bloch-point position (y = 50 nm).With decreasing cellsize -cell edge lengths: 5 nm (b), 2.5 nm (c), and 1 nm (d) -the difference between the integral S(x) and the corresponding integer values decreases.