---
abstract: 'Vector solitary waves are nonlinear waves of coupled polarizations that propagate with constant velocity and shape. In mechanics, they hold the potential to control locomotion, mitigate shocks and transfer information, among other functionalities. Recently, such elastic waves were numerically observed in compressible rubber-like laminates. Here, we conduct numerical experiments to characterize the possible vector solitary waves in these laminates, and expose a new type of waves whose amplitude and velocity oscillate periodically without dispersing in time. This oscillation is a manifestation of a periodic transfer of energy between the two wave polarizations, which we consider as internal mode of the solitary wave. We find that the vector solitary waves propagate faster at higher amplitudes, and determine a lower bound for their velocity. We describe a procedure for identifying which initial strains generate such vector solitary waves. This procedure also enables an additional classification between tensile and compressive solitary waves, according to the way that the axial strain changes as the waves propagate.'
author:
- |
    Ron Ziv and Gal Shmuel[^1]\
    [Faculty of Mechanical Engineering, TechnionIsrael Institute of Technology, Haifa 32000, Israel]{}\
title: Oscillating vector solitary waves in soft laminates
---

\#1[\#1]{} \#1[\#1]{} \#1[\#1\^]{} \#1[\_\#1]{} \#1[\#1\^]{} \#1[\#1\^[(m)]{}]{} \#1[\#1\^[(f)]{}]{} \#1[\#1\^[(p)]{}]{} \#1[\_\#1]{} \#1[\#1\^[(p)]{}]{} \#1[F\_[\#1]{}]{} \#1[C\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[b\_[\#1]{}]{} \#1[F\_[\#1]{}\^[-1]{}]{} \#1[\_[\#1]{}]{} \#1[d\_[\#1]{}]{} \#1[D\_[\#1]{}]{} \#1[e\_[\#1]{}]{} \#1[\#1\^]{} \#1[\#1]{} \#1\#2 \#1[I\_[\#1]{}]{} \#1 \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[h\_[\#1]{}]{} \#1[\#1\^]{} \#1[\#1\^[(p)]{}]{} \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[F\_[\#1]{}]{} \#1[C\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[b\_[\#1]{}]{} \#1[F\_[\#1]{}\^[-1]{}]{} \#1[\_[\#1]{}]{} \#1[d\_[\#1]{}]{} \#1[D\_[\#1]{}]{} \#1[e\_[\#1]{}]{} \#1[\#1\^]{} \#1[\#1]{} \#1\#2 \#1[(\#1)]{} \#1 \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[h\_[\#1]{}]{} \#1[\#1\^]{} \#1[\#1\^[(p)]{}]{} \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[\#1\^[(p)]{}]{} \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[\#1\^[(a)]{}]{} \#1[\#1\^[(b)]{}]{} \#1[\#1\^[(p)]{}]{} \#1[\#1\_[2,11]{}]{} \#1[\#1\_[2,1111]{}]{} \#1[\#1\_[2,1]{}]{} \#1[\#1\_[2,111]{}]{} \#1[\#1\_[,]{}]{} \#1[\#1\_[,]{}]{} \#1[\#1\_[,]{}]{} \#1[\#1\_[,]{}]{} \#1[\#1\_[(n)]{}]{} \#1[\^[(p)]{}]{} \#1[\#1\_[(n)m]{}\^[(p)]{}]{} \#1[u\_[\#1]{}]{} \#1[\#1]{} \#1\#2 \#1[\#1]{} \#1[I\_[\#1]{}]{} \#1 \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[\_[\#1]{}]{} \#1[h\_[\#1]{}]{} \#1[\#1\^]{} \#1[\#1\^[(p)]{}]{} \#1[\_[\#1]{}]{}\#1[\_[\#1]{}]{}\#1[\_[\#1]{}]{} \#1[P\_[\#1]{}]{} \#1\#2 \#1[I\_[\#1]{}]{} \#1[\#1]{} \#1[\#1]{} \#1[\#1\^[(a)]{}]{}\#1[\#1\^[(b)]{}]{}\#1[\#1\^[(n)]{}]{} \#1[\#1\^[(U)]{}]{}\#1[\#1\^[(L)]{}]{}

Introduction
============

Nonlinearities in wave mechanics are the source of fascinating phenomena, from instabilities [@Skipetrov200prl], harmonic generation [@Saltiel2008prl; @gonella17jmps; @khajehtourian2019nonlinear], extreme energy transfer [@ZHANG2018jmps] and non-reciprocal transmission [@Lepri2011prl], to the formation of shocks [@CHOCKALINGAM2020103746], rogue waves [@Baronio2012prl] and *solitary waves* [@SILLING2016jmps]. The latter are amplitude-dependent waves which propagate with a constant velocity and fixed shape, owing to a balance between dispersion and nonlinearity in the system [@dauxois2006physics]. Initially studied in fluid mechanics [@BOYD2015417], quantum mechanics [@Kasamatsu2006pra] and optics [@Stegeman1999science; @kivshar2003optical], solitary waves are recently gaining increased attention from the solid mechanics community [@nadkarni2014dynamics; @deng2018metamaterials; @DENG2020jmps; @Katz2019; @mo2019cnoidal]. This interest was triggered by the quest for mechanical *metamaterials*artificial composites with properties and functionalities not found in nature [@craster2012acoustic; @Christensen2015MRCComunications; @srivastava2015elastic; @bertoldi2017flexible; @Kadi2019nrp]; and its pertinent progress in additive manufacturing techniques for creating these architectured materials [@raney2015].

In addition to the mathematical significance of their analysis, the study of solitary waves in mechanical systems shows also technological potential in impact mitigation [@yasuda2019origami], mechanical logic gates [@raney16], nondestructive testing [@nasrollahi2017], wave focusing [@deng2019focusing], and robot locomotion [@Deng2020sciadv]. The framework of these works is almost exclusively of discrete models, with only few results on elastic *continua*. A specific class of continuum that exhibits nonlinearities is of soft materials, whose microscopic composition and ability to sustain large deformations [@ogden97book] are the source of constitutive and geometrical nonlinearities, respectively. Thus, soft materials constitute a barely explored platform to theoretically study solitary waves, and in turn experimentally realize using current 3D printing capabilities [@Bandyopadhyay2015MRS; @Truby2016; @GARCIA2019am].

To the best of our knowledge, the first numerical observation of one-dimensional solitary waves in elastic *periodic* continua was by @leveque2002finite and @bale2002siam, using a finite-volume method they developed, which was followed by the studies of @yong2003solitary [@andrianov2013dynamic; @andrianov2014numerical; @navarro2015nature; @hussein2018nonlinear] and the references therein. As explained in these works, the periodicity of the media causes dispersion, whose balance with geometrical and constitutive nonlinearities forms elastic solitary waves. Only recently, this principle for generating solitary waves was tested by @ziv2019b in *two-dimensional* continuum elastodynamics. There, the authors have developed a designated finite-volume method to simulate finite motions of soft laminates whose displacement field consists of two coupled components. By application of this method to a periodic repetition of compressible layers that differ in their mass density, the authors were able to observe the formation of *vector* solitary wavessolitary waves with (at least) two components and two polarizations that are coupled one with the otherin this case, when the axial and transverse components of the displacement field are coupled. To the best of our knowledge, that was the first report of vector solitary waves in an elastic continuum, within the framework of nonlinear elastodynamics [@ogden97book]. An earlier report of such waves in a *discrete* mechanical medium was given by @deng2017 [@deng2019prl], who conceived a model made of rigid squares that are connected using linear springs at their corners, thereby allowing for two coupled degrees of freedom for the squares: rotation and translation. While both the systems considered by @deng2017 [@deng2019prl] and the system considered here are two-dimensional and nonlinear, there are fundamental differences between them, as we explain in the sequel. These differences lead to opposite velocity-amplitude relation for the vector solitary waves in each system.

Here, we carry out a comprehensive study of the vector solitary waves that were observed by @ziv2019b, by analyzing a large set of numerical experiments that were generated using the method they developed in that work, with the following findings. First, we find that these vector solitary waves can be divided into two types which differ in their velocity and strain profile. We term the two types *quasi*-pressure and *quasi*-shear solitary waves, since in the limiting case when the heterogeneity vanishes and the strains are small, they reduce to the standard pressure and shear waves, respectively. We show that the profile of quasi-pressure vector solitary waves resembles the sech$^{2}$ function, which is the solitary wave solution of the well-known KdV equation [@korteweg1895xli]. We find that while quasi-pressure solitary waves maintain an identical wave profile as they pass between adjacent unit cells, the profile and velocity of quasi-shear solitary waves change between different cells in a periodic and permanent manner. This oscillation occurs through a periodic transfer of energy between the two polarizations of the vector solitary wave, which we consider as internal mode of the solitary wave. Similar long-lived oscillations were reported before in other fields [@Campbell1983ty; @Mantsyzov1995pra; @kivshar1998internal; @Szankowski2010prl], and were also reported recently in an acoustic system that is described using the continuum limit of a discrete model based on mechanicalelectrical analogies [@zhang2017bright]. Our results are the first report of oscillating vector solitary waves that are generated from the equations of continuum elastodynamics.

Next, we characterize the relation between the amplitude of the strain and the velocity of the solitary waves. This is done for a wide range of initial conditions, which were carefully chosen to generate in each experiment a single solitary wave. These experiments show that the velocity is a monotonically increasing function of the amplitude. This observation agrees with the results of @ziv2019b, who generated trains of solitary waves from a single set of initial conditions, and observed that the taller waves in the train propagate faster than the shorter waves. Notably, this velocity-amplitude relation is opposite to the relation for the vector solitary waves in the discrete mechanical model of @deng2017 [@deng2019prl], there the waves with smaller amplitudes are faster. We also determine the velocity of Bloch-Floquet waves in the laminate in the limit of low-frequency, long-wavelength linear elastodynamics [@santosa1991dispersive], to find they serve as a lower bound to the velocity of solitary waves in the nonlinear settings. Our results thus generalize the findings of @yong2003solitary and @andrianov2014numerical, who showed that one-dimensional solitary waves in nonlinear laminates are supersonic[^2], similarly to the feature of solitary wave solutions to the KdV equation [@dauxois2006physics].

Lastly, we describe a procedure for determining which initial deformations will generate vector solitary waves. This is done using contour maps of the characteristic velocities as functions of the strain, identifying the initial strain and loading path in these maps, and employing a certain criterion regarding the gradient of the characteristic velocity along the path, as explained in detail later. These maps are also useful in identifying how the axial strain in the laminate changes during the propagation of the solitary waves. Specifically, we identify two domains in the map for quasi-pressure solitary waves, namely, a domain of solitary waves that tend to increase the axial strain in the laminate and a separate domain of solitary waves that tend to decrease the axial strain in the laminate. In accordance with this tendency, we term the latter *compressive* solitary waves, notwithstanding the fact that the sign of axial strain in the laminate may be positive. Similarly, we term the former *tensile* solitary waves, notwithstanding the fact that the sign of axial strain may be negative. Interestingly, such a separation is absent from the map for quasi-shear solitary waves, since these waves can only decrease the axial strain, as explained later.

The rest of the paper is organized as follows. Sec.$\ $\[sec:Governing-equations\] contains a description of the relevant elastodynamics problem, together with its governing vectorial equations. Sec.$\ $\[sec:Method-of-solution\] first revisits solutions of the differential relations between the strain components in the benchmark problem of a homogeneous medium [@ZIV2019mom]. In the second part of this Sec., the computational study of the heterogeneous medium is carried out, where we also show that solutions for the homogeneous medium serve as an estimators for the differential relations in the main problem. A recap of our main results together with comments on future work close the paper in Sec.$\ $\[sec:Summary\].

Problem statement and governing equations\[sec:Governing-equations\]
====================================================================

We consider an infinite periodic repetition in the $X_{1}$ direction of two hyperelastic phases, namely, $a$ and $b$, governed by the same strain energy density function $\hyper$, and different initial mass density $\refden$. We set $X_{1}=0$ at the beginning of a certain $a$ phase, and denote its number by $n=0$, such that even and odd values of $n$ correspond to layers made of phase $a$ and $b$, respectively. At the initial state, the laminate is subjected to a combination of axial and transverse displacement fields as functions of $X_{1}$. Assuming that the laminate is uniform and infinite in $X_{2}$, the resultant motion $\Curx=\map\left(\Refx,t\right)$ is defined by

$$x_{1}=X_{1}+\dis 1\left(X_{1},t\right),\;x_{2}=X_{2}+\dis 2\left(X_{1},t\right),\label{eq:deformation}$$

and maps material points from the reference coordinate $\Refx$ to the current coordinate $\Curx$ (Fig.$\ $\[fig:Illustration\]). The problem amounts to determining the displacements $u_{i}\left(X_{1},t>0\right)$ for a given initial state. This is carried out by solving the balance equations for the corresponding first Piola-Kirchhoff stress tensor $\piola=\grad_{\defgrad}\hyper$, where $\defgrad=\grad_{\text{\ensuremath{\Refx}}}\map$ is the deformation gradient. When restricting attention to smooth waves[^3], the balance equations are

$$\Refdiv{\mathbf{\piola}}=\refden\deformation_{,tt},\label{eq:eom}$$

where we recall that $\refden$ is a periodic function of $X_{1}$, and it is tacitly assumed that $\map$ is continuous and twice differentiable, except at the interfaces between phases. For the given mapping , Eq.$\ $ is reduced to

$$\begin{aligned}\piolacomp{11,1}=\refden\derivativesec{u_{1}}t,\ \piolacomp{21,1}=\refden\derivativesec{u_{2}}t.\end{aligned}
\label{eq:motion}$$

In order to reformulate Eq.$\ $ in a more convenient form for the forthcoming analysis, we first introduce the Lagrangian strains and velocities of the axial and transverse displacements, namely, $$\pressurestrain:=\derivative{u_{1}}{X_{1}},\ \shearstrain:=\derivative{u_{2}}{X_{1}},\ \pressurevel:=\derivative{u_{1}}t,\ \shearvel:=\derivative{u_{2}}t.\label{eq:derivs}$$ In the limiting case of a simple shear $\left(u_{1}=0\right)$, the quantity $\shearstrain$ is called the amount of shear [@ogden97book], where in the limiting case of a uniaxial deformation $\left(u_{2}=0\right)$, the quantity $\pressurestrain$ is related to the axial stretch. We use Eq.$\ $ to carry out a reduction of order and rewrite Eq.$\ $ as $$\left(\begin{array}{c}
\pressurestrain\\
\shearstrain\\
\refden\pressurevel\\
\refden\shearvel
\end{array}\right)_{,t}+\left(\begin{array}{cccc}
0 & 0 & -\frac{1}{\refden} & 0\\
0 & 0 & 0 & -\frac{1}{\refden}\\
-\piolaonebyP & -\beta & 0 & 0\\
-\gamma & -\delta & 0 & 0
\end{array}\right)\left(\begin{array}{c}
\pressurestrain\\
\shearstrain\\
\refden\pressurevel\\
\refden\shearvel
\end{array}\right)_{,X_{1}}=\left(\begin{array}{c}
0\\
0\\
0\\
0
\end{array}\right),\label{eq:conservativeSystem}$$ where

$$\piolaonebyP=\derivative{\piolacomp{11}}{\pressurestrain},\;\piolaonebyQ=\derivative{\piolacomp{11}}{\shearstrain},\;\piolatwobyP=\derivative{\piolacomp{21}}{\pressurestrain},\;\piolatwobyQ=\derivative{\piolacomp{21}}{\shearstrain}.\label{eq:defineCoeff}$$

For hyperelastic materials we have that $\piolaonebyQ\equiv\piolatwobyP$ owing to the equality of mixed partials of $\hyper$.

The eigenvalues of Eq.$\ $ are the characteristic wave velocities at each material point, and are functions of the phase and deformation at that point, given by

$$c=\pm\sqrt{\frac{1}{2\refden}\left[\piolaonebyP+\piolatwobyQ\pm\sqrt{\left(\piolaonebyP-\piolatwobyQ\right)^{2}+4\piolaonebyQ\piolatwobyP}\right]}\eqqcolon\pm c_{\pm},\label{eq:wavespeeds}$$

where $c_{+}$ (resp.$\ $$c_{-}$) corresponds to the plus (resp.$\ $minus) sign of the inner square root. In the limit of linear elasticity, $c_{+}$ and $c_{-}$ are the velocities of pressure and shear waves, respectively [@DAVISON1966249]. Since in finite elasticity the corresponding modes involve both axial and transverse displacements, *i.e.*, there are no pure modes, we refer to those associated with $\shearwavespeed$ and $\pressurewavespeed$ as *quasi-*shear and *quasi-*pressure waves, respectively [@bland1965plane; @ZIV2019mom].

The eigenvectors of the matrix in Eq.$\ $ which are associated with the generated waves are

$$r_{1}=\left(\begin{array}{c}
1\\
\eta_{+}\\
\refden c_{+}\\
\eta_{+}\refden c_{+}
\end{array}\right),\;r_{2}=\left(\begin{array}{c}
-\eta_{-}\\
1\\
-\eta_{-}\refden c_{-}\\
\refden c_{-}
\end{array}\right),\ r_{3}=\left(\begin{array}{c}
-\eta_{-}\\
1\\
\eta_{-}\refden c_{-}\\
-\refden c_{-}
\end{array}\right),\;r_{4}=\left(\begin{array}{c}
1\\
\eta_{+}\\
-\refden c_{+}\\
-\eta_{+}\refden c_{+}
\end{array}\right),\label{eq:eigenvectors}$$

where $$\eta_{+}=\frac{\refden c_{+}^{2}-\piolaonebyP}{\beta},\quad\eta_{-}=\frac{\beta}{\alpha-\refden c_{-}^{2}},\label{eq:eta}$$ and the corresponding characteristic velocities are $-\pressurewavespeed,\,-\shearwavespeed,\,\shearwavespeed$ and $\pressurewavespeed$, respectively. Through algebraic manipulation and using the fact that $\alpha,\beta,$ and $\delta$ are derived from the same scalar potential, we find that $\text{\ensuremath{\eta_{+}}}$ is identical to $\eta_{-}$ and independent of $\rho_{L}$. The components of the eigenvectorsand specifically $\eta_{\pm}$are differential relations between the strain and velocity fields. When the initial conditions are piecewise-constant with a single jump discontinuity and the medium is homogeneous, the problem is called a Riemann problem, which can be solved by a variation of the method of characteristics. The method exploits the fact that the solution can be expressed as function of a single variablehere $X_{1}/t$in which case the corresponding waves are called simple waves [@davison2008book]. The resultant equations can be solved using standard methods such as RungeKutta methods. @ZIV2019mom have adapted such an approach to calculate smooth waves and shocks in semi-infinite compressible Gent media, when subjected to combined transverse and axial impacts. The challenge in solving Eq.$\ $ when the medium is heterogeneous and the initial data is not piecewise-uniform is associated with the interactions between waves that are repeatedly scattered at the interfaces. In this case, the ansatz of simple waves fails and the previous approach is no longer applicable. This brings a need for more sophisticated numerical solvers. The course taken by @ziv2019b was to develop a designated algorithm based on a finite-volume method[^4], by adapting the schemes of @LEVEQUE1997327 [@leveque2002finite] and @bale2002siam, to the type of problem considered here; this adaptation is used in the study described next.

Computational study\[sec:Method-of-solution\]
=============================================

*Homogeneous benchmark problem*.—Before we proceed to the study of the main problem, it is advantageous to present the strain relations of the homogeneous benchmark problem. As we will show in the sequel, the reason is that they serve as estimators to the laminated case. Thus, the differential relation between the axial and shear strains in quasi-shear and quasi-pressure waves traversing homogeneous media are given by [@DAVISON1966249; @ZIV2019mom]

\[integralCurves\] $$\begin{aligned}\textrm{quasi-pressure:\ \ } & \derivative{\shearstrain}{\pressurestrain}=\eta_{+},\end{aligned}
\label{eq:integralCurvesPressure}$$ \[etea\] $$\begin{aligned}\textrm{\;\;\;\;\ \;\;\;\,quasi-shear:\ \ } & \derivative{\shearstrain}{\pressurestrain}=-\eta_{-}^{-1},\end{aligned}
\label{eq:integralCurvesShear}$$

with the compatibility condition $$\shearstrain\left(\pressurestrain=\initial{\pressurestrain}\right)=\initial{\shearstrain},\label{eq:comp0}$$ where $\initial{\pressurestrain}$ and $\initial{\shearstrain}$ can be associated with axial and shear pre-strains, respectively. Eq.$\ $ is intimately related to the Poynting effect [@Poynting1909prsa], namely, the axial deformation that accompanies finite shear, hence occurs only at nonlinear deformations [@ogden97book]. Thus, this effect is prominent in soft materials, which are capable of undergoing large strains [@Cioroianu2013PRE; @Horgan2017prsa]. The exact manner in which this effect takes place depends on the constitutive response of the material [@Mihai2017nr]. Accordingly, the solution of Eq.$\ $ also depends on the constitutive response of the material, which in our study is modeled using the compressible model with the strain energy function [@olp2007] $$\hyper\left(\defgrad\right)=-\frac{\mu\jm}{2}\ln\left(1-\frac{\tr{\defgradT\defgrad}-3}{\jm}\right)-\mu\ln\det{\defgrad}+\left(\frac{\bulk}{2}-\frac{\mu}{3}-\frac{\mu}{\jm}\right)\left(\det{\defgrad}-1\right)^{2};\label{eq:Gent}$$ here, $\mu$ and $\kappa$ coincide with the shear and bulk moduli in the limit of small strains, respectively, and the dimensionless parameter $\jm$ models a stiffening of strain. The Gent model was originally developed to capture the stiffening of rubber materials . Recently, the model was also shown suitable for capturing nonlinear wave phenomena such as shear shocks and tensile-induced shocks, which were experimentally observed by @Catheline2003PRL and @espindola2017shear, and @NIEMCZURA2011442, respectively.

The first Piola-Kirchhoff stress tensor derived from Eq.$\ $ is $$\piola=\mu\left(1-\frac{\tr{\defgradT\defgrad}-3}{\jm}\right)^{-1}\defgrad+\left(\kappa-\frac{2\mu}{\jm}-\frac{2\mu}{3}\right)\det{\defgrad}\left(\det{\defgrad}-1\right)\defgradTi-\mu\defgradTi,\label{eq:gentstress}$$ and the corresponding components required for calculating Eq.$\ $ are $$\begin{aligned}\piolacomp{11} & =\mu\left(\pressurestrain+1\right)\left(1-\frac{\pressurestrain^{2}+\shearstrain^{2}+2\pressurestrain}{\jm}\right)^{-1}+\left(\kappa-\frac{2\mu}{\jm}-\frac{2\mu}{3}\right)\pressurestrain-\frac{\mu}{\pressurestrain+1},\\
\piolacomp{21} & =\mu\shearstrain\left(1-\frac{\pressurestrain^{2}+\shearstrain^{2}+2\pressurestrain}{\jm}\right)^{-1}.
\end{aligned}
\label{eq:pexplicit-1}$$ The resultant characteristic velocities are $$c_{\pm}^{2}=\frac{\vc_{3}}{2\refden}+\frac{\mu}{2\refden\vc_{1}^{2}}+\frac{\left(\vc_{1}^{2}+\shearstrain^{2}+\jm\right)\mu\jm}{2\refden\vc_{2}^{2}}\pm\frac{\sqrt{48\vc_{1}^{6}\mu^{2}\shearstrain^{2}J_{m}^{2}+\left(\vc_{1}^{2}\vc_{2}^{2}\vc_{3}+2\vc_{1}^{4}\mu J_{m}-2\vc_{1}^{2}\shearstrain^{2}\mu\jm+\vc_{2}^{2}\mu\right){}^{2}}}{2\sqrt{3}\refden\vc_{2}^{2}\vc_{1}^{2}},\label{eq:cmgent-1}$$ where $$\vc_{1}=\pressurestrain+1,\ \vc_{2}=\pressurestrain^{2}+2\pressurestrain+\shearstrain^{2}-\jm,\ \vc_{3}=\kappa-\frac{2\mu}{3}-\frac{2\mu}{\jm}.\label{eq:cg comp}$$ In the calculations to follow, we use representative moduli of rubber (see [@Marckmann2006yu; @Getz2017gsm], and the references therein), namely, $$\kappa=1\,\mathrm{MPa},\ \mu=200\,\mathrm{kPa},\ \jm=10.\label{eq:paremeters}$$ Numerical solutions of Eqs.$\ $- are shown in Fig.$\ $\[fig:integralCurves\]. We clarify that these solutions are not associated with a concrete problem, and are merely solutions of nonlinear differential equations with compatibility conditions. The resultant curves can be identified with the evolution of $\shearstrain$ and $\pressurestrain$ during the propagation of smooth waves which spread in a uniformly pre-strained homogeneous Gent medium, when subjected to impact [@ZIV2019mom]. Here, we present the shear strain $\shearstrain$ as function of the axial strain $\pressurestrain$, associated with the quasi-pressure (panel a) and quasi-shear (panel b) waves. Each curve is obtained using different values for the compatibility condition , which can be identified by the intersection of a certain curve with the vertical (panel a) and horizontal (panel b) axis. Accordingly, each curve in panel (a) corresponds to one of the conditions $$\shearstrain\left(\pressurestrain=0\right)=1.7,1.8,...,2.5,\label{eq:specificcomp}$$ where we note that the curves associated with the values $2,2.2$ and $2.4$ will we used in the sequel. Similarly, each curve in panel (b) corresponds to one of the conditions $$\shearstrain\left(\pressurestrain=-0.4,-0.2,...,0.6\right)=0.\label{eq:specificcomp-1}$$ We observe that for quasi-pressure waves $\shearstrain$ is a monotonically increasing function of $\pressurestrain$, while for quasi-shear waves $\shearstrain$ is a monotonically decreasing function of $\pressurestrain$. *It is important to note that the solutions to Eq.$\ $ are independent of the initial mass density $\refden$*. This property will be useful when solving the main problem, where the initial mass density of the medium varies in space.

$\ \ \ \ \ $

\
*Main problem: laminated media*.—Having at hand the solutions of the differential equations which are associated with a homogeneous medium, we return to the main problem of the laminated media. Specifically, consider a laminate whose layers are governed by Eqs.$\ $-, have an equal thickness of $1\,$cm, and the initial mass density of phases $a$ and $b$ is $500\,\mathrm{kg/m^{3}}$ and $4000\,\mathrm{kg/m^{3}}$, respectively. We specifically seek solutions to the response of the laminate when subjected to non-uniform initial strains. The initial strains we consider are in the form of a sum of a uniform part and a localized part about $X_{1}=0$, namely,

$$\epsilon_{i}\left(X_{1},t=0\right)=\begin{cases}
\initial{\epsilon_{i}}+\amp{\epsilon_{i}}\cos\frac{\pi X_{1}}{w}, & -\frac{w}{2}<X_{1}<\frac{w}{2},\\
\initial{\epsilon_{i}}, & \mathrm{otherwise},
\end{cases}\label{eq:initialState}$$

where $i=A,\,T$ denotes the axial and shear strains, respectively, and superscripts $\left(U\right)$ and $\left(L\right)$ correspond to the uniform and localized parts of the pre-strain, respectively.

\

As mentioned in Sec.$\ $\[sec:Governing-equations\], to solve the equations of motion we apply the finite-volume method of @ziv2019b, which was developed for analyzing motions comprising two coupled displacements in soft laminates.

We begin by showing in Fig.$\ $\[fig:pressureSoliton\](a) the distribution of the axial (left panel) and shear (right panel) strains as functions of space and time, when the laminate is subjected to a pre-strain with the parameters

$$\initial{\shearstrain}=2.2,\ \initial{\pressurestrain}=0,\ \amp{\shearstrain}=0.15,\ \amp{\pressurestrain}=0.3,\ w=0.08\,\mathrm{m}.\label{eq:quasi-pressure-strains}$$

We observe that a rightward propagating wave is generated. A movie of this propagation is also provided in the supplementary material online under the name S1. Remarkably, while the wave profile is different in phases $a$ and $b$, as it should, the respective profile in each phase is maintained in passing from one unit cell to another. Furthermore, by calculating the time of flight that takes the maximum of the strains to pass from a certain phase in one periodic cell to the same phase in adjacent cell, we find that its velocity is constant and equal for both phases; this wave is therefore a solitary wave. More specifically, since the solitary wave exhibits coupling between its axial and transverse components and polarizations, it is a *vector* solitary wave. To the best of our knowledge, these vector solitary waves were first observed by @ziv2019b using numerical experiments. Here, we carry out a comprehensive characterization of the possible vector solitary waves in the medium.

We first observe that the wave couples axial and shear strains of an identical sign, similarly to quasi-pressure waves in a homogeneous compressible Gent medium (Fig.$\ $\[fig:integralCurves\]a). Hence, we associate the solitary wave with a quasi-pressure mode, and term it a quasi-pressure solitary wave. To facilitate the identification of these features, we show in Fig.$\ $\[fig:pressureSoliton\](b) the corresponding axial (left panel) and shear (right panel) strains at the center of selected layers made of phase $a$, as functions of time. Arbitrarily, we chose each $6^{\mathrm{th}}$ layer from the range $64\leq n\leq88$; each curve corresponds to a different layer whose number appears on top of the curve. Evidently, the shape is maintained, and the equal spacing between the curves confirms that the velocity is constant too. Furthermore, we observe that the wave profile is reminiscent of the sech$^{2}$ profile, which, as mentioned, is the solitary wave profile associated with the well-known KdV equation [@korteweg1895xli]. This wave profile was also obtained for one-dimensional solitary waves in a laminate, using an asymptotic homogenization method [@andrianov2014numerical]. A comparison between the wave profile at $n=76$ and the sech$^{2}$ functions $$\pressurestrain=0.154\mathrm{sech}^{2}\left(1.33t\right),\ \shearstrain=2.2+0.108\mathrm{sech^{2}}\left(1.33t\right)\label{eq:sech-1}$$ (plotted in red dash), demonstrates the agreement between them.

\

Fig.$\ $\[fig:SolitonShear\](a) shows the distribution of the axial (left panel) and shear (right panel) strains as functions of space and time, when the laminate is pre-strained according to the parameters

$$\initial{\shearstrain}=2.2,\ \initial{\pressurestrain}=0,\ \amp{\shearstrain}=0.5,\ \amp{\pressurestrain}=-0.5,\ w=0.08\,\mathrm{m}.\label{eq:quasi-shear-strains}$$

We observe that a rightward propagating wave is generated, as also depicted in movie S2 in the supplementary material, with the following features. First, by examining the wave profile, we find that in contrast with quasi-pressure waves, here (*i*) the time dependency of the waves across each layer is not a sech$^{2}$ functionnor any other symmetric function; (*ii*) the wave profile changes not only between the different phases, but also between adjacent unit cells. The latter change in the profile occurs in a periodic manner without dispersing in time.

Second, by calculating the time of flight that takes the maximum of the strains to pass from a certain phase in one periodic cell to the same phase in adjacent cell, we find that the wave velocity oscillates about a constant value. These persistent oscillations of the amplitude and velocity occur through a periodic transfer of energy between the axial and transverse strains of the wave, which we consider as internal mode of the solitary wave. Similar long-lived oscillations were reported before in other fields [@Campbell1983ty; @Mantsyzov1995pra; @kivshar1998internal; @Szankowski2010prl], and recently were also reported in an acoustic system using the continuum limit of a discrete model based on mechanicalelectrical analogies [@zhang2017bright]. Our results are the first report of oscillating vector solitary waves that are generated from the equations of continuum elastodynamics.

Lastly, the axial and shear strains associated with the wave have an opposite sign, similarly to quasi-shear waves in a homogeneous compressible Gent medium (Fig.$\ $\[fig:integralCurves\]b). On account of these features, we term the wave *oscillating quasi-shear vector solitary wave*[^5]*.*

In order to characterize the dependency of the wave velocity and the relation between $\shearstrain$ and $\pressurestrain$ on the initial strain, we perform a set of numerical experiments with different initial states. Specifically, we set $\initial{\pressurestrain}=0$, consider three values of $\initial{\shearstrain}$, namely, 2, 2.2 and 2.4 (denoted in the following figures by blue, red and green markers, respectively), and for each value of $\initial{\shearstrain}$ we carry out calculations considering 6 combinations of $\amp{\pressurestrain},\ \amp{\shearstrain}$ and $w$, chosen from the ranges

$$-0.8<\amp{\pressurestrain}<0.8,\ \ \ 0.1<\amp{\shearstrain}<0.7,\ \ \ 0.07\,\mathrm{m}<w<0.08\,\mathrm{m}.$$

These combinations were selected by trial and error to yield a single vector solitary wave—either quasi-shear or quasi-pressure. By contrast, the initial conditions used by @ziv2019b delivered trains of solitary waves. The combinations for quasi-pressure solitary waves are ordered such that when we increase $\amp{\shearstrain}$, we also increase $\amp{\pressurestrain}$, where $\amp{\pressurestrain}$ is always positive. Conversely, the combinations for quasi-shear solitary waves are ordered such that when we increase $\amp{\shearstrain}$, we also decrease $\amp{\pressurestrain}$, where $\amp{\pressurestrain}$ is always negative. In both cases, the value of $w$ is decreased when $\amp{\shearstrain}$ is increased.

Our results for quasi-pressure solitary waves are shown in Fig.$\ $\[fig:Characterization\](a), in the following manner. From each simulation we extract the maximal state of strain over the time that it takes the generated wave to cross a unit cell. We specifically extract the maximal strain state at the middle of phases $a$ and $b$, respectively. The corresponding states are denoted in the $\left(\shearstrain,\pressurestrain\right)$ plane by disc (phase $a$) and diamond (phase $b$) markers.

In addition to the results from the full-wave simulations, we also present in Fig.$\ $\[fig:Characterization\](a) the solutions of Eq.$\ $ for different compatibility conditions , as previously depicted in Fig.$\ $\[fig:integralCurves\]. While these solutions are not associated with a concrete problem, they capture the way in which $\shearstrain$ and $\pressurestrain$ evolve during the propagation of smooth waves which spread in a uniformly pre-strained homogeneous Gent half-space. We specifically present using black curves the solutions to compatibility conditions that agree with the uniform part of the initial strain in the laminate ($\amp{\pressurestrain}=\amp{\shearstrain}=0$), namely, $$\shearstrain\left(\pressurestrain=0\right)=\initial{\shearstrain},\quad\initial{\shearstrain}=2,2.2,2.4.\label{eq:initial conditions}$$ We observe that the results extracted from the full-wave simulations of the laminate either coincide or fall very close to the curves that are associated with the impact response of the homogeneous medium. Specifically, we observe that the results for layers made of phase $a$ coincide with these curves for all values of $\initial{\shearstrain}$. We find it remarkable that these two seemingly unrelated problems are related in this way. This agreement between the full-wave simulations and the solution of Eq.$\ $ indicates that the analysis of the homogeneous medium can be used to predict quasi-pressure solitary waves in a soft laminate with a periodic modulation in its mass density. This result relies on the fact that Eq.$\ $ is independent of the mass density, and hence turns useful in the solution of a laminate whose only modulation is in that property.

Fig.$\ $\[fig:Characterization\](b) mirrors the presentation in panel a for quasi-shear solitary waves. Since in this case the profile changes between different cells, we present our results in the following way. First, we have calculated the maximal values of $\shearstrain$ and $\pressurestrain$ over the time that the solitary wave traverses a unit cell, as we previously did for the quasi-pressure waves. Since now different layers exhibit different maximal values owing to the internal mode, we here present the average value. Specifically, the average of the values for phase $a$ (resp.$\ $$b$) is taken over the number of layers within one period of the oscillation; the average is denoted by the disc (resp.$\ $diamond) markers. As in panel (a), to the results from the full-wave simulations we add in black curves the solutions of Eq.$\ $, when subjected to the compatibility conditions . Here again, the markers associated to the laminate either coincide or fall very close to these curves.

To highlight the oscillatory nature of the wave amplitude, we present in the inset of Fig.$\ $\[fig:Characterization\](b) the values over which the average was obtained for two exemplary markers; each point in the inset corresponds to a different layer that was used in the calculation of the average, where the number of three arbitrary layers is displayed. The slow spatial modulation in the wave profile that appears in Fig.$\ $\[fig:SolitonShear\] manifests itself here as closed ellipses, showing how the amplitude oscillates periodically about a mean. Interestingly, the orientations of the ellipses associated with each phase are opposite one to another.

The dependency of the velocity on the amplitude is shown in Fig.$\ $\[fig:velocities\](a) by plotting the velocity of quasi-pressure solitary waves as function of the maximal value of $\pressurestrain$ in phase $a$ over the time that the solitary wave traverses a unit cell. Here again, the cyan, red and green markers correspond to waves that were generated with $\initial{\shearstrain}=2,2.2$ and $2.4$, respectively. We observe that the highest and lowest clusters correspond to the largest and smallest $\initial{\shearstrain}$, respectively. We furthermore observe that within each cluster, markers at larger values of $\epsilon_{A}$ exhibit faster velocities. Accordingly, we conclude that the velocity is a monotonically increasing function of $\shearstrain$ and $\pressurestrain$. This is in agreement with the observation of @ziv2019b who generated trains of solitary waves, and found that the taller waves in the train propagate faster than shorter waves. This dependency is the opposite of the dependency of the vector solitary waves observed in the mechanical system conceived by @deng2017, as the two problems are fundamentally different. Notably, our system is a continuum solid which is constitutively and geometrically nonlinear, while the system of @deng2017 is captured by a model of linear springs and rigid squares. Moreover, the physical properties of our continuum are heterogeneous, as its initial mass density varies in space, while in the system of @deng2017 all the springs and squares share the same properties.

We show next that it is possible to bound from below the velocities $c_{+}$ and $c_{-}$ using the velocities in the linear limit (Bloch-Floquet waves) at the low-frequency, long-wavelength regime. Before we proceed, we note that @yong2003solitary have numerically found that the velocity of Bloch-Floquet compression waves in that regime is a lower bound for the velocity of one-dimensional compression waves in the nonlinear laminates they studied. @andrianov2014numerical later showed this result analytically for a different nonlinear model, using an asymptotic homogenization method [@andrianov2013dynamic]. Here, we generalize their conclusion to the present vectorial problem, where there are two coupled components of the displacement field and two coupled wave polarizations, using the limiting velocities of the quasi-shear and quasi-pressure waves. To this end, we follow @santosa1991dispersive, who rigorously showed that the velocity of the leading term in the Bloch expansion is simply the root square of the harmonic average of the stiffness over the weighted average of the mass density, and further showed that this result agrees with multi-scale asymptotics or homogenization. To carry out this calculation in the linear limit of our problem, we evaluate the instantaneous stiffness of the phases in the pre-strained configuration [@ogden97book; @LUSTIG2018jmps], and use the initial mass density. We denote the velocities that are calculated using the resultant average for the quasi-pressure and quasi-shear waves by $c_{+}^{\left(LL\right)}$ and $c_{-}^{\left(LL\right)}$, respectively. These values at each pre-strain are given by the dashed lines in Fig.$\ $\[fig:velocities\], confirming they are lower bounds. This observation implies that quasi-pressure and quasi-shear vector solitary waves are supersonic, thereby generalizing the results of @yong2003solitary and @andrianov2014numerical in the scalar case.

Fig.$\ $\[fig:velocities\](b) shows the velocity of quasi-shear solitary waves as function of the maximal values of $\left|\pressurestrain\right|$ (recall that here $\pressurestrain<0$) in phase $a$ over the time that the solitary wave traverses a unit cell, averaged over the number of layers within one period of the modulation. Here again, the velocity of quasi-shear solitary waves is a monotonically increasing function of the strain amplitude, and is bounded from below by $c_{-}^{\left(LL\right)}$.

We recall that for solitary waves to propagate there should be a balance between dispersion and nonlinearity. In our case of a finitely deformed laminate, this balance is achievable by pre-strain tuning: solitary waves can propagate only at specific initial deformations for which the material stiffens along the loading path [@hussein2018nonlinear]. The quantification of this stiffening is given by the gradient of the characteristic velocity along this path. Specifically in our case, this condition is satisfied when each characteristic velocity, *i.e.*, $\phaseA{c_{\pm}}$ and $\phaseB{c_{\pm}}$, is increasing monotonically along the loading path given in Eqs.$\ $ and for quasi-pressure and quasi-shear solitary waves, respectively. Since in our exemplary laminate the ratio between $\phaseA{c_{\pm}}$ and $\phaseB{c_{\pm}}$ is kept constant, it is sufficient to characterize the dependency of $\phaseA{c_{\pm}}$ on the pre-strain. This dependency is characterized in Fig.$\ $\[fig:quasi-pressure-solitary\], where contour plots of $\phaseA{\pressurewavespeed}$ (panel a) and $\phaseA{\shearwavespeed}$ (panel b) are shown as functions of $\left\{ \pressurestrain,\,\shearstrain\right\} $ pairs. We demonstrate how these maps are used to determine if a solitary wave will propagate using a representative example, by considering two different sets of uniform strains and two different sets of localized strains, namely,

$$\begin{aligned} & U_{C}=\left(\initial{\pressurestrain}=-0.35,\ \initial{\shearstrain}=1.98\right),\ \ \ U_{T}=\left(\initial{\pressurestrain}=0.2,\ \initial{\shearstrain}=2.15\right),\\
 & L_{1}=\left(\amp{\pressurestrain}=0.3,\ \amp{\shearstrain}=2.3\right),\ \ \ L_{2}=\left(\amp{\pressurestrain}=-0.6,\ \amp{\shearstrain}=1.92\right).
\end{aligned}$$

The dashed red curve is the solution of Eq.$\ $ with either $U_{C}$ or $U_{T}$ as the initial condition. As mentioned previously, the $\shearstrain-\pressurestrain$ relation of quasi-pressure solitary waves generated with inhomogeneous pre-strain whose uniform part is the same as $U_{C}$ (or $U_{T}$) approximately follows this curve. Furthermore, the relative location of the localized part of the strain with respect to the uniform strain determines the direction of the strain evolution along the red curve, such that $L_{1}$ and $L_{2}$ correspond to an increase and a decrease in $\pressurestrain$, respectively. When the initial strain in the laminate is the sum of the uniform strain $U_{T}$ and the localized strain $L_{1}$, we have that $\phaseA{\pressurewavespeed}$ is monotonically increasing along the dashed red curve, whose general direction is towards $L_{1}$. Therefore, this initial state supports a quasi-pressure solitary wave. When the initial strain is the sum of $U_{T}$ and $L_{2}$, we have that $\phaseA{\pressurewavespeed}$ is monotonically decreasing along the dashed red curve, whose general direction is towards $L_{2}$. Thus, this initial state does not support solitary waves and the generated wave will disperse. Similarly, when the uniform part of the initial strain is given by $U_{C}$ and the localized part is $L_{1}$, solitary waves cannot propagate, whereas when the localized part is $L_{2}$ such waves do propagate. We emphasize that without the observation that the $\shearstrain-\pressurestrain$ relations of the laminate approximately follow their counterparts in the homogeneous case, it would have not been possible to estimate *a priori* what should be the localized strains that generate a certain solitary wave.

These contour maps are also useful in predicting how the axial strain changes during the propagation of the solitary waves, *i.e.,* will it decrease increase. If the former occurs, we term the wave a compressive wave, owing to its tendency to decrease the axial strain in the laminate, notwithstanding the fact that the sign of axial strain in the laminate may be positive. Similarly, if the latter occurs, we term the wave a *tensile* solitary wave, notwithstanding the fact that the sign of axial strain may be negative. For compressive solitary waves to form, it is required that along the loading path of the solitary wave, the value of $\pressurestrain$ is decreasing and the phase velocity is monotonically increasing. Specifically for quasi-pressure waves, the domain for which this is possible is identified with an initial uniform strain that belongs to the domain left to the black curve in Fig.$\ $\[fig:quasi-pressure-solitary\](a). Conversely, the formation of tensile solitary waves requires that along the loading path of the solitary wave, the value of $\pressurestrain$ is increasing and the phase velocity is monotonically increasing. The domain of initial uniform strains that can satisfy this criterion is to the right of the curve. The condition for tensile waves cannot be satisfied for quasi-shear solitary waves where only compressive solitary waves are supported, independently of the initial strain.

Interestingly, the prescribed strains supporting the formation of solitary waves in laminates and shock waves in homogeneous media are similar, as they both require that the wave velocity will increase monotonically along the loading path. The latter problem of shocks in compressible Gent materials was studied by @ZIV2019mom when considering coupled shear-pressure shocks, and later by @CHOCKALINGAM2020103746 who restricted attention to shear shocks in the specialized incompressible case. The configuration considered by @ZIV2019mom consists of a pre-strained semi-infinite Gent material with the same model parameters as our laminate, namely, Eq.$\ $, subjected to an impact at the boundary. The analogy between the problems is manifested by the similarity of the prescribed strains supporting solitary waves and shock waves. Specifically, if the magnitude of the localized part of the initial strain is substituted as the magnitude of the impact at the boundary of the homogeneous half-space, then Fig.$\ $\[fig:quasi-pressure-solitary\] also predicts which impacts yield tensile shock waves and which impacts yield compressive shock waves. In other words, impacts that give rise to monotonically increasing phase velocities and in turn shocks in a homogeneous Gent material, also give rise to monotonically increasing phase velocities and in turn may form solitary waves in laminates made of the same Gent material whose mass density is modulated. This connection sheds light on the balance between the nonlinearity in the system and its dispersion: the constitutive nonlinearity causes the pulse to steepen, while the scattering at the interfaces disperses it and generates solitary waves [@hussein2018nonlinear], otherwise the steepening would have result in shock formation [@yong2003solitary].

Summary\[sec:Summary\]
======================

@ziv2019b have developed a finite-volume method to solve the equations governing finite-amplitude smooth waves with two coupled components in nonlinear compressible laminates. The application of this method to simulate the response of pre-strained compressible Gent two-phase laminates has revealed the generation of vector solitary waves, whose axial and transverse polarizations are coupled. Here, we have used the method to conduct a large set of numerical experiments with different initial conditions, in order to characterize these waves.

We have classified the possible vector solitary waves according to two types of profiles. The first classwhich was observed by @ziv2019bhas a profile that varies between the two phases, but remains fixed as the wave passes through the different periodic cells. We have termed this type as quasi-pressure, since it reduces to the standard pressure waves in the limit of small strains and vanishing heterogeneity, and showed that its components follow the profile of the sech$^{2}$ function.

The second class we have identified was not generated by @ziv2019b, and is more interesting; we termed it quasi-shear with similar reasoning. The profile and velocity of this class varies not only between different phases, but also between different cells in a periodic manner without dispersing in time. This long-lived oscillation is a manifestation of a permanent and periodic energy transfer between the axial and transverse polarizations, which we consider as internal mode of the vector solitary wave. To the best of our knowledge, this is the first report of such oscillating solitary waves in continuum elastodynamics.

Using initial conditions that generate a single solitary wave in each numerical experiment, we have studied the velocity-amplitude relation, and found that both types propagate faster at higher amplitudes. This observation agrees with the result of @ziv2019b, who generated a train of solitary waves of the first type, and found that the taller waves in the train are faster than the shorter waves. Notably, the vector solitary waves observed in the discrete mechanical systems of @deng2017 [@deng2019prl] exhibit an opposite behavior, namely, they are faster at smaller amplitudes. Since the systems considered by @deng2017 [@deng2019prl] and the system considered here are fundamentally different, this difference in the amplitude-dependent nature of the vector solitary waves is reasonable. In the limit of small strains, we have calculated velocity of low-frequency, long-wavelength Bloch-Floquet waves in the laminate, and found it serves as a lower bound the the velocity of the vector solitary waves.

In the last part of the paper, we have described a procedure for selecting initial strains that generate selected vector solitary waves. We emphasize that without the observation that the $\shearstrain-\pressurestrain$ relations of the laminate approximately follow their counterparts in the homogeneous case, it would have not been possible to estimate *a priori* what should be the localized strains that generate a certain solitary wave. We find this agreement remarkable, since it correlates between two seemingly unrelated problemsimpact of homogeneous media and localized initial strain in laminated media, and note it relies on the fact that Eq.$\ $ is independent of the mass density. The procedure for generating selected solitary waves also provides insights on the way that the axial strain in the laminate changes as the resultant waves propagate. This allows for an additional classification of the waves to tensile and compressive solitary waves, according to their tendency to increase or decrease that strain.

Future work can follow different paths from this point. One possibility is to attempt developing a homogenized system of equations [@yong2003solitary], based on which it may be possible to carry out further analytical investigation, *e.g.*, as done by @andrianov2013dynamic in the one-dimensional problem. This may be turn useful, for example, in identifying the solution function that describes the quasi-shear waves. Another possibility is to explore if additional types of vector solitary waves will form when using different parameters for the compressible Gent model than the set used here, or alternatively employ other types of constitutive models, such as the @ogden1972 and models. It would be also interesting to study the collision of vector solitary waves in the laminate, which may result with anomalous phenomena, as observed by @deng2019prl in their case of the discrete two-dimensional mechanical model. As pointed out earlier, the vector solitary waves reported here are different from those observed by @deng2019prl, and hence expect that their response to collisions will also be different.

Acknowledgments {#acknowledgments .unnumbered}
===============

We thank an anonymous reviewer whose constructive comments helped improving this paper. We acknowledge the support of the Israel Science Foundation, funded by the Israel Academy of Sciences and Humanities (Grant no.$\ $1912/15), the United States-Israel Binational Science Foundation (Grant no. 2014358), and Ministry of Science and Technology (grant no. 880011).

 {#section .unnumbered}

[71]{} \[1\][\#1]{} \[1\][`#1`]{} urlstyle \[1\][doi: \#1]{}

Igor V Andrianov, Vladyslav V Danishevskyy, Oleksandr I Ryzhkov, and Dieter Weichert. Dynamic homogenization and wave propagation in a nonlinear 1d composite material. *Wave Motion*, 500 (2):0 271–281, 2013.

Igor V Andrianov, Vladyslav V Danishevskyy, Oleksandr I Ryzhkov, and Dieter Weichert. Numerical study of formation of solitary strain waves in a nonlinear elastic layered composite material. *Wave Motion*, 510 (3):0 405–417, 2014.

E M Arruda and M C Boyce. . *J. Mech. Phys. Solids*, 41:0 389–412, 1993.

DS Bale, RJ Leveque, S Mitran, and JA Rossmanith. . *[SIAM JOURNAL ON SCIENTIFIC COMPUTING]{}*, [24]{}0 ([3]{}):0 [955–978]{}, [2002]{}. ISSN [1064-8275]{}. [doi: ]{}[[10.1137/S106482750139738X]{}]{}.

Amit Bandyopadhyay, Sahar Vahabzadeh, Anish Shivaram, and Susmita Bose. Three-dimensional printing of biomaterials and soft materials. *MRS Bulletin*, 400 (12):0 1162–1169, [DEC]{} 2015. ISSN [0883-7694]{}. [doi: ]{}[[10.1557/mrs.2015.274]{}]{}.

Fabio Baronio, Antonio Degasperis, Matteo Conforti, and Stefan Wabnitz. Solutions of the vector nonlinear schrödinger equations: Evidence for deterministic rogue waves. *Phys. Rev. Lett.*, 109:0 044102, Jul 2012. [doi: ]{}[10.1103/PhysRevLett.109.044102]{}. URL <https://link.aps.org/doi/10.1103/PhysRevLett.109.044102>.

Katia Bertoldi, Vincenzo Vitelli, Johan Christensen, and Martin van Hecke. Flexible mechanical metamaterials. *Nature Reviews Materials*, 20 (11):0 17066, 2017.

DR Bland. Plane isentropic large displacement simple waves in a compressible elastic solid. *Zeitschrift f[ü]{}r angewandte Mathematik und Physik ZAMP*, 160 (6):0 752–769, 1965.

J.P. Boyd. Dynamical meteorology | solitary waves. In Gerald R. North, John Pyle, and Fuqing Zhang, editors, *Encyclopedia of Atmospheric Sciences (Second Edition)*, pages 417 – 422. Academic Press, Oxford, second edition edition, 2015. ISBN 978-0-12-382225-3. [doi: ]{}[https://doi.org/10.1016/B978-0-12-382225-3.00374-1]{}. URL <http://www.sciencedirect.com/science/article/pii/B9780123822253003741>.

David K. Campbell, Jonathan F. Schonfeld, and Charles A. Wingate. Resonance structure in kink-antikink interactions in $\phi$4 theory. *Physica D: Nonlinear Phenomena*, 90 (1):0 1–32, 1983. [doi: ]{}[https://doi.org/10.1016/0167-2789(83)90289-0]{}. URL <http://www.sciencedirect.com/science/article/pii/0167278983902890>.

S. Catheline, J.-L. Gennisson, M. Tanter, and M. Fink. Observation of shock transverse waves in elastic media. *Phys. Rev. Lett.*, 91:0 164301, Oct 2003. [doi: ]{}[10.1103/PhysRevLett.91.164301]{}. URL <https://link.aps.org/doi/10.1103/PhysRevLett.91.164301>.

S. Chockalingam and T. Cohen. Shear shock evolution in incompressible soft solids. *Journal of the Mechanics and Physics of Solids*, 134:0 103746, 2020. ISSN 0022-5096. [doi: ]{}[https://doi.org/10.1016/j.jmps.2019.103746]{}. URL <http://www.sciencedirect.com/science/article/pii/S0022509619306684>.

Johan Christensen, Muamer Kadic, Oliver Kraft, and Martin Wegener. Vibrant times for mechanical metamaterials. *MRS Communications*, 5:0 453–462, 2015. [doi: ]{}[10.1557/mrc.2015.51]{}.

Adrian R. Cioroianu and Cornelis Storm. Normal stresses in elastic networks. *Phys. Rev. E*, 88:0 052601, Nov 2013. [doi: ]{}[10.1103/PhysRevE.88.052601]{}. URL <https://link.aps.org/doi/10.1103/PhysRevE.88.052601>.

Richard V Craster and S[é]{}bastien Guenneau. *[Acoustic metamaterials: Negative refraction, imaging, lensing and cloaking]{}*, volume 166. Springer Science [&]{} Business Media, 2012.

Thierry Dauxois and Michel Peyrard. *Physics of solitons*. Cambridge University Press, 2006.

L Davison. . *Journal of the Mechanics and Physics of Solids*, 140 (5):0 249–270, 1966. ISSN 0022-5096. [doi: ]{}[https://doi.org/10.1016/0022-5096(66)90022-6]{}. URL <http://www.sciencedirect.com/science/article/pii/0022509666900226>.

Lee Davison. *Fundamentals of Shock Wave Propagation in Solids*. Springer-Verlag Berlin Heidelberg, 2008.

B Deng, J R Raney, V Tournat, and K Bertoldi. . *Phys. Rev. Lett*, 118:0 204102, 2017.

Bolei Deng, Pai Wang, Qi He, Vincent Tournat, and Katia Bertoldi. Metamaterials with amplitude gaps for elastic solitons. *Nature communications*, 90 (1):0 3410, 2018.

Bolei Deng, Chengyang Mo, Vincent Tournat, Katia Bertoldi, and Jordan R Raney. Focusing and mode separation of elastic vector solitons in a 2d soft mechanical metamaterial. *Physical review letters*, 1230 (2):0 024101, 2019.

Bolei Deng, Vincent Tournat, Pai Wang, and Katia Bertoldi. Anomalous collisions of elastic vector solitons in mechanical metamaterials. *Physical Review Letters*, 122:0 044101, 2019.

Bolei Deng, Liyuan Chen, Donglai Wei, Vincent Tournat, and Katia Bertoldi. Pulse-driven robot: Motion via solitary waves. *Science Advances*, 60 (18), 2020. [doi: ]{}[10.1126/sciadv.aaz1166]{}. URL <https://advances.sciencemag.org/content/6/18/eaaz1166>.

Bolei Deng, Pai Wang, Vincent Tournat, and Katia Bertoldi. Nonlinear transition waves in free-standing bistable chains. *Journal of the Mechanics and Physics of Solids*, 136:0 103661, 2020. ISSN 0022-5096. [doi: ]{}[https://doi.org/10.1016/j.jmps.2019.07.004]{}. URL <http://www.sciencedirect.com/science/article/pii/S0022509619303321>. The Davide Bigoni 60th Anniversary Issue.

David Esp[í]{}ndola, Stephen Lee, and Gianmarco Pinton. Shear shock waves observed in the brain. *Physical Review Applied*, 80 (4):0 044024, 2017.

R Ganesh and S Gonella. . *J. Mech. Phys. Solids*, 99:0 272–288, 2017. ISSN 0022-5096. [doi: ]{}[http://doi.org/10.1016/j.jmps.2016.11.001]{}. URL <http://www.sciencedirect.com/science/article/pii/S0022509616305440>.

David Garcia, Ziling Wu, Jee Yun Kim, Hang Z. Yu, and Yunhui Zhu. Heterogeneous materials design in additive manufacturing: Model calibration and uncertainty-guided model selection. *Additive Manufacturing*, 27:0 61 – 71, 2019. ISSN 2214-8604. [doi: ]{}[https://doi.org/10.1016/j.addma.2019.02.014]{}. URL <http://www.sciencedirect.com/science/article/pii/S2214860418307711>.

A N Gent. . *Rubber Chem. Technol.*, 69:0 59–61, 1996.

R Getz, D M Kochmann, and G Shmuel. . *Int. J. Solids Struct.*, 113–114:0 24–36, 2017. ISSN 0020-7683. [doi: ]{}[http://doi.org/10.1016/j.ijsolstr.2016.10.002]{}. URL <http://www.sciencedirect.com/science/article/pii/S0020768316302931>.

C. O. Horgan and J. G. Murphy. Poynting and reverse poynting effects in soft materials. *Soft Matter*, 13:0 4916–4923, 2017. [doi: ]{}[10.1039/C7SM00992E]{}. URL <http://dx.doi.org/10.1039/C7SM00992E>.

MI Hussein and R Khajehtourian. Nonlinear bloch waves and balance between hardening and softening dispersion. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 4740 (2217):0 20180173, 2018.

Muamer Kadic, Graeme W. Milton, Martin van Hecke, and Martin Wegener. 3d metamaterials. *Nature Reviews Physics*, 10 (3):0 198–210, 2019. [doi: ]{}[10.1038/s42254-018-0018-y]{}. URL <https://doi.org/10.1038/s42254-018-0018-y>.

Kenichi Kasamatsu and Makoto Tsubota. Modulation instability and solitary-wave formation in two-component bose-einstein condensates. *Phys. Rev. A*, 74:0 013617, Jul 2006. [doi: ]{}[10.1103/PhysRevA.74.013617]{}. URL <https://link.aps.org/doi/10.1103/PhysRevA.74.013617>.

Shmuel Katz and Sefi Givli. Solitary waves in a nonintegrable chain with double-well potentials. *Phys. Rev. E*, 100:0 032209, Sep 2019. [doi: ]{}[10.1103/PhysRevE.100.032209]{}. URL <https://link.aps.org/doi/10.1103/PhysRevE.100.032209>.

Romik Khajehtourian and Mahmoud I. Hussein. Nonlinear dispersion relation predicts harmonic generation in wave motion, 2019.

Yuri S Kivshar and Govind Agrawal. *Optical solitons: from fibers to photonic crystals*. Academic press, 2003.

Yuri S Kivshar, Dmitry E Pelinovsky, Thierry Cretegny, and Michel Peyrard. Internal modes of solitary waves. *Physical review letters*, 800 (23):0 5032, 1998.

Diederik Johannes Korteweg and Gustav de Vries. Xli. on the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science*, 390 (240):0 422–443, 1895.

Stefano Lepri and Giulio Casati. Asymmetric wave propagation in nonlinear systems. *Phys. Rev. Lett.*, 106:0 164101, Apr 2011. [doi: ]{}[10.1103/PhysRevLett.106.164101]{}. URL <https://link.aps.org/doi/10.1103/PhysRevLett.106.164101>.

Randall J. LeVeque. Wave propagation algorithms for multidimensional hyperbolic systems. *Journal of Computational Physics*, 1310 (2):0 327 – 353, 1997. ISSN 0021-9991. [doi: ]{}[https://doi.org/10.1006/jcph.1996.5603]{}. URL <http://www.sciencedirect.com/science/article/pii/S002199919695603X>.

Randall J. LeVeque. *Finite Volume Methods for Hyperbolic Problems*. Cambridge University Press, Cambridge, 2002. ISBN 9780521009249. [doi: ]{}[DOI: 10.1017/CBO9780511791253]{}. URL <https://www.cambridge.org/core/books/finite-volume-methods-for-hyperbolic-problems/97D5D1ACB1926DA1D4D52EAD6909E2B9>.

Randall J LeVeque. Finite-volume methods for non-linear elasticity in heterogeneous media. *International journal for numerical methods in fluids*, 400 (1-2):0 93–104, 2002.

Randall J LeVeque and Darryl H Yong. Solitary waves in layered nonlinear media. *SIAM Journal on Applied Mathematics*, 630 (5):0 1539–1560, 2003.

O. Lopez-Pamies and P. Ponte Casta[ñ]{}eda. Homogenization-based constitutive models for porous elastomers and implications for macroscopic instabilities: Ii—results. *Journal of the Mechanics and Physics of Solids*, 550 (8):0 1702 – 1728, 2007. ISSN 0022-5096. [doi: ]{}[https://doi.org/10.1016/j.jmps.2007.01.008]{}. URL <http://www.sciencedirect.com/science/article/pii/S0022509607000269>.

Ben Lustig and Gal Shmuel. . *Journal of the Mechanics and Physics of Solids*, 117:0 37–53, 2018. ISSN 0022-5096. [doi: ]{}[https://doi.org/10.1016/j.jmps.2018.04.008]{}. URL <http://www.sciencedirect.com/science/article/pii/S0022509618302321>.

B. I. Mantsyzov. Gap 2$\pi$ pulse with an inhomogeneously broadened line and an oscillating solitary wave. *Phys. Rev. A*, 51:0 4939–4943, Jun 1995. [doi: ]{}[10.1103/PhysRevA.51.4939]{}. URL <https://link.aps.org/doi/10.1103/PhysRevA.51.4939>.

G. Marckmann and E. Verron. Comparison of hyperelastic models for rubber-like materials. *Rubber Chemistry and Technology*, 790 (5):0 835–858, 2020/02/16 2006. [doi: ]{}[10.5254/1.3547969]{}. URL <https://doi.org/10.5254/1.3547969>.

L. Angela Mihai and Alain Goriely. How to characterize a nonlinear elastic material? a review on nonlinear constitutive parameters in isotropic finite elasticity. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 4730 (2207):0 20170607, 2020/06/03 2017. [doi: ]{}[10.1098/rspa.2017.0607]{}. URL <https://doi.org/10.1098/rspa.2017.0607>.

Chengyang Mo, Jaspreet Singh, Jordan R Raney, and Prashant K Purohit. Cnoidal wave propagation in an elastic metamaterial. *Physical Review E*, 1000 (1):0 013001, 2019.

Neel Nadkarni, Chiara Daraio, and Dennis M Kochmann. Dynamics of periodic mechanical structures containing bistable elastic elements: From elastic to solitary wave propagation. *Physical Review E*, 900 (2):0 023204, 2014.

Amir Nasrollahi, Wen Deng, Piervincenzo Rizzo, Alex Vuotto, and Julie M. Vandenbossche. Nondestructive testing of concrete using highly nonlinear solitary waves. *Nondestructive Testing and Evaluation*, 320 (4):0 381–399, 2017. [doi: ]{}[10.1080/10589759.2016.1254212]{}. URL <https://doi.org/10.1080/10589759.2016.1254212>.

Pedro Franco Navarro, David J Benson, and Vitali F Nesterenko. Nature of short, high-amplitude compressive stress pulses in a periodic dissipative laminate. *Physical Review E*, 920 (6):0 062917, 2015.

J Niemczura and K Ravi-Chandar. . *Journal of the Mechanics and Physics of Solids*, 590 (2):0 442–456, 2011. ISSN 0022-5096. [doi: ]{}[https://doi.org/10.1016/j.jmps.2010.09.007]{}. URL <http://www.sciencedirect.com/science/article/pii/S0022509610001833>.

R W Ogden. *[Non-Linear Elastic Deformations]{}*. Dover Publications, New York, 1997.

Raymond William Ogden. Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids. *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences*, 3280 (1575):0 567–583, 1972. [doi: ]{}[10.1098/rspa.1972.0096]{}. URL <https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1972.0096>.

John Henry Poynting. On pressure perpendicular to the shear planes in finite pure shears , and on the lengthening of loaded wires when twisted. *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character*, 820 (557):0 546–559, 2020/06/03 1909. [doi: ]{}[10.1098/rspa.1909.0059]{}. URL <https://doi.org/10.1098/rspa.1909.0059>.

J Raney, N Nadkarni, C Daraio, D M Kochmann, J A Lewis, and K Bertoldi. . *Proc. Natl. Acad. Sci. U. S. A.*, 2016.

J R Raney and J A Lewis. . *MRS Bulletin*, 400 (11):0 943–950, 2015.

Solomon M. Saltiel, Dragomir N. Neshev, Robert Fischer, Wieslaw Krolikowski, Ady Arie, and Yuri S. Kivshar. Generation of second-harmonic conical waves via nonlinear bragg diffraction. *Phys. Rev. Lett.*, 100:0 103902, Mar 2008. [doi: ]{}[10.1103/PhysRevLett.100.103902]{}. URL <https://link.aps.org/doi/10.1103/PhysRevLett.100.103902>.

Fadil Santosa and William W Symes. A dispersive effective medium for wave propagation in periodic composites. *SIAM Journal on Applied Mathematics*, 510 (4):0 984–1005, 1991.

S.A. Silling. Solitary waves in a peridynamic elastic solid. *Journal of the Mechanics and Physics of Solids*, 96:0 121 – 132, 2016. ISSN 0022-5096. [doi: ]{}[https://doi.org/10.1016/j.jmps.2016.06.001]{}. URL <http://www.sciencedirect.com/science/article/pii/S0022509616302903>.

S. E. Skipetrov and R. Maynard. Instabilities of waves in nonlinear disordered media. *Phys. Rev. Lett.*, 85:0 736–739, Jul 2000. [doi: ]{}[10.1103/PhysRevLett.85.736]{}. URL <https://link.aps.org/doi/10.1103/PhysRevLett.85.736>.

Ankit Srivastava. . *International Journal of Smart and Nano Materials*, 60 (1):0 41–60, 2015.

George I. Stegeman and Mordechai Segev. Optical spatial solitons and their interactions: Universality and diversity. *Science*, 2860 (5444):0 1518–1523, 1999. ISSN 0036-8075. [doi: ]{}[10.1126/science.286.5444.1518]{}. URL <https://science.sciencemag.org/content/286/5444/1518>.

Piotr Szankowski, Marek Trippenbach, Eryk Infeld, and George Rowlands. Oscillating solitons in a three-component bose-einstein condensate. *Phys. Rev. Lett.*, 105:0 125302, Sep 2010. [doi: ]{}[10.1103/PhysRevLett.105.125302]{}. URL <https://link.aps.org/doi/10.1103/PhysRevLett.105.125302>.

R. L. Truby and J. A. Lewis. Printing soft matter in three dimensions. *Nature*, 540:0 371–378, 2016.

Hiromi Yasuda, Yasuhiro Miyazawa, Efstathios G Charalampidis, Christopher Chong, Panayotis G Kevrekidis, and Jinkyu Yang. Origami-based impact mitigation via rarefaction solitary wave creation. *Science advances*, 50 (5):0 eaau2835, 2019.

Jiangyi Zhang, Vicente Romero-Garc[í]{}a, Georgios Theocharis, Olivier Richoux, Vassos Achilleos, and Dimitrios J Frantzeskakis. Bright and gap solitons in membrane-type acoustic metamaterials. *Physical Review E*, 960 (2):0 022214, 2017.

Zhen Zhang, Leonid I. Manevitch, Valeri Smirnov, Lawrence A. Bergman, and Alexander F. Vakakis. Extreme nonlinear energy exchanges in a geometrically nonlinear lattice oscillating in the plane. *Journal of the Mechanics and Physics of Solids*, 110:0 1 – 20, 2018. ISSN 0022-5096. [doi: ]{}[https://doi.org/10.1016/j.jmps.2017.09.007]{}. URL <http://www.sciencedirect.com/science/article/pii/S0022509616308845>.

Ron Ziv and Gal Shmuel. Smooth waves and shocks of finite amplitude in soft materials. *Mechanics of Materials*, 135:0 67 – 76, 2019. ISSN 0167-6636. [doi: ]{}[https://doi.org/10.1016/j.mechmat.2019.05.002]{}. URL <http://www.sciencedirect.com/science/article/pii/S016766361930002X>.

Ron Ziv and Gal Shmuel. Observation of vector solitary waves in soft laminates using a finite-volume method. *International Journal of Non-Linear Mechanics*, 124:0 103502, 2020. ISSN 0020-7462. [doi: ]{}[https://doi.org/10.1016/j.ijnonlinmec.2020.103502]{}. URL <http://www.sciencedirect.com/science/article/pii/S0020746220301645>.

[^1]: Corresponding author. Tel.: +1 972 778871613. *E-mail address*: meshmuel@technion.ac.il (G. Shmuel).

[^2]: Solitary waves whose velocity exceeds the wave velocity in the linear quasi-static limit are termed *supersonic*.

[^3]: Roughly speaking, smooth waves propagate when the velocity at the tail of the wave does not exceed the velocity at its front, and changes monotonically in between these ends [@DAVISON1966249; @ZIV2019mom; @CHOCKALINGAM2020103746].

[^4]: For an excellent treatise on finite-volume methods, the reader is referred to the book of @LeVeque2002book.

[^5]: While the shape and velocity of the wave are not constant but rather oscillate about a mean, we retain the term *solitary*, in agreement with the terminology used in other fields [@Mantsyzov1995pra; @kivshar1998internal; @Szankowski2010prl]. This is due to the fact that the wave is stable, *i.e.*, it does not spread (or get squeezed) in the course of its propagation, in spite of the fact that it is governed by nonlinear equations. This also agrees with the interpretation of these oscillations as a result of internal mode of a solitary wave.
