{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction and Theory"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One of the biggest challenges of dynamical systems theory or nonlinear dynamics is the development of mathematical techniques that provide us with the capability of exploring the geometrical template of structures that governs transport in phase space. Since the early 1900, the idea of pursuing a qualitative description of the solutions of differential equations, which emerged from the pioneering work carried out by Henri Poincaré on the three body problem of celestial mechanics {% cite hp1890 --file LDs %}, has had a profound impact on our understanding of the nonlinear character of natural phenomena. Indeed, this powerful approach has now been widely embraced by the scientific community and its essence was nicely captured by Vladimir Arnold's statement that a complete description of classical mechanics boils down to the geometrical analysis of phase space. \n",
"\n",
"In this section we present the details of a mathematical tool whose potential brings us one step closer to fulfilling the long sought after dream envisioned by Poincaré. The method is known in the literature as Lagrangian descriptors (LDs) and has the capability of revealing the geometrical template of phase space structures that characterizes trajectories with qualitatively distinct dynamical behavior. As we will see, this method provides us with a systematic way of exploring phase space by means of looking at its dynamical skeleton using low-dimensional slices. This procedure allows for a complete reconstruction, that is, a *phase space tomography*, of the intricate geometry of underlying invariant manifolds that characterize transport mechanisms.\n",
"\n",
"Consider a general dynamical system:\n",
"\n",
"\\begin{equation}\n",
"\\dfrac{d\\mathbf{x}}{dt} = \\mathbf{v}(\\mathbf{x},t) \\;,\\quad \\mathbf{x} \\in \\mathbb{R}^{n} \\;,\\; t \\in \\mathbb{R} \\;,\n",
"\\label{eq:gtp_dynSys}\n",
"\\end{equation}\n",
"\n",
"where the vector field satifies $\\mathbf{v}(\\mathbf{x},t) \\in C^{r} (r \\geq 1)$ in $\\mathbf{x}$ and continuous in time. If the vector field does not depend on time, the system is called *autonomous*, and it is *non-autonomous* otherwise. For any initial condition $\\mathbf{x}(t_0) = \\mathbf{x}_0$ this system of differential equations has a unique solution with the same regularity as that of the vector field, which also depends continuously on the initial data {% cite coddington1984 --file LDs %}. The vector field can be specified as an analytical model, or it could also have been retrieved from numerical simulations as a discrete spatio-temporal dataset.\n",
"\n",
"The method of Lagrangian Descriptors was first introduced to analyze Lagrangian transport and mixing in geophysical flows {% cite madrid2009 mendoza2010 --file LDs %}. Initially, LDs relied on the arclength of the trajectories as they evolve from their initial conditions forward and backward in time {% cite mendoza2010 mancho2013lagrangian --file LDs %}. Since its proposal as a nonlinear dynamics tool, it was used to plan a transoceanic autonomous underwater vehicle missions {% cite ramos2018 --file LDs %}, manage marine oil spills {% cite gg2016 --file LDs %}, analyze the Stratospheric Polar Vortex {% cite alvaro1 alvaro2 curbelo2019a curbelo2019b --file LDs %}, and recently in chemical reaction dynamics {% cite craven2015lagrangian craven2016deconstructing craven2017lagrangian revuelta2019unveiling --file LDs %}.\n",
"\n",
"Our goal in this section is to give an introduction to the method of Lagrangian descriptors, which can reveal the geometry of phase space structures that determine transport in dynamical systems, such as in Eq. \\eqref{eq:gtp_dynSys}. The geometry of phase space structures is encoded in the trajectories of the system, which can be extracted using LDs from their initial conditions. The simple idea behind LDs is to seed a given phase space region with initial conditions and integrate a bounded and positive quantity (an intrinsic geometrical and/or physical property of the dynamical system under study) along trajectories for a finite time. This approach reveals the invariant phase space structures that make up the dynamical skeleton governing reaction dynamics. This is similar to the visualization techniques developed in fluid mechanics experiments in the laboratory to uncover patterns of flow structures by studying the evolution of particles in a moving fluid {% cite chien1986 --file LDs %}. A powerful analogy can be found in iron filings that align with the magnetic field lines of a magnet."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One of the biggest challenges that one is faced with when exploring the high-dimensional phase space of a dynamical system, such as those that occur in Hamiltonian systems, is that of the qualitative description of the behavior of ensembles of initial conditions, and recovering from their trajectory evolution the underlying template of geometrical phase space structures that governs the dynamical transport mechanisms of the flow. The problem that naturally arises in this context is that the trajectories of ensembles of initial conditions that start nearby might get ''lost'' with respect to each other very quickly, making the use of classical nonlinear dynamics techniques (Small Alignment Index, Lyapunov Exponents, etc...) that rely on tracking the location of neighboring trajectories computationally expensive and difficult to interpret. Furthermore, when Poincaré maps are applied for the dynamical analysis of high-dimensional systems, one might encounter the issue of trajectories not coming back to the selected surface of section, which would yield no relevant information whatsoever. \n",
"\n",
"The method of Lagrangian descriptors provides tremendous advantages in comparison to other methodologies to overcome these issues. For instance, it is a computationally inexpensive and straightforward to implement tool to explore nonlinear dynamics. But probably, the key and revolutionary idea behind the success of this technique is that it focuses on integrating a positive scalar function along trajectories of initial conditions of the system instead of tracking their phase space location. In this way, by emphasizing initial conditions, it directly targets the building blocks where the dynamical structure of phase space is encoded. The methodology offered by LDs has thus the capability of producing a complete and detailed geometrical *phase space tomography* in high dimensions by means of using low-dimensional phase space probes to extract the intersections of the phase space invariant manifolds with these slices {% cite demian2017 naik2019a naik2019b GG2019 --file LDs %}. Any phase space slice can be selected and sampled with a high-resolution grid of initial conditions, and no information regarding the dynamical skeleton af invariant manifolds at the given slice is lost as the trajectories evolve in time. Moreover, this analysis does not rely on trajectories coming back to the chosen slice, as is required for Ponicaré maps to work. In this respect, there is also another key point that needs to be highlighted which demonstrates the real potential of LDs with respect to other classical nonlinear dynamics techniques. Using LDs we can obtain *all* the invariant manifolds of the dynamical system *simultaneously*, the hyperbolic stable and unstable manifolds coming from *any* NHIM in phase space, and also the KAM tori. This provides an edge over the classical approach of computing stable and unstable manifolds that relies on locating first the NHIMs in phase space individually, and for every NHIM globalize the manifolds separately, for which a knowledge of the eigendirections is crucial. Consequently, the application of LDs offers the capability of recovering *all* the relevant phase space structures in one *shot* without having to study the local dynamics about equilibrium points of the dynamical system."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Formulations for Lagrangian Descriptors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Lagrangian descriptor is a scalar valued non-negative function $M$, that is integrated forward and backward for a fixed intergation time $\\tau$. Originally $M$ was defined using the arclength:\n",
"\n",
"\\begin{equation}\n",
"M(\\mathbf{x_0},t_0,\\tau) = \\int_{t_0-\\tau}^{t_0+\\tau} ||\\mathbf{v}(\\mathbf{x}(t;\\mathbf{x_0},t)|| \\; dt \\;,\n",
"\\label{eq:M_function}\n",
"\\end{equation}\n",
"\n",
"where $||\\cdot||$ is the Euclidean distance. $M$ can be naturally broken down into a forward ($M^f$) and backward ($M^b$) integral:\n",
"\n",
"\\begin{equation}\n",
"M(\\mathbf{x}_{0},t_0,\\tau) = M^b(\\mathbf{x}_{0},t_0,\\tau) + M^f(\\mathbf{x}_{0},t_0,\\tau) \\;,\n",
"\\end{equation}\n",
"\n",
"where we have that:\n",
"\n",
"\\begin{equation}\n",
"M^f(\\mathbf{x}_{0},t_0,\\tau) = \\int^{t_0+\\tau}_{t_0} ||\\mathbf{v}(\\mathbf{x}(t;\\mathbf{x}_0),t)|| \\; dt \\;\\;,\\;\\; \n",
"M^b(\\mathbf{x}_{0},t_0,\\tau) = \\int^{t_0}_{t_0-\\tau} ||\\mathbf{v}(\\mathbf{x}(t;\\mathbf{x}_0),t)|| \\; dt \\;.\n",
"\\end{equation}\n",
"\n",
"The advantage of this splitting is that $M^f$ highlights the stable manifolds of the dynamical system, and $M^b$ recovers the unstable manfifolds, while $M$ shows all the invariant manifolds simultaneously.\n",
"The intuitive idea why this tool works is that the influence of phase space structures on trajectories will result in differences (abrupt change) in arclength of nearby trajectories in the neighborhood of a phase space structure.\n",
"The method captures this distinct dynamical behavior separated by invariant phase space structures that results in abrupt changes of values of $M$. This detection of invariant manifolds has been mathematically quantified in terms of ''singular structures'' {% cite mancho2013lagrangian lopesino2017 --file LDs %}, where $M$ is non-differentiable. Once the manifolds are known, one can compute the NHIM at their intersection by means of a root search algorithm. An alternative method to recover the manifolds and their associated NHIM is by minimizing the function $M$ using a search optimization algorithm. This second procedure and some interesting variations are described in {% cite feldmaier2019 --file LDs %}. \n",
"\n",
"We remark that there is no general ''golden rule'' for selecting the value of $\\tau$ for exploring phase space. The *appropriate* (usually chosen by trial and error) value of $\\tau$ will unveil the relevant geometrical template of phase space structures. A very low value of $\\tau$ will not reveal any structures and very high value may lead to obscurity due to differences in magnitude of LD values. This means that $\\tau$ is intimately related to the time scales of the dynamical phenomena that occur in the system under consideration. One needs to bear in mind the compromise that exists between the complexity of the structures revealed by the method to explain a certain dynamical phenomenon."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An alternative definition of LDs relies on the $p$-norm of the vector field of the dynamical system, where $p \\in (0,1]$ and was first introduced in {% cite lopesino2017 --file LDs %} as:\n",
"\n",
"\\begin{equation}\n",
"M_p(\\mathbf{x}_{0},t_0,\\tau) = \\sum_{k=1}^{n} \\int^{t_0+\\tau}_{t_0-\\tau} |v_{k}(\\mathbf{x}(t;\\mathbf{x}_0),t)|^p \\; dt \\;, \n",
"\\label{eq:Mp_function}\n",
"\\end{equation}\n",
"\n",
"Although this alternative definition of LDs does not have such an intuitive physical interpretation as that of arclength, it has been shown allow for a rigorous analysis of the notion of ''singular structures'' and to establish the mathematical connection of this notion to invariant stable and unstable manifolds in phase space. Furthermore, forward integration will reveal stable manifolds and backward evolution unveils unstable manifolds, as in the arclegth version of LDs. Another important aspect of the $p$-norm of LDs is that, since in the definition all the vector field components contribute separately, one can naturally decompose the LD in a way that allows to isolate distinct dynamical effects such as hyperbolic and elliptic behavior. This was used in \\cite {demian2017,naik2019a} to show that the method can be used to successfully detect NHIMs and their stable and unstable manifolds in Hamiltonian systems. Furthermore, it is important to remark that with this definition of LDs one can mathematically prove that these phase space structures are detected as singularities of the $M_p$ scalar field, that is, at points where the function is non-differentiable and therefore its gradient takes very large values {% cite lopesino2017 demian2017 naik2019a --file LDs %}. Moreover, in this context it has also been shown that:\n",
"\n",
"\\begin{equation}\n",
"\\mathcal{W}^u(\\mathbf{x}_{0},t_0) = \\textrm{argmin } M_p^{(b)}(\\mathbf{x}_{0},t_0,\\tau) \\quad,\\quad \\mathcal{W}^s(\\mathbf{x}_{0},t_0) = \\textrm{argmin } M_p^{(f)}(\\mathbf{x}_{0},t_0,\\tau) \\;,\n",
"\\label{eq:min_LD_manifolds}\n",
"\\end{equation}\n",
"\n",
"where $\\mathcal{W}^u$ and $\\mathcal{W}^s$ are, respectively, the unstable and stable manifolds calculated at time $t_0$ and $\\textrm{argmin}(\\cdot)$ denotes the phase space coordinates $\\mathbf{x}_0$ that minimize the function $M_p$. In addition, NHIMs at time $t_0$ can be calculated as the intersection of the stable and unstable manifolds:\n",
"\n",
"\\begin{equation}\n",
"\\mathcal{N}(\\mathbf{x}_{0},t_0) = \\mathcal{W}^u(\\mathbf{x}_{0},t_0) \\cap \\mathcal{W}^s(\\mathbf{x}_{0},t_0) = \\textrm{argmin} M_p(\\mathbf{x}_{0},t_0,\\tau)\n",
"\\label{eq:min_NHIM_LD}\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this point, we would like to discuss the issues that arise from the definitions of LDs provided in Eqs. \\eqref{eq:M_function} and \\eqref{eq:Mp_function} when they are applied to analyze dynamics in open Hamiltonian systems. Notice that in both definitions all the initial conditions are integrated for the same time $\\tau$. Recent studies have revealed {% cite junginger2017chemical naik2019b GG2019 --file LDs %} that computing fixed-time LDs, that is, integrating all initial conditions chosen on a phase space surface for the same integration time $\\tau$, could give rise to issues related to the fact that some of the trajectories that escape from the potential energy surface can go to infinity in finite time or at an increasing rate. The trajectories that show this behavior will give NaN values in the LD scalar field, hiding some regions of the phase space, and therefore obscuring the detection of invariant manifolds. In order to circumvent this problem we explain here the approach that has been recently adopted in the literature {% cite junginger2017chemical naik2019b GG2019 --file LDs %} known as variable integration time Lagrangian descriptors. In this methodology, LDs are calculated, at any initial condition, for a fixed initial integration time $\\tau_0$ or until the trajectory of that initial condition leaves a certain phase space region $\\mathcal{R}$ that we call the {\\em interaction region}. Therefore, the total integration time in this strategy depends on the initial conditions themselves, that is $\\tau(\\mathbf{x}_0)$. In this variable-time formulation, given a fixed integration time $\\tau_0 > 0$, the $p$-norm definition of LDs with $p \\in (0,1]$ has the form:\n",
"\n",
"\\begin{equation}\n",
"M_p(\\mathbf{x}_{0},t_0,\\tau) = \\sum_{k=1}^{n} \\int^{t_0 + \\tau^{+}_{\\mathbf{x}_0}}_{t_0 - \\tau^{-}_{\\mathbf{x}_0}} |v_{k}(\\mathbf{x}(t;\\mathbf{x}_0),t)|^p \\; dt \\;,\n",
"\\label{eq:Mp_vt}\n",
"\\end{equation}\n",
"\n",
"and the total integration time is defined as:\n",
"\n",
"\\begin{equation}\n",
"\\tau^{\\pm}_{\\mathbf{x}_{0}} = \\min \\left\\lbrace \\tau_0 \\, , \\, |t^{\\pm}|_{\\big| \\mathbf{x}\\left(t^{\\pm}; \\, \\mathbf{x}_{0}\\right) \\notin \\mathcal{R}} \\right\\rbrace \\; ,\n",
"\\end{equation}\n",
"\n",
"where $t^{+}$ and $t^{-}$ are the times for which the trajectory leaves the interaction region $\\mathcal{R}$ in forward and backward time, respectively. It is important to highlighting that we select a large enough interaction region, the variable integration time LD definition given above in Eq. \\eqref{eq:Mp_vt} will approach the fixed-time LD definition in Eq. \\eqref{eq:Mp_function}. Thus, NHIMs and their stable and unstable manifolds will be captured by the phase space points for which the LD is non-differentiable and local minimum behavior given in Eqs. \\eqref{eq:min_LD_manifolds} and \\eqref{eq:min_NHIM_LD} is recovered. Moreover, KAM tori are also detected by contour values of the time-averaged LD. Therefore, the variable integration time LDs provides us with a suitable methodology to study the phase space structures that characterize escaping dynamics in open Hamiltonians, since it avoids the issue of trajectories escaping to infinity very fast."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Revealing the Phase Space Structure"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear System\n",
"\n",
"We describe this result for the two degrees-of-freedom system given by the linear quadratic Hamiltonian \n",
"associated to a index-1 saddle at the origin. This Hamiltonian and the equations of motion are \n",
"given by the expressions:\n",
"\n",
"\\begin{equation}\n",
"H(x,y,p_x,p_y) = \\dfrac{\\lambda}{2}\\left(p_x^2 - x^2\\right) + \\dfrac{\\omega}{2} \\left(p_y^2 + y^2 \\right) \\quad,\\quad \\begin{cases}\n",
"\\dot{x} = \\lambda \\, p_x \\\\\n",
"\\dot{p}_{x} = \\lambda \\, x \\\\\n",
"\\dot{y} = \\omega \\, p_y \\\\\n",
"\\dot{p}_{y} = -\\omega \\, y \n",
"\\end{cases}\n",
"\\label{eq:index1_Ham}\n",
"\\end{equation}\n",
"\n",
"Given the initial condition $\\mathbf{x}_0 = \\mathbf{x}(t_0) = \\left(x_0,y_0,p_x^0,p_y^0\\right) \\in \\mathbb{R}^4$, the general solution is:\n",
"\n",
"\\begin{equation}\n",
"\\begin{split}\n",
"x(t) = \\frac{1}{2} \\left[A e^{\\lambda t} + B e^{-\\lambda t}\\right] \\;\\;&,\\;\\;\n",
"p_x(t) = \\frac{1}{2} \\left[A e^{\\lambda t} - B e^{-\\lambda t}\\right] \\\\[.1cm]\n",
"y(t) = y_0 \\cos(\\omega t) + p^0_y \\sin(\\omega t) \\;\\;&,\\;\\;\n",
"p_y(t) = p^0_y \\cos(\\omega t) - y_0 \\sin(\\omega t) \n",
"\\end{split}\n",
"\\label{eq:gen_sol_linham}\n",
"\\end{equation}\n",
"\n",
"where $A = x_0 + p_x^0$ and $B = x_0 - p_x^0$. Since the system is autonomous, we can choose without loss of generality the initial time as $t_0 = 0$. Decomposing the $p$-norm LDs into the hyperbolic and elliptic components of the system we get:\n",
"\n",
"\\begin{equation}\n",
"M_{p}(\\mathbf{x}_{0},\\tau) = M^{h}_{p}(\\mathbf{x}_{0},\\tau) + M^{e}_{p}(\\mathbf{x}_{0},\\tau) \\;,\n",
"\\label{eq:hyp_comp_M}\n",
"\\end{equation}\n",
"\n",
"where the hyperbolic ($M^{h}$) and elliptic ($M^{e}$) parts are:\n",
"\n",
"\\begin{equation}\n",
"M^{h}_{p}(\\mathbf{x}_{0},\\tau) = \\int^{\\tau}_{-\\tau} |\\dot{x}(t;\\mathbf{x}_0)|^p + |\\dot{p}_{x}(t;\\mathbf{x}_0)|^p \\; dt \\;\\;,\\;\\; M^{e}_{p}(\\mathbf{x}_{0},\\tau) = \\int^{\\tau}_{-\\tau} |\\dot{y}(t;\\mathbf{x}_0)|^p + |\\dot{p}_{y}(t;\\mathbf{x}_0)|^p \\; dt\n",
"\\label{eq:ell_comp_M}\n",
"\\end{equation}\n",
"\n",
"By means of an asympotic analysis when the integration time is sufficiently large, i.e. $\\tau \\gg 1$, it is straightforward to show {% cite lopesino2017 demian2017 naik2019a --file LDs %} that the asymptotic behavior of the hyperbolic component is:\n",
"\n",
"\\begin{equation}\n",
"M_{p}^{h}\\left(\\mathbf{x}_0,\\tau\\right) \\sim \\left(\\dfrac{\\lambda}{2}\\right)^{p-1} \\dfrac{|A|^{p} + |B|^{p}}{p} \\, e^{p \\lambda \\tau}\n",
"\\label{eq:M_hyp_asymp}\n",
"\\end{equation}\n",
"\n",
"Therefore, this result shows that $M_{\\gamma}^{h}$ grows exponentially with $\\tau$ and also that the leading order singularities in $M_{p}^{h}$ occur when $|A| = 0$, that is, when $p^0_x = - x_0$, which corresponds to initial conditions on the stable manifold of the NHIM, or in the case where $|B| = 0$, that is, $p^0_x = x_0$, corresponding to initial conditions on the unstable manifold of the NHIM. Moreover, $M_{p}^{h}$ is non-differentiable at the NHIM, since it is given by the intersection of the stable and unstable manifolds. \n",
"\n",
"\\begin{equation}\n",
"\\lim_{\\tau \\to \\infty} \\dfrac{1}{2\\tau} M^{e}_{p}(\\mathbf{x}_{0},\\tau) = \\dfrac{2}{\\pi} \\left(\\omega R\\right)^{p} B\\left(\\dfrac{p+1}{2},\\dfrac{1}{2}\\right)\n",
"\\label{eq:lim_m_gamma}\n",
"\\end{equation}\n",
"\n",
"where $R$ is the radius of the circular periodic orbit described by the $y$ DoF in the center space $y-p_y$ . So the time average of the elliptic part of LDs converges to a value that depends on the energy of the $y$ DoF, and therefore in the limit its value is constant for all points of the periodic orbit. This result is of interest because it allows us to use LDs as a tool to recover phase space KAM tori, given that the dynamical system under study satisfies the conditions of the Ergodic Partition Theorem {% cite mezic1999 --file LDs %}. Therefore, if one analyzes the time average of LDs on a specific initial condition in phase space $\\mathbf{x}_0$, and its value converges, then this initial copndition would lie in an invariant phase space set cosisting of points that share the same time average value. Therefore, the contours of the time average of LDs, when it converges, identify invariant phase space sets. It is also important to highlight here that the limit value to which $M^{e}$ converges is directly related to the limit value to which the classical arclength definition of LDs tends when just applied to the elliptic component of the system. Indeed:\n",
"\n",
"\\begin{equation}\n",
"\\displaystyle{\\lim_{\\tau \\to \\infty} \\dfrac{1}{2\\tau} M^{e}_{p}} = \\dfrac{2}{\\pi} \\left(\\lim_{\\tau \\to \\infty} \\dfrac{1}{2\\tau} M^{e}\\right)^{p} B\\left(\\dfrac{p+1}{2},\\dfrac{1}{2}\\right)\n",
"\\end{equation}\n",
"\n",
"where $M^e$ represents the arclength LD in Eq. \\eqref{eq:M_function} applied only to the $y$ DoF. \n",
"\n",
"To conclude this theoretical discussion we will show how LDs attains a local minimum at the phase space points corresponding to the stable and unstable manifolds of the NHIM and a global minimum at the NHIM. Given an energy of the system, $H = H_0$, above that of the origin, we know from the Lyapunov Subcenter Manifold Theorem {% cite weinstein1973 moser1976 rabinowitz1982 --file LDs %} that a family of NHIM parametrized by the energy bifurcates from the index-1 saddle. The phase space points that lie on the stable (or unstable) manifold, that is $x_0 = -p_x^0$ ($x_0 = p_x^0$), contribute to the hyperbolic component of LDs described in Eq. \\eqref{eq:M_hyp_asymp} with $|A| = 0$ ($|B| = 0$) so that the Lagrangian descriptor has a local minimum at the manifold. Moreover, if the initial condition is on the NHIM, that is $x_0 = p_x^0 = 0$, then $|A| = |B| = 0$ and consequently $M^{h}_{p}(\\mathbf{x}_{0},\\tau) = 0$. Furthermore, all the energy of the system for these points concentrates on the $y$ DoF which evolves periodically, implying that $M^{e}_{p}(\\mathbf{x}_{0},\\tau) > 0$. As a result, LDs attain a global minimum at those points.\n",
"\n",
"In order to support the theretical argument presented above, we provide also a numerical computation of the $p$-norm LDs in the saddle space $x-p_x$ of the linear Hamiltonian given in Eq. \\eqref{eq:index1_Ham} using $p = 1/2$ and for an integration time $\\tau = 10$. In Fig. [fig:1](#fig:index1_lds) we illustrate how the method detects the stable (red) and unstable (blue) manifolds of the unstable periodic orbit at the origin, and these manifolds can be directly extracted from the ridges of the scalar field $||\\nabla M_p||$. Moreover, we show in Fig. [fig:1](#fig:index1_lds) C) that the LD scalar field is non-differentiable at the manifolds and also attains a local minimum on them, and we do so by looking at the values taken by $M_p$ along the line $p_x = 1/2$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"