Correlations between the leading Lyapunov vector and pattern defects for chaotic Rayleigh-Bénard convection

Many important challenges facing science and engineering today can be described as a spatially extended system that is driven far-from-equilibrium.¹ Examples include the complex dynamics of the weather and climate,² the patterns of interacting chemicals in biological systems,^3–5 chaotic electroconvection,^6,7 the striking dynamics of nonlinear optics,^8–10 and the formation of uniform crystal structures from a melt.¹¹ In many cases of interest, the systems are strongly driven, the dynamics are highly nonlinear, and a model description (if available) is very complex and difficult to evaluate. Furthermore, even if the model description is tractable, many of the important theoretical ideas and diagnostics are very difficult to compute.

One significant challenge is that the dimension of the dynamics of large, strongly driven systems is expected to be very large. For example, for many fluid systems under laboratory conditions, the dimension of the dynamics is expected to be on the order of hundreds or thousands.^12–14 If the equations describing the dynamics are known, a powerful approach to probe the dynamics is through the computation of Lyapunov vectors and exponents.^15,16 Using current algorithms and computing resources, it is now possible to compute Lyapunov diagnostics for laboratory scale systems.

Typically, the growth or decay of the Lyapunov vectors is used to determine the exponential growth or decay rates of infinitesimal perturbations in the tangent space which are the Lyapunov exponents. The spectrum of Lyapunov exponents represents the infinite-time average of the exponential growth or decay rate of the Lyapunov vectors. The leading Lyapunov exponent represents the largest exponential growth rate. Therefore, the leading Lyapunov vector represents the direction of a perturbation in a very high dimensional state-space that would grow exponentially with the maximum growth rate given by the leading Lyapunov exponent on average. From the spectrum of Lyapunov exponents, it is also possible to estimate the fractal dimension of the chaotic attractor in state space.

However, it is also possible to gain insight into the spatial and temporal features of the dynamics using the finite-time or instantaneous Lyapunov exponents and vectors. Rather than an infinite-time average, the Lyapunov exponents and vectors are computed for a finite duration of time. In this paper, we will focus only on the finite-time leading Lyapunov vector. The spatial features of the finite-time leading Lyapunov vector indicates the sensitivity of those regions to a perturbation during that period of time. In particular, regions where the magnitude of the leading Lyapunov are large indicate regions that are very sensitive to small perturbations. This can be used to quantitatively explore the correlations between the flow field topology and dynamics with the leading Lyapunov vector. This is precisely the approach we have taken here.

However, many important systems are beyond the current reach of numerics, yet it is possible to make detailed measurements of the dynamics. This often yields large amounts of data describing the spatiotemporal patterns of quantities of interest. For example, for a fluid system, this could be many detailed images of the time variation of the flow field patterns. For a chemical system, this could be images of the spatiotemporal variation of the concentration of reactants and products. In geophysical problems, one could imagine using detailed satellite imagery of plankton blooms in the oceans or detailed weather dynamics as sources of data.

It would be very insightful to have a fundamental understanding of the connection between the pattern dynamics that can be measured in the experiment and the deep insights that one gains from knowledge of the Lyapunov vectors and exponents. In this article, we discuss our efforts to address this difficult problem using the canonical pattern-forming system of Rayleigh–Bénard convection.

For several reasons, Rayleigh–Bénard convection presents an opportunity to study the fundamental connections between powerful ideas from dynamical systems theory with experimentally accessible flow field dynamics. Firstly, it is a pattern-forming system with a rich literature of deep insights from detailed experimental and theoretical studies.^1,17 Using modern algorithms and supercomputing resources, it is now possible to numerically simulate Rayleigh–Bénard convection for the precise conditions of the experiment.¹⁸ It is also possible to compute the Lyapunov vectors and exponents for the precise conditions of experiment.^14,19–21

There is a large body of research describing the rich dynamics of chaotic Rayleigh–Bénard convection.¹ The role of topological features in the dynamics has played an important part in building a physical understanding of the complex dynamics.^22–25 The leading Lyapunov vector has been computed and its spatiotemporal dynamics have been found to be highly localized in space and time.^12,19,21,26 The regions where the leading Lyapunov vector is large have been found to correlate with the presence of defect structures such as roll pinch-off events, targets, and spiral structures. The regions where the leading Lyapunov vector is large correspond to the regions in the flow field that are extremely sensitive to small perturbations. For example, near the defect structures identified by the leading Lyapunov vector, one would expect that small variations in the flow field properties would have a large effect on the dynamics such as pinch-off locations or reconnection events that could lead to large scale pattern changes later in time. Conversely, small perturbations to the flow field in a region where the Lyapunov vector has a small magnitude would be expected to quickly dissipate resulting in little effect on the pattern dynamics.

However, the picture that has emerged from these studies is more subtle than what was perhaps expected. Although the regions in the flow field where the magnitude of the Lyapunov vector is large and localized most often occur near a topological defect there are many instances where the topological feature is present and the magnitude of the leading Lyapunov vector is small at that location. In other words, the presence of a topological feature does not guarantee that the Lyapunov vector will be large. This was clearly demonstrated by Scheel and Cross¹⁹ in their study of the dynamics of Rayleigh–Bénard convection in a small cylindrical domain which showed that the time history of roll pinch-off events played an important role in their contribution to the long-time Lyapunov exponent.

In this paper, we present a careful and methodical exploration of this situation: namely, is it possible to build a quantitative link between the features of the flow field pattern and the spatial variation of the magnitude of the leading Lyapunov vector? First, there is the possibility that one simply has not used the correct topological feature to make this quantitative link. Second, and this is more likely to be the case, it is possible that one must include some explicit notion of the flow history in the diagnostic to improve this connection. We explore both of these ideas in depth.

A key element of our study is the availability of a large amount of accurate numerical data for chaotic three-dimensional convection for experimental conditions. It is essential that our approach be automated in the sense that thousands of images can be quantified without any subjective elements to the diagnostics. This will allow us to explore the statistical correlations between features of the flow field and the leading Lyapunov vector. We will use the ideas of precision and recall to probe the statistical efficacy of the different pattern diagnostics. Roughly speaking, the precision of a diagnostic indicates the probability that it will identify a region where the Lyapunov vector is large and the recall indicates what percentage of the regions with a large magnitude of the Lyapunov vector it is able to capture. We will find that the interplay between precision and recall for the different diagnostics we explore is able to provide interesting insights that would have been difficult to ascertain otherwise.

The article is organized as follows. In Sec. II, we briefly discuss Rayleigh–Bénard convection, the Boussinesq equations, and our numerical approach for computing the flow field and the Lyapunov vectors and exponents. In Sec. III, we describe our procedure for exploring the statistical correlations between the flow field topology and dynamics with the spatial variation of the leading Lyapunov vector. We use the ideas of precision and recall to quantify the usefulness of a wide variety of pattern diagnostics in terms of their ability to indicate where the leading Lyapunov vector magnitude will be significant. Lastly, in Sec. IV, we present some concluding remarks.

Rayleigh–Bénard convection occurs when a shallow layer of fluid is uniformly heated from below in a gravitational field. A schematic of the cylindrical convection domain we use in our study is shown in Fig. 1. The radius of the cylindrical domain is $r0$⁠, the depth of the fluid layer is $d$⁠, and the convention for the Cartesian coordinates $(x,y,z)$ are shown, where $z$ is opposing to the direction of gravity $g$⁠.

As the temperature difference between the bottom and top surfaces is increased, the quiescent fluid layer undergoes a supercritical bifurcation to convective motion. This temperature difference can be represented as the Rayleigh number $R$ which often acts as the control parameter. As the Rayleigh number is increased beyond its critical value $R>Rc$⁠, the convection rolls themselves become dynamic. Typically, the convection rolls become time dependent, chaotic, and eventually the rolls become unstable to smaller-scale structures, such as plumes, and the flow field becomes turbulent.

Rayleigh–Bénard convection contains all of the essential complexity and difficulties of spatially extended systems driven far-from-equilibrium, yet it is experimentally accessible. With today’s computing resources, it is also amenable to long-time numerical simulations for experimentally realistic conditions.

For moderate Rayleigh number, the Boussinesq equations describe the fluid motion. In the nondimensional form, these are

which represent the conservation of momentum, energy, and mass, respectively. In our notation, $u→(x,y,z,t)$ is the velocity vector, $p(x,y,z,t)$ is the pressure, $T(x,y,z,t)$ is the temperature, $t$ is time, $(x,y,z)$ are Cartesian coordinates, and $z^$ is a unit vector in the $z$ direction which opposes gravity. We have followed the typical convention and used the depth $d$ of the fluid layer as the length scale, the constant temperature difference $ΔT$ between the bottom and top surfaces as the temperature scale, and the vertical thermal diffusion time for heat $d2/α$ as the time scale, where $α$ is the thermal diffusivity of the fluid. The Prandtl number $σ$ is the ratio of the diffusivities of momentum and heat. Finally, the extent of the cylindrical domain is given by the aspect ratio $Γ=r0/d$⁠.

The boundary conditions are no-slip at all surfaces in contact with the fluid. The temperature of the hot bottom surface is $T(z=0)=1$⁠, and the temperature of the cold top surface is $T(z=1)=0$⁠. The lateral sidewalls of the domain are composed of a perfectly conducting material.

We also compute the spectrum of Lyapunov vectors and exponents. In the following, we provide only the essential ideas behind this computation and we refer the reader to Ref. 14 for a more detailed discussion. In order to compute the $Nλ$ largest Lyapunov vectors, we simultaneously evolve $Nλ$ copies of the tangent space equations. The tangent space equations are

\begin{aligngroup}

where $k=1,…,Nλ$⁠. The variables $δu→(x,y,z,t)(k),δp(x,y,z,t)(k)$⁠, and $δT(x,y,z,t)(k)$ are the $kth$ perturbations about the nonlinear orbit through state space given by $u→(x,y,z,t)$⁠, $p(x,y,z,t)$⁠, and $T(x,y,z,t)$⁠.

Using these perturbations, one can represent the $kth$ Lyapunov vector at time $t$ as the large column vector given by

where $(δu(k),δv(k),δw(k))$ are the $(x,y,z)$ components of the $kth$ perturbation velocity field $δu→k$⁠, and the superscript $T$ indicates a transpose. Equation (7) is meant to convey that the Lyapunov vector $v→g(k)(t)$ is represented as all of the values of $δu(k)(t)$ at time $t$ as one large row vector, followed by all of the components of $δv(k)(t)$ as one large row vector, and so on for $δw(k)(t)$ and $δT(k)(t)$⁠. The perturbation of the pressure field $δp(t)(k)$ is not included here since it is not an independent variable for incompressible flow.

In practice, the spectrum of Lyapunov vectors $v→g(k)$ is periodically reorthonormalized using a Gram–Schmidt procedure to avoid numerical errors associated with the tendency that all of the vectors will point toward the direction of fastest growth in the tangent space. The subscript $g$ on $v→g(k)$ indicates that these are the Gram–Schmidt orthogonalized vectors. As a result, the only Lyapunov vector pointing in a physically important direction will be the leading Lyapunov vector $v→g(1)(t)$⁠.

However, the Gram–Schmidt reorthonormalization is volume preserving in the tangent space, and, therefore, one can use these vectors to correctly compute the spectrum of Lyapunov exponents. Each time a reorthogonalization is performed the exponential growth or decay of each Lyapunov vector during that window of time is computed to yield the finite time or instantaneous Lyapunov exponent $λ~k(t)$⁠. The large fluctuations of the instantaneous Lyapunov exponents reflect the dynamics that has occurred for that period of time. The infinite time average of $λ~k(t)$ would yield the Lyapunov exponents $λk$⁠.

The Lyapunov exponents $λk$ are guaranteed to be in the order $λ1≥λ2≥⋯λNλ$⁠. Given the spectrum of Lyapunov exponents $λk$⁠, it is then possible to compute the fractal dimension $Dλ$ of the dynamics using the Kaplan–Yorke formula²⁷ 

where $K$ is the largest integer $n$ such that $∑i=1nλi>0$⁠. The fractal dimension $Dλ$ is the most common dimension computed for high-dimensional systems and provides an estimate for the lower bound of the number of chaotic degrees of freedom that would be necessary to describe the dynamics on average.²⁸ 

In order to numerically integrate Eqs. (1)–(3) and the $Nλ$ copies of Eqs. (4)–(6), we use a highly efficient parallel spectral element approach²⁹ and its implementation in the open source solver nek5000.³⁰ The approach is third-order accurate in time, exponentially convergent in space, and has been used broadly to study a wide range of problems in fluid dynamics.³¹ A more detailed description of the numerical approach as it pertains to explorations of Rayleigh–Bénard convection can be found in Refs. 14 and 18–21.

In this paper, we focus our attention on the spatiotemporal features of the leading Lyapunov vector $v→g(1)(t)$ and its statistical relationship to the flow field patterns and dynamics. In order to probe this quantitatively, we use the temperature perturbation field at the horizontal midplane $δT(x,y,z=1/2,t)(1)$ as a representation of the Lyapunov vector.²⁶ This allows us to identify the regions in physical space where the divergence of the Lyapunov vector is large which we can statistically correlate with an experimentally accessible quantity characterizing the flow field pattern such as the temperature field.

The approach we have taken is to explore ideas from pattern diagnostics and computational homology^32,33 for use on the temperature field in order to build an understanding of their statistical connection with the spatiotemporal variation of the leading Lyapunov vector. This understanding provides physical insights into the complex dynamics of high-dimensional systems driven far-from-equilibrium. For example, it may lead to a description of which topological structures in the flow field are contributing most to the dynamics or perhaps a suggestion of what a modal decomposition of the dynamics may look like. These are difficult and open challenges. In this paper, we have pushed these ideas further for the canonical problem of Rayleigh–Bénard convection.

We have chosen to explore convection in a cylindrical domain with an aspect ratio of $Γ=20.4$ containing a fluid with a Prandtl number $σ=0.84$ and for a Rayleigh number of $R=4000$ (see Fig. 1). These parameters are chosen in order to generate a high-dimensional chaotic state that is accessible to our numerical approach while also being in a regime that would be accessible to experiments using compressed gases.³³ 

The initial condition for our numerical simulation was a no-flow state with small, random perturbations to the temperature field. We then integrated Eqs. (1)–(3) for a long time to allow all initial transients to decay. This initial simulation was conducted for a time of $t=800$ which is approximately equal to the two horizontal diffusion times $τh$ for heat. The quantity $τh$ is the amount of time it would take for heat to diffuse from the center of the domain to the sidewall, where $τh=Γ2$⁠. The horizontal diffusion time for heat is expected to provide a useful estimate of the time required for initial transients to decay.^14,34

We next numerically integrated Eqs. (1)–(3) along with $Nλ=60$ copies of Eqs. (4)–(6) for another 100 time units. During this time, we performed periodic Gram–Schmidt reorthonormalizations and computed the Lyapunov vectors and exponents to generate the Lyapunov exponents and fractal dimension shown. We then continued the simulation for only the leading Lyapunov vector for another 350 time units. The time step in our simulation was $Δt=0.001$⁠, and we performed the reorthonormalizations every 10 time steps. This yielded 3500 images of the flow field and the leading Lyapunov vector that we have used in our analysis.

An image of a typical flow field pattern from our numerical simulation is shown in Fig. 2(a). The image shows a horizontal midplane slice (⁠$z=1/2$⁠) of the temperature field, where light regions indicate hot rising fluid and dark regions indicate cool falling fluid. The flow field pattern is composed of convection rolls and a variety of defect structures. Figure 2(b) shows the magnitude of the temperature perturbation field corresponding to the leading Lyapunov vector at the same horizontal midplane slice (⁠$z=1/2$⁠). We are interested in exploring the relationship between the high-magnitude regions of the leading Lyapunov vector and the topology and dynamics of the flow field pattern. However, one of our goals is to translate the results from this work from the numerical to experimental setting. Thus, the image shown in Fig. 2(a) is derived by projecting the numerical results onto a standard 8-bit gray scale filtration that corresponds to an experimentally accessible image. To simplify the presentation, this gray scale filtration is itself rescaled to the unit interval. In the analysis that follows we will find it useful to also morphologically dilate the leading Lyapunov vector magnitude function using a disk of one unit length as shown in Fig. 2(c). This will be discussed further in Sec. III A.

Figure 3(a) illustrates the spectrum of $Nλ=60$ Lyapunov exponents. The leading Lyapunov exponent is positive $λ1>0$⁠, which indicates that the dynamics are chaotic, as expected. Using Eq. (8), the fractal dimension $Dλ$ can be computed from the spectrum of Lyapunov exponents $λk$⁠. Figure 3(b) illustrates the time variation of $Dλ$ as it approaches an asymptotic value of $Dλ≈55$⁠. This indicates that the dynamics has approximately 55 active chaotic degrees of freedom, on average. The results shown in Figs. 3(a) and 3(b) have been computed for 100 time units and would be expected to continue this slow convergence if the computations were extended.

In Fig. 3(c), we show the time variation of the instantaneous leading Lyapunov exponent $λ~1(t)$ and our calculation of the leading Lyapunov exponent $λ1(t)$ for 100 time units. In our computations, $λ~1(t)$ is computed every 10 time steps and the large fluctuations in its value represents the growth or decay of the leading finite-time Lyapunov vector during that period of time. The leading Lyapunov exponent $λ1(t)$ is the running time average of $λ~1(t)$ which is converging toward a value of $λ1(t)≈0.57$ which is the value shown in Fig. 3(a) for $k=1$⁠. The leading Lyapunov vector that we analyze, as shown in Fig. 2(b), is a spatiotemporal representation of the fluctuations of $λ~1(t)$⁠.

Regions where the leading Lyapunov vector has high magnitude correspond to regions in the temperature field that exhibit a high degree of sensitivity to local perturbations. State-of-the-art methods for computing the leading Lyapunov vector apply only when one has full knowledge of the underlying dynamics, which we do not possess for experimental systems. Thus, our goal is to find experimentally accessible features that can be used to predict high-magnitude regions of the leading Lyapunov vector. For example, in Fig. 2(b), we would like to be able to predict the spatial locations where the Lyapunov vector magnitude is significant (not dark blue) using only information from the pattern flow field that has been computed from images such as what is shown in Fig. 2(a). There are many possible candidates to consider that may provide insight between the flow field patterns and the Lyapunov vector. We have studied a variety of candidates that we separate into two groups.

In the first group, which we will refer to as group I pattern diagnostics, we consider many of the pattern diagnostics that are typically used when studying convection patterns. These are the spatial variation of the local wavenumber of the convection rolls and the presence of topological defects. Furthermore, we quantify the role of several particular types of topological defects including spirals, targets, grain boundaries, and disclinations.

In the second group, or group II pattern diagnostics, we include diagnostics that contain additional information that may be important indicators of where the Lyapunov vector is significant. For example, the creation or annihilation of certain topological features as opposed to just their presence in the pattern. Specifically, we have explored temporal derivatives of the temperature field, the presence of roll pinch-off events, the emergence of target structures, and the dynamics of topological defects in the patterns.

Our approach to probe the ability of these different pattern diagnostics to indicate the spatial variation of the leading Lyapunov vector is the following. For each diagnostic, we define a pattern diagnostic function $F(x,y,t)$ that is based on the temperature field at the horizontal midplane $z=1/2$ at time $t$⁠. We use a threshold such that $F(x,y,t)$ is a binary valued function, where

The details behind each diagnostic function $F(x,y,t)$ is unique to the pattern diagnostic being studied, and we describe the specifics behind the individual computations in Secs. III B and III D. In all cases, the diagnostic function $F(x,y,t)$ is essentially an indicator function that identifies the presence or absence of a particular pattern feature for each pixel in the discretized domain.

We then seek to quantify the extent to which the pattern features indicated by $F(x,y,t)$ are correlated with the spatial variation of the leading Lyapunov vector at time $t$⁠. Figure 2(b) illustrates that the spatial variation of the Lyapunov vector contains both localized and long-range structures. In a typical image of the leading Lyapunov vector, there are several localized maxima present which are indicated by the red regions which occur on a length scale of a single convection roll (approximately a nondimensional distance of unity). However, there are also structures on a length scale of several convection rolls which are indicated by the light blue, green, and yellow regions. The interplay between these two length scales is dynamic and is not fully understood.²⁶ 

In order to quantify the spatial variation of the leading Lyapunov vector, we define another binary valued function $L(x,y,t;α),$ where

The parameter $α$ is a threshold that is used to determine how much of the Lyapunov vector we include in our analysis, where $α∈[0,1]$⁠. When $x,y,$ and $t$ are understood, we will write $L(α)$ instead of $L(x,y,t;α)$⁠. A large value of $α$ would only include the localized peak structures in the Lyapunov vector. As the value of $α$ is decreased, more of the larger-scale structure of the Lyapunov vector is also included. A value of $α=0$ would simply include the entire Lyapunov vector field.

We are interested is quantifying the correlations between the pattern diagnostic functions and the leading Lyapunov vector. Our approach is to build a statistical description using the binary valued functions representing $F(x,y,t)$ and $L(x,y,t;α)$ while varying the value of the threshold $α$⁠. One advantage of using binary valued functions for this purpose is that it can be efficiently applied to experimental results that have been digitized. In addition, the essential computation in building the statistics will become an efficient binary classification task. Overall, we anticipate that this approach should be straightforward to implement for a new problem and that it will readily scale to large amounts of data.

Since the proportion of the domain corresponding to high-magnitude regions of the Lyapunov vector is small, we perform this comparison of $F$ and $L$ using precision-recall curves,^35,36 where the variational parameter is the threshold $α$⁠. Since we do not expect to have a perfect predictor $F$ for the Lyapunov vector magnitude function $L$⁠, we also morphologically dilate $L$ by a disk with a diameter of one unit length which serves to increase the area accounted for by the local maxima³⁷ [see Fig. 2(c)]. Our use of morphological dilation is described in detail in  Appendix A.

Consider a chaotic flow field at time $t$ where we have computed $F(x,y,t)$ and $L(x,y,t;α)$⁠. The precision $Pα,t$ of the diagnostic function $F$ relative to the Lyapunov magnitude function $L$ with a threshold value of $α$ at time $t$ is defined as the conditional probability given by

In this definition, the probability $P$ is computed using the values of $F$ and $L(α)$ over all of the spatial points at time $t$⁠. The precision $Pα,t$ is the probability that the Lyapunov magnitude function $L(α)=1$ at some location in space given that the diagnostic indicator function $F=1$ at this location for a flow field at time $t$⁠.

The precision of a diagnostic function is a measure of its predictive value. For example, consider a diagnostic function $F$ that identifies $N$ locations in a flow field image where $F=1$⁠, while only $M$ of these locations also yield $L(α)=1$ for a particular value of the threshold $α$⁠. In this case, the precision of the diagnostic function is $Pα=M/N$⁠. Of the regions identified by $F$⁠, the precision is the probability that these identified regions correspond to a region where $L(α)=1$⁠.

Similarly, the recall $Rα,t$ of $F$ relative to $L$ at the threshold value $α$ is defined as

The recall $Rα,t$ is the probability that the diagnostic indicator function yields $F=1$ given that the Lyapunov magnitude function $L(α)=1$ at a particular point in space for a flow field at time $t$⁠. Consider again our example where a diagnostic function $F$ has identified $N$ locations where $F=1$ in a flow field image, where $M$ of these points also yield $L(α)=1$⁠. In this case, we also use the fact that there are a total of $NT$ locations identified such that $L(α)=1$⁠, where $NT≥M$⁠. The recall of the diagnostic function is then $Rα=M/NT$⁠. Therefore, the recall is the fraction of the total number of regions, where $L(α)=1$ that are identified by the diagnostic function $F$⁠.

A large value of the recall $Rα,t$ indicates that a large percentage of the regions where $L(α)=1$ are captured by the diagnostic function $F$⁠. However, one could imagine a diagnostic, where $F=1$ nearly everywhere on the pattern. This would yield a large value for the recall, yet such a diagnostic would clearly not be very useful and would, in our case, have low precision for $α>0$⁠. In practice, knowledge about both precision and recall of the diagnostic is what is insightful.

Ideally, one would like a diagnostic function with perfect precision and recall where $Pα,t(F,L)=Rα,t(F,L)=1$ for every time $t$⁠. This would only be achieved if regions where $F$ and $L(α)$ are unity exactly coincide for a flow field image at time $t$ using a particular value of $α$⁠. However, this is never achieved for the patterns we explore and, instead, we find a range of precision and recall values for a particular diagnostic function that vary as a function of $α$⁠.

We are interested in using these ideas to quantify a statistical relationship between a particular pattern diagnostic and the magnitude of the leading Lyapunov vector. In the following, we compute the variation of the precision and recall with the value of the threshold $α$ for many flow field images from the time series described above. We find it useful to compute the time average of the precision and recall over this time series which we denote as $P¯α(F,L)$ and $R¯α(F,L)$⁠, respectively. For our data, we have a total of 3500 flow field images separated in time by 0.1 time unit. We compute the time averaged precision and recall over the range of threshold $α∈[0,1]$ in increments of $Δα=0.05$⁠.

Since we are concerned with measuring how well a diagnostic predicts the regions of high magnitude in the Lyapunov vector, we do not include any time points $t$ in the average where that particular diagnostic returns $F=0$ for all spatial locations. We compare the efficacy of each diagnostic by studying the average precision-recall curves over the sampled time series by plotting the variation of $P¯α$ and $R¯α$ with the threshold $α$⁠.

The patterns of Rayleigh–Bénard convection are composed of counter-rotating convection rolls and defect structures.^1,38 We have developed an automated approach to compute the pattern diagnostic functions for topological defects. Our computational approach is briefly outlined as follows. We first identify several types of topological point defects present in the pattern using the algorithm of Bazen and Gerez.³⁹ In particular, we determine the presence of spiral structures, target structures, grain boundaries, dislocation pairs, convex disclinations, and concave disclinations. Sample images of these different topological defects from our numerical simulation are shown in Fig. 4. We accomplish this classification by clustering the computed topological point defects using single-linkage hierarchical clustering with a radius of $r=1$⁠. These clusters of topological point defects are grouped together based on the number of points in each cluster. For example, monopoles (singletons), dipoles (clusters with two points), etc. These groups are then analyzed for further classification.

Monopoles are classified as either convex or concave disclinations depending on the charge of the point defect (⁠$+1$ or $−1$⁠, respectively). Dipoles with oppositely charged point defects are classified as dislocations. A dipole with two positively charged point defects is classified as either a target or a spiral. In order to differentiate between targets and spirals, we use persistent homology to detect the presence (or absence) of a local extremum that is situated between the point defects in the case of a target (or spiral).³³ We classify clusters of point defects that are arranged in a near-linear sequence of alternating positively and negatively charged point defects as a grain boundary. Lastly, any cluster that does not fit into one of the above classification schemes, we refer to as an unclassified topological defect, although there are some exceptions to this rule which we describe in more detail in  Appendix B. We do not make any exceptions to the classification algorithm for topological defects located at or near the boundary of the domain.

Note that topological defects are point defects, while the actual pattern defects they correspond to take up a region in the temperature field. For example, visually, a concave disclination is the place where three convection rolls meet, and that the meeting region visually takes up about one roll-width. To account for this, we morphologically dilate the indicator function corresponding to each classified cluster of topological (point) defects by a disk of diameter unit length (see  Appendix A). This gives the classified defects in Fig. 5(b) (with the exception of the wavenumber defects) the shape of a union of balls, each centered at a topological defect.

In addition, we use the local wavenumber $q$ of the pattern as a pattern diagnostic. To compute the spatially varying local wavenumber, we use the approach of Egolf et al.^24,25 The wavenumber is an important parameter that is related to the linear stability of the convection rolls.¹ The linear stability of an infinite field of straight and parallel convection rolls was computed in detail by Busse.⁴⁰ The stability boundaries of straight and parallel convection rolls have proven to be insightful indicators of the pattern dynamics even for finite sized domains that are far from the critical threshold for convection. For example, when the local wavenumber of the pattern crosses a stability boundary this often results in a dynamic response by the pattern to return the wavenumber back within the bounds of linear stability. As a result of this interplay between the pattern dynamics and well known results from linear stability, knowledge of the local wavenumber is an insightful quantity.

For our investigation, where $R=4000$ and $σ=0.84$⁠, there are two particular linear stability boundaries of interest which occur at wavenumbers $q+$ and $q−$⁠. In our notation, $q+$ is the stability boundary of straight and parallel convection rolls at large wavenumber, $q−$ is the stability boundary at small wavenumber, and a local pattern wavenumber of $q$⁠, where $q−<q<q+$ would be in the stable region. A skew-varicose instability occurs for $q+≈2.7$ and a cross-roll instability occurs at $q−≈1.7$ (Ref. 17). We use these values of $q+$ and $q−$ as thresholds in defining two diagnostic functions. At spatial regions where the local pattern wavenumber $q>q+$⁠, we call this a high wavenumber location and where $q<q−$⁠, we call this region a low wavenumber location.

Figures 5(a)–5(b) show the pattern of convection rolls from Fig. 2(a) where we also use color contours to locate different diagnostic quantities of interest. The solid black lines indicate the pattern of the underlying convection rolls by plotting isocontours of $T=0$ at the horizontal midplane. Therefore, the black solid lines represent the center-line axis along the convection rolls where the fluid velocity in the vertical direction is zero. The Lyapunov magnitude function $L(α)$ is represented using the colors of gray and white where we have used a threshold value of $α=0.5$⁠. White indicates regions where the Lyapunov magnitude is above the threshold and gray indicates regions where it is below the threshold.

Figure 5(b) also shows the spatial variation of the eight different pattern diagnostic functions of group I where each particular diagnostic function is shown using a different color. The color bar indicates the particular color that is used to identify the location of each diagnostic function. If the diagnostic function is above its particular threshold, it is shown as its representative color at those regions. If the diagnostic function is below its particular threshold it is not shown in the figure. The different diagnostic functions shown in Fig. 5(b) represent regions of high local pattern wavenumber, regions of low local wavenumber, and the presence of a spiral structure, target structure, grain boundary, a dislocation pair, a convex disclination, and a concave disclination. More details about the specific thresholds used and the methods for computing the classifications is located in  Appendix B.

One important feature of the group I diagnostics is that one only needs an instantaneous image of the flow field pattern to compute these features. In other words, these diagnostics do not explicitly require information from previous times for their computation. We remove this restriction in Sec. III B, where we explore diagnostic functions that do require some knowledge of the flow field history for their computation.

Figure 5(b) conveys the complexity of the pattern and also the rich information contained in group I pattern diagnostic functions. There are several interesting features present which are typical of these patterns and dynamics. First, the diagnostic functions highlight a diverse range of pattern features with some regions containing overlapping diagnostics with other regions only highlighted by a single diagnostic. The connection between the diagnostic regions and the regions where the Lyapunov magnitude function is above its threshold value (represented as white regions) appears quite complicated. For example, several white regions are highlighted by different diagnostic functions, whereas several white regions are not highlighted by any of these diagnostics.

We next study the statistical connection between the different group I diagnostic functions shown in Fig. 5(b) and the Lyapunov magnitude function in detail. We accomplish this by quantifying 3500 images such as Fig. 5(b) using the large set of data generated by our long-time numerical simulation of the chaotic dynamics. We use the precision-recall approach described Sec. III A to generate a statistical description.

Figure 6(a) shows the average precision-recall curves for the group I diagnostics. Each curve represents the results for a different pattern diagnostic function as a function of the value of the threshold $α$ used for the Lyapunov magnitude function. Each different pattern diagnostic function is indicated using a different symbol. Each data point on a particular precision-recall curve is the average precision and average recall for that diagnostic using that value of $α$⁠, where the different values of $α$ are indicated in color.

In addition to the eight diagnostics for specific types of topological defects, we also include two new aggregate classifications in Fig. 6. We group together all of the topological defects in one diagnostic function which we refer to as the topological defect diagnostic function. In addition, in the process of analyzing the patterns, there were topological point defects identified which did not fit the signature of any of the particular defect structures that we have included. Any topological point defect that was not previously classified as a particular type of defect we have gathered together into the pattern diagnostic function that we refer to as unclassified.

The precision-recall curves shown in Fig. 6(a) exhibit some general trends. As the threshold $α$ is reduced, the recall decreases while the precision increases. For a large value of $α$⁠, such as $α≈0.9$⁠, this would indicate that only the spatial regions where the largest peaks of the Lyapunov magnitude function occur are included by setting $L(α)=1$ at these locations. As a result, there are fewer regions highlighted by $L(α)$ and these regions will be located near the largest peaks of the Lyapunov magnitude function. The average precision $P¯α$ of the different diagnostic functions will be low since they typically select numerous regions in the pattern yet now only a few regions have $L(α)=1$⁠. The result is that the pattern diagnostic functions will select regions that do not have $L(α)=1$ which results in a lower precision. In addition, the average recall $R¯α$ will be large since the pattern diagnostic functions will successfully locate more of the regions where $L(α)=1$⁠, on average.

Conversely, it is useful to also describe the results when the threshold $α$ is small, such as $α≈0.1$⁠. In this case, since the threshold for $L(α)$ is small, there are many regions in space where the Lyapunov magnitude function is larger than the threshold. As a result, the precision of the diagnostic functions will increase. This can be understood as follows: of the locations identified where the pattern diagnostic function yields $F=1$⁠, there will now be more locations where $L(α)=1$ which will yield a larger precision. However, for the smaller value of the threshold $α$ the recall will be smaller. This is because at the lower value of $α$ there are now many locations where $L(α)=1$ which will not also be captured by the pattern diagnostic function.

The only diagnostic that does not follow these trends are the regions where the local wavenumber is below the threshold. However, this diagnostic does not perform well for either precision or recall and we have not investigated this diagnostic further.

It is interesting to highlight that the pattern diagnostic function for target structures yields a large amount of precision although its recall is quite small. This indicates that when a target structure is present in the flow field, there is a significant probability that it will also be in a region where the leading Lyapunov magnitude function is above its threshold. However, there are typically very few target structures present in the pattern at any one time which results in a low value for the recall.

Of all of the group I diagnostics represented in Fig. 6(a), there are three which stand out in terms of achieving a significant amount of recall. These are the aggregate of all of the topological defects, the aggregate of all of the unclassified topological defects, and the regions where the local wavenumber is large. It is interesting that the individual topological defects we have quantified do not perform well yet, on average, the aggregate of all of the topological defects performs well. This suggests that there is not a particular topological defect structure, that we have identified in group I, that is a good indicator of where the Lyapunov magnitude function will be above threshold. Rather, this suggests that the different diagnostics capture features of the Lyapunov vector that, when added together, performs better.

In addition, the aggregate of all of the unclassified topological defects performs well in terms of both precision and recall. It is notable that the aggregate of the unclassified topological defects has a higher precision than the aggregate of all of the topological defects. This indicates that removing the specific topological defects we have identified from the rest of the topological defects selects specifically for regions that have a higher correspondence with the high-magnitude regions of the Lyapunov vector. This again suggests that the presence of a topological defect is what is important and not its classification as one of the typical canonical defects.

Lastly, the regions where the wavenumber is large also yields significant precision and recall. The connection between the Lyapunov vector and regions of large wavenumber has been discussed in the literature.^14,19,21,24 However, it is interesting to point out that the local wavenumber is outperformed as a diagnostic, in terms of both precision and recall, by the pattern diagnostic functions for the aggregate of all topological defects and for all unclassified topological defects.

Using the data presented in Fig. 6(a), we have computed the average precision and average recall of each diagnostic over all of the threshold values $α$⁠. Using these average values, we rank each of the diagnostics with respect to one another in terms of precision and recall which is shown in Fig. 6(b). The light shaded region around each symbol indicates the standard deviation about the mean. This allows one to discern between the different diagnostics and to compare their performance with one another on average.

The top three diagnostics in terms of the average recall are the aggregate of all topological defects, the aggregate of all of the unclassified topological defects, and the regions of large wavenumber. In particular, the aggregate of all of the topological defects is ranked first in average recall and is ranked third in average precision. It is also apparent from Fig. 6(b) that while the diagnostic function for target structures ranks first in precision it ranks fifth in terms of recall.

The leading Lyapunov vector quantifies the growth of perturbations in the tangent space, and, therefore, by definition, explicitly describes changes in the system over time. By contrast, the group I diagnostics, which are computed at a given instant on a single flow field image, can capture time-dependence only in an indirect, implicit manner. For example, when the local wavenumber pattern diagnostic is computed on a pattern snapshot, large or small values of $q$ in the pattern suggest that stability boundaries have been crossed and, as a consequence, spatially-localized, instability-driven time-dependence is imminent. In this section, we examine pattern diagnostic functions that include time-dependence explicitly with the goal of determining whether such diagnostics are more effective in locating regions of high magnitude of the leading Lyapunov vector. These diagnostics, which examine short sequences of pattern snapshots from flow image time series, probe topological point defects with a large velocity, the creation and annihilation of topological point defects, and temporal derivatives of the temperature field of the fluid at the horizontal midplane. Additionally, in this section, we introduce two pattern diagnostics for single flow field images that identify, in a novel way, roll pinch-off events and the emergence of target structures.

Figure 7 illustrates the relationship between a chaotic flow field image, the Lyapunov magnitude function, and the group II diagnostic functions at a particular instant of time. The solid black lines again indicate the pattern of convection rolls by plotting isocontours of $T=0$ at the horizontal midplane. The Lyapunov magnitude function is shown using a threshold of $α=0.5$⁠, where white indicates regions above threshold and gray indicates regions below threshold. Each pattern diagnostic is represented using a different color as indicated in the color bar. Figure 7 uses the same flow field image as that of Fig. 5 which allows for a direct comparison between the identified group I and II diagnostic functions for this flow field image.

A roll pinch-off event can be identified as a saddle point in the temperature field at the horizontal midplane. Persistent homology can be used to locate all of the saddle points in the flow field image. However, in a typical convection flow field there are many saddle points located throughout the pattern due to the complex sinusoidal variations of the pattern of convection rolls. The saddle points that are not associated with the roll pinch-off events can be filtered out. The spurious saddle nodes will not be contained in a range of values of the temperature field. We have found that saddle nodes that occur outside of the range $0.2≤T≤0.8$ are spurious, and those that occur within this range are roll pinch-off events.

Emerging targets are characterized by the existence of a local extremum that has a value of the temperature in the range $0.2≤T≤0.8$⁠. We use persistent homology to locate all local extrema in the temperature field and then, as with saddle points, filter their values. For reasons associated with discretization errors, we filter out any saddle points and local extrema that are tied to persistence points with a lifespan less than $δT=0.1$⁠. The lifespan is a measure of the vertical distance a persistence point is from the diagonal in a persistence diagram.

Figure 8 illustrates the use of persistence homology to identify the roll pinch-off events (blue and red circles) and emerging targets (white circles). Both of these defects can be seen as indicating the presence of a local minimum in the amplitude of the pattern. Figure 8 shows the chaotic flow field using a gray scale to indicate the spatial variation of the temperature. The red and blue circles indicate the location of a saddle point, and the white circle indicates a local minimum. In both cases, the values of the selected points occur in the range $0.2≤T≤0.8$⁠.

Figure 9 shows the corresponding regions in the persistence diagrams that are used to select the saddle points and local extrema for this annotation. A detailed description of the use of persistence diagrams for convection flow fields can be found in Ref. 33. For our purposes, we use the fact that each persistence point on a persistence diagram corresponds to a pair of critical points (saddle points or local extrema) in the underlying temperature field. The relationship between the type and dimension of the persistence diagram and the types of paired critical points is summarized in Table I. Thus, by selecting regions on each persistence plane, we are able to select particular critical points in the temperature field. More precisely, the red regions are used to locate saddle points corresponding to pinch-off events occurring in lower temperature ranges (⁠$0.2≤T≤0.5$⁠), the blue regions are used to locate saddle points corresponding to pinch-off events occurring in higher temperature ranges (⁠$0.5≤T≤0.8$⁠), and the gray regions are used to locate local extrema corresponding to emerging target structures. In all cases, persistence points that are a vertical distance of less than 0.04 temperature units from the diagonal are not selected. This ensures that only saddle points and local extrema associated to larger dynamical structures are selected.

We have also quantified the dynamics of the topological defects. In particular, we have used pattern diagnostic functions to locate regions containing topological defects with a large velocity to locate regions where topological defects are created, and to locate regions where topological defects are annihilated. This is done by matching topological point defects with the same charge from one frame in time to the next (see  Appendix C for details). We then use these results to determine which topological defects that have been created, annihilated, or moved farther than a threshold value of $d=0.2$⁠, where the duration between frames is 0.1 time units.

As a measure to identify regions where the pattern changes are rapid in time we have also computed the temporal derivative of the temperature field at the horizontal midplane $∂nT/∂tn$⁠, where $T(x,y,z=1/2,t)$ and $n$ is the order of the derivative. We have explored $n=1,2,3$ using first-order-accurate finite differencing to compute the temporal derivatives of the temperature field. To determine if the temporal derivative is changing rapidly, we have used the threshold $∂nT/∂tn≥0.16$⁠, which was determined by trial and error.

The precision-recall curves for the group II diagnostics are shown Fig. 10(a) and the average rankings of the diagnostics are shown in Fig. 10(b). It is clear that many of these diagnostics perform quite well in comparison with the group I pattern diagnostics shown in Fig. 6.

It is interesting that the roll pinch-off diagnostic function (up triangles) ranks highest in average recall. For example, the percentage of the area where the Lyapunov magnitude function is 0.95 or larger (⁠$α=0.95$⁠) that is also covered by the roll pinch-off diagnostic function is just larger than 60%. These results show that among the group II diagnostics, the roll pinch-off diagnostic is the most effective at locating regions where the Lyapunov magnitude function is above threshold. This is in agreement with the conclusions drawn from previous work on chaotic convection suggesting the significance of the these defect events in relation to the spatiotemporal dynamics of the leading Lyapunov vector.^12,14,19,26 However, the roll pinch-off diagnostic ranks last in precision. This indicates that this diagnostic also selects most regions, on average, that are not where the Lyapunov magnitude function is large.

The emerging target diagnostic function (stars) ranks first in precision, on average. In fact, the emerging target has the highest precision of all of the group I and II diagnostics. However, the average recall of the emerging target diagnostic is quite low and ranks seventh of the group II diagnostics. This indicates that if an emerging target structure is present in the flow field, there is a significant probability that the region occupied by the emerging target corresponds to a region where the Lyapunov magnitude function is above threshold. However, since the recall is low this indicates that the emerging target will not be capturing a significant portion of the area where the Lyapunov magnitude function is above threshold. In short, if an emerging target is present it is a good indicator of where the Lyapunov vector is large, on average.

We next discuss the diagnostic functions that are based on the dynamics of topological defects. The results from our analysis of the group I diagnostics suggest that the aggregate of all topological defects (shown in Fig. 6) is a good indicator in terms of both precision and recall. Our intention here is to explore the possibility of increasing the precision and recall performance of a diagnostic function based on the presence of a topological defect by including some additional measure of its variation with space and time.

Figure 10 yields that topological defects with a large velocity (pentagons) rank third in average recall and fourth in average precision. These results suggest that topological defects with a large velocity are relatively good indicators of where the Lyapunov magnitude function will be above threshold.

The precision and recall results for the annihilation of topological defects are plotted using left triangles in Fig. 10. The annihilation of a topological defect ranks second in precision and sixth in recall. The high value of precision indicates that the annihilation of a topological defect, if present, is a relatively good indicator of where the Lyapunov magnitude function is above threshold. However, since the recall is low the annihilation of a topological defect will not capture a large portion of the domain where the Lyapunov magnitude function is above threshold. Overall, the performance of the annihilation of a topological defect diagnostic function is quite similar to that of the emerging target.

The diagnostic function for the creation of a topological defect is shown in Fig. 10 using right triangles. The precision-recall curves of the creation and annihilation of topological defects shown in Fig. 10(a) are quite similar. However, there are some interesting differences that are indicated by the average diagnostic results shown in Fig. 10(b). In particular, the creation of a topological defect performs worse than the annihilation of a topological defect in terms of both recall and precision. In fact, the creation of a topological defect ranks last in recall and sixth in precision among all of the group II diagnostics. We currently do not have an understanding of why the annihilation of topological defect outperforms the creation of a topological defect.

We next discuss the results using temporal derivatives of the temperature field at the horizontal midplane as a pattern diagnostic function. We have explored the first order, second order, and third order time derivatives of the temperature field which are represented on Fig. 10 as circles, squares, and diamonds, respectively. The average recall of the first order derivative ranks second; however, its precision is poor and ranks seventh. As the order of the derivative increases, the rank of the average recall decreases while the rank of the average precision increases. This indicates that increasing the order of the derivative results in the diagnostic function covering less area where the Lyapunov magnitude function is above threshold while having a higher probability that the area it covers does coincide with regions where the Lyapunov magnitude function is above threshold.

In this section, we discuss some overall features of the diagnostic features we have explored. Figure 11 shows box and whisker plots for the distributions of nonzero probabilities over the entire time series for the pattern diagnostic regions, and Fig. 12 shows box and whisker plots for the distributions of the thresholded leading Lyapunov vector magnitude functions. The red lines in the center of each box denote the median value of the distribution, and the upper and lower extremes of each box give the first and third quartiles, respectively. The whiskers (vertical dashed lines) extend to the maximum and minimum values of the distribution, with the exception of outliers which are indicated by a “+” symbol.

Figure 11 shows that the probability distributions for the different pattern diagnostics are, in general, bounded away from one, typically taking up less than ten percent of the domain for the majority of diagnostics. This finding is supported by Figs. 5–7 which illustrate the regions selected by the different pattern diagnostic functions for a single flow field image. Thus, we can conclude that the recall [vertical axes in Figs. 6(a) and 10(a)] is not artificially inflated by our choice of diagnostics. Additionally, the low probability values for Fig. 12, specifically for larger thresholds, paired with the fairly low probability values in Fig. 11, indicates that the high recall values in Figs. 6(a) and 10(a) for higher magnitude thresholds are significant. This indicates that the center of the highest peaks are contained within the collection of diagnostics, for if not, the recall would decrease as the Lyapunov vector magnitude threshold increases.

The diagnostics of group I and II were selected using experience and intuition. There are connections between many of them and the diagnostics are not meant to be mutually disjoint with respect to one another. The relationships between the diagnostic functions are illustrated in Fig. 13 where we plot the conditional probability matrix $Pi,j$⁠. The conditional probability matrix is defined as $Pi,j=P(Fi|Fj),$ where $Fi$ and $Fj$ are pattern diagnostic functions and the indices $i$ and $j$ cycle through the 18 different diagnostic functions we have quantified. The magnitude of $Pi,j$ is plotted using the color scale shown. Red indicates a large value of $Pi,j$ which yields that pattern diagnostics $Fi$ and $Fj$ are related to one another in the sense that their conditional probability is large. Similarly, blue indicates a small value for the conditional probability which yields that the two diagnostics $Fi$ and $Fj$ are not related to one another in a statistical sense.

Therefore, the diagonal entries where $i=j$ yield $Pi,i=1$ since this compares a diagnostic function with itself which are represented as red. The off-diagonal entries indicate the average pairwise conditional probabilities of the different diagnostic functions. The first row of $Pi,j$ is for the diagnostic function that is the aggregate of all of the topological defects. The diagnostic functions for the different topological defects from group I are given in rows 2–8 as shown in Fig. 13. These topological defects are a subset of the diagnostic for the aggregate of all of the topological defects and as a result $P1,2:8=1$⁠. It is also clear that the diagnostic for the first order temporal derivative statistically contains the pattern diagnostics for emerging targets and for topological defects with a large velocity. Also, the different orders of the derivatives are statistically related with one another as expected. In addition, since the specific topological defects and the unclassified topological defects (⁠$F2,…,F8$⁠) partition the aggregate diagnostic for all of the topological defects, they are mutually disjoint. This is indicated by the submatrix $P2:8,2:8$ of Fig. 13 being approximately a $7×7$ identity matrix. The only notable exception is $P2,6$⁠, which gives the conditional probability of a concave disclination when restricting to target classifications. The reason for this anomaly is described in  Appendix B.

Our results suggest a complex relationship between the flow field structures and the spatiotemporal features of the leading Lyapunov vector of chaotic Rayleigh–Bénard convection. We have explored two different groups of pattern diagnostic functions. Group I diagnostics contain many of the conventionally used approaches for describing patterns such as the local wavenumber and different classifications of topological defects. Group I diagnostics can be determined from a single flow field image and include the time dependence of the dynamics only in an implicit sense. On the other hand, group II diagnostics have been chosen to explicitly include the time dependence of the dynamics and include, for example, the appearance or annihilation of a defect structure.

Generally speaking, we have found that group II diagnostics have better precision, but lower recall, than group I diagnostics. The increased precision is reasonable since group II diagnostics specifically include aspects of the time evolution of the patterns which were anticipated to be important. This insight regarding the importance of the time history of the flow is useful for researchers interested in building a description of high-dimensional dynamics using data driven approaches such as machine learning and different methods of modal decompositions.

We have found that powerful ideas from homology can be very useful in the quantification of the patterns of chaotic convection. Specifically, we have used ideas from homology to identify spirals, emerging target structures, roll pinch-off events, and to compute the velocity of topological point defects. Finally, we have used the ideas of precision and recall to build a statistical description of the ability of these diagnostics to describe the leading Lyapunov vector. The combination of high-resolution numerical data for three-dimensional chaotic convection, ideas from persistent homology to help generate pattern diagnostic functions that include the time history, and the statistical description gained from precision-recall curves has led to new insights. This approach is quite general and could be used on a wide range of problems where large amounts of highly resolved data is available.

One aspect that stands out from this study is the high precision of the diagnostic for the emerging target structure. This indicates that if an emerging target is present in the flow field it has the highest probability of all of the pattern diagnostics we have studied here to be in a region where the magnitude of the leading Lyapunov vector is above threshold. This indicates that the spatial region occupied by an emerging target will be highly sensitive to small perturbations to the flow field. This suggests that the active control of emerging target structures in some fashion could perhaps be used to stabilize or destabilize a flow field.

Our results also shed light on the influence of roll pinch-off events on the dynamics. The roll pinch-off has the highest recall yet lowest precision of all of the group II diagnostics. Also, the pattern diagnostic for the aggregate of all topological point defects outperforms all of the pattern diagnostics in terms of recall. It is interesting to highlight that we did not find a single type of topological point defect that significantly outperformed the other classified topological defects in terms of precision and recall. This suggests that if one is interested in identifying regions of the chaotic pattern where the leading Lyapunov vector magnitude is large, a useful first pass of the pattern is to identify the regions containing topological point defects.

For high-dimensional chaotic dynamics, there are many Lyapunov vectors that have positive Lyapunov exponents. For the dynamics we have studied here, there are 24 Lyapunov vectors that have positive Lyapunov exponents. By using only the leading Lyapunov vector, we are not using all of the rich information contained by the spectrum of Lyapunov vectors. It is possible that the connection between the patterns and the regions of rapid divergence in the tangent space may become clearer with the inclusion of more Lyapunov vectors.

In our study, we have used the Gram–Schmidt orthonormalized Lyapunov vectors and therefore only the spatiotemporal dynamics of the leading Lyapunov vector is physically relevant. However, it is now possible to compute the spectrum of covariant Lyapunov vectors^26,41–44 which are not orthonormalized and are expected to contain useful information in their spatiotemporal dynamics. As a result, one could imagine quantifying the correlations between the spectrum of covariant Lyapunov vectors and the pattern dynamics as we have done here for the leading Lyapunov vector. However, the computation of the covariant Lyapunov vectors in experimentally accessible convection domains is much more expensive^26,44 and is a topic of future interest.

The authors thank the following for many fruitful discussions and insights: Miro Kramar, Jeff Tithof, Balachandra Suri, Brett Tregoning, Logan Kageorge, and Saikat Mukherjee. The authors also thank Shaun Harker for his contributions to the code used to perform this study, in particular, the persistent homology computations and the implementation of the algorithm of Bazen and Gerez.³⁹ 

The research was supported by the DARPA MoDyL Program, DARPA Contract Nos. HR0011-17-1-0004 and HR0011-16-2-0033 and by National Science Foundation (NSF) Grant No. DMS-1622299. The computations were conducted using the resources of the Advanced Research Computing center at Virginia Tech.

Morphological dilation is a preprocessing technique used in image analysis.³⁷ For any binary image (i.e., with values in ${0,1}$⁠), it serves to increase the area corresponding to the nonzero pixels. The amount of area increased and the shape of the dilation is determined by what is called a structure element. In this paper, we always use a structure element in the shape of a disk since it is rotation-invariant, but it is common in the image analysis community to consider alternative shapes for structure elements, e.g., a square. In what follows, we give a mathematical description of the morphological dilation of a real-valued function.

Let $E⊆Rn$ be a subset of Euclidean space. For $x∈Rn$⁠, define the set

or the set $E$ translated by $x∈Rn$⁠. Let $1E$ be the indicator function on $Rn$ corresponding to the set $E$⁠, i.e., it takes the value $1$ on $E$ and $0$ otherwise. Given a real-valued function $G:D→R≥0$ on a subset $D⊆Rn$ of Euclidean space, we define the morphological dilation of $G$ by structure element $E⊆Rn$ as

In other words, given a subset $E⊆Rn$ and $x∈D$⁠, the value $(G⊕1E)(x)$ is defined to be the maximum value $G$ takes over the set $(x⊕E)∩D$⁠.

The use of the term dilation becomes apparent when one takes $E$ to be the unit ball centered at the origin in $Rn$⁠. To visualize, let $x∈Rn$ and set a function $F:=1{x}$⁠. Then $F⊕1E=1x⊕E$⁠, the indicator function on the unit ball centered at the point $x$⁠. Thus, the single nonzero point in $F$ takes up a region the size of $E$⁠.

Lastly, we note that given the above definition of a morphological dilation on real-valued functions, one can first threshold an image to binary values and then take its morphological dilation, or first take an image’s morphological dilation and then threshold, and the resulting binary-valued functions are the same. That is, morphological dilation commutes with thresholding. In the case of dilating diagnostic indicator functions, this does not apply, since the input functions are already binary-valued. However, we morphologically dilate the Lyapunov vector magnitude function by a disk of unit diameter [see Fig. 2(c)], and then threshold this dilated function in our analysis, which would yield the same results as first thresholding and then dilating.

As described in Sec. III B, topological defects are classified using the following approach. The topological point defects are clustered using single-linkage hierarchical clustering with a radius of $r=1$⁠, and then each cluster is classified as follows.

Singleton clusters, or monopoles, are classified according to their topological charge. Defects with positive charge (⁠$+1$⁠) are classified as concave disclinations and defects with negative charge (⁠$−1$⁠) are classified as convex disclinations. Any cluster with two topological point defects of opposite charge (⁠$+1,−1$⁠) is classified as a dislocation pair.

A cluster with three or more points that is arranged in a near-linear sequence of topological point defects of opposite charges is classified as a grain boundary. To detect this situation, we use basic ideas from principal components analysis (PCA). The set of topological point defects is treated as a point cloud in $R2$ and then PCA is run on this point cloud. The eigenvalue corresponding to the second principal component is then checked to see if its magnitude falls in the range $25r,32r$⁠, where $r=1$⁠, indicating that the points are arranged in a near-linear (but not completely linear) pattern. Finally, the points are orthogonally projected to the first principal component and checked for alternating charges in the resulting ordering. If both of these tests are true, then the cluster is classified as a grain boundary.

Dipole clusters with two positively charged topological defects are necessarily either targets or spirals. To distinguish the two, we use persistent homology to detect the existence of a local extremum that is positioned between the two topological defects through the method illustrated in Fig. 8.

Targets and spirals are often bounded by two or more concave disclinations, and if the target or spiral is small enough, the clustering radius may group together the target (or spiral) and its surrounding concave disclinations. For example, the target pattern identified in the 6 o’clock position in Fig. 5 is bounded by three concave disclinations. Any cluster of topological defects containing only two positively charged topological defects, regardless of how many defects are in the cluster, is checked for this scenario as follows. Let $r$ be the radius of the smallest circle circumscribing the two positively charged topological defects. If the circle of radius $1.25r$ centered at the midpoint between the two positively charged topological defects does not contain a negatively charged topological defect, then all negatively charged topological defects are marked as concave disclinations, and the two positively charged topological defects are classified using the method described above, as if they formed a dipole.

We use a simple algorithm for matching topological defects of like charge from frame to frame in the simulation. Our algorithm rests on the assumption that the vast majority of topological defects are either stationary or move only minimally as the pattern changes, and that the flow is sampled densely enough in time to resolve the majority of any ambiguities. We have found that this approach works well for the chaotic convection data we have explored.

Let $TDt(c)$ be the collection of topological defects of charge $c$ at time $t$⁠. For each $x∈TDt(c)$⁠, we compute the nearest defect $y∈TDt+Δt(c)$ and match $x$ to $y$⁠. We then reverse the roles of time points $t$ and $t+Δt$ and compute the reverse matchings. Any defects that were matched to each other in both rounds and have distance less than 0.7 from each other are consider matched to each other. All other points are considered unmatched. The distance filter is included to reduce the chance that a newly annihilated topological defect and a newly created topological defect are not erroneously matched to each other.