Directional Control of Rayleigh Wave Propagation in an Elastic Lattice by Gyroscopic Effects

We discuss the propagation of Rayleigh waves at the boundary of a semi-infinite elastic lattice connected to a system of gyroscopic spinners. We present the derivation of the analytical solution of the equations governing the system when the lattice is subjected to a force acting on the boundary. We show that the analytical results are in excellent agreement with the outcomes of independent finite element simulations. In addition, we investigate the influence of the load direction, frequency and gyroscopic properties of the model on the dynamic behavior of the micro-structured medium. The main result is that the response of the forced discrete system is not symmetric with respect to the point of application of the force when the effect of the gyroscopic spinners is taken into account. Accordingly, the gyroscopic lattice represents an important example of a non-reciprocal medium. Hence, it can be used in practical applications to split the energy coming from an external source into different contributions, propagating in different directions.


INTRODUCTION
According to their original definition, Rayleigh waves are a class of elastic waves that propagate on the surface of an infinite homogenous isotropic solid, and they are confined within a superficial region whose thickness is comparable with their wavelength (Rayleigh, 1885). These waves are well known by seismologists, since they are usually detected after the occurrence of an earthquake. They are also observed in other common urban activities, such as construction and demolition works and vehicular traffic, in addition to industrial processes and technologies, like mining exploration, nondestructive testing and design of electronic instruments. In the literature, they have been studied in depth especially with reference to continuous media (see, for instance, the classical treatizes by Viktorov (1967), Achenbach (1973) and Graff (1975)).
Surface waves traveling on the free boundaries of periodic media are usually referred to as Rayleigh-Bloch waves. They are of great importance in problems concerning the dynamic propagation of cracks in discrete systems, as discussed by Marder and Gross (1995) and Slepyan (2002) for uniform media and in successive works (Nieves et al., 2013) for non-uniform lattices. Similar localized phenomena may play a substantial role in non-uniform crack propagation as evidenced in (Piccolroaz et al., 2020), where lattice dissimilarity has been shown to promote or diminish localized deformations around the faces of the crack. Trapped modes associated with Rayleigh-Bloch waves in systems incorporating periodic gratings or periodic arrays of resonators were analyzed by Porter and Evans (1999), Porter and Evans (2005), Linton and McIver (2002) and Antonakakis et al. (2014) for scalar problems, by Colquitt et al. (2015) for vector systems and by Haslinger et al. (2017) and Morvaridi et al. (2018) for plates.
Recently, there has been an increasing interest in the design and fabrication of elastic media with unusual properties, often referred to as metamaterials (see the recent works by Filipov et al. (2015), Misseroni et al. (2016), Armanini et al. (2017), Bordiga et al. (2018), Bacigalupo et al. (2019), Al Ba'ba'a et al. (2020), Tallarico et al. (2020) and Wenzel et al. (2020) amongst others). In this paper, we focus the attention on an elastic metamaterial consisting of a triangular lattice connected to a system of gyroscopic spinners. The presence of gyroscopic spinners breaks the time-reversal symmetry of the system and makes the medium non-reciprocal, as proved by Nieves et al. (2020). We observe that the gyroscopic effect plays the role of magnetic bias (Bi et al., 2011) or angular momentum (Sounas et al., 2013) in linear non-reciprocal electromagnetic metamaterials and of circulating fluids in acoustic linear circulators (Fleury et al., 2014).
Firstly introduced by Brun et al. (2012) and later developed by Carta et al. (2014), a gyroscopic elastic lattice is a tunable system, whose dispersive properties can be varied by changing the spin and precession rates of the spinners. This type of metamaterial can be utilized to force waves to propagate along a line, whose direction is defined by the geometry of the medium (Carta et al., 2017). Gyroscopic spinners can also be employed to design topological insulators, where waves travel in one direction and are immune to backscattering (Nash et al., 2015;Süsstrunk and Huber, 2015;Wang et al., 2015;Garau et al., 2018;Lee et al., 2018;Mitchell et al., 2018;Carta et al., 2019;Garau et al., 2019;Carta et al., 2020). Furthermore, systems with gyroscopic spinners can be used to design coatings to hide the presence of objects in a continuous or discrete medium (Brun et al., 2012;Garau et al., 2019). In a recent work by Nieves et al. (2020), the dispersion analysis of Rayleigh waves in a semi-infinite triangular gyroscopic lattice has been carried out. The nonsymmetry of the eigenmodes of the lattice's particles at the free boundary for positive and negative values of the wave number has been linked to the non-symmetry of the system's response to an applied force, determined by means of finite element simulations. In addition, a comparison with an effective gyroscopic continuum discussed in that paper has corroborated the results for the discrete system when low values of the wave number are considered.
In this paper, we present for the first time the analytical derivation of the displacement field in a semi-infinite elastic lattice incorporating gyroscopic spinners, focusing the attention on Rayleigh waves. Conversely, the main objective of previous papers (Garau et al., 2018(Garau et al., , 2019 was related to the study of topologically-protected waveforms in lattices incorporating sub-domains with different values of the parameters of the gyroscopic spinners. The analytical results of the present paper are also verified with an independent finite element code. The displacement field produced by a point load acting on the boundary of the medium and the calculation of the energy flow demonstrate that the considered system is non-reciprocal, as the response of the forced lattice is not symmetric with respect to the point of application of the concentrated load. The influence of different physical quantities on the behavior of the gyroscopic system is also investigated through a detailed parametric analysis. While in  the wave field produced by a force on the boundary was determined numerically by using a finite element code, here the response of the system is calculated by means of a novel analytical formulation. The latter has many advantages. First, it does not require Adaptive Absorbing Layers (AAL) to prevent wave reflections at the boundaries of the computational domain; in fact, AAL are frequency-dependent and, as a consequence, their characteristic parameters need to be tuned manually every time the frequency is changed. The analytical formulation does not require the introduction of fictitious boundaries as in finite element codes. Second, the analytical formulation allows one to perform a parametric analysis quickly and efficiently, a task that is not straightforward in many finite element packages and that is performed here. Third, analytical results are necessary to validate the outcomes of finite element models. For these reasons, it is envisaged that the proposed analytical approach can be useful to the interested reader to tackle similar dynamic problems in discrete elastic systems.

MATERIALS AND METHODS
The material under consideration consists of a semi-infinite twodimensional triangular lattice of masses m, linked by elastic springs of stiffness γ, length L and whose inertia is negligible in comparison with m. The planar and three-dimensional representations of the lattice are shown in Figures 1A,C, respectively. In addition, each mass is attached to a gyroscopic spinner (see Figure 1B), whose configuration is described by the Euler angles θ, ϕ and ψ, denoted as the nutation, precession and spin angles, respectively. We assume that the nutation angle θ is small (θ ≪ 1), together with the displacement of the mass attached to the spinner. Each gyroscopic spinner has length l, moments of inertia I 1 about the axis of revolution and I 0 about the other two principal axes passing through the spinner's base. Firstly derived by Carta et al. (2018) and Nieves et al. (2018), the gyricity Ω is an independent parameter, representing the sum of the initial precession and spin rates, that remains constant throughout the motion: Ω _ ϕ + _ ψ Const. We also introduce the effective gyricity Ω * Ω[I 1 /(I 0 + ml 2 )] and the effective mass m * m[1 + I 0 /(ml 2 )], that both include the inertial contribution of the spinners.
In practice, such a lattice can be realized by constructing a triangular array of masses (represented, for example, by spheres) connected by thin elastic rods. At the bottom part of each mass a cylindrical hole can be drilled, where the tip of the gyroscopic spinner can be inserted. The connection needs to be frictionless, so that the spinning motion of the gyroscope is not transmitted to the mass, which can only move in the x 1 -or x 2 -direction without rotating. The spinning motion of the gyroscopic spinner can be applied by using an electric motor; in this way, its spin rate can remain constant during the motion. Each gyroscopic spinner is pinned at the base and its axis is parallel to the x 3 -direction in the initial configuration.
The boundary of the semi-infinite lattice is subjected to an oscillatory force of amplitude P and prescribed radian frequency ω 0 (see Figure 1A), which generates Rayleigh waves localized at the boundary and bulk waves propagating inside the medium. The displacement of a generic particle, whose position is identified by the multi-index n (n 1 , n 2 ) T , is denoted by u (n 1 ,n 2 ) (t) and is a function of time t.
We introduce the following normalizations: where the tilde symbol indicates a dimensionless quantity. In the following, the tilde is omitted for ease of notation and it is assumed that all the appearing quantities are dimensionless.

Governing Equations of the Forced Lattice in the Transient Regime
The equations of motion of a mass within the bulk of the lattice are given in normalized form as (Garau et al., 2019) € u where with Further, the dot denotes the time derivative, the vectors a (i) (i 1, 2, 3) are as shown in Figure 1A, and the matrix R is The governing equations for a mass belonging to the boundary of the semi-infinite lattice (n 2 0) are where is the Kronecker delta.
In addition, we assume that the lattice is at rest initially, namely:

Alternative Representations of the Governing Equations in the Transient Regime
We introduce the integer p 2 n 1 and the transient solution in the form where the time-dependent complex displacement amplitude U (p,n2) (t) U(p + n 2 , n 2 , t) satisfies the conditions of zero initial velocity and displacement. In addition, note that for α, β ∈ Z U ( p,n2 ) (t) n 1 n 1 + α n 2 n 2 + β U ( p+2α+β,n2+β ) (t).
In terms of these complex displacement amplitudes, Eqs. 2, 7 become: and where I is the identity matrix. In the last equation, the complex representation for the applied force was used.

Laplace and Fourier Transformed Equations for the Complex Displacement Amplitudes
Next, we apply the Laplace transform in time t to Eqs. 12,13 and we use the fact that the displacement amplitudes satisfy zero initial conditions. After this, we apply the discrete Fourier transform with respect to n 1 ∈ Z x 1 . In what follows, U LF n 2 denotes the Laplace and discrete Fourier transform of the solution, defined as Here, s and k 1 are the Laplace and discrete Fourier transform parameters, respectively. The Fourier transform parameter k 1 will also be referred to as the (normalized) wave number.
For the equations in the bulk, we obtain that the complex displacement amplitude satisfies for n 2 > 0, where the factor e (−n2/2) appearing throughout due to the application of the discrete Fourier transform has been cancelled (see Eq. 10). Similar conversions for Eq. 7 describing the forced boundary at n 2 0 lead to

Analysis of the Forced Problem in the Steady-State Regime
The transition to the steady-state regime (i.e. when t → + ∞) is made by multiplying Eqs. 15, 16 by the Laplace transform parameter s and taking the limit as s → + 0. In this limit, we define the discrete Fourier transform of the displacement amplitude in the steady-state regime as From Eqs. 15, 16, these transformed amplitudes then satisfy for n 2 > 0 and Frontiers in Materials | www.frontiersin.org January 2021 | Volume 7 | Article 602960 for n 2 0. Here, ω 0 − i0 lim s → +0 (ω 0 − is), and in what follows "i0" will represent a small imaginary term with positive imaginary part. The presence of this term enables one to determine how singular points of the Fourier transformed solution approach the real axis in the complex plane in passing to the steady-state regime. In this limit, the location of these singular points in the complex plane provides a way of ascertaining the direction of waves propagating in the medium relative to the load position. This is linked to the causality principle (see Section 3.3.2 in (Slepyan, 2002) for more details).

Solution to the Transformed Steady-State Forced Problem
The transformed amplitudes U F n2 are sought in the form Here, Λ is such that |Λ| < 1 for n 2 ≥ 1, since we are considering waves that decay into the bulk of the lattice.
The solution in Eq. 20 is inserted into Eq. 18 in order to find the eigensolutions for waves in the bulk; using Eq. 19 it is then possible to satisfy conditions at n 2 0 and find the response to the applied load. The full derivation is reported in the Supplementary Material, where it is shown that the solution has the form where U 0 , Λ and M 2 are the matrices specified in Supplementary Equations S8, S7, S11, respectively. Note, the right-hand side of Eq. 21 defines a 2π -periodic (4π -periodic) function for real k 1 when n 2 is even (odd).
If the gyricity is zero and the load acts in the horizontal (vertical) direction, the solution in Eq. 21 describes a vector function whose first component has even (odd) real and imaginary parts as functions of k 1 , whereas the second component has odd (even) real and imaginary parts. When gyricity is non-zero, these components have neither even nor odd real and imaginary parts. Such properties can be observed in Figure 2, where we report the horizontal and vertical components of the Fourier transform of the steady-state solution U F 0 , when Ω * 0 (parts (A) and (B)) and Ω * 1 (parts (C) and (D)) for a load acting in the horizontal direction. By further inspection of the analytical solution in Eq. 21, we pose the attention on the numerator in h j (j 1, 2), belonging to U 0 and M 2 (see Supplementary Equations S9, S11), where the presence of gyricity induces a competition between terms and the amplitude of the components depends on the sign of the wave number. This feature of the solution in Eq. 21 for Ω * ≠ 0 is responsible for the symmetry breaking in the considered problem (see also ).

The Displacements in the Forced Lattice and Associated Wave Phenomena
The displacements in the lattice can be found from inverting the discrete Fourier transform, taking into account the periodicity properties of Eq. 21. Then, the displacements of the nodes in the lattice are given by where Here, we note that U (n 1 ,n 2 ) lim t → +∞ U (p,n2) (t), where p 2 n 1 and the complex time-dependent amplitude U (p,n2) (t) were introduced in Section 2.2.

The Rayleigh Waves
Rayleigh waves carry some of the energy produced by the point load in both directions along the boundary. The lack of symmetry in the integral kernel of Eq. 23 with respect to the zero wave number results in a disparity between amplitudes of waves outgoing from the source to the left and to the right of the load.
The Rayleigh waves are defined by the degenerate values of the wave number of M 2 for a given frequency ω 0 , and these singular points are represented by the simple poles of Eq. 21. One can check that there exist no other singular points of this function for ω 0 ≠ Ω * (the special case ω 0 Ω * will be discussed later). On the dispersion diagram, these simple poles correspond to the intersections of the line ω ω 0 with the curve ω ω R , where ω R is given by and represents the dispersion curve of the system associated with Rayleigh waves. 1 The dispersion curve is plotted in Figure 3.
In Figure 3 we limit our attention to the intervals of periodicity for Eq. 21 mentioned in Section 2.4. Since ω 0 − i0 is present in Eq. 21, as discussed in Section 2.3, the degenerate wave numbers of Eq. 21 are perturbed and are slightly shifted from the real k 1 -axis in the complex plane. The causality principle enables the new locations of these points to be determined using information concerning the group velocity v g dω R /dk 1 (see also (Slepyan, 2002)). In particular, if at the previously mentioned intersection points v g > 0 (v g < 0) the frequency ω 0 and the corresponding wave number k 1 define waves propagating to the right (left) of the load. Additionally, when the frequency is ω 0 − i0, the associated degenerate value of Eq. 21 is located at k 1 − i0 (k 1 + i0). Consequently, the dispersion diagram in Figure 3 predicts that the degenerate values required later for the computation of outgoing waves from the source are at Using the above information, we can calculate the form of the Rayleigh waves produced by the point load by determining the simple poles of Eq. 23 and employing the residue theorem to compute Eq. 23 at a point located far from the application point of the oscillatory force. With this approach, it is possible to show that for n 2 0 where the term on the right-hand side is the Rayleigh wave produced to the right of the point load. On the other hand, we have where the right-hand side defines the Rayleigh wave traveling away from the load to the left. Here, It can be verified that when Ω * ≠ 0 Hence, the resulting dynamic response at the boundary of the semi-infinite gyro-elastic lattice to the far left and right of the point load is different. Moreover, these outwardly-propagating boundary waves cause the nodes to follow elliptical trajectories, as shown by the eigenmode analysis developed by Nieves et al. (2020).

Bulk Wave Radiation
Waves are also radiated along the rows of the lattice in the bulk. As n 2 → ∞, the amplitudes of such waves decrease. These waves are again attributed to the singular points of Eq. 21 associated with k 1 ± k 1R ∓ i0, ± (2π − k 1R ) ± i0, which also define Rayleigh waves along the boundary. The form of the waves depends on the row of the lattice considered and it is obtained following the same procedure outlined in Section 2.5.1. For n 2 > 0 with n 2 even, the displacements to the far right of the lattice behave as where the term on the right-hand side describes a wave traveling along the row defined by n 2 > 0 to the right. On the other hand, we have where the term on the right-hand side represents a wave propagating to the left along the row n 2 > 0. We note that for the wave numbers k 1 ± k 1R , ± (2π − k 1R ) and the frequency ω 0 − i0, the functions Λ j , j 1, 2, are complex with modulus less than unity. Hence, the matrix Λ n 2 in Eq. 21 ensures that these waves have a decreasing amplitude for increasing n 2 . For rows in the lattice defined by odd n 2 , we recall that Eq. 21 is a 4π -periodic function in k 1 . Hence, the inversion formula for the discrete Fourier transform is as presented in the second equation in Eq. 23. In this case, a slightly modified procedure is required to compute the waves radiated in the bulk, taking into account the singular points of the function Eq. 21 on the interval [−2π, 2π]. For odd n 2 , the displacements for n 1 → ∞ are whereas for n 1 → −∞ we have We point out that the amplitudes of the waves along the odd rows of the lattice, specified in Eqs. 32, 33, take into account the contributions from the wave numbers k 1 ± (2π − k 1R ) (see Figure 3).
There also exist preferential directions for energy radiation in the bulk. When the gyricity is zero, along these specific lines in the lattice, the nodal displacements decay slowly as O(n −1/3 2 ) for n 2 → ∞, as shown by Slepyan (2010). These directions can be identified by determining when Λ j , involved in Eq. 21, is complex with Λ j 1, j 1, 2, for the frequency equal to ω 0 . As discussed by Slepyan (2010), this effect is purely attributed to the lattice's micro-structure and is not found in the analogous continuous model. When the gyricity Ω * ≠ 0 is introduced, preferential directions for wave propagation can remain and increasing the gyricity causes the associated displacement amplitudes to decrease.

Resonant Modes
Next, we discuss the resonant case when the steady-state solution to the considered problem does not exist. We show this by investigating the derived solution in Eq. 21 along the boundary, where n 2 0.
When k 1 → 0, for either horizontal or vertical loading at the lattice boundary, the solution in Eq. 21 is bounded for ω 0 ≠ Ω * and admits the following asymptotic representation: where FIGURE 3 | Dispersion diagram for k 1 ∈ [−2π, 2π], relevant for Eq. 23 when n 2 is odd. The black (gray) curve indicates the function ω R for Ω * 1 (Ω * 0). The line ω ω 0 is shown for ω 0 0.6. The intersections of this line with the curve of ω R for Ω * 1 (Ω * 0) are represented by black (gray) crosses. The interval [−π, π] is instead relevant for Eq. 23 when n 2 is even. Wave numbers k 1 k 1R and 2π − k 1R , defining singular points of Eq. 21, are also shown.
Frontiers in Materials | www.frontiersin.org January 2021 | Volume 7 | Article 602960 and C 0 denotes a constant vector with non-zero entries depending on the gyricity Ω * , the frequency ω 0 and the load amplitude P. Clearly, D(Ω * , ω 0 ) 0 when ω 0 Ω * . In this degenerate case, the asymptote near k 1 → 0 takes the form that indicates the appearance of other singular points near k 1 0 in the solution in Eq. 21. In the above, if the load is horizontal (vertical) the vector S 0 (−i, −1) T (S 0 (1, −i) T ). We note that the singular points in this asymptote approach each other as the Laplace transform parameter s → + 0, with their limiting location being the real axis at k 1 0. An illustration showing the typical behavior of the first component of Eq. 21 in this case is given in Figure 4A. There, the kernel has a simple pole for positive wave numbers but not for negative wave numbers, which implies that the structure does not support Rayleigh waves propagating to the left of the load on the boundary. The second component of Eq. 21 possesses the same singular features. Further to Eq. 36, in this situation U F 0 can be written in the form: where W F 0 represents the term that is singular at k 1 k 1R − i0 of Eq. 21 and B F 0 is a bounded function for k 1 ∈ [−π, π]. Upon applying the inverse of the discrete Fourier transform, while the last two terms in the right-hand side of Eq. 37 give bounded contributions to the displacement, the first term leads to a function that is singular for x 1 n 1 + n 2 /2 0 and s → + 0. Hence, the displacement at the location of the load is unbounded for t → + ∞.
The asymptotic representation in Eq. 36 is also in agreement with the dispersion diagram for ω 0 Ω * , shown in Figure 4B.
The associated horizontal line on this diagram intersects the curves along which Λ j 1 (j 1, 2), represented by gray lines. These additional curves are useful in characterizing the main contribution of the integral contained in Eq. 23 in describing the lattice far-field and are connected with the appearance of waveforms and preferential directions in the bulk lattice (see also Slepyan (2010)). At k 1 0 on the dispersion diagram, the group velocity at the associated intersection point is zero (see Figure 4B, where at the point (k 1 , ω 0 ) (0, Ω * ) the group velocity is zero). Physically, such a point represents a resonance mode of the considered system, as discussed in (Slepyan, 2002, Section 3.3.5). The energy associated with this mode is unable to leave the location where the force is applied. Thus, the energy density at this location becomes unbounded as t → + ∞. In other words, in the case ω 0 Ω * the steady-state solution does not exist.

Determination of Energy Flow in the Steady-State Regime
In this section, we analyze the energy carried by the Rayleigh waves and by the waves radiated into the bulk of the gyroscopic lattice by calculating the energy flow through the boundary of a sufficiently large region of the medium, as shown in Figure 5. The considered region is the rectangle S enclosed by the half-plane boundary and by the segments zS i (i 1, . . . , 4), indicated in Figure 5.
In the steady-state regime, the rate at which energy is introduced into the system through the action of an oscillating point load applied at the node (n 1 , n 2 ) (0, 0), having the vector amplitude P and frequency ω 0 , can be computed using the formula: where the overline denotes the complex conjugate and U (n1,n2) is defined in Eq. 23. The rate at which energy flows through each segment zS i is given by the classical formula introduced by Brillouin (see Chapter V in Brillouin (1953)): for 1 ≤ i ≤ 4. In Eq. 39, F (n j ) represents the elastic force supplied by the nodes outside the considered region and connected to the node indicated by n j , positioned in the immediate proximity of zS i , within S. The set of the nodes with multi-indices n j located in the vicinity of zS i , inside S, is denoted by D i . For the conservation of energy law, the following equality is satisfied: In Section 3, the above formulae for the energy flow will be used to show quantitatively that the presence of gyroscopic spinners breaks the symmetry of energy propagation with respect to a vertical line passing through the application point of the load.

Finite Element Model
The results of the analytical formulation illustrated in the previous sections will be verified with a finite element model built in the commercial software Comsol Multiphysics (version 5.4).
The numerical model consists of massless truss elements and point masses inserted at the lattice's nodes. The effect of gyricity is simulated by imposing at each node a force that is proportional to the velocity, as in the first term on the right-hand side of the governing Eq. 2. Here the computational domain is finite, while the analytical treatment developed above is for a semi-infinite medium. In particular, the lattice has dimensions 100 × 46 3 √ . Adaptive Absorbing Layers (AAL) are introduced close to the top and vertical boundaries, in order to prevent waves from being reflected at these boundaries. AAL are created by assigning to the links located inside those regions a complex elastic modulus. A point load is applied to the middle point of the bottom boundary of the domain. The numerical simulations are performed in the time-harmonic regime.

The Displacement Field in the Semi-Infinite Gyro-Elastic Lattice Loaded on the Boundary
In this section, we show the response of the semi-infinite lattice with embedded gyroscopic spinners to a point force applied on the boundary. In particular, we investigate how the gyricity affects the behavior of Rayleigh waves traveling along the boundary and of the waves propagating into the bulk of the medium. Figure 6 shows the total displacement amplitude calculated at each node of the discrete system, produced by an oscillating force applied on the boundary. The results are based on the solution in Eq. 22, derived in Section 2. The force has unit amplitude and frequency ω 0 0.6, and it acts in the horizontal (vertical) direction in parts (A) and (C) (parts (B) and (D)).
In parts (A) and (B) of Figure 6 the effective gyricity is set equal to zero. It is apparent that when Ω * 0 the response of the system is symmetric with respect to the point of application of the force. This in agreement with the classical theory on Rayleigh wave propagation in elastic lattices without gyroscopic effects. We also observe that in the vicinity of the load there exist evanescent modes, which induce localized deformations in the neighboring links. In the bulk there are preferential directions along which waves are radiated from the source, as described in Section 2.5.2. The region between these rays show outwardly propagating oscillations that decay like O(n −1/2 2 ). Comparing parts (A) and (B), we note that in (B) there exist waves with larger amplitudes. Moreover, by looking at the wavelengths, we observe that the bulk waves generated by a horizontal force are of the shear type, while the vertical force induces bulk waves of the pressure type, as expected.
Figures 6C,D illustrate the total displacement amplitude field produced by a horizontal and vertical force, respectively, when the effective gyricity is Ω * 1. As in Figures 6A,B, the force has unit amplitude and frequency ω 0 0.6. The main difference between Figures 6A-D is that when the effective gyricity is non-zero the response of the system ceases to be symmetric with respect to the vertical line passing through the point of application of the force. In part (C) we note that the amplitude of Rayleigh waves propagating to the right of the force is much larger compared to that of the surface waves traveling to the left. In addition, the force appears to activate waves in its vicinity having significant amplitude and propagating at 120°to the positive horizontal axis. On the other hand, the displacement field produced by a vertical force (see part (D)) shows that the symmetry is still broken, but the amplitudes of the waves traveling to the right and to the left of the force along the boundary are now comparable.
We also point out that the non-symmetric displacement field in Figure 6 can be inverted by changing the sign of the effective gyricity Ω * .
In Figure 7 we present the total displacement amplitude fields computed by using the finite element model developed in Comsol Multiphysics and described in Section 2.7. The values of the parameters are the same as those considered in Figure 6. Comparing Figures 6, 7, we observe that the numerical and analytical results show an excellent agreement. This confirms the validity and accuracy of the analytical treatment discussed in Section 2.

The Energy Flow
Here, we report the analytical results concerning the energy flow furnished by the external force, referred to as W in , and the energy flow passing through each segment zS i in Figure 5, denoted as W out i (i 1, . . . , 4). The formulae for W in and W out i are given in Section 2.6.
As in Section 3.1, two different values of the effective gyricity are taken, namely Ω * 0 and Ω * 1. Moreover, two directions of the point load are considered, i.e. horizontal and vertical. The values of W in and W out i for all the examined cases are summarized in Table 1. By looking at Table 1, we notice that when the effective gyricity is zero W out 1 + W out 2 W out 3 + W out 4 , implying that the response of the system is symmetric with respect to the application point of the force. The above equality ceases to hold when Ω * 1. We also note that the energy partition between vertical and horizontal boundaries depends on the gyricity and on the direction of the force. Furthermore, for any of the four cases considered in Table 1, the energy balance is satisfied, since W in 4 i 1 W out i .

Parametric Analysis
In order to assess how the response of the micro-structured medium is affected by different physical quantities, we perform a parametric analysis where we vary the direction of the force, the radian frequency ω 0 of the force and the effective gyricity Ω * of the spinners. The behavior of the medium is evaluated quantitatively by calculating the percentages of the energy flows W out i passing through the boundaries zS i (i 1, . . . , 4) (see Figure 5) with respect to the energy input W in due to the external source.
In Figure 8 we show how the energy introduced into the system by the external oscillating force is divided into two parts, propagating in opposite directions relative to the position of the point force. In particular, the circles (squares) indicate the percentages of the energy flowing to the right (left) of the force with respect to the input energy flow. In Figure 8 it is assumed that the force acts in the horizontal direction. The five diagrams correspond to Ω * 0, 0.25, 0.5, 0.75, 1; in each diagram, several values of the frequency ω 0 of the point force are considered.
In Figure 9 the values of the effective gyricity and of the frequency of the external source are identical to those considered in Figure 8, but the outcomes are obtained by applying a concentrated oscillating load acting in the vertical direction. In both Figures 8, 9, when Ω * 0 the incoming energy is split into two equal contributions that propagate to the left and to the right of the force, both along the boundary and inside the bulk. Conversely, when Ω * ≠ 0 we FIGURE 7 | Same as in Figure 6, but obtained from numerical simulations performed in Comsol Multiphysics. observe a non-symmetrical distribution of energy in the system. We also notice that even when Ω * ≠ 0 the symmetry in the energy flow partition is retrieved for some specific values of the frequency, that change with the value of the effective gyricity. The crosses in both Figures 8, 9 represent the sums of the energy flow components traveling inside the system with respect to the input energy; it can be seen that, for every case considered, there is balance between input and output energy flows.
The insets in Figure 8 present the color maps of the displacement fields, computed at given gyricities and frequencies of the external force. In part (A), where Ω * 0, the wave pattern is clearly symmetric. When gyricity is introduced, as in part (B), the symmetry of the displacement field is broken, even if the frequency of the external force remains the same. In part (C), the wave pattern for a frequency near the intersection between the two sets of data is shown; in this case, the energy is split into two approximately equal parts traveling to the left and to the right of the load, but on each side the energy amounts propagating into the bulk and on the boundary are different. Here, it is apparent that the energy flowing into the bulk in the left-hand part of the lattice is approximately equal to the amount of energy carried in the right-hand part, where the energy is mainly concentrated on the boundary. FIGURE 8 | Percentages of energy flows propagating to the right (W out 1 + W out 2 ) and to the left (W out 3 + W out 4 ) of the point force with respect to the energy flow produced by the external load (W in ), calculated for several values of effective gyricity Ω * and frequency ω 0 . The results in this figure correspond to a point force acting in the horizontal direction. The insets illustrate the wave fields at specific frequencies ω 0 and gyricities Ω * .
In Figures 10, 11 we focus the attention on Rayleigh waves, showing how the contributions of the energy flows corresponding to surface waves traveling to the right (W out 1 ) and to the left (W out 4 ) of the point load vary with the effective gyricity of the spinners and with the frequency of the external force. In Figure 10 (Figure 11) the force acts in the horizontal (vertical) direction. It is apparent that the direction of the force, the gyricity and the frequency all influence the response of the system and, in particular, the relative amount of total and surface energy propagating to the left and to the right of the applied force.
The insets in Figure 10 illustrate, for different values of gyricity and frequency of the excitation, the amplitudes of the displacement components u (n 1 ,0) 1 and u (n 1 ,0) 2 of the lattice's nodes at the boundary and the displacement magnitude . The symmetry of the displacement field in the lattice without gyroscopic spinners, shown in Figure 10A, is broken when gyricity is incorporated into the system, as illustrated by the non-symmetric displacement profiles of Figures 10B,E. Now, we consider a special case, where the surface waves propagating to the left of the point load exhibit a negligibly small amplitude. This can be obtained by taking the effective gyricity Ω * 0.5, the frequency of the external force ω 0 0.35 and applying a horizontally-acting force (see also Figures 8, 10). The displacement field for this choice of the parameters, presented in Figure 12, clearly shows that the energy coming from the external source propagates into the bulk and practically only to the right of the FIGURE 9 | Same as in Figure 8, but for a vertically-acting point force.
Frontiers in Materials | www.frontiersin.org January 2021 | Volume 7 | Article 602960 force along the boundary. A similar effect could also be achieved in the resonant case ω 0 Ω * , but there the steady-state solution does not exist (see Section 2.5.3). The energy flow partition for the system in the neighborhood of ω 0 Ω * is investigated in Figure 8C, which shows that the energy flux distributions vary rapidly in the vicinity of ω 0 Ω * . Another interesting case is represented by the scenario where almost all the energy propagates along the boundary. This situation is shown in Figure 13, where Ω * 1, ω 0 0.9 and the point force acts in the vertical direction (see also Figures 9, 11). From the displacement field in this figure, it is apparent that most of the energy is localized at the boundary of the medium. Moreover, in this case the discrete system shows negligible preferential directionality (see also Figure 9E).

DISCUSSION
The diagrams of the displacement amplitude fields in Figures 6,  7, as well as the values of the energy flows in Table 1, show that the gyricity is capable of breaking the symmetry in the energy propagation of both Rayleigh and bulk waves propagating from the external source. This is a consequence of the non-reciprocity of the gyro-elastic lattice. 2 FIGURE 10 | Energy flow percentages associated with Rayleigh waves traveling to the right and to the left of the point force, indicated by W out 1 and W out 4 respectively, for different values of effective gyricity and frequency of the external source. In these computations, the force acts in the horizontal direction. The insets show the displacements of the points on the lattice boundary, calculated at the frequencies indicated by the arrows.
Examining the outcomes of the parametric analysis presented in Section 3.3, in particular Figures 8, 9, it is apparent that the distribution of energy in the system strongly depends on both the effective gyricity Ω * of the spinners and the radian frequency ω 0 of the external oscillating force. Generally, at low frequencies it is easier to partition the energy flow into two significantly different contributions, propagating in opposite directions. Nonetheless, for any given value of the effective gyricity, it is possible to find one or more values of ω 0 at which the energy is split into two equal parts. The lowest value of this frequency increases as the effective gyricity is increased. Comparing  Figures 8, 9, we observe that the energy symmetry breaking is more evident when the force acts in the horizontal direction. Of course, for an inclined direction of the force, the response of the system can be determined from the principle of superposition by summing the effects due to the horizontally-and vertically-acting forces, since the considered problem is linear.
The investigation of the individual components of the energy flows associated with Rayleigh waves, presented in Figures 10,  11 for a horizontally-and vertically-acting point force respectively, reveals that the propagation of energy along the boundary of the medium is affected significantly by the direction of the force. In the case when Ω * 0, the energy flow along the boundary is the same in both directions and it decreases (increases) as the frequency of the external source acting in the horizontal (vertical) direction is increased. On the other hand, when Ω * ≠ 0, the energy flowing to the right (left) is larger than in the opposite direction if the direction of the force is horizontal (vertical). Independently of the direction of the FIGURE 11 | Same as in Figure 10, but for a vertically-acting point load.
Frontiers in Materials | www.frontiersin.org January 2021 | Volume 7 | Article 602960 force, the difference between the energy flows in the two opposite directions generally decreases as the frequency is increased.
The parametric analysis discussed in Section 3.3 has been helpful in identifying special cases, where Rayleigh waves propagate only in one direction (see Figure 12) or where bulk waves have very small amplitudes compared with those of surface waves (see Figure 13).
The analytical formulation also made it possible to analyze resonant regimes where the steady-state solution cannot be reached (see Section 2.5.3).
The capability of the considered micro-structured system in creating preferential directionality in the propagation of waves can be exploited to design novel energy splitters, where the desired amount of energy propagating in a specific prescribed direction can be varied by changing the effective gyricity of the spinners. The tunability of the proposed model is an essential tool, that can be used in many practical applications where it is required to vary the output depending on the contingent needs. Important examples include electronic instruments converting mechanical energy into electric energy (and vice versa) and elastic filters that can be utilized both for protection and energy harvesting.