# Material Coherence from Trajectories via Burau Eigenanalysis of Braids

Melissa Yeung, David Cohen-Steiner, Mathieu Desbrun

<figure><img src="image_1.png" /></figure>

## To cite this version:

Melissa Yeung, David Cohen-Steiner, Mathieu Desbrun. Material Coherence from Trajectories via Burau Eigenanalysis of Braids. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2020, 30 (3), 10.1063/1.5128269. hal-02295987

HAL Id: hal-02295987

https://inria.hal.science/hal-02295987

Submitted on 24 Sep 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Material Coherence from Trajectories via Burau Eigenanalysis of Braids

Melissa Yeung,¹ David Cohen-Steiner,² and Mathieu Desbrun³

¹) Caltech

²) Inria - Méditerranée

³) Caltech/ShanghaiTech

(Dated: 24 September 2019)

In this paper, we provide a numerical tool to study material coherence from a set of 2D Lagrangian trajectories sampling a dynamical system, i.e., from the motion of passive tracers. We show that eigenvectors of the Burau representation of a topological braid derived from the trajectories have levelsets corresponding to components of the Nielsen-Thurston decomposition of the dynamical system. One can thus detect and identify clusters of space-time trajectories corresponding to coherent regions of the dynamical system by solving an eigenvalue problem. Unlike previous methods, the scalable computational complexity of our braid-based approach allows the analysis of large amounts of trajectories.

Studying two-dimensional flows and their induced transport and mixing properties is key to geophysical studies of atmospheric and oceanic processes. However, one often has only sparse tracer trajectories (e.g., positions of buoys in time) to infer the overall flow geometry. Fortunately, topological methods based on the theory of braid groups have recently been proposed to extract structures from such a sparse set of trajectories by measuring their entanglement. This braid viewpoint offers sound foundations for the definition of coherent structures. Yet, there has been only limited efforts in developing practical tools that can leverage topological properties for the efficient analysis of flow structures: handling a larger number of trajectories remains computationally challenging. We contribute a new and simple computational tool to extract Lagrangian structures from sparse trajectories by noting that the eigenstructure of the Burau matrix representation of a braid of particle trajectories can be used to reveal coherent regions of the flows. Detection of clusters of space-time trajectories corresponding to coherent regions of the dynamical system can thus be achieved by solving a simple eigenvalue problem. This paper establishes the theoretical foundations behind this braid eigenanalysis approach, along with numerical validations on various flows.

## I. Introduction

Studying two-dimensional flows and their induced transport is key to geophysical studies of atmospheric and oceanic processes. In particular, identifying regions of similar dynamical fate within experimental or numerical temporal data often leads to a better understanding of the overall flow geometry, from which material transport and mixing properties can be accurately quantified. Such an analysis can also be crucial in short-term forecasting of pattern evolution in complex 2D dynamical systems. While classic dynamical systems theory has established solid foundations (such as fixed points, orbits, and stable manifolds) to study structures in steady and time-periodic fluid flows, these tools become inadequate to capture the complexity of aperiodic large-scale flows: tracers in even a very simple Eulerian velocity field can undergo very different and intricate trajectories—a phenomenon often referred to as Lagrangian chaos, to distinguish it from its Eulerian counterpart (turbulence)¹. Consequently, numerous definitions of

the notion of *coherent regions* in flows have been formulated to help in this endeavor, from which algorithms can be developed to automatically generate a simplified skeleton of the overall dynamics of the system from the input data.

### A. Deformation-based coherent structures.

A large majority of approaches to material coherence discovery involve quantifying how the velocity field induces, in time, a deformation of the original space (see² and³ for two recent reviews). Instantaneous Eulerian diagnostics tools⁴,⁵ were first proposed to find coherent features in instantaneous velocity fields based on the eigenvalues of their gradient tensors. Later, the notion of *Lagrangian coherent structures* (a term originally coined by Haller and Yuan⁶) was introduced to describe material curves that delineate regions with qualitatively different tracer dynamics. These boundaries between coherent regions are, to a certain extent, finite-time analogs of invariant manifolds and can be derived from the finite-time Lyapunov exponent field of the flow⁷,⁸. Ensued a series of computational approaches with various theoretical guarantees and varying computational efficacy, each offering a segmentation of the flow into regions of coherent material transport⁹⁻¹². More recently, these separatrices were reformulated through local geometric properties of the flow map and its Cauchy-Green strain tensor field as material curves with locally maximal repelling, attracting, or shearing impact on neighboring fluid elements¹³,¹⁴. More global geometric definitions, such as least stretching material lines or stationary curves for the tangential stretching, were also shown to reveal further structures that include weakly but coherently stretching material vortices and Lagrangian jet cores¹⁵. Finally, authors have also classified directly coherent regions (as opposed to the separators between regions) via ergodicity-based¹⁶, observer-based¹⁷, and probabilistic¹⁸ approaches.

### B. Topology-based coherent structures.

While there is no universally agreed upon definition for coherent structures, most of the definitions and computational methods mentioned above assume a fine knowledge of the Eulerian velocity field in space and time to deduce a good approximation of the flow map. However, oceanic flows are often known only as a set of sparse particle trajectories in time (typically, buoys). Such a sparse sampling of the dynamical system does not lend itself well to a geometric analysis
of transport, rendering previous approaches unable to provide qualitative insights of the underlying flow. Fortunately, topological methods have recently been proposed to extract structures from a sparse set of trajectories by measuring their entanglement¹⁹. Topology and group-theoretic tools have been commonplace in the study of fluid dynamics since the early work of Leray²⁰ and Arnold²¹; in the case of particle trajectories, the theory of *braid groups*, a classical area of topology, was shown to be relevant to the analysis of mixing in two-dimensional flows²²—not only for periodic systems, but for aperiodic, finite-time systems as well. A set of trajectories (whether these are trajectories of particles in the flow or trajectories of stirring rods that induce a motion in the surrounding fluid) has a topological structure that can be represented using Artin's braid group²³. The resulting braid viewpoint offers a precise and versatile tool to distinguish different regimes of dynamics, to provide a lower bound on topological entropy, and in the case of positive topological entropy, to reveal the presence of topological chaos. Allshouse and Thiffeault²⁴ introduced such a topological method for detecting material coherence from a small set of particle trajectories. By examining the growth of loops that surround sets of particles trajectories, they identify coherent structures as loops that have negligible growth. Coherent regions are then defined as containing particles that possibly mix with other particles within the region itself but do not mix with particles outside the region; the set of trajectories arising from the particles within a coherent region forms a *coherent bundle*.

### C. Computational challenges.

Even if the use of braid groups offers sound foundations for the definition of coherent structures, there has been only limited efforts in developing practical computational tools for the efficient analysis of flow structures. A significant exception is the work of Thiffeault and Budišić²⁵,²⁶ which allows for the construction of an algebraic braid from trajectories by tracking the crossing of particles in space-time with respect to the projection onto an arbitrary line. Resulting braids can then be analyzed, allowing one to estimate the topological entropy of the flow, as well as compute the growth of loops under the action of the braid. However, finding coherent sets of trajectories that travel together but do not entangle other trajectories require an exhaustive search through loops enclosing pairs of particles to find the least growing ones. Drmatic improvement in efficiency was offered in²⁴ by limiting the exploration to “pair-loops” of nearby punctures, turning an exponential search into a quadratic one without significant artifacts—although it may end up merging multiple coherent regions into one single larger region. Yet, handling a larger number of trajectories remains computationally challenging. Allshouse and Thiffeault²⁴ examine up to 20 trajectories, while the more recent work of Francois et al.²⁷ examine up to 100 trajectories.

### D. Contributions.

In this paper, we present a scalable and efficient approach for detecting coherent sets from sparse tracer trajectories of a two-dimensional dynamical system. We adopt the trajectories-as-a-braid viewpoint pioneered by Thiffeault et

al., to which we contribute a new and simple computational tool to extract Lagrangian structures by noting that the eigenstructure of the *Burau matrix representation* of a braid of particle trajectories can be used to reveal coherent regions of the flows.²⁸ More precisely, we show that the *Burau matrix* of a reducible braid possesses eigenvectors that are almost piecewise constant over coherent sets. Thus, we deduce that one can detect and identify clusters of space-time trajectories corresponding to coherent regions of the dynamical system by solving an eigenvalue problem.

This paper both establishes the theoretical foundations behind this braid eigenanalysis approach and presents numerical validations on different flows. In particular, we apply our approach to study Aref's blinking vortex flow²⁹ and a modified Duffing oscillator²⁴. In the analysis of Aref's blinking vortex flow, our method distinguishes chaotic regions from KAM surfaces, over a range of flow strengths. In the analysis of the modified Duffing oscillator, we also detect two additional limit cycles in addition to detecting the two dominant regions of mixing detected by Allshouse and Thiffeault²⁴ due to our ability to handle a significantly larger amount of trajectories. We also discuss the computational complexity of our approach and contrast it with existing approaches which seek individual slow-growing loops to identify coherent regions.

## II. Braids and Dynamics

Throughout this paper, we assume that a flow on a two-dimensional domain $\mathcal{D} \subset \mathbb{R}^2$ is known via a set of $n$ disjoint particle trajectories $\{\mathbf{x}_i(t)\}_{i=1..n}$, where each trajectory $\mathbf{x}_i: [0, 1] \to \mathcal{D}$ records the path of the $i$th particle in time. Each of the trajectories can be viewed as a *strand* in space-time. Since these strands can be arbitrarily intertwined, the set of all trajectories resembles a “physical” space-time braid of entangled strands (see Fig. 1(a)). A topological analysis of this braid-like structure provides, remarkably, important dynamical information of the sampled flow.

### A. Braid terminology

We now define a few notions related to braids that will be relevant in our exposition for characterizing flows.

#### 1. Geometric and pure braid.

From a mathematical standpoint, we say that the collection $\{\mathbf{x}_i(t)\}_{i=1..n}$ of strands forms a *geometric braid* when the set of final positions $\mathbf{x}_i(1)$ are given by a permutation of the initial positions $\mathbf{x}_i(0)$. If this permutation is the identity, the braid is said to be a *pure braid*.

#### 2. Topological braid.

One can further discard geometric information (e.g., distances between particles) to form a more concise, purely topological description of a geometric braid: a *topological braid* is an *equivalence class* of geometric braids, where two geometric braids are considered equivalent when one can be deformed continuously into the other without having strands crossing in space-time. The product of two topological braids on $n$ strands can then be formed by concatenation, defining a group $B_n$ called the braid group.
<figure><img src="image_3.png" /><figcaption>FIG. 1. **From trajectories to topological braid.** (a) A set of space-time trajectories forms a physical braid as the trajectories entangle in time; (b) by tracking the order of the particles along a projection line as time goes by, one can construct a topological braid encoding the trajectories.</figcaption></figure>

## B. Encoding trajectories as a topological braid

A topological description of a set of particle trajectories can be determined by projecting the strands to a fixed space-time plane $\mathbb{R} \times [0, 1]$ and only retaining information about how strands pass over one another<sup>25</sup>. More specifically, we select an arbitrary projection line, and enumerate the strands $i = 1, \dots, n$ according to their ordering along the projection line (Fig. 1(b)). As particles move in time, and their strands cross one another, the projection of the strands onto the projection line will change, and the enumeration is updated accordingly. For $1 \le i < n$, we let $\sigma_i$ denote the braid consisting of a single crossing given by passing the $i$th strand behind the $(i+1)$st strand (Fig. 2(a)); one can thus think of the braid $\sigma_i$ as a simple clockwise half-twist interchanging strands $i$ and $i+1$. Conversely, $\sigma_i^{-1}$ denotes the single crossing given by passing the $i$th strand in front of the $(i+1)$st strand. Then the sequence of crossings in time can be recorded by a concatenation (product) of $\sigma_i$. (Following the standard practice in braid literature, we adopt the convention of composing our braids from left to right.) This sequence specifies the topological braid $\sigma$ corresponding to the intertwined input trajectories.

This composition of braids $\sigma_i$ is a convenient and well-studied description of a topological braid, namely, the $\sigma_i$ are the generators used in Artin's presentation<sup>23</sup> of the braid group $B_n$. An additional observation we will use is that if a set of particle trajectories forms a geometric braid, then a change in the projection line only changes the resulting topological braid by conjugation<sup>19,25,30</sup>.

## C. Braid dynamics

We now point out two important relationships between topological braids and the dynamical system they encode.

### 1. Action on closed curves.

For each braid on $n$ strands, we can associate an isotopy class of homeomorphisms on an $n$-punctured disk that fix the boundary pointwise<sup>31</sup>, where each generator $\sigma_i$ is associated to a homeomorphism from $D_n$ to $D_n$ that has support on a twice-punctured disk and that is described on this support by a half-twist. Simple closed curves in $D_n$ encompass-

ing two or more punctures may thus be subject to stretching and folding under the associated homeomorphism. Consequently, the braid group induces a (right) action on the set of isotopy classes of simple closed curves. This action can be exploited to study the flow: for instance, the growth rate of a loop around punctures under the action of a braid provides a lower bound on the topological entropy of the surrounding flow<sup>19</sup>.

### 2. Reducible braids.

Nielsen-Thurston theory<sup>31</sup> distinguishes a special class of braids, called *reducible braids*, which are particularly relevant to the study of dynamical systems. A reducible braid $\alpha$ on $n$ strands preserves a family $\mathcal{C}$ of disjoint simple closed curves, each enclosing more than one but fewer than $n$ punctures. The family $\mathcal{C}$ is referred to as a *reduction system* for $\alpha$.

In general, the reduction system of a reducible braid may be very complex geometrically, but every reducible braid is *conjugate* to a braid whose reduction system consists of a family of round curves (geometric ellipses)<sup>32–34</sup>, which we call a *round reduction system*. A reducible braid with a round reduction system can be described by a collection of tubes (i.e., tubular braids) and strands: each tube is traced out by the path of a simple closed curve in $\mathcal{C}$, and each tube may enclose other tubes or strands—see Fig. 2 for an illustration. We call particle trajectories enclosed by a tube a *coherent bundle*, and we call their initial positions (at time $t=0$) a *coherent region* (or coherent set), following Thiffeault’s nomenclature<sup>19</sup>.

<figure><img src="image_4.png" /><figcaption>FIG. 2. **Braids and dynamics.** (a) The topological braid of a set of trajectories can be encoded through a composition of Artin generators $\sigma_i$ (where trajectory $\mathbf{x}_i$ goes behind trajectory $\mathbf{x}_{i+1}$ as depicted) and their inverses; (b) a schematic example of 7 particles evolving in time (starting at the top) where the resulting braid is reducible as there are three tubular braids, formed by pairs or a triplet of particles.</figcaption></figure>

## D. Burau representation

A convenient tool for studying a braid, called the Burau representation<sup>35</sup>, is given by associating to each Artin generator
$\sigma_i$ the $n \times n$ block matrix

$$ \forall 1 \le i \le n, \quad \sigma_i \mapsto \begin{pmatrix} I_{i-1} & 0 & 0 & 0 \\ 0 & 1-s & s & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & I_{n-i-1} \end{pmatrix}, \quad (1) $$

parameterized by a (real or complex) variable $s$, where $I_k$ denotes the $k \times k$ identity matrix. For conciseness, we will sometimes express the Burau representation of $\sigma_i$ as

$$ \mathbf{B}[\sigma_i](s) := I_{i-1} \oplus \begin{pmatrix} 1-s & s \\ 1 & 0 \end{pmatrix} \oplus I_{n-i-1} $$

where we express the block structure of the matrix representation as a direct sum $\oplus$. A concatenation of generators corresponds to a sequence of multiplications of the aforementioned block matrices, resulting in an $n \times n$ Burau matrix whose elements are polynomials in $s$ with coefficients in $\mathbb{Z}$. This Burau representation is thus a homomorphism from the braid group $B_n$ to the general linear group of degree $n$. While the Burau representation is known to be not faithful for $n \ge 5^{36,37}$ (with a rather large kernel for $n \ge 6^{38}$), this matrix representation does nonetheless allow us to examine the dynamics of a motion of particles through a study of its eigenvectors.

The reduced Burau representation is a variant obtained by taking the quotient of the unreduced representation by the invariant subspace of the matrices $\mathbf{B}[\sigma_i](s)$ generated by the constant vector $(1, \dots, 1)$.

### III. Coherence from Burau eigenanalysis

We now state and prove the main theoretical result of this paper, which will be at the core of our computational approach to finding coherent structures from a set of trajectories.

**Main Result.** Let $\beta \in B_n$ be a reducible pure braid on $n$ strands with reduction system $\mathcal{C}_\beta$ consisting of $k$ simple closed curves. For each $\varepsilon > 0$, there exists $\delta > 0$ such that for $s \in \mathbb{C}$ with modulus 1 and argument at most $\delta$, the reduced Burau matrix $\bar{\mathbf{B}}[\beta](s)$ has $k$ eigenvectors whose span makes an angle at most $\varepsilon$ with the space of piecewise constant vectors on $(D_n, \mathcal{C}_\beta)$ modulo constant vectors.

In the remainder of this section, we prove this result in three main stages.

#### A. Decomposing a reducible braid

We begin with the claim that a reducible braid $\alpha \in B_n$ that preserves a family $\mathcal{C}$ of round curves can be written as a product of tubular braids with trivial interior braiding, followed by a product of interior braids, with trivial braiding between tubular braids.

<figure><img src="image_5.png" /></figure>

We follow the notation of Band and Boyland$^{39}$ and refer to the braid that moves the group of $\eta_1$ consecutive strands starting at strand $i$ behind the group of $\eta_2$ consecutive strands starting at strand $i + \eta_1$ (see inset) as:

$$ \sigma_{i, \eta_1, \eta_2} = (\sigma_{i+\eta_1-1} \cdots \sigma_{i+\eta_1+\eta_2-2}) (\sigma_{i+\eta_1-2} \cdots \sigma_{i+\eta_1+\eta_2-3}) \cdots (\sigma_i \cdots \sigma_{i+\eta_2-1}). \quad (2) $$

**Lemma 1.** Let $\alpha \in B_n$ be a reducible braid that preserves a finite non-empty family $\mathcal{C}$ of non-nested round curves. For each puncture $p_r$ of $D_n$ not enclosed by any curve of $\mathcal{C}$, extend $\mathcal{C}$ by adding a round curve $C_r$ enclosing the single puncture $p_r$. Let $k$ be the number of curves in the extended family $\mathcal{C}$.

Then there exists a finite sequence of tuples $(i(\ell), \eta_1(\ell), \eta_2(\ell))$ of non-negative integers and corresponding exponents $\varepsilon_\ell \in \{-1, +1\}$ for the generators of the tubular braids, as well as a finite sequence of (interior) braids $\beta_j$, $1 \le j \le k$, each supported on a punctured disk with boundary given by a reduction curve $C_j \in \mathcal{C}$, such that the braid $\alpha$ can be expressed as:

$$ \alpha = \prod_{\ell=1}^{L} \sigma_{i(\ell), \eta_1(\ell), \eta_2(\ell)}^{\varepsilon_{\ell}} \cdot \prod_{j=1}^{k} \beta_j. $$

*Proof.* Let $E_j$ denote the punctured disk enclosed by the curve $C_j \in \mathcal{C}$. Let

$$ \hat{D}_k = D_n \setminus \bigcup_{j=1}^{k} E_j $$

denote the $n$-punctured disk with $k$ disks (each containing punctures) removed. Then the braid $\sigma$ induces an automorphism $\hat{\sigma}$ on $\hat{D}_k$. The induced automorphism $\hat{\sigma}$ describes the tubular braid structure of $\sigma$ with respect to $\mathcal{C}$.

If we collapse each hole of $\hat{D}_k$ to a puncture, we can consider $\hat{D}_k$ as a $k$-punctured disk. Therefore, $\hat{\sigma}$ can be given by a sequence of Artin generators for the braid group $B_k$:

$$ \hat{\sigma} = \prod_{a_l} \hat{\sigma}_{a_l}^{\varepsilon_l} $$

where $1 \le a_l < k$ and $\varepsilon_l = \pm 1$ for all $\ell$. Viewed from above, each generator $\hat{\sigma}_{a_l}$ corresponds to a clockwise half-twist interchanging the holes obtained by removing $E_{a_l}$ and $E_{a_l+1}$. Thus, each $\hat{\sigma}_{a_l}^{\varepsilon_l}$ specifies a tuple $(\varepsilon_\ell, i(\ell), \eta_1(\ell), \eta_2(\ell))$, given by the direction of the half-twist, the minimum index of the points contained in $E_{a_l}$, the number of punctures in $E_{a_l}$, and the number of punctures in $E_{a_l+1}$, respectively. Defining the braids

$$ \beta_j = \prod_{\ell=L}^{1} \sigma_{i(\ell), \eta_1(\ell), \eta_2(\ell)}^{-\varepsilon_\ell} \cdot \sigma|_{E_j} $$

completes the decomposition. $\square$

#### B. Piecewise-constant vectors

Next, we show that the Burau representation $\mathbf{B}[\alpha](s)$ of a reducible braid $\alpha$ with a round reduction system maps a piecewise constant vector $v$ that respects the structure of $\alpha$ to a piecewise-constant vector. More precisely, we show that if $v$ is constant on the interior of every curve of a round reduction system $\mathcal{C}$, then the image $\mathbf{B}[\alpha](s)v$ is a piecewise constant vector, which is constant on the interior of every curve of $\alpha(\mathcal{C}) = \mathcal{C}$.
**Definition 2** (piecewise constant vectors). Let $\mathcal{C}$ be a collection of non-nested pairwise-disjoint simple closed curves in the n-punctured disk $D_n$. We say that a vector $\mathbf{v} = (v_1, \dots, v_n)$ is piecewise constant on $(D_n, \mathcal{C})$ if $\mathbf{v}_r = \mathbf{v}_{r'}$ whenever the punctures $p_r$ and $p_{r'}$ are enclosed by the same $C_{jr} \in \mathcal{C}$.

The tubular braid $\sigma_{i,\eta_1,\eta_2}$ described in Eq. (2) has an obvious round reduction system $\mathcal{C}$. The following observation is the basis of our method:

**Lemma 3.** If $\mathbf{v}$ is a vector that is piecewise constant on $(D_n, \mathcal{C})$, the image $\mathbf{B}[\sigma_{i,\eta_1,\eta_2}](s)\mathbf{v}$ is also piecewise constant on $(D_n, \sigma_{i,\eta_1,\eta_2}(\mathcal{C}))$ for all $s$.

*Proof.* Because the braid $\sigma_{i,\eta_1,\eta_2}$ has a grid structure with all its crossing of the same sign, the product of the Burau matrices of the involved generators is easily computed and yields that $\sigma_{i,\eta_1,\eta_2}$ has Burau matrix

$$I_{i-1} \oplus \left( \frac{\begin{array}{ccccccc} 1-s & s-s^2 & \cdots & s^{\eta_2-1}-s^{\eta_2} & s^{\eta_2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & & \vdots \\ 1-s & s-s^2 & \cdots & s^{\eta_2-1}-s^{\eta_2} & 0 & \cdots & s^{\eta_2} \\ \hline & & & & I_{\eta_2} & & \\ & & & & 0 & \cdots & s^{\eta_2} \end{array}}{I_{\eta_2}} \right) \oplus I_{n-i-\eta_1-\eta_2+1},$$

From the Burau matrix structure, it is clear that piecewise-constant vectors are sent to piecewise-constant vectors. $\square$

We can easily generalize Lemma 3 and show that the image of a piecewise constant vector under any reducible braid that preserves a finite non-empty family of round curves is also piecewise constant:

**Lemma 4** (invariance of piecewise-constant vectors). Let $\alpha$ be a reducible braid with a reduction system $\mathcal{C}$ consisting of $k$ round curves. Let $\mathbf{v}$ be a vector in $\mathbb{C}^n$ that is piecewise constant on $(D_n, \mathcal{C})$. Then the image $\mathbf{B}(\alpha)(s)\mathbf{v}$ is a vector that is piecewise constant on $(D_n, \mathcal{C}) = (D_n, \alpha(\mathcal{C}))$ for all $s$.

*Proof.* Write $\sigma$ in the form

$$\sigma = \prod_{\ell=1}^{m} \sigma_{i(\ell), \eta_1(\ell), \eta_2(\ell)}^{\epsilon_{\ell}} \cdot \prod_{j=1}^{k} \beta_j,$$

where $(i(\ell), \eta_1(\ell), \eta_2(\ell))$, $\epsilon_\ell$, and $\beta_j$ are given in the proof of Lemma 1. Since Burau matrices fix constant vectors, then for each $1 \le j \le k$ we have that $\mathbf{B}[\beta_j](s)\mathbf{v} = \mathbf{v}$ for all vectors $\mathbf{v}$ that are constant on the punctured disk $E_j$ enclosed by the curve $C_j$. Furthermore, by Lemma 1, the Burau matrices of the braids $\sigma_{i(\ell), \eta_1(\ell), \eta_2(\ell)}^{\epsilon_\ell}$ also send piecewise constant vectors to piecewise constant vectors. The claim follows. $\square$

In the case of pure braids, both spaces of piecewise constant vectors appearing in the above lemma coincide. We can thus expect the existence of piecewise constant eigenvectors. To do so, we use the fact$^{40}$ that reduced Burau matrices preserves a specific sesquilinear form, which is known to be positive definite when the Burau parameter $s$ has unit modulus and is close enough to 1 (more precisely, when the argument is less than $1/n$). The same is true for the reduced Burau matrix restricted to the space of piecewise constant vectors modulo constant, which is invariant by the above lemma. Hence this restriction is conjugate to a unitary matrix, implying that it is always diagonalizable. We thus deduce:

**Proposition 5.** Let $\alpha$ be a reducible pure braid with round reduction system $\mathcal{C}$. If $s$ has modulus 1 and argument less than $1/n$, the space of piecewise constant vectors on $(D_n, \mathcal{C})$ modulo constant vectors has a basis consisting of eigenvectors of the reduced Burau matrix $\bar{\mathbf{B}}[\alpha](s)$.

### C. Putting it all together

Now, recall that every reducible braid $\beta \in B_n$ is conjugate to a reducible braid $\alpha$ that preserves a family of round curves, and we can write $\beta = \gamma\alpha\gamma^{-1}$, for some $\gamma \in B_n$. Their corresponding reduced Burau matrices can thus be written as

$$\bar{\mathbf{B}}[\beta](s) = \bar{\mathbf{B}}[\gamma](s)\bar{\mathbf{B}}[\alpha](s)\bar{\mathbf{B}}[\gamma^{-1}](s).$$

When $s = 1$, $\bar{\mathbf{B}}[\gamma](s)$ is a permutation matrix. Since the Burau representation $\bar{\mathbf{B}}[\gamma](s)$ is continuous in $s$, then when $s$ is close to 1, the Burau matrix $\bar{\mathbf{B}}[\gamma](s \approx 1)$ is close to a permutation matrix. Thus, if $\bar{\mathbf{B}}[\alpha](s \approx 1)$ has an eigenvector that is piecewise constant on its components, then $\bar{\mathbf{B}}[\beta](s \approx 1)$ has an eigenvector that is almost piecewise constant on its components, yielding in the main result. More formally:

*Proof of the Main Result.* Let $\beta$ be a reducible pure braid on $n$ strands with reduction system $\mathcal{C}_\beta$. There exists a braid $\alpha$ conjugate to $\beta$ in $B_n$, with $\beta = \gamma\alpha\gamma^{-1}$ for some $\gamma \in B_n$, such that $\alpha$ has a round reduction system $\mathcal{C} = \gamma^{-1}(\mathcal{C}_\beta)$, whose curves are given by $\gamma^{-1}(C_j)$, for $C_j \in \mathcal{C}_\beta^{41}$. We consider $\gamma$ as an automorphism of the punctured disk as necessary.

Let their corresponding reduced Burau matrices be denoted

$$A(s) = \bar{\mathbf{B}}[\alpha](s)$$

$$B(s) = \bar{\mathbf{B}}[\beta](s)$$

$$C(s) = \bar{\mathbf{B}}[\gamma](s),$$

Then $B(s) = C(s)A(s)C(s)^{-1}$.

By Proposition 5, since $\alpha$ is a reducible pure braid with a round reduction system $\mathcal{C}$, there exists $\delta_A > 0$ such that for all $s$ with modulus 1 and satisfying $|s-1| < \delta_A$, the reduced Burau matrix $A(s)$ has a set of eigenvectors forming a basis of the space $E$ of piecewise-constant vectors on $(D_n, \mathcal{C})$.

If $\mathbf{v}(s)$ is such an eigenvector, $C(s)\mathbf{v}(s)$ is an eigenvector of $B(s)$ since

$$\begin{aligned} B(s)C(s)\mathbf{v}(s) &= C(s)A(s)C(s)^{-1}C(s)\mathbf{v}(s) \\ &= C(s)A(s)\mathbf{v}(s) \\ &\propto C(s)\mathbf{v}(s) \end{aligned}$$

Also because $C(1)$ is induced by the permutation matrix given by the action of $\gamma$ on the punctures, $C(1)E$ is the space of piecewise constant vectors on $(D_n, \gamma(\mathcal{C})) = (D_n, \mathcal{C}_\beta)$.

Now since $C(s)$ is continuous in $s$ at $s = 1$, there exists $\delta_C > 0$ such that for all $s$ satisfying $|s-1| < \delta_C$, we have that

$$\|C(s)C(1)^{-1} - I\| < \varepsilon. \quad (3)$$

Since

$$C(s)E = C(s)C(1)^{-1}C(1)E,$$

the claim follows. $\square$
# IV. Numerical implementation

Based on the previously discussed eigen-properties of the Burau representation of a pure reducible braid, we now describe an algorithmic approach to detecting coherent regions from a set of trajectories.

## A. From trajectories to braid

Given a 2D domain $\mathcal{D}$ and a set of $n$ particle trajectories $\mathbf{x}_j(t_k)$ in $\mathcal{D}$ for $1 \le j \le n$ and for time step indices $1 \le k < l$ we first convert the input into a topological braid by artificially closing each trajectories so that a pure braid can be constructed as explained in Sec. II.

## 1. Closure of trajectories.

So far, we have discussed the braid-theoretic interpretation of a set of disjoint particle trajectories that form a pure braid, whose properties reflect the underlying dynamics of a flow or a mechanical system. In a pure braid, the set of initial positions is the same as the set of final positions. However, in most physical systems, particle trajectories are not periodic, and the set of final positions are not typically the same as the set of initial positions. Consequently, the resulting algebraic braids formed by changing the projection line in Sec. II will not necessarily be conjugate$^{19,30}$. In order to fix this issue, we form a pure braid by appending the set of initial positions to the end of the set of particle trajectories; i.e., for each trajectory $(\mathbf{x}_i(t_0), \cdots, \mathbf{x}_i(t_{l-1}))$, we append the initial position to form the trajectory

$$ (\mathbf{x}_i(t_0), \cdots, \mathbf{x}_i(t_{l-1}), \mathbf{x}_i(t_l) \equiv \mathbf{x}_i(t_0)). $$

Note that this artificial closure may change the reducibility of the braid; but we will show in Sec. V that this theoretical shortcoming does not affect our ability to detect coherent sets in practice: when trajectories are sufficiently long, this small modification of the trajectories does not alter the overall dynamics significantly. Trajectory closures do, however, guarantee that the modified trajectories result in a pure braid, by construction.

## 2. Braid construction.

We use BRAIDLAB$^{25}$ to compute the algebraic braid $\{\sigma_i\}_{i=1..n}$ corresponding to the motion of the particles, with respect to some arbitrary projection line. We note that if our initial positions happen to be on a regular grid, one should perturb the initial positions and/or choose a projection line that is not aligned with the grid lines: this will help resolve some coincident particles on the projection line, allowing the topological braid to be computed. The resulting pure braid is then encoded as a sequence $\sigma$ of length $L$ of Artin generators $\sigma_i$ corresponding to the $L$ crossings of the piecewise-linear trajectories with respect to the projection line.

## B. From braid to Burau matrix

The sequence $\sigma$ of Artin generators specifies a Burau matrix corresponding to the braid of particle trajectories: for each generator $\sigma_i$ of the topological braid, its corresponding Burau matrix $\mathbf{B}[\sigma_i](s)$ is given by Eq. (1), and the product $\prod_{l} \mathbf{B}[\sigma_{at_l}](s)$ of these matrices in the order of the sequence $\sigma$ is the Burau matrix $\mathbf{B}$ of the braid.

## 1. Choice of Burau parameter.

Based on our results from Sec. III, we must pick a value $s$ such that $|s| = 1$ and $s \approx 1$ to construct the Burau matrix for the flow. We should also choose $s$ not to be too close to 1 so that $\mathbf{B}[\sigma](s)$ is not too numerically close to the identity matrix $\mathbf{B}[\sigma](1)$. Besides the guarantee of diagonalizability afforded by Proposition 5, one of the advantages of picking a value $s$ with a magnitude of one is that it prevents the rapid growth of the Burau matrix coefficients, which could grow exponentially over long durations otherwise. In practice, choices of $s$ with $|s-1|$ between $10^{-4}$ and $10^{-12}$ worked consistently well. In all our examples, we used $s = \exp(10^{-8}i)$.

## C. Finding coherent regions

We then proceed to an eigendecomposition of the Burau matrix $\mathbf{B}$. While the guarantees of Sec. III were proved for reduced Burau matrices, we found that using unreduced Burau matrices worked as well in practice. It is in fact not difficult to see that the theoretical guarantees can be extended to unreduced Burau matrices except in certain cases that are unlikely to happen; we omit further discussions on this point here, and focus instead of numerical details.

## 1. Eigensolver.

There are a number of numerical libraries (MATLAB, LAPACK, etc) that can be used to compute eigenvectors of a dense non-symmetric matrix like $\mathbf{B}$. In particular, one can exploit parallelism by finding on separate cores the eigenvectors corresponding to different search intervals of eigenvalues$^{42,43}$. These numerical methods thus allow us to handle braids with tens of thousands of strands on a personal computer. Note that one should not expect each relevant eigenvector to correspond to each coherent region: instead, a relevant eigenvector will typically be nearly piecewise constant over multiple coherent regions, see Figures 5 and 7 for examples.

## 2. Visualization.

We visualize eigenvectors of the Burau matrix as follows. We first compute the Voronoi diagram$^{44}$ of the initial positions $\{\mathbf{x}_i(0)\}_{1 \le i \le n}$, which tessellates the domain $\mathcal{D}$ into disjoint convex cells $\{V_i\}_{1 \le i \le n}$. For each eigenvector, we then simply color each cell based $V_j$ on the associated entry of the eigenvector of matrix $\mathbf{B}$. Note that when the initial positions are on a regular grid as we will use in Sec. V, the resulting visualization forms a pixelation of the domain (e.g., see Fig. 5).

## 3. Further processing.

We note that further processing can be done to help a user with finding material coherence. First, many of the eigenvectors found with our methods are quite easy to discard automatically through image processing; see Fig. 3 for examples of eigenvectors that do not depict meaningful material coherence. Second, we can also extract a series of potential coherent regions as follows: since each value of a relevant eigenvector $\mathbf{v}$ corresponds to a value at a particle at the first timestep $t_0$, we can compute the levelsets $L_\theta = \{i \mid |v_i - \theta| < \varepsilon\}$ consisting of all particles whose corresponding $i$th value in $\mathbf{v}$ is approximately $\theta$. These coherent sets can then be observed as
they evolve in time. Finally, boundaries of a coherent sets correspond to the material loops studied by Allshouse and Thiffeault²⁴. We can thus consider the isotopy class of loop(s) that separate $L_\theta$ from the rest of the domain: applying the method described by Allshouse and Thiffeault²⁴, we can track the growth of the material loop via Dynnikov coordinates and verify that the boundaries of the levelsets we found are indeed slow-growing to confirm coherence. Note that in Section V B, we will examine a modified Duffing oscillator and compare the coherent sets found in our method (Figure 7) with the coherent sets found by Allshouse and Thiffeault (Figure 8).

<figure><img src="image_6.png" /><figcaption>FIG. 3. Irrelevant eigenvectors. Eigenvectors of Aref's blinking vortex flow for $\mu = 50$. Most of the eigenvectors of a Burau matrix are nearly constant (left), or only show faint parts of coherent structures (middle/right); here, eigenvector 12, 426, and 574. Compare these results with Fig. 5(bottom right) showing an eigenvector that clearly reveals dynamical structures.</figcaption></figure>

## D. Discussion

Since our numerical probing of coherent sets only involves an eigenvalue problem, our exploration of material coherence from trajectories scales very well. In particular, we can handle braids arising from thousands or tens of thousands of particles on low-end personal computers, far beyond what was reported in Allshouse and Thiffeault²⁴ and Francois et al.²⁷.

Our approach is not without limitations. In particular, no strong guarantee is offered that this eigenanalysis will uncover all relevant coherent sets: our use of trajectory closures and the unfaithfulness of the Burau representation are two of the theoretical hurdles that may prevent a proper detection. Our approach should thus be used as a detection tool to narrow in on potential coherent sets. We demonstrate next that numerical experiments behave remarkably well on practical flows.

## V. Experiments on flows

In the following section, we demonstrate the relevance of our contribution on non-trivial flows. We focus on two different continuous dynamical systems: the blinking vortex flow, and the (modified) Duffing oscillator. Starting from a discretization of the domain $\mathcal{D}$ (a square in all our examples) through particles forming a regular grid, each particle is advected in time by the flow, generating an associated trajectory in space-time. The set of trajectories is then processed as described in Sec. IV, and compared to what is known in the literature.

### A. Blinking vortex flow

The blinking vortex flow was introduced by Aref as an idealization of stirring²³. The flow is given by a pair of vortices separated by a finite distance, blinking on and off periodically in an alternating fashion in an incompressible, inviscid fluid.

We consider a modified version of this flow in an unbounded domain (modeled on the complex plane). The velocity field due to a single point vortex located at $x = a$ on the x-axis is given by

$$ \begin{aligned} \dot{r} &= 0 \\ \dot{\theta} &= \frac{\Gamma}{2\pi r}, \end{aligned} $$

where $\Gamma$ is the strength of the vortex, and $r = \sqrt{(x-a)^2 + y^2}$ is the distance to the center of the vortex. The mapping, in dimensionless form⁴⁵, induced by identical two vortices at $\xi_i = \pm a$, each acting for time $T$, is given by the twist map

$$ \begin{pmatrix} x \\ y \end{pmatrix} \mapsto \begin{pmatrix} \xi_i + (x - \xi_i) \cos \Delta\theta - y \sin \Delta\theta \\ (x - \xi_i) \sin \Delta\theta + y \cos \Delta\theta \end{pmatrix}, $$

where $\Delta\theta = \frac{\mu}{r^2}$, with $\mu = \frac{\Gamma T}{2\pi a^2}$, and $r = \sqrt{(x - \xi_i)^2 + y^2}$. The parameter $\mu$ is the flow strength, and its value controls the behavior of the system. We make distances dimensionless with respect to $a$ and place the vortices at $\xi_i = \pm 1$.

<figure><img src="image_7.png" /><figcaption>FIG. 4. Poincaré sections for Aref's blinking vortex flow. From left to right, top to bottom: $\mu = 1, 5, 10, 20, 30, 35, 40, 50$.</figcaption></figure>

When both vortices act simultaneously ($T = 0, \mu = 0$), the system is integrable. We perturb the system by increasing $\mu$ from zero and check the Poincaré sections $t = kT, k \in \mathbb{Z}$. Chaotic regions appear for all $\mu > 0$⁴⁶. For small values of $\mu$, small chaotic regions exist near the elliptic and hyperbolic points. As $\mu$ increases, the size of the chaotic regions grow, destroying confining KAM surfaces as the chaotic regions merge. See a few Poincaré sections for various values of $\mu$ in Fig. 4. Using our numerical approach, we give eigenvectors for the blinking vortex flow with various values of $\mu$ in Fig. 5 for a set of 900 trajectories of particles, which were initially placed on a regular grid covering the domain. Coherent regions appear, indeed, as isolevels of the eigenvectors as expected.

### B. Modified Duffing oscillator

In this example, we study a modified Duffing oscillator, given by

$$ \begin{aligned} \dot{x} &= y + \alpha \cos(\omega t) \\ \dot{y} &= x(1 - x^2) - dy + \gamma \cos(\omega t) \end{aligned} \qquad (4) $$

with $\alpha = 0.1, \gamma = 0.14, d = 0.08, \omega = 1$. This compressible system is also studied by Allshouse and Thiffeault²⁴ with the same parameters as an example a system with two primary
<figure><img src="image_8.png" /><figcaption>FIG. 5. Burau Eigenvectors for Aref's blinking vortex flow. From left to right, and top to bottom: $\mu = 1, 5, 10, 20, 35, 50$. Nine hundred particles were used, started on a 30 × 30 regular grid.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6. Modified Duffing oscillator. Left: a three-coloring of initial particle positions by their position at time $t = 250$; Right: Phase portrait for trajectories belonging to each of the three regions.</figcaption></figure>

regions of mixing: (i) the limit cycle and its basin of attraction (colored in yellow in Fig. 6(left)) and (ii) the rest of the domain (blue/green).

Trajectories of the modified Duffing oscillator (4) belong to one of three types, each illustrated in Figure 6(right) with the corresponding colors. The yellow region corresponds to the yellow limit cycle and its basin of attraction. The green regions (eventually) attract to one of the two green limit cycles. The blue region attracts to neither the yellow limit cycle nor green limit cycles during the time period studied.

Using the computational method described in this paper, we obtain the Burau eigenvectors for the modified Duffing oscillator depicted in Figure 7. In addition to the two initial conditions found by Allshouse and Thiffeault²⁴, illustrated in Figure 8, the two additional limit cycles (green in Figure 6(left)) are also present.

We note that since the modified Duffing oscillator is a com-

pressible system, particle positions can in fact coincide. In order to form a well-defined braid, we select a sufficiently sparse sampling of the domain and a time window so that the dynamics are reasonably well-captured yet no two particles coincide at any time. Allshouse and Thiffeault²⁴ argue that the ability of the braid theoretic approach to detect coherent sets even in an incompressible flow is a testament to the wide applicability of the method.

<figure><img src="image_10.png" /><figcaption>FIG. 7. Two Burau eigenvectors for modified Duffing oscillator. Two of the eigenvectors of the Burau matrix computed from trajectories of the modified Duffing oscillator (1225 particles); the three main dynamical regions are present as expected.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 8. Previous work. Allshouse and Thiffeault²⁴ detect two types of initial conditions for the modified Duffing oscillator. The dots (numbered from left to right) are the initial conditions for the trajectories that are studied further by Allshouse and Thiffeault as representative trajectories for the two types of initial conditions.</figcaption></figure>

# VI. Conclusions

Our approach based on topological braids is especially advantageous when the input data are sparse, since it does not require nearby trajectories or derivatives of the velocity field. However, the braid approach is not without limitations. Accuracy is limited by the length of trajectory histories, and braid methods are only adapted to time intervals where particles never intersect at the observation scale. Fortunately, many geophysical flows, when observed at sufficiently large scale, are nearly incompressible, rendering this last limitation only mildly restrictive.

## 1. Faithfulness.

In this work, we have chosen to use the Burau representation for the analysis of spatially sparse particle trajectories. The Burau representation is not faithful but confers computational advantages (*n* vs. *n*²) in both space and time requirements compared to the faithful Lawrence-Krammer
representation⁴⁷. If a greater level of topological and dynamical detail is preferred, one may wish to consider the Lawrence-Krammer representation instead.

## 2. Parallelism.

The algorithm that we describe in this chapter lends itself naturally to parallelism. We do not discuss the details here, but we note that the matrix chain multiplication required by our algorithm is amenable to optimization. In particular, in addition to the usual considerations of the matrix chain multiplication problem⁴⁸,⁴⁹, we remark that by taking the sequence of braid generators into account, then depending on the flow, we can potentially partition the braid $\beta = \prod_{l}^{L} \sigma_{b_l}$ into sub-sequences of neighboring generators (e.g., each subsequence consists only of generators $\sigma_{i-j}$, $\sigma_{i-j+1}$, ..., $\sigma_{i+j}$, for small $j$). This effectively partitions the sequence of matrix multiplications into subsequences that each consist of (mostly) sparse matrix multiplications, thereby reducing the computational requirements of our analysis. Future work may wish to formalize these and other computational considerations: such a computational optimization may be crucial, in particular if a $(n^2 \times n^2)$ Lawrence-Krammer representation is used.

## Acknowledgments

MD gratefully acknowledges the Inria International Chair program during which he started this work, and all the members of the TITANE team for their support.

<footnote n="1">¹C. H. Amon, A. M. Guzmán, and B. Morel, “Lagrangian chaos, Eulerian chaos, and mixing enhancement in converging-diverging channel flows,” Physics of Fluids **8**, 1192–1206 (1996).</footnote>

<footnote n="2">²G. Haller, “Lagrangian coherent structures,” Annual Review of Fluid Mechanics **47** (2015), 10.1146/annurev-fluid-010313-141322.</footnote>

<footnote n="3">³M. R. Allshouse and T. Peacock, “Lagrangian-based methods for coherent structure detection,” Chaos: An Interdisciplinary Journal of Nonlinear Science **25**, 097617 (2015).</footnote>

<footnote n="4">⁴J. Weiss, “The dynamics of entrophy transfer in two-dimensional hydrodynamics,” Physica D **48**, 273–294 (1991).</footnote>

<footnote n="5">⁵H. Yang, “Chaotic transport and mixing by ocean gyre circulation,” in *Stochastic Modelling in Physical Oceanography*, edited by R. J. Adler, P. Müller, and B. L. Rozovskii (Birkhäuser, 1996) pp. 439–466.</footnote>

<footnote n="6">⁶G. Haller and G.-C. Yuan, “Lagrangian coherent structures and mixing in two-dimensional turbulence,” Physica D: Nonlinear Phenomena **147** (2000), 10.1016/S0167-2789(00)00142-1.</footnote>

<footnote n="7">⁷G. Haller, “Lagrangian coherent structures from approximate velocity data,” Physics of Fluids **14**, 1851–1861 (2002).</footnote>

<footnote n="8">⁸S. C. Shadden, F. Lekien, and J. E. Marsden, “Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional,” Physica D **212**, 271–304 (2005).</footnote>

<footnote n="9">⁹G. A. Voth, G. Haller, and J. P. Gollub, “Experimental measurements of stretching fields in fluid mixing,” Phys. Rev. Lett. **88**, 254501 (2002).</footnote>

<footnote n="10">¹⁰M. Mathur, G. Haller, T. Peacock, J. E. Ruppert-Felsot, and H. L. Swinney, “Uncovering the Lagrangian skeleton of turbulence,” Phys. Rev. Lett. **98**, 144502 (2007).</footnote>

<footnote n="11">¹¹J. Kasten, C. Petz, I. Hotz, H.-C. Hege, B. R. Noack, and G. Tadmor, “Lagrangian feature extraction of the cylinder wake,” Physics of Fluids **22**, 091108 (2010).</footnote>

<footnote n="12">¹²F. Lekien and S. D. Ross, “The computation of finite-time Lyapunov exponents on unstructured meshes and for non-euclidean manifolds,” Chaos **20**, 017505 (2010).</footnote>

<footnote n="13">¹³G. Haller, “A variational theory of hyperbolic Lagrangian coherent structures,” Physica D **240**, 574–598 (2011).</footnote>

<footnote n="14">¹⁴M. Farazmand and G. Haller, “Computing Lagrangian coherent structures from their variational theory,” Chaos **22**, 013128 (2012).</footnote>

<footnote n="15">¹⁵G. Haller and F. J. Beron-Vera, “Geodesic theory of transport barriers in two-dimensional flows,” Physica D: Nonlinear Phenomena **241**, 1680 – 1702 (2012).</footnote>

<footnote n="16">¹⁶M. Budišić and I. Mezić, “Geometry of the ergodic quotient reveals coherent structures in flows,” Physica D: Nonlinear Phenomena **241**, 1255–1269 (2012).</footnote>

<footnote n="17">¹⁷I. Mezić, “Analysis of fluid flows via spectral properties of the Koopman operator,” Annual Review of Fluid Mechanics **45**, 357–378 (2013).</footnote>

<footnote n="18">¹⁸G. Froyland and K. Padberg-Gehle, “Almost-invariant and finite-time coherent sets: Directionality, duration, and diffusion,” in *Ergodic Theory, Open Dynamics, and Coherent Structures*, edited by W. Bahsoun, C. Bose, and G. Froyland (Springer, 2014) pp. 171–216.</footnote>

<footnote n="19">¹⁹J.-L. Thiffeault, “Braids of entangled particle trajectories,” Chaos **20**, 017516 (2010).</footnote>

<footnote n="20">²⁰A. Borel, “Jean Leray and algebraic topology,” in *Selected Papers - Oeuvres Scientifiques*, Springer Collected Works in Mathematics, Vol. 1 (Springer-Verlag Berlin Heidelberg, 1998).</footnote>

<footnote n="21">²¹V. I. Arnold, “Sur la géométrie différentielle des groupes de lie de dimension infinie et ses applications à l’hydrodynamique des fluides parfaits,” Ann. Inst. Fourier **16**, 319–361 (1966).</footnote>

<footnote n="22">²²J.-L. Thiffeault, “Measuring topological chaos,” Phys. Rev. Lett. **94**, 084502 (2005).</footnote>

<footnote n="23">²³E. Artin, “Theory of braids,” Ann. of Math. **48**, 101–126 (1947).</footnote>

<footnote n="24">²⁴M. R. Allshouse and J.-L. Thiffeault, “Detecting coherent structures using braids,” Physica D: Nonlinear Phenomena **241**, 95–105 (2012).</footnote>

<footnote n="25">²⁵J.-L. Thiffeault and M. Budišić, “BRAIDLAB: a software package for braids and loops,” arXiv preprint arXiv:1410.0849 (2015).</footnote>

<footnote n="26">²⁶M. Budišić and J.-L. Thiffeault, “Finite-time braiding exponents,” Chaos: An Interdisciplinary Journal of Nonlinear Science **25**, 087407 (2015).</footnote>

<footnote n="27">²⁷N. Francois, H. Xia, H. Punzmann, B. Faber, and M. Shats, “Braid entropy of two-dimensional turbulence,” Sci. Rep. **5**, Art. 18564 (2015).</footnote>

<footnote n="28">²⁸Note that Burau matrices have been used before in, e.g.,³⁹, as a computational tool—but only to estimate braid entropy.</footnote>

<footnote n="29">²⁹H. Aref, “Stirring by chaotic advection,” J. Fluid Mech. **143**, 1–21 (1984).</footnote>

<footnote n="30">³⁰P. Boyland, “Topological methods in surface dynamics,” Topology Appl. **58**, 223–298 (1994).</footnote>

<footnote n="31">³¹B. Farb and D. Margalit, *A Primer on Mapping Class Groups*, Vol. PMS-49 (Princeton University Press, 2012).</footnote>

<footnote n="32">³²D. Bernardete, M. Gutierrez, and Z. Nitecki, “A combinatorial approach to reducibility of mapping classes,” in *Mapping Class Groups and Moduli Spaces of Riemann Surfaces*, Contemporary Mathematics, Vol. 150, edited by C.-F. Bodigheimer and R. M. Hain (Amer. Math. Soc., 1993) pp. 1–31.</footnote>

<footnote n="33">³³D. Bernardete, Z. Nitecki, and M. Gutierrez, “Braids and the nielsen-thurston classification,” J. Knot Theory Ramifications **4** (1995), 10.1142/S0218216595000259.</footnote>

<footnote n="34">³⁴J. S. Birman and T. E. Brendle, “Braids: A survey,” in *Handbook of Knot Theory*, edited by W. Menasco and M. Thistlethwaite (Elsevier, 2005).</footnote>

<footnote n="35">³⁵W. Burau, “Über zopfgruppen und gleichsinnig verdrillte verkettungen,” Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg **11**, 179–186 (1935).</footnote>

<footnote n="36">³⁶D. Long and M. Paton, “The Burau representation is not faithful for $n \ge 6$,” Topology **32**, 439–447 (1993).</footnote>

<footnote n="37">³⁷S. Bigelow, “The Burau representation is not faithful for $n = 5$,” Geometry & Topology **3**, 397–404 (1999).</footnote>

<footnote n="38">³⁸T. Church and B. Farb, “Infinite generation of the kernels of the Magnus and Burau representations,” Algebraic & Geometric Topology **10**, 837–851 (2010).</footnote>

<footnote n="39">³⁹G. Band and P. Boyland, “The Burau estimate for the entropy of a braid,” Algebraic & Geometric Topology **7**, 1345–1378 (2007).</footnote>

<footnote n="40">⁴⁰C. C. Squier, “The Burau representation is unitary,” Proc. Amer. Math. Soc. **90**, 199–202 (1984).</footnote>

<footnote n="41">⁴¹J. González-Meneses, “The nth root of a braid is unique up to conjugacy,” Algebraic & Geometric Topology **3**, 1103–1118 (2003).</footnote>

<footnote n="42">⁴²M. Gates, A. Haidar, and J. Dongarra, “Accelerating computation of eigen-vectors in the dense nonsymmetric eigenvalue problem,” in *High Performance Computing for Computational Science*, edited by M. e. a. Daydé (Springer International Publishing, 2015) pp. 182–191.</footnote>

<footnote n="43">⁴³J. Kestyn, E. Polizzi, and P. Tang, “FEAST eigensolver for non-Hermitian problems,” Tech. Rep. (arXiv:1506.04463, 2015).</footnote>
<footnote n="44">⁴⁴J.-D. Boissonnat and M. Yvinec, *Algorithmic geometry* (Cambridge University Press, 1998).</footnote>

<footnote n="45">⁴⁵J. Ottino, *The kinematics of mixing: stretching, chaos, and transport* (Cambridge University Press, 1989).</footnote>

<footnote n="46">⁴⁶M. F. Doherty and J. M. Ottino, “Chaos in deterministic systems: strange attractors, turbulence, and applications in chemical engineering,” *Chemical Engineering Science* **43**, 139–183 (1988).</footnote>

<footnote n="47">⁴⁷S. Bigelow, “The Lawrence-Krammer representation,” in *Topology and geometry of manifolds* (Proc. Sympos. Pure Math.), Vol. 71, edited by</footnote>

<footnote>G. Matić and C. McCrory (2003).</footnote>

<footnote n="48">⁴⁸T. C. Hu and M. T. Shing, “Computation of matrix chain products. part I,” *SIAM Journal on Computing* **11**, 362–373 (1982).</footnote>

<footnote n="49">⁴⁹T. C. Hu and M. T. Shing, “Computation of matrix chain products. part II,” *SIAM J. Comput.* **13**, 228–251 (1984).</footnote>