Breather stripes and radial breathers of the two-dimensional sine-Gordon equation

Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 Contents lists available at ScienceDirect Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns Research paper Breather stripes and radial breathers of the two-dimensional sine-Gordon equation P.G. Kevrekidis a , b , R. Carretero-González c , ∗, J. Cuevas-Maraver d , e , D.J. Frantzeskakis f , J.-G. Caputo g , B.A. Malomed h , i a Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003-4515, USA b Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UK c Nonlinear Dynamical System Group, Computational Science Research Center, and Department of Mathematics and Statistics, San Diego State University, San Diego, CA 92182-7720, USA d Grupo de Física No Lineal, Universidad de Sevilla. Departamento de Física Aplicada I, Escuela Politécnica Superior, C/ Virgen de África, 7, 41011 Sevilla, Spain e Instituto de Matemáticas de la Universidad de Sevilla, Avda. Reina Mercedes, s/n. Edificio Celestino Mutis, 41011 Sevilla, Spain f Department of Physics, National and Kapodistrian University of Athens, Panepistimiopolis, Zografos, Athens 15784, Greece g Laboratoire de Mathématiques, INSA de Rouen Normandie, Saint-Etienne du Rouvray 76801, France h Department of Physical Electronics, School of Electrical Engineering, Faculty of Engineering, Tel Aviv University, Tel Aviv 69978, Israel i Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica, Chile a r t i c l e i n f o Article history: Received 26 July 2020 Revised 19 October 2020 Accepted 30 October 2020 Available online 4 November 2020 Keywords: Nolinear waves sG equation Breathers Modulational instabiity a b s t r a c t We revisit the problem of transverse instability of a 2D breather stripe of the sine-Gordon (sG) equation. A numerically computed Floquet spectrum of the stripe is compared to analytical predictions developed by means of multiple-scale perturbation theory showing good agreement in the long-wavelength limit. By means of direct simulations, it is found that the instability leads to a breakup of the quasi-1D breather in a chain of interacting 2D radial breathers that appear to be fairly robust in the dynamics. The stability and dynamics of radial breathers in a finite domain are studied in detail by means of numerical methods. Different families of such solutions are identified. They develop small-amplitude spatially oscillating tails (“nanoptera”) through a resonance of higher-order breather’s harmonics with linear modes (“phonons”) belonging to the continuous spectrum. These results demonstrate the ability of the 2D sG model within our finite domain computations to localize energy in long-lived, self-trapped breathing excitations. ©2020 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction One of the most influential models in studies of solitary waves is the sine-Gordon (sG) equation, which has been extensively explored in numerous volumes [1–3] and reviews [4] . The one-dimensional (1D) version is the first nonlinear equation whose integrability was found (in the form of the Bäcklund transform, ca. 140 years ago [5] ). Later, its complete integrability ∗Corresponding author. E-mail address: [email protected] (R. Carretero-González). URL: http://nlds.sdsu.edu , http://www.csrc.sdsu.edu (R. Carretero-González) https://doi.org/10.1016/j.cnsns.2020.105596 1007-5704/© 2020 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 was systematically investigated by means of the inverse scattering transform (IST) [6] . Commonly known exact solutions to the sG equation are topological kink solitons and kink-antikink bound states in the form of breathers. Modes of the latter type are uncommon in continuum systems, but find broad realizations in their discrete counterparts [7,8] . While the 1D sG equation has become a textbook etalon of integrable models, far less is known about higher dimensional versions of the same equation, which are not integrable by means of the IST. In particular, some effort has been dedicated to the kinematics and dynamics of kinks in two and three dimensions (2D and 3D) [9–12] , including their ability to produce breathers as a result of collisions with boundaries [13] , and proneness to be stably pinned by local defects [14] . Nevertheless, despite the relevance of physical realizations of the sG model in higher dimensions (see, e.g., Ref. [15] for a recent example), the study of breathers in such settings has been scarce. While effective equations of motion for breathers were derived long ago [16] , these equations do not apply to steady isotropic breather profiles. Hence, the temporal dynamics of the latter states remains somewhat unclear. Moreover, there is a relevant related question concerning quasi-1D breather solutions in the higher-dimensional problem, i.e., breather stripes and planes in the 2D and 3D settings, respectively. The seminal work of Ref. [17] had predicted that such quasi-1D solutions are subject to transverse instabilities. This is a natural property due to the relation of the sG model, in the limit when it produces broad small-amplitude breathers, to the nonlinear Schrödinger (NLS) equation [1,2] , in which both bright and dark quasi-1D solitons are prone to transverse instabilities [18] —see, e.g., theoretical work of Ref. [19] and Refs. [20–22] for experimental realizations in atomic and optical physics. On the other hand, quasi-1D sG kinks are not vulnerable to the transverse instability. In this work, we aim to further explore the dynamics of quasi-1D sG breathers in the 2D geometry, including their transverse instability and its development into radially shaped breathing modes. We start the analysis by revisiting the key result of Ref. [17] for quasi-1D breathers in the 2D sG equation. We test this result for breathers with different frequencies ω b and demonstrate that their instability in the long-wavelength limit is accurately captured by the analytical treatment. Our numerical calculation of the respective Floquet multipliers (FMs) additionally permits the systematic identification of spectra of transverse-instability modes as a function of the corresponding perturbation wavenumber, k y . These results indicate that the quasi-1D sG breathers do not only maintain a main band of unstable wavenumbers, similarly to what is known for the NLS equation, where only a band of low wavenumbers is unstable. Importantly, they also display, for appropriate frequencies ω b , instability bubbles for higher wavenumbers. Near the NLS limit of ω b → 1 , we demonstrate how the NLS reduction, as well as a novel, to the best of our knowledge, variational approximation involving the transverse degree of freedom are able to capture instability and dynamical features, including the basic necking phenomenon occurring in the course of the breather stripe’s evolution. Studying the dynamics ensuing from the generic long-wavelength instabilities of the breather stripes, we identify spontaneous formation of localized breather waveforms. Further, “naively” initializing a quasi-1D radial breather profile, with x replaced by the radial variable, r, we find that the solutions adjust themselves and their frequencies in the course of a transient stage of the evolution, to form potentially very long-lived localized modes trapping the energy in a localized spatial region. This finding motivates us to study the radial time-periodic breathers, which we are able to find, as genuine periodic states, using fixed-point iteration techniques. Stability of these radial breathers is studied numerically in the framework of the Floquet theory. The manuscript is organized as follows. Section 2 reports the results concerning the breather stripes. In particular we expand the instability results for the stripes from Ref. [17] . The results, which are valid from the long-wavelength limit to larger wavenumbers, are obtained by connecting the sG dynamics to the NLS limit, using multiple-scale expansions, and also by dint of a variational approach, which is relevant for large ω b . We present numerical results demonstrating that the development of the breather-stripe instability nucleates long-lived radial breathers. The existence and stability of such radial breathers is the subject of Section 3 . Finally, in Section 4 we present a summary of the results, concluding remarks, and point out avenues for further research. 2. Breather stripes 2.1. Modulational instability of breather stripes The sG equation is a special case of the more general Klein-Gordon class of models which in 2D is written as u tt −∇ 2 u + V  (u ) = 0 , (1) where u (x, y, t) is a real field depending on coordinates (x, y ) and time t, subscripts with respect to the independent variables denote partial derivatives, and V (u ) is the potential that defines the particular model. The most studied among these models pertain to V (u ) = 1 −cos (u ) and V (u ) = (u 2 −1) 2 / 2 , which are referred to as the sG [3] and “φ4 ”[23] equations, respectively. From now on, we focus on the former one: u tt −∇ 2 u+ sin u = 0 . (2) The 1D version of the sG equation admits commonly known kink/antikink and breather solutions. Kinks are (generally speaking, traveling) 1D profiles that asymptote to the uniform background states, u ±∞ = { 0 , 2 π} , of the form u k = 4 arctan  exp s ( x −ct ) /  1 −c 2  , (3) 2 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 where s = ±1 is the topological charge, which distinguishes kinks and antikinks, and cis the velocity, which may take values | c| < 1 . On the other hand, breathers correspond to profiles which are oscillatory in time and localized in space. In particular, the exact 1D sG breather solution is u b (x, t) = 4 tan −1 β ω b sech ( βx ) cos ( ω b t ) , (4) where β≡ 1 −ω 2 b and the band of the breather frequencies is 0 ≤ω b < 1 . In the spirit of the work done in Ref. [14] for a kink stripe, we consider the transverse instability of the breather stripe u 2D (x, y, t) = u b (x, t) . (5) We then perturb the breather stripe according to u (x, y, t) = u b (x, t) + w (x, y, t) where the perturbation w is assumed small. Enforcing this solution to satisfy Eq. (2) yields, to first order, the following evolution equation for the perturbation: w tt −(w xx + w yy ) + cos (u b (x, t)) w = 0 . (6) This equation is a linear, inhomogeneous in x, wave equation so that we can study each transverse wave number k y separately and assume w (x, y, t) = ξ(x, t) exp (ik y y ) , (7) where ξ(x, t) describes the x -dependent shape (eigenmode) of the perturbation. Plugging (7) into (6) yields ξtt −ξxx + k 2 y + cos (u b (x, t) ξ= 0 . (8) This is a periodically forced wave equation and can be solved using Floquet theory, see Appendix for details. The final result is the linear operator M ξ(x, T ) ξt (x, T ) = M ξ(x, 0) ξt (x, 0) , (9) relating the solutions of Eq. (8) at t = 0 and t = T . If the maximal Floquet multiplier (obtained from the eigenvalues of M ) is such that | | > 1 , ξwill grow and the stripe will be unstable. In the case of instability (where at least one eigenvalue has modulus larger than 1) the eigenmode ψcorresponding to the maximum | | gives the solution ξ(x, t) and the modulation of the stripe through its perturbation w (x, y, t) . In the seminal work of Ref. [17] , Ablowitz and Kodama (AK) predicted that, in the long-wavelength limit, breather stripes are subject to modulational instability (MI). Their main result provides the instability growth rate λfor the stripe, as a function of the perturbation wavenumber, k y : λ=σk y , σ2 =  1 −ω 2 b ω b arcsin  1 −ω 2 b . (10) Given that the AK result is valid in the k y → 0 limit, we wish to investigate the instability for larger values of k y . To this end, one needs to solve Eq. (2) numerically, by means of a suitable fixed-point method (see details in the Appendix), taking into regard that a finite size (and a finite discretization parameter) of the solution domain induces gaps in the spectrum of the linear modes and, also that breathers may feature modified tails generated by “hybridization” (mixing) of the third and higher-order odd harmonics with linear modes (“phonons”), reminiscent of the formation of the so-called “nanoptera” in the φ4 equation [3] . The hybridization depends on the location of the linear modes, which, in turn, depends on the specifics of the domain size and boundary conditions. We have used periodic boundary conditions in domain involving x ∈ [ −L x , L x ) with L x = 100 and computed the periodic steady states corresponding to zero initial velocity u t (x, y, t = 0) = 0 . Finally, it is important to note that, addressing the dynamics with the underlying time-periodic solution, one needs to resort to Floquet analysis and the computation of FMs, in order to extract the (in)stability eigenvalues (see details in Appendix). The FMs, , are translated into (in)stability eigenvalues λ(in particular, for the comparison to the AK theory), via = exp (λT ) , where T = 2 π/ω b is the breather period. Fig. 1 depicts the dependence on k y of the MI growth rate for selected values of the breather stripe frequency, ω b . The eigenvalues have been rescaled by σ[see Eq. (10) ] to facilitate the comparison to the AK analytical prediction, given by Eq. (10) and shown by the dashed line in the figure. The MI spectrum indeed follows the AK prediction as k y → 0 , for all values of ω b , i.e., for long perturbation wavelengths, the MI eigenvalue is λ/σ≈k y , independent of ω b . Further, Fig. 1 demonstrates that, beyond the AK limit, λgrows with k y until reaching a maximum, and then decreases, falling to zero at some k y = ˜ k y . We refer to this first instability window ( 0 < k y < ˜ k y ) as the main MI “bubble”. The value ˜ k y , as well as the maximum of λ/σ, increase with ω b for ω b ≤0 . 7 , but this trend changes to a decrease for 0 . 7 < ω b < 1 . This result, which is in stark contrast with properties of MI for the NLS equation [24,25] , has its origin in increasing contributions 3 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 Fig. 1. (Color online) Modulational-instability eigenvalues for the breather stripe vs. wavenumber k y of the modulational perturbations, for the indicated values of the breather frequency. Solid curves depict the eigenvalues numerically determined by unstable Floquet multipliers (see the text). The dashed line represents the analytical Ablowitz-Kodama (AK) result for the long-wavelength limit, k y → 0 [17] , given by Eq. (10) . Fig. 2. (Color online) Snapshots of the breathers u b (x, t = 0) and corresponding eigenfunctions ψfor different MI windows/bubbles. The left and right sets of panels correspond to ω b = 0 . 6 and ω b = 0 . 4 , respectively. Snapshots of the breathers u b (x, t = 0) are displayed in the top panels. The subsequent ones, from top to bottom, display the u (left) and u t (right) components of the eigenfunctions ψcorresponding to the highest-MI wavenumbers in the bubbles of Fig. 1 . For ω b = 0 . 6 , two MI bubbles exist (with maxima located at k y ≈0 . 57 and 1.28), while for ω b = 0 . 4 there are three MI bubbles (with maxima located at k y ≈0 . 43 , 0.94 and 1.30). The computation domain is x ∈ [ −100 , 100) , half of which is shown here. from the third and higher harmonics of the breather’s oscillatory waveform when  1 −ω 2 b decreases and, consequently, truncating at the first harmonic is no longer a good approximation for expression (4) for ω b  0 . 7 . Furthermore, beyond the main MI bubble, Fig. 1 suggests another noteworthy feature: the existence of secondary MI bubbles for k y > ˜ k y . These secondary bubbles get narrower and more spaced with the increase of ω b ; in fact, they are not found for ω b > 0 . 7 . To understand MI in the secondary bubbles, and to compare it to the main MI window, k y < ˜ k y , we explore the shape of the corresponding eigenfunctions of the modulational perturbations. In particular, Fig. 2 depicts the breather profile and eigenfunctions corresponding to unstable eigenvalues in the main MI window and secondary bubbles for ω b = 0 . 6 and ω b = 0 . 4 . In the case of ω b = 0 . 6 (the set of left panels in the figure), MI takes place in the main window and in one secondary bubble, while in the case of ω b = 0 . 4 (the set of right panels) two secondary MI bubbles are present. All perturbation eigenfunctions in these MI bubbles are real as the corresponding eigenvalue responsible for the instability is also real (we depict the respective components u and u t in the left and right plots, respectively). Fig. 2 exhibits eigenfunctions that bear qualitatively similar shapes. Namely, the eigenfunctions in the main MI window are all localized 4 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 without oscillating tails, while in the secondary bubbles, tails are attached to the core of the eigenfunctions, getting stronger in higher-order bubbles. 2.2. Connection to NLS To better understand the MI spectrum for the breather stripes, we resort to the connection between the sG and the NLS equations, in the limit of ω b → 1 . Such a connection can be accurately established via the multiple-scales perturbation method [26] . Below we will briefly present the methodology and provide the main results (see also Refs. [1,2] for a detailed analysis). First we note that, in the limit of ω b → 1 , the sG breather of Eq. (4) has a small amplitude and, hence, one can expand the sG nonlinearity as sin u ≈u −u 3 / 6 . Thus, we seek solutions of the resulting φ4 -equation in the form of the asymptotic expansion: u (x, t) = N  n =1 n u n (x 0 , . . . , x N , t 0 , . . . , t N ) , (11) where u n are functions of the independent variables x j = j x and t j = j t( j = 0 , 1 , 2 , . . . N), while 0 < 1 is a formal small parameter. Substituting ansatz (11) in the equation, and collecting coefficients in front of different powers of yields a set of equations, the first three ones being O (1 ) : L 0 u 1 = 0 , (12) O (2 ) : L 0 u 2 + L 1 u 1 = 0 , (13) O (3 ) : L 0 u 3 + L 1 u 2 + L 2 u 1 −1 6 u 3 1 = 0 . (14) where L 0 ≡∂ 2 t 0 −∂ 2 x 0 + 1 , L 1 ≡2(∂ t 0 ∂ t 1 −∂ x 0 ∂ x 1 ) , L 2 ≡∂ 2 t 1 −∂ 2 x 1 + 2 ( ∂ t 0 ∂ t 2 −∂ x 0 ∂ x 2 ) . To leading order, O () , the solution to Eq. (12) is u 1 = A (x 1 , x 2 , . . . , t 1 , t 2 , . . . ) exp (i ) + c . c , (15) where A is a yet-to-be-determined function of the slow variables, c . c stands for the complex conjugate, and = kx 0 −ωt 0 , with wavenumber k and frequency ωobeying the linear dispersion relation, ω 2 = k 2 + 1 . (16) Substituting solution (15) in Eq. (13) , it is evident that the second term is secular and has to be removed. This yields equation A t 1 + v g A x 1 = 0 , where the v g = k/ωis the sG group velocity, as per Eq. (16) . This result suggests that the amplitude function A depends on variables x 1 and t 1 only through a traveling coordinate, X = x 1 −v g t 1 , (17) i.e., A = A (X, x 2 , x 3 , . . . , t 2 , t 3 , . . . ) . Furthermore, as concerns the solution of Eq. (13) [which is now reduced to L 0 u 2 = 0 ], we take the trivial solution u 2 = 0 , because a nontrivial one would be of the same form as Eq. (15) , hence its amplitude function can be included in field A . Finally, at the next order, O (3 ) , upon substituting u 1 from Eq. (15) and u 2 = 0 , Eq. (14) becomes: L 0 u 3 = 1 −v 2 g A XX + 2 iω ( A t 2 + v g A x 2 ) exp (iθ) + 1 6 A 3 exp (3 iθ) + 1 2 | A | 2 A exp (iθ) + c . c = 0 . (18) It is observed that secular terms ∝ exp (iθ) arise on the right-hand side of Eq. (18) [term ∝ exp (3 iθ) is not secular, because it is out of resonance with the uniform solution]. Removing the secular terms leads to the following NLS equation for complex amplitude A : iA t 2 + 2 ω 3 −1 A XX + ( 4 ω ) −1 | A | 2 A = 0 , (19) where we have removed the group-velocity term by redefining x 2 → x 2 −v g t 2 , and made use of identity 1 −v 2 g = 1 /ω 2 , resulting from Eq. (16) . The fact that the NLS Eq. (19) is valid in the band ω ≥1 , according to Eq. (16) , allows us to consider the limit of ω → 1 , i.e., k → 0 . In this case, with the vanishing group velocity v g → 0 , variable X, defined by Eq. (17) , carries over into x 1 , Eq. (19) reduces to the form of iA t 2 + 1 2 A x 1 x 1 + 1 4 | A | 2 A = 0 . (20) 5 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 Fig. 3. (Color online) Comparison between the (scaled) instability spectrum of the sG breather stripe (shown by colored dots) and the instability spectrum for a stationary NLS bright-soliton stripe (shown by the solid line, as per Refs. [24,25] ). The wavenumber and instability growth rate are rescaled as: ˜ k y ≡k y /  1 −ω 2 b and ˜ λ≡2 λ/ (1 −ω 2 b ) . The stationary soliton solution of the reduced NLS, expressed in terms of original variables x and t, is A = ηsech [(1 / 2)(ηx )] exp [(i/ 8) 2 η2 t] , (21) where ηis an arbitrary O (1) parameter. Thus, an approximate sG breather solution, valid up to order O () near the edge of the phonon band, is produced by Eqs. (15) and (21) : u (x, t) ≈2 η sech 1 2 ηx cos  1 −1 8 2 η2 t  . (22) Note that the frequency of the soliton solution (22) , ω s = 1 −(2 η2 / 8) , is smaller than the cutoff frequency, ω = 1 , which naturally means that this self-trapped excitation belongs to the phonon bandgap , as is known for the sG breather. This fact indicates a strong connection between the stationary NLS soliton and the sG breather in the limit of 0 < 1 −ω b 1 . (23) Indeed, in this limit, observing that the argument of tan −1 in solution (4) is small, it may be approximated by u b ≈ (4 β/ω b ) sech (βx ) cos (ω b t) . Then, letting η ≡2 β, which means that ω b ≈1 −(2 η2 / 8) ≡ω s , we see that the approximate NLS soliton solution, given by Eq. (22) , is identical to the sG breather in the limit case defined by Eq. (23) . Accordingly, it is expected that the breather stripe’s dynamics can be adequately approximated by the NLS Eq. (19) , where ωis identified as the breather stripe frequency, ω b . In particular, the connection between the NLS equation and the sG stripe allows us to approximate the stability spectrum of the sG breather stripe by that for the NLS bright-soliton stripes, which was studied in detail previously [24,25,27] . Specifically, in Fig. 3 we compare the rescaled numerically computed instability spectra of the sG breather stripe (dots) for different breather frequencies from region (23) with those for a stationary NLS bright-soliton stripe, taken from Refs. [24,25] . It is evident that, in the limit of ω b → 1 , the (scaled) instability spectrum of the sG stripe indeed smoothly approaches its NLS counterpart. 2.3. Variational approach for large frequencies In this Section we present a study of the evolution of sG breather stripes, based on the variational approximation, which is valid for relatively large breather frequencies. To this end, we note that the 2D sG Eq. (2) can be derived from the Lagrangian, L =   L dx dy, with density L = 1 2 u 2 t −u 2 x −u 2 y −2 sin 2 u 2 . (24) We rewrite the exact 1D sG breather solution (4) as u b = 4 tan −1 [ ( tan μ) sin ( t cos μ) sech ( x sin μ) ] , (25) where parameter μ≡arcsin  1 −ω 2 b takes values 0 ≤μ≤π/ 2 . We now focus on broad low-amplitude breathers in the limit of μπ/ 2 , approximating the exact 1D breather (25) by the following ansatz : u b = χ(y, t) sech ( μx ) . (26) 6 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 Here, χis a variational parameter accounting for the possibility of variations of the width of the stripe in the y direction, which may represent the so-called necking instability [27,28] . According to the VA, evaluating the Lagrangian with ansatz (26) substituted in Eq. (24) , and deriving the Euler-Lagrange equation from it [29] should yield an effective evolution equation for the width parameter χ( y, t ) . In the present form, the Lagrangian density (24) cannot be integrated analytically for ansatz (26) . Nonetheless, focusing on small-amplitude breather stripes, one may expand the nonlinear term to simplify the Lagrangian density: L = 1 2 u 2 t −u 2 x −u 2 y −1 2 u 2 + 1 24 u 4 . (27) Upon substituting ansatz (26) in Eq. (27) , the integration can be performed in the x direction, yielding the effective Lagrangian: μL eff =  + ∞ −∞ χ2 t −χ2 y + 1 −μ2 3 χ2 + 1 18 χ4 dy. (28) Finally, the variational equation for χ(y, t) , which governs the evolution of necking perturbations of the sG stripe under the condition (23) , is immediately derived from the Lagrangian (28) : χtt −χyy + 1 −μ2 3 χ−1 9 χ3 = 0 . (29) Before we compare the reduced dynamics, presented by Eq. (29) , to the evolution of the full sG breather stripe, we provide a generalization of the VA methodology to incorporate the full nonlinearity of the system. If instead of approximating sin (u ) up to third order as done above, one rewrites the nonlinearity in the effective Lagrangian as sin 2 (u/ 2) = (1 / 2) 1 − ∞ m =0 (−1) m u 2 m / (2 m )! and uses the identity,  + ∞ −∞ sech 2 m (x ) dx = 2 2 m −1 [ ( m −1 ) ! ] 2 ( 2 m −1 ) ! , (30) it is possible to evaluate the effective Lagrangian for the full sinusoidal nonlinearity with ansatz (26) . The respective extended Euler-Lagrange equation is χtt −χyy −μ2 3 χ− ∞  m =1 (−1) m 2 (m −1)! (2 m −1)! 2 (2 χ) 2 m −1 = 0 . (31) Although Eq. (31) generalizes Eq. (29) in that it takes the full nonlinearity into account, it does not provide an exact theory, because the underlying ansatz (26) is only valid for small amplitudes. Thus, the extended equation may provide for a more accurate VA model of the necking dynamics, but this model retains the approximate nature of the above analysis (based on the proposed ansatz). To gain insight into the validity of the VA models (29) and (31) , in Fig. 4 we compare the stability spectra provided by these approximations with the numerically found main instability window of the full sG model. Different panels correspond to replacing the infinite sum in Eq. (31) by the truncated sum  N t m =1 , that keeps the first N t terms in the expansion. In this notation, Eq. (31) reduces to Eq. (29) when N t = 2 . The results presented in Fig. 4 suggest that the MI is indeed adequately captured by the VA in the region (23) . This is true qualitatively in that region and even quantitatively as ω b → 1 . It is seen too that, as 1 −ω b increases, the match between the full sG model and its VA counterparts deteriorates. Nonetheless, the VA reproduces the correct overall trend and the shape of the instability spectra. As regards keeping more terms in the expansion of the nonlinearity in the Lagrangian, Fig. 4 shows that proceeding from N t = 2 to N t = 3 does not produce a significant improvement in the accuracy, and there is no discernible difference between N t = 4 and N t > 4 either. Therefore, from now on, we use the simplest reduced model (29) with N t = 2 for the comparison with the full sG model. 2.4. Dynamical evolution of breather stripes Having tested the reduced models in terms of the stability spectra, we now extend the comparison by numerically following the necking dynamics in the framework of the full sG model and the VA for a number of scenarios. The evolution was initiated by inputs u (x, y, 0) = u b (x, 0) 1 + ε cos n y π L y y , u t (x, y, 0) = 0 , (32) where n y is the perturbation’s integer mode number (associated with the wavenumber n y π/L y ) and εits relative strength, while u b (x, 0) is the exact 1D sG breather profile (4) . Results of the comparison are depicted in Fig. 5 . The dynamics is depicted by means of surface plots of u (x, y, t) at the indicated times. The input was a slightly perturbed exact sG breather stripe, to which one [(a)–(f)] or several [(g) and (h)] perturbation modes with wavenumbers k y were added as follows. 7 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 c Fig. 4. (Color online) Stability spectra of the sG breather (thick lines) and as produced by reduced model (31) (thin lines) for different breather frequencies. The top panel depicts the spectrum produced by Eq. (31) with only N t = 2 [i.e., by Eq. (29) ], instead of the infinite sum in the full sG equation. The second and third panels display results produced by Eq. (31) with N t = 3 and 4, respectively. There is no discernible difference in the spectra for N t ≥4 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) The perturbations were constructed by seeding the necking modulation as per Eq. (32) , with perturbation strength ε = 0 . 01 . Panels (a)–(d) depict the comparison for the perturbations with small wavenumbers n y = 1 , 3, 5, and 7, corresponding to the necking MI. Panels (e) and (f) correspond to large perturbation wavenumbers, with n y = 17 and n y = 19 , which are taken, respectively, just below and above the instability threshold. For the stable case (f) with n y = 19 , a stronger perturbation with ε = 0 . 1 was used to stress the stability of the sG breather against this perturbation wavenumber. Finally, panels (g) and (h) correspond to cases when a mix of perturbation modes with random strengths (chosen in interval −0 . 01 ≤ε ≤ 0 . 01 ) is introduced in the initial conditions. Case (g) contains perturbations with modes 1 ≤n y ≤5 , while in case (h) we used 1 ≤n y ≤22 . In all the cases, to match the underlying oscillation frequencies of the exact breather stripe and its VA counterpart, the time for the VA results is scaled by a fitting factor. This factor, for all the cases shown here, amounts to a reduction of the time between 8 and 10% . As seen in the figure, and as predicted by the stability spectra, perturbations with small wavenumbers always lead to necking MI. The perturbation eventually breaks the stripe into a chain of self-trapped “blobs” (localized modes). These time-oscillating “blobs” are identified as highly perturbed radial breathers, this conclusion being the second key aspect of our study (see the following Section 3 , and also Ref. [13] ). All the modes with perturbation indices 1 ≤n y ≤17 ( k y = n y π/L y ) [see Eq. (32) ] are unstable, and the VA is able to accurately capture the corresponding instability and ensuing dynamics, see panels (a)–(e). Naturally, as the necking perturbation grows, the VA becomes less accurate, as the underlying variational ansatz (26) is no longer an appropriate approximation for the profile of the deformed stripe. Nonetheless, for times up to those at which the stripe breaks into a chain of “blobs”, the VA evolves in tandem with the full sG dynamics. This is also true when, instead of using a single mode as the perturbation, one perturbs the stripe by a collection of modes, as shown in panels (g) and (h) where the perturbations contain, respectively, a combination of the lowest 5 modes and even as many as 22 ones. Finally, in panel (f) of Fig. 5 we present an example of the evolution with the perturbation index ( n y = 19 ) chosen beyond the instability boundary of the main MI window. In this stability case, 8 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 Fig. 5. (Color online) Dynamics ensuing from the destabilization of the unstable breather stripe under different perturbations. Different panels, depicting u (x, y, t) , compare the development of the modulational (necking) instability for the sG breather stripe, as produced by the full sG model (in top subpanels), and by the VA model given in Eq. (29) (in bottom subpanels) for ω b = 0 . 9 . The input was a slightly perturbed exact sG breather stripe, to which one [(a)– (f)] or several [(g) and (h)] perturbation modes with wavenumbers k y are added. The system was integrated in the domain | x | ≤40 , | y | ≤80 (only the segment | x | ≤15 is shown), with periodic boundary conditions along y and homogeneous Dirichlet boundary conditions along x . See the main text for more details. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) the simulations demonstrate that the full sG dynamics is closely followed by the VA for longer times, in comparison to the unstable cases. Fig. 6 expands on the results presented in Fig. 5 by considering random perturbation of the breather stripes and allowing the dynamics to evolve for longer times. We observe that, similarly to Fig. 5 , Fig. 6 shows the formation of long-lived spatially localized “blobs”. The top row of panels depicts snapshots of the field u (x, y, t) at the indicated moments of time. Since the solutions that we are inspecting are oscillatory in time, in the corresponding bottom row we plot snapshots of the corresponding 2D energy density, E(x, y, t) = 1 2 u 2 t + u 2 x + u 2 y + V ( u ) , (33) which allows us to partially eliminate the oscillations, and thus to better capture the emergence and dynamics of the “blobs” (localized breathers). The figure suggests that the breathers generated by the MI from the unstable stripes are generically robust, although some of them disappear, while others spontaneously emerge from condensation of the background random field. Note that the localized breathers travel and, as a result, collide or merge with neighboring ones. In Fig. 7 , we display 3D isocontour plots of the energy density in the (x, y, t) spatiotemporal continuum. This figure makes it possible to better follow the emergence, evolution, and interactions of the “blobs”. 3. Radial Breathers We now turn to the consideration of breathers with axial symmetry, i.e., u (x, y, t) = u (r, t) , with r =  x 2 + y 2 in a circular domain of radius R . Their existence is suggested by the findings reported in the previous Section, where necking MI splits the sG breather stripe into a chain of persistent localized modes (“blobs”). 3.1. Constructing radial breathers With the aim to construct initial conditions that directly lead to radial breathers, we take, as an input, the exact 1D breather (4) , with the x coordinate replaced by the radial one: u radial (r, 0) = u b (r, 0) . The results are summarized in Figs. 8 9 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 nearby ones. The apparent robustness of these modes has prompted us to study their existence and stability in more detail. It is important to note that, contrary to the 1D sG equation, its 2D counterpart is not integrable, hence genuine radial breathers cannot exist in infinite domains for infinitely long times due to resonances of multiple harmonics of the breather with the phonon (continuous) spectrum. Nonetheless, in the finite domain it is possible to find intraband breathers, similarly to how it was done for 1D phantom breathers that were discovered in Ref. [30] . In fact, our results show that such intraband radial breathers, possessing a “nanopteron” spatially oscillatory tail, do exist and may be stable in the finite domain. The full bifurcation structure of such radial breathers with arbitrary values of the oscillation frequency, ω b , is quite complex. Nonetheless, we have performed a detailed numerical study for the existence of broad small-amplitude radial breathers for ω b close to 1, in the vicinity of some of the underlying resonances. We have identified several solution branches and their stability, by means of numerical Floquet analysis, finding that such breathers may indeed be stable in the finite domain. There are numerous avenues deserving further study. For instance, a more complete characterization of the bifurcation scenarios for arbitrary values of ω b may be able to elucidate further nontrivial bifurcations and destabilization scenarios. Also, the dynamics and interactions of radial breathers clearly merit further study: for instance, a challenging objective is to investigate the outcome of interactions between radial breathers (including collisions between moving ones), as a function of their relative phase, velocity, and frequencies. Lastly, it is natural to generalize the consideration for 3D settings, and examine whether spherically symmetric radial breathers may survive in a finite domain, and/or whether they may be produced spontaneously via MI of quasi-2D breather planes embedded in 3D. These topics are presently under consideration, and results will be reported elsewhere. Declaration of Competing Interest No. CRediT authorship contribution statement P.G. Kevrekidis: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Validation, Writing - original draft. R. Carretero-González: Data curation, Funding acquisition, Investigation, Methodology, Project administration, Validation, Writing - original draft. J. Cuevas-Maraver: Investigation, Software, Formal analysis, Visualization. D.J. Frantzeskakis: Formal analysis, Methodology, Writing - review & editing. J.-G. Caputo: Formal analysis, Methodology, Writing - review & editing. B.A. Malomed: Conceptualization, Formal analysis, Methodology, Writing - original draft. Acknowledgments This material is based upon work supported by the US National Science Foundation under Grants PHY-1602994 and DMS-1809074 (P.G.K.). P.G.K. also acknowledges support from the Leverhulme Trust via a Visiting Fellowship and thanks the Mathematical Institute of the University of Oxford for its hospitality during this work. R.C.G. gratefully acknowledges support from the US National Science Foundation under Grant PHY-1603058 . J.C-M. thanks the Regional Government of Andalusia under the project P18-RT-3480 and MICINN, AEI and EU (FEDER program) under the project PID2019-110430GB- C21. B.A.M. appreciates support provided by Israel Science Foundation through grant Np. 1286/17 . Appendix A. Numerical implementation for periodic states and their Floquet stability analysis In this Appendix, we outline the Floquet analysis that was deployed for quasi-1D breathers embedded in rectangular domains and radial breathers in circular domains. These breathers are T -periodic solutions corresponding to steady states of the dynamics after time evolves from t = 0 to t = T . Furthermore, as these solutions only depend (non-trivially) on a single spatial coordinate —the x - or r-direction for stripe and radial breathers respectively—, they can be numerically computed as 1D solutions along this longitudinal direction. Nonetheless, it is crucial to note that the stability for these solutions has to be computed in the full 2D domain where these breather structures are embedded in. Namely, perturbations to these solutions need to be followed, not only along the longitudinal direction, but more importantly, along the transverse direction. To this end, as described below, we follow the growth of plane-wave transverse perturbations with wavenumbers k y and k θ, along the y - and θ-directions for, respectively, the breather stripes in the Cartesian plane and radial breathers in polar coordinates. Thus, each wavenumber, k y or k θ, yields a 1D perturbation equation determining the stability of the corresponding modulations over one period T that is treated, as detailed below, using Floquet analysis to obtain its corresponding eigenvalues (eigenfrequencies) and eigenvectors (eigenmodes). The first step in our numerical computations is to discretize in space the solution and its corresponding derivatives. To this end, we consider a finite-difference scheme with a (spatial) discretization parameter h = 0 . 1 , which transforms the PDE for u (x, t) into a set of Ncoupled ODEs for u n ≡u (x n , t) on the discrete grid { x n } . The discretized version of the quasi-1D sG equation (also known as the Frenkel-Kontorova model) reads: ¨ u n + sin (u n ) + 1 h 2 (u n +1 −2 u n + u n −1 ) = 0 , (37) 16 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 where n = −N/ 2 ... N/ 2 . Similarly, for the radial breather the corresponding discretized radial sG equation may be written as ¨ u n + sin (u n ) + 1 h 2 (u n +1 −2 u n + u n −1 ) + 1 2 nh 2 (u n −1 −u n +1 ) , n = 1 . . . N (38) To produce breathers in the numerical form, we have made use of two different techniques, based on the fact that the solutions are T -periodic, with T = 2 π/ω b : (i) a shooting method, based on the consideration of the map, Y (0) → Y (T ) , Y (t) = { u n (t) } { ˙ u n (t) } , (39) and (ii) a Fourier-transform implementation, based on expressing the solution of the discretized dynamical equations in the form of a truncated Fourier series: u n (t) = z 0 ,n + 2 k m  k =1 z k,n cos (kω b t) , (40) with k m being the maximum of the absolute value of the k in our truncation of the full Fourier series. In the numerics, k m = 11 was chosen. Note that because of the spatial evenness of the sG potential [ V (u ) = 1 −cos (u ) ], the Fourier coefficients with even k are zero and, consequently, only odd harmonics of ω b can resonate with phonons. After the introduction of Eq. (40) in the dynamical equations, one gets a set of N ×(k m + 1) nonlinear, coupled algebraic equations. For a detailed explanation of these methods, the reader is referred to Refs. [38,39] . In both methods, the continuation in the frequency is implemented via the path-following (Newton-Raphson) method. The Fourier-transformed methods have the advantage, among others, of providing an explicit analytical form of the Jacobian, which makes the calculations faster. However, for frequencies smaller than ∼0 . 6 , the Fourier coefficients decay slowly and a large number of coefficients should be kept to produce an accurate solution, which makes the Jacobian numerically expensive. Because of this, we developed the analysis based on method (i) to quasi-1D breathers in the range 0 . 1 ≤ω b < 1 , whereas the Fourier-transform methods (ii) were applied to the radial breathers, when we focused on frequencies around ω b = 0 . 9 . To study the spectral stability of the breathers, we introduce a small perturbation ξn (t) to a given solution u n, 0 (t) of the discretized dynamical equations as u n (t) = u n, 0 (t) + ξn (t) [8] . Then, defining the perturbation as ξn = ξx,n exp (ik y y ) in the quasi-1D breather stripe case yields ¨ ξx,n + cos (u n ) + k 2 y ξx,n + 1 h 2 (ξx,n +1 −2 ξx,n + ξx,n −1 ) = 0 . On the other hand, for the radial breather the corresponding dynamics for the perturbation may be written as ¨ ζr,n + cos (u n ) + k 2 θ n 2 h 2 ζr,n + 1 h 2 (ζr,n +1 −2 ζr,n + ζr,n −1 ) + 1 2 nh 2 (ζr,n −1 −ζr,n +1 ) = 0 , as ζn = ζr,n exp (ik θθ) . Due to the (temporal) periodicity of the solutions, Floquet analysis must be employed. In such a case, the stability properties are determined by the spectrum of the Floquet operator M (whose matrix representation is the monodromy), defined as: { ξn (T ) } { ˙ ξn (T ) } = M { ξn (0) } { ˙ ξn (0) } , (41) for the perturbation ξn in the quasi-1D breather stripe case and an identical equation for ζn for the radial breather case. The 2 N ×2 Nmonodromy eigenvalues ≡exp (i ) are dubbed the FMs (Floquet multipliers), with Floquet exponents (FEs) . A consequence of the fact that the Floquet operator is real is that, if is an FM, ∗is an FM too. Further, because of the symplecticity of the Floquet operator, 1 / is also an FM. In other words, FMs always come in quadruplets (, ∗, 1 / , 1 / ∗) if the monodromy eigenvalues are complex, and in pairs (, 1 / ) if the eigenvalues are real. As a consequence, a necessary and sufficient condition for a breather to be linearly stable is that must be real (i.e., that the corresponding FMs lie on the unit circle in the complex plane). Finally, note that FMs are related to eigenvalues λthrough relation = exp (λT ) . References [1] Dodd RK , Eilbeck JC , Gibbon JD , Morris HC . Solitons and nonlinear wave equations. London: Academic Press; 1982 . [2] Dauxois Th , Peyrard M . Physics of solitons. Cambridge: Cambridge University Press; 2006 . [3] Cuevas J, Kevrekidis PG, Williams FL, editors. The sine-Gordon model and its applications: from Pendula and Josephson junctions to gravity and high energy physics. Heidelberg: Springer-Verlag; 2014 . [4] Kivshar YuS , Malomed BA . Rev Mod Phys 1989;61:763 . [5] Rogers C , Schief WK . Bäcklund and darboux transformations. Cambridge: Cambridge University Press; 2002 . [6] Ablowitz MJ , Kaup DJ , Newell AC , Segur H . Phys Rev Lett 1973;30:1262 . [7] Flach S , Gorbach A . Discrete breathers – advances in theory and applications. Phys Rep 2008;467:1 . [8] Aubry S . Phys D 1997;103:201 . [9] Christiansen PL , Olsen OH . Phys Scr 1979;20:531 . Phys Lett A 1978;68:185; PS Lomdahl, OH Olsen, PL Christiansen, Phys Lett A 1980;78:125; PL Christiansen, PS Lomdahl, Phys D 1981;2:482. [10] Geicke J . Phys D 1982;4:197 . Phys Scr 1984;29:431; Phys Lett A 1983;98:147. 17 P.G. Kevrekidis, R. Carretero-González, J. Cuevas-Maraver et al. Commun Nonlinear Sci Numer Simulat 94 (2021) 105596 [11] Bogolubsky IL , Makhankov VG . JETP Lett 1976;24:12 . [12] Samuelsen MR . Phys Lett A 1979;74:21 . [13] Caputo J-G , Soerensen MP . Phys Rev E 2013;88:022915 . [14] Kevrekidis PG , Danaila I , Caputo J-G , Carretero-González R . Phys Rev E 2018;98:052217 . [15] Fedorov EG , Zhukov AV , Bouffanais R , Malomed BA , Leblond H , Mihalache D , Rosanov NN , Belonenko MB , George TF . Opt Express 2019;27:27592 . [16] Maslov EM . Physica D 1985;15:433 . [17] Kodama Y , Ablowitz MJ . J Math Phys 1980;21:928 . [18] Kivshar YuS , Pelinovsky DE . Phys Rep 20 0 0;331:117 . [19] Kuznetsov EA , Turitsyn SK . Zh Eksp Teor Fiz 1988;94:119 . [Sov Phys JETP 1988;67:1583]. [20] Anderson BP , Haljan PC , Regal CA , Feder DL , Collins LA , Clark CW , Cornell EA . Phys Rev Lett 2001;86:2926 . [21] Tikhonenko V , Christou J , Luther-Davies B , Kivshar YuS . Opt Lett 1996;21:1129 . [22] Gorza S-P , Deconinck B , Emplit Ph , Trogdon T , Haelterman M . Phys Rev Lett 2011;106:094101 . [23] Cuevas-Maraver J , Kevrekidis PG . A dynamical perspective on the φ4 model. Past, present and future. Cham: Springer; 2019 . [24] Cisneros-Ake LA , Carretero-González R , Kevrekidis PG . Phys Rev Res 2019;1:033043 . [25] Yang J . Nonlinear waves in integrable and nonintegrable systems. Society for industrial and applied mathematics; 2010 . [26] Jeffrey A , Kawahara T . Asymptotic methods in nonlinear wave theory. London: Pitman; 1982 . [27] Kuznetsov EA , Rubenchik AM , Zakharov VE . Phys Rep 1986;142:103 . [28] Gorza S-P , Kockaert P , Emplit Ph , Haelterman M . Phys Rev Lett 2009;102:134101 . [29] Malomed BA . Progr Optics 2002;43:71 . [30] Morgante AM , Johansson M , Aubry S , Kopidakis G . J Phys A 2002;35:4999 . [31] Fibich G . The nonlinear Schrödinger equation: singular solutions and optical collapse. Heidelberg: Springer; 2015 . [32] Malomed BA . Eur Phys J Spec Top 2016;s225:2507 . [33] Moll KD , Gaeta AL . Nonlinear optics: materials, fundamentals and applications. Sawchuk A, editor. Optical Society of America; 2002 . [34] Chen C.-A., Hung C.-L.. 2019. arXiv:1907.12550 . [35] Sakaguchi H , Malomed BA . Phys Rev E 2020;101:012211 . [36] Calvo MP , Sanz-Serna JM . SIAM J Sci Comp 1993;14:936 . [37] Kevrekidis PG , Cuevas-Maraver J , Pelinovsky DE . Phys Rev Lett 2016;119:094101 . [38] Archilla JFR , MacKay RS , Marín JL . Phys D 1999;134:406 . [39] Marín JL , Aubry S . Nonlinearity 1996;9:1501 . 18