Propagation of solitons in a two-dimensional nonlinear square lattice

We investigate the existence of solitary waves in a nonlinear square spring–mass lattice. In the lattice, the masses interact with their neighbors through linear springs, and are connected to the ground by a nonlinear spring whose force is expressed as a polynomial function of the masses out-of-plane displacement. The low-order Taylor series expansions of the discrete equations lead to a continuum representation that holds in the long wavelength limit. Under this assumption, solitary wave solutions are sought within the long wavelength approximation, and the subsequent application of multiple scales to the resulting nonlinear continuum equations. The study focuses on weak nonlinearities of the ground stiffness and reveals the existence of 3 types of solitons, namely a ‘bright’, a ‘dark’, and a ‘vortex’ soliton. These solitons result from the balance of dispersive and nonlinear effects in the lattice, setting aside other relevant phenomena in 2D waves such as diffraction that may lead to a field that does not change during propagation in nonlinear media. For equal constants of the in-plane springs, the governing equation reduces to the Klein–Gordon type, for which bright and dark solitons replicate solutions for one-dimensional lattices. However, unequal constants of the in-plane springs aligned with the two principal lattice directions lead to conditions in which the soliton propagation direction, defined by the group velocity, differs from the wave vector direction, which is unique to two-dimensional assemblies. Furthermore, vortex solitons are obtained for isotropic lattices, which shows similarities with results previously found in optics, thermal media and quantum plasmas. The paper describes the main parameters defining the existence of these solitary waves, and verifies the analytical predictions through numerical simulations. Results show the validity of obtained solutions and illustrate the main characteristics of the solitary waves found in the considered nonlinear mechanical lattice. The study provides an analysis of the physics of waves in nonlinear systems, and may lead to novel designs of devices that can be used for high-performance waveguides.


Introduction
The interest in waves propagating in nonlinear and dispersive media was prompted by Scott Russell [1] who observed in a water channel what he called wave of translation, a single wave in the surface of water that can travel undisturbed over very large distances.In 1895, Korteweg and de Vries developed a model (KdV) for waves on shallow water surfaces that could explain this observed phenomenon, since then called soliton or solitary waves [2].In a general sense, the term soliton is used to refer to any localized field that does not change during propagation due to a balance between dispersive and nonlinear effects.Soliton solutions have been also found in magnetic fields [3], electrical networks [4], or in biological systems [5,6], among others.In the field of nonlinear optics, the study of solitons has been particularly active [7].The balance between nonlinear effects and diffraction (spatial soliton), or between in periodic media is a very active research field.These systems can be widely used in mechanical, acoustic, electromagnetic, electronic and opto-mechanical applications.For one dimensional (1D) mechanical lattices, Zabusky and Kruskal [10] first showed that the continuum limit of Fermi-Pasta-Ulam lattice [11] was the KdV equation [2], and was therefore susceptible to admit soliton-like solutions.The Toda chain, with neighbor interaction [12] also constitutes a paradigmatic example of 1D lattice able to develop solitons; the Toda chain has motivated numerous works in this field [13][14][15].Also, Iizuka et al. [16] studied the wave propagation in random lattices with different nonlinearities, leading to soliton solutions.Recently, Nadkarni et al. [17] studied a 1D periodic chain of bistable elements; the continuum limit approach of this system in the long wave-length regime leads to a 1D nonlinear Klein-Gordon equation which, under certain conditions, exhibits the existence of bright and dark solitons.
In two dimensional (2D) lattice systems, it is necessary to cite the work by Toda [18], who studied, for the first time, the presence of solitons in (2D) mechanical systems.Later, Sreelatha and Joseph [19] presented a continuum treatment of soliton behavior in these kind of lattices.Interestingly, intense focussing effect and the generation of compact acoustic pulses (sound bullets) have been created by [20] by using a tunable, nonlinear acoustic lens, which consists of ordered arrays of granular chains.In contrast to optical bullets that are theorized to reside in optically dispersive/diffractive media, sound bullets arise when solitary waves created within the nonlinear granular medium coalesce in a linear nondispersive medium.Wang and Liu [21] investigated the dynamics of a two-dimensional system with nonlinear interactions, finding solitary wave solutions in certain directions for an hexagonal configuration.Deng et al. [22] found vector solitons in a structure comprising a network of squares connected by highly deformable ligaments.In any case, the number of studies devoted to solitons in 2D mechanical systems is comparatively smaller.
The paper presents a model for a 2D anisotropic square lattice resting on a weakly nonlinear substrate, with masses being submitted to out-of-plane displacement.The low-order Taylor series expansions of the discrete equations lead to a continuum representation that holds in the long wavelength limit.The continuum nonlinear equations are solved using a multiple-scale expansion [23].The leading order solution consists in a carrier wave modulated in amplitude by different solitons, which are derived from 1D or 2D nonlinear Schrödinger equations of the third order.These solitons result from the balance of dispersive and nonlinear effects in the lattice, setting aside other relevant phenomena in 2D waves such as diffraction that, as previously stated, may lead to a field that does not change during propagation in nonlinear media.The study shows that the propagation of bright and dark solitons is essentially one dimensional, although the lattice anisotropy affects their orientation with respect to the carrier wave.The existence of a third soliton type, namely the vortex soliton, previously identified in optics [24][25][26], thermal medium [27] or quantum plasmas [28], has specific bi-dimensional features.The propagation of these solitons in the lattice is verified by numerical solution of the discrete equations.The simulations provide the opportunity for the analysis of the physics of the considered nonlinear lattice, which may suggest exploiting of solitary waves for novel devices and waveguides.
The paper is organized as follows.After this introduction, Section 2 describes the lattice configuration, its dispersion properties and the continuum approximation adopted for the subsequent analytical solutions, which are presented in Section 3. Next, Section 4 presents the numerical simulation of the nonlinear lattice response, while Section 5 summarizes the main findings of the study and provides recommendations for followon investigations.

Lattice configuration and governing equations
We consider the propagation of waves in a spring mass lattice of square topology assembled in the ,  plane (Fig. 1(a)).The lattice consists of masses that undergo motion only in the direction perpendicular to the lattice plane, i.e. the  direction in Fig. 1.Linear springs with constants   and   connect neighboring masses in the  and  directions respectively, and provide a restoring force in direction , proportional to the relative vertical displacement between adjacent particles.In addition, a mass in the generic position ,  within the lattice is connected to the ground by a nonlinear spring.The behavior of this nonlinear spring is defined by a potential  (   ) , with   denoting the  direction displacement of the th mass.Thus, the force exerted by the nonlinear spring is generally expressed as: and the Hamiltonian of the lattice is given by: ] .
Accordingly, the momentum balance equation of the mass at location ,  is: We expand the force of the ground spring up to the third order: where  1 is the linear force coefficient while the coefficients of the nonlinear terms are obtained as: ,  3 = − 1 6 Based on this definition, the momentum balance equation can be rewritten as follows: where the following parameters and variables are introduced: with

Linear dispersion relation for the lattice
We first describe the linear dispersion relations of the lattice, which serve as reference and provide basic information on the propagation of plane waves.To this end, a linear approximation for the force in the ground spring is considered, i.e.: where, based on the definition of the force potential in Eq. ( 4), the relations f (0) = 0 and A plane wave solution is imposed in Eq. ( 7) of the form: where  0 is the magnitude, and  is the frequency associated with the wave vector  =    +   .Also,   ,   denote the position of the  As customary in the analysis of periodic systems, the assumed plane wave solution is rewritten in terms of a non-dimensional wave vector  =    +    =    +   .Substituting the resulting plane wave expression into the governing equation (Eq.( 7)) and further introducing a non-dimensional frequency: leads to the following dispersion relation for the linear lattice: A final parameter is introduced to account for the difference in the spring constants along the  and  directions.Namely, letting: gives: Examples of dispersion relations are presented in Fig. 2, where the surface defined in Eq. ( 15) is represented in the form of iso-frequency contours.The case of  = 1 (Fig. 2) is characterized by iso-frequency contours which in the long wavelength regimes appear as circular as a result of the isotropic nature of the lattice.In contrast, the case of  = 10 (Fig. 2) illustrates the difference in propagation which manifests itself even at low frequencies corresponding to the long wavelength regime.This frequency/wavenumber range, and the considerations herein are of importance in the context of the continuum approximation for the behavior of the lattice, and the analysis of the propagation of solitary waves based on such approximation, which is described in the next section.

Continuum limit and long wavelength approximation
The limiting case of long wavelengths is particularly interesting since the displacement amplitudes change only slowly from mass to mass.The lattice size has then no part to play and transition to a continuum is possible.Then, continualization may serve to provide a connection between parameters of discrete and equivalent continuum media, and to identify general properties of the discrete mechanical system.In the continuum limit, the displacement of first neighbors can be written using a Taylor expansion centered at the particle , , thus 4  12 4  12 where u = u(, ) denotes the continuous variable approximation of the out-of-plane motion of the lattice.The continuum counterpart of the governing Eq. ( 7) becomes 4  12 4  12 Assuming a plane wave solution in the linear regime, f (u) = −u, leads to the dispersion relation for the continuum which is given by: Upon truncation to the first Taylor expansion term in Eqs. ( 16), (17), which restricts the study to the long-wavelength, i.e.  = √  2  +  2  ≪ 1, Eq. ( 18) reduces to: The resulting equation also governs the dynamics of a membrane lying on a nonlinear foundation, which is submitted to different prestress in the directions  and .Thus, the results presented next are readily transferable to membrane-type systems consisting in ultra-thin sheets bonded or in contact to elastomeric or fluid substrates as found for example in MEMs applications [29,30].Else, it may guide the definition of an experimental set-up for the verification of the solutions found herein, which may be the object of future investigation.Within the longwavelength approximation, the coefficient  defined in Eq. ( 14) can be considered as an anisotropy index.Thus, in the isotropic case ( = 1), Eq. ( 20) is recognized as a 2D nonlinear Klein-Gordon type, which arises in various problems in science and engineering [31].
The linear approximation of the ground interaction leads to the corresponding dispersion relation which has the well-known form: Fig. 4 compares the approximate continuum dispersion relation in the long-wavelength regime, with the exact discrete one.The difference between both surfaces are small for |  | ≪ 1, |  | ≪ 1, which specifically illustrates the range of validity of the first order approximation, expressed in Eq. ( 20) and employed in the remainder of this paper.As expected, the higher order continuum approximation provides a better estimate of the dispersion relations, and therefore of the dynamic behavior of the lattice, over a broader range of wavenumbers and frequencies.

Multiple scale solution of the continuous equations
The continuous approximation truncated to the first order is employed as the basis for the analysis of the considered lattice in the presence of nonlinearities.To this end, we restrict the analysis to the case of weak nonlinearities so that a perturbation approach can be employed.
The considered continuum equation (Eq.( 20)) is further simplified by introducing the following set of non-dimensional variables: The proposed change of coordinates {, } constitutes a uniform scaling which preserves angles between lines.
For moderate displacement amplitudes the nonlinear force is approximated by a third-order Taylor expansion and the governing equation becomes with

Multiple-scale expansion
Eq. ( 23) is solved by considering a multiple-scale expansion [23], which leads to a solution process similar to the one previously employed for the 1D version of this problem in Nadkarni et al. [17].Consider the third-order ansatz: where  ≪ 1 is a book keeping parameter that enforces that wave motion occurs for small amplitude and therefore that the results obtained below hold in the weakly nonlinear regime.Also, the following scaled spatial and temporal scales are introduced, Accordingly, the unknown functions   in the solution expansion depend on the scaled variables, i.e.

Solution of the ordered equations
Substituting the ansatz (25) in Eq. ( 23) and equating terms of like order (up to  3 ), leads to the following set of ordered equations: • Order  2 : • Order  3 : where ℒ[⋅] is the following linear differential operator, corresponding to the linear governing equation expressed in Eq. ( 20): (31)

Solution for order 𝜀 1
The solution for the first order equation is in the form of a plane wave with unknown amplitude modulation, and is expressed in the following form: In Eq. ( 32) above and in the remainder of the paper, c.c. denotes the complex conjugate, while the parameter , which corresponds to the imposed phase modulation, is introduced for convenience: Here, ,   and   are related by the linear dispersion relation corresponding to the considered continuum description, which is expressed in Eq. ( 21).

Solution for order 𝜀 2
Following the solution for order  1 (Eq.( 32)), Eq. ( 29) can be rewritten as where the coefficients of  i in Eq. ( 28) are resonant and therefore lead to a secular term which is set equal to zero: Considering the ansatz for the non-homogeneous solution of Eq. ( 34) and equating oscillating and non-oscillating terms gives: and finally The non-oscillating term in  1 accounts for the potentially asymmetric response of the system due to an eventual non-nil value of the coefficient  2 .According to Eqs. ( 6) and ( 24),  2 represents the curvature of the force F (ū) at the origin, whose sign defines whether softening appears in compression or in tension.

Solution for order 𝜀 3
Substituting the solutions for  0 and  1 in Eq. ( 30) yields where the secular term in Eq. ( 35) is enforced, and where, again, the secular terms associated with the resonant term  i in Eq. ( 30) are set to zero, which gives: 2i Finally, the following ansatz: with is imposed in Eq. (39), which gives the following solution for order  3 :

Soliton solutions for the wave envelope
The solution of Eqs. ( 35) and (40) resulting from conditions on the secular terms provide the envelope function  (  1 ,  2 ,  1 ,  2 ,  1 ,  2 ) which modulates the plane wave  i in Eq. (32).Seeking for wave solutions of  propagating with the group velocity v , we will transform the Eqs.( 35) and (40) by first representing space position -at the different orders   -in a fixed coordinate system {  ,   } with the first axis aligned with v (see Fig. 5) and subsequently using the Galilean transformation of group velocity to the coordinate   aligned with v In terms of the variables   ,   and   , Eq. ( 35) reduces to thus  does not scale with  in time.Similarly, Eq. ( 40) becomes Therefore, the solution of Eq. (47) provides the shape of the envelope  in the frame of reference {  ,   }, which propagates with the group velocity v .

Bright and dark solitons
We first seek for a solution of Eq. (47) that does not vary in the direction  1 , therefore with constant amplitude along the direction normal to the group velocity.It is worth noting that these characteristics will remain in the reference {, } since the change of coordinates given by Eq. ( 22) preserves angles.
Then, the propagation of the soliton reduces to a one-dimensional problem.This leads to the 1D nonlinear Schrödinger (NLS) equation [32]: where in this case After some algebra thus  is positive.
It is worth noting that Eq. ( 48) corresponds to the one found and solved by Nadkarni et al. [17] for a weakly nonlinear 1D lattice containing bistable elastic elements.The dispersive and the nonlinear terms in Eq. ( 48) can balance and yield to solitons.Depending on the sign of the parameter in the nonlinear term, the NLS equation becomes focusing for  > 0, or defocusing for  < 0 [33].The focusing equation possesses soliton solutions decaying to a nil background state, called bright solitons as they are localized traveling 'humps'.The defocusing equation admits soliton solutions on a nontrivial background, called dark as they are localized traveling 'dips'.
Following the approach proposed by Remoissenet [34] for the case  > 0 (focusing medium), and back-substituting   2 t, leads to the expression defining a bright soliton, which consists in a harmonic carrier wave  i modulated by an envelope given by: where the subscript '' denotes the 'bright' soliton type, and where s and ω are respectively given by: and Thus, the solution to the leading order can be expressed as where   =   is the wave amplitude.Eq. (54) describes a solitary wave defined by an envelope that propagates in the direction and speed defined by the group velocity.The envelope modulates a carrier wave defined by a wave vector  that is in general not aligned with the group velocity v , and that is related to the frequency  through the dispersion relations of the medium in the considered form.The frequency  of the carrier wave is also modified by a (small) factor ω that depends on the amplitude   squared.Additionally, in the long wavelength regime the group velocity is expected to be smaller than the phase velocity since, according to Eq. ( 21), the slope of the dispersion surface tends to zero when  → 0. Therefore the envelope will be slower than the carrier wave.According to Shabat and Zakharov [35] the bulk energy is contained in the soliton, which propagates with practically permanent shape leaving behind a dispersive tail.Moreover, any initial localized disturbance eventually evolves into a bright soliton, as also stated by Nadkarni et al. [17].Fig. 6 shows the displacement solution (3D view and cross section along the ȳ = x line at a given instant of time) provided by Eq. (54).This case illustrates a situation where the group velocity and the wave vector are not aligned, which occurs in the case when the lattice is not isotropic, i.e. when  ≠ 1.Although the solution process follows closely the steps for the 1D soliton predictions [17] the misalignment presented above predicts a situation which is unique to the anisotropic 2D lattice considered herein and has no 1D counterpart.The validity of the solution in Eq. ( 54) for non-isotropic lattices will be tested through numerical simulations as described in Section 4. In the case  < 0 (defocusing medium), the solution of Eq. ( 48) is [34]: where subscript '' denotes a 'dark' soliton, with The solution to the leading order becomes where   =   is the wave amplitude.This solution, originally found by Shabat and Zakharov [35], predicts a line of nil value of the amplitude for s − v t = 0, and thus travels at the group velocity.Since the sign of the hyperbolic tangent is an antisymmetric function at zero (odd function), its sign swapping forces a phase inversion in the carrier wave along the line s − v t = 0. Fig. 7 shows the displacement solution provided by Eq. (57), again for the case of an anisotropic lattice.Specifically, the cross section along the ȳ = x line at a given instant of time illustrates the anti-symmetric shape of the wave enforced by the dark soliton.

Vortex soliton
A third soliton type is investigated for the case of the isotropic lattice ( = 1, or  0 =  0 =  0 ).This limiting condition is imposed because the equation corresponding to the secular terms given in Eq. (47) reduces to which, upon the following change of variables can be recognized as 2D NLS equation of the third-order [36], of the form: In polar coordinates, Eq. ( 60) becomes: For  < 0 we get the defocusing 2D NLS equation [33] which supports dark vortex solitons of the form where subscript '' denotes a 'vortex' soliton.The signed integer  (taken here as ±1 for vortex stability) is called the topological charge, and   is the frequency of the vortex.Substituting the previous ansatz in Eq. ( 61) leads to a nonlinear Bessel ordinary differential equation in terms of the vortex amplitude function  (): 1 2 The sign of the topological charge  does not modify the previous differential equation, but defines the chirality of the vortex.The unknown function  is obtained by imposing as boundary conditions a nil amplitude at the origin  (0) = 0, and a constant background amplitude far from the vortex core,  ( → ∞) =  ∞ .This value is evaluated by assuming that the amplitude of the envelope is stationary at infinity.This corresponds to letting  ′ ,  ′′ → 0 for  → ∞, which gives a background amplitude equal to: thus   must have the same sign as  (negative) for an envelope solution consistent with the singular boundary condition.A closed-form analytical expression for the boundary-valued problem in Eq. (63) does not exist, but a numerical solution can be found by transforming it into an initial value problem at  = 0 and using the shooting method [32].It is worth pointing out that, for  = ±1 , the numerical solution can be suitably fitted with a hyperbolic tangent , where  0 is a parameter that fits the numerical solution of Eq. (63).
Finally, the solution to the leading order becomes: with ω =    2 , and   =   ∞ .The polar coordinates  and  can be transformed back into cartesian coordinates x and ȳ with the expressions (59) and ) Fig. 8 shows the displacement solution provided by Eq. ( 64), where the vortex center is located at the origin and propagates with the group velocity, along with the corresponding cross section along the line ȳ = 0.

Discrete lattice configurations and integration scheme
The analytical predictions presented in the previous section are validated through comparison with the numerical simulation of the response of the discrete lattice.The nonlinear equations of motion for a finite lattice consisting of 1000 × 1000 masses are assembled and subsequently integrated in time using the Verlet integration scheme [37].The is considered in a free-free configuration, so none of the masses are constrained.
The perturbation to the lattice is applied in the form of initial conditions that correspond the various soliton solutions.Such solutions are applied to the entire extent of the upon the application of a spatial window that modulates to zero the initial displacements and velocities at the boundaries.This allows the observation of the propagation within the lattice of the initial perturbation for a sufficient time before the wave interacts with boundary reflections.Thus, a direct comparison between analytical results and numerical predictions can be performed in the absence of boundary effects.In general, the applied initial conditions are expressed as: where ,  = 1, … , ,  = 1000.Also in Eq. ( 67) the subscript  identifies the type of soliton considered: ''-bright, ''-dark, or ''-vortex according to the solutions presented in the previous section.Moreover, the coordinates ( 0 ,  0 ) define the center of the soliton at  = 0, while  denotes the applied spatial window, which in this case consists in a 2D generalized Gaussian window centered at  0 =  0  +  0 , which is defined as: ) where   =    +    is the position of the  particle,   defines the width of the domain along the direction , while  is the number of particles along , which in this case is either perpendicular to the soliton or parallel to it.The resulting window bounds the initial conditions within a box of defined size.Finally,  is the parameter describing the generalized Gaussian window, with  = 2 providing the normal distribution, while  > 2 provides increasingly sharper edges.
Based on the analytical derivations presented in the prior section, solitons travel at the group velocity   , whose components in the long wavelength approximation can be expressed as: where  = || is the magnitude of the wave vector of direction .Accordingly, the soliton propagates in the plane in a direction , which is given by: This expression illustrates the potential misalignment of the wave vector and soliton direction of propagation that occurs when the lattice is non isotropic, i.e. for  ≠ 1.The difference in propagation direction of the soliton and wave vector are illustrated schematically in Fig. 9.Note that based on Eq. ( 69), misalignment between wave vector and soliton occurs only when  ≠ 0,  2 , which correspond to propagation along the  and  directions respectively.In these cases, the soliton is always aligned with the wave vector regardless the level of anisotropy given the fact that for the type of lattice considered herein, propagation along the principal lattice directions appears unaffected by properties in the direction perpendicular to them.

Bright and dark solitons
Results and comparisons are first presented for the bright and dark solitons described in Section 3.3.1.In all simulations in this section and in the remainder of the paper, results are found for unitary mass  = 1, and assume  = 1 as inter-mass distance.Also for both sets of simulations considered in this section, the parameters of the lattice are set such that  1 = 10 and   = 20, and the imposed wave vector has magnitude  =  10 and direction  =  4 .While these parameters are chosen somewhat arbitrarily, the choice of wave vector amplitude is motivated by the need to ensure that the simulations are conducted within the limits of the long wavelength approximation.Simulations are then conducted for a variety of values of the anisotropy parameter , which allows evaluating the solution in the presence of directional misalignment between wave vector and soliton propagation.The bright soliton simulations, presented first, are conducted for parameters defining the nonlinear ground potential  2 = 5,  3 = 5, which lead to  > 0. Results are first presented for the case of  = 1 with a soliton originating at particle located at coordinates  0 =  0 = 500 (where 500 denotes the particle number along the ,  directions), which illustrates the propagation along the prescribed 45 • line given the isotropy of the lattice.The considered Gaussian window limits the extent of the soliton in the direction perpendicular to its propagation only, given that the bright nature of the soliton naturally limits its spatial extent along the propagation direction.The Gaussian window is defined by = 2 and   = 200, and it is centered at the maximum of the soliton.Fig. 10 shows snapshots at the beginning of simulations and at  = 0.3  , where   denotes the last time instant computed by the time integration scheme.This simulation time is chosen based on a estimated time or arrival to the boundary of the domain.It is computed as   =   ∕(  cos ), with   denoting the distance from the initial location position of the soliton ( 0 ,  0 ) to the nearest boundary.Given the chosen set of parameters, it is here equal to   ≈ 390 s.These surface plots illustrate that the shape of the soliton is maintained during the propagation, which as expected occurs along the 45 • direction.A different representation of the results is given in Fig. 11, where the colored contours of the propagating wavefield show the spatial extent of the soliton, along with the equal phase contours, which in this case of  = 1 are perpendicular to the propagation direction (Fig. 11).Such direction is also represented as a thin black line superimposed to the wavefield contours.The interpolated cross-section of the wavefield along this propagation line allows the direct comparison between numerical predictions and analytical solutions which are shown in Fig. 11.Specifically, the blue line corresponds to the results of the numerical simulations, while the red envelope is obtained from the analytical solution of Eq. ( 51).This solution is propagated in time according to the soliton group velocity   , computed from Eq. ( 68).The comparisons in Fig. 11(b), (d) demonstrate that the analytical procedure presented in this paper correctly predicts the existence of bright solitary waves, and evaluates its amplitude and speed of propagation, which match the predictions of numerical simulations.
A second set of results is here presented for a lattice characterized by   ≠   , and specifically for  = 10.The simulation results and the comparisons with analytical predictions are presented in Fig. 12, using the same formats previously described.In this case, one can observe the existence of a misalignment between wave vector, directed perpendicularly to the lines of constant phase in the wavefield, and the direction of propagation, which is the result of the lattice anisotropy along its principal directions as defined by  = 2.While the wave vector is imposed to be still directed along the  = 45 • direction, the soliton propagation occurs at an angle  ≈ 84 • , as given by Eq. ( 69).This misalignment does not have a 1D counterpart.The results from numerical simulations, of which only a representative subset is presented herein, confirm that analytical predictions hold also for nonisotropic 2D cases.
Next, the case of the dark soliton is investigated.The dark soliton solution exists for  < 0, which we impose by selecting a nonlinear force coefficient  2 = 0.5.All other parameters in the simulations are kept the same.The windowing function in this case is applied in both perpendicular and parallel direction to the soliton, as the harmonic spatial content would otherwise cover the entire spatial domain and reflections would begin as soon as simulations are started.To this end, a Gaussian window with  = 8 and   = 400 is applied in addition to the transverse Gaussian window previously considered for the case of the bright soliton.The propagation corresponding to imposed initial displacements and velocities corresponding to a dark soliton (Eq.( 57)) is first illustrated in the form of snapshots surfaces, which are shown at the initial time instant in Fig. 13(a), and at about  = 0.3  in Fig. 13(b).The surface plots show the propagation of the dark soliton along the expected direction.
As in the case of the bright soliton, the contour representation of the dark soliton propagation highlights the wavefield, the soliton direction of propagation, and shows the line along which the numerical solutions are interpolated for the subsequent comparison with the analytical envelope predictions.The plots for the isotropic case ( = 1) are presented in Fig. 14, while the results for a non-isotropic lattice, now with  = 4, are displayed in Fig. 15.Both configurations demonstrate a good agreement between the numerical results and the analytical envelope predictions provided by Eq. (55), which again hold even when the soliton propagates at approximately  ≈ 76 • for a wave vector along the 45 • .

Vortex soliton results
As a final set of tests, numerical simulations were conducted for the case of the vortex soliton predicted from the analytical solution in Eq. (64).As in the previous examples, initial conditions for displacements and velocities were applied, along with a prescribed wave vector.Given that the obtained solution holds only for isotropic lattices, the direction of propagation is aligned with the soliton, and propagation occurs equally in all directions.Thus, the case of wave vector direction  = 0 is considered among all possible choices.The initial conditions were again windowed with the chosen generalized Gaussian function, which was set with  = 8 along the propagation direction, and  = 2 in the direction to it.This choice allows the propagation of the soliton undisturbed by the boundary reflections, while minimizing the distortions coming from the sides.The width of the windows was set to 200 particles in both directions to minimize distortion of the center of the soliton located at the center point of the domain in the initial time.
The results for two vortex solutions are presented in what follows, specifically for the cases of   = −20 and   = −10.In both instances, the analytical solutions expressed by Eq. ( 64) and the associated time derivatives are applied as initial conditions.The corresponding values of the fitting parameter  0 are provided as obtained from the application of the shooting method mentioned in the theory section.These values are  0 = 0.28397 and  0 = 0.401599 respectively for   = −20 and   = −10.As in previous examples, the long wavelength approximation is enforced by letting   = ∕10, with   = 0 for the considered propagation direction.The value of the topological charge was chosen to be  = −1.The case  = +1 was also evaluated, and is here not presented as the results do not appear qualitatively different.The relative value of the topological charge is instead of importance in the study of the interaction between multiple solitons, whose combination may either increase or reduce the stability of propagation of individual solitons.The results for   = −20 are presented in Figs.16 and 17, which show surface plots corresponding to snapshots at two instants of time, and the direct comparison with analytical predictions along the middle line parallel to the direction of propagation.Similar results for   = −10 are presented in Figs.18 and 19.The results show a good correlation between analytical envelope predictions and numerical simulations, both in terms of amplitude and speed of propagation.The numerical results show distortions of the solitons along the direction perpendicular to the propagation direction, along with the back propagation of a dispersive wave, which is particularly noticeable for the case of   = −10 (see Fig. 18).The generation of this back propagating wave, its relation with the stability of the single soliton, and the influence of its characteristic parameters are currently under investigation.However, the presented results illustrate the existence of this form of solitary wave, which has no 1D counterpart, and which is well predicted by the analytical framework developed in the paper.

Conclusions
The paper investigates the existence of solitary waves in twodimensional nonlinear lattices.The lattices are composed of linear springs connecting nearest neighbors, and of nonlinear springs connecting masses to the ground, leading to forces expressed as polynomials in terms of the out-of-plane displacement.Solitary wave solutions are found from the solution of continuum equations derived in the long wavelength limit.Application of the method of multiple scales, within the assumption of weak nonlinearities, leads to an ordered set of equations and to corresponding secular terms that express differential equations in terms of the envelope modulating the propagation of plane waves.These equations admit solitary wave solutions of three kinds, corresponding to bright, dark and vortex solitons.In the case of the bright and dark solitons, the solution approach follows previous procedures used for the derivation of solitary waves in one-dimensional lattices, but predicts 2D wave effects such as the misalignment between wave vector and group velocity.This occurs when the lattice is non isotropic.The case of the vortex soliton is also derived from the solution of a NLS obtained from the secular terms in the perturbation solution, which holds for isotropic lattices.The validity of the analytical solutions and the existence of the solitary waves is verified through numerical simulations of the discrete lattice response, where analytical solitons displacement and velocity are imposed as initial conditions.Comparison between numerical results and analytical solutions confirm the existence of the solitary waves for a variety of lattice parameters and propagation directions.The presented analysis and numerical results illustrate the behavior of solitary waves of various kinds in a prototypical nonlinear mechanical lattice, and describe effects which are unique to a 2D domain.Further investigations will extend the study to consider the behavior of multi-vortices and their interaction, the effect of anisotropy, and consider potential application of these configurations for the nondispersive transfer of signals and information.

Fig. 2 .
Fig. 2. Dispersion relations for the linear lattice (with  0 = 4.5) represented as iso-frequency contours: isotropic lattice ( = 1) (a), and anisotropic lattice ( = 10) (b).Color bars correspond to values of the non-dimensional frequency  .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Dimensionless dispersion diagram for the linear lattice (with  0 = 4.5) along the edges of the irreducible Brillouin zone (IBZ) for two different values of anisotropy:  = 1 (blue, solid) and  = 10 (orange, dashed).Schematics of the first Brillouin zone is depicted below.

)Fig. 4 .
Fig. 4. Comparison of linear dispersion relations: discrete system (solid line), and continuum approximation (dotted lines) truncated to the first order, see Eq. (21) (a) and to the second order, see Eq. (19) (b).Color bars correspond to values of the non-dimensional frequency .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Schematic of directions of soliton propagation and its relation with the direction of the wave vector .Lines of constant phase for the wave are shown as thin black lines, while the soliton is represented by the thick red line, and its direction of propagation is defined by the thick red arrow perpendicular to the soliton line.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11 .
Fig. 11.Contour of numerically predicted bright soliton wavefields for  = 1 and direction of propagation considered for comparison with analytical results (black thin line): time  = 0(a),  = 0.3  (c).Numerical integration results (blue line) and analytical envelope evaluated from Eq. (51) (red line): time  = 0(b),  = 0.3  (d).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .Fig. 13 .
Fig. 12. Contour of numerically predicted bright soliton wavefields for  = 10 and direction of propagation considered for comparison with analytical results (black thin line): time  = 0 (a),  = 0.4  (c).Numerical integration results (blue line) and analytical envelope evaluated from Eq. (51) (red line): time  = 0 (b),  = 0.4  (d).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 14 .
Fig. 14.Contour of numerically predicted dark soliton wavefields for  = 1 and direction of propagation considered for comparison with analytical results (black thin line): time  = 0 (a),  = 0.3  (c).Numerical integration results (blue line) and analytical envelope evaluated from Eq. (55) (red line): time  = 0 (b),  = 0.3  (d).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 15 .Fig. 16 .
Fig. 15.Contour of numerically predicted dark soliton wavefields for  = 4 and direction of propagation considered for comparison with analytical results (black thin line): time  = 0 (a),  = 0.4  (c).Numerical integration results (blue line) and analytical envelope evaluated from Eq. (55) (red line): time  = 0 (b),  = 0.4  (d).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)