Brillouin-zone definition in non-reciprocal Willis monatomic lattices

Brillouin-zone (BZ) definition in a class of non-reciprocal Willis monatomic lattices (WMLs) is analytically quantified. It is shown that BZ boundaries only shift in response to non-reciprocity in one-dimensional WMLs, implying a constant BZ width, with asymmetric dispersion diagrams exhibiting unequal wavenumber ranges for forward and backward going waves. An extension to square WMLs is briefly discussed, analogously demonstrating the emergence of shifted and irregularly shaped BZs, which maintain constant areas regardless of non-reciprocity strength.


Introduction
Exploration of new wave phenomena has been a driving motivation in wave propagation research [1], with a growing interest in elastodynamic non-reciprocity [2].Within linear, time-invariant elastic media, a wave emanating from a source towards a receiver must not behave differently if the locations of the source and receiver are interchanged, thanks to elastodynamic reciprocity.Designs of artificial materials have enabled new avenues to manipulate waves beyond what is possible in natural materials, and breakage of wave reciprocity is no exception.Given its potential in enabling novel applications (e.g., one-way vibrational isolators [3]), researchers have introduced a variety of structural designs tailored for non-reciprocity, such as exploiting nonlinear instability with supra-transmission phenomenon [4].Comparably, granular crystals, with their nonlinear and bifurcating-chaotic behavior, have been exploited to asymmetrically transmit energy, with potential application in sensing and energy harvesting [5].Another choice to achieve non-reciprocity is by spatiotemporally modulating mechanical properties, with artificial momentum bias forcing waves to propagate differently in the forward and backward directions [3,[6][7][8][9][10].
When elastic media defy elastodynamic reciprocity, their dispersion diagram becomes skewed and asymmetric about a zero wavenumber, and conventional Brillouin-zone (BZ) definitions are rendered obsolete.Cassedy and Oliner are amongst the first contributors to study the skewed (non-reciprocal) nature of BZ in an attempt to generalize its definition for wave propagation in space-time modulated electrical circuits [11].Unlike spatial-only modulation with dispersion diagrams perfectly repeating along the wavenumber axis, the presence of time modulation creates dispersiondiagram periodicity along a line with a slope that is a function of such modulation (under small modulation amplitude assumption).However, the adequacy of traditional BZs to fully capture wave non-reciprocity remains an ongoing concern, especially for two-dimensional systems [8].
In this effort, BZ definition is unraveled for a class of non-reciprocal monatomic lattices (MLs), synthesized from a non-reciprocal wave equation governing longitudinal waves in an axially-moving elastic (continuum) rod as follows [12]: where ρ, E, and u(x, t) denote the density, modulus of elasticity, and the longitudinal displacement as a function of space x and time t of the rod, respectively.The rod's constant speed is defined as v 0 = βc, where c = E/ρ is the elasticmedium's sonic speed and β ∈ [−1, 1] (+ (−) sign indicates forward (backward) motion) is the non-dimensional relative moving velocity [12].As a matter of fact, the nature of non-reciprocity of Eq. ( 1) is intertwined with that of elastic rods with spatiotemporally modulated properties.At the low-frequency limit and slow modulation, the macroscopic behavior of an elastic rod with a spatiotemporal wave-like modulation is governed by a Willis-type equation [7], which is identical in form to that of Eq. ( 1) for moving rods [12].Not only that Eq. ( 1) is of Willis type, but it is also gyroscopic as it can be re-casted as [13]:  1).The lattice's masses m are evenly spaced with a lattice constant a and the displacement of the i th mass is donated as u i .As seen from the figure, two types of springs exist: (i) a conventional spring k arising from the elasticity of the rod and (ii) a negative spring −m(v 0 /a) 2 induced as a consequence of the rod's motion.Non-local coupling parameters that are proportional to the velocity of the next neighbors of an i th unit cell break the lattice's reciprocity (shown for the i th unit cell only). where x are the mass, gyral, and stiffness operators.Whether subjected to spatiotemporal modulation or motion at a constant speed, elastic rods obeying Eq. ( 1) exhibit linear momentum bias (inducing non-reciprocity), as evident from the mixed derivative term.The ramifications of the induced non-reciprocity on BZ shall be first established for the lattice version of Eq. ( 1) and then extended to a two-dimensional counterpart.

Mathematical model
To find an equivalent ML to that of the moving rod, henceforth referred to as Willis monatomic lattices (WMLs), the spatial derivatives in Eq. ( 1) are discretized in the spatial domain x using the central differencing scheme, resulting in where a is the finite difference spacing and constitutes the lattice's constant (Figure 1).The spring constant and mass of WML are defined as the effective stiffness k = EA/a and mass m = ρAa of a rod segment of length a and area A.
Using the introduced parameters, a few mathematical manipulations result in the equation of motion for the i th unit cell of WML: Examining Eq. ( 4) and Figure 1, it is first deduced that the motion of individual masses m is physically coupled through springs k, as in typical MLs [14,15].Having the momentum bias introduced by the rod's constant speed, the physics of the lattice's motion is influenced in two ways: 1.Besides the spring k, the i th unit-cell displacement is coupled with its next neighbors through negative springs of coupling coefficient −mv 2 0 /a 2 , which remains negative regardless of the sign of v 0 (or β).Consequently, the effective stiffness (k e = k − mv 2 0 /a 2 ) of the WML reduces as the rod's speed v 0 increases, jeopardizing dynamical stability when mv 2 0 /a 2 > k (equivalent to |β| > 1) as a consequence of its effective stiffness k e becoming negative.
2. Non-local coupling terms emerge, which are proportional to the velocity of next-neighboring masses of the i th unit cell.Such non-local interactions (related to Willis coupling [16]) arise from discretizing the mixed derivative in Eq. ( 1) and dictate the strength of lattice's non-reciprocity.While the physical motion of the elastic medium enables such non-local coupling, it may be alternatively achieved via feedback control [17][18][19].Lastly, it is noteworthy that such non-local (positive-negative) couplings are akin to skew-symmetric gyroscopic coupling studied in literature [20][21][22][23][24][25].

Dispersion relation
Assuming harmonic motion u i = ûi e iωt and applying Bloch boundary conditions u i±1 = e ±iq u i , where q is the non-dimensional wavenumber, ω is the excitation frequency, and i is the imaginary unit, a non-dimensional dispersion relation is derived from Eq. (4): Here, the definition of the non-dimensional frequency is Ω = ω/ω 0 , where ω 0 = √ k/m.Note that v 0 = βc ≡ aβω 0 .The dispersion relation in Eq. ( 5) has a single (positive-frequency) dispersion branch, which reads: The term sin(q) is what makes the dispersion relation in (6) non-reciprocal as its sign changes simultaneously with q.Corroborating our earlier observation regarding the lattice's dynamical stability with the speed v 0 , Eq. ( 6) indicates that |β| > 1 yields complex values of the frequency Ω, signaling dynamical instability.A positive (negative) β results in a positive (negative) shift ϕ, as seen from the bias towards the forward-going (backward-going) waves, while β = 0 recovers the dispersion diagram of a reciprocal ML.However, the free-wave dispersion does not capture the shifted nature of BZ, showing the driven-wave approach significance in unraveling the BZ shift ϕ.Observe that the attenuation (represented by the imaginary component of the wavenumber in the driven-wave dispersion) is smaller in WMLs relative to its reciprocal ML counterpart (β = 0) when Ω > 2, i.e., a frequency higher than the lattice's cutoff frequency Ω max = 2. Also, Ω > 2 results in a non-constant real component of the wavenumber in WML cases only.Note that, for β = ±0.5 shown here, the corresponding BZ shift is ϕ ≈ 0.3π.
(bottom) Spatiotemporal FFT for the time response of WML for an impulse excitation in the left and right ends, quantifying forwardand backward-going wavenumbers, respectively, and demonstrating excellent agreement with the new BZ definition in the middle panel.

Brillouin-zone definition
Conventionally, the dispersion relation is plotted by sweeping the non-dimensional wavenumber within the first BZ of the range q ∈ [−π, π] and solve for Ω in Eq. ( 6) for each value of q, a method known as free-wave dispersion relation.The results of this process are shown in the top row of Figure 2 for different values of β, where β ̸ = 0 induces dispersion asymmetry about zero wavenumber (i.e., q = 0).Note that β = 0 recovers the dispersion diagram of a classical (and reciprocal) ML [15] and results in a perfectly symmetric dispersion about the center of BZ (q = 0).However, the conventional BZ does not accurately capture the range of wavenumbers spanning the forward-going and backward-going waves for β ̸ = 0. To show that, Eq. ( 5) can be reformulated to solve for the wavenumber q by having the frequency Ω as an input, resulting in the driven-wave dispersion relation.After rewriting the trigonometric terms in their exponential form, one can show that: Observe that the first and last terms in Eq. ( 7) are complex conjugates, which can be rewritten as |z|e ±i(q−q s ) , where The emergent phase shift q s is a function of both Ω and β and shall dictate the new definition of BZ.Following the parametrization in Eq. ( 8), the dispersion relation in ( 7) is recast to: and the driven-wave formulation of the dispersion relation is obtained by solving for q: For β ̸ = 0, Eq. ( 10) output is no longer bounded by the traditional BZ boundaries q ∈ [−π, π], evincing the shifted nature of BZ due to the phase shift q s .Owing to the group velocity vanishing at BZ boundaries for periodic media [26], WML's group velocity must be first derived to precisely pinpoint its BZ boundaries.Then, the BZ shift shall be obtained by evaluating q s at the frequency corresponding to a zero group velocity.To do so, consider the (non-dimensional) group velocity c g for WML, derived by taking the derivative of Eq. ( 6) with respect to the wavenumber q (See Supplementary Note 1 for more discussion on group velocity and phase shift): By equating the group velocity in Eq. ( 11) to zero, the following wavenumber roots are found: Substituting Eq. ( 12) into Eq.( 6) reveals that the frequency at a zero group velocity is Ω max = 2, which is constant, independent of β, and equivalent to the cutoff frequency of a reciprocal ML.Finally, the BZ shift, denoted as ϕ and constitutes the key result in this study, is derived by evaluating q s at Ω max = 2: Equations ( 10) and (13) uncover that the BZ in its entirety is shifted and it is within the range q ∈ [−π + ϕ, π + ϕ], with its width remaining constant at 2π in analogy to a reciprocal ML.Depending on the sign of β (and subsequently ϕ), the BZ moves forward or backward, as seen in the middle row of Figure 2.Not only that non-reciprocity shifts the BZ, but it also forces the range of the wavenumber corresponding to the forward-going and backward-going waves to be unequal.To maintain a constant BZ width of 2π, the amount of shrinkage (or extension) in the forward-going wave wavenumber range (i.e., q ∈ [0, π + ϕ]) is compensated by a larger (smaller) wavenumber range of the backward-going waves (i.e., q ∈ [−π + ϕ, 0]), depending on the sign of ϕ.This shift ϕ disappears when the relative speed β is zeroed out as expected, which is verifiable from Eq. ( 13), and the traditional BZ range (q ∈ [−π, π]) is recovered.
A few additional observations from the driven-wave dispersion in Figure 2: (i) the attenuation at a frequency higher than the cutoff frequency Ω max = 2 is smaller in WMLs compared to a reciprocal ML. (ii) The real component of the wavenumber with Ω max > 2 in WMLs does not have a constant value and varies as the frequency increases, unlike its reciprocal counterpart with a constant real wavenumber of ±π within the attenuation zone (Ω max > 2).(iii) The BZ defined in MLs with β = 0 is symmetric about its center, allowing for an irreducible BZ to be defined within the range q ∈ [0, π].As WMLs have asymmetric dispersion branches, reducing the first BZ may be elusive.A zero wavenumber, however, remains the point that divides the forward-and backward-going waves for both lattices.Finally, the inequivalence of the forward and backward wavenumber ranges is numerically verified via a spatiotemporal fast-Fourier transform (FFT) of the time response of a finite WML (bottom row of Figure 2), exhibiting excellent agreement with the analytical results of the newly defined BZ definition in the middle row (See Supplementary Notes 2 for more details).

Extension to two-dimensional lattices
Next, the analysis for the one-dimensional WML is extended to a two-dimensional counterpart, using an elastic membrane with velocity modulations v x and v y in the x and y directions, respectively, and a proposed govering equation of: where u(x, y, t) symbolizes the membrane's transverse displacement.After discretizing Eq. ( 14) as a square grid, it can be shown that the dispersion relation for the resulting square WML is (See Supplementary Note 3): where q x,y and β x,y = v x,y /c are the wavenumber and relative velocity in the x (or y) direction, respectively.Note that β x,y ∈ [−1, 1] and the dispersion relation of a traditional square ML is recovered if β x = β y = 0, as expected.The dispersion surface Ω of the square WML in Eq. ( 15) is found by solving the quadratic equation for values of q x and q y within a newly defined (and shifted) BZ.To reveal such a BZ, an identical procedure to the one used in deriving Eq. ( 9) is followed here, yielding an alternative form of Eq. ( 15): where Numerically, it can be shown that the dispersion relation of the square WML has a constant cutoff frequency at Ω max = 2, regardless of β x and β y (See Supplementary Figure S3(a)).Analogous to Eq. ( 13), and after plugging in Ω max = 2 back into q x,y s in Eq. ( 17), it is concluded that the two-dimensional BZ corners are shifted with x and y direction shifts of: In a reciprocal square ML, the vertical (horizontal) boundaries of its perfectly square BZ correspond to a vanishing group velocity in the x-direction (y-direction) only.In extension, the same condition shall be imposed on the nonreciprocal WML case to find the horizontal and vertical boundaries of the shifted BZ.By doing so, it is found that the BZ for WML is no longer a square and has irregularly shaped boundaries.Yet, its area remains a constant value of 4π 2 , identical to its reciprocal case, regardless of the magnitude of β x and β y (See Figure 3(a,b) and Supplementary Figure S3(b)).
Another way to depict the two-dimensional dispersion relation is by using directivity plots and reducing the two variables q x and q y into a single (direction) angle equal to ψ = tan −1 (q y /q x ).In the reciprocal case, all four quadrants of the angle ψ show identical profiles of the dispersion surface, and the same cannot be said when β x ̸ = 0 and/or β y ̸ = 0 (Figure 3(c)).Complete dispersion surfaces in the directivity plot are generated if the wavenumber values within the non-reciprocal BZ boundaries (shown in Figure 3(b)) are substituted to the dispersion relation, which is not the case if the traditional BZ values of the wavenumbers q x and q y are used instead (See Figure 3(c) and Supplementary Figure S4).ψ = tan −1 (q y /q x ) (0,0) ψ (q x ,q y ) q y q x spatiotemporal FFT Finally, a spatiotemporal FFT of the time response of a square WML of finite size is shown in Figure 3(d) for β x = β y = 0 (reciprocal case) and β x = 0.5 and β y = 0.75 (non-reciprocal case).The FFT contours exhibit a close agreement with their analytical counterparts in Figure 3(c), substantiating the approach's validity and the completeness of the newly defined BZ (See Supplementary Note 2 for more details).

Concluding remarks
In summary, this study demonstrates that the Brillouin zone (BZ) remains constant in width/area with non-reciprocity in a class of Willis monatomic lattices (WMLs), yet its boundaries are shifted in a precisely quantifiable amount.Synthesized from a modified wave equation of longitudinal waves in moving elastic rod [12], the one-dimensional WMLs is studied analytically and its driven-wave dispersion relation is proven vital for the quantification of the BZ shift ϕ.It is shown that the proposed theory perfectly agrees with the numerical simulation, validating BZ shifting, its constant width, and inequivalence of wavenumber ranges occupying the forward-going and backward going waves.It is also established that an increase (decrease) in the forward-going wavenumber range is compensated by a shrinkage (enlargement) in the backward-going one.An extension to a square WML is also established, demonstrating the shifted nature of the BZ and its area being constant and unaffected by the degree of non-reciprocity, dictated by the velocity modulation in the x and y direction.Additionally, it is shown that the two-dimensional BZ is no longer a square, and its deformed vertical (horizontal) boundaries are found by setting the group velocity in the x (y) propagation direction to zero.The established clarification on the definition of BZs is envisioned to be a stepping stone for further investigations in BZ quantification for different types of modulations for wave non-reciprocity, as well as for periodic media with multi-mode dispersion relations.

Note 1. Group velocity and phase shift in one-dimensional Willis monatomic lattice
For completeness, the group velocity c g , the phase shift q s , and the Brillouin zone shift (corresponds to q s at the cutoff frequency of Ω max = 2) are shown in Figure S1 for different values of β.As proven earlier, the group velocity (shown here for the forward-going waves q ∈ [0, π + ϕ]) vanishes at the cutoff frequency Ω max = 2 (corresponding to q = π + ϕ) regardless of the magnitude of the relative velocity β.Worthy of note is that the initial slope of the dispersion branch (i.e., the group velocity at Ω = 0) is c 0 = 1 − β (c 0 = −1 − β) for dispersion branches pertaining to forward-going (backward-going) waves.For the phase shift q s , it is observed that its range is within [−π/2, π/2], corresponding to a range of relative velocity of β ∈ [−1, 1] and frequency Ω ∈ [0, 2].Being frequency-dependent is another clear characteristic of phase shift q s as seen from the series of curves in Figure S1 with β ̸ = 0, where q s appears to increase as the frequency increases.The Brillouin zone shift ϕ, on the other hand, does not depend on Ω since it corresponds to the phase shift q s at the constant cutoff frequency Ω max = 2.Note that the output of the group velocity c g and phase shift q s also depends on the excitation frequency Ω.The shift of the Brillouin zone ϕ signifies the values of phase shift q s at a frequency of Ω max = 2, and, thus is only a function of β.

Note 2. Numerical validation for the dispersion relations of Willis monatomic lattices
Numerical simulations are performed using a finite lattice of fixed boundaries and sufficient size to support the theoretical results and further validate the shift of the boundaries of the Brillouin zone (BZ).For the one-dimensional case, and to accurately quantifying the shift in BZ, the wavenumbers corresponding the forward-going (backward going) waves must be excited separately.Therefore, an impulse excitation is injected to a finite lattice of 151 masses at its left (right) end and the simulation time is chosen such that the impulse does not reach the other end and reflects (See Supplementary Videos S1-S6 for the finite lattice response for β = ±0.5 and β = 0).A spatiotemporal fast Fourier transform (FFT) is then performed to construct the numerical dispersion diagram from the impulse response of the lattice.The dispersion branch for the forward-going (backward-going) waves are constructed from the numerical response of the lattice for the impulse excitation on the left (right) end, which guarantees that the wavenumbers corresponding to forward-going and backward-going waves can be quantified separately.As seen in Figure S2, the wavenumber spanning the forward-going (backward-going) waves is smaller (larger) than the conventional range in a reciprocal media when β < 0 (β > 0), in agreement with the theoretical results (shown as white dashed-lines).It is now evident that the numerical simulation corroborates that BZ is no longer in the conventional range and its shift is precisely quantifiable based on the phase shift ϕ.
For the square lattices, an initial displacement was applied on the middle of a 40 × 40 lattice to excite all modes of the structure, and the simulation is terminated before reflections from the boundaries occur.Performing spatiotemporal FFT and projecting its amplitude on a directivity plot, the result is depicted in Figure 3(d) in the main manuscript.The wavenumbers used for this procedure are those for the newly proposed BZ, resulting in dispersion plots that are in excellent agreement with the analytically obtained one in Figure 3(c) in the main manuscript.For better results visualization, FFT's magnitude of 10% or smaller was filtered.

Note 3. Equation of motion and dispersion relation of square Willis monatomic lattices
The proposed equation of motion for a two-dimensional, non-reciprocal continuous elastic medium is: which correspond to the following mass, gyral, and stiffness operators of Eq. ( 2) in the main manuscript: Using central finite difference on the spatial domains, the discretization for the mixed derivatives, as well as the x and y second derivatives are as follows: where [˙] represents the derivative in time, t.Assuming a square grid, i.e., a x = a y = a, and using the discretization schemes shown in Eq. (S3), one can show that the equation of motion for the (i, j) th unit cell is given by: Knowing that ω 0 = c/a and using an equivalent mass of m per unit cell, the equation of motion (S4) can be re-written in a form incorporating the mass m and the stiffness parameter k = mω 2 0 : Upon the substitution of the Bloch theorem and making use of ω 0 = c/a and v x,y = cβ x,y , the equation of motion (S4) boils down to: üi,j + iω 0 β x sin(q x ) + β y sin(q y ) ui,j + 2ω 2 Assuming harmonic motion and normalizing the excitation frequency using ω 0 , the following non-dimensional dispersion relation for the square Willis monatomic lattice is obtained, which is identical to Eq. ( 15) in the main manuscript.The BZ corresponding to Eq. ( 15) is no longer a square as in reciprocal square MLs.It is, instead, of irregular deformed shape and its boundaries are defined based on setting the group velocity of the x-direction (y-direction) waves to zero to find the vertical (horizontal) boundaries, as proven in the main manuscript.Finally, it is important to re-iterate that: • The dispersion relation has a constant cutoff frequency of Ω max = 2, regardless of the magnitude of β x and β y , as shown in Figure S3(a).
• The BZ of square Willis monatomic lattices has an identical area to that of their reciprocal counterpart, regardless of the magnitude of β x and β y , which is numerically shown in Figure S3(b).The change in the BZ shape with different combinations of β x and β y is shown in Supplementary Movie S7.
• Using the traditional range of BZ for square lattices, i.e., q x ∈ [−π, π] and q y ∈ [−π, π], returns incomplete dispersion diagrams, which can be seen from the directivity plots shown in Figure S4 (See Figure 3 in the main manuscript for the complete dispersion diagrams for the same combinations of β x and β y ).
If q y = 0 (analogously q x = 0) in Eq. ( 15), the dispersion relation becomes very similar to that of the one-dimensional case presented in the main manuscript, i.e., Eq. (5).The shift at zero group velocity for a one-dimensional version of the square lattice is not equal to ϕ x,y , indicating that the BZ boundaries, indeed, do not form a square anymore.This is proven from analyzing the wave propagation in the square lattice along a single direction only, i.e., one of the components of the wavenumber has to be zero.As such, the dispersion relation reduces to (shown here for q y = 0, but identical procedure can be done for q x = 0): x ) sin 2 q x 2 = 0 (S7) Following identical procedure to that in the main manuscript to find the wavenumbers at which the cutoff frequency occurs (i.e., at a zero group velocity in the x-direction), the wavenumber solutions on BZ boundaries are: Substituting this wavenumber back into the dispersion relation in Eq. ( 15) and solving for the frequency Ω, one can show that the cutoff frequency is a function of the parameter β x and occurs at: which results in a phase shift q x s that is smaller than ϕ x , showing that the BZ of square Willis monatomic lattices should be of irregular shape.ψ = tan −1 (q y /q x ) (0,0) ψ (q x ,q y ) q y q x Figure S4: Contours of dispersion relations for square Willis monatomic lattices with a variety of combinations of β x and β y , presented via directivity plots, using the traditional definition of BZ for square lattices.As can be seen, the dispersion plots are not complete, in sharp contrast to the contours plotted using the newly proposed BZ definition shown in Figure 3 in the main manuscript.

Figure 1 :
Figure 1: Schematic of a non-reciprocal Willis monatomic lattice (WML), synthesized from discretizing Eq. (1).The lattice's masses m are evenly spaced with a lattice constant a and the displacement of the i th mass is donated as u i .As seen from the figure, two types of springs exist: (i) a conventional spring k arising from the elasticity of the rod and (ii) a negative spring −m(v 0 /a) 2 induced as a consequence of the rod's motion.Non-local coupling parameters that are proportional to the velocity of the next neighbors of an i th unit cell break the lattice's reciprocity (shown for the i th unit cell only).

Figure 2 :
Figure2: Dispersion diagram of WMLs, depicted using free-wave (top) and driven-wave (middle) methodologies.A positive (negative) β results in a positive (negative) shift ϕ, as seen from the bias towards the forward-going (backward-going) waves, while β = 0 recovers the dispersion diagram of a reciprocal ML.However, the free-wave dispersion does not capture the shifted nature of BZ, showing the driven-wave approach significance in unraveling the BZ shift ϕ.Observe that the attenuation (represented by the imaginary component of the wavenumber in the driven-wave dispersion) is smaller in WMLs relative to its reciprocal ML counterpart (β = 0) when Ω > 2, i.e., a frequency higher than the lattice's cutoff frequency Ω max = 2. Also, Ω > 2 results in a non-constant real component of the wavenumber in WML cases only.Note that, for β = ±0.5 shown here, the corresponding BZ shift is ϕ ≈ 0.3π.(bottom) Spatiotemporal FFT for the time response of WML for an impulse excitation in the left and right ends, quantifying forwardand backward-going wavenumbers, respectively, and demonstrating excellent agreement with the new BZ definition in the middle panel.

Figure 3 :
Figure 3: (a) Dispersion contours of a reciprocal square ML (β x = β y = 0) and a non-reciprocal square WML (β x = +0.5 and β y = +0.75),highlighting the definition of BZ in each case.For the WML case, the BZ area, while irregular in shape, is identical to that of its reciprocal counterpart.(b) A close-up of BZ for square WMLs, illustrating the shift in the x and y propagation directions (i.e., ϕ x and ϕ y , respectively), as well as the vanishing group velocity in the x-direction (y-direction) for the deformed vertical (horizontal) BZ boundaries.(c) Directivity plots for a variety of combinations of β x and β y , using the proposed definition of BZ for square WMLs.The depicted white lines correspond to the directions at which the cutoff frequency occurs.(d) Numerically-constructed directivity plots for β x = β y = 0 (reciprocal case) and β x = +0.5 and β y = +0.75(non-reciprocal case), achieved via spatiotemporal FFT, showing excellent agreement to their theoretical counterpart in sub-figure (c).

8 βFigure S1 :
Figure S1:(a) Group velocity c g for the forward-going waves q ∈ [0, π + ϕ], (b) Phase shift q s and (c) the Brillouin zone shift ϕ as a function of the relative velocity β.Note that the output of the group velocity c g and phase shift q s also depends on the excitation frequency Ω.The shift of the Brillouin zone ϕ signifies the values of phase shift q s at a frequency of Ω max = 2, and, thus is only a function of β.

Figure S2 :
Figure S2: Spatiotemporal FFT for the impulse response of the one-dimensional monatomic lattice, applied on the left and right ends for quantifying forward-and backward-going wavenumbers, respectively.Forward (backward) shift is observed with positive (negative) relative velocity β, which also correspond to a positive (negative) phase shift value of ϕ.A zero value of β returns the conventional dispersion relation of a reciprocal monatomic lattice with un-shifted BZ and symmetric dispersion branches.Results for non-reciprocal cases are for |β| = 0.5, as in Figure 2 in the main manuscript.

Figure S3 :
Figure S3: Numerical calculation of (a) the cutoff frequency and (b) the Brillouin zone area for the square Willis monatomic lattice over the entire range of β x and β y , showing their independence on these two parameters.