Traveling waves in 2D hexagonal granular crystal lattices

This study describes the dynamic response of a two-dimensional hexagonal packing of uncompressed stainless steel spheres excited by localized impulsive loadings. The dynamics of the system are modeled using the Hertzian normal contact law. After the initial impact strikes the system, a characteristic wave structure emerges and continuously decays as it propagates through the lattice. Using an extension of the binary collision approximation for one-dimensional chains, we predict its decay rate, which compares well with numerical simulations and experimental data. While the hexagonal lattice does not support constant speed traveling waves, we provide scaling relations that characterize the directional power law decay of the wave velocity for various angles of impact. Lastly, we discuss the effects of weak disorder on the directional amplitude decay rates.


Introduction
Granular crystals can be defined as ordered closely packed arrays of elastically interacting particles. The propagation of traveling waves in homogeneous one-dimensional (1D) granular chains has been studied extensively. These investigations have been summarized e.g. in the reviews [1][2][3]. In 1D, granular chains support the formation of highly localized traveling waves, which can be described by relatively simple numerical and analytical approaches [1,4,5]; a variety of other states including periodic traveling waves/nonlinear normal modes and an effective particle description thereof has also been proposed [6,7]. The highly nonlinear response of these crystals stems from the Hertzian contact interaction between compressed spheres, and the absence of tensile response between the grains [8]. In this setting, the presence or absence of static precompression has been used to control the intrinsic wave propagation properties of impulsive stresses [1,9]. Besides the features and applications associated with Fermi-Pasta-Ulam lattices (such as the existence of supersonic traveling waves), granular crystals have been proposed for different applications including shock and energy absorbing layers [9][10][11][12], actuating devices [13], acoustic lenses [14], acoustic diodes [15] and sound scramblers [16,17].
Although one-dimensional granular chains have been studied at considerable length, homogeneous higherdimensional crystals have only recently started to be explored in more detail. Examples of the topics examined include: the dynamic load transfer path in disordered two-dimensional systems [18][19][20][21], the stress-wave anisotropy in centeredsquare nonlinear chains with different materials [22][23][24], and the redistribution of energy in a square lattice with one or more interstitials [25]. While a large number of studies have investigated the dynamic response of basic two-or threedimensional granular hexagonal lattices , the vast majority focus on the weakly nonlinear and linear (i.e. compressed) or vibrational lattice response. Although the theme of excitation propagation is of fundamental interest in its own right, it may also be of significant practical interest. For instance, in applications such as sound absorption, the goal is often to find conditions that promote energy dispersion (or dissipation) and it is thus of interest to examine such properties of different types of bead arrangements in higherdimensions. In that light, we are interested in the dynamic response of impacting an initially at rest hexagonal lattice.
In one spatial dimension, traveling solitary waves are generated upon striking one end of an initially at rest chain of particles [1,47]. The properties of these solitary waves depend on the composition of the chain, on the geometry and material properties of the particles and on the degree of static precompression applied to the chain. The two main analytical approximations developed in order to study such traveling waves are based on a quasi-continuum approximation [1,4], where a partial differential equation is derived from the pertinent lattice model, and a binary collision approximation (BCA). The BCA relies on the assumption that at a given time, a significant portion of the energy of the traveling structure is centered over two lattice sites where the resulting equations can then be solved exactly [47]. We should note in passing that in addition to these methods and especially in 1D, the Fourier space, co-traveling frame analysis of [5] enables a computation of the exact traveling wave up to a prescribed numerical tolerance. Recently, an additional analytical technique was developed that can capture the modifications (and even decay) of traveling waves in granular chains when they are subjected to on-site perturbations [48]. In 2D configurations, one can study traveling structures by observing a row of beads along different radial directions from the impact bead. For a square packing, it has been shown that quasione-dimensional traveling solitary waves emerge upon striking the lattice [49], and are essentially described by the onedimensional theory. In other arrangements, such as hexagonal packings, we argue that such quasi-one-dimensional motion is impossible, since each adjacent row and column will be affected upon being struck, regardless of the striking angle. After initially impacting a single bead, the energy will gradually be spread over progressively more and more lattice sites. Since the energy continues to spread to an increasing number of beads, the (energy and) velocity magnitude will gradually decrease, and thus a perfect traveling solitary wave will be impossible to support. This is in contrast to the 1D (resp. 2D square) situation, where the amplitude of the velocity profile remains almost constant, and hence, supports a traveling solitary wave. Thus, hexagonally packed granular crystals are far better candidates for applications such as sound absorption.
In this paper, we use numerical, theoretical, and experimental tools in order to study how energy is spread throughout a hexagonally packed granular crystal lattice upon being struck at different "strike angles" and being observed at different "observation angles". In addition to confirming the absence of a traveling solitary wave, we identify scaling power law patterns for the decay of the velocity as the wave travels through multiple layers of beads within the structure. In a simple special case, we are able to generalize the BCA into a ternary collision approximation (TCA), which yields a reasonable numerical (and approximate semi-analytical) estimate for the velocity decay. Lastly, additional numerical simulations incorporating weak disorder in the particle lattice are performed in order to more realistically simulate and compare with the experimental results.

Experimental setup
Experiments were performed on a 2D hexagonal array of particles consisting of 11 spheres along each edge of the lattice (see Fig. 1). The spheres were stainless steel (type 316, purchased from mcmastercarr.com), with a manufacturer specified diameter d = 19.05 mm and tolerance tol = 0.0127 mm. For subsequent calculations and numerical simulations, the particles were assumed to have a Young's modulus E = 193 GPa, Poisson's ratio ν = 0.3, and density ρ = 8,000 kg/m 3 [49]. The hexagonal particle lattice was confined, but not compressed, by walls on all six sides, with a hole at the impact location along one edge. Since we are concerned with only the primary traveling excitation, a comparatively "soft" material, delrin, was chosen to line the con-  [49]) in order to delay reflections from the system edges. To aid in the lattice particle alignment, a slight tilt (<5 • ) was imposed in the x-direction of the experimental setup.
To excite the system, a striker sphere, identical to the spheres composing the 2D hexagonal granular crystal, impacted the central particle along one edge. The striker sphere rolled through a channel down the inclined experimental setup and its initial velocity was calculated from the drop height. To measure the motion of individual particles within the system we used custom fabricated sensor particles that consisted of a miniature triaxial accelerometer (PCB 356A01 with sensitivity 0.51 mV/ m s 2 ) embedded in a stainless steel sphere. A detailed description of the sensor particles can be found in [49]. Based on the data acquisition system used (NI BNC-2110 and NI PCI-6123 with simultaneous sampling rate at 500 kS/s), the number of sensors present in the array during a single experimental impact was limited to four (with data acquisition in both the x-and ydirections). Therefore, to better capture the system response, three different sensor configurations were used: (1) along the 0 • observation direction, (2) perpendicular to and near the impact, and (3) perpendicular to and further from the impact (see Fig. 1). The sensor in the lattice location closest to the impact was present in all three sensor configurations.
In an ideal hexagonal packing composed of perfectly spherical particles, each sphere would have an initial coordination number of 6. In experiments, however, slight deviations in particle sizes can result in local compressive forces or gaps in the initial contact lattice, thus potentially reducing the initial or at rest coordination number. Previous studies on the effects of weak disorder in tightly packed 2D granular arrays showed that the exact system response depends on the initial contact lattice, and is generally consistent between repeated impacts on an individual configuration [23]. In order to capture the variability caused by differences in initial contact lattices, 20 different particle packings were assembled and experimentally tested (for each sensor configuration, 1-3) with an initial striker velocity of V striker = 0.40 m/s. Each of the 20 initial assemblies was impacted repeatedly to calculate average wave front amplitude and arrival time values, which could differ due to slight variations in impact conditions, such as exact alignment and speed of the striker particle, and minor particle rearrangements. Additionally, experiments for 20 different initial particle packings with sensor configuration 1 were performed for striker velocities, V striker = 0.25 m/s and V striker = 0.70 m/s.
The impact conditions of the experimental assembly were chosen for experimental ease and consistency. In the numerical simulations and theoretical considerations below, the initial velocity V 0 was assigned to the n = 0 bead, that is the bead impacted by the striker sphere in experiments (see also Fig. 2). Numerical simulations, comparing the system response for an array with V striker = V 0 and V n=0 = V 0 , showed that the difference was negligible for the studied system. The experimental data is presented together with the results from the numerical simulations and theoretical predictions.

Model
In order to model the above experimental setup, we consider an ordered homogeneous lattice of uncompressed spherical beads arranged in a hexagonal configuration, neglecting, at least for the purposes of considerations herein, dissipative effects. The Hamiltonian of such a system has the form, where q m,n ∈ R 2 is the displacement of the bead with index {m, n} from its initial position (see Fig. 2 where A is the Hertzian contact stiffness, which depends on the elastic properties of the material and the geometric characteristics of the beads [1]. For the uniform packing of The bracket notation, defined by [r ] + = max(0, r ), is used to capture the lack of tensile strength between particles. All parameters can be rescaled to unity using the transformation, where the variables with the tilde represent the solutions in the physical scaling (i.e. the dimensional ones). We define q m,n = [x m,n , y m,n ] T andq m,n = p m,n = [u m,n , v n,m ] T . In these variables, the equations of motion have the form, This system was solved numerically using a standard 4th order Runge-Kutta method. In the numerical simulations, we have assumed that the walls consist of spheres of infinite radius. The numerical simulation results for the configuration corresponding to the experimental setup are shown in Fig. 3. Additionally, numerical simulations were performed for a larger array consisting of 45-spheres along the bounding hexagonal edge, compared to the 11-sphere per edge system tested in experiments. Figure 4 shows the evolution of the wavefront shape over the larger 2D hexagonal system. It is interesting to observe therein, in addition to the continuous decrease of the wavefront speed, its gradual deformation in what appears to be a semi-circular shape. We will return to this point later in the text (in the discussion after presenting our detailed results). Although this model ignores dissipation and disorder in the system, it has been found to provide a good comparison with the experimental observations, as we also detail below. It should also be noted that we are not concerned here with effects associated with reflections from the boundaries within a finite chain, although the latter problem is potentially of interest in its own right.

Ternary collision approximation: a theoretical approach for 30 • angle of observation
The 2D hexagonal, nonlinear mass-spring system represented by Eq. (4) is extremely complex. Indeed, many the- oretical studies of two-dimensional lattice equations are for simpler, single component systems [50][51][52][53], although traveling waves in 2D Hamiltonian square lattices of springs have also been examined; see e.g. [54]. To help reduce the complexity of the system we consider a setting which is amenable to a semi-analytical description. More specifically, we turn to an analog of the BCA that has been developed for 1D chains [47,[55][56][57][58]. We first assume that the beads travel approximately along fixed angles upon impact, based on symmetry considerations. From a theoretical point of view, it is more straightforward to develop the theory for a striking and observation angle of 30 • since the pattern does not alternate along this line of observation, in contrast to what is the case for a 0 • striking, see Fig. 2. Considering the time scale where only the beads adjacent to the struck bead in the direction of propagation are affected results in a so-called ternary collision approximation (TCA). There are four beads involved in the distribution of the energy at each collisional step, yet we use the symmetry to reduce the number of degrees of freedom to three. In the renormalized system, we have the following TCA for the system of the right panel of Fig. 2, where x 0 is the displacement along the impact direction, x 1 is the displacement along that same direction of the bead adjacent to the impacted bead, and x a is the displacement of the bead adjacent to the impacted bead at the angle θ . By symmetry, the bead adjacent along the −θ line will have the same contribution as the x a bead. Now define z = x 0 − x 1 and y = cos(θ )x 0 − x a . We have thenz where we used θ = π/3, which corresponds to a hexagonal packing. Although the set of Eq. (5) is remarkably simpler than the original system (4), it only lends itself to an exact (analytically obtainable) solution in the case where y(t) = ωz(t), where ω ∈ R. However, note that this prescription requires also that the initial conditions satisfy the same scaling. The definition of z and y necessitate (from the corresponding initial conditions) that ω = 1/2. However, using the above ansatz in system (5) yieldsz = −(2 + ω 3/2 )[z] 3/2 + , as well as the "compatibility condition" 2 + ω 3/2 = (1/2 + 3/2ω 3/2 )/ω, whose sole real solution is ω ≈ 0.382. From the above, it is clear that the compatible solution is close to (although not exactly) satisfying the initial conditions. In that light, we will also use for comparison the exact analytical solution of the form: where E in the relevant hypergeometric function expression stands for the constant "energy" of the oscillation associated with z(t) and can be computed as E = v 2 /2, where v is the initial velocity of the bead x 0 . As explained in [47], the main purpose of the BCA is to offer an approximation of the velocity of the traveling wave that propagates through the chain. In order to compute this "pulse velocity" one needs to define the time T n the pulse resides on bead n. Let t = 0 be the time when the velocities of beads n − 1 and n become equal. Then T n is the time when the velocities of beads n and n + 1 become equal and the pulse velocity observed at the nth bead is c n = 1/T n .  In order to approximate the residence time T n , the BCA combines two points of view. First, it is assumed that the interaction is not instantaneous, such that the velocity of the impacted bead will gradually decrease. However, the BCA is only valid if two beads are involved in the dynamics, which is clearly not the case over the entire residence time T n . Thus, strictly in terms of the BCA, it is not clear at all how to initialize the next step in the iteration. To circumvent this issue, it is supposed that the interaction is instantaneous such that the conservation of kinetic energy and momentum can be used to compute the velocity that bead n will emerge with after the collision with bead n − 1. For example, for a heterogeneous 1D chain, these conservations yield V n = 2V n−1 /(1+m n /m n−1 ), where m n is the mass of bead n [47]. Applying this relationship recursively yields Thus, the reduced equations in the BCA setting at the nth step are initialized with z(0) = 0,ż(0) = V n . Although it appears somewhat awkward to adopt these two viewpoints in order to make the approximation work, the results in 1D chains have been surprisingly good [47,[55][56][57][58]. For this reason we carry out a similar procedure to compute the pulse velocity of the traveling wave in the hexagonal lattice. Due to the complexity of the hexagonal system there is no clear way to obtain an analog of relationship (7) using conservation of kinetic energy and momentum. However, in both the 1D BCA and 2D TCA setting, we can compute the time it takes for the initially impacted bead to obtain a velocity equal to its adjacent partner. For a homogeneous 1D chain, the BCA predicts this to occur for a velocity equal to 50 % of the initial velocity. Thus, to understand what contribution the additional beads in the hexagonal configuration will absorb, we use the reduced Eq. (5) to compute when the velocities of beads 0 and 1 are equal (see Fig. 2 for the labeling conven-tion). The ratio of the velocity at this time compared to the 1D case is one way to quantify the fraction of momentum that is absorbed by the additional beads. Call this value c. For example, suppose the velocity at the time of intersection was predicted to be 0.45 in the TCA and 0.5 in the BCA. Then c = 0.9, since the bead only reached 90 % of the value it would have in the absence of the additional beads. If in addition, we assume that this ratio remains constant as the traveling structure propagates through the lattice, we expect that bead n emerges with the velocity V n = cV n−1 after collision with bead n − 1. Applying this relationship recursively yields, where V 0 is the impact velocity. Equation (8) will now play the role of (7) for the hexagonal system. In the rescaled system we found c ≈ 0.8452. Given the repetition of the fundamental building block of the right panel of Fig. 2, we now simply solve the TCA system (5) with initial values defined through (8) and compute the intersection time. Doing this for several bead locations yields a discrete set of approximations for the pulse velocity observed at those bead locations. See the left panel of Fig. 5 for an example with V 0 = 0.1. There, it can be seen that such an approximation of this average pulse velocity favorably compares with the full numerical results for the hexagonal lattice.
To amend the TCA to a 0 • strike (but still observing along the 30 • line) we need to understand how the energy is transferred among the first two beads that are in contact with the impact bead (see Fig. 2). Due to symmetry considerations, it is reasonable to conjecture that the velocity contribution along the 30 • (resp. −30 • ) will be half of the impact velocity. This would be the case if kinetic energy and momentum were conserved. This relation was observed to be fairly accurate (see the right panel of Fig. 7). At the 0 • strike and 30 • observation angle we have experimental data available for In experiments, the arrival times were calculated using a specified threshold, approximately 10 % of the maximum pulse amplitude. We should also note that the experiments use the arrival time at the first sensor position in Fig. 1 to define t = 0. Thus, in order to compare the experimental values to the TCA and numerical simulations we estimate the arrival time at the first sensor using the numerical simulation.

A numerical and experimental study: velocity decay rate for variable striking/observation angle combinations
In applications, not only is the time of arrival of the traveling wave important, but the magnitude of its amplitude as well. The TCA, although simple enough to afford analytical approximations, ignores the (weak) dependence of the decay rate c on V 0 and is restricted to an observation angle of 30 • . We wish to have a more accurate and robust description of the velocity decay rate and thus we now turn to a (strictly) numerical study of other combinations of striking/observation angles. We corroborate the obtained qualitative and quantitative results by means of experimental observations.

30 • strike
In panel (a) of Fig. 6, the magnitude of the velocity of each bead along the 30 • line is shown against time. Notice how the maximum velocity obtained for each bead decreases as we move further from the impact point, as expected. We fit these maximum points with a function of the form V 1 c n−1 where c is a decay parameter that will now depend on the initial strike velocity V 0 (in contrast to the situation above which had c independent of V 0 ). We start the fit at V 1 since the initial impact will only have a velocity component, and thus the dynamics of the collision of the n = 0 and n = 1 bead are inherently different than the dynamics of the n and n + 1 beads for n > 0. In the latter case, the bead will have a velocity and position component. In panel (b) of Fig. 6 we fit the decay rates c = c(V 0 ) with a linear function mV 0 + b with two free parameters m, b. We found that m ≈ 0.257 and b ≈ 0.832. In panel (c) of Fig. 6 we obtain the relationship between V 0 and V 1 , which has the form V 1 = αV 2 0 + βV 0 . We found α ≈ 0.206 and β ≈ 0.629, such that the decay between V 0 and V 1 is greater than that between V n and V n+1 for n > 0. We have then that, where V n is maximum velocity of the bead lying n beads away from the strike point on the 30 • line. Clearly this descrip- tion of the velocity decay is more accurate than (8) since the dependence on the initial impact velocity is taken into account. Nevertheless, we still observe that a power law decay of the maximum velocity of each bead provides an accurate description of the dynamics along the 30 • observation angle in the hexagonal chain. When observing along the zero degree line it turns out that the fitting is more accurate when starting at bead 2, see Fig. 2 for the labeling convention. Panel (a) of Fig. 7 shows a relavent example where the velocity fitting starts at bead 2, but for a 0 • strike (which is analyzed in Sect. 5.2). Since the relationship between V 1 and V 0 was linear, it is natural to suppose that the relationship between V 2 and V 0 is quadratic, i.e. V 2 ≈ αV 2 0 + βV 0 + γ [see panel (d) of Fig. 6 for example]. We have then that, where V n is the magnitude of the maximum velocity of the nth bead from the strike point on the 0 • line. Our fitting yielded m ≈ 0.3039, b ≈ 0.8756 and α ≈ −0.4788, β ≈ −0.0256, γ ≈ 0.1936.

0 • strike
Finally, we carry out the procedure for the case of striking at a 0 • angle. As mentioned above for a zero degree observation, the first bead behaves somewhat differently [see panel (a) of Fig. 7], and therefore we start the fitting again at bead 2. We found that our scaling relations continue to be valid for a 0 • strike, but now with m ≈ 0.1766, b ≈ 0.8878 and α ≈ −0.2643, β ≈ −0.0969, γ ≈ 0.2152. The optimal scaling along the 0 • observation angle has the same form as (10).
When observing at a 30 • angle, the formula that yielded the best fit had the form, with m ≈ 0.1484, b ≈ 0.8393 and α ≈ −0.1182, β ≈ 0.0203, γ ≈ 0.5192. However, for the range of initial velocities used, the relationship between V 1 and V 0 is approxi- Fig. 7] which yields the simplified formula,

Preparation scheme
Numerical simulations incorporating weak disorder were also performed for the 2D hexagonal system described in Sect. 2. Similar to [23], weak disorder was incorporated in the simulations by assigning each particle in the array a ran-dom diameter, based on a normal distribution with mean μ = d = 19.05 mm and standard deviation σ = 6 tol, where tol is the diameter tolerance specified in Sect. 2.
The initial resting positions of the particles in the weakly disordered lattice were found in two steps. First, the particles were given an initial hexagonal lattice spacing assuming d = d + tol, to avoid large repulsive forces between overlapping particles. To bring the particles in contact, a 5N force was applied to each of the edge particles and an artificial damping (of the type introduced in [59]) was used to settle the random particle motion. Secondly, the wall positions were found based on the slightly compressed settled configuration. Then, gravity was introduced along the x-direction (in agreement with the experimental tilt), the 5N force was removed, and the particles were again settled through a similar, but weaker damping process. Once the initial positions were obtained, the settled array was impacted with a striker sphere, V striker = 0.25, 0.40, or 0.70 m/s, along one edge of the system. This process was performed for 20 different initial lattice configurations. For the subsequent discussion, we will use the terms "ideal" and "weakly disordered" to distinguish between numerical simulations where the spheres are all assigned the same diameter, and those just described with a variable diameter.

Discussion of results
Based on the numerical simulations of the ideal hexagonal lattice (Figs. 3, 4), the wavefront shape can be initially described by a hexagonal pattern which gradually transitions into a more circular shape. This transition occurs as the scale of the structure changes from being comparable to the lattice "spacing" to being much wider than that. It is worthwhile to mention that this is an intriguing feature that perhaps a quasicontinuum or homogenization type approach (see e.g. [60] and references therein for recent work on the subject) may capture and is left as an interesting open problem. Figure 8 compares the initial wavefront shape observed in the ideal numerical simulations with those observed in the experiments and numerical simulations incorporating weak disorder. The average wave front arrival times from both experiments and numerical simulations incorporating weak disorder clearly show the initial hexagonal wave front propagating though the structure. The average arrival times in experiments are slightly longer than those predicted from the ideal numerical simulations. Physically, this can be understood since the dissipation present in experiments results in a decreased wave amplitude, consequently decreasing the wave front speed. Additionally, imperfections in the experimental contact lattices also result in redirection of wave amplitude and delays in the signal arrival time [23]. As Fig. 8 indicates, the relative arrival times in the weakly disordered simulations are notably longer than both the predicted values from the ideal simulations as well as the experiments. Since the simulations do not incorporate dissipation, the time delays observed in simulations of the weakly disorder system are entirely due to the effects of the particle misalignments. The discrepancy between the arrival times in experiments and weakly disordered numerical simulations suggests that the tolerance values used to simulate imperfections in the contact lattice are fairly conservative. Figure 9 compares the wave front amplitude decay along the 0 • and −30 • observation directions after a 0 • strike angle for experiments, ideal numerical simulations, weakly disordered numerical simulations, and the corresponding Eqs. (10) and (11). Overall the trends observed for all approaches are in good agreement. The presence of weak disorder in the numer-ical simulations results in decreased amplitude transmission along the 0 • observation direction and increased amplitude transmission along the −30 • observation direction, compared to the ideal simulations. This altered distribution of wave front amplitude can be seen in Fig. 9, where the mean amplitude from simulations incorporating weak disorder falls above the ideal simulation values along the ±30 • observation direction and below the ideal simulation values along the 0 • observation direction. The 30 • observation direction represents a line of spheres directly in contact, while the wave front must travel though a zig-zag of particle contacts along the 0 0 direction. Therefore, the number of contacts, or possible amplitude scattering points, is greater along the 0 • observation direction, which could help physically explain this amplitude redistribution phenomenon. Although the presence of dissipation in experiments results in lower average amplitudes compared to the weakly disordered numerical simulations, which neglect dissipation, we also observe this trend in average experimental amplitude values (see the middle panel of Fig. 9 for V 0 = 0.40 m/s).

Conclusions and future challenges
We presented a systematic study of the dynamic response of a 2D hexagonal, highly nonlinear lattice excited by an impulse. In this 2D hexagonal setting, because of the ever-increasing number of neighbors over which the energy is distributed, no genuine traveling wave excitations, i.e. with constant velocity, have been found to persist. The propagating pulse energy has been found to decay as a power law, both in our numerical computations and in our experimental observations. Detailed expressions were provided to describe these power laws at different angles of strike and of observation. In a special case (of 30 • strike and 30 • observation, according to our presented classification), a generalization of the binary collision approximation (dubbed the ternary collision approximation) was presented and utilized to give simple numerical and even approximate analytical expressions for the bead evolution. Lastly, the effects of weak disorder on the propagating wave structure were examined; the average spatial amplitude values from numerical simulations incorporating weak disorder were in good agreement with experimental values as well as the corresponding fitted equations derived from numerical simulations on the ideal hexagonal granular crystal. This agreement reveals that the level of disorder present in experiments does not cause significant deviation of the propagating wave structure from the predicted system response.
While this is a fundamental step towards an improved understanding of non-square lattices, there are numerous additional tasks to be considered. While invariant traveling solutions may not exist, the possibility of existence of selfsimilar decaying solutions of the discrete model or perhaps of a quasi-continuum approximation thereof cannot be excluded and should be further considered.
Additionally, while here we have concerned ourselves with the highly nonlinear limit of the model (i.e. no precompression), there has recently been a surge of activity in connection to linear or weakly nonlinear properties of nonsquare (such as honeycomb or hexagonal) lattices. This is due to some of the remarkable linear properties of these lattices, such as, e.g. conical diffraction and the existence of Dirac (diabolical) points examined in both the physical [61,62] and mathematical [63] communities. It would be particularly relevant to consider the possibility of such spectral and nonlinear features in the precompressed variant of the present chain. Some of these possibilities are currently under investigation and will be reported in future publications.