# Field line winding of braided vector fields in tubular subdomains

Christopher B Prior and Anthony R Yeates  Department of Mathematical Sciences,
Durham University, South Road, Durham, DH1 3LE, UK
<EMAIL_ADDRESS>
###### Abstract

Braided vector fields on spatial subdomains homeomorphic to the cylinder play
a crucial role in applications such as solar and plasma physics, relativistic
astrophysics, fluid and vortex dynamics, elasticity, and bio-elasticity. Often
the vector field’s topology – the entanglement of its field lines – is non-
trivial, and can play a significant role in the vector field’s evolution. We
present a complete topological characterisation of such vector fields (up to
isotopy) using a quantity called field line winding. This measures the
entanglement of each field line with all other field lines of the vector
field, and may be defined for an arbitrary tubular subdomain by prescribing a
minimally distorted coordinate system. We propose how to define such
coordinates, and prove that the resulting field line winding distribution
uniquely classifies the topology of a braided vector field. The field line
winding is similar to the field line helicity considered previously for
magnetic (solenoidal) fields, but is a more fundamental measure of the field
line topology because it does not conflate linking information with field
strength.

††: J. Phys. A: Math. Gen.

  * November 2018

##  1 Introduction

The entanglement of vector field integral curves (field lines) in tubular
subdomains homeomorphic to the cylinder has long been of wide interest. For
example, in stellar interiors, twisted bundles of magnetic field lines – known
as magnetic flux ropes – are potentially important for the generation of large
scale magnetic fields (?, ?, ?, ?). In stellar atmospheres such as the Sun’s
corona, twisted and braided magnetic fields play a crucial role in dynamic
phenomena such as coronal heating, jet formation, and coronal mass ejections
(?, ?, ?, ?, ?). Another area where the dynamic entanglement of magnetic field
is crucial is the formation of relativistic jets from black holes and neutron
stars (?, ?, ?, ?). In thin body elasticity, structures such as ropes, cables,
and biopolymers are treated as thin elastic tubes which can be internally
twisted and which often react to this twisting by forming looped and knotted
structures (?, ?, ?, ?, ?). A particularly well known example is plectoneme
formation, whereby an elastic tube submitted to increasing twisting becomes
unstable and loops at its centre, forming a self-contacting loop. As the twist
is further increased the number of loops increase and the tube becomes
supercoiled. DNA is an example of a biopolymer which exhibits this
supercoiling (and where twisting can be applied by optical tweezer
experiments, e.g., ?). Supercoiling acts as a mechanism for large DNA
molecules to compactify in order to fit into their local environment (?). In
all of these varied systems, a crucial aspect is their inability to
disentangle, at least completely. This often informs their physical behaviour,
and means that quantifying entanglement is important for understanding the
systems themselves.

Motivation for our work comes from the particular context of magnetised
plasmas, where the standard tool for quantifying entanglement is a volume
integral called magnetic helicity, $H$ (?, ?), analogous to a similar quantity
in fluid dynamics (?, ?, ?). For a “closed” magnetic field whose field line
curves are tangent to the boundary, the value of this integral is invariant
under an ideal evolution because field lines cannot reconnect and change their
topology. As a result, having non-zero $H$ constrains the amount of magnetic
energy that can be released (?). In such a magnetically-closed volume, it has
long been known that $H$ may be written as an average of the Gauss linking
integral between all pairs of magnetic field line curves (?, ?, ?). Recently,
? showed that the $H$ admits an analogous interpretation for so-called
“braided” magnetic fields, where magnetic field lines connect between two
planar boundaries rather than being tangent. In this case, $H$ is gauge
dependent, but the authors showed that there is a particular gauge, the
“winding gauge”, in which $H$ is the equal to the average winding number
between all pairs of field line curves. This extends the topological
interpretation of helicity beyond the restricted case of tangent fields, which
are not directly relevant for applications such as stellar atmospheres.

The main limitation of the winding-gauge helicity of ? is that it is defined
only for domains foliated by planes, for which the standard definition of
winding number can be used. However, in many contexts it would be useful to
consider tubular domains of more general shape. One example would be the
concentrated magnetic flux ropes that form during flux emergence into the
solar atmosphere. Such flux ropes typically occupy curved and distorted
domains that cannot be represented as foliations of parallel planes (?). Yet
to characterise the internal entanglement of such structures it is necessary
to isolate them from the surrounding magnetic field. Recent work has suggested
that these structures may be more stable if they have a complex, “braided”
internal field-line topology, rather than a uniformly twisted structure (?).
Even in stable loops, ? showed that different internal braiding patterns can
lead to significantly different internal heating. It is important to note that
the ? definition of winding-gauge helicity will fail not only if the tubular
domain as a whole has a curved shape, but also if its end “caps” are not flat.
Recently, ? have shown that accounting for such wrinkling of the Sun’s surface
can significantly affect the measured input of magnetic helicity into the
solar atmosphere. Observational studies of such helicity injection are
important because they can potentially be used try to predict the ensuing
behaviour of the atmospheric magnetic field (?, ?, ?), as well as providing
insight into the unobservable magnetic structure in the solar interior. Thus
there is a need to consider winding numbers on more general domains, so as to
quantify the entanglement of vector fields on those domains. This is the main
aim of this paper.

To fully quantify the entanglement of a braided vector field, it is not
sufficient to consider only a single, global integral like magnetic helicity,
$H$. To give an extreme example, since $H$ is a signed quantity, one can have
a topologically non-trivial vector field with zero total $H$, such as the
Borromean rings (?) or the equivalent magnetic braid (?). Moreover, in real
systems that are not perfectly ideal, it is often the case that there are
substantial local changes in entanglement, even when the total $H$ (the
average entanglement) is conserved to a good approximation (?, ?). A full
understanding of the field line topology therefore requires the study of
finer-grained invariants. As such, our main quantity of interest in this paper
is not the analogue of $H$ (which would be the average winding among all pairs
of field lines), but rather the analogue of field line helicity (?, ?, ?, ?,
?). This we call the field line winding, as it is defined for each field line
and measures the average winding of this field line with all others.

Our new measure of field line winding allows not only for more general shapes
of domain but also for more general classes of vector field. We still assume
that the vector field is “braided” (defined precisely in Section 2), but
unlike ? we do not require it to be solenoidal (divergence free). As will be
discussed in Section 6, this amounts to removing the field-strength
information from the winding measure. Magnetic helicity, for example, has
units of magnetic flux squared, and is effectively a confluence of both
topological and strength information. As shown by ?, a given set of field line
curves can have different helicity depending on the field strength. In this
paper, we focus purely on how to uniquely describe the field line topology,
without regard to the strength. This unweighted, purely topological measure
applies equally to non-solenoidal fields, whereas the magnetic helicity, for
example, is an ideal invariant only thanks to the solenoidal nature of the
magnetic field.

The paper is organised as follows. After stating our assumptions on the vector
field and domain in Section 2, we give the general definition of winding
numbers and field line winding in Section 3. These definitions rely on
establishing a least distorted vector field on the domain, whose field lines
have no mutual winding. This is addressed in Section 4. In Section 5, we prove
that the field line topology is completely determined by the distribution of
field line winding. This is an extension of an earlier completeness result
that applied only to a more restricted class of solenoidal vector fields (?,
?, ?). As well as extending the class of domains and vector fields considered,
we prove that the field line winding measure uniquely determines not only the
field line mapping (as in the earlier papers) but also whether or not the two
vector fields can be linked by an ideal evolution (an end-vanishing isotopy).
Thus we also strengthen the earlier result for magnetic fields. The relation
of our field line winding invariant to magnetic helicity is discussed in
Section 6, before concluding in Section 7.

##  2 Assumptions and notation

Figure 1: A tubular subdomain $\mathcal{M}$ is an embedding of the unit
cylinder in $\mathbb{R}^{3}$ given by a homeomorphism
$F:\mathcal{C}\to\mathbb{R}^{3}$. The discs $D_{z}$ of constant $z$ in the
cylinder $\mathcal{C}$ map to distorted surfaces $S_{z}$ that foliate
$\mathcal{M}$. Shown in green are a field line $F(\gamma)\in\mathcal{M}$ and
its preimage curve $\gamma\in\mathcal{C}$. A crucial part of this study will
be to determine a minimally-distorted choice for the map $F$, within a given
subdomain $\mathcal{M}$.  Figure 2: Three curves that could be field lines of
a braided vector field as defined in Section 2. The green curve would not be
monotone in $z$ but is still valid provided that the corresponding ${\bf v}$
can be deformed by isotopy to a field with field lines that are monotone in
$z$.

In this paper, we consider vector fields ${\bf v}$ defined on subdomains
$\mathcal{M}\subset\mathbb{R}^{3}$. Both are restricted, in the following
ways. The subdomain $\mathcal{M}$ is “tubular”, meaning that it may be thought
of as an embedding of the solid unit cylinder $\mathcal{C}$ in
$\mathbb{R}^{3}$, determined by some homeomorphism
$F:\mathcal{C}\to\mathbb{R}^{3}$, as in Figure 1. This homeomorphism $F$ is by
no means unique, and an important part of this paper will be to choose a
suitable $F$. Note that the boundary $\partial\mathcal{M}=S_{0}\cup S_{1}\cup
S_{s}$ is the union of two end caps $S_{0},S_{1}$ and a side boundary $S_{s}$,
and this composite boundary is assumed Lipschitz continuous. This means that
$\mathcal{M}$ can be foliated by a set of surfaces $S_{z},\,z\in[0,1]$. We
will choose $F$ so that these surfaces are the images of the discs $D_{z}$ of
constant $z$ in $\mathcal{C}$, where we define Cartesian coordinates on
$\mathcal{C}$ as $(x^{1},x^{2},z)$. This notation recognizes the special role
of the third coordinate in defining the foliation. In effect, $z$ becomes an
axial coordinate for $\mathcal{M}$. Throughout this paper we will denote
points in $\mathcal{C}$ by $x$ and points in $\mathcal{M}$ by $y$.

The vector field ${\bf v}$ is “braided”, meaning that it is Lipschitz
continuous, non-zero everywhere in $\mathcal{M}\cup\partial\mathcal{M}$, and
satisfies the boundary conditions

$\hat{\bf z}\cdot{\bf v}>0\quad\textrm{on $S_{0}$},\qquad\hat{\bf z}\cdot{\bf
v}>0\quad\textrm{on $S_{1}$},\qquad\hat{\bf n}\cdot{\bf v}=0\quad\textrm{on
$S_{s}$},$ (1)

where $\hat{\bf z}$ is a unit vector normal to $S_{z}$ pointing in the
direction of increasing $z$ and $\hat{\bf n}$ is the unit normal to $S_{s}$.
Thus the integral curves/field lines of ${\bf v}$ are all rectifiable curves
that connect from $S_{0}$ to $S_{1}$, as shown for the three example field
lines in Figure 2. The final assumption on ${\bf v}$ is that its field lines
can be deformed by isotopy to a vector field satisfying $\hat{\bf z}\cdot{\bf
v}>0$ on every $S_{z}$. By an isotopy we mean a continuous set of
homeomorphisms of the curves (which are automorphisms of the domain
$\mathcal{M}$). For such fields, we can define a mapping from $S_{0}$ to
$S_{z}$, or equivalently from $D_{0}$ to $D_{z}$ by following the field lines.
These mappings are diffeomorphisms thanks to the Picard-Lindelöf theorem (see,
e.g, ?). In particular, we define a diffeomorphism $f(1)$ from $D_{0}$ to
$D_{1}$ that will be used in Section 5.

##  3 Field line winding

Our goal in this paper is to describe the field-line topology of a braided
vector field ${\bf v}$ in $\mathcal{M}$. In other words, properties of the
field lines that remain unchanged under a continuous isotopy/deformation of
${\bf v}$. Since the field lines exit $\mathcal{M}$ through $S_{0}$ and
$S_{1}$, we must qualify this statement: we consider properties that remain
invariant under isotopies of ${\bf v}$ _that vanish on $S_{0}$ and $S_{1}$_.
These we call end-vanishing isotopies.

Our topological invariants are based on winding numbers between field lines.
In this paper, we take advantage of the embedding of $\mathcal{M}$ and define
the winding number between any two curves in $\mathcal{M}$ by calculating the
winding number between their two preimage curves in $\mathcal{C}$. An
alternative approach would be to work in $\mathcal{M}$ directly by choosing a
metric structure. The advantage of our approach is that it allows us to define
the angle in a simple and clear manner (i.e., on the cylinder), as well as
simplifying the proof in Section 5. Figure 3(a) shows field lines of an
example vector field on a tubular domain $\mathcal{M}$, and Figure 3(b) shows
the preimages of these curves in $\mathcal{C}$ for a particular choice of the
embedding map $F$. How best to choose $F$ will be addressed in Section 4.

(a)  (b)

(c)

.

Figure 3:  Example illustrating the mapping of field line curves from a
tubular domain $\mathcal{M}$ (a) to the cylinder $\mathcal{C}$ (b), along with
the resulting distribution of field-line winding $L_{\bf v}$ in $D_{0}$ (c).
The green curves in (a) and (b) show a particular field line $\gamma$ and its
preimage. The interpretation of $L_{\bf v}(\gamma)$ is its average winding
number with each of the other curves. Overall, this vector field shows
significant entanglement due to the intermixing of positive and negative
winding.

###  3.1 Pairwise winding number

Figure 4: Definition of the angle $\Theta$ between two curves $\gamma$,
$\widetilde{\gamma}$ on $\mathcal{C}$.

Given two preimage curves $\gamma$, $\widetilde{\gamma}$ in $\mathcal{C}$, the
most basic topological invariant is their pairwise winding number. In any
plane $D_{z}$, we can use the Cartesian components of the two curves $\gamma$
and $\widetilde{\gamma}$ to define the “angle” between them,

$\Theta(\gamma,\widetilde{\gamma},z)=\arctan\left(\frac{\gamma_{2}(z)-\widetilde{\gamma}_{2}(z)}{\gamma_{1}(z)-\widetilde{\gamma}_{1}(z)}\right),$
(2)

as shown in Figure 4. The net change in this angle as we follow the curves
from $z=0$ to $z=1$ is the pairwise winding number

$\displaystyle\mathcal{L}(\gamma,\widetilde{\gamma})$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{1}\frac{\mathrm{d}}{\mathrm{d}z}\Theta(\gamma,\widetilde{\gamma},z)\,\mathrm{d}z$
(3)
$\displaystyle=\frac{1}{2\pi}\int_{0}^{1}\frac{(\gamma_{1}-\widetilde{\gamma}_{1})\frac{\mathrm{d}}{\mathrm{d}z}(\gamma_{2}-\widetilde{\gamma}_{2})-(\gamma_{2}-\widetilde{\gamma}_{2})\frac{\mathrm{d}}{\mathrm{d}z}(\gamma_{1}-\widetilde{\gamma}_{1})}{(\gamma_{1}-\widetilde{\gamma}_{1})^{2}+(\gamma_{2}-\widetilde{\gamma}_{2})^{2}}\,\mathrm{d}z,$
(4)

where primes denote differentiation by $z$. Alternatively, we may write

$\mathcal{L}(\gamma,\widetilde{\gamma})=\frac{1}{2\pi}\Big{(}\Theta(\gamma,\widetilde{\gamma},1)-\Theta(\gamma,\widetilde{\gamma},0)\Big{)}+N(\gamma,\widetilde{\gamma}),$
(5)

where $N(\gamma,\widetilde{\gamma})$ counts the (signed) number of branch cut
crossings as we follow the curves in $z$. This makes clear that
$\mathcal{L}(\gamma,\widetilde{\gamma})$ is invariant under any isotopy of the
curves that does not move their end-points. In other words, it is a
topological invariant.

? showed that the definition of $\mathcal{L}(\gamma,\widetilde{\gamma})$ may
be generalised to allow for non-monotonic curves like the green curve in
Figure 2. Such a curve $\gamma$ is split into $n+1$ sections
$\gamma^{(0)},\ldots,\gamma^{(n)}$ using the $n$ turning points where
$\mathrm{d}\gamma_{z}/\mathrm{d}z=0$. For each section, we define the
indicator function

$\sigma^{(i)}=\left\\{\begin{array}[]{@{\kern 2.5pt}lL}\hfill 1&if
$\mathrm{d}\gamma^{(i)}_{z}/\mathrm{d}z>0$,\\\ \hfill-1&if
$\mathrm{d}\gamma^{(i)}_{z}/\mathrm{d}z<0$,\\\ \hfill
0&otherwise.\end{array}\right.$ (6)

We split $\widetilde{\gamma}$ and define $\widetilde{\sigma}^{(j)}$ in a
similar way. Then the pairwise winding number is the sum

$\mathcal{L}(\gamma,\widetilde{\gamma})=\sum_{i=0}^{n}\sum_{j=0}^{\widetilde{n}}\frac{\sigma^{(i)}\widetilde{\sigma}^{(j)}}{2\pi}\int_{z_{ij}^{\rm
min}}^{z_{ij}^{\rm
max}}\frac{\mathrm{d}}{\mathrm{d}z}\Theta(\gamma^{(i)},\widetilde{\gamma}^{(j)},z)\,\mathrm{d}z,$
(7)

where $[z_{ij}^{\rm min},z_{ij}^{\rm max}]$ is the mutual range of $z$ values
(if any) shared by the curve sections $\gamma^{(i)}$ and
$\widetilde{\gamma}^{(j)}$. Once again,
$\mathcal{L}(\gamma,\widetilde{\gamma})$ is invariant to any isotopy of the
curves that fixes their endpoints. It reduces to (3) if both curves have only
a single section stretching from $z=0$ to $z=1$. (Although not needed here, ?
also showed that when the curves are closed,
$\mathcal{L}(\gamma,\widetilde{\gamma})$ is equal to their Gauss linking
integral, hence the notation “$\mathcal{L}$”.)

(a)  (b)

Figure 5: How two different choices of mapping $F$ can lead to different
pairwise winding number $\mathcal{L}(\gamma,\widetilde{\gamma})$. Here the red
and green curves show $F(\gamma)$ and $F(\widetilde{\gamma})$, while yellow
lines indicate the surface that $F(\gamma)$ should belong to in each
coordinate system to give constant $\Theta(\gamma,\widetilde{\gamma},z)$. In
(a) $F(\gamma)$ lies in this surface, so
$\mathcal{L}(\gamma,\widetilde{\gamma})=0$, but in (b) it does not and
$\mathcal{L}(\gamma,\widetilde{\gamma})\neq 0$.

It is important to realize that the value of
$\mathcal{L}(\gamma,\widetilde{\gamma})$ depends on the choice of embedding
map $F$. This is illustrated by Figure 5. Nevertheless,
$\mathcal{L}(\gamma,\widetilde{\gamma})$ is a topological invariant for any
(fixed) choice of $F$, under deformations that vanish on $S_{0}$ and $S_{1}$.

###  3.2 Field line winding of a vector field

In principle, we could characterise the field line topology of ${\bf v}$ by
the set of all pairwise winding numbers
$\mathcal{L}(\gamma,\widetilde{\gamma})$ between all pairs of preimage field
lines. However, we will show in Section 5 that this would entail significant
redundancy. A more succinct description of the field line topology is given by
a quantity we call the field line winding. This is defined for each field line
$\gamma$, and is simply the average winding number of $\gamma$ with all other
field lines in $\mathcal{C}$, i.e.,

$L_{\bf
v}(\gamma)=\int_{D_{0}}\mathcal{L}\big{(}\gamma,\widetilde{\gamma}(x)\big{)}\,\mathrm{d}^{2}x.$
(8)

Here $\widetilde{\gamma}(x)$ denotes the field line starting from a point
$x\in D_{0}$. Since every field line $\gamma$ passes through $D_{0}$, $L_{\bf
v}$ defines a scalar distribution on $D_{0}$. An example of such a
distribution is shown in Figure 3(c).

By virtue of its definition, $L_{\bf v}$ is invariant under end vanishing
isotopies of ${\bf v}$ that vanish on $S_{0}$, $S_{1}$. Like
$\mathcal{L}(\gamma,\widetilde{\gamma})$, the field line winding $L_{\bf
v}(\gamma)$ is dependent on the choice of mapping from $\mathcal{C}$ to
$\mathcal{M}$; we will shortly fix this choice and hence the definition of
$\mathcal{L}$.

An alternative way to write (8) is to use Cartesian coordinates
$(x^{1},x^{2},z)$ on $\mathcal{C}$, then (4) gives

$L_{\bf
v}(\gamma)=\frac{1}{2\pi}\int_{0}^{1}\int_{D_{z}}\frac{(x^{1}-\gamma_{1})\frac{\mathrm{d}}{\mathrm{d}z}[\widetilde{\gamma}_{2}(x)-\gamma_{2}]-(x^{2}-\gamma_{2})\frac{\mathrm{d}}{\mathrm{d}z}[\widetilde{\gamma}_{1}(x)-\gamma_{1}]}{(x^{1}-\gamma_{1})^{2}+(x^{2}-\gamma_{2})^{2}}\,\mathrm{d}^{2}x\rmd
z,$

where $\widetilde{\gamma}(x)$ is the curve passing through the point
$x=(x^{1},x^{2},z)\in D_{z}$.

##  4 Choice of embedding map

For measuring field line winding on a fixed subdomain $\mathcal{M}$, any
choice of embedding map $F:\mathcal{C}\to\mathcal{M}$ will generate a
topological measure $L_{\bf v}$ that is invariant under end-vanishing
isotopies. However, as we saw in Figure 5, the actual value of $L_{\bf v}$ for
each field line depends on the chosen mapping. In order to allow for
comparison between vector fields on different $\mathcal{M}$, it is useful to
define $F$ uniquely and fix an absolute $L_{\bf v}$. In this section, we
propose to do this by identifying the “least distorted” vector field on
$\mathcal{M}$, that follows the shape of the tube as simply as possible. We
then define $F$ so that the field lines of this least distorted field are the
images of vertical lines on $\mathcal{C}$. This way, they will all have
$L_{\bf v}\equiv 0$. Vector fields with non-trivial field line topology will
then have $L_{\bf v}\neq 0$, at least for a subset of their field lines.

###  4.1 The least distorted field

A natural candidate for least-distorted braided vector field on $\mathcal{M}$
is a harmonic vector field ${\bf u}=\nabla\phi$ that satisfies

$\displaystyle\nabla^{2}\phi=0\quad\textrm{in $\mathcal{M}$},$ (9)
$\displaystyle\hat{\bf z}\times\nabla\phi=0\quad\textrm{on $S_{0}$, $S_{1}$},$
(10) $\displaystyle\hat{\bf n}\cdot\nabla\phi=0\quad\textrm{on $S_{s}$}.$ (11)

Here we reiterate that $\hat{\bf z}$ is the unit normal to the cross-sectional
surfaces $S_{z}$, and $\hat{\bf n}$ is the unit normal to the side boundary
$S_{s}$. Conditions (9)–(11) define ${\bf u}$ uniquely up to magnitude,
disregarding the trivial case of constant $\phi$. Condition (10) implies that
$\phi$ is constant on each of the end caps, so that the integral curves are
normal to these end caps. The actual values of these two constants control the
magnitude and sign of ${\bf u}$, but not the shape of its integral curves.
That this field can be assumed to exist on $\mathcal{M}$ is a result of
Theorem 1.1 of ?, which provides a mixed boundary condition Hodge
decomposition on Lipschitz domains. In particular the dimension of the
harmonic field space, given by (1.6) in that paper, is $1$ for a cylinder;
this dimension corresponds to the magnitude of the field ${\bf u}$.

To motivate the choice of condition (10), consider any other harmonic field
$\nabla\phi^{\prime}=\nabla(\phi+\psi)$ in $\mathcal{M}$ that satisfies (11)
and has the same flux as $\nabla\phi$ through $S_{0}$ and $S_{1}$, meaning
that

$\int_{S_{0}}\hat{\bf z}\cdot\nabla\psi\,\mathrm{d}S=\int_{S_{1}}\hat{\bf
z}\cdot\nabla\psi\,\mathrm{d}S=0.$ (12)

It follows from these boundary conditions, and from the constancy of $\phi$ on
$S_{0}$ and $S_{1}$, that the Dirichlet (stretching) energy of
$\nabla\phi^{\prime}$ is the orthogonal sum

$\int_{\mathcal{M}}|\nabla\phi^{\prime}|^{2}\,\mathrm{d}V=\int_{\mathcal{M}}|\nabla\phi|^{2}\,\mathrm{d}V+\int_{\mathcal{M}}|\nabla\psi|^{2}\,\mathrm{d}V.$
(13)

Therefore $\nabla\phi$ is the harmonic field in $\mathcal{M}$ that minimizes
this energy for a given axial flux. In this sense, it has the least distorted
integral curves. This is related to the notion of harmonic coordinates in
Riemannian geometry (?), which may be thought of as minimizing the Dirichlet
energy of a coordinate map from $\mathcal{M}$ to $\mathbb{R}^{3}$.

Figure 6: Example domains and computations of their minimally distorted fields
${\bf u}$.

Examples of this field ${\bf u}$ are shown in Figure 6 for various domains,
computed using the finite-element code IFEM2 in MATLAB (?). Note that the
curves of this least distorted field flow through the domain, contouring to
its shape and expanding (contracting) when the tube does. They are the natural
analogue of straight lines in a Cartesian domain.

###  4.2 Definition of the embedding map

It is most convenient to specify $F^{-1}$ rather than $F$ itself directly. Our
fundamental idea is that field lines of ${\bf u}$ should map to vertical lines
in $\mathcal{C}$, but this leaves considerable freedom in the choice of
$F^{-1}$, both in the $z$ coordinate (effectively distance along the field
lines) and in which field line of ${\bf u}$ maps to which vertical line in
$\mathcal{C}$.

To set the $z$ coordinate, note that ${\bf u}=\nabla\phi$ naturally defines a
foliation $\\{S_{z}\\}$ of $\mathcal{M}$ by taking each $S_{\phi}$ to be the
surface $\phi=\textrm{constant}$. Given the freedom in scaling the magnitude
of ${\bf u}$, we may always arrange that $\phi=0$ on $S_{0}$ and $\phi=1$ on
$S_{1}$. The fact that $\phi$ gives a valid foliation relies on the fact that
${\bf u}$ is non-zero everywhere on the interior of $\mathcal{M}$ (see A).

Choosing which field line of ${\bf u}$ maps to which vertical line is
equivalent to choosing $F^{-1}$ on one of the end caps, say $S_{0}$.
Equivalently, we must choose functions $(x^{1}(y),x^{2}(y))$ for $y\in S_{0}$.
To avoid measuring any spurious pairwise winding (cf. Figure 5), we require
each of the functions $x^{1}$ and $x^{2}$ to be harmonic in the two-
dimensional surface $S_{0}$, that is

$\nabla^{2}_{S_{0}}x^{1}|_{S_{0}}=\nabla^{2}_{S_{0}}x^{2}|_{S_{0}}=0.$ (14)

where $\nabla^{2}_{S_{0}}$ indicates the Laplacian operator in the surface
$S_{0}$. To specify a unique solution we also need to impose boundary
conditions on $\partial S_{0}$. In order to ensure that $F^{-1}$ (and hence
$F$) is one-to-one, these boundary conditions must take the form
$x^{1}=\cos(\theta(s)),x^{2}=\sin(\theta(s))$, where $\theta(s)$ is an angle
function of arclength $s$ along $\partial S_{0}$. An example is shown in
Figure 7, where we make the simple choice $\theta(s)=2\pi s/L$ for some
arbitrary point $s=0$, where $L$ is the perimeter length of $S_{0}$.

(a)  (b)

Figure 7: Example illustrating our proposed definition of the mapping $F$ on
the lower boundary, $S_{0}$. Black lines in (a) show the images under $F$ of
$x^{1}$ and $x^{2}$ contours, whereas those in (b) show the images of polar
coordinate $(r,\theta)$ contours under the same $F$.

##  5 Topological classification

In this section, we will show that the field-line winding $L_{\bf v}$
distribution for a braided vector field ${\bf v}$ completely determines the
topology of its field lines, in the following sense.

######  Theorem 1.

Let ${\bf v}$, ${\bf v}^{\prime}$ be two braided vector fields on the same
domain $\mathcal{M}$ whose field lines on $S_{s}$ are linked by an end-
vanishing isotopy. Then the field lines of ${\bf v}$ and those of ${\bf
v}^{\prime}$ within $\mathcal{M}$ can be linked by an end-vanishing isotopy if
and only if $L_{\bf v}=L_{{\bf v}^{\prime}}$ on all of $D_{0}$.

Notice that this is purely a result about the field line curves: the
magnitudes of ${\bf v}$ and ${\bf v}^{\prime}$ do not matter, because $L_{\bf
v}$ and $L_{{\bf v}^{\prime}}$ depend only on the geometry of the curves. The
boundary requirement on $S_{s}$ relates to a fundamental fact regarding the
classification group of parametrized diffeomorphisms of the unit disc (?, ?).
Essentially, if the corresponding field line mappings $f(1)$, $f^{\prime}(1)$
differ on the side boundary, then this difference could be compensated by an
opposite rotation on the interior, leading to both fields having the same
$L_{\bf v}$.

###### Proof of Theorem 1.

One direction is immediate: we know that $L_{\bf v}$ is invariant under any
isotopy of ${\bf v}$ that vanishes on $S_{0}$ and $S_{1}$. The other direction
is a deeper result. Our strategy will be to first show that equality of
$L_{\bf v}$ and $L_{{\bf v}^{\prime}}$ implies equality of the corresponding
field line mappings $f(1)$, $f(1)^{\prime}$ from $D_{0}$ to $D_{1}$. Then, we
will show that this implies the existence of an end vanishing isotopy between
the field lines of ${\bf v}$ and those of ${\bf v}^{\prime}$.

To show equality of the field line mappings $f(1)$ and $f^{\prime}(1)$,
consider the gradient of $L_{\bf v}(\gamma)$ with respect to the field line
startpoint $(r,\theta)\in D_{0}$. This is most succinctly expressed in the
language of differential forms as $\mathrm{d}_{\perp}L_{\bf v}$, and we claim
that

$\mathrm{d}_{\perp}L_{\bf
v}=f^{*}\alpha-\alpha,\qquad\textrm{where}\quad\alpha(r,\theta,z)=\frac{r^{2}}{2}\mathrm{d}\theta.$
(15)

To see this, note that we can write $L_{\bf v}$ in terms of the field line
mapping $f$ as

$L_{\bf
v}(\gamma)=\int_{D_{0}}\Big{[}\frac{1}{2\pi}\Big{(}f^{*}\Theta-\Theta\Big{)}+N(\gamma,\widetilde{\gamma}(x))\Big{]}\widetilde{r}\,\mathrm{d}^{2}x,$
(16)

where $\Theta(\gamma,\widetilde{\gamma}(x),z)$ is the angle made in the
surface $D_{z}$ by the curve $\gamma$ and another curve $\widetilde{\gamma}$
rooted at a point $x\in D_{0}$. Treating this as a function of the
$\gamma$-startpoint $(r,\theta)\in D_{0}$, we differentiate to find

$\displaystyle\mathrm{d}_{\perp}L_{\bf v}$
$\displaystyle=\frac{1}{2\pi}\int_{D_{0}}\Big{(}f^{*}\mathrm{d}_{\perp}\Theta-\mathrm{d}_{\perp}\Theta\Big{)}\,\mathrm{d}^{2}x.$
(17)

One may show by explicit calculation (see B) that

$\frac{1}{2\pi}\int_{D_{0}}\mathrm{d}_{\perp}\Theta\,\mathrm{d}^{2}x=\alpha,$
(18)

from which (15) follows.

Now consider the composite map $g:D_{0}\to D_{0}$, defined by
$g=f(1)^{-1}\circ f^{\prime}(1)$. Under our assumption $L_{\bf v}=L_{{\bf
v}^{\prime}}$, it follows from equation (15) that $g^{*}\alpha=\alpha$ (?).
Equivalently, in components,

$\displaystyle\frac{g_{r}^{2}}{2}\frac{\partial g_{\theta}}{\partial
r}=0,\qquad\frac{g_{r}^{2}}{2}\frac{\partial
g_{\theta}}{\partial\theta}=\frac{r^{2}}{2}.$ (19)

The first equation shows that $g_{\theta}(r,\theta)=G(\theta)$, but the fact
that $f(1)=f^{\prime}(1)$ on the side boundary (by our assumption on $S_{s}$)
means that $g_{\theta}(r,\theta)=\theta$. The second equation then shows that
$g_{r}(r,\theta)=r$, so $g$ is the identity and hence $f(1)=f^{\prime}(1)$.

This leaves us to show that there exists an isotopy between any two braided
vector fields with the same field line mapping $f(1)$. For this we utilise
some established results. The field lines $f(z)$ are parametrised curves in
the group of diffeomorphisms of the unit disc, $\mathrm{Diffeo}(D^{2})$. Since
$f(0)=f^{\prime}(0)=\mathrm{id}$ and $f(1)=f^{\prime}(1)$, the composite map

$f_{c}(z)=\left\\{\begin{array}[]{cc}f(z)&z\in[0,1)\\\
f^{\prime}(2-z)&z\in[1,2].\end{array}\right.$ (20)

forms a closed curve in this group which is anchored at the identity. If this
curve can be shrunk to the identity within the group by a homotopy then the
field lines can be linked by an isotopy. Thus we are interested in the
homotopy group $\pi_{1}(\mathrm{Diffeo}(D^{2}))$ of curves in
$\mathrm{Diffeo}(D^{2})$ (anchored at the identity). A number of established
facts related to this group allow us to answer the question. Firstly, there is
a restriction map

$\mathrm{Diffeo}(D^{2})\rightarrow\mathrm{Diffeo}(S^{1}),$ (21)

obtained by restricting each diffeomorphism in $D^{2}$ to the boundary
$\partial D^{2}=S^{1}$. This is a fibration with fiber
$\mathrm{Diffeo}_{\partial D^{2}}(D^{2})$, namely diffeomorphisms of the disc
that fix the boundary. This implies that the map
$\pi_{1}(\mathrm{Diffeo}((D^{2}))\rightarrow\pi_{1}(\mathrm{Diffeo}(S^{1}))$
is an isomorphism – see ? in conjunction with Theorem 4.41 in ?. Finally there
is a homotopy equivalence between $\mathrm{Diffeo}(S^{1})$ and the group
$O^{2}$ of rotations (see ?). Putting this information together implies that
the set of closed curves in the space $\mathrm{Diffeo}(D^{2})$, from which our
map $f_{c}$ is drawn, can be categorised up to homotopy by their behaviour
when restricted to the boundary of the cylinder. Then the fact that the
homotopy group has the structure of $O^{2}$ implies that closed curves in
$\mathrm{Diffeo}(D^{2})$ are homotopic to the identity if and only if they
describe no net integer rotation on the boundary. The boundary assumption in
the theorem ensures this for $f_{c}$. ∎

##  6 Relation to helicity

Since each pairwise winding number $\mathcal{L}(\gamma,\widetilde{\gamma})$ is
itself invariant under end-vanishing isotopies of ${\bf v}$, it is clear that
any functional

$W_{\bf
v}(\gamma)=\int_{D_{0}}w(x)\mathcal{L}(\gamma,\widetilde{\gamma}(x))\,\mathrm{d}^{2}x$
(22)

will be similarly invariant. The field line winding $L_{\bf v}(\gamma)$, with
$w\equiv 1$ is the simplest such invariant, and is already sufficient to
completely determine the field line topology of ${\bf v}$. Nevertheless, there
may be situations in which other weightings are of physical interest.

The most notable example is field line helicity, which is defined for the
subset of divergence-free braided vector fields. Typical examples are magnetic
fields (in magnetohydrodynamics) or vorticity fields. In this section, we
shall denote such fields as ${\bf b}$ rather than ${\bf v}$. The field line
helicity is the particular $W_{\bf v}$ invariant given by

$A_{\bf
b}(\gamma)=\int_{D_{0}}J_{0}(x)b_{z}(x)\,\mathcal{L}(\gamma,\widetilde{\gamma}(y))\,\mathrm{d}^{2}x,$
(23)

where $J_{0}$ is the Jacobian of the map $F^{-1}:S_{0}\to D_{0}$, and $b_{z}$
is the normal component of ${\bf b}$ on $S_{0}$. When
$\mathcal{M}=\mathcal{C}$, we have $J_{0}\equiv 1$ so

$A_{\bf
b}(\gamma)=\int_{D_{0}}b_{z}(x)\,\mathcal{L}(\gamma,\widetilde{\gamma}(x))\,\mathrm{d}^{2}x.$
(24)

For this case, ? showed that (24) is equivalent to the original definition of
field line helicity, $A_{\bf b}(\gamma)=\int_{\gamma}{\bf
a}\cdot\,\mathrm{d}{\bf l}$ (where ${\bf b}=\nabla\times{\bf a}$) provided
that the vector potential ${\bf a}$ is in the so-called winding gauge

${\bf a}(x^{1},x^{2},z)=\frac{1}{2\pi}\int_{D_{z}}\frac{{\bf
b}(\tilde{x}^{1},\tilde{x}^{2},z)\times(x^{1}-\tilde{x}^{1},x^{2}-\tilde{x}^{2},0)}{(x^{1}-\tilde{x}^{1})^{2}+(x^{2}-\tilde{x}^{2})^{2}}\,\mathrm{d}^{2}\tilde{x}.$
(25)

Indeed, they used the geometrical interpretation in terms of winding numbers
to motivate this particular choice of gauge.

For more general $\mathcal{M}$, we can similarly relate (23) to a particular
choice of vector potential. It is convenient to think of ${\bf b}$ as a
differential 2-form $\beta=\mathrm{i}_{\bf b}\mathrm{vol}^{3}$ on
$\mathcal{M}$, where $\mathrm{i}$ is the interior product and
$\mathrm{vol}^{3}$ the volume form (?). Then pushing forward the integral from
$D_{0}$ to $S_{0}$, and using equation (3.13) of ?, we have

$\displaystyle A_{\bf b}(\gamma)$
$\displaystyle=\int_{S_{0}}\beta_{z}(\tilde{y})\,\mathcal{L}(\gamma,\widetilde{\gamma}(\tilde{y}))\,\mathrm{d}\tilde{y}^{1}\wedge\mathrm{d}\tilde{y}^{2}$
(26)
$\displaystyle=\int_{0}^{1}\int_{S_{0}}\beta_{z}(\tilde{y})\frac{\mathrm{d}}{\mathrm{d}z}\Theta(\gamma,\widetilde{\gamma}(\tilde{y}),z)\,\mathrm{d}\tilde{y}^{1}\wedge\mathrm{d}\tilde{y}^{2}\,\mathrm{d}z$
(27)

Now, the divergence free property of ${\bf b}$ means that the form
$\beta_{z}\,\mathrm{d}y^{1}\wedge\mathrm{d}y^{2}$ is conserved along field
lines, so we can write

$A_{\bf
b}(\gamma)=\int_{0}^{1}\int_{S_{z}}\beta_{z}(\tilde{y})\frac{\mathrm{d}}{\mathrm{d}z}\Theta(\gamma,\widetilde{\gamma}(\tilde{y}),z)\,\mathrm{d}\tilde{y}^{1}\wedge\mathrm{d}\tilde{y}^{2}\,\mathrm{d}z.$
(28)

Using (4) and the fact that
$\mathrm{d}\gamma_{i}/\mathrm{d}z=\beta_{i}/\beta_{z}$, we find that

$A_{\bf b}(\gamma)=\int_{\gamma}\alpha,$ (29)

where the one-form
$\alpha=\alpha_{1}\mathrm{d}y^{1}+\alpha_{2}\mathrm{d}y^{2}+\alpha_{z}\mathrm{d}z$
has components

$\displaystyle\alpha_{1}(y)$
$\displaystyle=-\frac{\beta_{z}(\tilde{y})(y^{2}-\tilde{y}^{2})}{(y^{1}-\tilde{y}^{1})^{2}+(y^{2}-\tilde{y}^{2})^{2}}\,\mathrm{d}\tilde{y}^{1}\wedge\mathrm{d}\tilde{y}^{2},$
(30) $\displaystyle\alpha_{2}(y)$
$\displaystyle=\frac{\beta_{z}(\tilde{y})(y^{1}-\tilde{y}^{1})}{(y^{1}-\tilde{y}^{1})^{2}+(y^{2}-\tilde{y}^{2})^{2}}\,\mathrm{d}\tilde{y}^{1}\wedge\mathrm{d}\tilde{y}^{2},$
(31) $\displaystyle\alpha_{z}(y)$
$\displaystyle=\frac{\beta_{1}(\tilde{y})(y^{2}-\tilde{y}^{2})-\beta_{2}(y^{1}-\tilde{y}^{1})}{(y^{1}-\tilde{y}^{1})^{2}+(y^{2}-\tilde{y}^{2})^{2}}\,\mathrm{d}\tilde{y}^{1}\wedge\mathrm{d}\tilde{y}^{2}.$
(32)

This one-form corresponds to a vector potential ${\bf a}$ on $\mathcal{M}$
which is the generalisation of the winding-gauge vector potential (25).

Notice that $A_{\bf b}$ contains no additional topological information
compared to $L_{\bf b}$; the difference is that it weights each field line
differently according to $b_{z}$. This weighting has physical meaning because
$\beta_{z}\,\mathrm{d}y^{1}\wedge\mathrm{d}y^{2}$ is preserved along field
lines – thus the weighting is a property of the field line as a whole, not
just an arbitrary function of the endpoints on $S_{0}$. An interesting
consequence of the additional weighting is that, although $A_{\bf b}$ is
invariant under end-vanishing isotopies of ${\bf b}$ like $L_{\bf b}$, it may
behave differently from $L_{\bf b}$ if the domain $\mathcal{M}$ changes shape.
This remains to be investigated further.

##  7 Conclusions and comparison to other work

The primary development in this work (Theorem 1) is the proof of a topological
classification theorem for braided vector fields on tubular subdomains, which
are homeomorphic to the cylinder as defined in Section 4. Specifically, the
field lines of a given braided field may be continuously deformed into those
of another other field – without moving their end-points – if and only if the
two fields have the same distribution of field line winding, $L_{\bf v}$.
There is a small caveat: this invariant cannot distinguish overall full
(integer) rotations. These could be detected solely from the field line
mappings on the side boundary, and are unlikely to be problematic in practical
applications. Our result strengthens an earlier theorem of ? and extends it to
a broader class of vector fields; specifically to include those which have
divergence. The measure $L_{\bf v}$ has recently been applied to study
reconnection in a laboratory plasma experiment (?); here we have formalised
the underpinning result and extended it to tubular sub-domains that cannot be
foliated by parallel planes, including tubular domains foliated by non-
Euclidean surfaces.

The unique definition of $L_{\bf v}$ relies on the identification of a least
distorted vector field ${\bf u}$, defined in Section 4.1, whose integral
curves are mapped to straight lines on the cylinder (on which the field line
winding is defined). Our choice of ${\bf u}$ and its mapping to the cylinder
have three critical properties:

  1. 1.

The vector field ${\bf u}$ minimises the Dirichlet (distortion) energy amongst
all braided vector fields on a given tubular domain $\mathcal{M}$ with a
specified flux, i.e. it minimises

$E=\int_{\mathcal{M}}{\bf v}\cdot{\bf v}\,\rmd V,\quad\textrm{for specified
}\int_{S_{0}}\hat{\bf z}\cdot{\bf v}\,\rmd{S}.$ (33)

The latter constraint relates to the fact the energy varies with flux for a
fixed topology, whereas we are interested purely in finding the vector field
of least entangled topology.

  2. 2.

In this mapping, all field lines of ${\bf u}$ have zero pairwise winding
number, so the winding measure $L_{\bf v}$ is everywhere zero for ${\bf u}$
itself.

  3. 3.

To determine which field line of ${\bf u}$ maps to which vertical line on
$\mathcal{C}$, we use harmonic coordinates so as to avoid the erroneous
measurement of entanglement such as in Figure 5.

The combination of these three properties means that a non-trivial
distribution $L_{\bf v}$ necessarily indicates non-trivial internal
entanglement of the original vector field ${\bf v}$. This entanglement may be
seen as increasing the energy above the minimum level enforced by the
distorted global shape of the domain $\mathcal{M}$ itself. The only comparable
construction to this of which we are aware is the so-called Mermin-Ho basis,
which is defined on a pre-determined Riemannian manifold. Using the notion of
parallel transport with respect to a preferred vector field (not specifically
assigned in the theory) one can define a local basis on the manifold from
which any local deviations of an arbitrary vector from this field are
quantified (?, ?). Our work differs in that we define and justify the choice
of a preferred vector field. Moreover, this justification is based on a
globally defined notion (the winding rate
$\rmd\Theta(\gamma,\widetilde{\gamma},z)/\rmd z$), rather than a local notion
(parallel transport). The Godbillion-Vey invariant see e.g. (?), can be
related to the helicity of the integral curves of the foliation of a closed
domain so is not directly comparable here (see (?) for a discussion of its
applicability in MHD contexts).

The invariant $L_{\bf v}$ used in our classification theorem provides a
fundamental description of the field line topology, in the sense that it
depends only on the integral curves of the vector field and not on its
magnitude or flux. This is in contrast to the widely used (field line)
helicity invariants for divergence-free fields (?, ?, ?, ?, ?), which do
depend on the flux. Indeed, we showed in Section 6 that the field line
helicity – effectively, the helicity of an infinitesimal flux tube surrounding
a single field line – takes a similar integral form to $L_{\bf v}$ but with an
additional weighting of flux. As such, it mixes a purely topological quantity
(the pairwise winding) with a physical quantity (the flux). This indicated the
invariance of helicity under Euler flows (ideal motions) which vanish at the
domain boundaries requires not only the invariance of
$\mathcal{L}(\gamma,\widetilde{\gamma})$ but also the divergence-free
condition, which ensures preservation of flux. Since the flux can change
without altering the field line topology, $L_{\bf v}$ is revealed as a more
fundamental quantity than the (field line) helicity.

##  Appendix A Validity of the foliation defined by ${\bf u}$

(a)  (b)

Figure 8: Possible local structures of the constant-$\phi$ surfaces around
null points of harmonic fields ${\bf u}=\nabla\phi$. There are two
possibilities, depending on whether the null point is isolated (a) or part of
a null line (b).

In this appendix we show that the least-distorted harmonic field ${\bf u}$
defines a valid foliation of the domain $\mathcal{M}$. In particular, we will
show that ${\bf u}\neq{\bf 0}$ everywhere on the interior of $\mathcal{M}$, so
that there is a unique surface $S_{z}$ through all interior points and the
winding integrals are well defined. It is possible to have points with ${\bf
u}={\bf 0}$ on the boundary $\partial M$, but these do not prevent the
existence of the foliation.

To rule out null points of our harmonic field ${\bf u}=\nabla\phi$, the idea
is to consider the possible structures of constant-$\phi$ surfaces that could
result from such points, and show that such surfaces cannot exist. If $p$ is a
hypothetical point where ${\bf u}(p)={\bf 0}$, then the local linearisation of
${\bf u}$ around $p$ is constrained by the fact that ${\bf u}$ is harmonic.
Firstly, the tensor ${\bf u}$ is symmetric due to $\nabla\times{\bf u}={\bf
0}$, meaning that it has real eigenvalues $\lambda_{1}$, $\lambda_{2}$,
$\lambda_{3}$. Secondly, it is traceless due to $\nabla\cdot{\bf u}=0$, so
that $\lambda_{1}+\lambda_{2}+\lambda_{3}=0$. There are two possible
structures around $p$, shown in Figure 8. If all three eigenvalues are non-
zero, then in suitable local coordinates we may write

${\bf u}=\nabla\phi\approx\pm\left(x,y,-2z\right),\mbox{ so
}\phi\approx\pm\frac{1}{2}\left(x^{2}+y^{2}-2z^{2}\right).$ (34)

Then the surface $\phi=0$ has the topology of two cones whose tips meet at $p$
(Figure 8a), separating two-sheet hyperboloids with $\phi>0$ and one-sheet
hyperboloids with $\phi<0$. Alternatively, if one eigenvalue is $0$ then $p$
is part of a null line (Figure 8b) and, again in suitable local coordinates,
we may write

${\bf u}=\nabla\phi=\pm\left(x,0,-z\right),\mbox{ so
}\phi=\pm\left(x^{2}-z^{2}\right).$ (35)

The $\phi=0$ surface is then a pair of intersecting planes. The null line
could either be a closed loop or could intersect the boundary
$\partial\mathcal{M}$. It is not possible to have a null surface, nor for all
three eigenvalues to vanish at $p$ (unless ${\bf u}={\bf 0}$ on the whole
domain).

To rule out null points and lines, we must consider how the $\phi$ surfaces
extend outside the neighbourhoods of these structures. To do so, we use the
fact that $\phi$ defines a continuous foliation of $\mathcal{M}$ everywhere
outside the null set $\mathcal{N}=\\{p\in\mathcal{M}\,|\,{\bf u}(p)={\bf
0}\\}$. This follows from the Frobenius theorem. So the $\phi=0$ surfaces from
null points or lines extend smoothly away from these locations and can end
only on $\mathcal{N}$ or on the boundary $\partial\mathcal{M}$. In fact, they
can only intersect (transversely) the side boundary $S_{s}$, owing to the
boundary conditions of constant $\phi$ on $S_{0}$ and $S_{1}$.

(a)  (b)  (c)  (d)

Figure 9: Global extensions of the $\phi=0$ surfaces from (a) an isolated null
point, (b) two isolated null points, and (c) a null point and a null line. In
(a), the surfaces intersect the side boundary $S_{s}$, enclosing a volume
$\mathcal{V}$. In (d) the $\phi=0$ surface intersects $S_{0}$ and $S_{1}$
transversely, but this is forbidden by the boundary conditions
$\phi=\textrm{constant}$ on $S_{0}$ and $S_{1}$.

Consider first a single null point in the interior, as in Figure 9(a). Each of
the two $\phi=0$ surfaces must either be closed or end on $S_{s}$ (as shown in
Figure 9a). To show that neither is possible, let $\hat{\bf n}$ be the outward
normal from the indicated volume $\mathcal{V}$. Then observe that $\hat{\bf
n}\cdot{\bf u}$ is non-zero with the same sign everywhere on these two
surfaces where they bound $\mathcal{V}$, because $\phi<0$ on one side and
$\phi>0$ on the other. Since the two surfaces intersect $S_{s}$, and since
$\hat{\bf n}\cdot{\bf u}=0$ on $S_{s}$, there must be a net flux in/out of
$\mathcal{V}$, contradicting $\nabla\cdot{\bf u}=0$. If the surface were
closed instead of intersecting $S_{s}$, we would similarly obtain a
contradiction. So a single null point in the interior cannot exist.

This argument extends to rule out multiple null points or lines in the
interior (of which there can only be a countable number). Tracing the $\phi=0$
surfaces from some null point or line, we can extend these surfaces, possibly
incorporating further points in $\mathcal{N}$ – as depicted in Figure 9(b) or
(c) – until the composite surfaces either close on themselves or intersect
$S_{s}$. Either way we can find a finite volume enclosed with a non-zero flux
of ${\bf u}$ through its boundary, which contradicts $\nabla\cdot{\bf u}=0$.

Finally we point out that since the boundary is Lipschitz continuous it is
possible to embed the conic structure of a null surface $\phi=0$ in the
boundary. As an example this would be necessary on the bounding curves
$\partial S_{0/1}$ for an hourglass domain morphology, owing to the boundary
conditions on ${\bf u}$.

##  Appendix B Derivation of Equation (18)

To evaluate the left-hand side of (18) at the point $x_{0}=(r,\theta)$,
consider the components of $\mathrm{d}_{\perp}\Theta$ in polar coordinates
$(r,\theta)$. Differentiating the expression (2) for $\Theta(x_{0},x,z)$ gives

$\displaystyle\frac{\partial\Theta}{\partial r}$
$\displaystyle=\frac{(x_{0}^{1}-x^{1})\sin\theta-(x_{0}^{2}-x^{2})\cos\theta}{(x_{0}^{1}-x^{1})^{2}+(x_{0}^{2}-x^{2})^{2}},$
(36) $\displaystyle\frac{\partial\Theta}{\partial\theta}$
$\displaystyle=r\frac{(x_{0}^{1}-x^{1})\cos\theta+(x_{0}^{2}-x^{2})\sin\theta}{(x_{0}^{1}-x^{1})^{2}+(x_{0}^{2}-x^{2})^{2}}.$
(37)

To evaluate the integrals, change from $(x^{1},x^{2})$ to polar coordinates
centered at $x_{0}$,

$x^{1}=\rho\cos\varphi+x_{0}^{1},\qquad x^{2}=\rho\sin\varphi+x_{0}^{2}.$ (38)

In these coordinates,

$\displaystyle\frac{\partial\Theta}{\partial
r}=\frac{\sin(\varphi-\theta)}{\rho},\qquad\frac{\partial\Theta}{\partial\theta}=-r\frac{\cos(\varphi-\theta)}{\rho}.$
(39)

Integrating with respect to $(\rho,\varphi)$ we must account for the fact that
the $\rho$ limit of integration is a function of $\varphi$. Substituting
expressions (38) into the equation $(x^{1})^{2}+(x^{2})^{2}=1$ for the
boundary leads to the integration limit

$\displaystyle\rho(\varphi)=-r\cos(\varphi-\theta)+\sqrt{r^{2}\cos^{2}(\varphi-\theta)-r^{2}-1}.$
(40)

Using the phase behaviour of $\sin$ and $\cos$, we may then evaluate

$\displaystyle\frac{1}{2\pi}\int_{D_{0}}\frac{\partial\Theta}{\partial
r}\,\rho\mathrm{d}\rho\mathrm{d}\varphi=\frac{1}{2\pi}\int_{0}^{2\pi}\sin(\varphi-\theta)\rho(\phi)\,\mathrm{d}\varphi=0$
(41)

and

$\displaystyle\frac{1}{2\pi}\int_{D_{0}}\frac{\partial\Theta}{\partial\theta}\,\rho\mathrm{d}\rho\mathrm{d}\varphi=\frac{r}{2\pi}\int_{0}^{2\pi}\cos^{2}(\varphi-\theta)\,\mathrm{d}\varphi=\frac{r}{2}.$
(42)

## References

## References

  * [1]
  * [2] [] Aly J J 2018 Fluid Dynamics Research 50(1), 011408.
  * [3]
  * [4] [] Aref H 1984 Journal of fluid mechanics 143, 1–21.
  * [5]
  * [6] [] Arnold V I & Khesin B A 1992 Annual review of fluid mechanics 24(1), 145–166.
  * [7]
  * [8] [] Arnolʹd V I 2012 Geometrical methods in the theory of ordinary differential equations Vol. 250 Springer Science & Business Media.
  * [9]
  * [10] [] Bao J & Yang Q 2010 Nonlinear Dynamics 61(4), 773–781.
  * [11]
  * [12] [] Berger M 1988 Astronomy and Astrophysics 201, 355–361.
  * [13]
  * [14] [] Berger M A & Field G B 1984 Journal of Fluid Mechanics 147, 133–148.
  * [15]
  * [16] [] Berger M A & Prior C 2006 Journal of Physics A: Mathematical and General 39(26), 8321.
  * [17]
  * [18] [] Birman J S 2016 Braids, Links, and Mapping Class Groups.(AM-82) Vol. 82 Princeton University Press.
  * [19]
  * [20] [] Cantarella J & Parsley J 2010 Journal of Geometry and Physics 60(9), 1127–1155.
  * [21]
  * [22] [] Chen L 2009 University of California at Irvine .
  * [23]
  * [24] [] Childress S & Gilbert A D 1995 Stretch, twist, fold: the fast dynamo Vol. 37 Springer Science & Business Media.
  * [25]
  * [26] [] Contopoulos I, Christodoulou D M, Kazanas D & Gabuzda D C 2009 The Astrophysical Journal Letters 702(2), L148.
  * [27]
  * [28] [] Dalmasse K, Pariat E, Démoulin P & Aulanier G 2014 Solar Physics 289(1), 107–136.
  * [29]
  * [30] [] Del Sordo F, Candelaresi S & Brandenburg A 2010 Physical Review E 81(3), 036401.
  * [31]
  * [32] [] DeTurck D M & Kazdan J L 1981 in ‘Annales Scientifiques de l’École Normale Supérieure’ Vol. 14 Elsevier pp. 249–260.
  * [33]
  * [34] [] Forth S, Deufel C, Sheinin M Y, Daniels B, Sethna J P & Wang M D 2008 Physical Review Letters 100(14), 148301.
  * [35]
  * [36] [] Frankel T 2011 The geometry of physics: an introduction Cambridge University Press.
  * [37]
  * [38] [] Gol’dshtein V, Mitrea I & Mitrea M 2011 Journal of Mathematical Sciences 172(3), 347–400.
  * [39]
  * [40] [] Gopalswamy N, Akiyama S, Yashiro S & Xie H 2017 Proceedings of the International Astronomical Union 13(S335), 258–262.
  * [41]
  * [42] [] Goriely A & Tabor M 1998 Physical Review Letters 80(7), 1564.
  * [43]
  * [44] [] Grason G M 2009 Physical Review E 79(4), 041919.
  * [45]
  * [46] [] Hatcher A 2002 Algebraic topology Cambridge University Press.
  * [47]
  * [48] [] Heinz S & Begelman M C 2000 The Astrophysical Journal 535(1), 104.
  * [49]
  * [50] [] Kamien R D 2002 Reviews of Modern physics 74(4), 953.
  * [51]
  * [52] [] Komissarov S 1999 Monthly Notices of the Royal Astronomical Society 308(4), 1069–1076.
  * [53]
  * [54] [] Longcope D W & Malanushenko A 2008 The Astrophysical Journal 674(2), 1130–1143.
  * [55]
  * [56] [] Mann K 2013 New York J. Math 19, 583–596.
  * [57]
  * [58] [] Mermin N & Ho T L 1976 Physical Review Letters 36(11), 594.
  * [59]
  * [60] [] Moffatt H K 1969 Journal of Fluid Mechanics 35(1), 117–129.
  * [61]
  * [62] [] Moffatt H K & Proctor M R 1985 Journal of Fluid Mechanics 154, 493–507.
  * [63]
  * [64] [] Moraitis K, Pariat E, Valori G & Dalmasse K 2019 Astronomy & Astrophysics 624, A51.
  * [65]
  * [66] [] Moreau J 1961 C. R. Acad. Sci. Paris 252, 2810–2812.
  * [67]
  * [68] [] Mullinger A & Johnson R 1980 Journal of cell science 46(1), 61–86.
  * [69]
  * [70] [] Nordlund A, Brandenburg A, Jennings R L, Rieutord M, Ruokolainen J, Stein R F & Tuominen I 1992 The Astrophysical Journal 392, 647–652.
  * [71]
  * [72] [] Pohl W F 1968 Journal of Mathematics and Mechanics 17(10), 975–985.
  * [73]
  * [74] [] Prior C B & Neukirch S 2016 Journal of Physics A: Mathematical and Theoretical 49(21), 215201.
  * [75]
  * [76] [] Prior C & Gourgouliatos K N 2019 Astronomy & Astrophysics 622, A122.
  * [77]
  * [78] [] Prior C & MacTaggart D 2016 Geophysical & Astrophysical Fluid Dynamics 110(5), 432–457.
  * [79]
  * [80] [] Prior C & MacTaggart D 2019 Journal of Plasma Physics 85(2).
  * [81]
  * [82] [] Prior C & Yeates A 2014 The Astrophysical Journal 787(2), 100.
  * [83]
  * [84] [] Prior C & Yeates A 2016 Astronomy & Astrophysics 591, A16.
  * [85]
  * [86] [] Prior C & Yeates A R 2018 Physical Review E 98(1), 013204.
  * [87]
  * [88] [] Romano P, Pariat E, Sicari M & Zuccarello F 2011 Astronomy & Astrophysics 525, A13.
  * [89]
  * [90] [] Russell A J, Yeates A R, Hornig G & Wilmot-Smith A L 2015 Physics of Plasmas 22(3), 032106.
  * [91]
  * [92] [] Rust D M & Kumar A 1996 The Astrophysical Journal Letters 464(2), L199.
  * [93]
  * [94] [] Smale S 1959 Proceedings of the American Mathematical Society 10(4), 621–626.
  * [95]
  * [96] [] Starostin E & Van der Heijden G 2014 Journal of the Mechanics and Physics of Solids 64, 83–132.
  * [97]
  * [98] [] Taylor J B 1974 Physical Review Letters 33(19), 1139–1141.
  * [99]
  * [100] [] Thompson J T, van der Heijden G M & Neukirch S 2002 in ‘Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences’ Vol. 458 The Royal Society pp. 959–985.
  * [101]
  * [102] [] Török T & Kliem B 2005 The Astrophysical Journal Letters 630(1), L97.
  * [103]
  * [104] [] Webb G M, Dasgupta B, McKenzie J F, Hu Q & Zank G P 2014 Journal of Physics A: Mathematical and Theoretical 47(9), 095501.
  * [105]
  * [106] [] Wilmot-Smith A 2015 Phil. Trans. R. Soc. A 373(2042), 20140265.
  * [107]
  * [108] [] Wilmot-Smith A L, Pontin D I, Yeates A R & Hornig G 2011 Astronomy & Astrophysics 536, A67.
  * [109]
  * [110] [] Woltjer L 1958 PNAS 44(6), 489–491.
  * [111]
  * [112] [] Yeates A & Hornig G 2013 Physics of Plasmas 20(1), 012102.
  * [113]
  * [114] [] Yeates A & Hornig G 2014 in ‘Journal of Physics: Conference Series’ Vol. 544 IOP Publishing p. 012002.
  * [115]
  * [116] [] Yeates A & Hornig G 2016 Astronomy & Astrophysics 594, A98.
  * [117]
  * [118] [] Yeates A R & Page M H 2018 Journal of Plasma Physics 84(6).
  * [119]


# Comments on: ”Almost All Carbon/Oxygen White Dwarfs Can Support Double
Detonations”

Noam Soker  0000-0003-0375-8987 Department of Physics, Technion, Haifa,
3200003, Israel<EMAIL_ADDRESS>
###### Abstract

I critically review some claims in the paper “Almost All Carbon/Oxygen White
Dwarfs Can Support Double Detonations” (arXiv:2405.19417). The claim of that
paper that the community converges on a leading scenario of type Ia supernovae
(SNe Ia), the double detonation scenario, is wrong, as hundreds of papers in
recent years study five and more different SN Ia scenarios and their channels.
Moreover, the finding by that paper that the double detonation scenario with
the explosion of the secondary white dwarf (WD; the mass-donor WD) is common,
i.e., the triple-detonation channel and the quadruple-detonation sub-channel,
implies highly no-spherical explosions. The highly non-spherical explosions
contradict the morphologies of many SN Ia remnants. I find that the results of
that paper strengthen the claim that the double detonation scenario (with its
channels) might account for a non-negligible fraction of peculiar SNe Ia but
only for a very small fraction (or non at all) of normal SNe Ia.

Type Ia supernovae – Supernova remnants

##  1 Introduction

The community is far from a consensus on the leading scenarios for type Ia
supernovae (SNe Ia), as is evident from reviews and classifications in the
last decade (Maoz et al. 2014; Maeda, & Terada 2016; Hoeflich 2017; Livio &
Mazzali 2018; Soker 2018, 2019; Wang 2018; Jha et al. 2019; Ruiz-Lapuente
2019; Ruiter 2020; Aleo et al. 2023; Liu, Röpke, & Han 2023; Vinkó, Szalai, &
Könyves-Tóth 2023). There are five, six, or even more different SN Ia
scenarios (e.g., Liu, Röpke, & Han 2023; Soker 2024a), with even more channels
within some scenarios. The community does not converge on one-leading, two-
leading, or even three-leading scenarios. Hundreds of papers in the last five
years intensively discussed five scenarios and their channels for both normal
and peculiar SNe Ia. Liu, Röpke, & Han (2023) discuss in detail the single
degenerate (SD) scenario while Soker (2024a) argue for the core degenerate
(CD) scenario and the double degenerate scenario with the merger to explosion
delay (MED) time (the DD-MED scenario) for normal SNe Ia, namely, the group of
lonely white dwarf (WD) scenarios; other scenarios might account for a
fraction of peculiar SNe Ia. Some researchers prefer the double detonation
(DDet) scenario, where a WD accretes helium from a companion, and the ignition
of the helium (the first detonation) detonates the CO interior of the WD (the
second detonation).

Despite this obvious state of SN Ia scenarios and channels, Shen, Boos, &
Townsley (2024) open the Abstract of their paper “Almost All Carbon/Oxygen
White Dwarfs Can Support Double Detonations” with the statement “Double
detonations of sub-Chandrasekhar-mass white dwarfs (WDs) in unstably mass-
transferring double WD binaries have become a leading contender for explaining
most, if not all, Type Ia supernovae.” They then open their Introduction with
“The identity of Type Ia supernova (SN Ia) progenitors remains an unsolved
mystery (see, e.g., Liu, Röpke, & Han 2023 for a recent review), but recent
work appears to be converging towards a solution. There is now mounting
evidence that most, if not all, SNe Ia are the result of detonations of sub-
Chandrasekhar-mass white dwarfs (WDs) in double WD systems.” Concerning the
highly diverse view of the community, the above opening sentences deceive the
unfamiliar readers. It is fine to present one view and preference, but it
becomes unethical to present one preference as the view that the diverse
community converges on. More importantly, I find the new results that Shen,
Boos, & Townsley (2024) present to add problems to the DDet scenario as one of
the scenarios for normal SNe Ia.

##  2 The challenges of the double detonation scenario

There are two strong arguments against the DDet scenario (e.g., review by
Soker 2019) : (1) The DDet scenario predicts a highly non-spherical explosion,
contrary to the morphology of a large fraction of SNe Ia remnants (SNRs Ia)
that are globally spherical or axially symmetric. Some SNRs Ia do possess a
large departure from sphericity. (2) The surviving helium-donor companion is
not observed in SNRs Ia, and the observed number of surviving hypervelocity
WDs (the mass donor companions) is much below the number of normal SNe Ia
(Igoshev, Perets, & Hallakoun, 2023). By arguing that most explosions in the
DDet scenario are actually the triple detonation scenario (as coined and
simulated by Papish et al. 2015), Shen, Boos, & Townsley (2024) solve the
second problem and explain the nondetection of surviving companions. In the
triple detonation channel of the DDet scenario, helium detonation on the
surface of one WD triggers a CO detonation in the same WD. The exploding WD
ejecta triggers the explosion of the mass-donor WD companion, either the
detonation of helium of a helium WD or helium shell detonation that triggers a
CO detonation of the secondary WD in the quadruple detonation sub-channel of
the triple channel (a term coined by Tanikawa et al. 2019). Namely, both WDs
explode. (In the words of Shen, Boos, & Townsley 2024: “Our work suggests that
a majority of SNe Ia may arise from the double detonations of both WDs in a
merging double WD binary.”). The triple and quadruple detonation channels
result in an explosion similar to the DD scenario without a MED time.

As I now argue, solving the problem of the surviving companion with the triple
detonation channel of the DDet scenario or with the quadruple detonation sub-
channel, Shen, Boos, & Townsley (2024) make the predicted asymmetrical
explosion problem more severe. In the spirit of their paper, Shen, Boos, &
Townsley (2024) completely ignore the most severe problem of the DDet
scenario, namely, that the predicted asymmetrical explosion contradicts
observations of many SNRs Ia (see review by Soker 2019). Polarization
measurements of SNe Ia show that many have global spherical symmetry, although
there are non-spherical metal distributions (e.g., Cikota et al. 2019;
Hoeflich et al. 2023), which is also compatible with deflagration to the
detonation of WDs near Chandrasekhar mass, as in the single-degenerate,
double-degenerate with MED time, and the core degenerate.

The two-dimensional hydrodynamical simulations by Papish et al. (2015)
demonstrate the non-spherical ejecta of the DDet scenario, including the
triple detonation channel. In the case of the DDet scenario, when the
secondary (mass donor) WD does not explode, it casts an ‘ejecta shadow’ behind
it. This leads to a supernova remnant that has a clear departure from
spherical symmetry on a large scale. The SNR will appear spherical only if
observed along or near the axis connecting the two WDs. Even this is not
accurate, as the two WDs orbit each other, which will lead to deviation from a
symmetry axis. Papish et al. (2015) conducted two-dimensional simulations and,
therefore, could not follow the orbital motion of the two stars. In a simple,
non-hydrodynamical, calculation, Braudo & Soker (2024) demonstrate that the
orbital motion during the explosion adds some non-spherical structure to the
inner ejecta in the DDet scenario. The effect of the ejecta shadow is seen in
more detail and clearer in three-dimensional simulations (e.g., Tanikawa,
Nomoto, & Nakasato 2018; Tanikawa et al. 2019)

Papish et al. (2015) find that the triple-detonation channel of the DDet
scenario leads to an even larger departure of the ejecta from global spherical
symmetry than the DDet scenario. I demonstrate this in figure 1 of the
temperature and velocity maps in the meridional plane of a simulation by
Papish et al. (2015) of the triple-detonation channel. The simulation’s high-
temperature zone (red) is the material originating in the secondary WD. The
temperature and velocity maps show a large-scale departure from spherical
symmetry in the non-radial velocities and the unequal distribution of mass on
the two sides (one of the primary and one of the secondary WD). Tanikawa et
al. (2019), in their three-dimensional simulations, follow the outflow in much
more detail and to larger distances, one case of the triple-detonation channel
and one of the quadruple-detonation sub-channel. Their density and composition
maps present prominent large-scale departures from any symmetry, both in
density and composition.

Figure 1: Temperature and velocity maps in the meridional plane at six times;
figure built from a simulation by Papish et al. (2015) of the triple
detonation channel of the DDet scenario. In this case, there is a He WD of
$0.43M_{\odot}$, the blue circle at the upper-left panel, at an initial
distance of its center to the center of the explosion of $0.045R_{\odot}$.
Only one-quarter (not in full) of the meridional plane is shown. At
$t=2{~{}\rm s}$, the helium in the secondary WD is ignited (the high-
temperature zone), and the secondary WD is exploded. The velocity at each grid
point is proportional to the arrow length, with the inset showing a length
corresponding to a velocity of $10,000{~{}\rm km}{~{}\rm s}^{-1}$. Note the
highly non-spherical explosion at the last time (bottom-left panel).

The morphological expectations of large departure from axial symmetry and
spherical symmetry of the triple-detonation channel and the quadruple-
detonation sub-channel are in contradiction with the large fraction of SNRs Ia
that possess global spherical structure, e.g., the Tycho SNR and SN 1006, and
SNRs Ia that possess global axially-symmetric structures, e.g., SNR G1.9+0.3
(image by, e.g., Enokiya et al. 2023). There are some SNRs Ia with large and
global departures from any symmetry. I cannot rule out the DDet scenario for
these SNRs Ia. It is upon the supporters of the DDet scenario (with the
triple-detonation channel and quadruple-detonation sub-channel) to explain how
their high non-spherical explosion can account for at least one of the
observed SNR Ia. An example is the way the supporters of the core degenerate
scenario account for the point-symmetric morphology of SNR G1.9+0.3 as a
spherical explosion inside a planetary nebula (SNIP; Soker 2024b).

##  3 Summary

I conclude this comment on the recent paper by Shen, Boos, & Townsley (2024)
by reiterating my view that the DDet scenario and its triple-detonation
channel and quadruple-detonation sub-channel can account at most for a small
fraction of normal SNe Ia. The new results of Shen, Boos, & Townsley (2024)
reduce this fraction even more. Namely, I argue that the DDet scenario is not
a dominant scenario of normal SNe Ia. It likely contributes to peculiar SNe
Ia. My view is not in the consensus, as there is no consensus, and therefore,
the claim by Shen, Boos, & Townsley (2024) that the community converges (in
2024) on a leading scenario is wrong. A consensus on the SN Ia leading
scenario is for the future.

## References

  * Aleo et al. (2023) Aleo P. D., Malanchev K., Sharief S., Jones D. O., Narayan G., Foley R. J., Villar V. A., et al., 2023, ApJS, 266, 9.
  * Braudo & Soker (2024) Braudo J., Soker N., 2024, OJAp, 7, 7.
  * Cikota et al. (2019) Cikota A., Patat F., Wang L., Wheeler J. C., Bulla M., Baade D., Höflich P., et al., 2019, MNRAS, 490, 578.
  * Enokiya et al. (2023) Enokiya R., Sano H., Filipović M. D., Alsaberi R. Z. E., Inoue T., Oka T., 2023, PASJ, 75, 970. doi:10.1093/pasj/psad054
  * Hoeflich (2017) Hoeflich P., 2017, in Handbook of Supernovae, Springer International Publishing AG, 2017, p. 1151
  * Hoeflich et al. (2023) Hoeflich P., Yang Y., Baade D., Cikota A., Maund J. R., Mishra D., Patat F., et al., 2023, MNRAS, 520, 560.
  * Igoshev, Perets, & Hallakoun (2023) Igoshev A. P., Perets H., Hallakoun N., 2023, MNRAS, 518, 6223.
  * Jha et al. (2019) Jha, S. W., Maguire, K., & Sullivan, M. 2019, Nature Astronomy, 3, 706
  * Liu, Röpke, & Han (2023) Liu Z.-W., Röpke F. K., Han Z., 2023, RAA, 23, 082001.
  * Livio & Mazzali (2018) Livio, M., & Mazzali, P. 2018, Physics Reports, 736, 1
  * Maeda, & Terada (2016) Maeda, K., & Terada, Y. 2016, International Journal of Modern Physics D, 25, 1630024
  * Maoz et al. (2014) Maoz, D., Mannucci, F., & Nelemans, G. 2014, ARA&A, 52, 107
  * Papish et al. (2015) Papish O., Soker N., García-Berro E., Aznar-Siguán G., 2015, MNRAS, 449, 942. doi:10.1093/mnras/stv337
  * Ruiter (2020) Ruiter A. J., 2020, IAUS, 357, 1.
  * Ruiz-Lapuente (2019) Ruiz-Lapuente, P. 2019, New A Rev., 85, 101523
  * Shen, Boos, & Townsley (2024) Shen K. J., Boos S. J., Townsley D. M., 2024, arXiv, arXiv:2405.19417. doi:10.48550/arXiv.2405.19417
  * Soker (2018) Soker N., 2018, SCPMA, 61, 49502.
  * Soker (2019) Soker N., 2019, NewAR, 87, 101535.
  * Soker (2024a) Soker N., 2024a, OJAp, 7, 31.
  * Soker (2024b) Soker N., 2024b, RAA, 24, 015012.
  * Tanikawa, Nomoto, & Nakasato (2018) Tanikawa A., Nomoto K., Nakasato N., 2018, ApJ, 868, 90.
  * Tanikawa et al. (2019) Tanikawa A., Nomoto K., Nakasato N., Maeda K., 2019, ApJ, 885, 103.
  * Vinkó, Szalai, & Könyves-Tóth (2023) Vinkó J., Szalai T., Könyves-Tóth R., 2023, Univ, 9, 244.
  * Wang (2018) Wang B., 2018, RAA, 18, 049.


# Quantum search in sets with prior knowledge

Umut Çalıkyılmaz<EMAIL_ADDRESS>Sadi Turgut<EMAIL_ADDRESS>Middle
East Technical University, Physics Department, Dumlupınar Bulvarı 06800
Çankaya Ankara, Turkey

###### Abstract

Quantum Search Algorithm made a big impact by being able to solve the search
problem for a set with $N$ elements using only $O(\sqrt{N})$ steps.
Unfortunately, it is impossible to reduce the order of the complexity of this
problem, however, it is possible to make improvements by a constant factor. In
this paper we pursued such improvements for search problem in sets with known
probability distributions. We have shown that by using a modified version of
quantum search algorithm, it is possible to decrease the expected number of
iterations for such sets.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

As in any other advancement in science, the idea of quantum computing has
emerged from need. Classical computing is insufficient for simulating complex
quantum systems mainly because the memory needed to store the state of a
quantum system as classical information goes exponentially with the system
size. As a better way of simulating this kind of systems, Richard Feynman
proposed the use of quantum computers, computers employing quantum systems to
store and process data [1, 2]. Shortly after, other advantages of using such
information processing method are noticed. First, the superiority of quantum
computing over classical computing is shown using some specially designed
problems [3, 4, 5, 6]. Then Shor proved the possibility of solving the age-old
factorization problem in polynomial time by using a quantum computer [7].

A few years after that, Grover showed that another classical problem, namely
the search problem, can be solved in a shorter time using a quantum algorithm
[8, 9]. The search problem is a problem of finding an element, satisfying some
condition, inside an unordered set. It is solved classically by trying every
element in the set until a solution is found. Grover’s algorithm works by
applying successive rotations on some initial quantum state until it is
transformed into a desired state. The initial state is the superposition of
all elements in the set with equal coefficients, and the desired final state
is the superposition of only the solutions.

Grover’s algorithm is not able to solve the search problem in a polynomial
time, but it reduces the number of required trials dramatically. In case of a
single solution, which is the only case that we have examined in this paper,
searching a set of $N$ elements classically takes $O(N)$ trials. Grover’s
algorithm finds the solution for the same problem in only $O(\sqrt{N})$
iterations of the Grover operator.

It has been shown that the order of the number of iterations required to solve
this problem using a quantum algorithm cannot be smaller than $O(\sqrt{N})$
[10]. However, some improvement to the Grover’s algorithm by a constant factor
is still possible. The exact Grover algorithm requires the repetition of
Grover operator $\pi/4\sqrt{N}$ times. In Ref. 11, a different method, where
the Grover operator is repeated by a smaller number of times, is introduced.
As finding a solution is not guaranteed at the measurement stage, it might be
necessary to repeat the attempt several times until the solution is found.
With this approach, the average time to find the solution can be reduced by
$12\,\%$.

The aim of this paper is achieving a similar reduction in the case that the
elements of the search space has differing probabilities of being a solution.
It is assumed that for each element $i$, there is a known probability $p_{i}$,
which denotes the probability of element $i$ being the solution. To achieve a
better result, different initial states than that of Grover’s method are
utilized. In these states, coefficients of different indices can be different
from each other. Here, the cost of preparing these states are neglected. The
method using these generic initial states is called as the generic search.

In practical applications of search algorithms, it is sometimes possible to
assume differing probabilities for the elements of the search set. For
example, when brute force is used to crack a password, some combinations are
assumed to be more probable and being tried before others. Quantum search can
be utilized for such a task, and generic search, or any similar method using
generic initial states, can be used to speed up the process.

In the next section of this paper, the procedure and the parameters of generic
search is defined and various equations are derived. The third section is
reserved for finding the optimality conditions for the parameters of the
method. In the fourth section, the estimated values for the optimum values of
the parameters and the expected number of steps are presented. The fifth
section concludes the paper by summarizing our study and giving some ideas
about the implementation of the generic search.

##  II Definition of the Method

###  II.1 The Procedure and the Parameters of the Method

Generic search uses consecutive quantum searches until a solution to the
search problem is found, such as the method given in [11]. The difference of
generic search is that it employs varying initial quantum states for the
searches. Thus, we expect that the optimal number of applied Grover iterations
would also be different for each step.

The $j^{th}$ quantum search in this method is named as the $j^{th}$ step of
the algorithm. The initial state used for the $j^{th}$ step is denoted as
$\left|\psi^{(j)}\right\rangle$ and defined as

$\left|\psi^{(j)}\right\rangle=\sum_{i=0}^{N}c_{i}^{(j)}\left|i\right\rangle~{},$
(1)

where $c_{i}$ are assumed to be real and positive for convenience. At the
$j^{th}$ step, $m^{(j)}$ iterations are applied to this initial state before
measuring the final state. As one might expect, $m^{(j)}<\pi/4\sqrt{N}$ for
each $j$, where $N$ is the number of elements in the search space.

We followed a convention similar to the one used in [12] to define the
parameters of the algorithm. Let $s$ denote the index of the solution that we
sought after. Then, the state representing the solution is
$\left|s\right\rangle$. The non-solution state for step-$j$
$\left|ns^{(j)}\right\rangle$ is defined as

$\left|ns^{(j)}\right\rangle=\frac{1}{\sqrt{1-\left|c^{(j)}_{s}\right|^{2}}}\sum_{i\neq
s}c^{(j)}_{i}|i\rangle~{}.$ (2)

Similar to the original Grover algorithm, the oracle operator $O$ reflects a
state along $|s\rangle$,

$O=\mathds{1}-2|s\rangle\langle s|~{},$ (3)

and the succeeding reflection operator for the $j$th step is

$R^{(j)}=\mathds{1}-2\left|\psi^{(j)}\right\rangle\left\langle\psi^{(j)}\right|~{}.$
(4)

Note that if the state $|\psi^{(j)}\rangle$ can be prepared easily on a
quantum computer, then the operation $R^{(j)}$ can also be implemented
easily111If $G_{j}$ is the unitary map corresponding to the gates that
prepares the state $|\psi^{(j)}\rangle$ from the standard state,
$|\psi^{(j)}\rangle=G_{j}|0|$, then
$R^{(j)}=G_{j}(\mathds{1}-2|0\rangle\langle 0|)G_{j}^{\dagger}$.. We will
assume this for all $|\psi^{(j)}\rangle$ and ignore the complexity of state
preparation.

Let $\left|\xi\right\rangle$ be a quantum state lying in the plane spanned by
$\left|s\right\rangle$ and $\left|ns^{(j)}\right\rangle$. As it was defined in
the Grover algorithm, the Oracle operator reflects a state with respect to the
subspace spanned by all the vectors of non-solution indices. When applied to
$\left|\xi\right\rangle$, $O$ reflects this state with respect to
$\left|ns^{(j)}\right\rangle$ and the reflected state
$\left|\xi^{\prime}\right\rangle$ would still be in the plane. Then, if
$R^{(j)}$ is applied to $\left|\xi^{\prime}\right\rangle$, the resulting state
$\left|\xi^{\prime\prime}\right\rangle$ is also in the plane. The reason is
that $\left|\psi^{(j)}\right\rangle$ is a superposition of
$\left|s\right\rangle$ and $\left|ns^{(j)}\right\rangle$, which can be seen
from Eq. (1), so it lies in the plane spanned by these two states. Moreover,
being a product of two reflections, the operator $G^{(j)}=R^{(j)}O$ acts as a
rotation for $\left|\xi\right\rangle$. The angle of rotation is
$2\phi^{(j)}_{s}$ where $\phi^{(j)}_{s}$ is the angle between
$\left|\psi^{(j)}\right\rangle$ and $\left|ns^{(j)}\right\rangle$,

$\phi^{(j)}_{s}=\arcsin\left|\left\langle\psi^{(j)}|s\right\rangle\right|=\arcsin
c_{s}^{(j)}~{}.$ (5)

In step $j$, the operator $G^{(j)}$ is applied $m^{(j)}$ times to the initial
state $\left|\psi^{(j)}\right\rangle$. Then, the angle between the final state
of step $j$ and $\left|ns^{(j)}\right\rangle$ is

$\theta^{(j)}_{s}=(2m^{(j)}+1)\phi^{(j)}_{s}~{}.$ (6)

One important point is that the index of the solution is not known prior to
the search. Thus, defining the final angle depending on $s$ is not useful.
Instead, the hypothetical final angle in the case that index $i$ is the
solution must be found, for each $i$. This parameter can be defined as

$\theta^{(j)}_{i}=(2m^{(j)}+1)\arcsin c_{i}^{(j)}~{}.$ (7)

Let $P^{(j)}(\mathrm{succ}|i)$ denote the probability that the correct
solution would be identified successfully if $i$ were the solution (if $s=i$).
Then

$\displaystyle P^{(j)}(\mathrm{succ}|i)$ $\displaystyle=\left|\left\langle
i\left|\left(G^{(j)}\right)^{m^{(j)}}\right|\psi_{j}\right\rangle\right|^{2}$
(8) $\displaystyle=\sin^{2}\theta^{(j)}_{i}~{},$

and

$\displaystyle P^{(j)}(\mathrm{fail}|i)$
$\displaystyle=1-P^{(j)}(\mathrm{succ}|i)$ (9)
$\displaystyle=\cos^{2}\theta^{(j)}_{i}~{}.$

In the rest of this paper, the final angles $\theta^{(j)}_{i}$ are used as
primary variables for convenience. Once the $\theta^{(j)}_{i}$ values are
known, the initial angles and the coefficients can be calculated easily using
Eq. (7).

###  II.2 Expected Number of Iterations

The purpose of generic search is to reduce the expected number of Grover
iterations required to solve the search problem. To find the minimum possible
value for this variable, first, its mathematical expression must be
constructed. As a first step, assume that the solution is the index $i$ and
let $E_{i}$ be the expected number of iterations in this case. It is defined
as

$\displaystyle E_{i}$
$\displaystyle=m^{(1)}+\cos^{2}\theta_{i}^{(1)}m^{(2)}+...$ (10)
$\displaystyle=\sum_{j\geq
1}m^{(j)}\prod_{k=1}^{j-1}\cos^{2}\theta_{i}^{(k)}~{},$

where we follow the convention that when $j=1$, the product in the expression
above is equal to $1$, i.e.,

$\prod_{k=1}^{0}\cos^{2}\theta_{i}^{(k)}=1~{}.$ (11)

In Eq. (10), each term in the sum is the number of iterations applied in the
$j^{th}$ step times the probability of ever carrying out the $j^{th}$ step.
For instance, the first quantum search is applied definitely, so the first
term of the sum is $m^{(1)}$. Likewise, the second search is applied in the
case that the solution is not found in the first step, with a probability of
$\cos^{2}\theta_{i}^{(1)}$, so the second term in the sum is
$\cos^{2}\theta_{i}^{(1)}m^{(2)}$.

Since the solution is not known initially, writing the expected number of
iterations in a more general manner is required. This can be done using
$E_{i}$ values defined in Eq. (10) and the probabilities $p_{i}$ as

$E=\sum_{i=0}^{N-1}p_{i}E_{i}=\sum_{i=0}^{N-1}p_{i}\sum_{j\geq
1}m^{(j)}\prod_{k=1}^{j-1}\cos^{2}\theta_{i}^{(k)}~{}$ (12)

Minimizing $E$ requires finding the optimal values of the parameters $m^{(j)}$
and $\theta_{i}^{(j)}$. However, $E$ depends on infinitely many parameters and
optimizing them using numerical methods is impossible. A modified version of
generic search with a finite number of parameters can be used to estimate the
minimum of expected number of iterations.

###  II.3 Modified Method

The procedure is slightly changed to keep the number of parameters finite. In
this modified version, maximum number of steps is set to be $n$. If the
solution is not found in $n-1$ steps, in the $n^{th}$ step,
$m^{(n)}=(\pi/4)\sqrt{N}$ Grover iterations are applied to the initial state
of the original Grover algorithm (which ensures finding the solution).

The expected number of iterations $E$ for this method is very similar to the
original one. The main difference is the number of steps. The equation for $E$
for the modified method is

$E=\sum_{i=0}^{N-1}p_{i}\sum_{j=1}^{n}m^{(j)}\prod_{k=1}^{j-1}\cos^{2}\theta_{i}^{(k)}~{}.$
(13)

The parameters of the last step are predetermined so, it leaves $n-1$ steps to
optimize.

###  II.4 Constraint Equations

The parameters of the method are not independent from each other as it can be
seen in Subsection II.1. A relation between $m^{(j)}$ and $\theta_{i}^{(j)}$
can be found using Eq. (7). For $N\gg 1$, it can be assumed that $c_{i}^{(j)}$
are very small and that $\arcsin c_{i}^{(j)}\approx c_{i}^{(j)}$. By also
assuming that $m^{(j)}$ are $O(\sqrt{N})$, Eq. (7) can be approximated as

$\theta_{i}^{(j)}\approx 2m^{(j)}c_{i}^{(j)}~{}.$ (14)

By taking the square of both sides in this equation, summing with respect to
the index $i$ and using the normalization condition, we obtain the constraint
equation below

$\sum_{i=0}^{N-1}\left[\theta_{i}^{(j)}\right]^{2}=4\left[m^{(j)}\right]^{2}~{}.$
(15)

For generic search, there are infinitely many constraint equations and for the
modified method, there are $n-1$ of them.

##  III Optimality Conditions

The optimum values for the modified method can be found, since there are only
a finite number of parameters. To do this, equations of optimality are derived
using Lagrange multipliers method. To define a function for minimizing the
expected number of iterations, the definition of $E$ Eq. (13), and the
constraint equations Eq. (15) can be used. The resulting function is

$\mathcal{L}=\sum_{i=0}^{N-1}p_{i}\sum_{j=1}^{n}m_{j}\prod_{k=1}^{j-1}\cos^{2}\theta_{i}^{(k)}-\sum_{j=1}^{n-1}\lambda^{(j)}\Gamma^{(j)}~{},$
(16)

where parameter $\lambda^{(j)}$ is the Lagrange multiplier for the $j^{th}$
constraint and

$\Gamma^{(j)}=\sum_{i=0}^{N-1}\left[\theta_{i}^{(j)}\right]^{2}-4\left[m^{(j)}\right]^{2}~{}.$
(17)

At the optimal point, the variation of $\mathcal{L}$ is equal to zero. For its
variation to be equal to zero, the partial derivative of $\mathcal{L}$ with
respect to each parameter must be equal to zero. This can be shown as

$\delta\mathcal{L}=\frac{d\mathcal{L}}{dm^{(j)}}\delta
m^{(j)}+\frac{d\mathcal{L}}{d\theta_{i}^{(j)}}\delta\theta_{i}^{(j)}+\frac{d\mathcal{L}}{d\lambda^{(j)}}\delta\lambda^{(j)}~{}.$
(18)

By equating the partial derivatives to zero, three optimality conditions are
found. These are equations of the optimal values of the parameters $m^{(j)}$,
$\theta_{i}^{(j)}$ and $\lambda^{(j)}$; which are denoted as $m^{(j)*}$,
$\theta_{i}^{(j)*}$ and $\lambda^{(j)*}$ respectively. These optimality
conditions are written as

$\sum_{i=0}^{N-1}p_{i}\prod_{k=1}^{j-1}\cos^{2}\theta_{i}^{(j)*}=8\lambda^{(j)*}m^{(j)*}~{},$
(19)
$\sum_{i=0}^{N-1}\left[\theta_{i}^{(j)*}\right]^{2}=4\left[m^{(j)*}\right]^{2}~{},$
(20)
$p_{i}\cos\theta_{i}^{(\ell)*}\sin\theta_{i}^{(\ell)*}\sum_{j=\ell+1}^{n}m^{(j)*}\prod\limits_{\begin{subarray}{c}k<j\\\
k\neq\ell\end{subarray}}\cos^{2}\theta_{i}^{(k)*}=\lambda^{(\ell)*}\theta_{i}^{(\ell)*}~{}.$
(21)

Solving these equations gives the optimal values for the parameters. Using
these values, the minimum of the expected number of iterations can be found.
Unfortunately, finding analytical solutions to this set of solutions does not
seem possible. Instead, we estimated the optimum values of the parameters
using numerical methods. The procedures of these methods and the results of
the analyses are given in the next section.

##  IV Numerical Analysis

###  IV.1 Procedure

In the numerical analyses, the optimality equations for the modified method
are used to estimate the parameters. In each analysis, number of steps is set
to $n=10$ and the size of the search space is set to $N=10^{6}$. These numbers
are thought to be sufficient to find a close estimation to the original method
for big search spaces. In any case, applying a method with infinitely many
parameters is not applicable in reality. Thus, finding the minimum for the
original method would only determine the lower bound for the expected number
of iterations. The values found in the analyses are considered as estimations
of these limit values.

Optimality conditions given in Eq. (19), Eq. (20) and Eq. (21) are put to use
for the analyses. Initial values for $\lambda^{(j)}$ are set arbitrarily, and
then enhanced iteratively. In this procedure, Newton-Raphson method is
employed to find the inverse of some transcendental function which is derived
from the optimality equations. After updating $\lambda^{(j)}$ values for all
$j$ values, which means at the end of each iteration $m^{(j)}$,
$\theta_{i}^{(j)}$ and $E$ values are calculated. Iteration is halted after
the difference between two successive $E$ values is less than a predetermined
precision value. The details of this numerical analysis can be seen in [14].

Analyses are made for eight probability distributions depicting power law,
exponential and Gaussian behavior. To obtain a discrete probability
distribution from a continuous one, it is assumed that the continuous
distribution function $f(x)$ is a probability distribution on $0\leq x\leq 1$
interval and we defined

$p_{i}=\int_{i/N}^{(i+1)/N}f(x)dx~{},\quad(\textrm{for }i=0,1,\ldots,N-1)~{}.$
(22)

For power law distributions, $f(x)=(n+1)x^{n}$, ($n=0,1,\ldots,5)$. For
exponential distribution we have $f(x)=Ae^{-cx}$ where $c$ is a large value
(we have used $c=30$) and $A$ is for normalization. For half-normal
distribution we have $f(x)=Ae^{-cx^{2}}$ where we have taken $c$ to be $18$.
Since the algorithm does not depend on how the elements are ordered, the fact
that whether $f(x)$ is increasing or decreasing is irrelevant. However, there
is a relation between the standard deviation and $E$ for monotonic
distributions. Thus, we gave the results for monotonic distributions in the
next section. This is also the reason why only half of the normal distribution
is used. If a normal distribution with a mean inside the $(0,1)$ interval were
used, identical results would have been obtained with the half-normal ones. We
also made the same analyses for the randomized versions of the chosen
distributions to verify that ordering of the elements does not change the
expected number of iterations. These randomized versions are random
permutations of the given ordered distributions.

###  IV.2 Results

We estimated the optimal values of the parameters for each selected
probability distribution. Using these values, we calculated the minimum of the
expected number of Grover iterations. These $E$ values are of order
$O(\sqrt{N})$. We also calculated the standard deviations $\sigma$ which are
$O(N)$. To make the results of our study independent of the size of a search
set, we also calculated $E/\sqrt{\sigma}$ values in an effort to find a
constant of the method. These results and the improvements against
$\pi/4\sqrt{N}$, for all the distributions are given in Table 1.

Table 1: The results of the numerical analysis for chosen probability distributions $f(x)$ | $E/\sqrt{N}$ | $E/\sqrt{\sigma}$ | Improvement
---|---|---|---
$1$ | $0.690$ | $1.285$ | $12\,\%$
$2x$ | $0.627$ | $1.292$ | $20\,\%$
$3x^{2}$ | $0.559$ | $1.270$ | $29\,\%$
$4x^{3}$ | $0.507$ | $1.254$ | $35\,\%$
$5x^{4}$ | $0.466$ | $1.241$ | $41\,\%$
$6x^{5}$ | $0.433$ | $1.232$ | $45\,\%$
exp |  $0.213$222This result is valid for exponential distribution with c=30, but not for every exponential distribution. | $1.168$ |  $73\,\%$333The amount of improvement is calculated for the exponential distribution with c=30.
hnorm |  $0.406$444This result is valid for the half-normal distribution with c=18, but not for every half-normal distribution | $1.281$ |  $48\,\%$555The amount of improvement is calculated for the half-normal distribution with c=18.

The first entry in the table belongs to the uniform distribution. For this
distribution, there is no prior knowledge, thus, no improvement on the results
of [11] are obtained as expected. However, for the distributions with lower
standard deviation values, lower $E$ values are found. For these distributions
the decrease in the expected number iterations is more than $12\,\%$. The
highest improvement, which is $73\,\%$, is achieved for the exponential
distribution.

These results suggests that generic search can find the solution in a set with
prior knowledge faster than its predecessors. They also suggest that, for
monotonic probability distributions, there might be a linear relationship
between the square root of the standard deviation and the expected number
iterations to be used in the generic search. Using the results given in 1 and
the least square errors method, this linear relationship is estimated as

$E\approx 1.267\sqrt{\sigma}~{}.$ (23)

How well our results fit in this relation is shown in Figure 1.

$1.267$$1.150$$1.300$$1$$2x$$3x^{2}$$4x^{3}$$5x^{4}$$6x^{5}$$exp$$hnorm$$E/\sqrt{\sigma}$
Figure 1: The fitting of the minimums of the expected number of iterations

The analyses conducted on the randomized versions of the given probability
distributions have verified that the expected number of Grover iterations does
not depend on the order of the elements. For each of the distributions, the
estimated minimum of $E$ is the same for the ordered and randomized versions.
It also shows that Eq. (23) can be used to estimate the minimum of $E$ for any
distribution, not only the monotonic ones. However, for non-monotonic
distributions, the standard deviation of the distribution itself should not be
inserted into the equation. Instead, the standard deviation of the ordered
version of the distribution must be used.

##  V Conclusion

We have shown that generic search method provides an improvement in the
expected number of iterations for search spaces with non-uniform probability
distributions. We also derived an equation that can be used to find an
approximation to the expected number of iterations for any distribution. We
also show that generic search is suitable to be used for non-monotonic
distributions, aside from the monotonic ones.

In this paper, we assumed that any initial state can be prepared without much
complexity. However, some of the initial states to be used to achieve a
minimum expected number of iterations might be hard-to-prepare. Computational
cost of preparing such quantum states might outweigh that of applying Grover
iterations. In that case, generic search does not provide any advantage in the
practical applications. For practical use, some other methods that are similar
to the one introduced here but use easy-to-prepare quantum states can be
designed. Examples of such methods are given in [14]. Obviously, these methods
does not ensure improvements as big as the one ensured by generic search, but
some other methods giving closer results to it can still be developed. Thus,
the results that are given for the generic search shows the theoretical
minimum of the expected number of iterations for any method employing generic
initial states.

## References

  * Feynman [1982] R. P. Feynman, Simulating physics with computers, International Journal of Theoretical Physics 21, 467 (1982).
  * Feynman [1985] R. P. Feynman, Quantum mechanical computers, Optics News 11, 11 (1985).
  * Deutsch [1985] D. Deutsch, Quantum theory, the church-turing principle and the universal quantum computer, Proceedings of the Royal Society A 400, 97 (1985).
  * Deutsch and Jozsa [1992] D. Deutsch and R. Jozsa, Rapid solution of problems by quantum computation, Proceedings of the Royal Society A 439, 553 (1992).
  * Simon [1994] D. R. Simon, On the power of quantum computation, in  _Proceedings of the 35th Annual Symposium on Foundations of Computer Science_ (IEEE, 1994).
  * Bernstein and Vazirani [1994] E. Bernstein and U. Vazirani, Quantum complexity theory, in  _Proceedings of the 25th Annual ACM Symposium on Theory of Computing_ (ACM, 1994).
  * Shor [1994] P. W. Shor, Polynomial-time algorithms for quantum computation: Discrete logarithms and factoring, in  _Proceedings of the 35th Annual Symposium on Foundations of Computer Science_ (IEEE, 1994).
  * Grover [1996] L. K. Grover, A fast quantum mechanical algorithm for database search, in  _Proceedings of the 28th Annual ACM Symposium on Theory of Computing_ (ACM, 1996).
  * Grover [1997] L. K. Grover, Quantum mechanics helps in searching for a needle in a haystack, Physical Review Letters 79, 325 (1997).
  * Bennett  _et al._ [1997] C. H. Bennett, E. Bernstein, G. Brassard, and U. Vazirani, Strenghts and weaknesses of quantum computing, SIAM Journal on Computing 26, 1510 (1997).
  * Boyer  _et al._ [1998] M. Boyer, G. Brassard, P. Høyer, and A. Tapp, Tight bounds on quantum searching, Fortschritte der Physik 46, 493 (1998).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information (10th aniversary edition) (Cambridge University Press, 2010) Chap. Quantum Search Algorithm, pp. 248–276.
  * Note [1] If $G_{j}$ is the unitary map corresponding to the gates that prepares the state $|\psi^{(j)}\rangle$ from the standard state, $|\psi^{(j)}\rangle=G_{j}|0|$, then $R^{(j)}=G_{j}(\mathds{1}-2|0\rangle\langle 0|)G_{j}^{\dagger}$.
  * Çalıkyılmaz [2021] U. Çalıkyılmaz,  _Quantum Search in Sets with Prior Knowledge_ , Master’s thesis, Middle East Technical University (2021).