Extraction of Distinguished Hyperbolic Trajectories for 2D Time‐Dependent Vector Field Topology

This paper does two main contributions to 2D time‐dependent vector field topology. First, we present a technique for robust, accurate, and efficient extraction of distinguished hyperbolic trajectories (DHT), the generative structures of 2D time‐dependent vector field topology. It is based on refinement of initial candidate curves. In contrast to previous approaches, it is robust because the refinement converges for reasonably close initial candidates, it is accurate due to its adaptive scheme, and it is efficient due to its high convergence speed. Second, we provide a detailed evaluation and discussion of previous approaches for the extraction of DHTs and time‐dependent vector field topology in general. We demonstrate the utility of our approach using analytical flows, as well as data from computational fluid dynamics.


Introduction
Vector field topology (VFT) is a widely used concept for the visualization of steady vector fields [HH89]. It is defined by isolated stationary points and their separatrices, which are (instantaneous) streamlines that converge to saddle-type critical points in forward or reverse direction and partition the domain into regions of similar flow behavior. While this concept can be applied to time-dependent vector fields by freezing an instance in time, it is unable to capture the true Lagrangian motion of particles over time, even if an appropriate frame of reference is chosen. Instead, the notion of Lagrangian coherent structures (LCS) is used for the definition of vector field topology in time-dependent flow. In 2D, LCS represent material lines [SLM05], which act as transport barriers and are locally the most attracting or repelling. They are included in the set of ridges of the finite-time Lyapunov exponent (FTLE) The main challenge posed by time-dependent vector field topology is the extraction of the DHTs. In order to obtain intersections of ridges in forward and reverse FTLE fields, a dense grid of particles would need to be integrated in forward and backward time. To ensure robust extraction of the ridges, this would need to be done at very high resolution (or in an adaptive manner). However, performing this costly computation at a single instance of time would not be enough: numerical integration started at such an intersection point exponentially accumulates errors along the repelling or attracting manifolds in forward or backward time, and additionally DHTs can appear and disappear over time. So far, two local extraction methods have been proposed to solve this problem. Machado et al. [MBES16] locally track critical points in a Galilean-invariant reference frame, and locally refine the solution toward the closest pathline. Secondly, Rojo and Günther [RG20] solve a least-squares optimization problem to obtain more accurate frames of reference, in which critical points follow pathlines more closely. While the local refinement scheme of Machado et al. does not need to converge to a DHT or at all, paths of critical points in an optimal frame of reference, as computed by Rojo and Günther, do not need to be DHTs or even pathlines, as we are going to show in Section 6.8.
Both approaches have in common, that they are mostly local in nature. In this work, we bridge the gap between local and global integration-based techniques. First, we locally extract initial candidates by tracking critical points in appropriately chosen reference frames. Second, we globally compute the dynamics in the localized flow along each initial candidate, using numerical integration. This is similar to computing the localized FTLE [KPH * 09] along a single pathline, but with the path given by the initial candidate instead of pathline integration. Finally, closely following the works of Ide et al. [ISW02] as well as Branicki and Wiggins [BW09], we compute a time-dependent coordinate transformation, which decouples repelling and attracting directions into one-dimensional systems, which allows for separate integration along the repelling direction in backward time and along the attracting direction in forward time, thus refining the initial candidate toward the DHT. Having obtained a set of DHTs, we compute streak manifolds in the 3D space-time domain, seeded along the DHTs at an offset determined by the coordinate transformation. The streaklines obtained in this way represent LCS, and thus the topology of time-dependent flow.
Our contributions include: • We combine existing local DHT extraction techniques with global numerical approaches [ISW02,BW09] and refinement. • We obtain consistent offsets for seeding streak manifolds.
• We evaluate different approaches to obtain DHT candidates.
• We compare our technique with previous approaches, and validate our results using the FTLE. • We show that no purely local technique is able to accurately extract LCS-based features in general time-dependent flow.

Related Work
Haller [Hal01] and Shadden et al. [SLM05] proposed to compute the largest finite-time Lyapunov exponent on a grid and obtain Lagrangian coherent structures as ridges in these fields. This concept has proven successful for the visualization of unsteady vector fields [GLT * 09, SP09]. Instead of using a computational grid to obtain the necessary flow map gradients, Kasten et al. [KPH * 09] used the localized flow along trajectories. Generally, computing gradients of the flow map, which maps an initial position of a massless particle to its position after a given integration time, requires very fine sampling especially near LCS. To counter this, adaptive methods have been proposed [GGTH07,SP07,SRP11]. Sadlo and Weiskopf [SW10] adaptively computed intersection points of ridges of the forward and backward FTLE fields. The authors have shown, that such points lie on DHTs, and that generalized streak lines [WTS * 07], seeded along DHTs, represent LCS. In this way, steady vector field topology is extended to unsteady vector fields by replacing streamlines with generalized streaklines. This method has been extended to 3D by Üffinger et al. [USE13]. Both methods require numerical integration of hyperbolic trajectories, which leads to exponential error growth due to their attracting and repelling subspaces. A computationally expensive, but accurate alternative is to compute the FTLE in dense time slices instead [BSDW12,SBDW13]. Besides computational cost, using the FTLE for time-dependent flow analysis also has the disadvantage, that a poor choice of its integration parameter can lead to inconsistencies with LCS [SLM05,BW10]. Also based on the flow map, Bujack et al. [BDZG19] defined time-dependent equivalents of saddles, sinks, and sources, which are also included in the DHT definition by Ide et al. [ISW02]. In this work, we extract LCS explicitly in geometric representation, avoiding computation of discretized flow maps and error-prone extraction of LCS by ridges therefrom.
Steady vector field topology has been introduced by Helman and Hesselink [HH89] for the visualization of steady vector fields. Since it is defined by streamlines and not pathlines, it is not suitable for visualization of general time-dependent vector fields. Theisel et al. [TWHS04] considered a pathline topology, which splits the space-time domain into regions of instantaneous hyperbolicity. The extraction of vortex core lines in time-dependent vector fields is closely related. The criterion by Sujudi and Haimes [SH95] is widely used for their extraction in steady 3D vector fields, which can be formulated using the parallel vectors operator [PR99], and requires other streamlines to swirl around the vortex core line. Weinkauf et al. [WSTH07] extract pathline cores by applying the parallel vectors operator to the extended space-time flow. Due to the special structure of this flow, the problem is reduced to finding critical points in 2D or extracting 3D vortex core lines from a modified flow, where the feature flow field [TS03] is subtracted. Günther et al. [GST16] noted, that subtracting the feature flow field can be regarded as observing the flow in a Galilean-invariant frame of reference. Using linear optimization, Günther et al. [GGT17] obtained locally objective frames of reference. This has been extended to affine invariance [GT20]. Hadwiger et al. [HMTR18] formulated the problem of finding optimal reference frames as a global optimization problem. Rojo and Günther [RG20] extended the local reference frame optimization framework to displacement transformations, which includes the preceding methods as special cases by allowing for arbitrary spatial deformations. They define a timedependent vector field topology in terms of the steady vector field topology in the steady frame of reference. Closely related to this strictly local approach is that of Machado et al. [MBES16]. These authors extract bifurcation lines using the parallel vectors operator from the space-time domain, which is equivalent to tracking critical points in the Galilean-invariant frame of reference defined by the feature flow field. In a second step, the authors use a local refinement technique [MSE13] to correct the initial line toward the nearest pathline. The time-dependent vector field topology is finally obtained as streak lines, as proposed by Sadlo and Weiskopf [SW10]. Since Rojo and Günther compute separatrices as streamlines in a steady vector field instead, they are not limited by the time domain of the dataset. However, as we are going to show, this approach does not yield LCS in general.
Haller [Hal00] defined hyperbolic trajectories (here denoted DHT), as those pathlines that stay in a hyperbolic region for locally the longest time. His notion of hyperbolicity is based on the determinant of the instantaneous Jacobian being negative, and is employed by Sadlo and Weiskopf [SW10], Üffinger et al. [USE13], and Machado et al. [MBES16]. A similar notion based on exponential dichotomy in dynamical systems [Cop78] is that of a distinguished hyperbolic trajectory. It was introduced by Ide et al. [ISW02], who also formulated a numerical algorithm based on Fourier transform, which computes a DHT from an initially obtained path of critical points in the lab frame of a dataset. Since the Fourier transform imposes time-periodicity on the dataset, an approach geared toward more general datasets was developed [MSW04]. For solving the integral equation associated with a DHT, both approaches use a linear approximation of the flow. For the extension to 3D, Branicki and Wiggins [BW09] use a fixed-point iteration instead, which was originally used for proving existence and uniqueness theorems of DHTs [JSW03]. In our work, we closely follow these approaches in order to correct an initial candidate line toward a DHT, while for identifying initial candidate lines, we rely on Haller's definitions. In this way, we obtain a 2D time-dependent vector field topology consistent with LCS as obtained by FTLE ridges.

Fundamentals
We consider a 2D time-dependent vector field u(x,t) ∈ R 2 , with x ∈ Ω ⊂ R 2 , defined over a finite time interval t ∈ [t 0 ,t N ]. The Lagrangian motion of massless particles in such a field is described by pathlines, which are tangent curves given an initial seed point x(t 0 ) = x 0 .

FTLE and Lagrangian Coherent Structures
The flow of a time-dependent vector field is given by its flow map φ φ φ T t0 (x), which maps a seed point x at time t 0 to its final position after advection time T . Exponential separation, the main organizing structure of the flow, can be obtained from the FTLE field, ∇u ∇u ∇u Figure 2: The localized flow along a pathline x(t) (green) describes the evolution of infinitesimal perturbations δ δ δe 1 (t 0 ) = e 1 , δ δ δe 2 (t 0 ) = e 2 (left ellipse) around it over a time interval [t 0 ,t N ]. In first-order approximation, the evolution is obtained from the Jacobian ∇u along the pathline. The integrated perturbations δ δ δe 1 (t), δ δ δe 2 (t) form the columns of the fundamental solution matrix Xt 0 (t).
where σmax(·) denotes the largest singular value. For an ndimensional flow, the FTLE is a spectrum of n exponents, but the largest is typically used as a synonym. Ridges in ς T t0 (x) for negative advection time represent attracting LCS, while they represent repelling LCS for positive advection time [Hal01]. Generally, particles separate exponentially in forward time, when they are close to repelling LCS, and they separate exponentially in backward time, when they are close to attracting LCS.

Streak-Based Topology
Since LCS are material lines [SLM05], the intersection of repelling and attracting LCS represents pathlines. These include special pathlines, from which streaklines seeded at an offset represent the corresponding LCS. Interpreting these pathlines as saddle-type degenerate streaklines motivates the definition of streak-based topology [SW10]. The generating pathlines of this topology are the DHTs, which are those pathlines, that reside in hyperbolic regions (where det ∇u < 0), for locally the longest time [Hal00].

Separation and Attachment at No-Slip Boundaries
In time-dependent flow with no-slip boundaries, the LCS defined by DHTs need to be extended by streak manifolds seeded along space-time separation and attachment lines at the no-slip boundaries [MBES16]. Separation and attachment lines can be obtained [Ken98] using the parallel vectors operator [PR99].

Localized Finite-Time Lyapunov Exponents
While the FTLE field conceptually describes separation of neighboring particles, we now consider the (localized) finite-time Lyapunov exponents of single pathlines. Local separation and attraction in an infinitesimal neighborhood ( Figure 2) of a pathline x(t), with x(t 0 ) = x 0 , over a time interval [t 0 ,t N ], is described by the localized flow This ODE describes the evolution of an infinitesimal neighborhood δ δ δx(t) along x(t), which is captured entirely by its solution with initial condition Xt 0 (t 0 ) = I. We call this solution the fundamental solution matrix Xt 0 (t). In previous work, it has been used with orthogonal matrices Bt 0 (t), Rt 0 (t), and a diagonal ma- , the localized finitetime Lyapunov exponent spectrum is obtained as Note, that the finite-time Lyapunov exponent is often refered to as the largest one, i.e., ς t−t0 t0 = λ 1 t0 (t), which measures the exponential separation from time t 0 to t. However, since integration can be reversed, we have Xt (t 0 ) = Xt 0 (t) −1 . Thus, we also obtain exponential attraction from the smallest Lyapunov exponent λ 2 t0 (t), i.e., exponential separation in backward direction from time t to t 0 . This fact can also be exploited for computing both FTLE fields using a grid of pathlines from a single computation [HS11].

Distinguished Hyperbolic Trajectories
If a pathline x(t) is structurally stable, the localized flow (Equation 3) can be transformed into the linear system d dt which is defined by the diagonal, time-independent matrix This is achieved by the time-dependent coordinate transformation A structurally stable pathline is characterized by none of the diagonal entries of Dt 0 (t N ) being zero. It is attracting if all entries of Dt 0 (t N ) are negative, and repelling if all entries are positive.
A hyperbolic trajectory has both positive and negative entries in Dt 0 (t N ) in this notation. Since with this notion, almost all trajectories are hyperbolic, Ide et al. [ISW02] define a DHT as a hyperbolic trajectory, which remains in a bounded neighborhood B for all time, while all other trajectories starting in B leave B at exponential rate in forward and backward time, and, in addition, the trajectory it is not an intersection of attracting and repelling manifolds of other DHTs. The latter condition is important (Section 6.3), since such non-distinguished trajectories manifest themselves as falsepositive ridge intersections of the forward and backward FTLE fields, but they are not topological generators of the time-dependent flow. This definition of a DHT is more general than the notion of Haller We now use the preceding discussion to obtain a refinement scheme, which corrects an initial candidate linex(t) toward a nearby DHT. Using the same definitions as above, we change into the coordinate frame w(t) = T(t)(x(t) −x(t)). The flow (Equation 1) then takes the form where the nonlinear part h(w(t),t) is given by [ISW02] h In Equation 9, attracting and repelling behavior is decoupled into one-dimensional systems. Thus, we are able to integrate along attracting directions in forward time, and along repelling directions in backward time. An approximate DHT for a given finite-time in- where w i and h i denote the ith component of w and h, and d i the diagonal entries of Dt 0 (t N ). Since integration starts and stops at the time boundaries of the initial candidate, a solution w(t) takes the special values The decoupled integration in Equation 11 reverses the repelling behavior of the DHT and thus corrects toward it in both forward and backward time (see Figure 3). This also explains the initial values in Equation 12: the DHT does not exist across these time boundaries, and integration beyond them, even if the data would permit it, would not undergo the same hyperbolic behavior and thus not yield a correction. Finally, we obtain a refinement of the initial candidatex(t) toward an approximate DHT, An approximate DHTx DHT obtained in this way, converges to the actual DHT for infinite integration time, in cases, where a DHT exists for all time [ISW02] and the initial candidate is within a bounded region, depending on its Lyapunov exponents [JSW03].

Method
Based on the previous discussion, we extract the time-dependent vector field topology of a vector field u(x,t) ∈ R 2 given on a time interval [t 0 ,t N ] using the following steps: 1. Locally extract a set of initial candidate linesx(t). 2. Refine the initial candidates toward DHTs. 3. Seed streaklines along the obtained DHTs.

Obtaining Initial Candidates for DHTs
Candidatesx(t) for DHTs are obtained by tracking saddle-type critical points in a suitable frame of reference using the parallel vectors operator in the space-time domain. Given a decomposition into observer motion f(x,t) and an observed vector field w(x,t), we obtain initial polylines as paths of critical points in w, which we filter by instantaneous hyperbolicity det ∇u < −τ h with a large positive threshold τ h . These short but robust solution lines are then extended by seeding at their endpoints streamlines in the observer motion field f in forward and reverse direction, until integration leaves the hyperbolic region (det ∇u ≥ 0), or a domain boundary is reached (e.g., Figure 12o). Notice, that f is not a feature flow field of w in general, since it neglects possible motion of observed critical points in w relative to the observer motion f. In Section 6.8, we evaluate using candidates in the lab frame [ISW02], the Galilean-invariant frame of reference defined by the feature flow field [MBES16], and an optimal frame of reference [RG20].

Refinement toward DHTs
The initial linex(t) is given as a polyline (x 0 ,t 0 ), . . . , (x N ,t N ) in space-time. The Jacobian matrix is computed at each of these space-time locations, . From this, the singular value decomposition (Equation 4) of the fundamental solution matrix X(t,t 0 ) is computed by numerical integration of the initial value problem where we linearly interpolate J i between the discrete time steps t 0 , . . . ,t N to obtain J(t). To avoid numerical issues, that occur for strong hyperbolicity or long integration times, we use the continuous SVD method (see Section 1 in the supplemental material) to obtain the matrices Bt On each of these non-overlapping intervals, as well as on intervals shifted by T /2, we compute a refined DHT. The overlapping results are averaged using cosine weights, thus ensuring a uniform precision along the DHT. The refinement toward the DHT is obtained as the solution of the implicit integral Equations 11, which we solve using a fixed-point iteration. Starting with initial guess w (0) (t) = 0, we compute h(w (0) (t),t) according to Equation 10, and evaluate the integrals using trapezoidal rule, yielding an approximate solution w (1) (t) of Equation 11. This process is iterated until w ( j+1) − w ( j) drops below a predefined threshold τ f or a maximum number of iterations is reached. A refined DHT is then obtained asx DHT (t) from Equation 13. Since the localization of the initial pathx(t) is fixed, we usex DHT (t) as new initial candidate, and repeat the entire process. The iteration terminates when x ( j+1) DHT reaches a predefined threshold τ i .

Seeding Streak Manifolds
The attracting and repelling manifolds of the DHTs, i.e., the LCS, are extracted by computing streamsurfaces in space-time (streaklines in space), offset in the perturbation directions that belong to the respective Lyapunov exponents. These directions yield linear approximations of the corresponding LCS and are obtained as the columns of T(t) −1 . Previous work [USE13,MBES16] used the real eigenvectors of the Jacobian as an estimate, which, due to their instantaneous nature, are in general not well aligned with the manifolds (Figure 4). In general, DHTs only exist over finite time intervals, which are typically shorter than the time domain. As shown by Üffinger et al. [USE13], streak manifold integration has to continue across the entire space-time domain, however.

Hyperbolicity Strength
The localized view on the distinguished hyperbolic trajectories in Section 3 gives a quantitative measure for exponential separation along the trajectory. It is given by the diagonal matrix Dt 0 (t N ), i.e., nearby trajectories separate over time as exp(tDt 0 (t N )). Therefore, we may use ϑ = max(t N − t 0 )Dt 0 (t N ) to measure the hyperbolicity of a trajectory, which not only depends on the finite-time Lyapunov exponents, but also on the duration in time.
Using the method described in Section 4.1, the same DHT may be extracted multiple times. In order to obtain a set of unique DHTs, we filter DHTs that are too close to each other, according to minimum point-wise distance, and choose those among the corresponding ones, that have the largest hyperbolicity strength ϑ. By choosing the distance threshold smaller than the cell size of the grid defining u(x,t), the interpolation makes it unlikely, that DHTs are filtered as false negatives in this process. A thorough threshold based on the interpolation scheme is subject of future work.

Implementation
We use an embedded Runge-Kutta 4/5 scheme with adaptive step size and dense output to numerically obtain the fundamental solutions of the localized systems. For our results, we use a relative tolerance of 10 −3 and absolute tolerance of 10 −6 . The maximum step size is chosen as the mean step between discrete time steps t i+1 −t i of the candidate line, and the initial step size as a 10th of this. We use the dense output to obtain solutions at the predefined discrete time steps t i , which may not be reached exactly due to the adaptive step size. The iterations in the DHT refinement (Section 4.2) are stopped if two consecutive steps differ less than τ f = τ i = 10 −10 , or after at most 100 steps. A C++ prototype of our DHT refinement is provided in the supplemental material.

Results
In the following, we first discuss time-dependent vector field topology at simple analytical examples, to build intuition about related issues, and compare the efficiency and accuracy of our DHT refinement method to the local refinement approach by Machado et al. [MSE13]. Finally, we evaluate the different approaches to obtain initial candidate lines (Section 4.1) at two numerical flow simulation datasets, and compare the approaches by Rojo and Günther [RG20] as well as Machado et al. [MBES16] to our DHT refinement approach on these datasets. As ground truth, we compare against ridges in the FTLE fields in all our examples. We implemented the local bifurcation line refinement as described by the authors [MSE13], and used the provided prototype for computing the optimal reference frames [RG20], where we chose second-order optimization with a neighborhood of 41 2 nodes.

Skewing Oscillating Gyre-Saddle
Based on the Skewing Gyre-Saddle and Oscillating Gyre-Saddle examples proposed by Sadlo and Weiskopf [SW10], we construct a model for a hyperbolic region, which oscillates sinusoidally between (0.2, 0.8) ⊤ and (−0.2, −0.8) ⊤ with a period of 4 s, with the saddle directions skewing at a period of 1 s. Figures 3 and 4 show this dataset over the time interval [0, 4]. While this dataset is timeperiodic and thus would allow to compute the FTLE beyond the time boundaries, we choose not to do so. Since we only extract the DHT and its attracting and repelling manifolds over this time interval, the resulting topology will only be consistent with ridges in the FTLE fields constrained to the same time interval. For example, the forward FTLE at t=0 s (see Figure 3, front slice) can thus only constrain the DHT to a line at this instant of time, while the forward and backward FTLE fields at the center of the time interval (t=2 s) constrains it to the point, where forward and backward FTLE ridges intersect. This is perfectly consistent with the computation of streak manifolds, as shown in Figure 3.

The Beads Problem
The Beads problem is an example for an attractor, around which pathlines exhibit swirling motion. In the following, we use the analytic model and ground truth previously used by Weinkauf and Theisel [WT10]. While traditional feature extraction methods fail to find this attractor, it can be extracted from particle density estimation [WCW * 11], using streakline cores [WT10], or as rotational invariant vortex core [GST16]. This attractor can be regarded as a "sink-type DHT" [BDZG19], as it has only negative Lyapunov exponents. Note, that the definition of a DHT due to Ide et al. [ISW02] also includes this case, since their notion of hyperbolicity means non-zero Lyapunov exponents rather than saddle-like behavior. We therefore may refine the erroneous parallel vectors solution toward the attractor using the algorithm described in Section 4.2. In this case, integrating the decoupled system (Equation 9) is equivalent to integration of the original flow. The initial candidate is thus increasingly refined toward the ground truth with increasing time (to the right in Figure 5), while the point at the left time boundary is left unchanged. Integration-based methods are able to extract this attractor accurately near the beginning of the time domain, but require a dense computation of pathlines. Our method is only accurate at later time steps, but only requires a localized integration.

Unsteady Saddle Connectors
In 3D steady vector fields, the separatrices of saddle-type critical points can intersect in streamlines, which form connections between two different (heteroclinic) or one (homoclinic) saddle-type critical point. These are also called saddle connectors [TWHS03]. Similar connections can be formed by hyperbolic trajectories in unsteady vector fields [MW98]. They are typically formed at a time t i by manifolds of DHTs, which have come into existence at a time in the past (t < t i ) or in the future (t > t i ). For sufficiently long advection time, they can be observed as intersections of ridges in the forward and backward FTLE fields, which, however, do not represent DHTs. For illustration, we construct an analytical example, from  two stationary saddle-type linear fields, with instantaneous critical points at (0, 0) ⊤ and (1, 1) ⊤ . Both are made nonlinear with Gaussian window functions. The saddle at (1, 1) ⊤ is faded into and out of existence over a short time period, resulting in a short DHT, and one DHT that exists for all times. The DHTs and their manifolds are shown in Figure 6. The intersection curves of their manifolds are visible in the forward and backward FTLE fields computed at a time slice, where the shorter DHT does not exist. These unsteady equivalents of saddle connectors can be hyperbolic trajectories, but they do not need to be. However, they are not DHTs, since they are not generators of the time-dependent topology of the flow. Since they arise from non-local mechanisms, however, local extraction methods (such as ours) do usually not yield candidate lines for these false-positives (cf. Figure 8).

Convergence of Refined DHTs
We use a simple analytical model, for which the ground truth is known, to measure convergence of our DHT refinement, and compare it to the refinement scheme of Machado et al. [MSE13]. The vector field is defined componentwise as u i (x,t) = d i x i + A i sin(ω i t). Since it is periodic in time, the DHT for all times is obtained [ISW02] as Locations of vanishing acceleration, and thus the parallel vectors solution, are obtained, by straightforward calculation, as To avoid a possible influence of the numerically evaluated parallel vectors operator, we use this analytical representation instead. We fix the parameters d 1 = 3, d 2 = −3, ω 1 = 2, ω 2 = 3, A 1 = A 2 = 1, and sample the analytical field on a regular grid with 100 3 nodes over the domain pointwise mean square error to the ground truth (Figure 7a). For the computation of our DHT refinement, 100 fixed point iterations, together with the continuous SVD, took 28 ms on average, while 100 iterations of the bifurcation line refinement [MSE13] took 19 ms. However, our DHT refinement converges quickly after about 20 iterations, while the bifurcation line refinement takes 2000 iterations until a good solution is obtained. The MSE for both methods stays above a rather large value. With our DHT refinement, the reason for this are the initial values at the ends of the candidate line (Equation 12), while the bifurcation line refinement scheme cannot distinguish between the DHT and other path lines converging toward it at the end points. The convergence in the case of a nonlinear field is shown in Section 6.7. Next, we measure the influence of the time length of the initial candidate on the two algorithms (Figure 7b). For different amounts s of time, we sample the PV line over the time interval [−s/2, s/2], while keeping the sampling density equal to the previous experiment. Fixing the number of fixed point iterations at 100 for the DHT refinement, and the number of iterations for the bifurcation line refinement at 2000, we measure the distance of the center vertex to the ground truth, since both methods tend to be inaccurate at the endpoints. Since our algorithm is integration-based, the accuracy increases smoothly with the amount of time available along the candidate line. The bifurcation line refinement scheme, on the other hand, is generally unstable. Finally, we measure the impact of splitting the time interval into intervals of T = (t N −t 0 )s/ϑ for different s. This was introduced with s = 15 in Section 4.2 to avoid numerical errors. With fluctuations due to numerical errors, accuracy increases exponentially with s ( Figure 7c).

Convergence of FTLE Ridges
A small area of interest (I in Figure 13e) of the streak topology computed in the Cylinder Flow (see Section 6.8 for a detailed discussion) is shown enlarged in Figure 8. At moderate resolution, the ridge locations in the forward and backward FTLE fields cannot be reliably determined (Figure 8a). Only with much increased resolution (Figure 8c), the five ridges in the backward FTLE field (blue) become apparent. Increasing resolution even further (Figure 8d), leads to aliasing artifacts due to the Runge-Kutta 4 integrator with fixed step size used in our implementation. We also note, that only one of the five ridge intersections actually belongs to a DHT (Section 6.3), as detected by our streak topology extraction. Except near the seeding locations of the streak manifolds, where numer- ical integration time is zero, extracted ridges in the FTLE fields approach our streak manifolds for increasing resolutions of FTLE. Figures 8e-8h show distances between streak manifolds and ridges.

Stability under Perturbation
We analyze the stability of our DHT refinement when applied on a perturbed initial candidate, which investigates robustness against inaccurate candidate extraction, and applied on the perturbed ground truth DHT, which investigates how far away a candidate is allowed to be from the DHT for convergence. We employ two kinds of perturbations: symmetric perturbation, zero at the center, linearly increasing outward, and asymmetric perturbation with zero at the beginning, linearly increasing in forward time.
First, we use the same analytical model as in Section 6.4 over the time span [−2, 2], where the ground truth is known (Figure 9). The ground truth has Lyapunov exponents 3, −3, and a spatial range of 1.0 × 1.0. We introduce perturbations, that linearly vary between zero and (0.6, 0.6) ⊤ . In all cases, our DHT refinement is able to compensate for the perturbations only in directions where sufficient integration time is available. The errors of the streak manifolds are partially corrected by streak integration. Second, we perform the same analysis for the Cylinder Flow (Section 6.8), where we use the path of the instantaneous critical point behind the cylinder in the time interval [0, 2.5] as initial candidate, and use the DHT computed from it as ground truth. We computed its Lyapunov exponents as approximately 5.56663, −5.47087, and a spatial range of approximately 0.5 × 0.02. We impose, relative to the spatial scales, similar perturbations varying between zero and (0.06, 0.06) ⊤ . The same results as in the analytical case can be observed here (Figure 10). However, larger perturbations have caused divergence of the DHT computation in our experiments. While there are theoretical results on this convergence [JSW03], we leave further investigation at numerical datasets for future work.
The vector fields are defined on the spatial domain [0, 0] × [2, 2]. They contain a saddle-type critical point at approximately (1 + a (i) (t), 1 + a (i) (t)) ⊤ . Starting from t = 0, in u (1) , the critical point moves from (1.5, 1.5) ⊤ to (1, 1) ⊤ over 5 s, while the critical point in u (2) moves to (0.5, 0.5) ⊤ over 10 s. Both fields are symmetric in time, and they coincide over the time interval [−5 s, 5 s]. Thus, any local algorithm, that at most considers data on this interval, would yield the same result at t = 0 s. However, as Figure 11 shows, the location of the DHT, which starts to become visible using the FTLE over the time interval [−10 s, 10 s], differs tremendously in the two flows. We therefore conclude, that hyperbolic trajectories are not defined by local properties of the flow, but by the global dynamics. The convergence of the DHT refinement is shown in Figures 11d and 11e. After the first fixed-point iteration (i to ii, white to brown), the exact location is not yet reached, since the localized flow along the initial candidate (orange sphere) was used. The result of this is used to recompute the localized flow and start a second fixed-point iteration in our scheme (ii to iii), which quickly reaches the intersection of FTLE ridges.

Cylinder Flow
We now evaluate our method at a CFD simulation of a flow behind a cylinder. It was computed using the Gerris flow solver [Pop04] and is provided by Günther et al. [GGT17]. We compute DHTs and the time-dependent vector field topology over the first ten sec-onds of the dataset. We compare the different extraction methods for obtaining initial candidates, and compare the computation of DHTs with the bifurcation line refinement method by Machado et al. [MSE13]. The results are evaluated using the distance to the nearest intersection of ridges in the FTLE fields. Since extracting ridge lines or ridge surfaces in the space-time domain [BSDW12] requires a prohibitively large FTLE resolution (see Section 6.5), we extract those local maxima in the product of the two FTLE fields, that have a persistence greater than 0.5. The topological simplification was performed using TTK [TFL * 17]. Using an FTLE resolution of 2000 × 1000, this approach yields reasonably noise-free results (see Figures 12k-12m), except where the FTLE field does not exhibit sharp ridges or ridges are too close to each other.
Extracting critical points from the lab frame of reference, only results in short candidate lines (Figure 12a), because the saddles and nodes, that are periodically generated behind the cylinder, cancel each other out after a short amount of time in this frame of reference. Both our DHT refinement (Figure 12c) and the bifurcation line refinement [MSE13] (Figure 12b) are only able to make minor corrections, since both methods rely on sufficiently long candidate lines. Computing streak manifolds from these short segments misses most of the repelling LCS, but is already able to obtain large parts of the attracting LCS (Figure 12u).
Extracting critical points in the Galilean-invariant frame of reference defined by the feature flow field [MBES16], and in an optimal frame of reference [RG20] both result in missing line segments due to numerical noise in the area behind the cylinder. Therefore, we obtain initial robust segments by filtering, where det ∇u > −10, i.e., τ h =10, and integrate from their ends along the respective observer motions (Figures 12o and 12p), as long as integration stays in a hyperbolic region. While in the Galilean-invariant reference frame, integration for some candidate lines leaves the hyperbolic region before the domain boundary is reached, using an optimal reference frame, the domain boundary is always reached. An exception is the first hyperbolic trajectory, which is created at the beginning of the simulation. This hyperbolic trajectory stops existing before t=5 s, since the backward FTLE ridge belonging to its attracting manifold does not intersect with a forward FTLE ridge at this instance of time (II in Figure 13e). In the accompanying video this event can be observed at around t=4 s of the dataset.
While the candidate lines in the optimal frame of reference are more accurate, our DHT refinement reaches similar results in both cases, which very closely follow FTLE ridge intersections (Figures 13a-13d), except near the right domain boundary, where the ends of the initial candidates are reached. The bifurcation line refinement is very unstable near the cylinder, where the refined lines oscillate. Furthermore, it deviates more from the FTLE ridge intersections near the right domain boundary than our DHT refinement. Streak manifolds computed from the DHTs obtained from both approaches capture most of the attracting and repelling LCS (Figures 12v and 12w). Since initial candidates obtained using the optimal reference frame [RG20] are longest, they yield larger streak manifolds. Thus, this is the most accurate option, and we use this to obtain candidates for DHT refinement in our streak-based topology. our method, for the time steps t=5 s and t=7.5 s. While our streakbased topology closely matches ridges in the FTLE fields, critical points in the steady VFT only are accurate in regions, where the LCS are moving at constant speed. Separatrices in the steady VFT only follow FTLE ridges for short amounts of time, if at all. Since streamlines of a steady vector field cannot intersect, intersections of the repelling and attracting LCS cannot be captured by this concept (see Section 6.3 and Figure 8).

Convective Flow
We now consider a CFD simulation of bouyant air flow with two obstacles. This flow differs from the flow behind a cylinder, because it is confined to a closed container on a spatial domain of 0.1 m 2 with no-slip boundaries. We consider the time interval [49 s, 52 s] of this dataset. Much of its unsteady topology is generated by a slow moving DHT near the center of the domain, for which an initial candidate can be obtained by tracking critical points in the lab frame. As in the previous section, initial candidates obtained from the optimal frame of reference [RG20] yield the most accurate results (see Section 2 in the supplemental material). The streak manifolds in the space-time domain are shown in Figure 1. In Figure 14, we compare the streak topology obtained in this way with the steady VFT proposed by Rojo and Günther [RG20]. Again, in the steady VFT, neither do the critical points resemble LCS intersections, nor do separatrices follow LCS for longer integration times. On the other hand, in this dataset, many FTLE ridges are not captured by our approach. Some of these cases are investigated in Figure 15, where we seed pathlines across some of those missed FTLE ridges. We have found, that these either correspond to separation induced by shear, or other weak separation, and are thus false-positives in the FTLE field.

Performance
Computing LCS directly from the FTLE over the entire time domain of a dataset requires ridge extraction from densely evaluated FTLE fields for different starting times [BSDW12]. Our method, on the other hand, computes LCS as streamsurfaces in the spacetime domain seeded at an offset from locally extracted DHTs, resulting in a geometric representation of LCS that varies smoothly over time and is not affected by FTLE resolution. Unless reference frame optimization [RG20] is used for obtaining initial candidates, our method differs from that of Machado et al. [MBES16] only in the DHT refinement step, which has a comparable computational cost (see Section 6.4). As the authors have shown, their approach is about two orders of magnitude faster than a dense computation of FTLE fields, and, thus, so is our method.

Limitations
Since the DHT refinement is integration-based, it requires sufficiently long initial candidates, depending on the dataset. The sharpness of ridges in the FTLE fields computed within the time interval of an initial candidate can give an indication whether the candidate is sufficiently long. In turbulent flows, the extraction methods for obtaining initial candidates can miss features or only yield short segments. This could be partially overcome by enforcing temporal coherence in the local extraction methods, i.e., considering a larger time interval at increased computational cost.

Conclusion
We presented an approach to 2D time-dependent vector field topology, that similarly to the previous work of Machado et al. [MBES16] relies on local extraction of candidate lines for hyperbolic trajectories, from which distinguished hyperbolic trajectories are refined. Our algorithm for refining DHTs closely follows the work of Ide et al. [ISW02], as well as Branicki and Wiggins [BW09]. We have shown, that this approach is faster and more reliable than the previous method [MSE13], and that the resulting streak-based topology is more accurate than vector field topology in a steady frame of reference [RG20]. While the extension of the notion of a DHT to 3D is straightforward [BW09], they should be replaced by hyperbolic path surfaces in the 4D space-time domain [USE13]. We therefore want to treat this case in future work.