High level domain specific languages for the finite element method underpin high productivity programming environments for simulations based on partial differential equations (PDE)  while employing automatic code generation to achieve high performance. However, a limitation of this approach is that it does not support operators that are not directly expressible in the vector calculus. This is critical in  applications where PDEs are not enough to accurately describe the physical problem of interest. The use of deep learning techniques have become increasingly popular in filling this knowledge gap, for example to include features not represented in the differential equations, or closures for unresolved spatiotemporal scales. We introduce an interface within the Firedrake finite element system that enables a seamless interface with deep learning models. This new feature composes with the automatic differentiation capabilities of Firedrake, enabling the automated solution of inverse problems. Our implementation interfaces with PyTorch and can be extended to other machine learning libraries. The resulting framework supports complex models coupling PDEs and deep learning whilst maintaining separation of concerns between application scientists and software experts.

§ INTRODUCTION

The design of efficient and composable software relying on high-level languages and automatic differentiation (AD) tools is a major and rapidly growing aspect of modern scientific computing and is of great interest for the machine learning (ML) community. The development of such software requires a vast range of knowledge spanning several disciplines, ranging from applications expertise to mathematical analysis to high performance computing and low-level code optimisation. Software projects relying on automatic code generation have grown in prominence, as their design enables a separation of concerns, increases productivity, and facilitates the collaboration between scientists with different specialisations. Examples of such projects include: Firedrake <cit.>, FEniCS <cit.> and FreeFEM++ <cit.>, in the domain of finite element methods (FEM); and the ML frameworks PyTorch <cit.>, TensorFlow <cit.>, Theano <cit.> and MXNet <cit.>. The Unified Form Language (UFL) <cit.> is a domain specific language (DSL) embedded in Python for the finite element method forming the mathematical programming interface of the Firedrake and FEniCS finite element systems. UFL equips Firedrake and FEniCS with a highly expressive interface to specify the variational forms of PDEs and discrete function spaces, providing the abstractions needed for code generation. Finite element problems are solved using gradient-based methods to minimise the problem residual. Consequently, UFL is a fully differentiable language, with inbuilt automatic differentiation.

Domain-specific languages for PDEs are by definition very specific to partial differential equations. The models of real phenomena of interest to scientists and engineers are, regrettably, seldom so straightforward. Alongside the fundamental physical laws expressed as PDEs are various empirical parametrisations, closures, and regularisation terms which represent aspects of the system for which a more fundamental model is either not known or is practically infeasible for some reason.
This scenario is ubiquitous in science and engineering, ranging from geoscience <cit.> to structural mechanics <cit.> to name but two fields. In order to extend the high productivity, high performance capabilities of UFL and Firedrake to these scenarios, we introduce a foreign function interface to UFL in the form of external operators, by which we mean any operator which is not directly expressible in the vector calculus notation of the existing UFL language.

Machine learning (ML) is a natural tool for creating empirical model components from observed data, and provides the motivating example for the external operator concept. ML also demonstrates the criticality of a differentiable programming approach, since backpropagating the neural network involves the differentiation of both the neural net itself and the PDE to which it is coupled. The work presented here allows deep learning frameworks to be coupled directly into Firedrake as external operators, creating common environment for developing PDE models with deep learning components. The external operator feature composes seamlessly with the dolfin-adjoint library <cit.> enabling automatic differentiation. It interfaces with the PyTorch library and its automatic differentiation engine (torch.autograd), and interfaces to other ML frameworks would be straightforward.

§ THE UNIFIED FORM LANGUAGE (UFL)

The Unified Form Language (UFL) <cit.> is an embedded domain-specific language in Python which provides symbolic representations of finite element simulations. UFL forms are compiled by a domain-specific compiler, which takes the high-level description of the weak form of PDEs provided by UFL and translates this representation into low-level code that assembles the sparse matrices and vectors of the finite element problem.

§.§ Forms

UFL is organised around representing finite element variational forms, and in particular multilinear forms. A multi-linear form (or linear k-forms) is a map from the product of a given sequence $\left\lbrace V_{j} \right\rbrace_{j=0}^{k-1}$ of function spaces to a space $K$:

\begin{equation}
\label{def:multilinear_form}
V_{k-1} \times V_{k-2} \times \cdots \times V_{0} \rightarrow K
\end{equation}

that is linear in each argument. The arity of a form $k$ is the number of argument spaces. Variational forms with arity $k = 0,\ 1,\ $ and $2$ are respectively named functionals, linear forms and bilinear forms. These can be assembled to produce a scalar, a vector and a matrix, respectively. The UFL variational forms can be parametrised by coefficient functions. In this case, the form is expressed as a mapping from a product of a sequence $\left\lbrace W_{j} \right\rbrace_{j=1}^{n}$ of coefficient spaces and the argument spaces:

\begin{equation}
\label{def:multilinear_form}
W_{1} \times W_{2} \times \cdots \times W_{n} \times V_{k-1} \times V_{k-2} \times \cdots \times V_{0} \rightarrow K.
\end{equation}

We refer to (<ref>) as a linear $k$-form with $n$ coefficients. While the multilinear form is necessarily linear in each argument, it can be nonlinear in each coefficient function. An argument is an unknown function in a finite element space, while a coefficient is a known function in a finite element space.

The multilinear forms represented in UFL constitute the weak form of PDEs and comprise a sum of integrals over subspaces of the problem domain.

Let $V$ be a suitable function space and $f$ a given function in $V$. The modified Helmholtz problem is given by: find $u \in V$ such that:

\begin{equation}
\label{example:var_form}
\int_{\Omega} u\, v + \nabla u \cdot \nabla v\ \diff x = \int_{\Omega} f\, v \diff x \quad \forall v \in V.
\end{equation}

The left-hand side of (<ref>) is a 2-form, with $u$ and $v$ as the arguments, and the right-hand side is a 1-form with one coefficient $f$ and one argument $v$.

The symbolic representation of forms in UFL is sufficient to generate the code for matrix and vector assembly, however the execution of that code also requires the simulation data: mesh topology and coordinates, and the value of each coefficient at each finite element node. This is provided by finite element frameworks such as Firedrake and FEniCS, which subclass UFL objects and attach problem data to them, as well as orchestrating the (parallel) assembly and solution operations using the assembly code generated from UFL.

listingFiredrake code defining the variational forms in (<ref>) using UFL. Note the similarity between lines 9 and 10, and (<ref>).

from firedrake import *
mesh = |UnitSquareMesh|(10, 10)
V = |FunctionSpace|(mesh, "Lagrange", 1)
u = |TrialFunction|(V)
v = |TestFunction|(V)

f = |Function|(V)

a = (u * v + |inner|(|grad|(u), |grad|(v))) * dx
L = f * v * dx

§.§ The limited scope of UFL

As a domain specific language, UFL is specialised to the particular application domain it was built for: the symbolic representation of finite element forms. The vector calculus syntax of UFL can represent the weak form of essentially any partial differential equation. However, realistic applications often augment PDEs with terms not readily expressible in vector calculus. In particular, where the physical basis for part of a model is not well understood it may be advantageous to represent this using a neural network trained on suitable data. These terms could not hitherto be represented in UFL.

§ EXTERNAL OPERATORS: OPENING A CLOSED LANGUAGE

We present an expressive, flexible and powerful interface for incorporating arbitrary operators in UFL, and providing their implementation to Firedrake. This is achieved by defining symbolic external operator objects at the UFL level, and equipping them with the appropriate mathematical properties, such as the ability to be differentiated, so that they can be incorporated directly in UFL expressions, and hence be incorporated directly into the operator assembly and problem solution operations provided by Firedrake.

Let $(V_{i})_{1\le i\le l}$, $(W_{i})_{1\le i\le k}$ and $X$ be finite element spaces. An external operator $N$ mapping $k$ operands to $X$ is defined as

\begin{equation}
\label{def_external_operator_OPERATOR}
\begin{aligned}
N \vcentcolon  W_1 \times \cdots \times W_k \times V_l \times \cdots \times V_{1} &\longmapsto  X \\
u_1,\ \ldots,\ u_k,\ v_l,\ \ldots,\ v_1\  \quad &\longrightarrow N(u_1, \ldots, u_k; v_l, \ldots, v_1)
\end{aligned}
\end{equation}

or equivalently,

\begin{equation}
\label{def_external_operator_FORM}
\begin{aligned}
N \vcentcolon  W_1 \times \cdots \times W_k \times  V_l \times \cdots V_1 \times X^{*} &\longmapsto \reels \\
u_1,\ \ldots,\ u_k,\ v_l,\ \ldots,\ v_1,\ v^{*} \quad &\longrightarrow N(u_1, \ldots, u_k; v_l, \ldots, v_1, v^{*})
\end{aligned}
\end{equation}

where $X^{*}$ is the dual space to $X$. That is, the space of bounded linear functionals over $X$. $N$ can be nonlinear with respect to its operands $u_1, \ldots, u_k$ but it is linear with respect to its arguments $v_l, \ldots, v_1, v^{*}$. The numerical evaluation of $N$ is not specified in the UFL language but left to the specific implementation of the external operator $N$.

The equivalence of  (<ref>) and (<ref>) follows from the reflexivity of $X$. That is,
the canonical injection from $X$ into $X^{**}$ is surjective. As a consequence, we can identify $X$ and $X^{**}$, the space of bounded linear functionals on $X^{*}$. In other words, $X \equiv X^{*} \longmapsto \reels$.

Equation (<ref>) suggests that external operators behave like operators between finite element spaces whereas definition (<ref>) suggests that they act as forms in the sense that they are multilinear and scalar-valued. There is no contradiction, both points of view are complementary: they are two different ways to see external operator that are reflected in the implementation. Using (<ref>), an external operator can be used in a form anywhere where a coefficient or argument could be used, while using (<ref>) it is possible to have whole form terms which are simply external operators.

We define a symbolic UFL object representing external operators, the ExternalOperator, whose implementation is left to be specified in Firedrake by the user. In other words, the user can build their own external operator by subclassing the external operator class and providing an implementation of the evaluation of the operator and its derivatives. In the case of neural networks, that can be achieved by calling the evaluation and backpropagation of the neural network using the machine learning framework considered by the user (e.g. PyTorch or TensorFlow). Since Firedrake is embedded in Python, the coupling to these frameworks is straightforward. The $m$ hyperparameters of the neural net can be represented in UFL as a coefficient in the very simple finite element space $\reels^m$, which is constant in space over the problem domain.

§ DIFFERENTIATION

Solving a partial differential equation and training a neural network with a PDE constraint are problems that require the evaluation of the gradient of a form. In UFL, the derivative of a form is based on the Gâteaux derivative. Because external operators can be used in variational forms, we need to extend UFL automatic differentiation rules to handle external operators.

Let $N$ be an external operator of the following form:

\begin{equation}
\label{eq:def_N_star}
\begin{aligned}
N \vcentcolon  V \times V \times  X^{*} &\longmapsto \reels \\
u,\ m, \ v^{*} \quad &\longrightarrow N(u, m; v^{*})
\end{aligned}
\end{equation}

and consider a PDE, defined by the residual form $F$ linear with respect to $v$, such that:

\begin{equation}
\label{residual_form_F}
F(u, m, N(u,m; v^{*});v) = 0 \quad \forall v\in V.
\end{equation}

In order to solve (<ref>) we need to compute the Jacobian of the residual form, i.e. we need to take the Gâteaux derivative in the direction $\hat{u} \in V$:

\begin{equation}
\label{eq:jacobian_F_formula}
\frac{\diff F(u, m, N(u,m; v^{*}); \hat{u}, v)}{\diff u} = \frac{\partial F(u, m, N; \hat{u}, v)}{\partial u} + \frac{\partial F(u, m, N; \widehat{N}, v)}{\partial N} \cdot \frac{\partial N(u, m; \hat{u}, v^{*})}{\partial u}
\end{equation}

where $\hat{N}\in X$ is an argument (unknown function). In UFL, $\frac{\partial N(u, m; \hat{u}, v^{*})}{\partial u}$ is a new external operator representing the Gâteaux derivative of $N$. It is evaluated using the derivative evaluation implementation provided by the user. Listing <ref> shows the UFL code defining and differentiating an external operator and a form containing that operator.

listingCompute $\frac{\partial F}{\partial u}$ and $\frac{\partial N}{\partial u}$ using UFL automatic differentiation

V = |FunctionSpace|(...)
u = |Coefficient|(V)
m = |Coefficient|(V)
v = |TestFunction|(V)
uhat = |TrialFunction|(V)

# |$N : V \times V \times X^{*} \longmapsto \mathbb{R}$|
#      |$u,\ m,\ v^{*} \ \ \ \longrightarrow N(u, m; v^{*})$|
N = |ExternalOperator|(u, m, function_space=X)

# Define a given form F
F = u * N * v * dx

# Symbolically compute the derivative |$\frac{\partial N(u, m; \hat{u}, v^{*})}{\partial u}$|
dNdu = |derivative|(N, u, uhat)

# Symbolically compute the derivative |$\frac{\diff F}{\diff u}$|
dFdu = |derivative|(F, u, uhat)

Automatic differentiation is also employed to compute the gradient in the task of training a neural network constrained by a PDE. In this case, we would have an external operator of the form $N(u, m; v^{*})$ where $u$ would be the input feeding the model and $m$ the parameters. The training would be performed with respect to a cost functional computed from $u$, the solution to the PDE,  (<ref>) in the cost function. The computation of the gradient of the functional with respect to the neural net parameters $m$ requires the solution of the adjoint PDE. This in turn requires taking the adjoint of $\frac{\partial F}{\partial m}$, which by chain rule necessitates taking the adjoint of $\frac{\partial N}{\partial m}$. This last stage is implemented using backpropagation. In fact, evaluating the adjoint of the gradient of the external operator with respect to the parameters will call the evaluation implementation from the external operator subclass representing the model, which itself delegates evaluation to  the machine learning framework in which the model is implemented. At the UFL level, this can be simply achieved via the code in Listing <ref>.

listingAssembling $\frac{\partial N}{\partial m}^{*} \cdot y$ in Firedrake

# Symbolically compute the derivative |$\frac{\partial N}{\partial m}$|
dNdm = |derivative|(N, m)

# Symbolic operation: |$\frac{\partial N}{\partial m}^{*}\cdot y$|
backprop_y = |action|(|adjoint|(dNdm), y)

# Evaluation

§ EXAMPLE: SEISMIC INVERSION

Seismic inversion is the use of seismic reflection data to infer the material properties of the Earth's subsurface. In this process, waves are used to interrogate a medium, and the medium's response is recorded and processed to obtain a description of the subsurface. The problem arises predominantly in exploration geophysics (e.g. oil and gas prospection) <cit.> and geotechnical site characterisation <cit.>. This problem can be formulated as an inverse wave problem. This problem is invariably ill-posed: various solution artifacts result in indistinguishable functionals but nonphysical solutions. A simple model problem can be formulated as follows:

\begin{equation}
\label{regularizer_example}
\min_{c \in P}\ \ \frac{1}{2}\norm{\varphi(c) - \varphi^{obs}}_{V}^{2} + \alpha \mathcal{R}(c)
\end{equation}

where $\varphi^{obs}$ are the observed data, $c$ the scalar wave speed, $\alpha$ the regularisation factor and $\varphi \in V$ is the wave displacement of the medium such that:

\begin{equation}
\label{regularizer_wave_Example}
\begin{aligned}
- \frac{\partial^{2} \varphi}{\partial t^{2}}  - c^{2} \Delta \varphi &= f &\textrm{in}\ \Omega\\
\varphi(t=0) = \varphi_0.
\end{aligned}
\end{equation}

Equation (<ref>) is referred to as the forward problem. In practice, one may use a more complex formulation of this equation to take into account additional physical effects such as elasiticy, acousticity or anisotropy. More realistic applications may also have complex boundary conditions.

Equation (<ref>) is the function to optimise. The first term accounts for the mismatch error between the solution obtained by the forward problem (the PDE) and the observed data. We refer to this term as the fidelity term. $\mathcal{R}$ is a regularisation operator which attempts to counter from noisy observations and to help making the problem well-posed. The PDE is based on fundamental physical laws but the regulariser is more heuristic in nature. A number of recent works have employed deep learning to build regularisers for inverse problems, especially for seismic inversion <cit.> and medical imaging <cit.>. External operators in UFL provide a high-level programming abstraction for this problem. The cost function simply contains an external operator representing the neural network model, i.e. $\mathcal{R}(c) = \frac{1}{2}\norm{N(c)}^{2}$, with $N$ an external operator. Listing <ref> illustrates this in Firedrake.

listingOutline of seismic inversion using a neural net regulariser implemented as an ExternalOperator in Firedrake.

from firedrake import *
from firedrake_adjoint import *

# Get a pre-trained PyTorch model
model = ...
# Define the external operator from the model
pytorch_op = |neuralnet|(model, function_space=...)
N = pytorch_op(vel)

# Solve the forward problem defined by equation |(<ref>)|
phi = F(c)
# Assemble the cost function: |$J = \frac{1}{2} \norm{\varphi(c) - \varphi^{obs}}^{2} + \frac{\alpha}{2} \norm{N(c)}^{2} $|
J = |assemble|(0.5*(|inner|(phi-phi_obs, phi-phi_obs) +
alpha*|inner|(N, N))*dx)

# Optimise the problem
Jhat = |ReducedFunctional|(J, |Control|(c))
c_opt = |minimize|(Jhat, method="L-BFGS-B", tol=1.0e-7,
options="disp": True, "maxiter" : 20)

Recovered wave speed $c$ as a function of position ($x,z$) obtained for the model (<ref>)-(<ref>): exact velocity (upper left), without regularisation (upper right), Tikhonov regulariser (lower left), neural network-based regulariser (lower right).

The exact version of Firedrake (including UFL) used in this paper is archived in <cit.> while the code used to run the example and generate the output figure is archived in <cit.>

§.§ Acknowledgement

This work was funded by a President’s PhD Scholarship at Imperial College London.


# Screw dislocations in cubic chiral magnets

Maria Azhar  Institut für Theoretische Festkörperphysik, Karlsruhe Institute
of Technology, 76131 Karlsruhe, Germany  Volodymyr P. Kravchuk  Institut für
Theoretische Festkörperphysik, Karlsruhe Institute of Technology, 76131
Karlsruhe, Germany  Bogolyubov Institute for Theoretical Physics of National
Academy of Sciences of Ukraine, 03680 Kyiv, Ukraine  Markus Garst  Institut
für Theoretische Festkörperphysik, Karlsruhe Institute of Technology, 76131
Karlsruhe, Germany  Institute for Quantum Materials and Technology, Karlsruhe
Institute of Technology, 76131 Karlsruhe, Germany

###### Abstract

Helimagnets realize an effective lamellar ordering that supports disclination
and dislocation defects. Here, we investigate the micromagnetic structure of
screw dislocation lines in cubic chiral magnets using analytical and numerical
methods. The far field of these dislocations is universal and classified by an
integer strength $\nu$ that characterizes the winding of magnetic moments. We
demonstrate that a rich variety of dislocation-core structures can be realized
even for the same strength $\nu$. In particular, the magnetization at the core
can be either smooth or singular. We present a specific example with $\nu=1$
for which the core is composed of a chain of singular Bloch points. In
general, screw dislocations carry a non-integer but finite skyrmion charge so
that they can be efficiently manipulated by spin currents.

Introduction – The Dyzaloshinskii-Moriya interaction (DMI) in cubic chiral
magnets like MnSi, FeGe, or Cu2OSeO3 favours helimagnetic long-range order in
a large region of the phase diagram [1, 2, 3, 4, 5, 6, 7]. A finite field
$\boldsymbol{H}$ aligns the helix axis and, in addition, tilts the magnetic
moments towards the field direction giving rise to a conical magnetic helix,
see Fig. 1(a). This helimagnetic ordering realizes a one-dimensional periodic
texture that shares many similarities with other emerging lamellar structures
found, e.g., in various soft matter systems [8, 9, 10].

In particular, in the limit of weak spin-orbit coupling (SOC) the phase
transition from the paramagnetic to the helimagnetic phase at
$\boldsymbol{H}=0$ is a fluctuation-driven first-order transition similar to
the ones in certain cholesteric liquid crystals or diblock co-polymers [11,
12, 13]. The correlation length above the critical temperature $T_{c}$
possesses a temperature dependence that is well described by weak
crystallization theory [14]. This indicates that the paramagnetic regime just
above $T_{c}$ is characterized by strong correlations that are maintained by
the large density of states of paramagnons [15]. It is still an important open
issue whether these pronounced magnetic correlations are also at the origin of
the non-Fermi liquid behavior observed in MnSi and FeGe upon suppressing the
critical temperature towards zero with pressure [16, 17, 18, 19, 20].

Figure 1: (a) Conical magnetic helix with pitch $\lambda_{h}$ and cone angle
$\theta_{0}$ enclosed by the magnetic moments and the direction of the applied
magnetic field $\boldsymbol{H}$. (b) Helimagnetic order is characterized by
equidistant isosurfaces, the hallmark of lamellar order, where, e.g., the
x-component of magnetization, $n_{x}$, assumes a particular value. (c)-(e)
Example of a screw dislocation with strength $\nu=1$ illustrated by different
vertical cross-sections. The core of this specific example contains a chain of
magnetic Bloch points.

In this context, the question arises as to whether this intriguing
paramagnetic-helimagnetic phase transition can also be understood from the
dual perspective as a defect-mediated melting transition [21, 22, 23]. A
prerequisite to address this question is an understanding of the elementary
defects of helimagnetic order. It is well-known that defects of lamellar
structures in general consist of disclinations and dislocations, which are
line defects in case the lamellae are embedded in three dimensional space [8].
For helimagnets, such defects were discussed on a phenomenological level by
Kléman [24]. Recently, it was shown both theoretically and experimentally that
domain walls of helimagnetic order might consist of an arrangement of
disclinations and edge dislocations [25, 26, 27] very similar to domain walls
in cholesteric liquid crystals [28]. Moreover, at helimagnetic twist grain
boundaries screw dislocations are expected to occur [8, 29]. It was also
demonstrated that the motion of edge dislocations is an important relaxation
process for disordered helimagnets possibly accounting for the large
relaxation times observed experimentally [30, 31, 32, 33].

However, there exist additional line excitations within the conical helix
phase that are distinct from dislocations and disclinations. In particular,
chiral magnets are famously known to host skyrmions, i.e., topological two-
dimensional magnetic textures [34]. In bulk magnets, the skyrmion textures
extend along the third direction forming skyrmion strings that either condense
into a lattice or exist as metastable excitations of the field-polarized phase
[35, 36, 7], for a recent review see Ref. [37]. It has been demonstrated in
[38] that such metastable skyrmion configurations also persist within the
conical helix phase upon decreasing the magnetic field below the critical
field $H_{c2}$. Such skyrmion strings attract each other and can form clusters
or even networks [39, 40, 41]. Nevertheless, these skyrmion strings within the
conical helix phase possess an exponentially decaying far field and, in
contrast to dislocation lines, are characterized by a vanishing Burgers
vector. In addition, the conical helix can also support localized large-
amplitude excitations like bound pairs of hedghog defects, i.e., Bloch points
[42] and even Hopfions [43].

In the present work, we theoretically investigate in detail screw dislocations
in cubic chiral magnets. The helimagnetic order defines equally spaced
isosurfaces where, e.g., the $x$-component of the magnetization assumes the
same value, see Fig. 1(b). The deviation of isosurfaces from their equilibrium
configuration is described by the displacement field $\boldsymbol{u}$ 111The
displacement field $\bm{u}$ describes small deviations of the fronts of
constant phases $\phi(\bm{r})=\text{const}$ of the helimagnetic order where
$\phi(\bm{r})=2\pi(\bm{r}+\bm{u})\hat{\bm{z}}/\lambda_{h}$.. The integral
along a loop enclosing a dislocation line, $\oint
d\boldsymbol{u}=\boldsymbol{b}$, is finite and given by the Burgers vector
$\boldsymbol{b}$ indicating that $\boldsymbol{u}$ is singular at the
dislocation core. For a screw dislocation,
$\boldsymbol{b}=\lambda_{h}\nu\hat{\bm{z}}$ is aligned with the helix axis
$\hat{\bm{z}}$ and its size is an integer multiple of the helix pitch
$\lambda_{h}$, where $\nu\in\mathbb{Z}\setminus\\{0\\}$ characterizes the
strength of the screw dislocation (sdν).

Using analytical arguments and numerical simulations we determine the
micromagnetic structure of screw dislocations. In the limit of small SOC, when
the influence of magnetocrystalline anisotropies is negligible, we find that
they possess in the far field the expected universal form of lamellar
structures with a displacement vector $\boldsymbol{u}=u_{z}(x,y)\hat{\bm{z}}$
where [8, 9, 10]

$\displaystyle u_{z}(x,y)=\frac{\lambda_{h}}{2\pi}\nu\chi.$ (1)

Here, $\chi$ is the polar angle of cylindrical real-space coordinates
$(\rho,\chi,z)$. Moreover, we show that the magnetization texture at the core
of screw dislocations can either be smooth or might contain Bloch points,
which we illustrate explicitly for screw dislocations with $|\nu|=1$.

Theory of cubic chiral magnets – The magnetic energy functional $E=\int
d\bm{r}\mathcal{E}$ of cubic chiral magnets possesses a density that reads in
leading order in SOC

$\displaystyle\mathcal{E}=A(\partial_{i}\bm{n})^{2}+D\bm{n}(\bm{\nabla}\times\bm{n})-M_{s}\mu_{0}Hn_{z}.$
(2)

Here, $\bm{n}$ is a unit vector specifying the orientation of the local
magnetization and the magnetic field is applied in the $z$-direction,
$\boldsymbol{H}=H\hat{\bm{z}}$; $A$ is the exchange constant, $D>0$ is the DMI
assuming a right-handed chiral magnetic system, $M_{s}$ is the saturation
magnetization and $\mu_{0}$ is the magnetic constant. Importantly, for zero
field $H=0$ this density is isotropic with respect to a combined rotation of
spin and real space. This rotational symmetry is explicitly broken by
magnetocrystalline anisotropies that are however weak in the limit of small
SOC and will be mostly neglected in the following. We also neglect for
simplicity the magnetic dipolar interaction.

It is convenient to consider the representation of
$\bm{n}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$ in terms of polar
angle $\theta$ and azimuthal angle $\phi$. For large fields
$H>H_{c2}=D^{2}/(2A\mu_{0}M_{s})$ the ground state of Eq. (2) is field-
polarized with $\theta=0$. The ground state for small fields $0\leq H<H_{c2}$
is the conical helix with a position dependent polar angle $\phi=2\pi
z/\lambda_{h}$ where the helix pitch $\lambda_{h}=4\pi A/D$, and the cone
angle $\theta=\theta_{0}$ with $\cos\theta_{0}=H/H_{c2}$.

Far field of screw dislocations – Performing an asymptotic analysis of the
Euler-Lagrange equations of Eq. (2), for details see the supplementary
material [45], we find that the conical helix supports screw dislocations with
an asymptotic behavior for large distances from their core $\rho\to\infty$,

$\displaystyle\theta$
$\displaystyle=\theta_{0}+\frac{2\nu\sin^{2}\theta_{0}}{1+\sin^{2}\theta_{0}}\frac{\lambda_{h}}{2\pi\rho}\sin\big{(}(\nu-1)\chi+\frac{2\pi
z}{\lambda_{h}}\big{)}+\mathcal{O}(\rho^{-2}),$ (3) $\displaystyle\phi$
$\displaystyle=\frac{2\pi}{\lambda_{z}}\Big{(}z+u_{z}(x,y)\Big{)}+\mathcal{O}(\rho^{-2}),$
(4)

where $u_{z}$ is the universal displacement field of Eq. (1). The dependence
of $\phi$ indicates that the $(x,y)$-components of the magnetization form a
vortex within each plane perpendicular to the applied field, i.e., for each
value of $z$, and the winding number is just given by the disclination
strength $\nu$. The structure of the vortex changes from plane to plane as a
function of $z$ due to the linear dependence of $\phi$ on $z$, see Fig. 1.

The far-field configuration allows to determine the topological skyrmion
charge for a screw dislocation within each $z$-plane, $N_{\rm top}(z)=\int
dxdy\rho_{\rm top}$ where $\rho_{\rm
top}=\frac{1}{4\pi}\bm{n}(\partial_{x}\bm{n}\times\partial_{y}\bm{n})$. Note
that, in contrast to skyrmions, here $N_{\rm top}(z)$ is not an integer. At
infinity $\rho\to\infty$, the magnetization $\bm{n}$ encircles as a function
of real-space angle $\chi$ $\nu$-times the $\hat{\bm{z}}$ axis with a fixed
polar angle $\theta_{0}$. Assuming that the magnetization texture is smooth in
a given $z$-plane, we obtain for the charge

$\displaystyle N_{\rm top}(z)=\frac{\nu}{2}\left(1-h\right)+n(z),$ (5)

with the reduced field $h=H/H_{c2}$, and $n(z)\in\mathbb{Z}$ is an integer
that depends on the magnetization at the core of the screw dislocation. As we
will see below, in case that $n(z)$ varies with $z$ the core is singular and
contains Bloch points.

Figure 2:  Core of the screw dislocation sd${}^{+}_{-1}$ with strength
$\nu=-1$ obtained by micromagnetic simulations for $H=0$. The magnetic moments
represented by arrows form within a given $z$-plane an antivortex-like
structure. The orange isosurface in the main panel is defined by $n_{z}=1/2$.
The magnetization at the core is preferentially pointing in $\hat{\bm{z}}$
direction; the configuration with opposite core magnetization,
sd${}^{-}_{-1}$, is degenerate at $H=0$ (not shown). Note that the periodicity
of the core structure along $z$ is characterized by a wavelength
$2\lambda_{h}$.

Plugging the asymptotics of Eqs. (3) and (4) into Eq. (2) we obtain for the
energy of a screw dislocation line per length

$\displaystyle\varepsilon_{\rm sd}=\varepsilon^{\rm core}_{\rm
sd}+A(2\pi\nu)^{2}\frac{h^{2}(1-h^{2})}{2-h^{2}}\log\frac{R}{\rho_{\rm
core}},$ (6)

where $h=H/H_{c2}$. The length scale $R$ specifies the extension of the system
in radial direction, and $\rho_{\rm core}$ is the linear size of the
dislocation core with the associated core energy $\varepsilon^{\rm core}_{\rm
sd}$. The far-field tail of the screw dislocation gives rise to a contribution
to the energy that in general diverges logarithmically with the radial system
size $R$. This logarithmic contribution vanishes at the transition $H=H_{c2}$
to the field-polarized phase because the dislocation ceases to be defined when
the cone angle vanishes, $\theta_{0}=0$. Moreover, it also vanishes in zero
field $H=0$ due to the rotational symmetry of Eq. (2). It is well-known that
lamellar structures emerging in an isotropic environment are characterized by
particularly soft small-amplitude, i.e., phonon excitations that possess the
Landau-Peierls form [8]. In this case the contribution to the energy from the
far field of screw dislocations, which can be captured in terms of a static
phonon field, vanishes [8, 9, 10, 46] in agreement with Eq. (6). Technically,
this is here due to a cancellation of exchange and DMI energies. The prefactor
of the logarithmic contribution in Eq. (6) remains finite however at $H=0$ if
magnetocrystalline anisotropies are taken into account that explicitly break
the continuous rotational symmetry of the theory [45].

Core structure of screw dislocations – Having established the far field of
screw dislocations we now turn to the discussion of their core structures. We
employ micromagnetic simulations [45] in order to determine the core and its
energy. First, we focus on the case with strength $\nu=-1$, see Fig. 2. The
magnetization can here be continuously extrapolated from the far field towards
the core resulting in a smooth texture. There exist two energetically
degenerate configurations at zero field where the core magnetization is either
aligned or anti-aligned with the field, respectively denoted by
sd${}^{+}_{-1}$ and sd${}^{-}_{-1}$ in the following. In order to decrease DMI
energy, the core deforms elliptically such that it mimics a small Bloch-domain
wall with preferred chirality. In addition, we found that close to zero field
the ellipse is further deformed into a banana-shape structure, see cross-
sections in Fig. 2, that leads to a periodicity of the core along the $z$-axis
with an enhanced wavelength $2\lambda_{h}$. For finite field $H$, the
configuration with the aligned core magnetization sd${}^{+}_{-1}$ is
energetically favoured, see Fig. 3. The precise $H$-dependence of the
dislocation energy depends on the system size $R$ but it vanishes for $H\to
H_{c2}$.

Figure 3: Energy per length of the various screw dislocation lines shown in
Figs. 2 and 4, as a function of $H$ with a comparison to the skyrmion string
energy (gray lines) whose core magnetization is aligned (Sk+) or anti-aligned
(Sk-) with the field. The energies were obtained by micromagnetic simulations
for a cylinder-shaped system with radius $R=5\lambda_{h}$. Dashed line shows
the logarithmic contribution of the far field in Eq. (6) assuming
$\rho_{\text{core}}=\lambda_{h}/2$ for illustration (for simplicity we assume
that $\rho_{\text{core}}$ does not depend on $H$). The energy of dislocations
sd${}^{+}_{-1}$ and sd${}^{+}_{1}$ vanishes at $H_{c2}$ where they can be
identified as vortex lines of the XY-order parameter. The energy of
sd${}^{{\rm Sk}-}_{1}$ merges at $H_{c2}$ with that of the skyrmion string
Sk-. Points $A$-$D$ mark field values where the corresponding configurations
became unstable in numerical simulations due to the lattice discreteness.
Figure 4:  Core of various screw dislocations with strength $\nu=1$ obtained
by micromagnetic simulations. Orange and blue isosurfaces are respectively
defined by $n_{z}=1/2$ and $n_{z}=-1/2$ except in panel (a) where $n_{z}=1/2$
and $n_{z}=-1/4$. (a) Dislocation sd${}^{+}_{1}$ with core magnetization
aligned with the field for $H=0.21H_{c2}$. (b) Dislocation sd${}^{{\rm
Bp}}_{1}$ for $H=0$ with an alternating core magnetization separated by Bloch
points (yellow spheres) with alternating topological charges $\pm 1$. (c)
Dislocation sd${}^{{\rm Sk}-}_{1}$ at $H=0$ with an anti-aligned core
magnetization; it smoothly connects to a skyrmion string configuration for
$H\to H_{c2}$.

There exists in fact a screw dislocation for each strength $\nu$ whose energy
vanishes at $H_{c2}$. The phase transition at $H_{c2}$ corresponds to a magnon
condensation [47] that is in the XY-universality class. Employing a standard
Holstein-Primakoff expansion around the field-polarized state
$\bm{n}=\hat{\bm{z}}+(\psi\,e^{-i2\pi
z/\lambda_{h}}(\hat{\bm{x}}+i\hat{\bm{y}})+c.c.)+\mathcal{O}(|\psi|^{2})$ for
$H>0$ the complex spin wave function $\psi$ can be identified with the
corresponding XY-order parameter. When it condenses for $H\lesssim H_{c2}$
with non-zero constant $\psi$, long-range conical order emerges. The $U(1)$
symmetry of the complex wave function $\psi$ then also supports vortex line
solutions given by $\psi\sim e^{-i\nu\chi}$ with an amplitude vanishing at the
core. The size of the vortex core is determined by the correlation length of
the condensate, $\xi\sim 1/\sqrt{1-H/H_{c2}}$ [48], and is thus independent of
$z$. This implies a cylindrical core of the vortex line close to $H_{c2}$.
Hence, these vortex lines can be identified with a special type of screw
dislocations of the helimagnetic order that possess a smooth magnetization at
the core that is aligned with the applied field.

Let us consider the corresponding solution with strength $\nu=1$ denoted by
sd${}^{+}_{1}$ in Fig. 4(a). Here, the vortex structure of the
$(x,y)$-components of magnetization in the far field continuously alternate
from a divergenceless to a rotationless configuration as a function of $z$
with a topological charge $N_{\rm top}=\frac{1}{2}(1-H/H_{c2})$ that is
independent of $z$. It is instructive to focus on the planes
$z=\frac{1}{2}(m\pm\frac{1}{2})\lambda_{h}$ in Eq. (4) with $m\in\mathbb{Z}$
where a divergenceless configuration is realized. As the core is approached,
the magnetization smoothly rotates in a right-handed or left-handed manner
until it is aligned with the field at the center, see lower panels in Fig.
4(a). Whereas the former texture is favoured by the DMI, the latter is
disfavoured. As a consequence, the core, that is cylindrical close to
$H_{c2}$, becomes undulated along the $z$-axis for smaller fields with a
contraction on the $z$-planes housing the disfavoured textures. The energetic
cost of the latter also leads to an increase of the dislocation energy for
decreasing $H$. In the simulations we found that this particular screw
dislocation structure cannot be maintained for lowest magnetic fields.
Instead, a first order transition to a different structure occurs at point $B$
in Fig. 3, to which we turn next.

The energy cost of the disfavoured configurations shown in the lower right
panel of Fig. 4(a) can be avoided by switching the core magnetization within
these planes. An alternating core magnetization along the $z$-axis is indeed
characteristic for the screw dislocation sd${}^{\rm Bp}_{1}$ shown in Fig.
4(b). At zero field, it realizes a Bloch-like meron structure on the planes
$z=\frac{1}{2}(m\pm\frac{1}{2})\lambda_{h}$ that are both favoured by the DMI.
These merons possess alternating skyrmion charges $N_{\rm top}=\pm\frac{1}{2}$
that implies the presence of Bloch points with alternating topological charges
$\pm 1$ 222The topological charge $c$ of a Bloch point is defined by
$c=\int_{V}\mathrm{d}\bm{r}\bm{\nabla}\cdot\bm{\Omega}$ where $V$ is a small
volume around the Bloch point and
$\Omega_{i}=\varepsilon_{ijk}\frac{1}{8\pi}\bm{n}(\partial_{j}\bm{n}\times\partial_{k}\bm{n})$.
positioned in the core on the intermediate planes $z=m\lambda_{h}/2$ with
$m\in\mathbb{Z}$. This screw dislocation sd${}^{{\rm Bp}}_{1}$ with a chain of
Bloch points at its core is the most stable configuration for $\nu=1$ close to
zero field but it is energetically more costly than the dislocations
sd${}^{\pm}_{-1}$ with $\nu=-1$. It can be maintained for a large field range
but it becomes unstable in the simulations at point $C$ in Fig. 3 where
oppositely charged Bloch points annihilate before reaching the critical field
$H_{c2}$.

There exist a third screw dislocation with strength $\nu=1$ that again
possesses a smooth core texture without singularities. Its core magnetization
is either fully aligned or anti-aligned with the applied field. Both
configurations, sd${}^{{\rm Sk}+}_{1}$ and sd${}^{{\rm Sk}-}_{1}$,
respectively, are degenerate at zero field, but the anti-aligned core is
energetically favoured at finite $H$, see Fig. 3. In the simulations, the
configuration with the aligned core can only be stabilized for small field
values up to point $D$. In order to elucidate the core structure, we focus in
Fig. 4(c) on sd${}^{{\rm Sk}-}_{1}$ at $H=0$ and consider again the planes
$z=\frac{1}{2}(m\pm\frac{1}{2})\lambda_{h}$ with a divergenceless
configuration of the magnetization in the far field. As the core is approached
within these planes, the magnetization smoothly rotates in a Bloch-like
fashion that is favoured by the DMI. In half of these planes it is sufficient
to rotate the magnetization by $\pi/2$ but in the complementary planes a
rotation by $3\pi/2$ is required in order to reach a uniformly magnetized
core. The topological charge within each plane is given by $N_{\rm
top}=-\frac{1}{2}(1+H/H_{c2})$. As the field increases, the in-plane texture
transforms from a meron with $N_{\rm top}=-\frac{1}{2}$ to a skyrmion with
$N_{\rm top}=-1$. In fact, as the critical field $H_{c2}$ is approached this
screw dislocation smoothly converts into a skyrmion configuration of the
field-polarized state. The screw dislocation sd${}^{{\rm Sk}}_{1}$ close to
$H_{c2}$ can thus be viewed as a bound state of a skyrmion string with a
vortex of the XY-order parameter $\psi$.

Summary – Screw dislocations belong to a diversity of topological textures in
chiral magnets. They are line excitations of the conical helix phase and
possess a universal far field that is characterized by their strength
$\nu\in\mathbb{Z}$. Their micromagnetic core structure can be multifarious.
Focusing on $\nu=\pm 1$ we showed that it can be either smooth or singular
comprising a chain of Bloch points. The core of screw dislocations with larger
strength $|\nu|$ can be even richer that will be shown elsewhere. Moreover,
screw dislocations carry a finite non-integer skyrmion charge $N_{\rm top}$
and can thus be efficiently manipulated by spin currents [37].

M.G. is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) in the framework of SFB 1143 (project A07; project-id 247310070),
project-id 270344603 and SPP 2137 (project-id 324327023).

## References

  * Bak and Jensen [1980] P. Bak and M. H. Jensen, Journal of Physics C: Solid State Physics 13, L881 (1980).
  * Nakanishi  _et al._ [1980] O. Nakanishi, A. Yanase, A. Hasegawa, and M. Kataoka, Solid State Communications 35, 995 (1980).
  * Ishikawa  _et al._ [1976] Y. Ishikawa, K. Tajima, D. Bloch, and M. Roth, Solid State Communications 19, 525 (1976).
  * Lebech  _et al._ [1989] B. Lebech, J. Bernhard, and T. Freltoft, Journal of Physics: Condensed Matter 1, 6105 (1989).
  * Uchida  _et al._ [2006] M. Uchida, Y. Onose, Y. Matsui, and Y. Tokura, Science 311, 359 (2006).
  * Adams  _et al._ [2012] T. Adams, A. Chacon, M. Wagner, A. Bauer, G. Brandl, B. Pedersen, H. Berger, P. Lemmens, and C. Pfleiderer, Physical Review Letters 108, 237204 (2012).
  * Seki  _et al._ [2012] S. Seki, X. Z. Yu, S. Ishiwata, and Y. Tokura, Science 336, 198 (2012).
  * Chaikin and Lubensky [2000] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ , 1st ed. (Cambridge University Press, 2000).
  * de Gennes and Prost [1995] P. G. de Gennes and J. Prost,  _The Physics of Liquid Crystals_ (Oxford University Press, 1995).
  * Kleman and Lavrentovich [2004] M. Kleman and O. D. Lavrentovich, eds.,  _Soft Matter Physics: An Introduction_ (Springer New York, 2004).
  * Janoschek  _et al._ [2013] M. Janoschek, M. Garst, A. Bauer, P. Krautscheid, R. Georgii, P. Böni, and C. Pfleiderer, Physical Review B 87, 134407 (2013).
  * Buhrandt and Fritz [2013] S. Buhrandt and L. Fritz, Physical Review B 88, 195137 (2013).
  * Živković  _et al._ [2014] I. Živković, J. S. White, H. M. Rønnow, K. Prša, and H. Berger, Phys. Rev. B 89, 060401 (2014).
  * Brazovskii  _et al._ [1987] S. A. Brazovskii, I. E. Dzyaloshinskii, and A. R. Muratov, Sov. Phys. JETP , 625 (1987).
  * Kindervater  _et al._ [2019] J. Kindervater, I. Stasinopoulos, A. Bauer, F. X. Haslbeck, F. Rucker, A. Chacon, S. Mühlbauer, C. Franz, M. Garst, D. Grundler, and C. Pfleiderer, Physical Review X 9, 041059 (2019).
  * Pfleiderer  _et al._ [2001] C. Pfleiderer, S. R. Julian, and G. G. Lonzarich, Nature 414, 427 (2001).
  * Pfleiderer  _et al._ [2007] C. Pfleiderer, P. Boni, T. Keller, U. K. Rossler, and A. Rosch, Science 316, 1871 (2007).
  * Pedrazzini  _et al._ [2007] P. Pedrazzini, H. Wilhelm, D. Jaccard, T. Jarlborg, M. Schmidt, M. Hanfland, L. Akselrud, H. Q. Yuan, U. Schwarz, Y. Grin, and F. Steglich, Physical Review Letters 98, 047204 (2007).
  * Ritz  _et al._ [2013] R. Ritz, M. Halder, M. Wagner, C. Franz, A. Bauer, and C. Pfleiderer, Nature 497, 231 (2013).
  * Schmalian and Turlakov [2004] J. Schmalian and M. Turlakov, Physical Review Letters 93, 036405 (2004).
  * Toner [1982] J. Toner, Physical Review B 26, 462 (1982).
  * Grinstein  _et al._ [1986] G. Grinstein, T. C. Lubensky, and J. Toner, Physical Review B 33, 3306 (1986).
  * Zhai and Radzihovsky [2021] Z. Zhai and L. Radzihovsky, Annals of Physics 10.1016/j.aop.2021.168509 (2021).
  * Kléman [1970] M. Kléman, Philosophical Magazine 22, 739 (1970).
  * Li  _et al._ [2012] F. Li, T. Nattermann, and V. L. Pokrovsky, Phys. Rev. Lett. 108, 107203 (2012).
  * Nattermann and Pokrovsky [2018] T. Nattermann and V. L. Pokrovsky, Journal of Experimental and Theoretical Physics 127, 922 (2018).
  * Schoenherr  _et al._ [2018] P. Schoenherr, J. Müller, L. Köhler, A. Rosch, N. Kanazawa, Y. Tokura, M. Garst, and D. Meier, Nature Physics 14, 465 (2018).
  * Bouligand [1983] Y. Bouligand, Dislocations in solids (North-Holland Publishing Company, New York, 1983) Chap. 23.
  * Martin  _et al._ [2017] N. Martin, M. Deutsch, G. Chaboussant, F. Damay, P. Bonville, L. N. Fomicheva, A. V. Tsvyashchenko, U. K. Rössler, and I. Mirebeau, Phys. Rev. B 96, 020413 (2017).
  * Bauer  _et al._ [2017] A. Bauer, A. Chacon, M. Wagner, M. Halder, R. Georgii, A. Rosch, C. Pfleiderer, and M. Garst, Physical Review B 95, 024429 (2017).
  * Milde  _et al._ [2020] P. Milde, L. Köhler, E. Neuber, P. Ritzinger, M. Garst, A. Bauer, C. Pfleiderer, H. Berger, and L. M. Eng, Physical Review B 102, 024426 (2020).
  * Dussaux  _et al._ [2016] A. Dussaux, P. Schoenherr, K. Koumpouras, J. Chico, K. Chang, L. Lorenzelli, N. Kanazawa, Y. Tokura, M. Garst, A. Bergman, C. L. Degen, and D. Meier, Nature Communications 7, 12430 (2016).
  * Schoenherr  _et al._ [2021] P. Schoenherr, M. Stepanova, E. N. Lysne, N. Kanazawa, Y. Tokura, A. Bergman, and D. Meier, (2021), arXiv:2105.00658 [cond-mat.mtrl-sci] .
  * Bogdanov and Hubert [1994] A. Bogdanov and A. Hubert, Journal of Magnetism and Magnetic Materials 138, 255 (1994).
  * Mühlbauer  _et al._ [2009] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Science 323, 915 (2009).
  * Yu  _et al._ [2010] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Nature 465, 901 (2010).
  * Back  _et al._ [2020] C. Back, V. Cros, H. Ebert, K. Everschor-Sitte, A. Fert, M. Garst, T. Ma, S. Mankovsky, T. L. Monchesky, M. Mostovoy, N. Nagaosa, S. S. P. Parkin, C. Pfleiderer, N. Reyren, A. Rosch, Y. Taguchi, Y. Tokura, K. von Bergmann, and J. Zang, Journal of Physics D: Applied Physics 53, 363001 (2020).
  * Leonov  _et al._ [2016] A. O. Leonov, T. L. Monchesky, J. C. Loudon, and A. N. Bogdanov, Journal of Physics: Condensed Matter 28, 35LT01 (2016).
  * Du  _et al._ [2018] H. Du, X. Zhao, F. N. Rybakov, A. B. Borisov, S. Wang, J. Tang, C. Jin, C. Wang, W. Wei, N. S. Kiselev, Y. Zhang, R. Che, S. Blügel, and M. Tian, Phys. Rev. Lett. 120, 197203 (2018).
  * Sohn  _et al._ [2019] H. R. O. Sohn, S. M. Vlasov, V. M. Uzdin, A. O. Leonov, and I. I. Smalyukh, Phys. Rev. B 100, 104401 (2019).
  * Leonov  _et al._ [2021] A. O. Leonov, C. Pappas, and I. I. Smalyukh, (2021), arXiv:2103.12950 .
  * Müller  _et al._ [2020] G. P. Müller, F. N. Rybakov, H. Jónsson, S. Blügel, and N. S. Kiselev, Phys. Rev. B 101, 184405 (2020).
  * Voinescu  _et al._ [2020] R. Voinescu, J.-S. B. Tai, and I. I. Smalyukh, Phys. Rev. Lett. 125, 057201 (2020).
  * Note [1] The displacement field $\bm{u}$ describes small deviations of the fronts of constant phases $\phi(\bm{r})=\text{const}$ of the helimagnetic order where $\phi(\bm{r})=2\pi(\bm{r}+\bm{u})\hat{\bm{z}}/\lambda_{h}$.
  * [45] Supplementary information.
  * Santangelo [2006] C. D. Santangelo, Liquid Crystals Today 15, 11 (2006).
  * Giamarchi  _et al._ [2008] T. Giamarchi, C. Rüegg, and O. Tchernyshyov, Nature Physics 4, 198 (2008).
  * Ivanov and Sheka [1995] B. A. Ivanov and D. D. Sheka, Low Temperature Physics 21, 881 (1995).
  * Note [2] The topological charge $c$ of a Bloch point is defined by $c=\intop\ilimits@_{V}\mathrm{d}\bm{r}\bm{\nabla}\cdot\bm{\Omega}$ where $V$ is a small volume around the Bloch point and $\Omega_{i}=\varepsilon_{ijk}\frac{1}{8\pi}\bm{n}(\partial_{j}\bm{n}\times\partial_{k}\bm{n})$.
  * Vansteenkiste  _et al._ [2014] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, AIP Advances 4, 107133 (2014).
  * [51] The Object Oriented MicroMagnetic Framework, developed by M. J. Donahue and D. Porter mainly, from NIST. We used the 3D version of the 1.2b3 release.
  * Cortés-Ortuño  _et al._ [2018] D. Cortés-Ortuño, M. Beg, V. Nehruji, R. A. Pepper, and H. Fangohr, OOMMF extension: Dzyaloshinskii-Moriya interaction (DMI) for crystallographic classes T and O (2018).
  * Karhu  _et al._ [2012] E. A. Karhu, U. K. Rößler, A. N. Bogdanov, S. Kahwaji, B. J. Kirby, H. Fritzsche, M. D. Robertson, C. F. Majkrzak, and T. L. Monchesky, Physical Review B 85, 094429 (2012).
  * Wilson  _et al._ [2013] M. N. Wilson, E. A. Karhu, D. P. Lake, A. S. Quigley, S. Meynell, A. N. Bogdanov, H. Fritzsche, U. K. Rößler, and T. L. Monchesky, Phys. Rev. B 88, 214420 (2013).
  * Rohart and Thiaville [2013] S. Rohart and A. Thiaville, Phys. Rev. B 88, 184422 (2013).
  * Dzyaloshinskii and Ivanov [1979] I. E. Dzyaloshinskii and B. A. Ivanov, JETP Letters 29, 592 (1979).
  * Kosevich  _et al._ [1990] A. M. Kosevich, B. A. Ivanov, and A. S. Kovalev, Physics Reports 194, 117 (1990).
  * Ackerman and Smalyukh [2017] P. J. Ackerman and I. I. Smalyukh, Physical Review X 7, 011006 (2017).
  * Tai and Smalyukh [2018] J.-S. B. Tai and I. I. Smalyukh, Physical Review Letters 121, 187201 (2018).
  * Sutcliffe [2018] P. Sutcliffe, Journal of Physics A: Mathematical and Theoretical 51, 375401 (2018).
  * Hopf [1931] H. Hopf, Mathematische Annalen 104, 637 (1931).
  * Kravchuk and Sheka [2007] V. Kravchuk and D. Sheka, Physics of the Solid State 49, 1923 (2007).
  * Malozemoff and Slonzewski [1979] A. P. Malozemoff and J. C. Slonzewski,  _Magnetic domain walls in bubble materials_ (Academic Press, New York, 1979).
  * Döring [1968] W. Döring, Journal of Applied Physics 39, 1006 (1968).
  * Pylypovskyi  _et al._ [2012] O. V. Pylypovskyi, D. D. Sheka, and Y. Gaididei, Physical Review B 85, 224401 (2012).
  * Charilaou [2020] M. Charilaou, Physical Review B 102, 014430 (2020).
  * Li  _et al._ [2020] Y. Li, L. Pierobon, M. Charilaou, H.-B. Braun, N. R. Walet, J. F. Löffler, J. J. Miles, and C. Moutafis, Physical Review Research 2, 033006 (2020).

SUPPLEMENTARY MATERIALS

In section I we provide details of the micromagnetic simulations. In section
II we present the derivation of the universal far field configuration of screw
dislocations in chiral magnets and compare with micromagnetic simulations. In
addition, we discuss the influence of rotational symmetry breaking on the
screw dislocation energy at zero field. In section III the micromagnetic core
structures of screw dislocations is further analyzed. In particular, we
demonstrate the use of linking numbers to distinguish between the
configurations sd${}_{1}^{\text{Sk}\pm}$ and sd${}_{1}^{\pm}$, and we further
compare the skyrmion tube Sk+ with the screw dislocation
sd${}_{1}^{\text{Sk}+}$. Finally, we discuss the energy and positions of Bloch
points contained in the screw dislocation sd${}^{\rm Bp}_{1}$.

##  I Micromagnetic simulations of screw dislocations in chiral magnets

For the micromagnetic simulations we used both MuMax [50] and OOMMF [51, 52]
with a discrete cubic lattice with an elementary unit cell with size $\Delta
x=\Delta y=\Delta z=0.5$ nm. We simulated a cylindrical sample with periodic
boundary conditions along the cylinder axis, i.e., $z$-axis that coincides
with the magnetic field direction. On the edge of the cylinder we applied free
boundary conditions. For concreteness, we used parameters appropriate fro
MnSi, i.e., $A=0.4$ pJ/m, $\mu_{0}M_{s}=0.203$ T, $D=0.28$ mJ/m2 [53],
implying $\lambda_{h}\approx 18$ nm and $\mu_{0}H_{c2}=0.6$ T. The size of the
simulated volume was $L_{z}=m\lambda_{h}$ along the $z$-axis with an integer
$m$ between $1<m\leq 4$. The radius of the cylinder was varied between 5 and
11 times the helix wavelength $\lambda_{h}$.

##  II Far field of screw dislocations in chiral magnets

###  II.1 Asymptotic expansion of the Euler-Lagrange equations

The Euler-Lagrange equations for the energy density functional of Eq. (2) in
the main text can be written in the form

$\displaystyle-\bm{\nabla}^{2}\theta+\sin\theta\cos\theta\left[(\bm{\nabla}\phi)^{2}-2\partial_{\zeta}\phi\right]+h\sin\theta-2\sin^{2}\theta(\bm{\varepsilon}\cdot\bm{\nabla}\phi)$
$\displaystyle=0,$ (7a)
$\displaystyle-\bm{\nabla}\cdot\left(\sin^{2}\theta\bm{\nabla}\phi\right)+2\sin^{2}\theta\left(\bm{\nabla}\theta\cdot\bm{\varepsilon}\right)+2\sin\theta\cos\theta\partial_{\zeta}\theta$
$\displaystyle=0.$ (7b)

Here, we used the angular parameterization of the unit vector
$\bm{n}=\sin\theta\bm{\varepsilon}+\cos\theta\hat{\bm{z}}$ with
$\bm{\varepsilon}=\hat{\bm{x}}\cos\phi+\hat{\bm{y}}\sin\phi$, and the
derivatives are taken with respect to the dimensionless coordinate
$2\pi\bm{r}/\lambda_{h}$ with $\zeta=2\pi z/\lambda_{h}$ being the
dimensionless coordinate in direction of the field $h=H/H_{c2}$. Here, we
denote the dimensionless cylindrical coordinates by $\\{\rho,\chi,\zeta\\}$.

Let us first consider the case $h=0$ where the ground state is given by the
helix with $\theta=\pi/2$, $\phi=\zeta$. Substituting into (7) a vortex-like
solution $\theta=\pi/2$, $\phi=\nu\chi+\zeta$, one can see that Eq. (7b) is
fulfilled but Eq. (7a) has a non-vanishing left-hand-side (LHS) of order
$\mathcal{O}(\rho^{-1})$. In order to compensate this contribution, we look
instead for solutions of the form
$\theta=\pi/2+\vartheta_{1}(\chi,\zeta)/\rho$ and $\phi=\nu\chi+\zeta$. In the
limit of large $\rho$ Eq. (7a) reduces in this case to

$\frac{\partial_{\zeta\zeta}^{2}\vartheta_{1}-\vartheta_{1}+2\nu\sin[(\nu-1)\chi+\zeta]}{\rho}+\mathcal{O}\left(\frac{1}{\rho^{3}}\right)=0,$
(8)

and the LHS of Eq. (7b) is only of order $\mathcal{O}\left(\rho^{-2}\right)$.
At this order, the equations are fulfilled by choosing
$\vartheta_{1}=\nu\sin[(\nu-1)\chi+\zeta]$ so that the numerator in (8)
vanishes. The accuracy can be iteratively improved in a next step by looking
for a solution with the form $\theta=\pi/2+\vartheta_{1}/\rho$ and
$\phi=\nu\chi+\zeta+\varphi_{2}(\chi,\zeta)/\rho^{2}$, where $\vartheta_{1}$
was found in the previous step. The LHS of Eq. (7a) is still on the order
$\mathcal{O}\left(\rho^{-3}\right)$ and Eq. (7b) obtains the form

$\frac{\partial_{\zeta\zeta}^{2}\varphi_{2}+\nu(2-\nu)\sin[2(\nu-1)\chi+2\zeta]}{\rho^{2}}+\mathcal{O}\left(\frac{1}{\rho^{4}}\right)=0.$
(9)

The function $\varphi_{2}$ is determined by demanding that this equation
vanishes up to order $\mathcal{O}\left(\frac{1}{\rho^{4}}\right)$ yielding
$\varphi_{2}=\frac{1}{4}\nu(2-\nu)\sin[2(\nu-1)\chi+2\zeta]$.

Repeating this procedure we obtain the asymptotic solution for the screw
dislocation in the form of a series up to any required accuracy. Up to terms
of order $\mathcal{O}\left(\rho^{-5}\right)$ this solution is

$\displaystyle\theta=\frac{\pi}{2}+\frac{a_{11}\sin\psi}{\rho}+\frac{a_{21}\sin\psi+a_{22}\sin
3\psi}{\rho^{3}}+\mathcal{O}\left(\frac{1}{\rho^{5}}\right)$ (10a)
$\displaystyle\phi=\nu\chi+\zeta+\frac{b_{11}\sin
2\psi}{\rho^{2}}+\frac{b_{21}\sin 2\psi+b_{22}\sin
4\psi}{\rho^{4}}+\mathcal{O}\left(\frac{1}{\rho^{6}}\right),$ (10b)

where $\psi=(\nu-1)\chi+\zeta$ and coefficients $a_{ij}$ and $b_{ij}$ are
given by

$\begin{split}&a_{11}=\nu,\quad a_{21}=\frac{3}{8}\nu^{2}(2-\nu),\quad
a_{22}=-\nu\left(\frac{1}{5}-\frac{\nu}{4}+\frac{\nu^{2}}{24}\right);\\\
&b_{11}=\frac{\nu(2-\nu)}{4},\quad
b_{21}=\nu^{2}\left(\frac{8}{5}-\frac{5\nu}{4}+\frac{\nu^{2}}{8}\right),\quad
b_{22}=-\nu(2-\nu)\left(\frac{3}{80}-\frac{\nu}{16}+\frac{\nu^{2}}{32}\right).\end{split}$
(11)

Thus, one can conjecture that the deviation of the screw dislocation of the
helimagnetic order from the vortex-like solution can be expressed in odd
powers of $1/\rho$ for $\theta$ and in even powers for $\phi$.

However, this symmetry is broken by the magnetic field when the helix is
becoming conical. Applying the procedure described above for the case $h\neq
0$ we obtain

$\displaystyle\theta=$
$\displaystyle\theta_{0}+\frac{2\nu\sin^{2}\theta_{0}}{1+\sin^{2}\theta_{0}}\frac{\sin\psi}{\rho}-\frac{\nu^{2}h}{\sqrt{1-h^{2}}}\frac{\cos^{4}\theta_{0}-2\sin^{4}\theta_{0}\sin^{2}\psi}{\rho^{2}\left(1+\sin^{2}\theta_{0}\right)^{2}}+\mathcal{O}\left(\frac{1}{\rho^{3}}\right),$
(12a) $\displaystyle\phi=$
$\displaystyle\nu\chi+\zeta+\frac{\nu(2-\nu)}{2}\frac{\sin^{2}\theta_{0}}{1+\sin^{2}\theta_{0}}\frac{\sin
2\psi}{\rho^{2}}-\frac{4h\nu^{2}}{\sqrt{1-h^{2}}}\frac{[1-(1+\nu)\sin^{2}\theta_{0}]\cos\psi}{\rho^{3}(1+\sin^{2}\theta_{0})^{2}}+\mathcal{O}\left(\frac{1}{\rho^{4}}\right),$
(12b)

where $\cos\theta_{0}=h$. Based on (12) we conclude that the asymptotic
expansion (12) is valid for $\rho\gg\max\\{1,|\nu|h/\sqrt{1-h^{2}}\\}$. As the
critical field is approached $h\to 1^{-}$, the size of the core of the screw
dislocation is governed by the correlation length $\xi\sim 1/\sqrt{1-h}$ that
diverges for $h\to 1^{-}$, see discussion in the main text.

Generally, the asymptotic solution is obtained with the help of the Ansatz

$\theta=\theta_{0}+\sum\limits_{m=1}^{\infty}\frac{\vartheta_{m}(\chi,\zeta)}{\rho^{m}},\quad\phi=\nu\chi+\zeta+\sum\limits_{m=1}^{\infty}\frac{\varphi_{m}(\chi,\zeta)}{\rho^{m}}.$
(13)

Substituting (13) into (7) and collecting terms at the same powers $\rho^{-m}$
one obtains a cascade of sets of equations with two equations in each set.
Solving the first set ($m=1$) results in $\vartheta_{1}$ and $\varphi_{1}$.
Solving the second set ($m=2$) using the solution for $\vartheta_{1}$ and
$\varphi_{1}$ results in $\vartheta_{2}$ and $\varphi_{2}$. The iterations are
repeated until the required accuracy is achieved. In our particular case, it
turns out that the each set consists of two independent equations (this is
true at least for $m\leq 4$).

###  II.2 Comparison of the far field with micromagnetic simulations

Figure S1: Far field of three distinct screw dislocations, sd${}_{-1}^{+}$,
sd${}_{1}^{\rm Bp}$ and sd${}_{1}^{+}$. The angles
$\vartheta=\theta-\theta_{0}$ and $\varphi=\phi-\nu\chi-\zeta$ extracted from
micromagnetic simulations (solid line) are compared with the asymptotic
expansion of Eq. (12) (dashed line). For all cases $h=0.5$, $\chi=\pi/4$, and
the dimensionless coordinate $\zeta=2.7$, 1.7 and 4.1 for sd${}_{-1}^{+}$,
sd${}_{1}^{\text{Bp}}$ and sd${}_{1}^{+}$, respectively. The simulations were
carried out for a cylinder-shaped sample with height $\lambda_{h}$ and radius
$11.1\lambda_{h}$. The dimensionless radial coordinate is thus bounded between
$0<\rho<2\pi 11.1\approx 69.7$. The strong deviations close to the edge of the
cylinder are due to surface twist effects. Good agreement between numerics and
analytics is obtained in the intermediate range where Eq. (12) is applicable,
see also Fig. S2.  Figure S2: Further comparison of the far field between
numerics and analytics for screw dislocations sd${}_{1}^{\rm Bp}$ and
sd${}_{-1}^{+}$. The quantity $\vartheta\rho$ obtained from the micromagnetic
simulations (solid lines) is compared to its analytically expected asymptotic
value $2\nu\sin^{2}\theta_{0}[1+\sin^{2}\theta_{0}]^{-1}\sin\psi$ (dashed
lines) for the case $h=0.5$, $\chi=\pi/4$ for different values of $\zeta$. The
green shaded area indicates the region of the dislocation core. The surface
effect is quite pronounced hampering a quantitative comparison with the
analytical asymptotic solution at larger distances $\rho$.

The comparison of the asymptotic solutions (12) in the presence of a finite
magnetic field with the micromagnetic simulations are presented in Fig. S1.
Good agreement is found in the far field except close to the sample edge where
the DMI-induced surface twist [54, 55] leads to significant distortions not
accounted for by the analytical treatment. This is further illustrated in Fig.
S2 that shows $\vartheta\rho=(\theta-\theta_{0})\rho$ which is expected to
approach a constant value in the far field (dashed line). The approach to this
value is hampered by the surface twist that appears to possess a large
penetration depth in the present case; this requires further study. Fig. S2
also allows to estimate the size of the screw dislocation core that is
indicated by the green shaded region. For the chosen parameters in this
example, the core is on the order of a single helix wavelength $\lambda_{h}$,
i.e., $\rho\sim 2\pi$.

###  II.3 Influence of uniaxial anisotropy on the screw dislocation energy

The energy per length of a screw dislocation line in general depends
logarithmically on the system size, see Eq. (6) of the main text. This
logarithmic dependence disappears in zero field for the theory of Eq. (2) of
the main text where the theory possesses rotational symmetry. It was argued in
the main text that a logarithmic dependence remains even at $H=0$ when
magnetocrystalline anisotropies are taken into account that explicitly break
the rotation symmetry.

In order to demonstrate this point, we consider in the following a modified
Langrange density that contains an additional uniaxial anisotropy $K$,

$\mathcal{E}=A(\partial_{i}\bm{n})^{2}+D\bm{n}(\bm{\nabla}\times\bm{n})-M_{s}\mu_{0}Hn_{z}+Kn_{z}^{2}.$
(14)

Whereas an uniaxial anisotropy does not comply with the symmetry of cubic
chiral magnets, it is easier to treat as compared to cubic magnetocrystalline
anisotropies and it serves our purpose.

One of the Euler-Lagrange equations still coincides with Eq. (7b), whereas the
other in dimensionless units is now given by

$\begin{split}-\bm{\nabla}^{2}\theta+\sin\theta\cos\theta[(\bm{\nabla}\phi)^{2}&-2\partial_{\zeta}\phi-\kappa]+h\sin\theta-2\sin^{2}\theta(\bm{\varepsilon}\cdot\bm{\nabla}\phi)=0,\end{split}$
(15)

where $\kappa=4KA/D^{2}$ is the dimensionless anisotropy coefficient. The
conical state solution $\phi=\zeta$ and $\theta=\theta_{0}$ with
$\cos\theta_{0}=h/(1+\kappa)$ exists if $\kappa>-1$ and $|h|<1+\kappa$. For
the far field of the screw dislocation we now obtain

$\displaystyle\theta=\theta_{0}+\frac{2\nu\sin^{2}\theta_{0}}{1+(1+\kappa)\sin^{2}\theta_{0}}\frac{\sin\psi}{\rho}+\mathcal{O}(\rho^{-2}),$
(16a)
$\displaystyle\phi=\nu\chi+\zeta+\frac{\nu(2-\nu)\sin^{2}\theta_{0}}{2[1+(1+\kappa)\sin^{2}\theta_{0}]}\frac{\sin
2\psi}{\rho^{2}}+\mathcal{O}(\rho^{-3}).$ (16b)

The corresponding energy per length is

$\varepsilon_{\rm sd}=\varepsilon^{\rm core}_{\rm
sd}+A(2\pi\nu)^{2}\sin^{2}\theta_{0}\frac{1-(1-\kappa)\sin^{2}\theta_{0}}{1+(1+\kappa)\sin^{2}\theta_{0}}\log\frac{R}{\rho_{\rm
core}},$ (17)

where the field dependence is incorporated into $\theta_{0}$. For the case
$\kappa=0$ one has $\sin^{2}\theta_{0}=1-h^{2}$, and the result of Eq. (6) in
the main text is reproduced. For $h=0$, i.e., $\sin^{2}\theta_{0}=1$ we obtain
instead

$\varepsilon_{\rm sd}=\varepsilon^{\rm core}_{\rm
sd}+A(2\pi\nu)^{2}\frac{\kappa}{2+\kappa}\log\frac{R}{\rho_{\rm core}}.$ (18)

Thus, for a non-vanishing uniaxial anisotropy $\kappa$ the screw dislocation
energy is logarithmically divergent even at zero magnetic field.

##  III Further discussions on the magnetic microstructure of screw
dislocations

###  III.1 Linking numbers for screw dislocations sd${}_{1}^{\text{Sk}+}$ and
sd${}_{1}^{+}$

Figure S3: Field evolution of the skyrmion number calculated for an arbitrary
horizontal plane $z=\text{const}$ for screw dislocations
sd${}_{1}^{\text{Sk}+}$ and sd${}_{1}^{+}$ as well as for the skyrmion tube
Sk+ as obtained from micromagnetic simulations. This explicitly confirms that
the skyrmion number of both dislocations is $N_{\rm
top}=\frac{1}{2}(1-H/H_{c2})$ whereas $N_{\rm top}=1$ for Sk+. Note that due
to the slow algebraic decay of the dislocations’ far field, an accurate
determination of $N_{\text{top}}$ requires significantly larger system size
within the $(x,y)$-plane as compared to the exponentially localized skyrmion
tubes. Alternatively, for the particular case of $\nu=1$ one can calculate
$N_{\text{top}}$ for planes $\zeta=n\pi$ (corresponding to
$z=\pi\lambda_{h}/2$), where the leading corrections to the asymptotic
solutions (12) vanish.  Figure S4: Isosurfaces of screw dislocations
sd${}_{1}^{\text{Sk}+}$ (panel (a) and (b)) and sd${}_{1}^{+}$ (panel (c) and
(d)) as obtained by micromagnetic simulations. The magnetic field is
$H=0.21H_{c2}$ corresponding to points $B$ and $D$ in Fig. 3 of the main text.
The coloring represents the angle $\phi=\arg(n_{x}+in_{y})$. Panels (b) and
(d) correspond to vertical cross-sections of (a) and (c), respectively, and
also show isolines on the isosurfaces characterized by constants $\phi=0$ and
$\phi=\pi$ corresponding to the dark-red and dark-cyan lines, respectively.

Whereas the stability regime of the screw dislocations sd${}_{1}^{\text{Sk}+}$
and sd${}_{1}^{+}$ as a function of magnetic field is quite different, see
Fig. 3 of the main text, they nevertheless share common features. They both
possess the same strength $\nu=1$, a smooth core, their core magnetization is
aligned with the $z$-axis, and the skyrmion number in each plane is for both
given by $N_{\rm top}=\frac{1}{2}(1-H/H_{c2})$, see Fig. S3 (such that
sd${}_{1}^{\text{Sk}+}$ smoothly connects to a skyrmion string with $N_{\rm
top}=1$ at the negative critical field $H=-H_{c2}$). Coincidentally, both
become also unstable at a similar magnetic field and transform into sd${}^{\rm
Bp}_{1}$ by the nucleation of singular Bloch points, see points $B$ and $D$ in
Fig. 3 of the main text.

Nevertheless, both screw dislocations sd${}_{1}^{\text{Sk}+}$ and
sd${}_{1}^{+}$ can be distinguished topologically by linking numbers. Both
dislocations possess toroidal isosurfaces with a negative magnetization
surrounding the core line, cf. Fig.4(a) and (c) of the main text. This suggest
a certain similarity with Hopfions [56, 57, 58, 59, 60]. However in contrast
to Hopfions, the screw dislocations sd${}_{1}^{\text{Sk}+}$ and sd${}_{1}^{+}$
are extended line defects. They are not spatially localized and, as a
consequence, they can not be characterized by a Hopf invariant [61]. Note that
Fig. 4 of the main text shows only a single period along the $z$-axis each
containing a toroidal isosurface. These isosurfaces are again shown in Fig.
S4(a) and (c) with a coloring that is reflecting the phase
$\phi=\arg(n_{x}+in_{y})$.

The nature of the toroidal isosurfaces in both cases is quite distinct. This
is best illustrated by considering the plane $z=\frac{1}{4}\lambda_{h}$ shown
in the lower panels of Fig. 4(a) and (c). When approaching the core from the
far field, the magnetization first twist right-handedly for sd${}_{1}^{+}$ but
then reverses and twists left-handedly in order to align the magnetization
with the $z$-axis at the center. In total the magnetization rotates by $\pi/2$
but due to the twist reversal a toroidal isosurface is created. In contrast,
the magnetization always twists in a right-handed fashion upon approaching the
core for sd${}_{1}^{\text{Sk}+}$ with a total rotation angle of $3\pi/2$. The
additional rotation by $\pi$ results in the toroidal isosurface for
sd${}_{1}^{\text{Sk}+}$.

This distinct nature of toroidal isosurfaces is reflected in the topological
properties of isolines where the magnetization points in a specific direction.
Examples of isolines on the toroidal isosurfaces are shown in Fig. S4(b) and
(d). Whereas pairs of isolines on the toroidal isosurfaces of
sd${}_{1}^{\text{Sk}+}$ are interlinked and thus characterized by a finite
linking number $N_{\rm link}=1$, the corresponding pairs for sd${}_{1}^{+}$
are not linked and thus $N_{\rm link}=0$. At the same time, pairs of isolines
defined on isosurfaces close to the core are linked for sd${}_{1}^{+}$ with
$N_{\rm link}=1$ and not linked for sd${}_{1}^{\text{Sk}+}$ with $N_{\rm
link}=0$. This particular situation, that different pairs of isolines possess
different linking numbers for the same screw dislocation, is here possible
because these line defects are extended and not spatially localized.

Figure S5: Comparison between the skyrmion tube Sk+ in panel I and the screw
dislocation sd${}_{1}^{\text{Sk}+}$ in panel II at zero field. All data were
obtained by micromagnetic simulations. Each panel (a) illustrates the
micromagnetic structure within a single period $\lambda_{h}$ along the
magnetic field axis where orange and blue isosurfaces are defined by
$n_{z}=1/2$ and $-1/2$, respectively. Panels (b) show cross-sections of the
magnetization. Panel (c) displays the distribution of the gyrofield
$\bm{\Omega}=\hat{\bm{x}}_{i}\frac{\epsilon_{ijk}}{8\pi}\bm{n}\cdot\left[\partial_{j}\bm{n}\times\partial_{k}\bm{n}\right]$;
here, the magenta surface is defined by
$|\bm{\Omega}|=0.5\Omega_{\text{max}}$. Panel (d) shows cylindrical components
of $\bm{\Omega}$ in units of $\lambda_{h}^{-2}$.

###  III.2 Comparison of the screw dislocation sd${}_{1}^{\text{Sk}+}$ and the
skyrmion tube Sk+ embedded in the conical helix

The aim of this section is to compare the core structures of the screw
dislocation sd${}_{1}^{\text{Sk}+}$ and the skyrmion tube Sk+. Furthermore, we
will illustrate the merger of both as the negative critical field $H\to-
H_{c2}$ is approached from above.

The key difference between the two configurations is the different behavior of
their far field. In the limit $\rho\to\infty$ far away from the core, the
magnetization within each $z$-plane is uniform and vortex-like for the
skyrmion tube and the screw dislocation sd${}_{1}^{\text{Sk}+}$, respectively.
This is illustrated in Fig. S5. As the far field of the skyrmion tube Sk+ is
uniform the topological skyrmion number evaluated for each $z$-plane is an
integer $N_{\text{top}}=1$ that is independent of the applied magnetic field,
see Fig. S3. For details on the field evolution of the skyrmion tube within
the conical state see Ref. [41].

In contrast, in zero magnetic field $H=0$ the screw dislocation
sd${}_{1}^{\text{Sk}+}$ forms a meron within each $z$-plane with the skyrmion
number $N_{\rm top}=1/2$. This is consistent with the formula $N_{\rm
top}=p\nu/2$ for vortices in two-dimensional ferromagnets [57] where $\nu$ is
the strength and $p=n_{z}(\rho=0)$ is the polarity, here $p=1$. In the
presence of a finite magnetic field, the mapping of the order parameter from
any $z$-plane covers more or less than half of the Bloch sphere for negative
and positive $H$, respectively. This is reflected by the field-dependent
skyrmion number, $N_{\rm top}=(1-H/H_{c2})/2$ for the screw dislocation, see
Fig. S3.

Figure S6: Evolution of $\text{sd}_{1}^{\text{Sk}+}$ and Sk+ in a magnetic
field. As the negative critical field is approach $h\to-1$ the two merge into
the same skyrmion tube within the field-polarized state. Data is obtained by
means of micromagnetic simulations for a cylinder-shaped sample with radius
$R_{0}=5\lambda_{h}$. Color of the isosurfaces represents the value of the
$z$-component $n_{z}$, see legend. In the insets, the positions on the axes
are given in units of $\lambda_{h}$.

As the negative critical field is approached $H\to-H_{c2}$ the screw
dislocation sd${}_{1}^{\text{Sk}+}$ and the skyrmion tube Sk+ merge into the
same structure, i.e., a skyrmion tube within the field-polarized state, see
the illustration in Fig. S6. As discussed in the main text,
sd${}_{1}^{\text{Sk}+}$ can be considered close to $-H_{c2}$ as a bound state
of a skyrmion tube with a vortex of the XY-order parameter of the magnon
condensation transition at $H_{c2}$. The latter endows sd${}_{1}^{\text{Sk}+}$
with the vortex-like far field.

The smooth evolution of sd${}_{1}^{\text{Sk}+}$ as a function of field $H\to-
H_{c2}$ or, equivalently, of sd${}_{1}^{\text{Sk}-}$ for $H\to H_{c2}$ is to
be contrasted with sd${}^{+}_{-1}$ that switches its core magnetization before
reaching the critical field $H_{c2}$, see Fig. 3 of the main text. The latter
behavior is similar to vortices in two-dimensional ferromagnets that can also
show a switching of their core polarity [62].

###  III.3 Screw dislocation sd${}^{\rm Bp}_{1}$ and the energy of its Bloch
points

The aim of this section is to provide some quantitative explanations about the
equilibrium positions of the Bloch points (BP) along the core of
sd${}_{1}^{\text{Bp}}$ at zero field. Moreover, we discuss the question as to
why antivortex-like dislocations sd${}_{-1}^{\pm}$ do not possess BPs.

Let us consider a BP $\bm{n}_{\textsc{bp}}=c\,\mathfrak{R}\,\bm{r}/|\bm{r}|$
of charge $c$ with $\mathfrak{R}$ being the matrix of spatial rotations [63].
The Bloch point charge $c=\pm 1$ is defined by the relation [63]

$\bm{\nabla}\cdot\bm{\Omega}=c\delta(\bm{r}-\bm{r}_{\textsc{bp}}),\qquad\bm{\Omega}=\hat{\bm{x}}_{i}\frac{\varepsilon_{ijk}}{8\pi}\bm{n}\cdot[\partial_{j}\bm{n}\times\partial_{k}\bm{n}],$
(19)

where ${\bm{r}_{\textsc{bp}}}$ denotes the Bloch point position. If
$\mathfrak{R}=\mathfrak{R}_{\hat{\bm{z}}}(\alpha)$ describes the rotation by a
constant angle $\alpha$ around the $z$-axis, then the BP has a vortex
structure within the $(x,y)$-plane similarly to the BPs shown in Fig. 4(b) of
the main text. In zero magnetic field, the exchange energy density of the BP
$\mathscr{E}_{\text{ex}}=A(\partial_{i}\bm{n})^{2}=2A/r^{2}$ and the DMI
energy density is
$\mathscr{E}_{\textsc{dmi}}=D\,\bm{n}\left[\bm{\nabla}\times\bm{n}\right]=2D\sin\alpha\frac{z}{r^{2}}$.
The DMI contribution vanishes for the integrated energy $E_{\textsc{bp}}$ of a
sphere with radius $R_{\textsc{bp}}$ centered on the BP. Finally, we reproduce
the well known [64, 65] result for the BP energy $E_{\textsc{bp}}^{0}=8\pi
AR_{\textsc{bp}}$.

An important distinction in the present case is that the angle $\alpha$ is not
a constant because the BP is immersed into the helical state, namely
$\alpha=\frac{2\pi}{\lambda_{h}}(z-z_{\textsc{bp}})+\frac{\pi}{2}(1-c)$, where
$z_{\textsc{bp}}$ is the BP position. The last summand appears by demanding
that at $z=0$ the asymptotic value
$\bm{n}(\rho\to\infty,\chi=0,z=0)=\hat{\bm{x}}$. For such a position dependent
angle $\alpha(z)$, the BP energy has the form

$E_{\textsc{bp}}=E_{\textsc{bp}}^{0}\left[1-\left(\frac{2\pi}{3}\right)^{2}\frac{R_{\textsc{bp}}^{2}}{\lambda_{h}^{2}}+c\,\mathcal{A}\cos\left(2\pi\frac{z_{\textsc{bp}}}{\lambda_{h}}\right)\right],$
(20)

with $\mathcal{A}=\frac{\lambda_{h}}{\pi R_{\textsc{bp}}}\left[\text{Si}(2\pi
R_{\textsc{bp}}/\lambda_{h})-\sin(2\pi
R_{\textsc{bp}}/\lambda_{h})\right]\approx
2\left(\frac{2\pi}{3}\right)^{2}R_{\textsc{bp}}^{2}/\lambda_{h}^{2}$ where
$\text{Si}(x)$ is the sine integral and the approximation is valid for
$R_{\textsc{bp}}\ll\lambda_{h}$. Minimization of $E_{\textsc{bp}}$ with
respect to $z_{\textsc{bp}}$ results in the equilibrium BP positions

$z_{\textsc{bp}}^{c}=\lambda_{h}\left[\frac{c+1}{4}+m\right],\qquad
m\in\mathbb{Z}.$ (21)

These equilibrium BP positions coincide with the ones obtained by means of
micromagnetic simulations for sd${}^{\rm Bp}_{1}$, see Fig. 4(b) of the main
text. The resulting BP energy at equilibrium for
$R_{\textsc{bp}}\ll\lambda_{h}$ is then given by $E_{\textsc{bp}}\approx
E_{\textsc{bp}}^{0}\left[1-\frac{(2\pi)^{2}}{3}R_{\textsc{bp}}^{2}/\lambda_{h}^{2}\right]$.

Let us now estimate the energy of a BP which has an anti-vortex structure in
its horizontal cross-section, i.e.
$\bm{n}_{\textsc{bp}}=c\,\mathfrak{R}_{\hat{\bm{z}}}(\alpha)\mathfrak{R}_{\hat{\bm{x}}}(\pi)\bm{r}/|\bm{r}|$.
A direct calculation shows that the corresponding Bloch point energy
$E_{\textsc{bp}}^{\textsc{av}}=E_{\textsc{bp}}^{0}\left[1-\left(\frac{2\pi}{3}\right)^{2}\frac{R_{\textsc{bp}}^{2}}{\lambda_{h}^{2}}\right]$
is independent on the BP position, and it is higher than the energy of the BP
with the vortex-like horizontal cross-section. The latter supports the
numerically obtained fact that the anti-vortex screw dislocations
sd${}_{-1}^{\pm}$ do not contain Bloch points.

We note that the statics and dynamics of a single BP with a vortex cross-
section in a cubic helimagnet was recently studied numerically [66] as well as
experimentally [67].