1. Introduction
1.1. Vortex rings
Vortex rings are ubiquitous flow phenomena in both applied and theoretical settings, with applications including sound generation, transport, mixing and vortex interactions (Shariff & Leonard Reference Shariff and Leonard1992). In geophysical settings, vortex rings can be used to model entrainment and dispersion in particle clouds (Bush, Thurber & Blanchette Reference Bush, Thurber and Blanchette2003). They play important roles in the initial jets of volcanic eruptions (Taddeucci et al. Reference Taddeucci, Alatorre-Ibarguengoitia, Palladino, Scarlato and Camaldo2015) and the transport of contaminated sediments disposed of in open-water settings (Ruggaber Reference Ruggaber2000). In biomechanical settings, vortex rings have been observed in the motions of blood in the human heart (Arvidsson et al. Reference Arvidsson, Kovács, Töger, Borgquist, Heiberg, Carlsson and Arheden2016) and in the propulsive motion of oblate medusan jellyfish (Dabiri Reference Dabiri2005). Remarkably, separated vortex rings augment dandelion seed dispersal by prolonging flight through drag enhancement (Cummins et al. Reference Cummins, Seale, Macente, Certini, Mastropaolo, Viola and Nakayama2018). In aerodynamic settings, vortex rings are responsible for the so-called vortex ring state, which negatively impacts lift in helicopters (Johnson Reference Johnson2005) and the performance of offshore wind turbines (Kyle, Lee & Früh Reference Kyle, Lee and Früh2020). In experimental and numerical settings, the formation and pinch-off of vortex rings are of particular interest in jet flows involving nozzles and orifices (Gharib, Rambod & Shariff Reference Gharib, Rambod and Shariff1998; Mohseni, Ran & Colonius Reference Mohseni, Ran and Colonius2001; Krueger & Gharib Reference Krueger and Gharib2003; O'Farrell & Dabiri Reference O'Farrell and Dabiri2014; Limbourg & Nedić Reference Limbourg and Nedić2021).
Vortex rings are also associated with complex instabilities and dynamics that relate more generally to the sustenance of turbulence. Flow instabilities in vortex rings depend primarily on the core vorticity distribution, the circulation Reynolds number 
$( Re_\varGamma = \varGamma / \nu )$ and the slenderness ratio 
$( \delta = a / R )$ (Balakrishna, Mathew & Samanta Reference Balakrishna, Mathew and Samanta2020). Here, 
$\varGamma$ is the circulation, 
$\nu$ is the kinematic viscosity, 
$a$ is the core radius and 
$R$ is the ring radius. We focus on the evolution of thin-cored vortex rings with Gaussian core vorticity profiles, no swirl and centroids (
$Z$) that propagate along the 
$z$-axis. In cylindrical coordinates (
$r, \theta, z$), this initial vorticity profile is written as
\begin{equation} \omega_\theta (r, z; t = 0) ={\pm} \frac{\varGamma_0}{{\rm \pi} a_0^2} \exp \left(-\frac{(z - Z_0)^2 + (r - R_0)^2}{a_0^2}\right)\!, \end{equation}
where subscripts 
$({\cdot })_0$ denote parameter values at 
$t = 0$ and the sign of 
$\omega _\theta$ dictates the propagation direction. Since Gaussian vortex rings only satisfy the governing equations with infinitesimal core thickness, they initially undergo a rapid period of equilibration in which vorticity is redistributed throughout the core (Shariff, Verzicco & Orlandi Reference Shariff, Verzicco and Orlandi1994; Archer, Thomas & Coleman Reference Archer, Thomas and Coleman2008; Balakrishna et al. Reference Balakrishna, Mathew and Samanta2020). Following instability growth, transition is often marked by the development of secondary vorticity in a halo around the core vorticity (Dazin, Dupont & Stanislas Reference Dazin, Dupont and Stanislas2006; Bergdorf, Koumoutsakos & Leonard Reference Bergdorf, Koumoutsakos and Leonard2007; Archer et al. Reference Archer, Thomas and Coleman2008). During turbulent decay, the shedding of secondary vortex structures to the wake can result in a stepwise decay in circulation (Weigand & Gharib Reference Weigand and Gharib1994; Bergdorf et al. Reference Bergdorf, Koumoutsakos and Leonard2007).
Stability analyses of thin vortex rings are often (classically) formulated in terms of asymptotic expansions in 
$\delta$ (Widnall, Bliss & Tsai Reference Widnall, Bliss and Tsai1974; Widnall & Tsai Reference Widnall and Tsai1977; Fukumoto & Hattori Reference Fukumoto and Hattori2005). Infinitesimally thin vortex rings (
$\delta \to 0$) are neutrally stable (Shariff & Leonard Reference Shariff and Leonard1992). For rings with finite thickness (
$\delta > 0$), the curvature instability occurs at first order in 
$\delta$ and the elliptic instability occurs at second order in 
$\delta$. The curvature and elliptic instabilities occur at short wavelengths and arise due to parametric resonance between Kelvin waves with core azimuthal wavenumbers separated by one and two, respectively (Fukumoto & Hattori Reference Fukumoto and Hattori2005; Hattori, Blanco-Rodríguez & Le Dizès Reference Hattori, Blanco-Rodríguez and Le Dizès2019). The curvature instability is attributed to a dipole field produced by the vortex ring curvature (Fukumoto & Hattori Reference Fukumoto and Hattori2005; Blanco-Rodríguez et al. Reference Blanco-Rodríguez, Le Dizès, Selçuk, Delbende and Rossi2015; Blanco-Rodríguez & Le Dizès Reference Blanco-Rodríguez and Le Dizès2017). By contrast, the elliptic instability is attributed to a quadrupole field generated by straining induced by the ring or some external source (Fukumoto & Hattori Reference Fukumoto and Hattori2005; Blanco-Rodríguez et al. Reference Blanco-Rodríguez, Le Dizès, Selçuk, Delbende and Rossi2015; Blanco-Rodríguez & Le Dizès Reference Blanco-Rodríguez and Le Dizès2016).
This elliptic instability acts to break up elliptic streamlines and is key to the development of three-dimensional transitional and turbulent flows (Kerswell Reference Kerswell2002). In the context of vortex rings (or, more generally, strained vortices), it is sometimes called the Moore–Saffman–Tsai–Widnall instability (Fukumoto & Hattori Reference Fukumoto and Hattori2005; Chang & Llewellyn Smith Reference Chang and Llewellyn Smith2021) based on the initial investigations of Moore & Saffman (Reference Moore and Saffman1975) and Tsai & Widnall (Reference Tsai and Widnall1976). The elliptic instability dominates the curvature instability for thin Gaussian vortex rings without swirl. However, the curvature instability becomes increasingly important for vortex rings with increasing 
$Re_\varGamma$ and decreasing 
$\delta$, as well as in vortex rings with swirl (Blanco-Rodríguez & Le Dizès Reference Blanco-Rodríguez and Le Dizès2017; Hattori et al. Reference Hattori, Blanco-Rodríguez and Le Dizès2019).
While interesting in their own right, thin vortex rings often form canonical building blocks of more complex turbulent flows. Modified vortex geometries, such as elliptic vortex rings (Cheng, Lou & Lim Reference Cheng, Lou and Lim2016, Reference Cheng, Lou and Lim2019) and trefoil knots (Yao, Yang & Hussain Reference Yao, Yang and Hussain2021; Zhao et al. Reference Zhao, Yu, Chapelier and Scalo2021), provide alternative means of probing vortex dynamics and interactions. Collisions between vortex rings and other vortex rings, walls and free surfaces are also commonly studied to investigate mechanisms underlying the turbulent cascade and the generation of small scales (see Mishra, Pumir & Ostilla-Mónico (Reference Mishra, Pumir and Ostilla-Mónico2021) for a review). These mechanisms can be characterized using a variety of collision geometries, including head-on collisions (Cheng, Lou & Lim Reference Cheng, Lou and Lim2018; McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2018, Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021), inclined collisions (Kida, Takaoka & Hussain Reference Kida, Takaoka and Hussain1991; Yao & Hussain Reference Yao and Hussain2020a,Reference Yao and Hussainc) and axis-offset collisions (Zawadzki & Aref Reference Zawadzki and Aref1991; Smith & Wei Reference Smith and Wei1994; Nguyen et al. Reference Nguyen, Phan, Duong and Le2021), among others. Boundary layers play an important role in vortex–wall interactions (e.g. by causing rebounding events) (Walker et al. Reference Walker, Smith, Cerra and Doligalski1987) and interactions with free surfaces can often be understood in terms of image vortices (Archer, Thomas & Coleman Reference Archer, Thomas and Coleman2010).
Here, we focus on head-on collisions between identical vortex rings of opposite circulation, which have been noted for their rapid enstrophy production (Lu & Doering Reference Lu and Doering2008; Ayala & Protas Reference Ayala and Protas2017; Kang, Yun & Protas Reference Kang, Yun and Protas2020). They have been classically studied in the contexts of the formation of smaller rings through vortex reconnection and the formation of turbulent clouds at high 
$Re_\varGamma$ (Oshima Reference Oshima1978; Lim & Nickels Reference Lim and Nickels1992; Chu et al. Reference Chu, Wang, Chang, Chang and Chang1995). Many recent investigations have focused particularly on the mechanisms (e.g. instabilities) underlying these transitional and turbulent processes (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2018, Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021).
For the head-on vortex ring collisions under consideration, the elliptic instability competes and interacts with the longer-wavelength Crow instability. The Crow instability (Crow Reference Crow1970) is associated with the mutual interaction of perturbed counter-rotating vortices, which, in the linear regime, locally displaces the vortices without modifying their core structures (Leweke, Le Dizès & Williamson Reference Leweke, Le Dizès and Williamson2016). Mishra et al. (Reference Mishra, Pumir and Ostilla-Mónico2021) provides a focused review of vortex ring collisions in the context of the these instabilities. For collisions at relatively low Reynolds numbers, the Crow instability can lead to the pinch-off of secondary vortex rings via local reconnections. At higher Reynolds numbers, the elliptic instability favours rapid disintegration of the vortex rings into a turbulent cloud.
McKeown et al. (Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020) proposed that iterative elliptic instabilities between successive generations of antiparallel vortices can mediate the turbulent cascade in head-on vortex ring collisions. Mishra et al. (Reference Mishra, Pumir and Ostilla-Mónico2021) also observed that the elliptic instability tends to dominate at high 
$Re_\varGamma$, although this behaviour is also sensitive to the slenderness ratio and vorticity distribution. In a different configuration involving symmetrically perturbed antiparallel vortices, Yao & Hussain (Reference Yao and Hussain2020b) attributed the turbulent cascade at high 
$Re_\varGamma$ to an avalanche of successive vortex reconnections. In general, Ostilla-Mónico et al. (Reference Ostilla-Mónico, McKeown, Brenner, Rubinstein and Pumir2021) found that collisions between counter-rotating vortices are indeed highly sensitive to the geometry of their configuration. They particularly found that the mechanisms mediating the cascade bear resemblance to the reconnection scenario (Yao & Hussain Reference Yao and Hussain2020b) when the vortices are nearly perpendicular, whereas they are more reminiscent of the iterative elliptic instability scenario (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020) when the vortices are more acutely aligned. These recent works share two common themes, (i) that the mode of transition and the formation of a cascade are sensitive to the details of the initial flow configuration and (ii) that the interplay between relevant instabilities is simultaneously important to the flow physics and difficult to capture.
1.2. Velocity gradients and vortices
The elliptic instability, which typically dominates head-on collisions between the vortex rings of interest at high 
$Re_\varGamma$ (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021), is associated with elliptic streamlines (Kerswell Reference Kerswell2002). This generic feature of strained vortical flows can be used to characterize the elliptic instability, which is typically difficult to discern in the complex interactions of multiscale vortices (Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021; Ostilla-Mónico et al. Reference Ostilla-Mónico, McKeown, Brenner, Rubinstein and Pumir2021). Given the inherent complexity of turbulent flows, the geometry of local streamlines provides a relatively simple and interpretable means for characterizing flow features (e.g. vortices).
The instantaneous trajectory of a materially advecting fluid particle follows the streamlines, which are frame dependent. At a critical point, e.g. in a frame advecting with the particle, the velocity gradient tensor (VGT), 
$\boldsymbol{\mathsf{A}} = \boldsymbol {\nabla } \boldsymbol {u}$, determines, to linear order, the local structure of streamlines (Perry & Fairlie Reference Perry and Fairlie1975; Perry & Chong Reference Perry and Chong1987; Chong, Perry & Cantwell Reference Chong, Perry and Cantwell1990). The scale-invariant shape of local streamlines is captured by normalizing the VGT as 
$\tilde{\!\boldsymbol{\mathsf{A}}} = \boldsymbol{\mathsf{A}}/A$ (Girimaji & Speziale Reference Girimaji and Speziale1995; Das & Girimaji Reference Das and Girimaji2019), where 
$A = \left\Vert \boldsymbol{\mathsf{A}} \right\Vert_{\rm F} = {\rm tr}( \boldsymbol{\mathsf{A}}^{\rm T} \boldsymbol{\mathsf{A}})^{1/2}$ is the Frobenius norm of the VGT, 
$(\cdot)^{\rm T}$ represents the transpose and, unless otherwise stated, non-bold versions of bold tensor quantities represent their Frobenius norms. This normalized VGT has been used to investigate the scalings, forcings and non-local features of the VGT dynamics (Das & Girimaji Reference Das and Girimaji2019, Reference Das and Girimaji2020a, Reference Das and Girimaji2022) and a similar analysis of vorticity gradients has been used to classify the geometry of local vortex lines (Sharma, Das & Girimaji Reference Sharma, Das and Girimaji2021).
The principal invariants of 
$\tilde {\!\boldsymbol{\mathsf{A}}}$ instantaneously characterize local streamline topologies and geometries (Chong et al. Reference Chong, Perry and Cantwell1990; Das & Girimaji Reference Das and Girimaji2019, Reference Das and Girimaji2020a,Reference Das and Girimajib). They are given by
\begin{equation} p_{\rule{0pt}{1.6ex} \mkern-3mu A} =-{\rm tr}(\tilde{\!\boldsymbol{\mathsf{A}}}),\quad q_{\rule{0pt}{1.6ex} \mkern-3mu A} = \tfrac{1}{2}({\rm tr}(\tilde{\!\boldsymbol{\mathsf{A}}})^{2} - {\rm tr}(\tilde{\!\boldsymbol{\mathsf{A}}}^2)),\quad r_{\rule{0pt}{1.6ex} \mkern-4mu A} =-{\rm det}(\tilde{\!\boldsymbol{\mathsf{A}}}), \end{equation}
where 
${\rm tr}({\cdot })$ and 
${\rm det}( {\cdot } )$ represent the trace and determinant, respectively. For incompressible flows (
$p_{\rule {0pt}{1.6ex} \mkern -3mu A} = 0$), four classes of local streamline topologies are separated by degenerate geometries in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane. Using the invariants of 
$\tilde {\!\boldsymbol{\mathsf{A}}}$ is advantageous compared with using the invariants of 
$\boldsymbol{\mathsf{A}}$ since the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane is a bounded phase space and it provides a more complete representation of streamline geometries (Das & Girimaji Reference Das and Girimaji2019, Reference Das and Girimaji2020a,Reference Das and Girimajib). For example, the aspect ratio of purely elliptic local streamlines (
$q_{\rule {0pt}{1.6ex} \mkern -3mu A} > 0$ and 
$r_{\rule {0pt}{1.6ex} \mkern -4mu A} = 0$) is completely characterized by 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$, but not by 
$Q = A^2 q_{\rule {0pt}{1.6ex} \mkern -3mu A}$. However, while the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane efficiently characterizes local streamline geometries at critical points, additional parameters are required to fully describe all geometries (Das & Girimaji Reference Das and Girimaji2020a).
Following Das & Girimaji (Reference Das and Girimaji2020b), we consider the local streamline geometry in the context of the modes of deformation of a fluid parcel: extensional straining, (symmetric and antisymmetric) shearing and rigid rotation. The well-known Cauchy–Stokes decomposition of the VGT, 
$\tilde {\!\boldsymbol{\mathsf{A}}} = \tilde {\!\boldsymbol{\mathsf{S}}} + \tilde {\!\boldsymbol{\mathsf{W}}}$, disambiguates contributions from the symmetric strain rate tensor, 
$\tilde {\!\boldsymbol{\mathsf{S}}} = ( \tilde {\!\boldsymbol{\mathsf{A}}} + \tilde {\!\boldsymbol{\mathsf{A}}}^{\rm T})/2$, and the antisymmetric vorticity tensor, 
$\tilde {\!\boldsymbol{\mathsf{W}}} = ( \tilde {\!\boldsymbol{\mathsf{A}}} - \tilde {\!\boldsymbol{\mathsf{A}}}^{\rm T})/2$. It has enabled insightful characterizations of the VGT dynamics from the perspective of the strain rate eigenframe (Tom, Carbone & Bragg Reference Tom, Carbone and Bragg2021). However, it does not disambiguate symmetric shearing from extensional straining in 
$\tilde {\!\boldsymbol{\mathsf{S}}}$ or antisymmetric shearing from rigid rotation in 
$\tilde {\!\boldsymbol{\mathsf{W}}}$. This limitation motivated the development of the triple decomposition of the VGT (Kolář Reference Kolář2007), which disambiguates all three fundamental modes of deformation.
Kolář (Reference Kolář2004, Reference Kolář2007) originally formulated the triple decomposition of the VGT by identifying a ‘basic’ reference frame in which motions associated with elongation, rigid rotation and pure shearing can be isolated. However, identifying a basic reference frame requires a challenging pointwise optimization problem, the solution of which is typically approximated over a finite number of frames (Kolář Reference Kolář2007; Nagata et al. Reference Nagata, Watanabe, Nagata and da Silva2020). More recently, Gao & Liu (Reference Gao and Liu2018, Reference Gao and Liu2019) introduced a unique triple decomposition, based on a related vorticity tensor decomposition (Liu et al. Reference Liu, Gao, Tian and Dong2018; Gao et al. Reference Gao, Yu, Liu and Liu2019), that is more computationally practical than that of Kolář (Reference Kolář2004, Reference Kolář2007). This triple decomposition is formally performed in a local ‘principal’ coordinate system (
$x^*, y^*, z^*$), which is related to the global coordinate system (
$x, y, z$) by an orthogonal transformation. In this principal frame, denoted by 
$( {\cdot } )^*$, the triple decomposition is given in normalized form by
\begin{equation} \tilde{\!\boldsymbol{\mathsf{A}}}^* = \underbrace{\begin{bmatrix} \dot{\epsilon}_{x^*} & 0 & 0 \\ 0 & \dot{\epsilon}_{y^*} & 0 \\ 0 & 0 & \dot{\epsilon}_{z^*} \end{bmatrix}}_{\boldsymbol{\dot{\epsilon}}^*} + \underbrace{\begin{bmatrix} 0 & 0 & 0 \\ \dot{\gamma}_{z^*} & 0 & 0 \\ \dot{\gamma}_{y^*} & \dot{\gamma}_{x^*} & 0 \end{bmatrix}}_{\boldsymbol{\dot{\gamma}}^*} + \underbrace{\begin{bmatrix} 0 & -\dot{\varphi_{z^*}} & 0 \\ \dot{\varphi_{z^*}} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}}_{\boldsymbol{\dot{\varphi}}^*}. \end{equation}
Here, 
$\boldsymbol {\dot {\epsilon }}^*$, 
$\boldsymbol {\dot {\gamma }}^*$ and 
$\boldsymbol {\dot {\varphi }}^*$ represent the normal straining, pure shearing and rigid body rotation tensors, respectively. Their constituents can be directly identified from the components of the VGT in the principal frame (Gao & Liu Reference Gao and Liu2018, Reference Gao and Liu2019; Das & Girimaji Reference Das and Girimaji2020b). Their representations in the global coordinates (
$\boldsymbol {\dot {\epsilon }}$, 
$\boldsymbol {\dot {\gamma }}$ and 
$\boldsymbol {\dot {\varphi }}$) can subsequently be recovered by inverting (i.e. transposing) the original orthogonal transformation (Gao et al. Reference Gao, Yu, Liu and Liu2019).
The components of the normal straining tensor represent the real parts of the eigenvalues of 
$\tilde {\!\boldsymbol{\mathsf{A}}}$, which are identical to those of 
$\tilde {\!\boldsymbol{\mathsf{A}}}^*$. For points with rotational local streamlines, 
$\tilde {\!\boldsymbol{\mathsf{A}}}$ has a pair of complex eigenvalues and the real eigenvector defines the local rotation axis. In this case, the transformation to the principal frame is identified by (i) using a real Schur decomposition to align the 
$z^*$-axis with the real eigenvector of the VGT and (ii) orienting the 
$x^*$–
$y^*$ plane to minimize the local rotational speed (Liu et al. Reference Liu, Gao, Tian and Dong2018; Das & Girimaji Reference Das and Girimaji2020b). One advantage of (1.3) is that it provides representations of the strength (
$2 \dot {\varphi }_{z^*}$) and the axis (
$z^*$) of rigid rotation that are Galilean invariant (Wang, Gao & Liu Reference Wang, Gao and Liu2018). Unlike the rotational case, the VGT has only real eigenvalues when the local streamline geometry is non-rotational (
$\boldsymbol {\dot {\varphi }}^* = \boldsymbol {0}$). In this case, the principal frame is identified by using a Schur decomposition to transform the VGT into a triangular tensor. The modes of deformation are then isolated by decomposing this transformed tensor into a normal, diagonal tensor representing normal straining and a non-normal, strictly triangular tensor representing pure shearing (Keylock Reference Keylock2018; Das & Girimaji Reference Das and Girimaji2020b).
The triple decomposition enables refined analyses of the influences of fundamental constituents of the VGT. For example, the original triple decomposition (Kolář Reference Kolář2007) has been used to show that lifetimes of fundamental flow structures at macroscopic scales (where viscosity can be neglected) can be related to stability of rigid rotation, linear instability of pure shearing and exponential instability of irrotational straining (Hoffman Reference Hoffman2021). At small scales, the more recent triple decomposition (Gao & Liu Reference Gao and Liu2018, Reference Gao and Liu2019) has been used to show that pure shearing is typically the dominant contributor to energy dissipation (Wu et al. Reference Wu, Zhang, Wang, Zou and Chen2020) and intermittency (Das & Girimaji Reference Das and Girimaji2020b) in turbulent flows. Further, the symmetric and antisymmetric components of 
$\boldsymbol {\dot {\gamma }}^*$ are given by 
$\boldsymbol {\dot {\gamma }}^*_{\rule {0pt}{-2ex} \mkern -2mu S} = ( \boldsymbol {\dot {\gamma }}^* + \boldsymbol {\dot {\gamma }}^*{}^{\rm T} )/2$ and 
$\boldsymbol {\dot {\gamma }}^*_{\rule {0pt}{-2ex} \mkern -2mu W} = ( \boldsymbol {\dot {\gamma }}^* - \boldsymbol {\dot {\gamma }}^*{}^{\rm T} )/2$, respectively. In this manner, the triple decomposition is more refined than the Cauchy–Stokes decomposition since 
$\tilde {\!\boldsymbol{\mathsf{S}}}^* = \boldsymbol {\dot {\epsilon }}^* + \boldsymbol {\dot {\gamma }}^*_{\rule {0pt}{-2ex} \mkern -2mu S}$ and 
${\tilde {\!\boldsymbol{\mathsf{W}}}^* = \boldsymbol {\dot {\varphi }}^* + \boldsymbol {\dot {\gamma }}^*_{\rule {0pt}{-2ex} \mkern -2mu W}}$ (Gao & Liu Reference Gao and Liu2019; Das & Girimaji Reference Das and Girimaji2020b). As described in detail by Das & Girimaji (Reference Das and Girimaji2020b), the triple decomposition also enables a natural characterization of local streamline topologies and geometries. Similar topological analyses of vortical flow features have also been proposed (Nakayama Reference Nakayama2017), but we focus on the triple decomposition for the advantages outlined herein.
The ability of the triple decomposition to capture local streamline structure in terms of fundamental modes of deformation has guided efforts to define improved vortex criteria. There are an abundance of criteria to identify vortices that are based on various features (e.g. eigenvalues) of the VGT and that adopt various philosophies of what constitutes a vortex (Chakraborty, Balachandar & Adrian Reference Chakraborty, Balachandar and Adrian2005; Epps Reference Epps2017; Günther & Theisel Reference Günther and Theisel2018; Liu et al. Reference Liu, Gao, Dong, Wang, Liu, Zhang, Cai and Gui2019a; Haller Reference Haller2021). Debates surrounding these criteria primarily involve their (i) philosophical underpinnings, (ii) threshold sensitivities and (iii) observational invariances.
Regarding (i), the Cauchy–Stokes decomposition underlies many common symmetry- based vortex criteria, including the 
$Q$ (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988) and 
$\lambda _2$ (Jeong & Hussain Reference Jeong and Hussain1995) criteria. Local streamline topology underlies many common geometry-based vortex criteria, including the 
$\varDelta$ (Chong et al. Reference Chong, Perry and Cantwell1990) and 
$\lambda _{ci}$ (Chakraborty et al. Reference Chakraborty, Balachandar and Adrian2005) criteria. Like the geometry-based methods, and unlike the symmetry-based methods, the rigid vorticity criterion (
$\dot {\varphi }_{z^*} > 0$) (Tian et al. Reference Tian, Gao, Dong and Liu2018) captures all rotational local streamline geometries under the assumption that rigid rotation is an essential ingredient of a vortex (Liu et al. Reference Liu, Gao, Dong, Wang, Liu, Zhang, Cai and Gui2019a; Das & Girimaji Reference Das and Girimaji2020b). The philosophical distinction between symmetry-based and geometry-based criteria also underlies the so-called ‘disappearing vortex problem’ in which, fixing the VGT configuration and strain rate, increasing only the vorticity magnitude can remove a geometry-based vortex from the flow (Chakraborty et al. Reference Chakraborty, Balachandar and Adrian2005; Kolář & Šístek Reference Kolář and Šístek2020, Reference Kolář and Šístek2022). However, we here adopt the geometry-based viewpoint since, unlike vorticity, rigid rotation persistently underlies rotational local streamline topologies in all inertial frames. This interpretation in terms of local streamline topology has the potential to elucidate connections to related (e.g. elliptic) instabilities.
Regarding (ii), the Omega (
$\varOmega$) class of vortex criteria (Liu et al. Reference Liu, Wang, Yang and Duan2016; Dong et al. Reference Dong, Wang, Chen, Dong, Zhang and Liu2018; Dong, Gao & Liu Reference Dong, Gao and Liu2019; Liu & Liu Reference Liu and Liu2019) is advantageous since it uses quantities that are bounded and less threshold sensitive than the aforementioned methods. Regarding (iii), whereas most common vortex criteria are Galilean invariant, they are typically not objective since they are not preserved in rotating reference frames (Epps Reference Epps2017; Günther & Theisel Reference Günther and Theisel2018). However, the objectivized 
$\dot {\varphi }_{z^*}$ (Liu, Gao & Liu Reference Liu, Gao and Liu2019b) and objectivized 
$\varOmega$ (Liu et al. Reference Liu, Gao, Wang and Liu2019c) criteria, which are formulated by replacing 
$\tilde {\!\boldsymbol{\mathsf{W}}}$ with its deviation from its global spatial mean, remain invariant in these reference frames. Moreover, they are among the only compatible (i.e. self-consistent) objectivized vortex criteria out of the modifications commonly associated with the vortex criteria we have discussed (Haller Reference Haller2021). This advantage enhances the experimental verifiability and clarifies the physical significance of visualizations of the corresponding vortex structures.
Synthesizing the advantages of the geometry-based vortex definitions and the 
$\varOmega$ class of vortex criteria, we identify vortices using the 
$\varOmega _r$ method in the present investigation. This criterion is formulated in terms of the quantity
\begin{equation} \varOmega_r = \frac{(\boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{e}_{z^*})^2}{2 (\boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{e}_{z^*})^2 - 4 \lambda_{ci}{\vphantom{\left( \boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{e}_{z^*} \right)}}^{\mkern-14mu 2} + 4\varepsilon_{vort}}, \end{equation}
where 
$\lambda _{ci}$ is the imaginary part of the complex eigenvalues of 
$\boldsymbol{\mathsf{A}}$, 
$\boldsymbol {e}_{z^*}$ is the unit vector along the 
$z^*$-axis and 
$\varepsilon _{vort}$ is a numerical threshold used to prevent division by zero. Vortices are theoretically identified as spatially connected regions satisfying 
$\varOmega _r > 0.5$ when 
$\varepsilon _{vort} = 0$. In practice, however, vortices are identified using a small 
$\varepsilon _{vort} > 0$ and 
$\varOmega _r \geq 0.52$ (Liu et al. Reference Liu, Gao, Dong, Wang, Liu, Zhang, Cai and Gui2019a) to e.g. remove weak vortices.
1.3. Contributions
In this paper, we utilize the advantageous properties of geometry-based analyses of the VGT to efficiently characterize turbulence initiated by a vortex ring collision. We use the adaptive, multiresolution computational techniques discussed in § 2 to perform a direct numerical simulation of this flow at 
$Re_{\varGamma _0} = 4000$. In § 3, we establish the fidelity of our simulation and we visualize and discuss the various regimes of its evolution. In § 4, we analyse the partitioning of the velocity gradients to characterize these regimes in terms of the modes of deformation. In § 5, we introduce a geometry-based phase space that characterizes the action of the elliptic instability and its interplay with other mechanisms driving the turbulent flow. Our analyses reveal statistical features of the VGT that are similar to those of previous simulations. They also provide tools with the potential to help disentangle mechanisms underlying vortex interactions during transition and turbulent decay. Finally, we summarize our results in the context of previous works and highlight promising future research prospects in § 6.
2. Methods
2.1. Computational method
To efficiently simulate a turbulent vortex ring collision, we adopt a recently developed multiresolution solver for viscous, incompressible flows on unbounded domains (Liska & Colonius Reference Liska and Colonius2016; Dorschner et al. Reference Dorschner, Yu, Mengaldo and Colonius2020; Yu Reference Yu2021; Yu, Dorschner & Colonius Reference Yu, Dorschner and Colonius2022). Yu et al. (Reference Yu, Dorschner and Colonius2022) provide a detailed discussion of the formulation, properties and performance of the method. We summarize the key advantages of the solver here and expound the computational formulation in Appendix A. The advantages we discuss allow us to simulate a relatively high Reynolds number vortex ring collision at a relatively low computational cost.
The Navier–Stokes equations (NSE) are spatially discretized onto a staggered Cartesian grid using a second-order-accurate finite-volume scheme that endows discrete operators with useful properties (i) (Liska & Colonius Reference Liska and Colonius2016). Discrete differential operators are constructed to mimic the symmetry, orthogonality and integration properties of their continuous counterparts. They also commute with the Laplacian and integrating factor operators, as defined in Appendix A. Furthermore, the discretization of the nonlinear term in the momentum equations preserves relevant (e.g. energy) conservation properties in the absence of viscosity. Together, the mimesis, commutativity and conservation properties of the discretization scheme facilitate fast, stable, high-fidelity simulations of turbulent flows.
The computational methods we employ also have high parallel efficiency (ii) and linear algorithmic complexity (iii). The computational efficiency of the flow solver is primarily centred around solving the discrete pressure Poisson equation on a formally unbounded grid using the lattice Green's function (LGF) (Liska & Colonius Reference Liska and Colonius2014, Reference Liska and Colonius2016; Liska Reference Liska2016). The Poisson equation is obtained by taking the divergence of the NSE in rotational form, such that the source term is 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {r}$, where 
$\boldsymbol {r} = \boldsymbol {u} \times \boldsymbol {\omega }$ is the Lamb vector. By considering flows with at least exponentially decaying far-field vorticity, the approximate support of this source field can be captured using a finite computational domain. Given a source cutoff threshold, the finite domain is adaptively truncated to capture only the regions relevant to the Poisson problem. Solving the Poisson problem over this domain involves the convolution of the LGF with the source field. The flow solver achieves (ii) and (iii) by efficiently evaluating this convolution via a fast multipole method (Liska & Colonius Reference Liska and Colonius2014) that compresses the kernel using polynomial interpolation. This method is accelerated by exploiting the efficiency of fast Fourier transforms on a block-structured Cartesian grid.
In addition to spatially adapting the extent of the computational domain, adaptive multiresolution discretization (iv) is achieved by using adaptive mesh refinement (AMR) to reduce the number of degrees of freedom required for solutions. As discussed previously (Dorschner et al. Reference Dorschner, Yu, Mengaldo and Colonius2020; Yu et al. Reference Yu, Dorschner and Colonius2022), the present AMR framework is carefully constructed to preserve the desirable operator properties (i) and augment the efficiencies (ii, iii) associated with the uniform-grid framework (Liska & Colonius Reference Liska and Colonius2014, Reference Liska and Colonius2016). In the AMR framework, the computational grid is partitioned into multiple levels, each with double the resolution in each direction as the previous level. The spatial regions associated with each level are non-overlapping, except for extended regions that are used to compute a combined source term that includes a correction induced by the difference between the coarse-grid and fine-grid partial solutions. As formulated in Appendix A, a region is refined when its combined source exceeds a threshold and it is coarsened when its combined source falls below a smaller threshold. As shown in § 3.1 (see table 1), this AMR formulation drastically reduces the number of computational cells required to capture a head-on vortex ring collision compared with a fixed-resolution scheme.
Table 1. Reference times used to analyse the initial (red), transitional (green) and turbulent (blue) regimes of the present vortex ring collision, where 
$t^* = 14.77$ is the time of maximum dissipation. Here, 
$\bar {R}_p$ represents the mean vortex ring radius (see Appendix B), 
$V$ is the volume of the computational domain and 
$\phi _k$ is the fraction of 
$V$ occupied by level 
$k$ of the AMR grid, which has grid spacing 
$\Delta x_k = \Delta x_{base}/2^{k}$ with 
$\Delta x_{base} = 0.04$.
2.2. Vortex ring collision simulation
As depicted in figure 1, we consider a flow configuration in which the vortex rings are initialized with opposing circulations such that they propagate toward one another along the 
$z$-axis and meet at the collision plane at 
$z = 0$. The rings are initialized a distance 
$L_z = 2.5 R_0$ apart, which is sufficiently large to mitigate their mutual influence during the most vigorous period of equilibration. Both rings are initialized with Gaussian vorticity distributions (1.1) such that 
$Re_{\varGamma _0} = 4000$ and 
$\delta _0 = 0.2$. Unless otherwise stated, we use the initial circulation, 
$\varGamma _0 = 1$, and radius, 
$R_0 = 1$, of each ring to non-dimensionalize all variables. To excite transition, we randomly perturb the radii of the vortex rings using the first 32 Fourier modes in 
$\theta$, which are prescribed random phases and uniform magnitudes, 
$R_{pert} = 5 \times 10^{-4}$. Consistent with previous tests (Yu et al. Reference Yu, Dorschner and Colonius2022), these initial perturbations are sufficiently large to dominate perturbations incurred by discretization errors.
Figure 1. Initial geometry of the flow configuration used to simulate the head-on collision between vortex rings. The shading of the vortex cores reflects their Gaussian vorticity profiles.
The computational mesh we use has 
$N_{level} = 2$ levels of refinement beyond the base level such that the ratio of the coarsest-grid spacing to the finest-grid spacing is 
$\Delta x_{base} / \Delta x_{fine} = 4$. Based on preliminary simulations of turbulent vortex rings (Liska & Colonius Reference Liska and Colonius2016) and vortex ring collisions (Yu et al. Reference Yu, Dorschner and Colonius2022), we select 
$a_0/\Delta x_{base} = 5$ and 
$\Delta t/ \Delta x_{fine} = 0.35$ to ensure the flow is well resolved throughout the simulation. Finally, parameters controlling the spatial and mesh refinement thresholds are chosen as discussed in Appendix A.
2.3. Simulation integral metrics
We track the evolution and fidelity of the simulation using integral metrics associated with incompressible flows (Liska & Colonius Reference Liska and Colonius2016). Particularly, we compute the hydrodynamic impulse, helicity, vortical kinetic energy and enstrophy of the flow, which are denoted by 
$\boldsymbol {I}_{V}$, 
$H$, 
$K_V$ and 
$E$, respectively. These integrals are formally evaluated on an unbounded domain, but we evaluate them using the finite AMR grid as
\begin{equation} \left.\begin{gathered} \boldsymbol{I}_{V}(t) =\int_{V(t)} (\boldsymbol{x} \times \boldsymbol{\omega})\,{\rm d} V,\quad H(t) = \int_{V(t)} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\omega} )\,{\rm d} V, \\ K_{V}(t) = \int_{V(t)} \boldsymbol{u} \boldsymbol{\cdot} (\boldsymbol{x} \times \boldsymbol{\omega})\,{\rm d} V,\quad E(t) = \tfrac{1}{2} \int_{V(t)} \lvert \boldsymbol{\omega} \rvert^2 \,{\rm d} V, \end{gathered}\right\} \end{equation}
where 
$V(t)$ is the time-varying AMR grid. The impulse is the appropriate measure of momentum since it converges for flows on unbounded domains with compact vorticity. The vortical kinetic energy can also be expressed as 
$K_V = K + K_{\partial V}$, where 
$K$ represents the kinetic energy and 
$K_{\partial V}$ is a correction term based on the flow at the boundary of the grid, 
$\partial V(t)$. These metrics can be expressed as
\begin{equation} K(t) = \tfrac{1}{2} \int_{V(t)} \lvert \boldsymbol{u} \rvert^2\,{\rm d} V,\quad K_{\partial V}(t) = \int_{\partial V(t)} \boldsymbol{x} \boldsymbol{\cdot} \big((\boldsymbol{u} \boldsymbol{u} ) \boldsymbol{\cdot} \boldsymbol{n} - \tfrac{1}{2} \lvert \boldsymbol{u} \rvert^2 \boldsymbol{n}\big)\,{\rm d} S, \end{equation}
where 
$\boldsymbol {n}$ is the normal vector of 
$\partial V$ (Wu, Ma & Zhou Reference Wu, Ma and Zhou2015). For vanishing far-field velocity, 
$K_{\partial V}$ vanishes on unbounded domains and, for the present grid, we make use of the smallness of 
$K_{\partial V}$ when analysing dissipation.
In the absence of non-conservative external body forces, the hydrodynamic impulse is conserved for incompressible flows on unbounded domains (Saffman Reference Saffman1993). The helicity would also be conserved in the absence of viscosity, and it is useful for assessing simulation fidelity as the vortex rings initially approach the collision plane since the evolution of the flow is dominated by inviscid effects. These integral metrics initially evaluate to 
$\boldsymbol{I}_{V}(0) = \boldsymbol{0}$ and 
$H(0) = 0$ due to the spatial symmetries of the initial flow configuration. These initial symmetries hold to the extent that the vorticity is well captured and the contributions of the random perturbations used to excite instability growth are negligible. As the flow evolves, subsequent deviations from these initial values reflect the degree to which the corresponding initial symmetries are broken.
The enstrophy and the kinetic energy provide a more detailed picture of simulation fidelity during transition and turbulent decay, when viscous dissipation at small scales becomes relevant. For unsteady, incompressible flows on unbounded domains, the dissipation governs the decay rate of kinetic energy and can be expressed in terms of the enstrophy. Comparing these integrals is useful for characterizing (i) the degree to which small-scale features are resolved during peak dissipation and (ii) the flux of kinetic energy out of the finite computational domain (Archer et al. Reference Archer, Thomas and Coleman2008). We therefore introduce effective Reynolds numbers, which are given by
\begin{equation} \frac{Re^{eff}_S (t)}{Re_{\varGamma_0}} =-\frac{\varPhi_S (t)}{{\rm d} K/{\rm d} t},\quad \frac{Re^{eff}_W (t)}{Re_{\varGamma_0}} =-\frac{\varPhi_W (t)}{{\rm d} K/{\rm d} t}, \end{equation}
where 
$\varPhi _S$ is the volume-integrated dissipation and 
$\varPhi _W = 2 E/Re_{\varGamma _0}$ is its enstrophy-based counterpart (Serrin Reference Serrin1959). Here, we differentiate 
$K$ instead of 
$K_V$ when computing the effective Reynolds numbers to prevent amplification of the noise associated with adaptations in the computational domain, to which 
$K_V$ is more sensitive. This is justified since 
$K$ and 
$K_V$ are nearly identical throughout the present simulations (see figure 2). The ratio 
$Re^{eff}_S$ is useful for assessing spatial resolution since the dissipation can vary significantly during transition and turbulent decay. The corresponding Kolmogorov scale, 
$\eta = (\nu ^3/\varPhi _S)^{1/4}$, can also be used to validate the selected grid spacings. The difference between 
$Re^{eff}_S$ and 
$Re^{eff}_W$ reflects the relative significance of the acceleration of the flow on 
$\partial V$ through the boundary integral in the Bobyleff–Forsyth formula (Serrin Reference Serrin1959). Together, the error metrics defined in this section comprehensively characterize the fidelity of the simulation as its flow structures evolve and the computational domain adapts accordingly.
Figure 2. Temporal evolution of the integral metrics defined in § 2.3 over the course of the simulation. The vertical lines correspond to the reference times in table 1 and they are coloured accordingly. The horizontal lines in the 
$Re^{eff}$ panel represent 
$Re_{\varGamma _0}$ (solid) with 10 % margins (dashed). The impulse magnitude is normalized by that of each vortex ring in isolation, 
$\lvert \boldsymbol {I}_{V1} \rvert \approx 1.02{\rm \pi} \approx 3.204$. The enstrophies 
$E_{\mathcal {E}}$ and 
$E_{\mathcal {C}}$ are computed using vorticities located at the edges and centres, respectively, of the computational cells.
3. Evolution of integral metrics and vortical structures
3.1. Evolution of integral metrics
Figure 2 shows the evolution of the integral metrics from § 2.3 over the course of the simulation. In the subsequent analysis, we reference the various regimes of flow development with respect to the time, 
$t^* = 14.77$, at which maximum dissipation is attained. Table 1 qualitatively characterizes the state of the simulation at each reference time we consider for the initial, transitional and turbulent regimes of the simulation.
The initial evolution of the flow involves a rapid period of equilibration (
$t \lesssim 0.25 t^*$) and the propagation of the equilibrated rings towards the collision plane (
$0.25 t^* \lesssim t \lesssim 0.50 t^*$). The interaction of the rings accelerates their radial expansion (
$0.50 t^* \lesssim t \lesssim 0.75 t^*$) and the elliptic instability eventually emerges along the expanding rings (
$0.75 t^* \lesssim t \lesssim 0.90 t^*$). Appendix B supports the importance of the elliptic instability during the early stages of transition. Subsequently, the flow transitions to turbulence (
$0.90 t^* \lesssim t \lesssim t^*$) and rapidly produces small-scale flow structures. Following transition, the flow undergoes turbulent decay for the remainder of the simulation (i.e. for 
$t \gtrsim t^*$). See § 3.2 for visualizations of the flow at the reference times from table 1 associated with each of these regimes of evolution.
As the vortex rings initially propagate towards the collision plane (
$t \lesssim 0.50 t^*$), the kinetic energy decays slowly and the enstrophy and dissipation are relatively small. The effective Reynolds numbers rapidly adjust to the value of 
$Re_{\varGamma _0}$ during the initial equilibration period (
$t \lesssim 0.25 t^*$) and remain roughly constant as the equilibrated rings approach the collision plane (
$0.25 t^* \lesssim t \lesssim 0.50 t^*$). The helicity is well conserved in this regime since the flow evolves in a nearly inviscid fashion. Further, the impulse is initially small and grows relatively slowly during this period. These results suggest that the symmetries associated with the handedness and momentum distribution of the flow are well preserved in the initial regime of evolution.
As the rings expand radially at the collision plane (
$0.50 t^* \lesssim t \lesssim 0.75 t^*$) and the elliptic instability emerges (
$0.75 t^* \lesssim t \lesssim 0.90 t^*$), the kinetic energy decays more rapidly and the dissipation grows. During these periods, the helicity symmetry remains well preserved and the effective Reynolds numbers remain relatively constant near 
$Re_{\varGamma _0}$, suggesting that the flow is well resolved. However, the impulse magnitude varies more rapidly in time due to the rapid radial expansion of the rings. In following the expanding vortical flow at the collision plane, the adaptations of the domain break the symmetry associated with impulse integral more significantly than during the initial evolution of the rings. The resulting growth in 
$\lvert \boldsymbol {I}_V \rvert$ is primarily attributed to its component in the 
$z$ direction, along which the domain is compressed as the flow concentrates about the collision plane.
As the flow transitions to turbulence (
$0.90 t^* \lesssim t \lesssim t^*$), the kinetic energy decays even more rapidly and the dissipation approaches its maximum value. Due to the proliferation of small-scale flow structures during this period, the effective Reynolds numbers drop to their minimum values at 
$t \approx t^*$, when the flow is most difficult to resolve. The increased difference between 
$Re^{eff}_S$ and 
$Re^{eff}_W$ reflects that the acceleration of the flow near 
$\partial V$ is more relevant at this time. The rapid generation of small-scale flow structures also implies that viscosity plays a more important role in this regime. Correspondingly, the helicity begins to vary in time in this regime, reaching its maximum rate of change at the time of peak dissipation. Its variations reflect that vortex lines in the flow undergo rapid topological changes during transition. By contrast, the impulse magnitude decays to a roughly constant value as the radial expansion of the rings slows in this transitional regime.
During the turbulent decay of the flow (
$t \gtrsim t^*$), the kinetic energy becomes small and the dissipation decays rapidly, eventually falling below its initial value (at 
$t \approx 1.62 t^*$). The dissipation matches the kinetic energy decay rate more closely for this regime than for transition. Further, as the turbulence develops, the effective Reynolds numbers agree well with one another and, to a lesser extent, with 
$Re_{\varGamma _0}$. These features reflect, respectively, that the acceleration near 
$\partial V$ is less significant and that the small scales are relatively well resolved, especially with respect to the transitional period. The helicity variations in this turbulent regime also eventually slow relative to those observed during transition. Similarly, the impulse remains roughly constant around its value at 
$t = t^*$. Whereas the 
$z$ component dominates the impulse magnitude during the radial expansion of the rings, all impulse components have similar magnitudes in this turbulent regime.
The evolution of the integral metrics characterizes the various regimes of the flow and supports the fidelity of our simulation. For example, the helicity is well conserved during the nearly inviscid evolution of the vortex rings and its subsequent variations are relatively small in magnitude. Further, the variations in the impulse magnitude throughout the simulation remain less than 5 % of the impulse associated with each vortex ring in isolation. These results indicate that the symmetries associated with the handedness and momentum distribution of the flow also remain well preserved in the appropriate regimes of the simulation.
Our dissipation analysis also suggests that the small-scale flow structures remain reasonably well resolved throughout the entire simulation. After equilibration, the maximum relative deviation in 
$Re^{eff}_S$ from 
$Re_{\varGamma _0}$ is roughly 6.5 % and it occurs around 
$t \approx t^*$. This relative deviation is similar to that of a previous simulation of a single vortex ring at 
$Re_{\varGamma _0} = 7500$ using a finite computational grid (Archer et al. Reference Archer, Thomas and Coleman2008). Moreover, it is considerably smaller during the approach and radial expansion of the rings and, to a lesser extent, during turbulent decay. Even during peak dissipation, when the Kolmogorov scale is smallest, the finest grid has acceptable resolution since 
$\Delta x_{fine} / \eta \approx 3.42$. Altogether, these results suggest that our simulation is well resolved and support our analysis of the mechanisms underlying transition and turbulent decay.
3.2. Evolution of vortical flow structures
For the present simulation, we identify vortices using the 
$\varOmega _r$ criterion (Liu et al. Reference Liu, Gao, Dong, Wang, Liu, Zhang, Cai and Gui2019a) with a numerical threshold of 
$\varepsilon _{vort}$ = 0.04. This criterion provides connections to the triple decomposition of the VGT and the structure of local streamlines. Due to the well-preserved symmetries of the flow, the global spatial mean of the vorticity tensor is nearly zero and, hence, the 
$\varOmega _r$ criterion is nearly objective (Liu et al. Reference Liu, Gao, Wang and Liu2019c) for the present simulation. We specifically visualize the flow using 
$\varOmega _r = 0.52$ and 
$\varOmega _r = 0.93$ to investigate the structures of the vortex boundaries and the vortex cores, respectively. We colour these structures using 
$\cos \theta ^*$, where 
$\theta ^*$ is the angle between the 
$z^*$-axis and the 
$z$-axis. This colour scheme enables the identification of antiparallel vortices along the 
$z$-axis, which play an important role in mediating transition and generating small-scale flow structures in the present vortex ring collision. In figure 3, we visualize the vortical structures in the flow at each reference time from table 1.
Figure 3. Visualizations of the vortex boundaries (
$\varOmega _r = 0.52$, left side) and vortex cores (
$\varOmega _r = 0.93$, right side), coloured by 
$\cos \theta ^*$, for each reference time from table 1. A movie depicting the evolution of the vortex boundaries from the auxiliary viewpoint (leftmost column) is provided as supplementary material available at https://doi.org/10.1017/jfm.2024.90.
During the initial evolution of the flow, the equilibrated vortex rings approach the collision plane and expand radially due to their mutual interaction. In this regime, the thinning of the vortex boundaries and cores illustrates the mechanisms driving the shift from a rigid-rotation-dominated regime to a shearing-dominated regime. Further, the visualizations at 
$t = 0.75 t^*$ depict the emergence of the short-wave elliptic instability, which is consistent with previous vortex ring collision simulations in similar parameter regimes (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021).
The transitional regime of the flow is marked by the development of secondary vortex filaments and the subsequent generation of small-scale vortical flow structures. At 
$t = 0.90 t^*$, the visualizations show the development of secondary vortical structures around the circumference of the collision. These structures consist of antiparallel vortex filament pairs that arise in regions where the elliptic instability drives local interactions between the rings. This behaviour supports the notion that the elliptic instability mediates the initial transition of the rings, leading to the development of secondary vortical structures. The antiparallel secondary filaments become increasingly densely packed as transition progresses and they mediate the proliferation of small-scale vortical flow structures, e.g. as observed at 
$t = t^*$.
In Appendix B, we decompose the flow into azimuthal Fourier modes to characterize the wavenumbers of the perturbations that dominate transition. Our analysis confirms that the short-wave elliptic instability, with wavelength of the order of the core radius, mediates the initial stages of transition. We further show that, at 
$t = t^*$, the most prominent antiparallel vortices occur at the second harmonic of an originally dominant perturbation. Taken together, our visualizations and perturbation analysis are consistent with the initial stages of the iterative elliptic instability pathway, which is driven by subsequent generations of antiparallel vortex filaments (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020).
During the turbulent decay of the flow, the geometric features of the vortex boundaries remain similar at each reference time. However, as energy is dissipated, the smallest-scale vortices are progressively destroyed and the vortical flow structures grow larger in time. The structures of the vortex cores and boundaries reinforce the importance of the interactions between the secondary vortex filaments in mediating the evolution of the turbulent flow. The vortex boundaries also show the formation and ejection of vortex rings from the turbulent cloud resulting from the collision. These ejections, which are a hallmark of the Crow instability, often occur in regions where antiparallel vortex filaments interact and are of similar size to those filaments. This observation provides further evidence of the interplay between the elliptic and Crow instabilities driven by interacting vortex filaments around the turbulent cloud (Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021; Ostilla-Mónico et al. Reference Ostilla-Mónico, McKeown, Brenner, Rubinstein and Pumir2021). In what follows, we develop machinery to probe these mechanisms in the context of features of the velocity gradients, with a particular emphasis on characterizing the action of the elliptic instability among other mechanisms.
4. Partitioning of velocity gradients
Here, we investigate the partitioning of the velocity gradients to characterize the evolution of the flow in the context of the fundamental modes of deformation. We first consider volumetric weighted averages of the relative contributions of various constituents of 
$\tilde {\!\boldsymbol{\mathsf{A}}}$ to the strength of the velocity gradients. These averages may be expressed as
\begin{equation} \langle {\xi}\rangle_{\mkern-3mu A^2} = \frac{\displaystyle\int_{V(t)} A^2 \xi \,{\rm d} V}{\displaystyle\int_{V(t)} A^2\,{\rm d} V}, \end{equation}
where 
$\xi \in \{\tilde {S}^2, \tilde {W}^2\}$ for the Cauchy–Stokes decomposition and 
$\xi \in \{\dot {\epsilon }^2, \dot {\gamma }^2, \dot {\varphi }^2, 2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} \}$ for the triple decomposition. Here, we have used that the Frobenius inner product, denoted by :, of a symmetric tensor with an antisymmetric tensor is zero and that 
$\dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu W} = \dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu S} = \dot {\gamma }^2/2$. The shear–rotation correlation term, 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} = 2 \dot {\varphi }_{z^*} \dot {\gamma }_{z^*} > 0$, reflects the presence of shearing in the plane of rigid rotation. All of the relative contributions discussed are unitarily invariant and, thus, they apply to both the principal coordinates and the global coordinates. Figure 4 shows how these relative contributions evolve during the simulation for both decompositions of 
$\tilde {\!\boldsymbol{\mathsf{A}}}$.
Figure 4. Relative contributions of the constituents of the Cauchy–Stokes decomposition and the triple decomposition to 
$A^2$. The vertical lines correspond to the reference times in table 1 and they are coloured accordingly.
Consistent with the equivalence of 
$\varPhi _S$ and 
$\varPhi _W$ for incompressible flows on unbounded domains (Serrin Reference Serrin1959), 
$\langle {\tilde {S}^2}\rangle _{\mkern -3mu A^2} \approx \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2} \approx 0.50$ for the present simulations. The largest deviations from this balance occur during equilibration (
$t \lesssim 0.25 t^*$) and around the time of peak dissipation (
$t \approx t^*$). These deviations are consistent with the behaviour of the effective Reynolds numbers in figure 2 and their smallness further validates the ability of the finite computational grid to approximate a formally unbounded flow. However, since 
$\langle {\tilde {S}^2}\rangle _{\mkern -3mu A^2}$ and 
$\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$ remain relatively constant throughout the simulation, they provide limited information about the nature of the velocity gradients as the flow progresses through its initial, transitional and turbulent regimes.
Compared with the constituents of the Cauchy–Stokes decomposition, the constituents of the triple decomposition show more pronounced variations associated with the different regimes of evolution. For the initial Gaussian vorticity profiles, the contribution of rigid rotation to the enstrophy dominates the contribution of antisymmetric shearing. During equilibration (
$t \lesssim 0.25 t^*$), the fluctuations in all contributions of the triple decomposition constituents reflect the redistribution of velocity gradients in the cores of the vortex rings. As the equilibrated rings approach the collision plane and spread (
$0.25 t^* \lesssim t \lesssim 0.75 t^*$), 
$\langle {\dot {\epsilon }^2}\rangle _{\mkern -3mu A^2}$ and 
$\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2}$ decrease and 
$\langle {\dot {\gamma }^2}\rangle _{\mkern -3mu A^2}$ and 
$2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ increase. As the elliptic instability emerges, these contributions level off in a regime where antisymmetric shearing dominates rigid rotation and shear–rotation correlations are enhanced. The subsequent development of the elliptic instability (
$0.75 t^* \lesssim t \lesssim 0.90t^*$) is marked by slight rebounds in the contributions of 
$\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2}$ and 
$\langle {\dot {\gamma }^2}\rangle _{\mkern -3mu A^2}$. These rebounds are associated with the emergence of secondary vortex filaments and, hence, the nonlinear evolution of the elliptic instability. The transition to turbulence (
$0.90 t^* \lesssim t \lesssim t^*$), which is associated with the generation of small scales and enhanced dissipation, is marked by a decrease in the contribution of 
$2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$.
Remarkably, even though the flow is not stationary during turbulent decay, the relative contributions of the constituents of the triple decomposition to the strength of the velocity gradients remain roughly constant after transition. Further, as summarized in table 2, these ‘equilibrium’ relative contributions are similar to those computed by Das & Girimaji (Reference Das and Girimaji2020b) for forced isotropic turbulence at high Taylor-scale Reynolds numbers. This agreement suggests that the velocity gradient partitioning may encode a relatively common balance in unbounded, incompressible turbulence with appropriate symmetries. In this balance, shearing makes the largest contribution to the velocity gradients and rigid rotation makes the smallest contribution.
Table 2. Comparison of the equilibrium partitioning of the velocity gradients for the present vortex ring collision with the partitioning computed for forced isotropic turbulence (Das & Girimaji Reference Das and Girimaji2020b). Here, the equilibrium partitioning is computed as the mean over the turbulent decay regime (
$t \gtrsim t^*$) and it is insensitive to the length of the averaging interval.
Beyond the strength of velocity gradients, it is also useful to examine the interplay between the modes of deformation in the context of vortical flow structures. Here, we introduce a new phase space defined by the relative contributions of 
$\langle {\dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu W}}\rangle _{\mkern -3mu A^2}$, 
$\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2}$ and (implicitly) 
$2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ to 
$\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$. The upper bound of this phase space is found by maximizing 
$\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2} / \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$, which occurs when 
$2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2} / \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2} = 0$. The bottom boundary is found by minimizing 
$\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2} / \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$, which occurs when 
$2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2} / \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2} = 2 \sqrt{\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2} \langle {\dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu W}}\rangle _{\mkern -3mu A^2}}/\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$. Correspondingly, these boundaries may be expressed as
\begin{align}
&\left(1-\sqrt{\frac{\langle
{\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu
W}}\rangle_{\mkern-3mu A^2}}{\langle
{\tilde{W}^2}\rangle_{\mkern-3mu A^2}}}\right)^2 \leq
\frac{\langle {\dot{\varphi}^2}\rangle_{\mkern-3mu
A^2}}{\langle {\tilde{W}^2}\rangle_{\mkern-3mu A^2}}\ \, \, \leq 1
- \frac{\langle {\dot{\gamma}^2_{\rule{0pt}{-2ex}
\mkern-2mu W}}\rangle_{\mkern-3mu A^2}}{\langle
{\tilde{W}^2}\rangle_{\mkern-3mu A^2}},
\end{align}
\begin{align} &\Updownarrow
\qquad\qquad \Updownarrow \nonumber\\ &-2 \left(\frac{\langle
{\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu
W}}\rangle_{\mkern-3mu A^2}}{\langle
{\tilde{W}^2}\rangle_{\mkern-3mu A^2}} -
\sqrt{\frac{\langle {\dot{\gamma}^2_{\rule{0pt}{-2ex}
\mkern-2mu W}}\rangle_{\mkern-3mu A^2}}{\langle
{\tilde{W}^2}\rangle_{\mkern-3mu A^2}}}\right) \geq
\frac{2 \langle {\boldsymbol{\dot{\varphi}} \mkern-5mu :
\mkern-5mu \boldsymbol{\dot{\gamma}}}\rangle_{\mkern-3mu
A^2}}{\langle {\tilde{W}^2}\rangle_{\mkern-3mu A^2}}\!\geq 0,
\end{align}
highlighting that lower and upper bounds of the rigid rotation contribution correspond to the upper and lower bounds, respectively, of the shear–rotation correlation contribution. The maximum value of 
$2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}/\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$ varies along the lower boundary of this phase space to ensure that the relative contributions sum to unity. The global maximum occurs when 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} / \tilde {W}^2 = 0.5$ or, equivalently, 
$\dot {\varphi }^2 / \tilde {W}^2 = \dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu W} / \tilde {W}^2 = 0.25$ at all points, which corresponds to the (pointwise) maximum of 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} = ( \sqrt {2} + 1 )^{-1}$ reported by Das & Girimaji (Reference Das and Girimaji2020b). This maximum corresponds to local streamlines in the principal frame for which shearing occurs exclusively in the plane of rigid rotation.
Figure 5 depicts the trajectory of the flow in this shear–rotation phase space and elucidates how the relationships between the constituents of enstrophy associated with the triple decomposition evolve in time. Following rapid variations during equilibration, the trajectory returns to a position in phase space similar to that of the initial condition at 
$t \approx 0.25 t^*$. The trajectory then undergoes a shift across the phase space as the equilibrated rings approach the collision plane and spread radially (
$0.25 t^* \lesssim t \lesssim 0.75 t^*$). This shift from a rigid-rotation-dominated regime to a shearing-dominated regime is associated with an enhanced contribution from 
$2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}/\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$. The development of the elliptic instability (
$0.75 t^* \lesssim t \lesssim 0.90 t^*$) and transition (
$0.90 t^* \lesssim t \lesssim t^*$) are associated with a shift in the direction of the trajectory towards smaller contributions of 
$2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}/\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$. During turbulent decay (
$t \gtrsim t^*$), the trajectory remains roughly fixed in the phase space and very close to its position at 
$t = t^*$.
Figure 5. Shear–rotation phase space trajectory of the flow, highlighting the evolution during equilibration (left inset) and transition and turbulent decay (right inset). The contours represent increments of 0.05 in the colour scale. The white circle marks the initial condition and the fill colours of the black circles correspond to the reference times they represent from table 1.
Considering this phase space trajectory in the context of the dissipation (see figure 2) reveals that, while the initial growth in dissipation is associated with enhanced shear–rotation correlations, its subsequent enhancement during transition is associated with a reduction in shear–rotation correlations. In Appendix C, we reexamine the visualizations from figure 3 in terms of the shear–rotation correlations in the flow to highlight their relationship to the vortical flow structures. In § 5, we interpret the effects of these shear–rotation correlations using a new, related phase space, based on local streamline geometry, to characterize the elliptic instability and other mechanisms.
5. Statistical geometry of local streamlines
5.1. Phase space transformations
The elliptic instability is associated with the resonance of the vortical flow with the underlying strain field and acts to break up elliptic streamlines. Consistent with this picture, we introduce a new geometry-based phase space that captures local flow features that (i) are conducive to the elliptic instability and (ii) characterize its action.
To address (i), we consider the angle, 
$\theta _\omega$, between the vorticity vectors associated with antisymmetric shearing and rigid rotation, which is given by
\begin{equation} \theta_\omega = \cos^{-1}\left(\sqrt{\frac{(\boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}})^2}{\dot{\varphi}^2 \dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}}}\right)\!. \end{equation}
Our focus on shear straining is consistent with the classical models of strained vortices used to characterize the elliptic instability (Kerswell Reference Kerswell2002). Since decreasing 
$\theta _\omega$ corresponds to increasing the alignment between shearing and rigid rotation, it can be associated with conditions conducive to the elliptic instability.
To address (ii), we consider the aspect ratio, 
$\zeta$, of the elliptic component of rotational local streamlines in the plane of rigid rotation, which is given by
\begin{equation} \zeta = \sqrt{1-e^2} = \sqrt{\frac{\dot{\varphi}^2}{\dot{\varphi}^2 + 2 \boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}}}}, \end{equation}
where 
$e$ represents the eccentricity of an ellipse with aspect ratio 
$\zeta$. This aspect ratio characterizes the scale-invariant geometry of the local streamlines in the plane of rigid rotation. As such, it can be used alongside 
$\theta _\omega$ to characterize the action of the elliptic instability by identifying how alignment between shearing and rigid rotation affects local streamline geometry.
The present shear–rotation and geometry-based phase spaces can be understood through nonlinear transformations of the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space. The transformations we derive express the relative contributions of the triple decomposition to the velocity gradients using 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$ and 
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$, and they represent the inverse transformations to those presented by Das & Girimaji (Reference Das and Girimaji2020b). However, an additional parameter, 
$\theta _\omega$, is generally required to evaluate our transformations. This extra parameter demonstrates that the invariants of 
$\tilde {\!\boldsymbol{\mathsf{A}}}$ alone cannot generally characterize the relative contributions of the constituents of the triple decomposition to the velocity gradients.
For rotational local streamlines, the transformations are given by
\begin{align} \left.\begin{gathered} \dot{\epsilon}^2 = \frac{(-3 \times 2^{2/3} q_{\rule{0pt}{1.6ex} \mkern-3mu A} + 3^{2/3} (2 \sqrt{3 \varDelta} + 9 \lvert r_{\rule{0pt}{1.6ex} \mkern-4mu A} \rvert)^{2/3})^2}{6^{5/3} (2 \sqrt{3 \varDelta} + 9 \lvert r_{\rule{0pt}{1.6ex} \mkern-4mu A} \rvert)^{2/3}},\quad \dot{\gamma}^2 = 2 \dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W} = 1 - 2 \dot{\epsilon}^2 - 2 q_{\rule{0pt}{1.6ex} \mkern-3mu A}, \\ \dot{\varphi}^2 = \left(-\sqrt{\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}} \cos \theta_{\omega} + \sqrt{1 - \dot{\epsilon}^2 - \dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W} (2 - \cos^2 \theta_{\omega})}\right)^2,\quad 2 \boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}} = 2 \sqrt{\dot{\varphi}^2 \dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}} \cos \theta_{\omega}, \end{gathered}\right\} \end{align}
where 
$\varDelta = q_{\rule {0pt}{1.6ex} \mkern -3mu A}{\mkern -6mu}^3 + \frac {27}{4} r_{\rule {0pt}{1.6ex} \mkern -4mu A}{\mkern -6mu}^2$ is proportional to the discriminant of the characteristic equation of 
$\tilde {\!\boldsymbol{\mathsf{A}}}$. In this case, 
$\dot {\epsilon }^2$ and 
$\dot {\gamma }^2$ can be determined directly from 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$ and 
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ and 
$\dot {\varphi }^2$, 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and 
$\zeta$ can be determined if 
$\theta _\omega$ is known. For non-rotational local streamlines, the transformations are given by
\begin{equation} \dot{\epsilon}^2 =-2 q_{\rule{0pt}{1.6ex} \mkern-3mu A},\quad \dot{\gamma}^2 = 1- \dot{\epsilon}^2,\quad \dot{\varphi}^2 = 0, \end{equation}
which can be determined directly from a single parameter, 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$. The rotational and non-rotational transformations are continuous with one another at their boundary (
$\varDelta = 0$) and they are both symmetric about the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$-axis. However, the aspect ratio is only well defined for 
$\varDelta > 0$, consistent with our focus on rotational local streamlines.
In figure 6, we illustrate how 
$\dot {\epsilon }^2$, 
$\dot {\gamma }^2$, 
$\dot {\varphi }^2$, 
$2\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and 
$\zeta$ vary within the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space. In this phase space, the rotational geometries are externally bounded by 
$3\sqrt {3}\lvert r_{\rule {0pt}{1.6ex} \mkern -4mu A} \rvert = ( 1 + q_{\rule {0pt}{1.6ex} \mkern -3mu A} ) ( 1 - 2 q_{\rule {0pt}{1.6ex} \mkern -3mu A} )^{1/2}$ and the non-rotational geometries are externally bounded by 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A} = -\frac {1}{2}$. For the rotational geometries, we display 
$\dot {\varphi }^2$, 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and 
$\zeta$ for 
$\theta _\omega = 43.57^\circ$, which corresponds to the mean value at 
$t = 1.25 t^*$ (see figure 8) and approximates the equilibrium value during turbulent decay. We document how each of these quantities varies with 
$\theta _\omega$ in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space in Appendix D.
Figure 6. The 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space (a) and the variations in 
$\dot {\epsilon }^2$, 
$\dot {\gamma }^2$, 
$\dot {\varphi }^2$, 
$2\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and 
$\zeta$ in this space when 
$\theta _\omega = 43.57^\circ$ is held constant. Panel (a) shows the boundaries corresponding to the symmetry-based (dotted green) and the geometry-based (solid blue) vortex criteria (see § 1.2). It also labels the four classes of non-degenerate local streamline topologies: stable-focus–stretching (SFS), unstable-focus–compression (UFC), stable-node–saddle-saddle (SNSS) and unstable-node–saddle-saddle (UNSS).
Figure 6 shows that, generically, the contribution of pure shearing (
$\dot {\gamma }^2$) tends to dominate the velocity gradients near the origin of the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane. The contribution of normal straining (
$\dot {\epsilon }^2$) grows large when the velocity gradients are dominated by the strain rate tensor. By contrast, for rotational geometries, the contribution of rigid rotation (
$\dot {\varphi }^2$) grows large near the external boundary in regions where the vorticity tensor dominates. The contribution of shear–rotation correlations (
$2\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$) grows largest in the intermediate region of the phase space and it decays near the discriminant line and the external boundary. The aspect ratio 
$\zeta$ is unity along the external boundary and it decays to zero at the discriminant line.
The elliptic instability, which is relevant in strained vortical flows, is expected to be most active for an intermediate range of aspect ratios. Further, vortex stretching and squeezing, which are known to play important roles in turbulent flows, can most readily be associated with the SFS and UFC streamline topologies, respectively (Lopez & Bulbeck Reference Lopez and Bulbeck1993). In the corresponding regions in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane, the transformations in figure 6 suggest that the interplay between shearing and rotation is pertinent to these fundamental turbulent processes. Within this context, the transformations have the potential to characterize the role of the elliptic instability, among other mechanisms, in mediating such processes. In what follows, we investigate this premise by analysing the evolution of the velocity gradient distributions in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$, shear–rotation and 
$\zeta$–
$\theta _\omega$ phase spaces.
5.2. Phase space distributions
The joint probability density functions (p.d.f.s) of the normalized velocity gradients in the phase spaces we investigate encode information about the local streamline geometries. Although these phase spaces are related to one another, the choice of phase space plays an important role in interpreting the statistical distributions of the flow. In the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space, incompressible turbulent flows often follow a near-universal teardrop-like distribution about the origin (Das & Girimaji Reference Das and Girimaji2019, Reference Das and Girimaji2020a). The shear–rotation phase space highlights the distribution of rotational streamline geometries and it characterizes the interplay of rigid rotation and antisymmetric shearing. The 
$\zeta$–
$\theta _\omega$ phase space also considers rotational streamline geometries and it characterizes the flow in terms of geometric features of local streamlines that are associated with the elliptic instability.
Figure 7 shows the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ and the shear–rotation phase space distributions at reference times, selected from table 1, that pertain to the development of the elliptic instability, transition and turbulent decay. As the elliptic instability emerges (
$t \approx 0.75 t^*$), the velocity gradients are concentrated near the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$-axis. Since 
$2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ is relatively large around this time, the rotational regions of the flow are concentrated near the bottom boundary of the shear–rotation phase space, particularly near the location where 
$2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ is maximized.
Figure 7. Joint p.d.f.s of the velocity gradients satisfying 
$A^2/\textrm {max}(A^2) \geq 0.1\,\%$ in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space (a) and 
$W^2/\textrm {max}(W^2) \geq 0.1\,\%$ in the shear–rotation phase space (b) at times 
$t = 0.75 t^*$, 
$0.90 t^*$, 
$1.00 t^*$ and 
$1.25 t^*$ (from top to bottom). The blue triangles represent the centroids of the distributions and the magenta contours represent the p.d.f. levels for which 90 % of the flow (by volume) resides at higher p.d.f. levels. These contours are smoothed by using coarser p.d.f. bins to ensure that they roughly enclose the regions with higher p.d.f. levels.
As the flow transitions (
$t \approx 0.90t^*$), the distribution remains centred about the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$-axis but it expands towards larger 
$\lvert r_{\rule {0pt}{1.6ex} \mkern -4mu A} \rvert$ as it begins to fill the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space. Similarly, the p.d.f. in the shear–rotation phase space begins extending away from the bottom boundary of the phase space, although its bulk remains concentrated at the boundary. The broadening p.d.f.s represent the generation of more diverse local streamline topologies (see figure 6), which is consistent with the formation of more complex vortical structures (see figure 3).
At the time of peak dissipation (
$t \approx t^*$), the distributions populate nearly all of the area in both the phase spaces. During the subsequent turbulent decay (as shown for 
$t \approx 1.25 t^*$), the flow approaches its equilibrium distributions in both phase spaces, which remain similar to the distributions for 
$t \approx t^*$. In the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space, the equilibrium p.d.f. above the 
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$-axis is concentrated slightly left of the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$-axis. This slight preference of the SFS topology is consistent with the typical presence of positive vortex stretching in regions of turbulent flows with rotational geometries. Below the 
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$-axis, the p.d.f. is concentrated along the discriminant line for 
$r_{\rule {0pt}{1.6ex} \mkern -4mu A} > 0$. The equilibrium distribution of our vortex ring collision in this phase space is similar to the near-universal teardrop-like shapes reported previously for forced isotropic turbulence (Das & Girimaji Reference Das and Girimaji2019, Reference Das and Girimaji2020a). This similarity suggests that, in addition to the velocity gradient partitioning (see table 2), the teardrop-like distribution may be more broadly applicable to incompressible flows with appropriate symmetries.
The p.d.f.s in the shear–rotation phase space evolve similarly to those in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space (e.g. by broadening) and specifically highlight rotational geometries. However, the difference between vortex stretching and squeezing, which is encoded in the sign of the real eigenvalue of the VGT, cannot be distinguished in this phase space since its constituents are non-negative. The high concentration of the p.d.f.s near the lower corner of this phase space, representing the origin of the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space, shifts the centroids accordingly and highlights the importance of shearing in the generation and evolution of turbulent flows. The equilibrium distribution remains relatively concentrated about the lower boundary of this phase space, including regions where 2
$\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ is relatively large. This behaviour highlights the potential for the elliptic instability to be active during turbulent decay, but the breadth of the distribution suggests it may not be a completely dominant mechanism driving turbulent flow in rotational regions.
To further investigate the evolution of the flow in the context of the elliptic instability, we again consider the relationship between the 
$\zeta$–
$\theta _\omega$ and shear–rotation phase spaces. Using 
$\theta _\omega$ and 
$\zeta$, the constituents of the shear–rotation phase space and the corresponding shear–rotation correlation term and vortex identification criterion are given by
\begin{equation} \left.\begin{gathered} \frac{\dot{\varphi}^2}{\tilde{W}^2} = \frac{4 \zeta^4 \cos^2 \theta_\omega}{1 + (-2 + 4\cos^2 \theta_\omega)\zeta^2 + \zeta^4},\quad \frac{\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}}{\tilde{W}^2} = \frac{(1 - \zeta^2)^2}{1 + (-2 + 4\cos^2 \theta_\omega)\zeta^2 + \zeta^4},\\ \frac{2 \boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}}}{\tilde{W}^2} = \frac{4 \zeta^2 \cos^2 \theta_\omega (1 - \zeta^2)}{1 + (-2 + 4\cos^2 \theta_\omega) \zeta^2 + \zeta^4}, \quad \varOmega_r = \left(1 + \left( 1 - \sqrt{\frac{\dot{\varphi}^2}{\tilde{W}^2}} \right)^2\right)^{-1}. \end{gathered}\right\} \end{equation}
Figure 8 shows the p.d.f.s of the velocity gradients alongside the distributions of 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}/\tilde {W}^2$ and 
$\varOmega _r$ in the geometry-based phase space at the same times as those in figure 7.
Figure 8. Joint p.d.f.s (a–d) of the velocity gradients in the 
$\zeta$–
$\theta _\omega$ phase space at the same times, at the same points and in the same style as those in figure 7. The bottom plots superimpose the 90 % contours and centroids for 
$t = 0.75 t^*$ (red) and 
$t = 1.25 t^*$ (black) on 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} / \tilde {W}^2$ (e) and on 
$\varOmega _r$ (f), as given by (5.5).
As the elliptic instability emerges (
$t \approx 0.75 t^*$), the p.d.f. is highly concentrated at small values of 
$\theta _\omega$ over a broad range of 
$\zeta$. This distribution is consistent with the enhancement of 
$2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ around this time (see figure 4) and reflects conditions conducive to the elliptic instability. The centroid of the distribution is located around 
$\zeta \approx 0.5$, which suggests that rotational streamlines with this aspect ratio may be particularly susceptible to these conditions.
During transition (
$t \approx 0.90 t^*$), the distribution broadens to a much larger range of 
$\theta _\omega$ and, at the high end of this range, 
$\zeta$ becomes increasingly correlated with 
$\theta _\omega$. This broadening reflects the diversification of rotational local streamline topologies to include those for which shearing and rigid rotation are not well aligned. However, despite this broadening, the p.d.f. is still concentrated at small 
$\theta _\omega$. This behaviour is consistent with the importance of the elliptic instability during transition (see Appendix B).
During turbulent decay (
$t \gtrsim t^*$), the features of the p.d.f. are relatively invariant in time. In the equilibrium distribution (approximated at 
$t \approx 1.25 t^*$), the centroid is located around 
$\theta _\omega \approx 44^\circ$, which is considerably larger than its value (
$\theta _\omega \approx 3^\circ$) at 
$t = 0.75 t^*$. Further, the 90 % contour of the p.d.f. spans nearly the entire range of 
$\theta _\omega$ and highlights a well-defined sharpening in the correlation of 
$\zeta$ with 
$\theta _\omega$ with increasing 
$\theta _\omega$. The upper limit 
$(\theta _\omega,\zeta ) = (90^\circ,1)$ corresponds to circular local streamlines in the plane of rigid rotation subject to out-of-plane shearing.
The shifts in the centroid and the 90 % contour support the hypothesis that, while the elliptic instability still plays a role in the turbulent regime, other mechanisms also contribute significantly to the local flow structure. Specifically, they suggest that, in addition to the breakup of elliptic streamlines via the elliptic instability, the deformation of vortices, e.g. through the action of the Crow instability, may play an important role during turbulent decay. As discussed in Appendix C, the partitioning of shear–rotation correlations between the cores and boundaries of vortices in the flow (see figure 12) provides an interesting opportunity for analysing these mechanisms. Given current challenges in disentangling the elliptic and Crow instabilities in turbulent flows (Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021; Ostilla-Mónico et al. Reference Ostilla-Mónico, McKeown, Brenner, Rubinstein and Pumir2021), the present geometry-based (
$\zeta$–
$\theta _\omega$) phase space has the potential to help distinguish flow features associated with these ubiquitous mechanisms.
Altogether, the results in this section reinforce the notion that the elliptic instability is the dominant mechanism mediating the transition of the present vortex ring collision. They also support the notion that, after transition, the elliptic instability is no longer a strictly dominant mechanism underlying the turbulent decay of the flow. The results point to increased contributions from rotational geometries with out-of-plane shearing, which may reflect interactions associated with mechanisms like the Crow instability. The results also highlight the ability of the new shear–rotation and geometry-based phases spaces to characterize the relative contributions of different modes of deformation and rotational features of local streamlines, respectively.
6. Conclusions
We use a recently developed adaptive, multiresolution numerical scheme based on the LGF to efficiently simulate the head-on collision between two vortex rings at a relatively high Reynolds number (
$Re_{\varGamma _0} = 4000$). The fidelity of this simulation is confirmed using various integral metrics that reflect the symmetries, conservation properties and discretization errors of the flow. We provide a detailed analysis of the initial evolution, transition and turbulent decay of the flow to elucidate flow features that are pertinent to the mechanisms driving its evolution, e.g. the elliptic instability.
Our visualizations of vortex structures enable qualitative characterizations of the various regimes through which the flow evolves. They depict the short-wave elliptic instability as the mechanism driving the initial transition of the rings as they merge at the collision plane. Consistent with previous studies (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021), late transition and (to a lesser extent) turbulent decay are mediated by antiparallel secondary vortex filaments that arise from local interactions associated with the elliptic instability. We confirm that the elliptic instability dominates transition by analysing the scales of dominant wave-like perturbations in that regime. During turbulent decay, we observe local ejections of vortex rings in regions where antiparallel vortex filaments interact. This observation supports the notion that interplay between the elliptic and Crow instabilities can impact vortex interactions, consistent with previous findings (Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021; Ostilla-Mónico et al. Reference Ostilla-Mónico, McKeown, Brenner, Rubinstein and Pumir2021).
Our analysis of the flow centres around using the triple decomposition of the VGT to characterize the contributions of axial straining, shearing, rigid rotation and shear–rotation correlations to the velocity gradients. The mutual interaction of the rings is marked by the development of shearing-dominated vorticity and enhanced shear–rotation correlations, reflecting conditions conducive to the elliptic instability. These conditions are consistent with the initial elliptic instability observed in our visualizations and previously in similar configurations (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021). During turbulent decay, the relative contributions of the different modes of deformation to the velocity gradient strength (which is not stationary) are roughly invariant in time, suggesting an equilibrium partitioning of the VGT. This equilibrium partitioning is remarkably similar to the partitioning observed for forced isotropic turbulence (Das & Girimaji Reference Das and Girimaji2020b), suggesting that it may provide a broadly applicable avenue for modelling incompressible flows with appropriate symmetries.
During the transition and turbulent decay of the flow, we also consider instantaneous distributions of the velocity gradients in various phase spaces. The broadening of the phase space distributions in these regimes reflects the generation of more diverse local streamline topologies. The distributions in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space show that the present vortex ring collision produces velocity gradients that follow the near-universal teardrop-like distribution observed previously for forced isotropic turbulence (Das & Girimaji Reference Das and Girimaji2019, Reference Das and Girimaji2020a). In addition to the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space, we introduce the shear–rotation phase space to characterize the interplay of shearing and rigid rotation in rotational settings and highlight the role of their correlations during transition and turbulent decay.
Finally, we introduce a geometry-based (
$\zeta$–
$\theta _\omega$) phase space to further characterize the action of the elliptic instability (and other mechanisms) during transition and turbulent decay. As the rings interact, the emergence of the elliptic instability spurring transition is associated with the alignment of shearing and rigid rotation (
$\theta _\omega \approx 3^\circ$). In this regime, the elliptic local streamlines in the plane of rigid rotation have aspect ratios centred about 
$\zeta \approx 0.5$. During late transition and turbulent decay, the generation and interaction of secondary vortical structures broadens the distribution to include larger 
$\theta _\omega$, and the equilibrium distribution is ultimately centred near 
$\theta \approx 44^\circ$. In this regime, regions with high 
$\theta _\omega$ and high 
$\zeta$ become increasingly correlated as they approach (
$\theta _\omega,\zeta$) = (
$90^\circ,1$). In conjunction with our visualizations, these results suggest that proximity to vortex cores and boundaries may be a useful tool for modelling the interplay between mechanisms such as the elliptic and Crow instabilities. As a whole, the geometry-based phase space we introduce has the potential to help distinguish effects associated with the elliptic instability (small 
$\theta _\omega$) and other mechanisms, which is an ongoing challenge for turbulent flows driven by interacting vortex filaments (Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021; Ostilla-Mónico et al. Reference Ostilla-Mónico, McKeown, Brenner, Rubinstein and Pumir2021).
Moving forward, the VGT phase spaces we introduce may provide a useful setting for analysing a broad class of turbulent flows. For vortex ring collisions, analysing regimes where the Crow instability dominates the elliptic instability would clarify the extent to which the phase spaces can disambiguate these mechanisms. More generally, it would also be useful to identify the conditions under which (i) the equilibrium partitioning of the VGT (Das & Girimaji Reference Das and Girimaji2020b) and (ii) the teardrop-like distribution in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space (Das & Girimaji Reference Das and Girimaji2019, Reference Das and Girimaji2020a) are applicable. The present VGT analyses are limited by the local, instantaneous nature of the streamline geometries under consideration. It would be interesting to generalize these analyses to capture features that are non-local and that persist in time.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.90.
Acknowledgements
The authors gratefully acknowledge B. Dorschner and K. Yu for their extensive guidance on the computations.
Funding
R.A. was supported by the Department of Defense (DoD) through the National Defense Science and Engineering Graduate (NDSEG) Fellowship Program. This work used the Extreme Science and Engineering Discovery Environment (XSEDE, Towns et al. Reference Towns2014), which is supported by National Science Foundation grant number ACI-1548562. Specifically, it used the Bridges-2 system, which is supported by NSF award number ACI-1928147, at the Pittsburgh Supercomputing Center (PSC).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Computational formulation
Here, we briefly document the adaptive computational framework. We refer to Yu et al. (Reference Yu, Dorschner and Colonius2022) for a detailed description of the framework and a discussion of its novel aspects.
The non-dimensional, incompressible NSE are given by
\begin{equation} \partial_t \boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} =-\boldsymbol{\nabla} p + \frac{1}{Re} \nabla^2 \boldsymbol{u},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}
where 
$\boldsymbol {u} = (u, v, w)$ is the velocity, 
$p$ is the pressure, 
$t$ denotes time and 
$Re$ is the Reynolds number. We focus particularly on the class of unbounded flows obeying the following far-field boundary conditions: 
$\boldsymbol {u}(\boldsymbol {x},t) \rightarrow \boldsymbol {0}$, 
$p(\boldsymbol {x},t) \rightarrow p_\infty$ and 
$\boldsymbol {\omega }(\boldsymbol {x},t) \rightarrow \boldsymbol {0}$ (exponentially) as 
$\lvert \boldsymbol {x} \rvert \rightarrow \infty$. These boundary conditions differ slightly from the more generic (time-varying) free-stream conditions considered by Liska & Colonius (Reference Liska and Colonius2016). For the present simulations, variables are non-dimensionalized using the initial radius and circulation of each vortex ring (
$R_0$ and 
$\varGamma _0$, respectively) and 
$Re$ is given by the initial circulation Reynolds number (
$Re_{\varGamma _0}$).
The NSE are spatially discretized on the composite grid, which contains a series of uniform staggered Cartesian meshes with increasing resolution. Figure 9 depicts the locations of various vector and scalar flow variables on the cells of these meshes. We use 
$\mathcal {Q} \in \{\mathcal {C}, \mathcal {F}, \mathcal {E}, \mathcal {V}\}$ to denote operations that are constrained to the corresponding locations on the cells. The semi-discrete NSE on the composite grid are given by
\begin{equation} \frac{{\rm d} \boldsymbol{\mathsf{u}}}{{\rm d} t} - \boldsymbol{\mathsf{N}}(\boldsymbol{\mathsf{u}}) =-\boldsymbol{\mathsf{G}} \textit{p}_{tot} + \frac{1}{Re} \boldsymbol{\mathsf{L}}_\mathcal{F} \boldsymbol{\mathsf{u}},\quad \boldsymbol{\mathsf{D}} \boldsymbol{\mathsf{u}} = 0, \end{equation}
where we represent discretized variables and operators using sans-serif symbols, which are bold in vector settings and non-bold in scalar settings. Here, 
$\boldsymbol{\mathsf{G}}$, 
$\boldsymbol{\mathsf{D}}$, 
$\boldsymbol{\mathsf{L}}$ and 
$\boldsymbol{\mathsf{N}}$ represent the discrete forms of the gradient, divergence, Laplace and nonlinear operators, respectively. We have used the rotational form of the convective term in (A2a) such that 
$\textit {p}_{tot} = \textit {p} + \frac {1}{2}|\boldsymbol{\mathsf{u}}|^2$ discretely represents the total pressure perturbation and 
$\boldsymbol{\mathsf{N}}(\boldsymbol{\mathsf{u}}) = \boldsymbol{\mathsf{r}}$ discretely represents the Lamb vector, 
$\boldsymbol {r} = \boldsymbol {u} \times \boldsymbol {\omega }$.
Figure 9. Unit cells of the staggered Cartesian grid at the base and refinement levels of the mesh, showing the locations of relevant flow variables.
The semi-discrete momentum equations, subject to the continuity constraint, are integrated in time using the IF-HERK method (Liska & Colonius Reference Liska and Colonius2016; Yu et al. Reference Yu, Dorschner and Colonius2022). This method combines an integrating factor (IF) technique for the viscous term with a half-explicit Runge–Kutta (HERK) technique for the convective term. In the HERK time-stepping scheme (Brasey & Hairer Reference Brasey and Hairer1993; Liska & Colonius Reference Liska and Colonius2016; Yu et al. Reference Yu, Dorschner and Colonius2022), the task of integrating (A2a,b) at each time step is subdivided into 
$N_{stage}$ stages. Using a block lower–upper (LU) decomposition and the mimesis and commutativity properties of the relevant operators, the subproblem associated with stage 
$i$ of the HERK scheme is formulated on the composite grid as
\begin{equation} \textit{L}_\mathcal{C} \textit{p}_{tot}^i = \textit{f}^{\,i} = \boldsymbol{\mathsf{D}} \boldsymbol{\mathsf{r}}^i,\quad \boldsymbol{\mathsf{u}}^i = \boldsymbol{\mathsf{H}}_\mathcal{F}^i (\boldsymbol{\mathsf{r}}^i - \boldsymbol{\mathsf{G}} \textit{p}_{tot}^i). \end{equation}
Here, 
$\boldsymbol{\mathsf{H}}$ represents the IF operator and 
$\textit {f} = \boldsymbol{\mathsf{D}} \boldsymbol{\mathsf{r}}$ approximately represents the divergence of the Lamb vector. For brevity, we omit the exact dependencies of 
$\boldsymbol{\mathsf{r}}^i$ on various flow variables from stages 1 to 
$i$ of the HERK scheme. We refer to the formulation in § 2.4 of Yu et al. (Reference Yu, Dorschner and Colonius2022) for these details and for the corresponding Butcher tableau.
While the discrete operators in (A2a,b) and (A3a,b) are formally defined on the unbounded composite grid, they are practically applied to the finite subset representing the AMR grid. The operator 
$\textit {R}_\mathcal {Q}$ restricts variables from the composite grid 
$( {\cdot } )$ to the AMR grid as 
$\widehat { ( {\cdot } ) } = \textit {R}_\mathcal {Q} ( {\cdot } )$. In the other direction, the operator 
$\textit {P}_\mathcal {Q}$ approximates variables on the composite grid using the values on the AMR grid as 
$( {\cdot } ) \approx \textit {P}_\mathcal {Q} \widehat { ( {\cdot } )}$. Using these operators, solutions to the subproblems associated with each stage of the HERK formulation in (A3a,b) can be approximated on the AMR grid.
The two steps of each subproblem in (A3a,b) involve (i) solving the discrete pressure Poisson equation and (ii) applying the IF to recover the velocity. The solution to the pressure Poisson equation on the AMR grid can be expressed as
\begin{equation}
\widehat{\textit{p}_{tot}^i} =
\textit{R}_\mathcal{C} \mathcal{G}_\mathcal{C} *
\textit{f}^{\,i} = \textit{R}_\mathcal{C}
\textit{L}^{-1}_\mathcal{C}
\textit{f}^{\,i} \approx
\textit{R}_\mathcal{C}
\textit{L}^{-1}_\mathcal{C}
\textit{P}_\mathcal{C}
\widehat{\textit{f}^{\,i}},
\end{equation}
where 
$\mathcal {G}_\mathcal {C}$ is the LGF and 
$*$ represents the discrete convolution. We efficiently evaluate (A4) using a fast multipole method (Liska & Colonius Reference Liska and Colonius2014; Dorschner et al. Reference Dorschner, Yu, Mengaldo and Colonius2020) that accelerates solutions by incorporating summation techniques based on the fast Fourier transform. This method is key to enabling the linear algorithmic complexity and high parallel efficiency of the flow solver. The application of the IF operator can similarly be expressed as
\begin{equation} \widehat{\boldsymbol{\mathsf{u}}^i} = \boldsymbol{\mathsf{R}}_\mathcal{F} \boldsymbol{\mathsf{H}}_\mathcal{F}^i (\boldsymbol{\mathsf{r}}^i - \boldsymbol{\mathsf{G}} \textit{p}_{tot}^i) \approx \boldsymbol{\mathsf{R}}_\mathcal{F} \boldsymbol{\mathsf{H}}_\mathcal{F}^i \boldsymbol{\mathsf{P}}_\mathcal{F} (\widehat{\boldsymbol{\mathsf{r}}^i} - \widehat{\boldsymbol{\mathsf{G}}\textit{p}_{tot}^i}), \end{equation}
where 
$\widehat {\boldsymbol{\mathsf{G}} \textit {p}_{tot}^i } = \boldsymbol{\mathsf{R}}_\mathcal {C} \boldsymbol{\mathsf{G}} \textit {p}_{tot}^i$. The application of the IF operator represents a convolution with an exponentially decaying kernel and it can also be evaluated using fast LGF techniques (Liska & Colonius Reference Liska and Colonius2016; Yu et al. Reference Yu, Dorschner and Colonius2022).
At each time step, the simulation adapts the extent of the AMR grid and adaptively refines regions within the AMR grid according to the spatial adaptivity and mesh refinement criteria, respectively. The spatial adaptivity criterion sets the boundaries of the AMR grid to capture regions where the source of the pressure Poisson equation exceeds a threshold, 
$\varepsilon _{adapt}$, relative to its maximum value in the domain. In other words, the AMR grid is adaptively truncated to capture the subset of the unbounded domain satisfying
\begin{equation} \lvert \textit{f} (\boldsymbol{\mathsf{x}}, \textit{t}) \rvert < \varepsilon_{adapt} \lVert \textit{f} \rVert_\infty (\textit{t}), \end{equation}
where 
$\varepsilon _{adapt} = 10^{-6}$ for the present simulation. One caveat is that the IF convolution involves a velocity source that decays slower than vorticity. Correspondingly, its evaluation requires the velocity field in a slightly extended domain based on a cutoff distance that is selected to capture the IF kernel with high accuracy. For the present simulation, the initial rectangular domain is large enough to contain as a subset the domain satisfying the adaptivity criterion (A6).
At each level 
$k$ of the AMR grid, the mesh refinement criteria are formulated in terms of a combined source, 
$\textit {f}_k(\textit {t})$, which includes the source of the pressure Poisson equation and a correction term. The correction term accounts for the differences between the partial solutions on the coarse and fine grids and it is evaluated using an extended region that can overlap with neighbouring levels. We refer to Yu et al. (Reference Yu, Dorschner and Colonius2022) for the details of its formulation and implementation, which we omit for brevity. Using the combined source, a region is refined or coarsened when
\begin{equation} \textit{f}_k (\boldsymbol{\mathsf{x}}, \textit{t}) > \alpha^{N_{level}-k} \textit{f}_{max}(\textit{t}) \quad {\rm or}\quad \textit{f}_k (\boldsymbol{\mathsf{x}}, \textit{t}) < \beta \alpha^{N_{level}-k} \textit{f}_{max}(\textit{t}), \end{equation}
respectively, where 
$\alpha \in (0,1)$ and 
$\beta \in (0,1)$ and we select 
$\alpha = 0.125$ and 
$\beta = 0.875$ for the present simulation. In these criteria, the combined source is evaluated relative to its maximum blockwise root-mean-square (BRMS) value computed over all blocks and previous times, which is expressed as
\begin{equation} \textit{f}_{max}(\textit{t}) = \max_{\tau < \textit{t}}{\rm BRMS}(\textit{f}_k(\boldsymbol{\mathsf{x}},\tau)). \end{equation}Appendix B. Instability development during transition
Whereas the short-wave elliptic instability has a wavelength of the order of the vortex core radius, 
$a(t)$, the long-wave Crow instability occurs at wavelengths much larger than 
$a(t)$ (Leweke et al. Reference Leweke, Le Dizès and Williamson2016; McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021). Here, we track the development of wave-like instabilities around the azimuth of the ring at the reference times from table 1 associated with the transition to turbulence. Figure 10 shows closeups of the vortical flow structures at these times from the auxiliary viewpoint in figure 3.
Figure 10. Magnified versions of the auxiliary viewpoints in figure 3 at the reference times associated with instability growth and transition.
The dominant scales of the wave-like perturbations are identified by decomposing the flow into azimuthal Fourier modes, which are denoted using 
$\hat {({\cdot })}$ and have corresponding wavenumbers 
$m$. To obtain these Fourier modes, we linearly interpolate the flow at the finest level of the AMR grid to a uniform cylindrical grid. This uniform grid is discretized into 
$N_\theta = 1001$ points in the azimuthal direction and it has spacings, 
$\Delta r_{uni} = \Delta z_{uni} = 0.01$, that are consistent with that of the finest level, 
$\Delta x_{fine} = 0.01$. In the following analysis, we limit our consideration to 
$(r,z)$ pairs for which 
$\langle \omega _z^2 \rangle _\theta$ exceeds 10 % of 
$\max _{r,z} \langle \omega _z^2 \rangle _\theta$, where 
$\langle {\cdot } \rangle _\theta$ denotes azimuthal averaging.
As depicted in figure 11, we characterize instability development using two flow variables. First, following previous studies (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021), we approximate the positions, 
$\boldsymbol {X}^+$ and 
$\boldsymbol {X}^-$, of the vortex cores by identifying the locations of the pressure minima for 
$z > 0$ and 
$z < 0$, respectively. The radial perturbations of these cores, 
$R_p^+(\theta )$ and 
$R_p^-(\theta )$, about their mean radial positions, 
$\bar {R}_p^+$ and 
$\bar {R}_p^-$, characterize the interactions between the rings. We measure the strength of the vortex core perturbations using 
$| \hat {R}_p | = [( | \hat {R}_p^+ |^2 + | \hat {R}_p^- |^2 )/2]^{1/2}$. Second, we consider 
$\omega _z$ since it captures the antiparallel vortex structures that develop around the ring during transition. The Fourier coefficient amplitudes of the corresponding perturbations, 
$| \hat {\omega }_z |$, are coloured according to their displacement from average vortex ring radius, 
$\bar {R}_p = (\bar {R}_p^+ + \bar {R}_p^-)/2$.
Figure 11. Fourier coefficient amplitudes: 
$| \hat {\omega }_z |$ (coloured by radial displacement from 
$\bar {R}_p$) and 
$| \hat {R}_p |$ (black). At each of the reference times shown, 
$\delta$ represents the average slenderness ratio of the vortex rings. The values of 
$| \hat {\omega }_z |$ at wavenumbers pertinent to instability development are marked by symbols and the shaded regions represent wavenumbers for which 
$C > 10$.
To characterize the wavelengths of the dominant perturbations in terms of the elliptic and Crow instabilities, we estimate the core radii, 
$a^+$ and 
$a^-$, of the vortex rings at each reference time. Following McKeown et al. (Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020), we fit a two-dimensional Gaussian function to the vortex core profiles, 
$\langle \omega _\theta \rangle _\theta$, for each ring. However, whereas McKeown et al. (Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020) estimated the core radii by averaging the standard deviations of the Gaussian fits, 
${\sigma ^\pm }_{\mkern -12mu {\rm min}}$ and 
${\sigma ^\pm }_{\mkern -12mu {\rm max}}$, we instead estimate the core radii as 
$a^\pm = [({\sigma ^\pm }_{\mkern -12mu {\rm min}})^2 + ({\sigma ^\pm }_{\mkern -12mu {\rm max}})^2]^{1/2}$. Our estimates produce core radii that are consistent with the definition in (1.1) for Gaussian vortex rings. Hence, we correctly identify 
$a_0 = 0.2$ for the initial condition, which is larger than the alternate core radius definition (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Ostilla-Mónico et al. Reference Ostilla-Mónico, McKeown, Brenner, Rubinstein and Pumir2021) by a factor of 
$\sqrt {2}$. We compute the average slenderness ratio of the vortex rings as 
$\delta = (\delta ^+ + \delta ^-)/2$, where 
$\delta ^\pm = a^\pm /\bar {R}_p^\pm$. This definition is consistent with the Gaussian fits we consider since the radial locations of their centroids coincide with 
$\bar {R}_p^\pm$ to within 0.38 %. The average ratio of the perturbation wavelength to the core thickness can be expressed as 
$C = (C^+ + C^-)/2$, where 
$C^\pm = 2{\rm \pi} /m \delta ^\pm$. Here, we loosely associate 
$C > 10$ and 
$C < 10$ with the Crow and elliptic instabilities, respectively.
At 
$t = 0.75 t^*$, the wavenumber (
$m = 36$) of the largest 
$\hat {\omega }_z$ perturbation corresponds to 
$C = 5.08$. A similar value, 
$C = 5.38$, is obtained at 
$t = 0.90 t^*$ for the 
$m = 52$ perturbation. However, 
$m = 76$ represents the largest perturbation at both 
$t = 0.90 t^*$ and 
$t = 0.95 t^*$, for which 
$C = 3.68$ and 
$3.97$, respectively. At each of these times, the dominant short-wave perturbations are of the order of the core thickness. This result suggests that the corresponding development of secondary antiparallel vortex filaments (see figure 10) can be associated with the development of the elliptic instability.
The vortex core perturbations further support the notion that the elliptic instability is the dominant mechanism in the transitional regime. The perturbations at 
$t = 0.75 t^*$ are too small, relative to 
$\Delta r_{uni} = 0.01$, to resolve. The dominant core perturbations are resolved for 
$t = 0.90 t^*$ and 
$t = 0.95 t^*$ and, consistent with the vorticity perturbations, they are largest at 
$m = 52$ and 
$m = 76$. For both 
$\hat {\omega }_z$ and 
$\hat {R}_p$, we also observe prominent perturbations at 
$m = 80$, but we do not speculate on their source.
For a similar vortex ring collision with 
$Re_{\varGamma _0} = 3500$ and 
$\delta _0 = 0.1$, Mishra et al. (Reference Mishra, Pumir and Ostilla-Mónico2021) attributed the growth of the 
$m = 40$ mode to the elliptic instability in a regime where 
$2 {\rm \pi}R_p(t) / m \approx 0.2\unicode{x2013} 0.4$. Employing a crude volume-conserving approximation for the vortex cores (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021), the core radii can be modelled as 
$a^\pm = a_0 \sqrt {R_0/R_p^\pm }$, which suggests that 
$C \approx 2.26\unicode{x2013} 6.38$ in that regime. Hence, the scales of the dominant perturbations relative to the core thickness in the present case are consistent with those previously attributed to the elliptic instability in a similar collision.
At 
$t = t^*$, we remarkably observe that the dominant 
$\hat {\omega }_z$ perturbation occurs at 
$m = 152$, which is the second harmonic of the 
$m = 76$ perturbation that governed the generation of secondary vorticity. This observation suggests that the elliptic instability retains an important role in mediating the production of subsequent generations of vortical structures at progressively smaller scales. It thus qualitatively supports the initial stages of iterative elliptic instability scenario leading to the generation of turbulence (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020). Identifying the later stages of this pathway would require a more refined analysis of the orientation of each generation of vortices relative to previous generations.
Although perturbations for which 
$C > 10$ are non-negligible, their signatures in 
$\hat {\omega }_z$ are not as prominent as those for which 
$C < 10$. This observation supports the notion that the Crow instability plays a secondary role to the elliptic instability in the transitional regime, consistent with previous studies in similar configurations (McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Mishra et al. Reference Mishra, Pumir and Ostilla-Mónico2021). Nevertheless, especially as the flow becomes turbulent, the broadening range of active scales obscures the interplay between these mechanisms. As this occurs, the vortex core perturbations gain significant energy at lower wavenumbers (
$C > 10$), indicating that long-wave mechanisms like the Crow instability may become important. Altogether, while relatively limited, the present analysis of instability development confirms the pre-eminence of the elliptic instability during transition and supports our interpretation of the corresponding velocity gradients.
Appendix C. Shear–rotation correlations and vortical flow structures
The visualizations in figure 3 help identify antiparallel vortex filaments and interactions between vortices, but the comparisons between the vortex boundary (
$\varOmega _r = 0.52$) and core (
$\varOmega _r = 0.93$) structures provide relatively little information. In figure 12, we visualize the same vortex structures but instead colour them using 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} / \tilde {W}^2$ to probe how conditions conducive to the elliptic instability are structured throughout the vortices in the flow.
Figure 12. Visualizations of the vortex boundaries (
$\varOmega _r = 0.52$, left side) and vortex cores (
$\varOmega _r = 0.93$, right side), coloured by 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} / \tilde {W}^2$, for each reference time from table 1. A movie depicting the evolution of the vortex boundaries from the auxiliary viewpoint (leftmost column) is provided as supplementary material.
As the vortex boundaries merge and expand radially, the shear–rotation correlations are relatively large at the collision plane and the outer boundaries in 
$z$ and they are relatively small at the inner and outer boundaries in the radial direction. This structuring illustrates how shear–rotation correlations are especially enhanced in regions where the vortex boundaries become thinner, corresponding to the shift from a rigid-rotation-dominated regime to a shearing-dominated regime. During transition, the secondary vortex filaments are initially associated with relatively high and low shear–rotation correlations near their boundaries and cores, respectively. As the turbulence develops, this structuring of 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} / \tilde {W}^2$ within the vortices remains similar to that of the secondary vortices mediating transition.
This persistent partitioning opens up an interesting possibility of analysing the action of various mechanisms (e.g. the elliptic and Crow instabilities) in turbulent flows based on their proximity to vortex cores. For example, the phase space transformations in § 5 can be used to characterize local streamline geometries throughout vortices using the structure of the shear–rotation correlations. Consistent with the transformations depicted in figure 8, our results suggest that local streamlines are more elliptic near vortex boundaries and more circular near vortex cores. This conceptual picture is consistent with the notion that the breakup and displacement of vortex core structures can be loosely associated with the elliptic and Crow instabilities, respectively.
Appendix D. Effect of shear–rotation alignment
Here, we characterize the effect of the alignment between shearing and rigid rotation, as measured by 
$\theta _\omega$, on the phase space transformations associated with rotational local streamline geometries. Figure 13 depicts how the corresponding transformations vary with 
$\theta _\omega$ in the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space. When 
$\theta _\omega = 90^\circ$, shearing and rigid rotation occur in orthogonal planes. In this case, 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} = 0$, 
$\zeta = 1$, and the region where 
$\dot {\varphi }^2$ dominates 
$A^2$ extends the furthest from the external boundary of the phase space. When 
$\theta _\omega = 45^\circ$, the regions where 
$\dot {\varphi }^2$ and 
$\zeta$ are large concentrate more sharply near the external boundary and 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ grows in the intermediate region between the boundaries of the rotational geometries. The concentration of 
$\dot {\varphi }^2$ and 
$\zeta$ and the amplification of 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ are most extreme when 
$\theta _\omega = 0^\circ$. In this case, the peak contribution of 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ is 
$( \sqrt {2} + 1 )^{-1}$ (Das & Girimaji Reference Das and Girimaji2020b) and, for all 
$\theta _\omega < 90^\circ$, it occurs when 
$r_{\rule {0pt}{1.6ex} \mkern -4mu A} = 0$. The location of this maximum approaches 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A} \rightarrow \frac {1}{4}$ as 
$\theta _\omega \rightarrow 90^\circ$ and 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A} \rightarrow \frac {1}{2\sqrt {2}}$ as 
$\theta _\omega \rightarrow 0^\circ$. The qualitative features of the distributions vary more significantly from 
$\theta _\omega = 90^\circ$ to 
$\theta _\omega = 45^\circ$ than they do from 
$\theta _\omega = 45^\circ$ to 
$\theta _\omega = 0^\circ$.
Figure 13. Transformations to 
$\dot {\varphi }^2$, 
$2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and 
$\zeta$ from the 
$q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–
$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space for various alignment angles, 
$\theta _\omega$. The plots are in the same style as those in figure 6.
