{'fraction_non_alphanumeric': 0.04021341764647762, 'fraction_numerical': 0.017834206106115996, 'mean_word_length': 4.2070739549839224, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hyperbolic and Elliptic Transport Barriers in Three-Dimensional Unsteady Flows

We develop a general theory of transport barriers for three-dimensional unsteady flows with arbitrary time-dependence. The barriers are obtained as two-dimensional Lagrangian Coherent Structures (LCSs) that create locally maximal deformation. Along hyperbolic LCSs, this deformation is induced by locally maximal normal repulsion or attraction. Along shear LCSs, the deformation is created by locally maximal tangential shear. Hyperbolic LCSs, therefore, play the role of generalized stable and unstable manifolds, while closed shear LCSs (elliptic LCSs) act as generalized KAM tori or KAM-type cylinders. All these barriers can be computed from our theory as explicitly parametrized surfaces. We illustrate our results by visualizing two-dimensional hyperbolic and elliptic barriers in steady and unsteady versions of the ABC flow.


Introduction
Detecting transport barriers is important in a number of areas, including geophysical flows [29], plasma fusion [25], reactive flows [24], and molecular dynamics [28]. For steady and temporally periodic flow models in these areas, classical dynamical systems theory identifies key invariant manifolds acting as phase space barriers [22]. Even in this well-understood setting, however, only specific examples of de facto barriers have been identified. Indeed, even for steady flows, no general approach to defining and locating multi-dimensional transport barriers has been available. A commonly used informal definition of barriers as surfaces with zero transverse flux is easily seen to be inadequate. Indeed, any randomly chosen surface of trajectories (material surface) admits zero normal flux [16].
A number of heuristic flow diagnostics have nevertheless been employed to infer transport barriers indirectly, mostly targeting two-dimensional flows (see [4]and [23] for reviews). These diagnostics can be highly effective for simple flows, especially if one tunes their threshold parameters to match expectations. However, when it comes to discovering more complex flows, the lack of an exact mathematical foundation renders these approaches problematic, manifested by frame-dependence, false positives and false negatives [2,18,14].
For exploration, decision making and forecasting, one would ideally require a general approach with solid mathematical foundations. Such a general approach, the geodesic theory of transport barriers, has recently emerged for two-dimensional unsteady flows [16,6,17]. Specifically, Haller and Beron-Vera [16] constructs transport barriers as curves most closely shadowed by least-stretching geodesics of the Cauchy-Green strain tensor derived from the flow map. As a further improvement, F. J. Beron-Vera and Haller [6], Haller and Beron-Vera [17] constructs transport barriers as distinguished geodesics of the Green-Lagrange strain tensors derived from the flow map. The objective of the present work is to extend these ideas to unsteady flows in three dimensions.
We are unaware of other approaches that would directly target transport barrier surfaces in multi-dimensional unsteady flows. Scalar fields associated with the flow map such as the finite-time Lyapunov exponents (FTLE) and finite-size Lyapunov exponents have been used as indicators of hyperbolic coherent structures [13,21,3,27,11,26]. Under certain conditions, select ridges of these fields can be rigorously related to repelling and attracting transport barriers [15,8,7,19].
Related methods exist for multi-dimensional coherence detection, including an ergodicity-based approach to visualizing elliptic structures in steady and temporally periodic flows [5], as well as a probabilistic approach to locating almost invariant sets in phase space [10]. Both methods seek coherent domains (ergodic components or almost invariant sets) via a modal truncation of an infinitedimensional operator (the Koopman or the Peron-Frobenius operator, respectively). This process then yields scalar fields whose topology is expected to reveal coherent sets. Specifically, in [5], lowindex eigenfunctions are proposed as indicators of dynamically distinct regions of the phase space. In [10], the left and right eigenfunctions of the second largest singular value of the Peron-Frobenius operator are thresholded to maximize the coherence ratio of pair of sets.
Our focus here is a direct variational construction of transport barriers as parametrized surfaces. Instead of starting with a particular mathematical quantity and arguing for its relevance in barrier detection, we start with a physically motivated question: What objective property makes transport barriers observable in physical and numerical experiments?
We put forward the same answer that has been well-tested in two-dimensional flows. Specifically, for a time-evolving surface M(t) to be an observed transport barrier, the following two properties should hold: T1 M(t) must be a material surface, i.e., a two-dimensional invariant manifold in the extended phase space of positions and time. This implies that barriers locally divide the phase space and have zero flux between their two sides.
T2 M(t) must impose locally extreme deformation on nearby sets of initial conditions. This is achieved either by locally maximal normal repulsion or attraction (hyperbolic barrier), or locally maximal tangential shear (shear barrier).
Properties T1-T2 provide an extension of the concept of a multi-dimensional Lagrangian Coherent Structure (LCSs) from the purely hyperbolic case treated in [15] to the general case. Solving the extremum problem described in T2 leads to vector fields of admissible unit normals for hyperbolic and shear barriers. It turns out that surfaces orthogonal to these admissible normal fields can only exist at locations where the helicity of these normal fields vanishes. Strain-and shear-helicity generically vanish on computable two-dimensional surfaces; these zero sets necessarily contain the transport barriers we seek. The intersection of transport barriers with select two-dimensional reference surfaces turns out to satisfy ordinary differential equations (strain and shear ODEs). These ODEs can be solved numerically, yielding parametrized reduced strainlines and reduced shearlines on the reference surfaces. Open reduced shearlines of zero helicity signal generalized jets (parabolic barriers), while closed reduced shearlines mark invariant tubes or invariant tori (elliptic barriers). Extracting such parametrized curves over a parametrized family of reference surfaces leads to explicitly parametrized two-dimensional transport barriers.
This construction applies to any three-dimensional flow with general time-dependence, and uncovers key barriers that shape tracer patterns over a finite time of observation. This time can be arbitrarily short or long: our approach, by construction, will locate barriers that best explain tracer patterns developing over the observational period chosen. Over longer time intervals, the same approach yields increasingly accurate approximations for classic hyperbolic and elliptic invariant manifolds, should those exist in the given flow.
We first illustrate these results on the steady and time-periodic ABC flows, which have welldefined steady and time-periodic transport barriers given by invariant manifolds. Even in these flows, we obtain new, explicit barrier surfaces that were previously only inferred from numerical images. Next, we consider a chaotically forced version of the ABC flow over a finite time interval. For this flow, transport barriers can only be constructed as temporally aperiodic material surfaces in the extended phase space. Remarkably, we obtain that select hyperbolic barriers and torus-type shear barriers continue to exist even in this fully aperiodic setting. The latter tori tori deform aperiodically in time, yet continue to provide sharp boundaries for coherent Lagrangian vortices. Indeed, they exhibit minimal deformation while nearby material elements in their exteriors stretch exponentially.

Set-up and notation
Consider the dynamical systeṁ with a smooth vector field v(x, t) defined over a time interval of length T , for locations x in a compact set D. We assume a finite time time interval in (1) since data sets obtained from physical measurements or numerical methods are only known over such intervals. Moreover, coherent structures of physical interest are typically transient in nature (eddies, hurricanes, etc). A trajectory of (1) starting from x 0 at time t 0 is denoted by x(t; t 0 , x 0 ). The flow map of (1) is then defined as which is as smooth in x 0 as the function v(x, t) in x.
Assuming that v is of C r with r ≥ 1, the flow gradient ∇F t t0 (x 0 ) can be computed. This linear mapping advects small initial perturbations ξ 0 to x 0 along the trajectory x(t, t 0 , x 0 ) to the evolved denotes the Cauchy-Green strain tensor, and < · , · > is the classic Euclidean inner product.
We will be interested in stationary values of total perturbation growth (3) over the time interval [t 0 , t 0 + T ]. These values are precisely the eigenvalues of the symmetric, positive definite matrix C t0+T t0 (x 0 ). The eigenvalues λ i (x 0 ) and their corresponding orthonormal eigenvectors ξ i (x 0 ) satisfy From now on, we restrict our discussion to an open set U of initial conditions where the eigenvalues of C t0+T t0 (x 0 ) are disjoint:

Three-dimensional transport barriers
Here we give a formal definition of transport barriers building on the properties T1-T2 described in the Introduction. According to T1, a time-dependent transport barrier M(t) must be a material surface, i.e., an invariant manifold in the extended phase space of the variables (x, t). This necessarily implies for any time t ∈ [t 0 , t 0 + T ]. As long as M(t 0 ) is a smooth surface, so is the surface M(t) for any fixed time t. The family M(t) is equally smooth in t by our smoothness assumption on (1). At an initial point x 0 ∈ M(t 0 ), let n 0 denote a unit normal to M(t 0 ). Then, as discussed in [15], a smoothly varying unit normal to M(t) along the trajectory x(t, t 0 , x 0 ) is given by For any initial point x 0 ∈ M(t 0 ) and initial unit normal n 0 to M(t 0 ) at x 0 , we define the normal repulsion ρ(x 0, n 0 ) of M(t) along the trajectory x(t; t 0 , x 0 ) as the normal component of the growth of n 0 under the linearized flow between times t 0 and t 0 + T [15]. Specifically, we have with the geometry illustrated in Fig. 1. Similarly, we define the tangential shear σ(x 0 , n 0 ) as the tangential component of the growth of n 0 under the linearized flow along the trajectory x(t; t 0 , x 0 ) between times t 0 and t 0 + T [16]. Specifically, we have with the geometry shown in Fig. 1.
We seek transport barriers as material surfaces that maximize normal repulsion or tangential shear with respect to perturbations to their tangent spaces. We do not insist on this maximizing property under all perturbations to the material surface: we only consider perturbations to their tangent spaces. This is because we seek a well-defined local directionality for the transport barrier, while in principle allowing for it to have a finite thickness. In other words, the barrier may a priori be a thin set of nearby, parallel surfaces.
if for any point x 0 ∈ M(t 0 ) and for any other material surfaceM(t 0 ) with x 0 ∈M(t 0 ) and with unit normaln 0 n 0 at x 0 , we have (ii) A material surface M(t) ⊂ R 3 is called an attracting hyperbolic LCS over the interval [t 0 , t 0 +T ], if for any point x 0 ∈ M(t 0 ) and for any other material surfaceM(t 0 ) with x 0 ∈M(t 0 ) and with unit normaln 0 n 0 at x 0 , we have (iii) A material surface M(t) ⊂ R 3 is called a shear LCS over the interval [t 0 , t 0 + T ], if for any point x 0 ∈ M(t 0 ) and for for any other material surfaceM(t 0 ) with x 0 ∈M(t 0 ) and with unit normaln 0 n 0 at x 0 , we have withn 0 denoting a unit normal toM(t 0 ) at the point x 0 .
(iv) A material surface M(t) ⊂ R 3 is called a transport barrier over the interval [t 0 , t 0 + T ], if it is either a hyperbolic or a shear LCS over [t 0 , t].
Remark 1. By (i) of Definition 1, hyperbolic LCS exhibit strictly maximal repulsion relative to all perturbations to their tangent spaces. This is to guarantee that any potential stretching in directions tangent to the LCS is smaller than in its normal direction. This condition guarantees both the observability and the robustness of a repelling hyperbolic barrier [15]. The same observation applies to attracting hyperbolic barriers in backward time.
Remark 2. By (ii) of Definition 1, the shear σ(x 0 , n 0 ) along a shear barrier M(t 0 ) is not required to be strictly maximal among all possible perturbations to the tangent space of M(t 0 ). As it turns out below, there are always two choices of the normal n 0 yielding the exact same, locally largest value of σ(x 0 , n 0 ) at any point x 0 . Shear barriers, therefore, exhibit strictly maximal shear only with respect to small enough perturbations of their normals. There will always exist a unique, finite perturbation to their normal yielding exactly the same locally maximal shear at x 0 . The two shear extrema at x 0 represent maximal shear with two different signs, which is disguised by the absolute value appearing in the definition of σ(x 0 , n 0 ).
The following theorem shows that for a material surface M(t) to be a transport barrier over [t 0 , t 0 + T ], its initial position must be orthogonal to a direction characterizing maximal repulsion or maximal shear. We use the notation T x0 M(t 0 ) for the tangent space of M(t 0 ) at a point x 0 .
Proof. See Appendix C.
Remark 3. The above necessary conditions for hyperbolic LCS have previously been obtained from slightly different considerations, along with examples illustrating their meaning in [15,7,9]. For an explicit example of how shear LCS can be found in three-dimensional, unsteady parallel shear flows, we refer to Appendix A.
Remark 4. Unlike in the two-dimensional case [16], the shear LCSs obtained in (iii) of Theorem 1 generally do not preserve their surface area under an incompressible flow map F t0+T t0 , even though they still preserve their enclosed volume (cf. Appendix B). This enables their use in detecting material footprints of commonly observed toroidal vortices, such as growing smoke rings.
Remark 5. A related recent paper [12] shows how quasi-invariant hyperbolic LCS can be used to compute a specific family of hyperbolic barriers (normally hyperbolic invariant manifolds) in steady flows of arbitrary dimension. Theorem 1 requires the initial position M(t 0 ) of a transport barrier to be orthogonal to ξ 3 (hyperbolic barrier) or to n ± (shear barrier). In general, if a two-dimensional surface is orthogonal to a three-dimensional vector field ρ(x), then any local parametrization p(s 1 , s 2 ) : U ⊂ R 2 → R 3 of the surface must satisfy the first-order quasi-linear system of PDEs This system of PDEs will only have a smooth solution through a given point x 0 if this point is contained in a transport barrier. To locate such barrier points, we now give computable necessary conditions for transport barrier locations. In stating these conditions, we will use the helicity H ρ (x) of a three-dimensional vector field ρ(x), defined as with × denoting the cross product. (i) Suppose that M(t) is a repelling hyperbolic LCS. Then at all points x 0 ∈ M(t 0 ), we must have (iii) Suppose that M(t) is a shear LCS. Consider the two vector fields Then at all points x 0 ∈ M(t 0 ), we must have for one choice of the sign in ±.
Proof. See Appendix D.
Remark 6. The problem of finding surfaces orthogonal to vector fields is locally equivalent to finding surfaces tangent to two smooth vector fields. For the existence of such tangent surfaces, the Frobenius Integrability Theorem provides a necessary condition. This can be shown equivalent to the zero helicity conditions described above (cf. Appendix D).
Theorem 3 provides specific scalar equations of the form (11) and (13) that a transport barrier M(t 0 ) must satisfy. Rather than solving these equations numerically, we locate the intersection curves of all potential transport barriers with a family of two-dimensional reference surfaces. Out of all these intersection curves, we then select the ones on which the appropriate helicity condition in Theorem 2 vanishes. This leads to the following result:

Theorem 3. [Necessary condition for intersections of transport barriers with reference surfaces]
Let M(t) ⊂ R 3 be a material surface over the time interval [t 0 , t 0 + T ]. Also, let Π(s 1 ) denote a smooth, one-parameter family of two-dimensional orientable surfaces with smooth normal vector fields n Π(s1) (x).
(i) Suppose that M(t 0 ) is a repelling hyperbolic LCS. Then the intersection curve γ s1 = Π(s 1 ) ∩ M(t 0 ) is a trajectory of a two-dimensional differential equation on Π(s 1 ), given by This trajectory must also satisfy the condition (ii) Suppose that M(t 0 ) is an attracting hyperbolic LCS. Then the intersection curve γ s1 = Π(s 1 ) ∩ M(t 0 ) is a trajectory of a two-dimensional differential equation on Π(s 1 ), given by This trajectory must also satisfy the condition is a trajectory of a two-dimensional differential equation on Π(s 1 ), given by for some choice of the sign in ±. This trajectory must also satisfy the condition with the same choice of the sign. If the trajectory γ s1 (s 2 ) is a closed orbit, then M(t) is an elliptic barrier.
Proof. See Appendix E.

Computation of transport barriers
Theorem 3 provides a practical algorithm for the computation of transport barriers as different types of LCSs in three-dimensional flows. The barriers can be reconstructed from their intersections with a family of orientable hypersurfaces.
In the simplest case, these hypersurfaces are just two-dimensional planes. For this case, we summarize below the extraction of hyperbolic LCS (generalized stable and unstable manifolds) and elliptic LCS (invariant cylinders and generalized KAM tori). Further hints on the numerical implementation of these algorithmic steps can be found in Appendix D.
H3 Integrate such reduced strainlines as long as the running average of |H ξ3 | stays below 0 .
H4 Filter the reduced strainline segments so obtained to find the ones that approximate the zero sets of H ξ3 most closely. Specifically, if the Hausdorff distance x − y two strainline segments γ 1 andγ s1 are smaller than a small threshold value d 0 , then discard either γ 1 orγ s1 .
H5 Vary the parameter s 1 in the definition of the plane family Π(s 1 ) to obtain uniform coverage of the domain of interest. Repeat H1-H4 for each s 1 . Obtain repelling hyperbolic LCSs by fitting a surface to the parametrized curve family γ s1 (s 2 ), as described in Remark 4.
SH3 Integrate such reduced shearlines as long as the running average of H n± stays below 0 .
SH4 Keep only reduced shearlines that form limit cycles.
SH5 Vary the parameter s 1 in the definition of the plane family Π(s 1 ) to obtain uniform coverage of the three-dimensional domain of interest. Repeat SH1-SH4 for each s 1 . Starting from a closed shearline γ s1min (s 2 ) on the plane Π(s 1min ) corresponding to the lowest value of the parameter s 1 , obtain a discretized approximation γ s1 (s 2 ) to a closed shear barrier by always selecting the closest closed reduced shearline in the planes Π(s 1 ) under increasing s 1 . Obtain elliptic LCS by fitting a surface to the curve-family γ s1 (s 2 ), as described in Remark 4.

Steady ABC flow
As a first example, we consider the steady ABC floẇ an exact solution of Euler's equation. We select the parameter values A = √ 3, B = √ 2 and C = 1.0. This well-studied set of parameter values yields the Poincare map shown in Fig. 2. We first use the theory developed here to construct the full two-dimensional transport barriers suggested by this Poincare map. Because the flow is steady, the transport barriers we seek are also invariant manifolds in the phase space, not just in the extended phase space.
We therefore only need to carry out the computational steps H1-H4 and SH1-SH4 of Section 4 to obtain intersection curves between barriers and a single reference plane Π. We then advect these intersection curves under the flow map to obtain the full two-dimensional barriers.
In Sections 5.2 and 5.3, we consider temporally periodic and aperiodic versions of (20) where this simple approach will no longer suffice. The present steady case is only considered here to provide a consistency check on a well-studied steady flow.

Elliptic LCSs in the steady ABC flow
We first perform the computational steps SH1-SH4 of the previous section in one of the vortical regions seen in the Poincare map plot of Fig. 2. The lower panels of the same figure show orbits the Poincare map in blue (color only in the online version), as well as closed reduced shearlines (green) obtained from the computational steps SH1-SH4 on the plane or two different integration times. In both cases, a uniform grid G 0 of 1000 × 1000 initial points was used to compute the Cauchy-Green strain tensor. The tolerance parameter in the computational step SH3 is chosen to be 0 = 10 −2 .
While the plots in Fig. 2 show curves in the (x, y) plane, the analysis is inherently threedimensional. Indeed, computing the vector n ± in equations (18)- (19) requires fully three-dimensional trajectory integration. By the periodic nature of the phase space, the shear LCS obtained from the advection of closed, reduced shearlines are two-dimensional tori. To bring out the toroidal nature of these barriers, we introduce new coordinates with the help of the approximate spatial core (x 0 (z), y 0 (z), z) obtained by advecting the vortical center point of Fig. 2. Using this center curve, we introduce the toroidal coordinate systemx where R i are positive constants. A nested family of invariant tori obtained from this transformation is shown in Fig. 3.

Repelling hyperbolic LCSs in the steady ABC flow
We compute the Cauchy-Green strain tensor C 3 0 over a 500 × 500 grid on the plane z = 0. The tolerance parameter in the computational step H3 is chosen to be 0 = 10 −4 . Under this tolerance level, intersections or repelling hyperbolic LCSs with the z = 0 plane are shown in the left panel of Fig. 5, obtained as parametrized curves. Also shown is the vertical line y = 0.5, along which we compute the strain helicity H ξ3 (cf. the right panel of Fig. 5). The latter figure illustrates that the reduced strainlines shown in the left panel are indeed intersections of the z = 0 plane with repelling hyperbolic LCSs in the sense of Definition 1.
To illustrate the dynamical impact of the barrier surface emanating from the reduced strainlines, we select one of these strainlines (shown in green in the left panel of Fig. 6). We perturb this reduced strainline segment in the x direction by ±0.01 to obtain the blue and red curves shown in the same panel. We then advect all three parameterized curves from t 0 = 0 to t 0 + T = 3 to obtain the surfaces shown in the left panel of Fig. 6. Note that the blue and red curves have noticeable upward and downward z-drifts, respectively, while the surface evolving from the green reduced strainline has no z-drift.

Time-periodic ABC flow
We now consider a temporally periodic version of the ABC flow, given by The first return map to the plane z = 0 is now a non-autonomous map. Therefore, to gain insight into the flow from classical tools, only a fully three-dimensional temporal Poincare map can be used. This would result in spatially scattered points, as opposed to the sharply defined shear and hyperbolic barriers that we will continue to obtain form our approach.

Elliptic LCSs in the time-periodic ABC flow
We compute the Cauchy-Green strain tensor C 30π 0 over a 500 × 500 grid in the plane z = 0. This integration length is equal to 15 iterations of the temporal Poincare map F 2π 0 . The tolerance parameter in the computational step SH3 is chosen to be 0 = 10 −2 . Fig. 7 shows the closed reduced shearlines we find as limit cycles of the equation (18). The shear-helicity zero distribution along these curves is similar to that in the steady case, and hence is omitted here for brevity.
Next, we iterate the outermost closed reduced shearline under the Poincare map F 2π 0 . The result is a two-dimensional invariant torus for F 2π 0 , shown in Fig. 7 under the embedding (21). This torus is an intersection of a three-dimensional invariant torus of the full, spatially and temporally periodic flow (defined over the toroidal phase space T 4 ) with the t = 0 hyperplane.
We now illustrate the barrier property of the three-dimensional torus represented by the elliptic LCS of Fig. 7. To this end, we advect two initial conditions from the interior of the two-dimensional   torus starting from the z = 0 reference plane, and two other initial conditions from the exterior of this torus within the same plane. These four initial conditions are placed on the grey circle shown in the left panel of Fig. 7. The center of this circle is on the outermost torus barrier, and is advected as a blue trajectory. As seen in the right panel of Fig. 7, the blue trajectory indeed remains confined to a quasi-periodically deforming transport barrier in phase space. This barrier keeps both the red and the yellow initial conditions from spreading. In contrast, the green and black initial conditions launched from outside the outermost torus barrier develop large excursions over time.

Repelling hyperbolic LCSs in the time-periodic ABC flow
We now compute repelling hyperbolic barriers for the time-periodic ABC flow using the slicing method described in steps H1-H5 of section 4.1. We select the discrete family of 21 planes Π(s 1 ) = (x, y, z) ∈ [0, 2π] 3 : z = s 1 , s 1 = 0.00, 0.005, 0.01, . . . , 0.1, with a 500 × 500 uniform grid G 0 within each such plane. Over the time interval ranging from t 0 = 0 to t 0 + T = 4.0, we carried out the procedure outlined in steps H1-H5 of section 4. In step H2, the grid G 1 was chosen as 600 × 10 in (x, y), and the helicity parameter was chosen as 0 = 10 −4 . The filtered reduced strainlines obtained from H1-H4 on the Π(0) plane are shown in the left panel of Fig. 9. The right panel of the same figure shows the reconstructed barrier surface by performing step H5 across the plane family Π(s 1 ) and interpolating smooth surfaces over the resulting reduced strainline segments.

Chaotically forced ABC flow
Here we consider a temporally aperiodic version of the ABC flow, given by the equations x = (A + F (t)) sin z + C cos y, with F (t) representing a chaotic signal. The signal is generated by a trajectory close to the strange attractor of a periodically forced and damped Duffing oscillator (see Fig. 10). The temporally aperiodic flow (23) admits neither a well-defined spatial nor a well-defined temporal autonomous first return map. Therefore, the simplified barrier visualization methods used for elliptic barriers in the steady (Fig. 3) and time-periodic (Fig. 7) ABC flows are no longer applicable.

Elliptic LCSs in the chaotically forced ABC flow
We compute the Cauchy-Green strain tensor C 100 0 over a 500 × 500 grid in each member of the plane family The forthcoming computations were carried out in a parallelized fashion over the 150 s 1 -slices defined in The closed reduced strainlines obtained form SH1-SH4 on the Π(0) plane are shown in the upper left panel of Fig. 11. The tolerance parameter in the computational step SH3 is chosen to be coherence property is remarkable, as illustrated by Fig. 12. In this figure, a circle of one million initial conditions is selected as a perturbation to the smaller diameter of the torus barrier. Just after an advection time of t = 13.0, the ring quickly loses all its coherence, stretching and folding by a large amount in a visibly chaotic fashion. Fig. 13 shows the same type of verification of the optimality of the barrier that we employed in Fig. 8 for the time-periodic ABC flow. Again, tracers launched inside the barrier remain confined to the interior of the barrier, while tracers launched slightly outside the barrier exhibit large excursions.

Repelling hyperbolic barriers in the chaotically forced ABC flow
Finally, we compute repelling hyperbolic barriers for the chaotically forced ABC flow using steps H1-H5 of section 4. The Cauchy-Green strain tensor C 5 0 is computed over the same plane family used in section 5.2.2 for the time-periodic case. The grids G 0 and G 1 , as well as the admissible upper bound 0 on the helicity norm, are also selected the same as in section 5.2.2. Fig. 14 shows the final result, the set of extracted repelling hyperbolic barriers in the chaotically forced case.

Conclusions
We have developed a unified theory of transport barriers for three-dimensional unsteady flows. The barriers are attracting, repelling or shear LCS, which are constructed to prevail as skeletons of material deformation over a given finite time interval of observation. Out of general shear LCSs, elliptic LCSs are distinguished by their tubular or toroidal topology. Our approach renders all these LCSs as explicitly parametrized surfaces with precisely understood impact on neighboring material blobs. This is to be contrasted with alternative approaches that seek the domains separated by transport barriers from various indicator functions, without yielding specific dynamical information or a parametrization for the barrier itself.
Our approach closely reproduces known transport barriers in the steady ABC flow, and provides similar results for time-periodic and time-aperiodic version of the same flow. Remarkably, generalized   KAM tori (Lagrangian vortex rings) continue to exist in the general aperiodic case, providing sharp boundaries for coherent toroidal islands in an otherwise chaotic flow. Although similar notions of multi-dimensional hyperbolic LCSs have been used earlier [15], the theory of elliptic LCS as well as the computational methodology developed here for all types of LCSs is new. Our notion of an ellkiptic LCS extends that arising in the two-dimensional theory of shear barriers developed in [16]. This two-dimensional theory has identified highly coherent Lagrangian eddies from satellite-mapped surface velocities in the Agulhas leakage of the Southern Ocean [6]. The direct analogy with the two-dimensional theory promises similar results in the application of the present techniques to three-dimensional numerical or experimental flow data.
Appendix A. Shear LCSs in unsteady, three-dimensional parallel shear flows Consider the flowẋ (t) = u(z, t), where the dependence of u, v, and w on their arguments is smooth but otherwise arbitrary. Trajectories of (A.1) satisfy We introduce the functions a(z 0 , t 0 , T ) and b(z 0 , t 0 , T ) as Suppressing the arguments of a and b, we obtain the Cauchy-Green strain tensor in the form We now show that the planes z = k = const. are shear LCSs in the sense of Definition 1, as obtained from an application of Theorem 1. To do this, we use an expression for the angle φ that the vector v ± := ξ 2 × n ± encloses with the vertical planes z = k for general 3D flows.

(A.5)
Proof. The two sides of equation (A.5) represent two different ways of computing v ± , C t0+T The left-hand side is computed using the matrix elements of C t0+T t0 . To compute the same quantity and arrive at the quantity on the right-hand side, recall first that n ± = αξ 1 ± βξ 3 , where Hence which proves the lemma.
Since ξ 2 is always orthogonal to n ± , to show that z = k is a shear LCS, it remains to argue that v ± is also tangent to z = k. Since the flow is incompressible, we conclude that λ 1 λ 3 = 1. As a result, the right-hand side of the angle formula in Lemma 1 is one. Using our specific form of the Cauchy-Green strain tensor (A.4), the angle formula (A.5) becomes cos φ a 2 + b 2 cos φ + 2 sin φ (a cos θ + b sin θ) = 1, (A.7) which has φ = 0 as a solution. Therefore, we conclude that both ξ 2 and v ± = n ± × ξ 2 are tangent to the plane z 0 = k, which is therefore a shear LCS provided that a(z 0 , t 0 , T ) = 0 or b(z 0 , t 0 , T ) = 0.

Appendix B. Evolution of LCS surface area
We consider how the surface area of an LCS changes under the flow map. We have the following general result where p si := ∂p ∂si .
Proof. Since F t0+T t0 • p is a parameterization of M(t 0 + T ), the vectors ∇F t0+T t0 p si are tangent to F t0+T t0 (U ). The advected surface area, by definition, is then (ii) Suppose that M(t) is an attracting hyperbolic LCS. Then we have (iii) Suppose that M(t) is a shear LCS. Then we have In the special case of a volume-preserving flow, we have Proof. We proof the result for (iii), as cases (i) and (ii) are similar. The shear vector field n ± has unit length, and hence p s1 × p s2 , p s1 × p s2 = p s1 × p s2 2 . Observe that Using the definition of n ± , one sees that n ± , C t0+T Proposition 1 shows that the final surface area along a shear LCS in incompressible flow is obtained by integrating the initial surface element p s1 × p s2 weighted by 4 √ λ 2 . In Appendix A, we showed that λ 2 = 1 golds globally in space and time, and hence the corresponding shear LCS surface area is conserved. For the steady ABC flow, we find that λ 2 computed over elliptic LCSs oscillates around one. Fig. Appendix B shows this along a specific closed reduced shearline.

Appendix C. Proof of Theorem 1
We first need the following proposition which allows us to compute the normal repulsion and tangential shear in terms of C t0+T t0 and n 0 . Proposition 2. The quantities ρ and σ can be expressed as Proof. The proof of the first identity can be found in [15]. As for the second identity, we make use of the following formula for the unit normal of a multi-dimensional material surface [15]: We then obtain which proves the proposition.
We are now ready to prove Theorem 1. Let S be a hyperbolic LCS and x 0 a point on S. Thus, by definition, for any other surfaceS passing through x 0 such that T x0S = T x0 S, the normal repulsion rate ρ(x 0 , nS 0 ) alongS at x 0 is smaller than the normal repulsion rate ρ(x 0 , n S 0 ) along S at x 0 . Thus, at each point x 0 on S, the quantity ρ(x 0 , n 0 ) is maximized with respect to changes in n 0 . Thus, we want to solve the following constrained optimization problem: maximize ρ(x 0 , n 0 ) with respect to n 0 under the constraint that n 0 = 1. To this end, expand n 0 in the eigen-basis {ξ 1, , ξ 2 , ξ 3 } of the Cauchy-Green strain tensor C t0+T t0 (x 0 ) : We then have and setting the gradient of ρ(x 0 , n 0 ) proportional to the gradient of the constraint n 0 = 1 yields − 2ρ 2 λ 2 λ 3 n 1 = γn 1 , − 2ρ 2 λ 1 λ 3 n 2 = γn 2 .
− 2ρ 2 λ 1 λ 2 n n = γn 3 , where γ is an appropriate constant. Thus, assuming that λ 1 > 0, it follows that two of the coordinates n 1 , n 2 , n 3 must be zero. Therefore, the only extremum directions of the normal repulsion rate are the eigenvectors ξ i of the Cauchy-Green strain tensor. Since ρ(x 0 , ξ i ) = √ λ i , it follows that n 0 = ξ 3 is the global maximum of the normal repulsion ρ with respect to changes in n 0 . Therefore, a repelling hyperbolic LCS is necessarily orthogonal to ξ 3 . The proof of statement (ii) is analogous, but involves the global minimum of the normal repulsion rate.
The prove statement (iii), we need to find the maximizing normal directions n 0 of the tangential shear σ(x 0 , n 0 ) under the constraint that n 0 = 1. We again represent n 0 in the eigen-basis {ξ 1, , ξ 2 , ξ 3 } of C t0+T t0 (x 0 ) as and seek to maximize σ(n 1, , n 2 , n 3 ) := σ(x 0, n 0 ) subject to the constraint Setting the gradient of σ(n 1 , n 2 , n 3 ) proportional to the gradient of the constraint n 2 i = 1 leads to the equations 1 2σ which can also be rewritten as or, equivalently, Assume now that for some index i, we have n i = 0, i.e., the i th component of the unknown normal vector n 0 of the shear LCS is nonzero. In that case, we must have where the constant K > 0 is the same for any choice of i. Taking the square of equation (C.6) gives and hence there can be at most two distinct λ i values for which (C.6) holds. Also note that there has to be at least two distinct λ i values for which (C.6) holds, otherwise n 0 would be an eigenvector, and hence a local minimizer of the tangential shear . We conclude that there exist precisely two eigenvalues, λ k and λ l > λ k , of C t0+T t0 that satisfy (C.6). All eigenvalues of C t0+T t0 (x 0 ) are simple by assumption. Therefore, by our argument above, n 0 must be of the form n 0 = n k ξ k + n l ξ l , n 2 k + n 2 l = 1. (C.7) Substituting the expression (C.7) into (C.6) with i = k and i = l, eliminating the common constant 2σγ from the resulting two equations, and using the notation we obtain the system of equations for the unknowns a 2 and b 2 . The solution of these equations is given by Proposition 3. Let v be a smooth vector field in R 3 and S a surface orthogonal to v. Then for any x ∈ S, the helicity of vmust vanish, i.e., Consider an open neighborhood D ⊂ S of x in S . By Stokes' Theorem, we have that Since v is orthogonal to S, the integral on the right-hand side of (D.3) is zero. Thus, since v = v, n n, we haveˆD Since D was arbitrary, H v must vanish on S.
Theorem 2 then follows directly from Theorem 1 and Proposition 3.

Appendix D.2. Relation to Frobenius Integrability
We can rephrase the problem of computing a surface orthogonal to ξ 3 for repelling hyperbolic, ξ 1 for attracting hyperbolic, and n ± for shear LCS as finding surfaces tangent to {ξ 1 , ξ 2 }, {ξ 2 , ξ 3 } and {ξ 2 , n ± × ξ 2 }, respectively. The problem of finding surfaces tangent to two specified vector fields is then related to the Frobenius Integrability Theorem [1].
As a special case, this theorem states that if X and Y are two vector fields in R 3 , then necessary and sufficient conditions for the existence of a foliation of R 3 by surfaces tangent X and Y is that In our context, we do not seek to find global foliations of R 3 , but only individual, isolated surfaces. Nevertheless, as we show in this section, the Frobenius condition (D.5) is still a necessary, albeit not sufficient condition for the existence of such surfaces.
If X, Y and Z are smooth vector fields in R 3 , consider both the Frobenius and helicity conditions We show in Proposition 4 that the zero sets of F and H Z coincide, provided that X, Y and Z are pairwise orthogonal vector fields. Combined with Theorem 2, Proposition 4 then shows that (D.5) is a necessary condition. We make use of the following lemma, which states that the zero sets of F and H Z are invariant under a nonlinear rescaling of the three vector fields involved.
Proof. By definition, we have that However, X T Y = Y T X = 0 by our orthogonality assumption, and hence As for the claim on the helicity, note that Finally, using the identity (a × b) · c = (b × c) · a with a = ∇ × Z, b = Y , and c = X, we obtain that (D.11) is equivalent to ∇ × Z, Z = 0 (D. 16) as claimed.
Remark 9. Frobenius Integrability Theorem applied to the existence of tangent foliations provides a necessary and sufficient condition. By contrast, the zero helicity condition (and its equivalent Frobenius condition) are only necessary conditions for the existence of isolated surfaces normal to a vector field v. For example, let v(x, y, z) = (y, z, x). Then H v (x, y, z) = ∇ × v, v = −y − z − x, which has a plane as its zero set, but this plane is not orthogonal to v. Thus H v (x, y, z) = 0 is not sufficient for the existence of a surface orthogonal to v.
In the 3D setting of the present paper, we cannot find a linear transformation Ω that would globally fix the orientation of ξ 3 relative to ξ 1 . One therefore does not have simply two globally defined vector fields to integrate. Rather, one starts the integration and has to consider at each point along a reduced shearline all four possible directions, resulting in four possible choices of n 0 at that point. The correct vector to select is the one that results in a smooth reduced shearline (as the transverse intersection of a shear LCS with the reference plane). The method used in this paper to select the correct vectors at each time step is: 1) At the initial condition r(0) in the reference plane Π(s 1 ), compute the vectors ξ 1 (r(0)) and ξ 3 (r(0)).