Supertransmission channel for an intrinsic localized mode in a one-dimensional nonlinear physical lattice

The combination of nonlinearity and discreteness makes possible localized vibrational modes in homogeneous lattices. Such an excitation, called an intrinsic localized mode (ILM),¹ or discrete breather (DB),^2,3 can exist in a lattice of any dimension. It was shown early on that for an ILM frequency close to the plane wave normal mode spectrum mobility was possible but when far removed from the band it was stationary stable and pinned at a lattice point.⁴ A Peierls-Nabarro (PN) potential barrier is often introduced to characterize the fact that an ILM located halfway between two lattice sites has a different energy than one centered on a lattice site. Starting with the work on kink dynamics⁵ in a discrete sine-Gordon system, effort has gone into trying to remove the PN barrier for localized excitations in various discrete 1D lattices both for practical^6–8 and theoretical applications.^9–12 For a zero velocity excitation the PN barrier has been successfully eliminated with specific parameters for 1D discrete equations of the sine-Gordon,⁵ Klein-Gordon (KG),^9–11 and nonlinear Schrödinger^13–15 forms. Slowly running ILMs have been shown to exist, in principle, for nonlinear intersite interactions, and for Klein-Gordon and nonlinear Schrödinger systems with specific potentials.^10,16 This has been followed by a broad effort showing that radiationless traveling localized modes can be produced with nonlinear Schrödinger equation forms^17–23 using specific conditions, such as stability inversion,²¹ saturable nonlinearity,^19,20,22 and a gain-loss mechanism.¹⁸ 

Another direction demonstrating localized excitation mobility makes use of the concept of moving embedded solitons (ESs), where the soliton internal frequency is inside a spectrum of linear modes. Although originally identified for continuous nonlinear optical model systems that support analytical solutions,^17,24–26 there has even been a family of stable ES identified where the analytical solutions have not been found.²⁷ More recently, localized excitations have been shown to exist in some discrete nonlinear models where they are referred to as embedded lattice solitons (ELS).^17,28 The general result is that for either an ES or an ELS to exist there can be no resonance between it and the small amplitude linear waves that are also present. Although much less is known about ELS, a specific example is the demonstration of a two-parameter family of an exact ELS for a discrete version of a complex modified Korteweg-de Vries equation that includes next-nearest-neighbor (NN) coupling.²⁸ 

The recent focus on the resonance interaction in discrete systems between the moving nonlinear localized excitation and the plane wave spectrum as the key element,^29,30 rather than the PN barrier itself, has led us to address the question how best to identify the resonance for strongly localized ILMs in discrete nonlinear physical systems where analytic solutions are unknown. We have found a practical procedure, whereby the resonance strength can be readily identified and then, depending on the relative contributions of the different nonlinearities, eliminated. This is demonstrated for an ILM traveling in a driven and damped 1D Newtonian lattice that contains both onsite and NN intersite nonlinear forces. The approach is to focus attention on the FT of the nonlinear force presented in a 2D $(k,ω)$ map. Although the driven-damped condition no longer characterizes a mathematically ideal lattice, this complication is compensated by precise measurements obtained for the traveling ILM, including its shape. Using these measured shapes simulations have been carried out for traveling ILMs under the same parameter conditions but without the driver or damper. The end result for both driven and undriven lattices is the discovery of simple specific conditions for a supertransmission channel for traveling ILMs. Such discrete equations describe a variety of physical systems ranging from electrical transmission lines,³¹ to micromechanical arrays,^32,33 to 1D antiferromagnets,³⁴ to nonlinear photonic crystal waveguides,³⁵ and to other systems treated within a tight binding formulation.³⁶ 

This paper is organized as follows. In Sec. II, we develop the expected features of the 2D $(k,ω)$ map, first for a simple linear traveling pulse in the continuous limit and then for two well known discrete soliton examples, namely, the Toda lattice and the Ablowitz-Ladik (AL) model. Section III begins with mobile ILMs in the driven and damped 1D physical lattice that contains both onsite and NN intersite nonlinear forces. By varying the relative strengths of these two contributions, the Klein-Gordon onsite nonlinear force model, the NN intersite nonlinear model, and the supertransmission model for both hard and soft nonlinearities can be readily studied. A summary and conclusions are presented in Section IV. The  Appendix follows.

For a display of the 2D FT map of interest consider the following 1D pulse traveling in the continuum limit:

where ω_c and k_c are the frequency and wavenumber of a carrier wave, and $Ã(x)$ is the envelope. The Fourier transform in space and time of this excitation is

Thus in $(k,ω)$ space the amplitude components of the traveling ILM, $| ψ̃(k,w) |$⁠, are along a line $ω=v(k−kc)+ωc$ centered at k_c with slope v.³⁷ 

For a discrete lattice, as displayed in Fig. 1, a localized moving $| ψ̃(k,w) |$ in the $(k,ω)$ space representation becomes a line centered at the carrier wavenumber k_c (identified by the vertical bar) and at the carrier frequency ω_c (the intersection of the bar and the line). Two different types of mobile ILMs are shown. Superimposed in this figure is the dispersion curve for a small amplitude KG lattice, represented by the dotted curve. For a traveling local mode that is sharply peaked, the FT components of $| ψ̃(k,w) |$ in the $(k,ω)$ range over a large extent in k-space and intersect the linear dispersion curve where the dashed and dots cross. We shall show that the nonlinear force component at this crossing accounts for the resonance between the plane wave modes and the mobile localized excitation.

For this special integrable lattice mobile localized excitations, first called “lattice solitons” by Toda,^38,39 have been well studied and do not interact with the small amplitude plane wave spectrum. The purpose here is to demonstrate in the $(k,ω)$ space representation that the FT of the nonlinear force of such an excitation does not intersect the plane wave dispersion curve. The equations of motion are characterized by

where n is the lattice site, x_n is the displacement, and a and b are nonlinear potential parameters. The linear dispersion curve $ω=2ab sin | k/2 |=ωBW sin | k/2 |$ is represented by the dotted curve in Fig. 2(a), where $ωBW=2ab$ is the bandwidth. The single soliton solution is

where $ω=ab sinh μ$⁠.

To quantify the nonlinear force for the Toda equation, we linearize the force term of Eq. (3), namely, $−ab{ 2rn−rn−1−rn+1 }$ and then subtract it from Eq. (3). The result is

where the right hand side now describes the nonlinear force. The next step is to carry out the FT on the time dependent displacement $xq(tp)$ of the moving lattice soliton

where the displacement (real quantity) is first transformed over the time axis into a one-side $(0∼ωmax)$ complex spectrum. Then, the FT of the resulting complex values along the space axis gives the two sided $(k=−π∼π)$ spectrum. The result is the $| x̃(k,ω) |$⁠, the lattice soliton amplitude map. To obtain the corresponding spectrum for the nonlinear force, we insert the same numerical data, $xq(tp)$⁠, into the right hand side of Eq. (5) to obtain the time dependent nonlinear force $fq(tp)$⁠. Equation (6) is then used to produce $| f̃(km,ωn) |$⁠.

Figure 2 presents the numerical illustration of these findings with the following parameters. The lattice dispersion curve is represented by the dots in Fig. 2(a) with $ab=ωBW2/4$ and for the soliton $μ=0.5$⁠. The soliton displacement was calculated for a lattice of N = 500 particles for a time duration $T=500/ωBW$ with a time step of $0.5/ωBW$ and resulting FT amplitude $| x̃(km,ωn) |$ is displayed using a gray (log) scale. The maximum FT frequency $ωFTmax=πωBW$ is well beyond the maximum linear frequency $ω=ωBW$⁠. The dark straight bands represent the 2D FT of the single soliton solution (darker means larger amplitude). Because the solution is a pulse, not an envelope, the line feature starts from origin $(k,ω)=(0,0)$⁠. The line is tangent to the dispersion curve at the origin so it is the only point where the soliton line and the linear dispersion curve coincide. Panel (b) shows the FT amplitude $| x̃(k,ω) |$ along that linear band. The maximum occurs at the origin. The 2D map for the nonlinear force $| f̃(k,ω) |$ is displayed in Fig. 2(c). A similar band structure to that shown in panel (a) appears for the nonlinear force but there is no longer amplitude at the origin. This feature is brought out more clearly in panel (d) where again the amplitude along the band vs wavenumber is plotted. The conclusion is that the nonlinear force spectrum is zero where the soliton spectrum and the linear dispersion curve coincide so there is no possibility of a resonance between the two kinds of excitations at finite frequency.

There is added value in considering the same sort of nonlinear force spectral approach for the AL lattice, which supports a different kind of linear dispersion curve and a different kind of traveling mode solution for the single lattice soliton.^40,41 The equations of motion of AL model are

where K determines the bandwidth of the linear spectrum and λ is the nonlinear parameter. The dispersion curve is $ω=2K(cos k−1)$ and thus $ω<0$ for K > 0. We consider $λ>0$ to provide solitons. The single soliton solution is

where k_c and ω_c are the wave number and frequency of the carrier part of the solution, and μ and v are width parameter and velocity, respectively. To analyze the solution, first we choose k_c and ω_c, then μ and v are calculated as follows:

Inserting these into Eq. (8), gives $xn(t)$ over the fixed time interval with the appropriate time step as was just discussed for the Toda lattice soliton.

To illustrate the soliton and nonlinear force spectral results, a specific set of parameters is now chosen; namely, $kc=−1250π, ωc=−0.0278ωBW$⁠, where the bandwidth in this case is $ωBW=4K$⁠. The dotted curve in Fig. 2(e) identifies the dispersion curve. The wavenumber region $k=−2π∼0$ is shown instead of $k=−π∼π$⁠, in order to display the entire structure of the traveling soliton. (In this case the frequency spreads both positively and negatively, because of the complex displacement in Eq. (8).) The FT of the traveling mode is again mapped using a gray scale, where darker means larger amplitude. The dark narrow band tangent to the dispersion curve in Fig. 2(e) identifies the soliton spectrum. Figure 2(f) shows the FT amplitude along that tangent line. The peak amplitude is at the carrier wave-number and it decreases smoothly away from that location. The right hand side of Eq. (7) is the nonlinear force. When the same procedure is carried out here as was described for the Toda lattice one finds the nonlinear force spectrum given in Fig. 2(g). Now the dark band fades away where it crosses the linear dispersion curve. Perhaps more convincing is the plot of the nonlinear force amplitude along this band shown in panel (h) where the sharp dent implies zero nonlinear amplitude at the crossing point.

The particular physical lattice system we have in mind is the 1D micromechanical array, which has both onsite and intersite forces. It has been demonstrated experimentally that both components can have nonlinear contributions and that ILMs can be readily produced.³² As mentioned in the introduction, there are a number of different kinds of physical systems that obey this kind of nonlinear, nonintegrable equation so its choice is not at all restrictive. To describe the general case appropriate to such a 1D physical system periodic boundary conditions are introduced for the nonlinear equations of motion for a mono-element lattice containing both onsite and NN intersite force terms

With Λ = 0, the linear dispersion curve is $ω2=ωO2+2ωI2(1−cos k)$⁠, already described by the dotted curve in Fig. 1. In this example, the onsite resonance frequency $ωO=1.45ωBW$⁠, the bandwidth frequency $ωBW=ωO2+4ωI2−ωO$⁠, and the normalized intersite frequency $ωI=0.99ωBW$ all have similar values. The $(±)$ sign of 5th and 6th terms of Eq. (11) describe positive (hard) and negative (soft) nonlinearities, respectively. Here g is the fractional strength of the nonlinear intersite force and $(1−g)$⁠, the onsite nonlinear strength. The factor η is used to balance the strength of the two nonlinearities at a particular k-value, e.g., the intersite nonlinearity is more effective than the onsite one when the carrier wave is close to the zone boundary $k=π$⁠. This factor is determined from $η=∑ xn4/∑ (xn−xn−1)4$ where $xn=x0 cos kcn$ is evaluated at the carrier wave number k_c. The right hand side of Eq. (11) is for a propagating wave driver, which excites only the target carrier wave mode. Three parameters are borrowed from an actual cantilever array study:³² the damping term $τ=8.75(ms)=2770ωBW$⁠, giving a quality factor Q = 6400; the nonlinear coefficient $Λ=6.5×1020(s−2m−2)$⁠; and the driver amplitude $α=1000 (m/s2)$⁠. The lattice size is N = 400.

Since there is no analytic solution for an ILM in this lattice, we need a different procedure than was used earlier to generate $xn(t)$ for the lattice soliton examples. This procedure is to produce a running ILM in the numerical simulations first by choosing k_c then finding the normal mode frequency ω_n at k_c. Next slowly increase the driver frequency from ω_n, which automatically increases the ILM amplitude: the larger the frequency difference from the normal mode, the larger the ILM amplitude and the narrower it is in real space. In steady state one obtains the time dependent displacement $xn(tp)$ for Eq. (6). It is first transformed over the time axis into one-side $(0∼ωmax)$ complex spectra. A long time interval $T=6280/ωBW$ is used for fine frequency resolution, and a small time step for a large maximum frequency $(ωmax=20ωBW)$⁠. Then, the FT of the resulting complex values along space axis gives the two sided $(k=−π∼π)$ spectrum.

The FT displacement spectrum along the tangent line can be calculated from the envelope with a 1D-FT as indicated by Eqs. (1) and (2). Similarly, the nonlinear force components on the 2D-FT map along the tangent line can be calculated from only the envelope by using a 1D-FT. The FT components of the onsite and intersite nonlinear forces represent real spectra as shown in the  Appendix and can cancel against one another when they have opposite sign.

The ILM example presented in Fig. 3(a) shows the resulting FT amplitude $| x̃(km,ωn) |$ as a function of $(k,ω)$ for the KG model with onsite nonlinear coupling (g = 0). Here $kc=144π/200=2.2619$ is one of available k values for the N = 400 lattice. As expected, the FT amplitude band is centered on a line tangent to the linear dispersion curve. The arrow identifies where a resonance occurs between the FT amplitude components and normal modes of linear dispersion curve. Figure 3(b) brings out more clearly the FT amplitude along the center of the tangent band. A strong resonance is apparent. The traveling ILM looses energy by generating waves in the normal mode spectrum. The broad amplitude peak on the right hand side centered at $kc=2.2619$ is associated with the traveling ILM. The sharp depression centered at k_c (arrow) is due to negative response of the lattice away from the ILM spatial region. Such canceling has been reported for a stationary ILM.⁴² (The noisy structure for k < 0 is produced by the large envelope signal at higher frequency, which is continued from $k=π$⁠, and due to the finite FT time.) Figure 3(c) shows the 2D FT amplitude spectrum of the nonlinear force $−Λxn3$ calculated from x_n in the equations of motion. The amplitude spectrum of the nonlinear force intersects the normal mode spectrum but frame (d) shows that only a monotonic dependence occurs in this region.

For the pure NN intersite nonlinear lattice (g = 1.0), the 2D FT amplitude results for the ILM presented in Figs. 3(e)–3(h) are somewhat different. Again a resonance occurs in frame (e) where the FT amplitude along the tangent band intersects the plane wave dispersion curve but nearby there is a hole in this amplitude band. The FT amplitude along the center of the tangent band displayed in frame (f) shows that both a resonance and antiresonance are evident. We now turn our attention to the properties of the antiresonance. Frame (g) presents the 2D FT amplitude spectrum of the nonlinear force $−Λη{ (xi−xi+1)3+(xi−xi−1)3 }$ calculated from x_n in the equations of motion. The arrow identifies the location of the missing amplitude components. Again looking at the FT amplitude along the center of the tangent band, frame (h), it is evident that the strong antiresonance is a feature of the nonlinear force.

Qualitatively new results are obtained if both NN intersite and onsite nonlinear forces are included. Frames (i) and (j) demonstrate that by tuning the mixing ratio between the two nonlinear force terms in Eq. (11) to g = 0.924 the amplitude resonance shown in frames (a) and (b) and the resonance and antiresonance shown in frames (e) and (f) are eliminated. How this occurs is illustrated in frames (k) and (l). Changing the mixing ratio of the two kinds of nonlinear force terms shifts the hole position of the nonlinear force components displayed in frame (g) to the position shown in frame (h). For this specific set of parameters in Eq. (11), a supertransmission channel has appeared since there is no emission of plane waves produced by the propagating ILM.

Supertransmission parameters for the ILMs with different sets of k_c and ω_c, covering ILMs located both outside the plane wave spectrum (⁠$1.45ωBW∼2.45ωBW$⁠) as well as inside, are given in Table I. Because of k_c and ω_c dependencies of the resonant crossing point, the mixing ratio g must be tuned for these positive (hard) nonlinear cases $ωc>ωn$⁠. Examining this table shows that a smaller g is required for a larger k_c for this positive nonlinear case. Tuning is accomplished by adjusting the onsite term, because of the small η (high effectiveness of intersite nonlinear term) for these k_c. In all the cases, the mode shape in real space can be described by a carrier wave times an envelope. The last column in Table I presents the real space-time average full width at half maximum (FWHM) for each case. As expected, the larger the difference in frequency between the ILM carrier and the plane wave normal mode from which it originated, the sharper the real space pattern.

To further test the specific design parameters, which produce zero nonlinear force components at the resonant crossing point and hence greatly reduced radiation loss by the moving ILM, additional simulations have been carried out for traveling ILMs with a variation of Eq. (11), now without the driver or damper. Figure 4(a) presents ILMs velocities versus time for several k_c conditions, each with two different values of the positive nonlinear onsite and intersite strengths: g = 1 (pure intersite nonlinear model) and the corresponding supertransmission tuned g value. Three different k_c values are shown for the positive nonlinear lattice with N = 3200. Each t = 0 initial condition is obtained from the steady state driven-damped simulations with N = 400 in Eq. (11). The other 2800 lattice points are filled with zeros. The deceleration of the ILMs observed for each of the g = 1 cases is due to the radiation loss related to the presence of the resonance while for tuned g cases no deceleration is observed. This figure shows that by tuning the mixing ratio of the onsite and intersite nonlinearities to a specific combination for an ILM with a given $(kc,ωc)$⁠, one can greatly reduce the radiation loss. Since we obtained fairly constant velocity by launching an N = 400 lattice ILM into an N = 3200 lattice, these results demonstrate that this super transmission channel doesn't require background excitations, which may occur in the initial condition for the ILM obtained in the driven-damped simulations.

Figures 5(a) and 5(b) show the real space energy density plot for an ILM traveling either in the supertransmission channel (a) or the corresponding lossy case (b) for a specific k-value. A gray scale with a magnification of 100 is used to bring out the differences between the two cases. For the supertransmission lattice, the ILM travels at a constant speed without emitting small amplitude plane waves while for the nonlinear NN lattice, plane wave emission decreases the ILM speed. Note the plane wave emission evident in (b) is in the opposite direction from that of the ILM and the accumulated waves re-form into a weak localized backward directed wave packet represented by the relatively sharp gray band.

For the negative (soft) nonlinear case (g = 0) describes an ILM below the bottom of the plane wave spectrum (⁠$ωc<ωn$⁠) for an untuned, regular KG lattice. Here tuning is accomplished by increasing the intersite term, because of the very low effectiveness of the intersite nonlinear component. Figure 4(b) shows the velocity versus time results for a soft nonlinear lattice with $kc=0.2513$⁠, g = 0.0038. The parameters are given in the last row of Table I. Now the onsite nonlinearity is much more effective than the intersite nonlinearity since k_c is close to zero. In comparing the unturned case with the tuned case, note that the initial velocity for the regular lattice is larger than for the tuned case because the added negative nonlinear intersite term of the tuned one reduces its ILM velocity. The much smaller deceleration observed here for the untuned case than for the other cases shown in Fig. 4(a) is partially because of a smaller ILM amplitude, and partially because of the relatively remote crossing point at $Δk∼4.7$ compared to the cases shown in Fig. 4(a) where $Δk=2.3∼3.8$⁠. It appears the resonance effect is more important for a fast moving ILM like that shown in Fig. 1(a) than a slow moving one like in Fig. 1(b) because for the former the intersection point of the nonlinear force Fourier component is closer to k_c, the carrier point. It is also true that if the carrier frequency of the moving ILM is far enough from the band, the intersection point will occur in a much higher order Brillouin zone, producing a much weaker resonance effect so that the Peierls-Nabarro barrier would now play an important role in determining ILM mobility.

Recipe for production of supertransmission channel. Consider the positive nonlinear case. (1) With driver and damping in the equations of motion (Eq. (11)), choose k_c and find the normal mode frequency ω_n. (Initial conditions do not matter.) (2) Increase driver frequency from ω_n to produce locked ILM. Use the minimum α necessary to reach the high amplitude ILM state and pass the low amplitude chaotic state quickly. In such a driven-damped condition, the nonlinear parameter Λ can have any value >0. (3) Start from g = 1 (nonlinear NN model). (4) After the driven mode reaches steady state (several 10 thousand periods) take the double FT of the time dependent displacement giving the two-sided $(k,ω)$ spectrum. (5) Calculate the FT displacement spectrum along the tangent line. (6) Calculate the nonlinear force FT components along the tangent line (remember the  Appendix). (7) Check spectra in intersection region for that g value. (8) Change g and do simulation again. (9) Repeat until zero nonlinear force is achieved.

As long as mobile ILMs can be generated and the resonance intersection occurs at one location on the $(k,ω)$ map the range of parameter tunability to produce supertransmission extends over most of the plane wave dispersion curve shown in Fig. 3(i). For this particular dispersion curve producing supermobile ILMs in or above the top of the plane wave spectrum is straightforward but generating those below the bottom require special care. For cases near $kc∼π$ or $kc∼0$ supermobility is not possible because the velocity is too small and multiple resonances occur in k space. The key element for supertransmission is the presence of at least two nonlinearities, which have opposite sign where the intersection occurs in the $(k,ω)$ map, while the sum of these nonlinearities is large enough at k = k_c to produce the traveling ILM. This occurred for the physical equations considered here since the FT components of the nonlinear onsite force always have the same sign while a region of those of the NN nonlinear intersite force have opposite sign. [Note: A somewhat less physical combination we explored briefly that produced a supermobile ILM involved NN nonlinear forces together with 3rd NN nonlinear forces with a strength ratio of 1.0:0.5.] It should be noted that most of the physical parameters used in this paper were taken from the existing microcantilever arrays but the supertransmission mixing ratios displayed here would require a larger onsite nonlinearity; hence for an experimental demonstration a redesign would be required, such as widening the head of each cantilever.

By simulating traveling ILMs in the nonintegrable 1D physical lattices with driver and damping combinations, and either hard or soft onsite and intersite nonlinearities, we have found a way to suppress the resonant interaction between the nonlinear force and the plane wave modes. Examining a running ILM, with internal frequency located either outside or inside the plane wave spectrum, in a $(k,ω)$ map can be used to display the ILM amplitude and nonlinear force components centered on a band tangent to the normal mode dispersion curve. For both KG and NN nonlinear models, we have demonstrated that the signature of a retarded ILM is characterized by a resonance in $(k,ω)$ space where the plane wave dispersion curve and the FT amplitude components of the ILM intersect. It was discovered that for the NN case nearby the resonance location there is an antiresonance in this FT amplitude spectrum. We have shown that for situations where both onsite and intersite nonlinearities occur, the position of this hole in the nonlinear force spectrum depends on the specific mix of the two nonlinearities. By changing this mix, one can move the nonlinear force hole to the amplitude intersection point, thereby removing the retardation effect for a traveling ILM characterized by a specific frequency and k-vector. Thus for both driven and un-driven 1D nonlinear physical lattices with coupling of the kind shown in Eq. (11) simple specific conditions have been found for a supertransmission channel for a traveling ILM. Similar supertransmission channel features are to be anticipated for running ILMs in the 1D physical systems with more complex combinations of onsite and intersite nonlinearities. Because discrete equations of the form studied here characterize a variety of nonlinear transmission lines and other physical systems treated within a tight binding approximation such an engineered, intrinsic, low loss transmission channel is expected to be a very useful property.

M.S. was supported by JSPS-Grant-in-Aid for Scientific Research No. 25400394. A.J.S. was supported by Grant No. NSF-DMR-0906491 and he acknowledges the hospitality of the Department of Physics and Astronomy, University of Denver, where much of this work was completed.

The zero nonlinear force contribution at the intersection with the linear dispersion curve is achieved by summing two force components, the nonlinear onsite and intersite terms. To see how this follows, it is not necessary to consider the entire nonlinear force spectrum but only to examine the 2D complex amplitude components along the tangent line in the $(k,ω)$ map. Both nonlinear forces are generated from the three times multiplication of x_n, which is the convolution in reciprocal space. The signal of interest is centered around k_c, not $3kc$ so the FT of the onsite term on the tangent where $ω=v(k−kc)+ωc$ is

where $x̃(k)$ is the FT of the envelope part of x_n and $*$ indicates the convolution. In our case $x̃(k)$ only has a real part because of the symmetry of the envelope so that Eq. (A1) is also real.

A corresponding expression for the intersite term is

where we have made use of the fact that

The end result is a complex spectrum. The other intersite term is

and when Eq. (A2) is added to Eq. (A4) the total intersite spectrum is real. Thus, both onsite and intersite nonlinear force spectra are real.