Breathers in lattices with alternating strain-hardening and strain-softening interactions

Nonlinear lattices play a pivotal role in the modeling of a wide range of physical, biological, chemical, and electrical systems FPUreview ; pgk:2011 ; moti ; Morsch ; sievers . One of the prototypical examples, and most relevant for the present study, is the Fermi-Pasta-Ulam-Tsingou (FPUT) chain, which is simply a mass-spring system with nonlinear springs FPU55 . While the types of behavior that are possible in FPUT chains span a wide range, the structure we will focus on here is the so-called discrete breather, which is a solution that is localized in space and periodic in time. Discrete breathers (or just breathers) have been studied in a variety of physical systems, including, optical waveguide arrays and photorefractive crystals moti , micromechanical oscillator arrays sievers2 , Josephson-junction ladders alex ; alex2 , electrical lattices with nonlinear elements remoissenet , layered antiferromagnetic crystals lars3 ; lars4 , halide-bridged transition metal complexes swanson , dynamical models of the DNA double strand Peybi , Bose–Einstein condensates in optical lattices Morsch , and magnetic lattices moleron ; Mehrem2017 ; Marc2017 ; magBreathers ; Chong_2021 ; cantilevers among others.

One mechanical realization of an FPUT lattice that has attracted significant attention is the so-called granular chain, in which case the nonlinear interaction is described by a power-law with nonlinearity exponent $p=3/2>1$. Granular chains consist of closely packed arrays of particles that interact elastically upon compression Nester2001 ; granularBook ; yuli_book ; Lindenberg_review ; gc_review ; sen08 . The contact force can be tuned to yield responses that range from near-linear to purely nonlinear, and the effective stiffness properties can also be easily changed by modifying the material, geometry, or contact angle of the elements in contact Nester2001 , although the stiffness will be of the strain-hardening type (implying an increase of elastic modulus with strain). Due to the strain-hardening nature of the granular chain, the monomer one (where all particles are identical) does not have genuine breather solutions JAMES201339 (although so-called dark breathers are possible dark ; dark2 ). Breathers are possible in granular chains with defects Theocharis2009 ; Nature11 or in mass-dimer chains Boechler2010 ; Theocharis10 ; hooge12 and have been studied in great detail.

More recently, strain-softening materials (implying a decrease of elastic modulus with strain) have been emerging as a new playground for the formation of nonlinear waves Deng2017 ; Raney , and the dynamics here will be distinct from their strain-hardening counterparts. For instance, rarefaction waves form in lattices with strain-softening interactions shock_trans_granular ; yasuda ; HEC_DSW instead of the classic solitary waves in the case of hardening lattices Nester2001 . Examples of lattices that have been argued to exhibit strain-softening interaction include tensegrity carpe , origami origami ; origami2 , stainless steel cylinders separated by polymer foams RarefactionNester , and hollow elliptical cylinders HEC_DSW . These strain-softening lattices can be modeled with FPUT type lattices with a power-law interaction with nonlinear exponent $0<p<1$. Breathers in such strain-softening lattices are far less studied than their strain-hardening counterparts.

More novel still, to the best of our knowledge, is the study of breathers in lattices that feature both strain-hardening and strain-softening interactions, the focus of the present study. Such a lattice could be possible in an origami setting where each origami unit cell is tailored to have strain-hardening or strain-softening behavior. This is possible, in principle, given the highly tunable nature of origami by virtue of changes in crease-angle (see Bistable_origami for an example of an origami unit cell that can be tuned to exhibit both softening and hardening behavior, depending, e.g., on the folding angle for Tachi-Miura polyhedra origami ). Another possible realization is a chain of stainless steel cylinders in which only a subset of contacts is separated with a polymer foam in order to induce the strain-softening behavior.

The paper is structured as follows. The model equations and basic linear analysis of such a periodically alternating nonlinearity are given in Sec. II. Our main aim herein is to explore this model and especially its nonlinear time-periodic solutions. The numerical study of such breather (exponentially localized in space and periodic in time) waveforms of the damped-driven model (being more realistic for the potential experimental realization of such structures) is detailed in Sec. III. The ideal case of the Hamiltonian lattice is considered in Sec. IV, where the Nonlinear Schrödinger (NLS) equation is derived using a multiple-scale analysis in order to construct an analytical prediction of both acoustic and optical breathers. Section V concludes the paper and offers suggestions for future studies that build on the present one.

We now turn to the study of time-periodic solutions to the nonlinear problem with damping and driving. We will use primarily numerical computations for their study, but we discuss an analytical approximation in the next section using a reduction to the NLS equation in the Hamiltonian limit of the problem.

We compute time-periodic orbits, and their spectral stability, of Eq. (1) with period $T_{d}=1/f_{d}$ with high precision using a standard Newton-type procedure, see Appendix A for details. To investigate the dynamical stability of the obtained states, a Floquet analysis is used to compute the multipliers associated with the solutions. If a solution has all Floquet multipliers on or inside the unit circle, the solution is called (spectrally) stable. An instability that results from a multiplier on the (positive) real line is called a real (exponential in nature) instability. However, there can also be oscillatory instabilities, which correspond to complex-conjugate pairs of Floquet multipliers lying outside the unit circle (in the complex plane). In the bifurcation diagrams that are to follow in this paper, solid blue segments will correspond to stable parametric regions while the dashed-dotted red segments correspond to real unstable regions and green dashed-dotted segments correspond to oscillatorily unstable regions. The bifurcation diagrams are obtained using pseudo-arclength continuation as it is implemented in the software package AUTO Doedel . In all of the bifurcation diagrams that follow below, the monitored quantity that will be used is the average energy over one period $T_{d}$ defined as:

where

is the total energy, and

is the (discrete) energy density. Note that $V_{n}$ in Eq. (18) is given by Eq. (6).

For the numerical computations, we consider a chain with $N=42$ particles, and set $p=1$, $\delta=0.3$, and $\gamma=8\times 10^{-3}$ (unless explicitly stated otherwise). The drive amplitude $a_{d}$ and frequency $f_{d}$ will be used as continuation parameters. Let us start our presentation by considering Figs. 2-5, which correspond to frequency continuations. In particular, we consider the following cases: linear monomer chain ($p_{n}\equiv p$) in Fig. 2, a strain-hardening monomer chain ($p_{n}=p+\delta$) in Fig. 3, a strain-softening monomer chain ($p_{n}=p-\delta$) in Fig. 4, and a dimer chain with alternating stiffness ($p_{n}$ given by Eq. (2)) in Fig. 5. Figures 2-5 showcase the average energy [cf. Eq. (16)] as a function of the driving frequency $f_{d}$. When $p_{n}\equiv 1\,\forall n=1,\dots,N$, i.e., for the case corresponding to a linear lattice, it can be discerned from Fig. 2 that the resonant peaks match perfectly with the resonant frequencies of the respective linear problem shown with vertical dashed-dotted black lines (and solved numerically) in this case, as expected. Adding nonlinearity to the system will cause these linear resonant peaks to deform, or “bend”. For example, supposing that all contacts are strain-hardening, the resonant peaks will bend toward smaller frequencies, as shown in Fig. 3 where $p_{n}=p+\delta>1$. If all contacts are strain-softening, the resonant peaks will bend toward larger frequencies as shown in Fig. 4 where $p_{n}=p-\delta<1$. The observations in Figs. 3 and 4 are expected, and are inline with the established behavior for oscillators with hardening or softening nonlinearities Mook . In the bottom panels of Figs. 3 and 4, the resonant frequencies are depicted by vertical dashed-dotted lines (similar to Fig. 2) in order to ease the visualization of the bending of the resonant curves (to the left or right).

Having investigated the baseline cases of pure hardening and pure softening stiffness, we now turn to the aspect of alternating stiffness, i.e., Eq. (1) with nonlinearity exponents $p_{n}$ given by Eq. (2). The frequency continuation for this case is shown in Fig. 5 for $a_{d}=0.15$. Notice the resonant peaks close to the top edge of the acoustic band and the bottom edge of the optical band all bend toward the gap due to the interplay of softening and hardening. This finding is shown in detail in the bottom panel of Fig. 5, which is a zoom of the top panel. Note that the eigenfrequency modes (shown with vertical dashed-dotted lines) help highlight this effect as well as the existence of a frequency gap. Note that in traditional mass-dimer lattices, resonant peaks from either the top of the acoustic band or the bottom of the optical band will bend into the frequency gap, but not both, like in the present case granularBook . The nature of the bending in the lattice with alternating stiffness can be understood in the small amplitude limit by inspection of the strain of the eigenvectors of the linear problem. In particular, within the acoustic band, the larger amplitude strains occur for odd indices, see the bottom left panel of Fig. 1. According to Eq. (2), the odd numbered interactions are of the softening type. Hence, one would expect softening-type dynamics, involving resonant peaks bending toward the right, as seen in Fig. 5. The situation in the optical band is reversed. The even numbered interactions have higher amplitude in the strain (see the bottom right panel of Fig. 1) indicating that hardening-type behavior should be observed, involving resonant peaks bending to the left, like in Fig. 5.

A discussion about the stability characteristics of the pertinent branches of Figs. 2-5 is now in order. In the linear lattice of Fig. 2, we observe that the entire branch is spectrally stable. In the nonlinear cases of Figs. 3-5, the branches are mostly stable, with exceptions occurring typically between two turning points. These segments, in line with what is known for classical oscillator systems Mook , are found to be unstable. For example, in Fig. 3 and for $f_{d}\approx[0.29,0.35]$, the alternation of stable and unstable segments between turning points can be discerned. The instability in this case corresponds to an exponential one (and is shown with dashed-dotted red segments). A similar observation is found in Fig. 4 for drive frequencies $f_{d}\approx[0.12,0.25]$. Interestingly, in some portions of the relevant branches, and in particular those of Figs. 3 and 5 (see the bottom panels therein), we find oscillatory instabilities. In these cases, the branches have been colored green to reflect this modified instability feature and the associated expected dynamics.

Besides the solutions that live on the main resonant curve shown in Fig. 5, there are also “isolas” of solutions, namely closed curves of solutions in the bifurcation diagram Isolas ; see Fig. 6 for example (and isola_strategy on how it was obtained). This isola lives entirely within the spectral gap. Periodic solutions with frequency that lie within a spectral gap are interesting since they are (necessarily) spatially localized. Such solutions are called nonlinear localized modes, or breathers Flach2007 . The case of damped-driven lattices where there is a driving source has attracted considerable attention in the present context. This is because of the intriguing interplay of localization due to the “defect” (i.e., the source) and that of the bulk breathers (or traveling waves) potentially available in the medium. This interplay can be used to produce hysteresis, metastability, resonances/antiresonances and related hallmarks of nonlinear dynamics that have been explored in a number of publications Nature11 ; hooge12 ; PhysRevE.92.063203 . It is interesting to highlight that there is a portion of the relevant isola branch that is found to be spectrally stable. The starting point of our subsequent study in this work will be to explore the time-periodic solutions of this system that have a frequency that lies in the spectral gap in the presence of damping and driving. Having examined those, we will subsequently turn to the case of the breathers in the Hamiltonian lattice.

We now fix the frequency (within the gap), and vary the amplitude of the drive, as shown in Fig. 7. We show the dependence of $\overline{E}$ as a function of $a_{d}$ for a number of cases. The left and middle panels of Fig. 7 correspond to the hardening and softening cases, respectively, both with $f_{d}=0.23$ where the well-known hysteretic bifurcation diagram is obtained; see, e.g., Nature11 ; hooge12 . It should be noted in passing that the numerical results depicted in these panels are shown for amplitudes of the drive of up to $a_{d}=0.25$ and $0.5$, respectively. The stability of the hysteretic bifurcation diagram is in line with earlier works Nature11 ; hooge12 . The fold bifurcations in the left panel of the figure occur at $a_{d}\approx 0.1961$ and $a_{d}\approx 0.1329$. In the middle panel the bifurcation values are $a_{d}\approx 0.2726$ and $a_{d}\approx 0.1767$. Notice, however, that there are distinctive features such as the short interval of oscillatory instability (for $a_{d}\approx[0.1768,0.1865])$) preceded by a very narrow segment (for $a_{d}\approx[0.1767,0.1768]$) of stable breathers in the middle panel of the softening case very near the bifurcation point. The right panel of the figure is a representative example for the alternating stiffness case (i.e., with nonlinearity exponents given by Eq. (2)) with $f_{d}=0.25$ lying in the frequency gap. The bifurcation curve (that is shown for drive amplitudes of up to $a_{d}=0.7$) becomes far more elaborate. It features a cascade of turning points, the first of which (i.e., at $a_{d}\approx 0.6935$) represents a fold bifurcation. Note that the labels (a)-(c) in the figure are connected with the numerical results of Fig. 8. Importantly, many of the numerous intermediate branches are found to be unstable through real instabilities, although the branch eventually features an oscillatory instability for sufficiently larger amplitude (as before). This additional complexity of the bifurcation diagram is introduced by means of the interplay between strain-softening and strain-hardening dynamics.

In Fig. 8 we show select examples of the dynamics of breather solutions (with no perturbation added on top of the respective initial conditions) for $f_{d}=0.25$, which is in the spectral gap. Those solutions correspond to values of $a_{d}=0.4$ and are associated with the labels of the right panel of Fig. 7. The left column of Fig. 8 shows the strain profiles, i.e., $y_{n}=u_{n-1}-u_{n}$ where the insets depict the Floquet spectra. The solution in the top row is spectrally stable (and remains accordingly robust, although the energy quantity we measure oscillates due to the presence of damping and driving) whereas the ones in the middle and bottom are unstable. The one in the middle row is an exponentially unstable waveform since the instability is characterized by a real Floquet multiplier pair (with a dominant unstable mode of $\lambda_{r}\approx 1.63576$), and the bottom one is oscillatory unstable corresponding to a complex multiplier quartet (with a dominant instability of $\lambda\approx 1.07981\pm 0.265122$). The stability results are corroborated in the respective panels of the middle column of the figure which present the spatio-temporal evolution of the energy density [cf. Eq. (18)]. It can be seen in the middle panels of the second and third rows of the figure that after an unstable breather gets destroyed, a wavepacket is created which propagates with decreasing amplitude toward the bulk of the lattice. An additional part of the relevant energy is localized at the left edge of the computational domain and appears to localize forming a smaller amplitude breathing structure therein. At longer times, these structures gradually asymptote to the stable solution of the top left panel of the figure. This finding can be further explored by looking at the right column of the figure. In particular, we present the total energy $E$ of Eq. (17) as a function of time in the respective cases. In the top right panel of the figure, the total energy of the stable solution oscillates, as discussed above, but maintains a nearly constant average of about 0.245 for $1000$ periods. On the other hand, in the middle and bottom panels, the energy experiences large deviations, highlighting the onset of the instability. Both unstable solutions eventually approach the stable one shown in the top row. This can be seen in the energy plots, where the average energy approaches 0.245 at $t\gtrapprox 1600$ for the unstable breather with label (b) and at $t\gtrapprox 1800$ for the unstable breather with label (c).

In all the numerical results that we have presented so far, $\gamma=8\times 10^{-3}$ was considered in Eq. (1). Next, we will study the effect of varying the dissipation parameter $\gamma$, where we gradually decrease it in order to approach the Hamiltonian limit. Moreover, this study will be beneficial while connecting the respective results presented herein with the ones of Sec. IV for $\gamma=0$. Let us start our discussion by considering Fig. 9 which corresponds to the value of driving frequency of $f_{d}\approx 0.2566$ being itself the cut-off frequency of the bottom optical edge of the optical band (see the right panel of Fig. 1). From the left to the right panels of the figure, the dependence of $\overline{E}$ is shown as a function of $a_{d}$ for $\gamma=8\times 10^{-3}$, $8\times 10^{-4}$, and $\gamma=8\times 10^{-5}$, respectively. We mention in passing the existence of a cascade of turning points in the panels, with the first four happening (as we trace the bifurcation curve from $a_{d}=0$ with increasing arclength) at $a_{d}\approx(0.33847,0.28153,0.50280,0.26376)$ (left), $a_{d}\approx(0.12421,0.03476,0.32090,0.03386)$ (middle), and $a_{d}\approx(0.12325,0.00351,0.32084,0.00339)$ (right). A movie illustrating how the solution profile changes as one moves along the bifurcation curve of the right panel is included in the supplement.

As the value of the dissipation parameter decreases, when $a_{d}\rightarrow 0$, Eq. (1) approaches its Hamiltonian sibling, i.e., when $(\gamma,a_{d})\rightarrow(0,0)$. Indeed, the solutions at the turning points in the right panel with $a_{d}\ll 1$ (and $\gamma\ll 1$), aside from the lowest one introduced by the drive and vanishing in this limit, are expected to be proximal to Hamiltonian breathers, since they are large amplitude (energy) states that exist in the lattice with nearly zero drive and zero damping. Recall that in the above diagrams we have fixed the frequency at the lower optical band edge. Bifurcation diagrams for frequencies deeper in the gap generally follow the same trend, namely that the turning points move toward smaller drive amplitude as the damping is decreased. However, the overall bifurcation structure is more intricate, as we explore next. Moreover, we would like to report an intriguing feature that happens at the turning points of the above mentioned bifurcation diagrams. As we approach a turning point, two eigenvalues are moving (one from the first and the other one from the fourth quadrant inside the unit circle) toward the real axis, until they collide exactly at the turning point, and past that, a pair of real multipliers emerges. We note that for the Hamiltonian problem at hand, such a collision happens at $1+0\mathrm{i}$ but for the present damped-driven problem, the collisions happen slightly inside the unit circle. More concretely, in the middle panel of Fig. 9, the eigenvalues at the first turning point (i.e., $a_{d}\approx 0.1242$) are $\lambda_{1}\approx 0.9992$ and $\lambda_{2}\approx 0.9976$. Interestingly and in line with similar earlier observations in jesus , they satisfy the relation $\lambda_{1}\lambda_{2}=e^{-\gamma T_{d}}$ which for the Hamiltonian case, i.e., $\gamma=0$, it gives $\lambda_{1}\lambda_{2}=1$. Indeed, we have $\lambda_{1}\lambda_{2}\approx 0.9968$ and $e^{-\gamma T_{d}}\approx 0.9968$. We further checked the relation $\lambda_{1}\lambda_{2}=e^{-\gamma T_{d}}$ at the other turning points and cases in $\gamma$ (right panel of the figure), and it holds indeed. It should also be noted that this implies that the stability change does not arise at the turning point but only nearby, given that the collision happens inside of the unit circle, and hence further parametric tuning is needed for one of the (now real) multipliers to cross the $1+0\mathrm{i}$ point.

The connection between the damped-driven and Hamiltonian breathers is explored in Fig. 10. We chose a frequency that is in the gap, namely $f_{d}=0.25$, and a small damping parameter $\gamma=8\times 10^{-4}$. At this frequency, the amplitude continuation (analogous to the one shown in Fig. 9) yields a more intricate bifurcation diagram (see the top left and middle panels in the figure). It is expected that a solution lying at a turning point with small driving amplitude (see the orange square marker in the inset of the top middle panel) will be proximal to an associated Hamiltonian breather, i.e., an intrinsic localized mode of Eq. (1) with $\gamma=0$. The velocity profile of the damped-driven breather (corresponding to the orange square maker of the top middle panel) is shown in blue in the top right panel of the figure. The solution itself compares fairly well with the Hamiltonian breather, which is shown in red.

In order to construct Hamiltonian breathers of Eq. (1) numerically, we consider a lattice with $N=201$ beads and employ zero Dirichlet and free boundary conditions at the endpoints of the chain, respectively, i.e., $u_{0}=0$ and $u_{N+1}=u_{N}$. The larger lattice size was needed in order to minimize the effect of the boundary conditions (since the strain of the breather solution is nearly zero at the boundaries for large lattices). The initial guess to our solvers was provided by the associated linear mode at the bottom edge of the optical band of the Hamiltonian problem (i.e., when $\gamma=0$). The Hamiltonian breathers found in this way are centered at the middle of the lattice. Thus, in order to compare the obtained breather to the damped-driven one (which is localized at the source of the drive, namely the left boundary) we translate the Hamiltonian breather so that the largest amplitude is on the boundary. This is how the comparison shown the the top right panel of Fig. 10 was made, showcasing the proximity between the two.

The insets of the top right panel of Fig. 10 indicate that the breathers of the damped-driven model and its Hamiltonian version are both exponentially unstable (with dominant unstable modes of $\lambda_{r}\approx 1.23721$ and $\lambda_{r}\approx 1.19449$, respectively), in addition to being extremely proximal to each other profile-wise. It should be reminded here that the difference in the density of Floquet multipliers populating the unit circle between the damped-driven (left) and Hamiltonian (right) insets stems from the different corresponding size lattices. This stability analysis is corroborated in the bottom left and middle panels of Fig. 10 where the spatio-temporal evolution of the energy density of (randomly) perturbed breathers is shown therein. It can be discerned from the figure that the instability in the damped-driven case manifests itself slightly earlier than in the Hamiltonian one, in line with our observation of the slightly larger (instability) growth rate. Moreover, the bottom right panel showcases the total energy as a function of time in both cases. Although the (averaged) energy in the damped-driven case (shown in blue therein) is bounded, it starts deviating at $t\approx 60$ from its original value signaling the onset of the instability. On the other hand, and as per the Hamiltonian nature of the breather, the energy is conserved (and shown in red), i.e., in this case, the energy is not a suitable diagnostic for discerning the instability, yet its conservation is indicative of the accuracy of the numerical simulation. On the other hand, if we probe the total energy of the first 5 sites, one can see that this local energy portion starts deviating around $t\approx 150$ from its initial (unstable equilibrium) value, thus signaling the onset of the instability, in line with the dynamical evolution of the left panel of Fig. 11. In both cases, over suitably long times, we observe that the resulting dynamical outcome is rather similar although the Hamiltonian system features longer-lived (and total-energy-preserving) dispersive radiation. In the dissipative case, the relevant energy appears to decay away due to the damped nature of the lattice bulk. While here we compared a Hamiltonian breather with damped-driven breather located at a single turning point in the bifurcation diagram, there is in fact a zoo of such solutions (one at each turning point in the bifurcation diagram), see the middle and right panels of Fig. 11. A more systematic study of all the relevant waveforms may be an interesting topic for further study.

We have just shown that when the damping is small, the Hamiltonian breathers are close to breather solutions of the damped-driven lattice. In the former case we can analytically approximate breather solutions by deriving an NLS equation to describe slow modulations in time and space of an underlying carrier wave. The NLS equation has an exact soliton solution, and thus one can obtain an analytical approximation in this way. The derivation of the NLS equation we present is similar to the one performed in Huang2 for a mass-dimer, but differs in a few key ways, which is why we include the derivation here.

At first, we re-write the equation of motion of Eq. (1) with $\gamma=0$ in terms of even ($v_{n}\coloneqq u_{2n}$) and odd ($w_{n}\coloneqq u_{2n+1}$) particles as

Recall that $p_{1}=p-\delta$, while $p_{2}=p+\delta$ above and in what follows. Upon assuming $|w_{n-1}-v_{n}|\ll 1$ and $|w_{n}-v_{n}|\ll 1$, Eqs. (19a)-(19b) reduce into

where $M_{n}$ ($n=1,2$) is defined by

with coefficients $J^{(1)}=-1$, $J^{(2)}_{n}=p_{n}$, $J^{(3)}_{n}=-p_{n}\left(p_{n}-1\right)/2$, and $J^{(4)}_{n}=p_{n}\left(p_{n}-1\right)\left(p_{n}-2\right)/6$.

Next, solutions in the form of fast oscillating yet small in amplitude patterns modulated by envelopes (that vary slowly in space and time) are sought by the following multiple scale Ansätze:

where $\varepsilon\ll 1$ as well as $X=\varepsilon\left(n-ct\right)$, $T=\varepsilon^{2}t$, and $E_{n}(t)=e^{i\left(kn+\omega t\right)}$. Note that $U_{s,j},W_{s,j}\in\mathbb{C}$, and also $U_{s,-j}=\overline{U_{s,j}}$ and $W_{s,-j}=\overline{W_{s,j}}$ with the bar standing for complex conjugation. Based on Eqs. (21a)-(21b), we also have that

together with the following Taylor expansions:

Since we are interested in standing breather solutions, we now set $k=\pi$ and hence $c=0$. Thus, the two possible frequencies from which we perturb are $\omega_{+}(\pi)$, corresponding to the the lower optical cut-off frequency, and $\omega_{-}(\pi)$, corresponding to the upper acoustic cut-off frequency.

We now plug Eqs. (21a)-(21b) into Eqs. (20a)-(20b), and collect coefficients of each $\mathcal{O}\left(\varepsilon^{s}E_{n}^{j}(t)\right)$ term. In particular, Eq. (20a) (after suppressing the $t$ dependence for $E_{n}$) gives, up to the first two orders:

Similarly, Eq. (20b) gives:

From Eqs. (24b) and (25b), we can express the underlying system as $\mathcal{M}\left[U_{1,1}\,\,W_{1,1}\right]^{T}=0$ with

which has a non-trivial solution provided that $\mathcal{M}$ has a vanishing determinant. It should be noted in passing that the matrix $\mathcal{M}$ is the same as the one of Eq. (12) with $\gamma=0$. The vanishing determinant of $\mathcal{M}$ yields the dispersion relation

where the “$+$” subscript corresponds to the lower cut-off of the optical band, and the “$-$” subscript corresponds to the upper cut-off of the acoustic band. In addition, from Eq. (25b) (same argument can be deduced by using Eq. (24b)), we obtain:

based on the dispersion relation.

Upon utilizing Eqs. (24a) and (25a), Eqs. (24d) and (25d) yield:

which (with the aid of Eq. (28)) can be cast into

Since the matrix $\mathcal{M}^{\left(\pm\right)}$ is singular (note that the superscripts emanate from $\omega_{\pm}(\pi))$, it follows that $\left(\mathcal{M}^{\left(\pm\right)}\right)^{T}\left(=\mathcal{M}^{\left(\pm\right)}\right)$ (where $T$ stands for its transpose) has a non-trivial one-dimensional kernel. Indeed, it is then a direct calculation to identify:

such that $\left(\mathcal{M}^{\left(\pm\right)}\right)^{T}\,\mathcal{W}^{\ast\left(\pm\right)}=\mathbf{0}$ (upon using Eq. (27) too). In addition, and upon considering the Fredholm alternative for Eq. (30), we can show that the latter indeed has a solution if and only if the right hand side of Eq. (30) is orthogonal to the co-kernel of $\mathcal{M}^{\left(\pm\right)}$ (that is, the kernel of $\left(\mathcal{M}^{\left(\pm\right)}\right)^{T}$).

We now consider Eqs. (24e) and (25e). With the aid of Eqs. (27) and (28), we arrive at

Next we consider the $\varepsilon^{2}\,E^{0}$ equations, and in particular Eq. (24c) which is solved with respect to $W_{2,0}$ (the exact same solution is obtained from Eq. (25c)), thus yielding:

Next, we focus on the $\varepsilon^{3}\,E^{0}$ equations for $k=\pi$ (shown in the Appendix B) from which we will obtain an equation for $U_{1,0}$ that depends on $U_{1,1}$ and its (complex) conjugate. Upon utilizing Eqs. (24a), (25a), and (28), one can find the $W_{3,0}^{\left(\pm\right)}$ component. Then, if one plugs the expression for $W_{2,0}^{\left(\pm\right)}$ [cf. Eq. (33)] as well as $W_{3,0}^{\left(\pm\right)}$ into the $\varepsilon^{3}\,E^{0}$ equation for $v_{n}$, and solves subsequently with respect to $\partial_{X}^{2}U_{1,0}^{\left(\pm\right)}$, one obtains:

Note that Eq. (34) can be integrated once with respect to $X$, thus yielding

where $f(T)$ is an arbitrary function of $T$. Considering our interest in structures that asymptotically vanish at $\infty$ we set $f(T)\equiv 0$. It should be noted in passing that when $\delta\equiv 0$ (i.e., the monomer case), the prefactor $\nu_{0}^{\left(\pm\right)}=4\left(p-1\right)$ is exactly the same as the one in granularBook (see, Eq. (A19) therein).

Finally, we focus on the $\varepsilon^{3}\,E^{1}$ equations for both $v$ and $w$, where we will find the NLS equation. At first, we express $W_{2,1}^{\left(\pm\right)}$ in terms of $U_{2,1}^{\left(\pm\right)}$:

which emanates from the (singular) system of Eq. (30). Next, we utilize Eqs. (24a), (25a), (28), (32), and Eq. (33). This way, we obtain the following system at $\mathcal{O}\left(\varepsilon^{3}\,E^{1}\right)$: