{'fraction_non_alphanumeric': 0.027753940653087034, 'fraction_numerical': 0.00993631723951041, 'mean_word_length': 4.4928057553956835, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Painting the Phase Space of Dissipative Systems with Lagrangian Descriptors

In this paper we apply the method of Lagrangian descriptors to explore the geometrical structures in phase space that govern the dynamics of dissipative systems. We demonstrate through many classical examples taken from the nonlinear dynamics literature that this tool can provide valuable information and insights to develop a more general and detailed understanding of the global behavior and underlying geometry of these systems. In order to achieve this goal, we analyze systems that display dynamical features such as hyperbolic points with different expansion and contraction rates, limit cycles, slow manifolds and strange attractors. Furthermore, we study how this technique can be used to detect transition ellipsoids that arise in Hamiltonian systems subject to dissipative forces, and which play a crucial role in characterizing trajectories that evolve across an index-1 saddle point of the underlying potential energy surface.


I. INTRODUCTION
Dissipative dynamical systems are ubiquitous in Nature. For example, in mechanical systems they occur when energy is pumped into the system through an external driving force acting on it, or energy is drained by means of mechanisms such as friction. The process by which dissipation manifests itself in phase space is through the contraction and/or expansion of volumes of ensembles of initial conditions as they evolve in time. These systems display a rich variety of dynamical behavior governed by objects such as limit cycles, attractors, slow manifolds, etc., and their understanding is of paramount importance to unravel the underlying dynamics at play [14,22,29,35,38].
The goal we pursue in this work is to illustrate how the method of Lagrangian descriptors (LDs) can be used to explore the phase space of dynamical systems subject dissipation. We believe that this technique would be of great help to gain insight into the dynamical mechanisms that take place in such systems, and provide some guidance in order to carry out a detailed analysis of any system of interest. The advantage of applying LDs is that it is a scalar-based trajectory diagnostic tool that is very simple to implement computationally. This technique was first developed a decade ago for the study of Lagrangian transport and mixing in geophysical flows [26,30], and in this context it has been applied for instance to the real-time assessment of oil spills [11,13] and also to the path planning of transoceanic missions of autonomous underwater vehicle missions [37]. The identification of hyperbolic trajectories and their stable and unstable manifolds in the time-dependent vector fields defined by ocean currents plays a crucial role in the analysis of these problems, and it was the main focus of these studies. However, other geometrical objects (limit cycles, attractors and slow manifolds) that are relevant for dissipative systems have not been analyzed in the literature with this tool. The contribution of this work is aimed towards filling this gap.
The capability of LDs to reveal phase space structure is addressed here through many classical examples from nonlinear dynamics. The contents of this paper are outlined as follows. In Section II we provide a brief introduction on the method of Lagrangian descriptors, and describe the numerical setup we have used for this tool throughout this work when exploring the phase sapce of dissipative systems. Section III is devoted to presenting the results obtained from our analysis of several well-known model dynamical systems by means of LDs. Firstly, in Sec. III A we look at simple saddle systems that display different expansion and contraction rates, and thus their dynamics are characterized by two distinct time scales. Then, in Sec. III B we move on to illustrate how LDs can be used to detect bifurcations, and in particular, Andronov-Hopf bifurcations that give rise to the birth of limit cycles. Moreover, we demonstrate that LDs are capable of revealing the limit cycle structures in the classical van der Pol oscillator. Section III C is focused on the analysis of systems with slow manifolds, and Sec. III D describes how this tools succeeds in highlighting the intricate and fractal structure of attractors and strange attractors in the Duffing oscillator system. We extend our study of dissipative dynamics further in Sec. III E, where we illustrate the potential that LDs brings for the analysis of transition ellipsoids in Hamiltonian systems subject to dissipative forces. The development of techniques such as LDs for identifying these geometrical objects in phase space is of paramount importance, since they greatly facilitate the dynamical study of trajectories evolving across an index-1 saddle of the underlying potential energy surface that describes the model system. Finally, in Sec. IV we present the conclusions of this work and outlook for future research.

II. LAGRANGIAN DESCRIPTORS
This section briefly introduces the method of Lagrangian descriptors (LDs) that we have used in this work to analyze the phase space structures that govern the dynamics of dissipative systems. This technique is a scalar-based diagnostic tool that was originally introduced in the context of fluid mechanics to study transport and mixing in geophysical flows [26,27,30], but in the past years it has been found to provide useful information for the analysis of the high dimensional phase space of chemical systems, see e.g. [1,2] and references therein. Although the method was initially developed to analyze continuous-time dynamical systems [24,27] with general time dependence, it has also been applied to maps [9,23], stochastic dynamical systems [3] and holomorphic dynamics [10].
The idea behind the method of Lagrangian descriptors is very simple. Take any initial condition and accumulate along its trajectory the values attained by a positive scalar function that depends on the phase space variables. This calculation is carried out both forward and backward in time as the system evolves. Once this computation is done for a grid of initial conditions, the scalar output obtained from the method will highlight the location of the invariant stable and unstable manifolds intersecting this slice. These manifolds will appear as 'singular features' of the scalar field, that is, they are detected at points where the values of LDs display an abrupt change, which indicates distinct dynamical behavior. Forward integration of trajectories detects stable manifolds while backward evolution does the same for unstable manifolds. It is also important to remark here that this technique also reveals the structure of KAM tori, but in this case the tori regions correspond to smooth values of the LD function. In fact, tori can be visualized by computing long-term time averages of LDs as discussed in [9,24].
Consider a general time-dependent dynamical system given by the evolution law: where the vector field defining flow satisfies f (x, t) ∈ C r (R n ) × C(R), with r ≥ 1. The method of LDs was originally defined in [26] as the computation of the arclength of trajectories starting at the initial condition x(t 0 ) when they evolve forward and backward for a fixed integration time τ > 0. In this work we will use the definition of LDs based on the p-norm, which was first presented in [24] to provide a rigorous mathematical foundation for this technique. Given a value p ∈ (0, 1], fix forward and backward integration times, τ f > 0 and τ b > 0 respectively, the diagnostic is given by: where x 0 = x(t 0 ) is any inital condition. Notice that it can be decomposed into its backward and forward components: Although this alternative definition of LDs does not have such an intuitive physical interpretation as that of arclength, it has been shown to provide many advantages. For example, it allows for a rigorous mathematical analysis of the notion of "singular structures" and to establish a mathematical connection of this notion to invariant stable and unstable manifolds in phase space. In [24] it was shown that forward integration reveals stable manifolds, while backward evolution unveils unstable manifolds, at points where the scalar field L p is non-differentiable. Therefore its gradient displays jumps [7,24,33]. This property that the method detects manifolds where the function is nondifferentiable, allows for an easy extraction of these geometrical objects by means of edge detection algorithms similar to those used in image processing, such as Sobel filters and many others. Here, we have applied for this purpose the simple approach of computing the norm of the gradient vector of the scalar field, ||∇L p ||, and also its laplacian ∆L p .
We finish this overview of the method by pointing out that there exists many systems where initial conditions may escape to infinity at a very fast rate, or even in finite time (blow up) when we integrate them forward or backward.
In these situations, one must adapt LDs so that the underlying geometry of the phase space can still be recovered when applying this tool. To do so, it suffices to implement the strategy of accumulating the value of LDs along a trajectory for the whole integration time chosen, or until the trajectory leaves a fixed region of the phase space, whatever happens first. This procedure has been proved to successfully circumvent the issue caused by escaping trajectories on the values attained by the LD scalar field [9,12,17,18].

III. RESULTS
A. Dissipative Saddle Systems

The Linear Saddle
Let us begin our study of the phase space of dissipative systems by analyzing the simplest case of a linear saddle with different expansion and contraction rates. Consider the dynamical system: where λ > 0 is the parameter that measures expansion about the origin, and µ > 0 determines the contraction rate. This system is dissipative, since the divergence of the vector field is λ − µ, which is not zero when λ = µ. The analytical solution of this system has the expression: where (x 0 , y 0 ) is the initial condition of the trajectory at time t 0 = 0. For this system, the origin is a hyperbolic equilibrium point and its stable and unstable manifolds are given by: If we apply the p-norm definition of LDs to Eq. (4) it is straightforward to show that the manifolds of the system are revealed at points where the scalar field L p is non-differentiable. This was shown for the first time in [24]. Consider a grid of initial conditions (x 0 , y 0 ) in the plane at ime t 0 = 0 and integrate them forward and backward for a time interval τ f and τ b respectively. Then, LDs yields: and this expression is non-differentiable at points (0, y 0 ), which correspond to the stable manifold, and also at (x 0 , 0), which is the unstable manifold. Recall that the forward component of LDs highlights the stable manifold whereas the backward contribution reveals the unstable manifold. Notice however that both coordinate axis are singularities, even if we only integrate forward or backward: If the integration times are large enough, i.e. τ f 1 and τ b 1 we get: and this result illustrates that, in forward integration, the singularity that dominates is that corresponding to the stable manifold, while backward integration enhances the unstable manifold. However, an issue arises in this system when it comes to visualizing both manifolds simultaneously with LDs in the same picture, since there are two different timescales at play when λ = µ. There is a competition between the exponential growth, described by e λt , and the decay rate identified with e −µt , and this can obscure phase space structure. We illustrate this behavior in Fig. 1, where we have calculated LDs for p = 1/2 and the system parameters λ = 1 and µ = 2, using an integration time τ f = τ b = 8. The values of LDs are normalized by subtracting the minimum and dividing by the difference between the maximum and minimum. In panels A) and B) we depict L respectively, and we can see how the method nicely highlights the location of the stable and unstable manifolds respectively. However, when we plot L p = L (f ) p in C), only the unstable manifold is visible and we would like to get both manifolds in the same picture. This happens because λ < µ, and therefore, when we compute LD backward in time for an initial condition (0, z 0 ) on the stable manifold using τ b = τ , and integrate a symmetric initial condition (z 0 , 0) on the unstable manifold forward for τ f = τ , this yields: If we divide these two quantities we can compare their magnitude: and, whenever τ 1, this gives: Therefore, if λ < µ, it is clear that L becomes the dominating term. This explains why in Fig. 1 C) the values of LDs on the unstable manifold dominate those attained on the stable manifold, and hence we cannot see the stable manifold when we plot the quantity L p . Nevertheless, we would like to point out that we can easily extract the location of the points corresponding to the stable and unstable manifolds by means of computing the gradient of L ( f ) and L ( b) respectively, and overlay both results in the same plot. This is shown in Fig. 1 D) where we have used the standard coloring scheme widely applied in the literature, that is, red for the unstable manifold and blue for the stable manifold.
This visualization problem in the total LD scalar field was first discussed in detail in [24], where it was shown that one can circumvent this issue by changing the value of p that defines the p-norm used to compute LDs. The argument that was developed in that paper assumed that τ f = τ b , that is, all initial conditions are integrated for the same time forwards and backwards. However, it seems more natural to proceed as follows. Fix a value p ∈ (0, 1] to define the p-norm of LDs, and then adjust the forward and backward integration times τ f and τ b so that the time scales of the problem compensate and the contributions of the forward and backward components of LDs become comparable. We will determine here a condition for τ f and τ b that solves the issue of depicting all the phase space structures simultaneously in the same simulation. Take the initial conditions (z 0 , 0) and (0, z 0 ) on the unstable and stable manifolds respectively, and integrate the one on W u forward in time for τ f , and the one located on W s backward in time for τ b . Impose that the values of LDs obtained are comparable, that is: Taking τ f , τ b 1 and simplifying, this gives: Therefore, we have found a simple formula that determines how to set the forward and backward integration times in terms of the Lyapunov exponents of the saddle system. In particular, for the case λ = µ the formula indicates that forward and backward integration times should be chosen to be equal. In order to validate that this approach helps improving the visualization of the total LD field, we calculate LDs for the model parameters λ = 1, µ = 2 and use a forward integration time of τ f = 8, and a corresponding backward integration time of τ b = 4.346 determined from the formula in Eq. (16). We display the total LD function in Fig. 2 C), and this confirms that we acn adjust the forward and backward integration accordingly to obtain a complete picture of the phase space of the system, despite the different expansion and contraction rates, and hence the timescale difference is captured by the method.

A Nonlinear Saddle Example
Next, we move on to analyze the nonlinear system with a saddle point at the origin: and focus our analysis in the case where λ = −2 and µ = 1. It is a simple exercise to show that this system can be reduced to a non-autonomous (in x) first order linear ODE that can be solved by means of an integrating factor. In fact, the solution is: where (x 0 , y 0 ) is the initial condition of the trajectory at time t 0 = 0. This system has an unstable equilibrium point at the origin and its stable and unstable manifolds are given by: We will show that the method of LDs detects the stable manifold when we integrate forward, whereas the unstable manifold is highlighted by the backward integration. For this computation, note that we can not obtain a closed-form expression for the integrals. Hence, we need to perform an asymptotic approximation similar to the one explained in [24]. Consider first the forward component of LDs: As we can see, this expression is non-differentiable when x 0 = 0, which coincides with the stable manifold of the system. Similarly, the backward component of LDs can be approximated by: Since τ b 1, the leading order singularity in the gradient of the backward component appears in the term that contains |y 0 − x 2 0 /2|. Therefore, it reveals the equation of the unstable manifold of the system. This analytic derivation demonstrates that the method successfully highlights the stable and unstable manifolds of the saddle point at the origin. We illustrate numerically in Fig. 3 that the manifolds are detected. Indeed, panel A) displays the total LD values, and in B) we overlay the LD field with with the stable (blue) and unstable (red) manifolds described in Eq. (19). As a validation, we extract the manifold locations by calculating the gradient of the LD function, ||∇L p ||, which allows us to detect the singular features (image edges) of the scalar output provided by LDs.

Andronov-Hopf Bifurcations
Our next goal is to show how Lagrangian descriptors can be used to detect the presence of limit cycles. We look first at the system described by the 2D normal form for a supercritical Andronov-Hopf bifurcation [28], which is given by the following equations: We will consider here the case where σ > 0. If one writes Eq. (22) in polar coordinates, it is straightforward to show that the system becomes: and therefore, it is clear that when β > 0 the system has a stable limit cycle at the circle of radius β/σ, while for β < 0 the limit cycle disappears, and the origin becomes an asymptotically stable equilibrium point. One needs to be careful about the integration of initial conditions, because in this case, many trajectories will blow up in finite time. This is nicely exemplified if we solve analytically the system in Eq. (23) for the situation where β = 0. For this parameter value, the analytical solution of the system (23) is: where (r 0 , θ 0 ) is the initial condition of the trajectory at time t 0 = 0. We can clearly see that this solution blows up at time t = −1/ 2σr 2 0 < 0, so that one would run into trouble when computing trajectories backward in time. This issue is also present in the solutions to the dynamical system in Eq. (23) for β = 0, and can be easily checked by integrating the ODE using separation of variables.
Our goal now is to apply LDs to the system in Eq. (22) in order to reveal its phase space. Since there is an issue of trajectories blowing up in finite time we will use the strategy of stopping those trajectories that escape a circular region of radius R = 4 centered at the origin. In this way, we avoid problems when computing LDs. In order to reveal the phase space of this system we compute LDs for τ f = τ b = 8 using variable time integration, so that the initial conditions that escape the region of the phase space delimited by a circle of radius R = 4 centered about the origin are stopped. We run experiments for the mode parameter values σ = 1, and β = −0.5, 0, 0.5 and integrate trajectories forwards and backwards for a time τ f = τ b = 8 or until they leave the circle of radius R = 4, what happens first. We present the results of this simulations in Fig. 4. In panel A), which corresponds to β = −0.5, the origin is a stable focus, and the LD values highlight that dynamical behavior. In B), for β = 0, this case is the critical value of the parameter where the system undergoes a supercritical Andronov-Hopf bifurcation. The origin is still a weakly asymptotically stable focus and a limit cycle is about to emerge. However, the LD scalar field seems to be indicating the existence of a limit cycle at points where its values are very large. This is not correct, and we will comment on this issue later. Lastly, in panel C) we show how LDs nicely capture the location of the limit cycle present in the phase space of the system for β = 0.5.
At this point we would like to discuss the false positive given by LDs for the case where β = 0. First, we calculate the forward component of LDs: Thus, it is clear that there are no singularities in the forward integration. However, if we integrate backward, we obtain that: and these expressions are non-differentiable when r 0 = 1/ √ 2στ b . It is important to observe that this is a false positive, since we know that the system in Eq. (23) does not have a limit cycle for β = 0. Notice that this false positive shrinks as the integration time gets large, and it will tend to disappear when τ b goes to infinity. This situation could be interpreted in the sense that the method is identifying the critical value of β for which the system's phase space undergoes a bifurcation, and a limit cycle is about to be born. This potential issue of a false positive shows the importance of carefully checking the output provided by the method when analyzing the dynamics of a system. In summary, if a structure tends to disappear as the integration time gets large, it is reasonable to expect that this structure does not have any dynamical significance, although this should be always verified by other means before drawing any further conclusions.

The van der Pol Oscillator
We focus our attention now on the classical van der Pol oscillator system, which is a paradigmatic example of a system displaying limit cycle behavior [29,38,42,43]. This dynamical system is described by the following second order differential equation:ẍ where µ ≥ 0. As a system of first order ODEs, this can be rewritten as: In order to compute LDs for this system we will use again the strategy of stopping trajectories if they escape a certain phase space region, since backward integration also has blow up issues in finite time. For this purpose we set a circle of radius R = 20 about the origin, and compute LDs forward and backward for τ f = τ b = 50, or until the trajectory leaves this circle, what happens first. To illustrate how the limit cycle of the system gets larger and its shape is distorted as µ is increased, we have simulated the system for µ = 0.1, 0.5, 1.5, 3. Results of the total LD field are shown in Fig. 5, and the panels from left to right correspond to the cases from lower to higher values of µ. It is clear from these plots that the method successfully recovers the geometry and location of the limit cycle for the van der Pol oscillator, and how it deforms as the model parameter value is increased.

C. Analysis of Systems with Slow Manifolds
In this subsection we look at the capability of LDs to reveal slow manifolds in systems that display a separation of time scales [22]. Slow manifolds are relevant in many applications, such as those that involve the dynamical study of gliding animals [34]. We focus our attention first on the dynamical system we introduced in Eq. (17) which has been widely studied in the literature as a basic model for slow-manifold dynamics [4,25]. It is well known that for λ µ < 0 this system has a slow manifold at the curve: We will study here the case where µ = −0.05 and λ = −1. For this situation, the dynamics in the y-coordinate rapidly evolves towards makingẏ very small, and at that point, trajectories follow the asymptotically attracting slow manifold so that dynamics on x become dominant. To reveal the location of the slow manifold we apply LDs and integrate trajectories for τ f = τ b = 5. In Fig. 6 we display from left to right the total, forward and backward LD scalar field respectively, and for comparison we overlay in magenta the location of the slow manifold given by Eq. (29). This test illustrates that the method is correctly highlighting the slow manifold location in the phase space of the system. In order to provide further evidence of the capability of this tool for revealing slow manifolds, we focus next on two classical model dynamical systems that have been widely studied in the literature. First we analyze the classical mechanics problem of a bead in a rotating hoop. A detailed analysis of the stability of this system can be found in [38]. We are interested in solving the dynamical system which is given by: where f (φ) = (µ cos(φ) − 1) sin(φ). We will look at the case where µ > 1 for which the system has two unstable equilibrium points at φ = 0 and φ = π, and two stable equilibria at φ = arccos (1/µ). we will study the slow manifold that exists in the system for ε = 0.02 and µ = 2.3. We compute LDs by integrating trajectories forward and backward for the same time τ f = τ b = 10. The output of these numerical experiments is presented in Fig. 7. As before, we have superimposed in the picture the slow manifold of the system in magenta, which is given by the curve Ω = f (φ), in order to validate the results obtained. We have also depicted with yellow dots the stable equilibria, and magenta dots represent unstable equilibria. Notice that the method also reproduces for this problem the slow manifold, and it is also highlighting the stable and unstable manifolds of the saddle equilibrium points of the system. We check this fact by extracting the geometrical features of the phase space from the LD scalar field by means of applying the laplacian operator to it. Observe that with this approach we can also obtain an accurate approximation of the slow manifold from the ridges of the laplacian field ∆L p . We finish our analysis of the detection of slow manifolds with Lagrangian descriptors by revisiting the van der Pol oscillator problem in Eq. (28) where µ 1. In this case, in order to visualize the phase space of the system we will carry out a Liénard transformation. Consider the following definitions: then we can write Eq. (28) in the new variables x and w as: This system has a slow manifold at the curve w = F (x) = x 3 /3 − x, and in order to test how LDs perform when it comes to detecting this structure, we run a simulation using µ = 10 and we integrate trajectories for τ f = τ b = 50 or until they leave a circle of radius R = 6 centered at the origin. Results are shown in Fig. 8, and we have overlaid in magenta the curve that represents the true slow manifold of the system. This comparison between the output of LDs and the location of the slow manifold underlines the success of the method in highlighting this object. Furthermore, we can extract this information from the gradient of LDs, and obtain an approximation of the slow manifold from the scalar field values of the LD output. Notice that this analysis is also revealing the slow and fast branches of the underlying limit cycle of the system.

D. Identification of Attractors
The next dynamical topic that we would like to address is the capability of the method of Lagrangian descriptors to detect attractors and their basins of attraction. To illustrate this fact, we have chosen the Duffing equation as the model problem to study [8,14,20,21]. The dynamics of the periodically forced and damped Duffing oscillator is described by Newton's second law:ẍ which can be rewritten as a system of first order ODEs in the form: We will consider in our simulations that β = 1. We begin by illustrating how in the conservative unforced case, that is, γ = δ = 0 we get the classical phase space picture of a double well system with a saddle point at the origin and two homoclinic orbits surrounding the center equilibria. To do so, we calculate LDs by using equal integration times forward and backward τ f = τ b = 20 and depict the results in Fig. 9. In panel A) we present the total LD field, and in B) we have extracted the stable and unstable manifolds as ridges of the LD output by means of the gradient. The location of the centers is marked with yellow dots, and the magenta dot at the origin corresponds to the saddle point. Notice how LDs nicely highlight the figure eight shape formed by the homoclinic orbits that enclose the center equilibria of the system. In the next experiment we continue working with the unforced Duffing oscillator (γ = 0) but switch on the damping in the system in order to show that LDs are capable of detecting the boundaries of the basins of attraction of both centers, that now turn into stable foci. We set the damping coefficient to δ = 0.3, and for this parameter value, the eigenvalues at the origin are λ 1 ≈ 0.8612 and λ 2 ≈ −1.1612. Notice that this would imply that we have two different time scales at play in the neighborhood of the origin, and this could difficult the simultaneous visualization of the manifolds when we plot the scalar output of the total LD (forward plus backward), as we discussed for the example of the dissipative saddle. This becomes evident in Fig. 10 A) when we compute the total LDs (forward plus backward contributions) for τ f = τ b = 25. Panels B) and C) show the forward and backward components of LD respectively, and in them we can see that the method reveals the location of the stable and unstable manifolds of the system. These manifolds form the boundaries of the basins of attraction of the stable foci depicted with yellow dots in panel D). Despite the fact that it is difficult to visualize both manifold structures in the total LD field because of the different time scales at play, we can easily extract their locations by analyzing the gradient of the forward and backward components of LD separately and plot these information together as we illustrate in D).
We continue the analysis of attractors with LDs by demonstrating that this tool can be used to reconstruct the intricate geometry of strange attractors. To do so, we explore two different cases for the model parameters and introduce forcing into the system, making it nonautonomous. We look at two different situations, one where α = 1, δ = 0.3, γ = 0.5 and ω = 1.2 (displayed in the left column of Fig. 11), and for the other we select α = 0, δ = 0.05, γ = 7.5 and ω = 1 (right column of Fig. 11). The second case corresponds to the setting for the Ueda attractor [39,40]. The existence and properties of these strange attractors for the Duffing equation have been widely studied in the literature, see e.g. [14,15,21,31] and references therein. In order to reveal the strange attractor and its stable and unstable manifolds at time t = 0 we compute LDs with an integration time τ f = τ b = 20. Since the system is nonautonomous, the attractor is going to evolve and change shape with time, so our analysis is going to provide a snapshot of it at the initial time t = 0. If we want to obtain a picture of the attractor at any other time, say t = t 1 , one just needs to compute LDs using that instant t = t 1 as the initial time. We compare our results with the classical technique used to reconstruct strange attractors that consists in calculating a Poincaré section by strobing the trajectory of an initial condition at integer multiples of the forcing period. For this test we have chosen (1, 0) as the initial condition at time t = 0, and we have recorded its location every period of the forcing term for 15000 periods. In Fig. 11 A) and B) we depict the forward component of LDs, while in the panels below, C) and D), the geometry of the stable (blue) and unstable (red) manifolds has been recovered from the LD values with a gradient-type filter. Finally, and as a validation of the results obtained with LDs, we present in E) and F) the corresponding Poincaré map in purple.
It is important to remark here that the reconstruction of the attractor by means of a Poincaré section approach requires a extremely large integration time, which would be prone to the accumulation of numerical error in the simulation. On the other hand, with the two examples we have discussed above, LDs provide a clear advantage for recovering the structure and location of the attractor, since our experiments were carried out with a small integration time of τ f = τ b = 20, which roughly corresponds to a couple of periods of the forcing, a number that by no means is comparable to the large time interval needed for the Poincaré map, Moreover, Poincaré maps can be only constructed for periodically forced systems, whereas LDs can be applied in a straightforward manner to both aperiodic and periodic forcing scenarios.

E. Transition Ellipsoids in Hamiltonian systems with two degrees of freedom subject to dissipation
We would like to finish this work by exploring the capability that the method of Lagrangian descriptors brings for the dynamical analysis of Hamiltonian systems with two degrees of freedom (DoF) subject to damping effects. In many applications, such as astrodynamics [16,19], structural mechanics [5,49], capsize ship motion [32], chemical reactions [41], etc., it is of crucial importance to address transition events across an index-1 saddle critical point of a potential energy surface (PES). The phase space objects responsible for controlling this transport mechanism are the stable and unstable manifolds of the normally hyperbolic invariant manifold (NHIM) that exists in the vicinity of the index-1 saddle [44]. For 2 DoF Hamiltonian systems, the NHIM is an unstable periodic orbit, that is, it has the geometry of S 1 , and its stable and unstable manifolds have the topology of cylinders (S 1 × R). These 'tube manifolds' are known in the literature as spherical cylinders [6,36] and act as conduits, i.e. a transportation network, that trajectories moving across the index-1 saddle follow along their evolution. Moreover, they are codimension-1 objects in the energy hypersurface of the Hamiltonian system, and thus, they are barriers to transport that trajectories can not cross [45]. These cylinders partition the energy surface into 'reactive' and 'non-reactive' trajectories. All trajectories starting from an initial condition inside the tube manifolds are termed as 'reactive', which means that they will cross the index-1 saddle region of the PES, whereas the initial conditions outside the tubes will not move across the index-1 saddle. However, when damping is added to the Hamiltonian system, it has been recently shown that the topology of the manifolds of the NHIM changes significantly and transforms from a cylindrical geometry into an ellipsoid [46][47][48].
In particular, the stable cylinder of a conservative Hamiltonian system, which characterizes the trajectories crossing from one side to the other of the index-1 saddle in forward time, becomes an ellipsoid, known as a transition ellipsoid. This phase space structure determines all transition events in forward time in the dissipative setup. Furthermore, the unstable periodic orbit associated to the index-1 saddle equilibrium point that existed in the conservative system disappears and converts into a focus-type asymptotic orbit tending to the equilibrium point itself.
In order to demonstrate that LDs can be used as a simple mathematical tool to detect the location and geometry of this ellipsoid, we will analyze an uncoupled Hamiltonian system with 2 DoF defined from a double well potential in the x DoF, and a harmonic oscillator in y. To this model, we will add a linear damping into Hamilton's equations of motion. Consider the Hamiltonian function: where m 1 and m 2 are the masses of the x and y DoF respectively, and ω is the angular frequency of the harmonic oscillator. For our analysis we will set m 1 = m 2 = a = b = ω = 1. We can easily include the effect of dissipation in the system by writing Hamilton's equations in the form: where γ x , γ y > 0 measure the dissipation strength in each DoF. We will look in this work at the case where dissipation has the same influence on both DoF, that is, γ x = γ y = γ, and we fix for all our calculations m 1 = m 2 = a = b = ω = 1.
In order to analyze the dynamics of the system in Eq. (36), we define the following Poincaré surfaces of section (PSOS): The main reasons for choosing these phase space slices are the following. First, Σ 1 serves to obtain a global picture of the basins of attraction of the stable equilibrium points located at the wells of the system. Second, we will use the section Σ 2 to probe the ellipsoidal geometry of the stable manifolds, since it cuts the ellipsoid longitudinally across its equator. Finally, Σ 3 is a transverse plane to the ellipsoid that would yield a circular cross-section, which gives us the classical picture known as a reactive islands [36]. This name comes from the fact that all initial conditions chosen inside the closed curve obtained in this cut are reactive trajectories that cross the index-1 saddle region on forward time, and any other initial condition outside it never crosses the saddle along its evolution.
We start looking at the effect of increasing the strength of the damping in the geometry of the invariant manifolds. It is straightforward to show that the conservative Hamiltonian system (γ = 0) has an equilibrium point of saddle×center stability type at the origin. When dissipation is included, its stability changes, and converts into a saddle×focus. We consider here three different values for the damping coefficient, γ = 0.1, 0.25, 1. If we compute LDs in the PSOS Σ 1 , see Fig. 12, the basins of attraction of the wells in the system are revealed. The spiraling of the unstable manifold that emanates from the equilibrium point at the origin and ends in the equilibrium point in the left and right wells reduces as the damping increases, and this explains why it takes less time for trajectories to approach the stable equilibria. Notice also how the stable manifolds get closer to the origin for large values of the friction, and this represents that the amount of trajectories that can cross from well to well across the saddle point of the PES reduces significantly with large damping values. This implies that the transition ellipsoid becomes smaller in size. Notice that the stable manifold in the total LD scalar field displayed in the left column of Fig. 12 becomes harder to locate as dissipation is increased, but this is expected in the output of LDs, as we have discussed for other examples in this work, since the time scales of the system separate. If we would like to fix this visualization issue, one could integrate initial conditions forward and backward for different times. However, this is not necessary because the manifolds can be directly extracted from the forward and backward components of LDs by means of applying the gradient or laplacian of the LD scalar field, or other edge detection algorithms used in image processing.
In order to clearly visualize the ellipsoidal geometry of the stable manifold we calculate LDs on the configuration space PSOS Σ 2 . We depict in Fig. 13 the LD values (left column) and the stable (blue) and unstable (red) manifolds obtained from the laplacian operator of the LD function (right column). In this phase space slice, the ellipsoidal geometry is nicely captured by LDs, confirming its existence in the system. Also, it shrinks in size as the damping gets large. To finish, we analyze the system further to validate that the structure highlighted by LDs is in fact the stable manifold, and that this transition ellipsoid controls the transport of trajectories across the index-1 saddle. To do so, we calculate LDs on the PSOS Σ 3 using a dissipation of γ = 0.25, and show the results of the forward LD in Fig. 14 A). This slice intersects the ellipsoid transversely into a circle. To test if this closed curve corresponds to the ellipsoid boundary, we select three different initial conditions, the red circle outside the ellipsoid, the blue diamond inside the manifold, and the magenta square on the boundary of the stable manifold: We evolve them forward in time and plot their trajectories in 3D in panel B) and projected onto configuration space in C). This numerical experiment convincingly shows that LDs is correctly detecting the boundary of the transition ellipsoid, since the red The initial energy of the system is H0 = 0.05. Panels in each row correspond, from top to bottom, to dissipation values of γ = 0.1, 0.25, 1, respectively. We have depicted the initial energy boundary of the system as a magenta curve. The first column displays the total LD values, and the column on the right shows the sable (blue) and unstable (red) manifolds extracted from the LD scalar field by means of applying the laplacian operator.
initial condition does not cross the index-1 saddle and is attracted by the stable equilibrium point of the left well. On the other hand, the blue trajectory moves across the saddle point, and then it ends its evolution at the right well equilibrium point. Finally, the magenta trajectory which starts exactly at the boundary of the ellipsoid, behaves as expected, and asymptotically tends towards the saddle point at the origin. It does so following a spiral-like trajectory that evidences the saddle-focus stability nature of the saddle equilibrium.

IV. CONCLUSIONS
In this work we have demonstrated how the method of Lagrangian descriptors can be successfully implemented to reveal the relevant phase space structures of dissipative dynamical systems. This study substantially broadens the applicability spectrum of this tool and exemplifies how a simple scalar diagnostic technique can provide relevant insights in order to understand dynamical phenomena in real problems.
Our analysis of different classical models from nonlinear dynamics has shown that this technique has the capability to detect the presence of limit cycles, basins of attraction, strange attractors and slow manifolds. Moreover, by including the possibility of integrating trajectories forward and backward for different times we have generalized the definition of LDs, and this allows us to account for the separation of time scales displayed in these type of problems. On the other hand, in the context of Hamiltonian systems subject to dissipative forces, where the topology of the stable and unstable manifolds changes due to the presence of friction, we have explored this effect on a double well system with two degrees of freedom. The study of this problem by means of LDs convincingly illustrates the advantage that this tool brings for the location of transition ellipsoids which characterize the phase space trajectories that are allowed to move across the index-1 saddle separating two well regions of the underlying potential energy surface.
In the future, our work will focus on applying Lagrangian descriptors to develop a more detailed understanding on the impact of dissipation on the geometry of invariant manifolds and their interaction, and the role that these phase space structures play on governing the dynamics in many Hamiltonian systems subject to friction that arise in real problems and applications. The initial energy of the system is H0 = 0.05. Panels in each row correspond, from top to bottom, to dissipation values of γ = 0.1, 0.25, 1, respectively. We have depicted the initial energy boundary of the system as a magenta curve. The first column displays the total LD values, and the column on the right shows the sable (blue) and unstable (red) manifolds extracted from the LD scalar field by means of applying the laplacian operator. We have marked three different initial conditions, the red circle is outside the transition ellipsoid, the blue diamond is inside, and the magenta square is on its boundary. The magenta curve depicts the initial energy boundary of the system. B) Three dimensional dynamics in forward time of the initial conditions selected in panel A) as they evolve from t0 = 0 to t = 20. We have superposed LDs calculated on the phase space slices Σ2 and Σ3, and the green surface represents the initial energy shell. C) Configuration space projection of the trajectories in panel B), superposed with LDs calculated on the section Σ2. Potential energy contours are depicted as black cuves and the outer grey area corresponds the forbidden region for the initial energy of the system.