{'fraction_non_alphanumeric': 0.03973050496964846, 'fraction_numerical': 0.01383496764725502, 'mean_word_length': 4.497726272553909, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Incremental constitutive tensors and strain localization for prestressed elastic lattices: Part II -- incremental dynamics

Floquet-Bloch wave asymptotics is used to homogenize the in-plane mechanical response of a lattice grid of elastic Rayleigh rods, possessing a distributed mass density together with a rotational inertia and subject to incremental dynamics, superimposed to a given state of axial loading of arbitrary magnitude. In contrast to the quasi-static energy match (addressed in Part I of the present study), the vibrational properties of the lattice are directly represented by the acoustic tensor of the equivalent solid, without passing through the constitutive tensor which is instead arrived via energy equivalency. The acoustic tensor is shown to result independent of the rods' rotational inertia and allows directly the analysis of strong ellipticity of the equivalent continuum, evidencing coincidence with macro-bifurcation in the lattice. On the other hand, the dynamics of micro-bifurcation corresponds to a vanishing frequency, at finite wavelength, of the lowest dispersion branch of the lattice. Dynamic homogenization reveals the structure of the acoustic branches close to ellipticity loss, and of forced vibrations (both in physical space and Fourier space) showing low-frequency wave localizations. A comparison between the strain localization occurring near ellipticity loss and the forced vibration of the lattice, both corresponding to the application of a concentrated pulsating force, shows that the homogenization technique allows an almost perfect representation of the lattice. Therefore, the presented results pave the way for the design of architected cellular materials to be used in applications where extreme deformation are involved.


Introduction
A lattice grid of axially-preloaded elastic rods, subject to in-plane incremental kinematics, has been shown (in Part I of this study, via quasi-static energy match with an equivalent continuum), to provide a way to the realization of a material evidencing macroscopic (strain localization) and microscopic (lattice buckling) instabilities within the elastic range. Such lattice represents an example of architected cellular material to be used for large strain applications where instabilities become a major issue. The latters are investigated in the present article within the dynamic range, thus employing an asymptotic homogenization based on the Floquet-Bloch wave technique, exploited to a level of generality never achieved so far. The vibrational properties of a cellular material are deeply affected by the emergence and development of localized signals, edge waves, and topologically protected modes [1][2][3][4][5][6], an example being that reported in Fig. 1, where the dynamic emergence and propagation of discontinuity wavefronts (rectilinear and curvilinear) is shown in the so-called 'pinscreen', a material (made up of a perforated plate having each hole filled with a movable pin) on the verge of ellipticity loss.
Homogenization techniques based on the asymptotic analysis of wave solutions dates back to Brillouin [7] and Born [8], and has received significant contributions in recent years when the case of random and periodic media has been considered [9][10][11][12][13][14][15][16][17] and also extended to the high-frequency regime [18][19][20]. With the exception of [17], these developments have been so far produced for the analysis of wave propagation in continuous materials, not in structures, so that their practical implementation required the systematic use of numerical techniques (typically finite elements). It is shown in the present article that low-frequency effective properties can be derived analytically for lattices composed of rods (incrementally loaded in-plane) through a direct computation of the wave asymptotics 1 . Recent results on beam lattices [17] are here extended to the case of elastic lattices, axially stressed up to an arbitrary amount, whose incremental dynamics is derived without restrictions on the rods' constitutive law and without neglecting the rotational inertia of the rods' cross section. The low-frequency asymptotics of the lattice waves is shown to be governed by the spectral properties of the acoustic tensor associated to the effective solid and capturing the incremental response of the system.
While the final results coincide with those obtained via energy match, it is shown that Floquet-Bloch asymptotics leads directly to the acoustic tensor, not to the constitutive operator, of an equivalent, prestressed, solid. In this way, the correspondence between loss of ellipticity in the continuum and degeneracy of the acoustic properties in the lattice is directly demonstrated. Moreover, the elastic rods are characterized by an axial mass density, equipped with rotational inertia (the Rayleigh model [21][22][23]), so that it is possible to prove that the latter does not influence the vibrational properties, expressed by the acoustic tensor, of the equivalent continuum. The homogenization, obtained analytically, is exploited to investigate the lattice response near macro-instability threshold (coincident with the failure of strong ellipticity of the effective continuum), therefore unveiling microscopic features of lattice dynamics loaded up to the verge of shear band formation.
The response to the application of a pulsating concentrated force is finally analyzed in the spirit of [24] and [25], when applied both to the lattice (in the physical and Fourier spaces) and to the equivalent solid, with special detail on low-frequency wave localizations. The time-harmonic infinite-body Green's function (see [25]) of the effective solid is compared to the forced response of lattice obtained through the finite elements methodology, in which the application of PML boundary conditions to a finite grid simulates an infinite domain. The acoustic branches of the lattice structure and its continuum approximation are shown to be in almost perfect agreement, thus confirming that the lattice is an excellent candidate for the realization of architected cellular materials to be used to harness micro-buckling and strain localization, two features which can control the dynamic response of a structure with tunable functionality [26,27]. This paper is organized as follows. The mathematical setting for incremental wave propagation is developed in Section 2 for an axially pre-loaded lattice of elastic rods, organized in an arbitrary periodic geometry. The asymptotic analysis of lattices waves is derived in Section 3, leading to the homogenization result that provides the acoustic tensor associated to the incremental effective Cauchy continuum. The developed homogenization scheme is applied in Section 4 to the same grid of elastic rods considered in Part I, so that failure of ellipticity and micro-buckling are investigated now through the analysis of the dispersion characteristics of the lattice. In Section 5, the emergence and characterization of dynamic strain localization is presented by means of a perturbative approach (applied to the lattice and the effective continuum), which demonstrates how the macroscopic localization direction, as well as the number of localization bands, can effectively be designed through appropriate modifications of the underling lattice structure.

Incremental dynamics of preloaded lattices: governing equations
In this section, the governing equations for incremental wave propagation in an axially-preloaded lattice of elastic rods are formulated. These are obtained (i) by solving for time-harmonic vibrations the incremental dynamics of a single rod (derived in Appendix A), (ii) by using this solution to formulate the equations of motions of a single unit cell, and finally (iii) by applying the Bloch theorem to obtain the equations governing the incremental dynamics of the infinite lattice.
An infinite two-dimensional lattice structure is considered, composed of nonlinear elastic rods which are axially preloaded (or prestretched) from an unloaded reference configuration B 0 to a preloaded configuration B used as reference in an updated Lagrangean formulation of incremental dynamics (Fig. 2).

Fig. 2.
A periodic two-dimensional lattice of (axially and flexurally deformable) elastic rods is considered, preloaded from the stress-free configuration B 0 (upper part, on the left) by means of a pure axial loading state, transforming B 0 to the current, preloaded configuration B (upper part, on the right). The latter configuration, used as reference in an updated Lagrangean description, can be represented as the tessellation of a single unit cell along the vectors of the direct basis {a 1 , a 2 }. The incremental dynamic response (lower part where an incremental deformation is shown at three different instants of time) is defined on B by the incremental displacement field of each rod u(s, t), here decomposed in an axial and transverse component, u(s, t) and v(s, t).
By introducing a local coordinate s k for each rod of a given unit cell, the incremental kinematics is described by means of the following fields u k (s k , t) = {u k (s k , t), v k (s k , t)} , ∀k ∈ {1, ..., where N b is the number of rods in the unit cell, and the two in-plane displacement components, respectively axial and transverse, are denoted by u k (s k , t) and v k (s k , t), while the rotation of the cross-section θ(s, t) is assumed to satisfy the unshearability condition θ(s, t) = v (s, t). 2 Time-harmonic solutions are sought, so that, by introducing the circular frequency ω, the dependence on time t can be represented as ,v k (s k )} are functions of the coordinate s k only.

Bloch's theorem
Wave propagation in an infinite periodic elastic medium can effectively be analyzed through the application of Bloch theorem. Essentially, the theorem states that the time-harmonic solutions of the equations of motion possess a modulation in space having the same periodicity of the medium, a condition expressed by the following requirement 3 u( where k is the Bloch vector and the modulation ϕ(x) is periodic with respect to the direct basis, so that it satisfies ϕ(x + n j a j ) = ϕ(x) ∀{n 1 , n 2 } ∈ Z 2 ∀x ∈ R 2 .
Eq. (2a) can be equivalently expressed as Note that, in the case of a lattice made up of rods, the waveform ϕ(x) (as well as the field u(x, t)) is described by the displacement of the rods constituting the unit cell, and thus it is defined only for x corresponding to the location of the structural elements.
The importance of the Bloch theorem lies in the fact that it allows to obtain the response of an infinite structure by complementing the equations of motion of a single unit cell with suitable boundary conditions. These so-called 'Floquet-Bloch' conditions enforce the periodicity of the lattice, by relating nodal displacements and forces on the boundary of the unit cell, according to the property expressed by Eq. (2b).

Time-harmonic solution of a preloaded elastic rod
The analytic solution is obtained for time-harmonic vibrations of a Rayleigh rod. In the framework of a linearized theory, the equations of motion governing the incremental dynamics of an axially pre-stretched Rayleigh rod are the following where a superimposed dot denotes time differentiation, γ(λ 0 ) is the current linear mass density, γ r (λ 0 ) is the current rotational inertia, λ 0 is the axial pre-stretch and P (λ 0 ) the corresponding axial preload (assumed positive in tension), while A(λ 0 ) and B(λ 0 ) are, respectively, the current axial and bending stiffness. The analytic derivation of equations (3) is reported in Appendix A, where the formulation provided in Part I of this work is extended to the dynamic case. Moreover, the derivation of the current stiffnesses A(λ 0 ) and B(λ 0 ) from strain-energy functions, as well as the their identification for rods made up of an incompressible elastic material can be found in Appendix A of Part I. In the following, the parameters γ(λ 0 ), γ r (λ 0 ), A(λ 0 ), and B(λ 0 ) will simply be denoted as, γ, γ r , A, and B, and treated as independent quantities for generality. The substitution of Eq. (1) into Eq. (3) leads to a system of linear ODEs for the functionsû(s) andv(s). As the system is fully decoupled, the solution is easily obtained in the form where {C u 1 , C u 2 , C v 1 , ..., C v 4 } are 6 arbitrary complex constants and the characteristic roots β u j and β v j are given by with l being the current length of the rod,ω = ω l γ/A the non-dimensional angular frequency, p = P l 2 /B the non-dimensional preload, Λ = l/ B/A the slenderness of the rod, and r = γ r A/(γB) is the dimensionless rotational inertia.

Exact time-harmonic shape functions, mass and stiffness matrices
To facilitate the asymptotic expansion, it is instrumental to identify the 6 constants {C u 1 , C u 2 , C v 1 , ..., C v 4 } , with the degrees of freedom at the rod's ends (i.e. its nodal displacements). This allows a dimensional reduction through a direct application of the compatibility conditions at the joints.
For any given rod of length l, the following notation for the nodal quantities is introduced so that collecting the degrees of freedom at the two ends of the rod in the vectorq = {u 1 , v 1 , θ 1 , u 2 , v 2 , θ 2 } yields the solution of system (5) in the form which is now a linear function of the nodal displacements q(t) =q e −i ω t . The 2-by-6 matrix N (s; ω, P ) acts as a matrix of frequency-dependent and preload-dependent 'shape functions' which is the exact functional basis in which the time-harmonic response of the rod can be represented. The representation (6) can also be considered as the definition of a 'finite element' endowed with shape functions built from the exact solution. Moreover, these time-harmonic shape functions reduce to the quasi-static solution when ω → 0 [so that for vanishing preload, in the limit lim ω→0 N (s; ω, 0) the usual shape functions for beam elements, employed for instance in [28], are recovered]. Note that, in the following, N (s; ω, P ) will simply be denoted as N (s; ω), to simplify notation.
By employing Eq. (6), the exact mass and stiffness matrices of a rod under time-harmonic vibration can be computed. For the k-th rod the kinetic energy and the elastic strain energy are given by where J k and E k are matrices collecting the inertia and stiffness terms, respectively, while B k (s k ; ω) is the strain-displacement matrix, which are defined as follows is a row vector containing the derivative of the shape functions corresponding to the transverse displacement v. Note that the kinetic energy (7a) accounts for translational as well as rotational inertia of the rod. In addition, the contribution of the axial preload can be included in the potential energy as (details are provided in Appendix A) By combining Eqs. (7b) and (8), the potential energy of the k-th rod is denoted as From Eqs. (7), (8) and (9) the frequency-dependent mass and stiffness matrices are naturally defined as The matrices for the quasi-static case can be obtained by evaluating the limit ω → 0. In particular, the quasi-static stiffness matrix lim ω→0 K k (ω) is found to be identical to that reported in Part I, where it is used to formulate the incremental equilibrium and analyze the bifurcation of preloaded lattices.

Equations of motion for the isolated unit cell
As expressions (7)-(9) govern the incremental dynamics of a single rod, the kinetic and potential energies of a single unit cell can be obtained through a summation of contributions from each rod and therefore the equations of motion for an isolated unit cell can be derived from Hamilton's principle applied to the Lagrangian where q is the vector collecting all the degrees of freedom of the unit cell and W(q) = f · q is the incremental work done by the internal actions f acting on the boundary nodes.
The stationarity of the action associated to the Lagrangian (11) leads to the following Euler-Lagrange equations where M (ω) and K(ω) are, respectively, the mass and stiffness matrices of the unit cell. The global matrices can be easily obtained by assembling the local matrices expressed by Eqs. (10) for all the rods in the unit cell.
The assumption that the external action on the isolated unit cell is time-harmonic, f (t) =f e −i ω t , allows a substitution of the time-harmonic response q(t) =q e −i ω t , so to obtain the governing equation for the isolated unit cell with the definition A(ω) = −ω 2 M (ω) + K(ω). Note that the dimension of the linear system (13) is 3N j , where N j is the number of nodes of the unit cell.

Application of the Floquet-Bloch conditions
In order to describe the response of the infinite lattice, the governing equations of the isolated unit cell (13) have to be complemented with the Floquet-Bloch conditions. The application of these conditions to periodic beam lattices is well-known [28,29] and is briefly summarized in the following.
(a) (b) Fig. 3. In order to impose the appropriate Floquet-Bloch conditions required by Eq. (2b), vector q collecting the degrees of freedom of the unit cell is conveniently partitioned by distinguishing the sets of inner q i and boundary nodes, both located at corners {q lb , q lt , q rb , q rt } and on the edges {q l , q b , q r , q t } (a). The corresponding force vector f is partitioned in the same way (b) and the forces acting on the boundary nodes have to satisfy equilibrium as well as the Floquet-Bloch conditions (the shift factor z j = e i k·a j is simply a consequence of Eq. (2b)).
(a) As Eq. (2b) requiresq for all pairs of nodes {p, q} such that x q − x p is an integer, linear combination of the lattice vectors {a 1 , a 2 }, the relations to be imposed on the degrees of freedom of the unit cell derive directly from Eq. (2b), evaluated at x ∈ ∂C and for n j ∈ {0, 1}, hence obtaininĝ succinctly written asq where Z(k) andq * are defined according to Eq. (15a), and z j = e i k·aj , with j = 1, 2. In Eq. (15a) the vector q has been partitioned to denote the inner and boundary nodes according to the notation sketched in Fig. 3. The same partitioning is also introduced for the force vectorf .
A substitution of Eqs. (15) into Eq. (13) provides so that a left multiplication by Z(k) H , where the superscript H denotes the complex conjugate transpose operation 4 , leads to the reduced system where the following definition is introduced Note that the dimension of system (16) is smaller than the dimension of system (13). In fact, the imposition of the Floquet-Bloch conditions has permitted to express the equations of motion of the lattice only in terms of the reduced variablesq * andf * . If external loads are not present in the infinite lattice (so that the boundary forces shown in Fig. 3b are purely internal actions), Eq. (2b) impliesf * = 0 and therefore the following system of equations is obtained, governing Floquet-Bloch wave propagation within the lattice where the matrix of the reduced system has been defined as A * (ω, k) = Z(k) H A(ω)Z(k). Note that this matrix is non-symmetric, but for a conservative system is always Hermitian, so that A * (ω, k) = A * (ω, k) H and ω 2 ∈ R.

Generalized eigenvalue problem for the lattice's dynamics
Eq. (17) defines a homogeneous linear system for the unknown vectorq * , in which the angular frequency ω and the wave vector k are for the moment undetermined.
It is important to note that the matrix-valued complex function A * (ω, k) depends on all the mechanical and geometrical parameters of the unit cell. Since Eq. (17) is homogeneous, all non-trivial solutions are obtained by imposing the condition that defines, although implicitly, the dispersion relation, relating ω and k as ω(k). For every given wave vector k and for each of the corresponding roots ω(k) obtained from Eq. (18), the non-trivial solutions of Eq. (17) provide the modes of the Floquet-Bloch waves propagating through the lattice at each frequency. This means that the vectorq * is an implicit function of ω and therefore of k, so that the dependence on k can be made explicit so that Eq. 17 is rewritten as 3 Incremental dynamics of preloaded lattices: asymptotic analysis The scope of this section is the analysis of the low-frequency/long-wavelength asymptotic behaviour of the lattice dynamics, which will enable the identification of the effective continuum material.

Wave propagation in a prestressed Cauchy continuum
Before introducing the homogenization technique, it is essential to recall here the fundamental equations governing incremental wave propagation in a prestressed hyperelastic continuum. An appropriate form of the equations for the prestressed continuum has to be selected, to result compatible with the formulation of the lattice dynamics introduced in Section 2. Specifically, it is observed that the equations of motion for the lattice are (i) obtained in the context of a linearized theory, and (ii) referred to a preloaded reference configuration, therefore, the dynamics of the unknown 'equivalent' continuum has to be formulated in the context of the incremental theory of nonlinear elasticity by means of a relative Lagrangian description [30]. This is based on incremental constitutive laws relating the increment of the first Piola-Kirchhoff stressṠ to the incremental deformation gradient L = grad uṠ through the elasticity tensor C where T is the Cauchy stress, defining here the prestress, and E is a fourth-order tensor endowed with the left and right minor symmetries and the major symmetry. Moreover, Eq. (21) implies that the number of unknown components of C is at most 9 for a 2D material and 27 for the 3D case.
In the absence of body forces, the incremental equations of motion for the continuum can be written in the usual form where u is the incremental displacement field and ρ h > 0 the mass density. Assuming the usual plane wave representation for incremental displacement, Eq. (22) leads to the following eigenvalue problem governing the wave propagation in a homogeneous Cauchy material whose acoustic tensor A (C) (k) is defined with reference to the unit vector k as A (C) pr (k) = k q C pqrs k s . The eigenvectors a represent the wave amplitudes, while the eigenvalues ω 2 h are the roots of the characteristic equation Note that the subscript h differentiates ω h (k) from ω(k), so that the former defines the angular frequency of a wave propagating through the equivalent homogeneous Cauchy continuum, while the latter the dispersion relation of the lattice.

Asymptotic expansion of Floquet-Bloch waves
A perturbation method is now developed for the equations governing wave propagation in a lattice made of prestressed elastic rods, through a generalization of the technique proposed by Born [8] for lattices involving only point-like mass interactions. A rigorous link is established between the low-frequency solutions of Eq. (19) and spectral characteristics of the equivalent continuum governed by Eq. (23). The important point is that a wave propagating in an homogeneous Cauchy continuum can be seen as a special case of a Floquet-Bloch wave that is characterized by: • a linear dispersion relation, so that writing k = n (with n being a unit vector), the following condition is obtained where ω (1) n depends only on the direction n; • a spatially uniform modulation In general a periodic structured medium does not satisfy these properties for every wave vector k, but its response can be analyzed about the point k = 0 and compared with the constitutive equations of the continuum, Eqs. (25)- (26). In fact, the linear relation (25) can be considered as the first-order term in the asymptotic expansion of the dispersion relation ω( n) centered at = 0 (so that {ω, k} = 0) and along the direction n of the k-space. This asymptotic expansion, truncated at the N -th term is In the above expansion the O( 0 ) term vanishes because the point {ω, k} = 0 certainly satisfies the dispersion equation (18). This follows from the fact that, setting {ω, k} = 0, Eq. (19) becomes so that, since the matrix Z(0) prescribes equal displacements on corresponding sides of the unit cell (preventing rigid rotations), the only possible solutionsq * (0) correspond to 2 rigid-body translations, which represent the nullspace of K(0) without rigid rotations. These rigid translations are essentially the limit of the eigenmodes ϕ k (x) for k → 0. Their derivation requires first the construction of the asymptotic expansion of ϕ k (x), in complete analogy to Eq. (27), and then the computation of the limit → 0. As a result, the zeroth-order term of the waveform ϕ (0) n is indeed uniform in space (independent of x), and therefore the acoustic properties of the equivalent Cauchy continuum have to satisfy conditions which define an 'acoustic equivalence' between the lattice and the continuum. An effective method to obtain the series expansions (27) and (29) is outlined in the following. As the waveform ϕ k (x) for a lattice made up of rods is governed by the vector of degrees of freedomq * , solution of the eigenvalue problem (19), the expansion ofq * (ω(k), k) is performed along an arbitrary direction n in the so that the first termq * (0) n can be used to identify the left-hand side of Eq. (30b). To this end, the matrix so that the eigenvalue problem (19) is rewritten through a substitution of the series representations (31) and (32) as Since the dispersion relation ω(k) is formally inserted in the above expansions, Eq. (33) has to be satisfied for every value of , which means that the left-hand side has to vanish at every order in . Thus the following sequence of linear systems is obtained which has to be solved for the unknown vectorsq * (j) n . It is clear that the computation of these vectors starts from the solution of the zeroth-order equation and then, sequentially, the higher-order terms are to be obtained. At the j-th order, the matrix of the linear system is A * (0) n and multiplies the unknown vectorq * (j) n , so that the constant term (not involving the unknownq * (j) n ) contains all the previously determined vectors {q * (0) n , ...,q * (j−1) n }. Moreover, it is important to observe that the terms A * (j) n in the expansion (32) can be computed explicitly once the series S N n (ω)( ) has been determined. It is recalled that, as shown by Eq. (28), the matrix of each linear system A * (0) n is singular and it has a two-dimensional nullspace spanned by two linearly independent vectors, t 1 and t 2 , which represent the two in-plane rigid-body translations. 5 Thus, every linear combination in the form is a solution of the zeroth-order equation in (34). This implies that the matrix A * (0) n is not invertible, so that the solvability of the j-th linear system depends on the form of its right-hand side, which has to satisfy the following condition, known as the Fredholm alternative theorem or, equivalently, using Eq. (35), the condition In principle, Eqs. (34) and (36) are sufficient to compute the series representations (31) and (27), thus making conditions (30) explicit.

The acoustic tensor for a lattice of elastic rods
The perturbation method outlined in Section 3.2 is general enough to provide, up to the desired order, the series representation of the acoustic properties of a preloaded lattice subject to incremental dynamics.
It will be proved in the following that it is always possible to employ the above-described perturbation technique to construct an eigenvalue problem governing the propagation of waves in a lattice (where elements are subject to both axial and flexural deformation) in the low-frequency and long-wavelength regime. In particular, this eigenvalue problem will possess the following properties: (i) the eigenvalues identify the first-order term ω (1) n of both acoustic branches of the dispersion relation; (ii) the eigenvectors govern the zeroth-order term of the Floquet-Bloch waveform for both acoustic waves, through coefficients {α 1 , α 2 } in the linear combination (35); (iii) the algebraic structure of the problem is exactly equivalent to that governing wave propagation in a Cauchy material, Eq. (23).
The construction of the above eigenvalue problem allows the rigorous definition of the 'acoustic tensor for a lattice of elastic rods' and from the latter the identification of the elasticity tensor representing a material equivalent to the lattice. In fact, this eigenvalue problem defines eigenvalues and eigenvectors satisfying the conditions of acoustic equivalence, Eq. (30).
In order to construct the eigenvalue problem, the solution of the sequence of the linear systems (34) is pursued up to the order O( 2 ). The equations involve the following terms of the series (32) where a series expansion has been introduced for the matrices K(ω( n)), M (ω( n)) and Z( n) as → 0.
It is important to note that: , only the zeroth-order terms of the matrices K and M (which correspond to the quasi-static limit, (ii) the zeroth-order matrix Z (0) , and consequently A * (0) , is independent of the direction n (owing to continuity); while A * (1) n is linear in n, A * (2) n is quadratic in n; (iii) the linear term ω and its solvabilty condition requires two conditions which are always satisfied. In fact, a use of Eq. (37) 2 yields a scalar product which vanishes because Z (0) t i is a rigid-body translation, so that it cannot produce any stress, hence K (0) Z (0) t i = 0. Since Eq. (38) is always solvable, all its solutions can be expressed in the form where t 1 and t 2 are the solutions of the following two linear systems Note that t 1 and t 2 are defined up to an arbitrary rigid-body translation. By employing Eqs. (35) and (39), the linear system of order O( 2 ) reads as which (because A * (0) is singular) admits a solution if and only if the right-hand side is orthogonal to both t 1 and t 2 , namely that in matrix form can be written as 6 Up to order O( 1 ) the coefficients {α 1 , α 2 } and the linear term ω (1) n remain completely arbitrary, but now they have to satisfy the system (41) in order to make Eq. (40) solvable. In fact, the homogeneous system (41) represents an eigenvalue problem with eigenvectors {α 1 , α 2 } and eigenvalues ω with the following definitions Eq. (42) is an eigenvalue problem, and the following properties can be deduced: n andM (0) are real and symmetric, also the matrices Ξ and Γ are real and symmetric, hence the eigenvalues ω (ii) Since Z (0) , t 1 and t 2 are independent of n and Z (1) n , t 1 (n) and t 2 (n) are all linear in n, each component of the 2-by-2 matrix Ξ is a quadratic form in n; 6 Note that vectors t j (n) may contain an arbitrary rigid-body translation. This would apparently lead to a non-uniqueness in the form of Eq. (41), because the terms t j · A * (1) n t i (n) are present. This lack of uniqueness is only apparent, because t j · A * (1) n t i = 0. 7 The resulting expression has been simplified using again the property K (0) Z (0) t i = 0.
(iii) The components of the 2-by-2 matrix Γ admit the following simplifications k=1 γ k l k and |C| are, respectively, the average mass density of the lattice and the area of the unit cell. Note that the matrix Γ contains only terms in the form t i ·M (0) t j , where t i are rigid translations, and therefore the rotational inertia of the rods plays no role (recall the definition (10a)).
Finally the eigenvalue problem (42) can be written in the standard form where α = {α 1 , α 2 } and tensor A (L) (n) reads It is important to note at this stage, that the eigenvalue problem (43) has exactly the same structure of Eq. (23). Furthermore, tensor A (L) (n), the 'acoustic tensor of the lattice', uniquely defines the eigenvalues and eigenvectors appearing on the right-hand side of the equivalence conditions (30), so that ω h = ω (1) n and a = ϕ (0) n = α. This implies that the acoustic equivalence holds if and only if the 'acoustic tensor of the lattice' coincides with the acoustic tensor of the Cauchy material. Therefore, the equivalent Cauchy continuum has to satisfy the acoustic equivalence condition (valid for every unit vector n) It is important to note that the equivalence condition has been obtained without introducing restrictive assumptions on the lattice structure, so that the homogenization method is completely general and includes a generic state of axial preload acting on the lattice. Moreover, the presented technique can easily be extended to three-dimensional lattices.

Identification of the incremental constitutive tensor equivalent to a preloaded lattice
The perturbation method developed in Section 3 leads to the determination of the acoustic tensor of an effective Cauchy continuum, equivalent to the low-frequency response of a preloaded lattice of rods. As the method is entirely based on the dynamics of the periodic medium, the acoustic tensor is obtained directly, without any prior computation of the effective constitutive tensor, which is instead traditional in standard energy-based homogenization techniques [31][32][33][34][35].
In this section the steps for retrieving the incremental (or 'tangent') constitutive tensor C are outlined from the acoustic tensor given by Eq. (45).
As the condition (45) has to hold for an arbitrary direction of propagation, it can equivalently be expressed by applying the Hessian with respect to n on both sides of Eq. (45) to obtain where the right-hand side can be regarded as a data defined by the lattice structure, namely, the Hessian of tensor (44). By considering the symmetry with respect to the {k, l} indices, Eq. (46) provides a linear system of 54 equations in a three-dimensional setting or 12 equations in a two-dimensional setting, while the rank of the system is found to be 26 or 8, respectively. By recalling that the unknown tensor C has the form (21), it is clear that, if the system is solvable, all but one of the unknown components of C can be determined as these are 27 for a three-dimensional lattice and 9 for a two-dimensional. In order to solve for the identification, results obtained by Max Born [8] can now be generalized to prove that: (i) the system is solvable when the equations of motions of the lattice satisfy the rotational invariance and (ii) the solution is unique, except for the spherical part of the prestress (i.e. tr T /3 in 3D and tr T /2 in 2D) which remains undetermined for the system (46), but can be determined by matching the first-order incremental work introduced in Part I.

Bifurcation and loss of ellipticity in a lattice of preloaded elastic rods
In order to demonstrate the effectiveness of the homogenization method developed in Section 3, the condition of loss of ellipticity for a preloaded two-dimensional grid lattice of elastic rods is derived and analyzed. The geometry of the reference configuration is sketched in Fig. 4 and is composed of two sets of rods, inclined at an angle α and characterized by an axial stiffness A 1 = A 2 = A and slenderness ratios Λ 1 = l A/B 1 , Λ 2 = l A/B 2 , where the subscripts 1 and 2 are relative to the horizontal and inclined rods, respectively. For simplicity, the linear mass density is assumed to be the same for both rods γ 1 = γ 2 = γ, while the rotational inertia contribution has been shown in Section 3.3 to be negligible for the homogenization result. The resulting grid of rods is also considered stiffened by diagonal linear springs (Fig. 4), whose stiffness Fig. 4. Current configuration of a rhombic lattice of preloaded elastic rods, with the associated unit cell C highlighted. The direct basis for the lattice is denoted by the pair of vectors {a 1 , a 2 }. Labels 1, 2, and S denote the horizontal rods, the inclined rods, and the diagonal springs, respectively. The stiffness of inclined springs, the axial and flexural rigidity of the rods, the preloads P 1 and P 2 , as well as the grid angle α can all be varied to investigate different incremental responses. k S is conveniently made dimensionless by introducing the parameter κ = k S l/A. In this configuration, the lattice is subject to a preload state defined by the axial forces P 1 and P 2 , made dimensionless respectively as p 1 = P 1 l 2 /B 1 and p 2 = P 2 l 2 /B 2 . Hence, a lattice configuration is completely defined by the parameter set {p 1 , p 2 , Λ 1 , Λ 2 , κ, α}. Note also that the lattice structure considered includes the simpler case of a rectangular grid that has been analyzed in [36].

Macro and micro bifurcations as degeneracies of the dispersion relation
The effect of the diagonal reinforcement (springs labeled with S in Fig. 4) on the bifurcation of the lattice has been systematically investigated in Part I of this study, where it has been shown to play a fundamental role in determining the wavelength critical for bifurcation. Specifically, it has been demonstrated that an increase in the spring stiffness induces a transition of the critical bifurcation from macroscopic to microscopic, and in particular the bifurcation is characterized by an infinite wavelength when κ = 0.
On the other hand, it can be argued that the dynamical setting formulated here may provide a more direct mechanical interpretation of the significant difference between macroscopic and microscopic bifurcation, compared to a static analysis. In fact, as the homogenization scheme introduced in Section 3 proves that the long-wavelength asymptotics for waves propagating in the lattice is governed by the acoustic tensor of the effective medium, Eq. (44), it becomes now clear that a macro-bifurcation in the lattice has to be equivalent to failure of ellipticity in its equivalent continuum. Hence, a macro-bifurcation occurs when the velocity of the acoustic long waves of the lattice vanishes along some directions. Moreover, a clear interpretation of short-wavelength bifurcations (micro-bifurcations) is also provided by the analysis of the dispersion relation of the lattice (18), interpreted now as a function of the preload state. The latter can be used to identify the condition of buckling in the lattice as the 'propagation' of a Bloch wave at vanishing frequency. In fact, regardless of the critical wavelength, macro and micro bifurcations can be visualized by plotting the evolution of the dispersion surfaces along a loading path up to loss of stability.
The essential difference between the two kinds of bifurcation is exemplified in Fig. 5 for two square grids (α = π/2 and Λ 1 = Λ 2 = 10), one without diagonal springs (upper row in the figure) and the other with κ = 0.4 (lower row in the figure), subject to equibiaxial compression (p 1 = p 2 ) of increasing magnitude (from left to right in the figure). The dispersion surfaces (plotted in the non-dimensional space {k 1 l, k 2 l, Ω} with Ω = ω l γ/A and γ r = 0) show that the macro-bifurcation in the grid without springs occurs with the progressive lowering, and eventually vanishing, of the slope of the acoustic branches at the origin, while the dispersion surface attains non-null frequency for every other wave vector. On the contrary, the micro-bifurcation occurring in the grid reinforced with springs is characterized by non-vanishing slope of the acoustic branches at the origin, but instead the preload-induced lowering of the dispersion surface causes the generation of a zero-frequency wave with non-null wave vector (corresponding to a finite wavelength buckling). These dispersion surfaces can be considered the dynamic counterpart of the bifurcation surfaces presented in Part I of this study.
Failure of ellipticity and the consequent emergence of strain localization in the lattice is now investigated following the following steps: (i) the acoustic tensor A (C) (n) of the homogenized continuum is analytically calculated as an explicit function of the parameter set {p 1 , p 2 , Λ 1 , Λ 2 , κ, α}; (ii) loss of ellipticity is analyzed for cubic, orthotropic, and fully anisotropic lattices by identifying the prestress states leading to a vanishing eigenvalue of the acoustic tensor, while the corresponding eigenvector identifies the localization mode; (iii) the ellipticity domain in the {p 1 , p 2 }-space and its dependence on lattice parameters {Λ 1 , Λ 2 , α} is determined; (iv) the predicted directions and modes of strain localization are compared with the maps of displacements resulting in the equivalent elastic material from the application of a concentrated pulsating force (the time-harmonic Green's function), thus following the so-called 'perturbative approach to localization' introduced in [37]; (v) the behaviour of the equivalent elastic solid is validated through numerical simulations (with finite elements) of the forced low-frequency response of the corresponding preloaded lattice.

Acoustic tensor, eigenvalues, eigenvectors, and ellipticity domain
With reference to the basis {e 1 , e 2 }, the acoustic tensor A (C) (n) for the continuum equivalent to the homogenized response of the lattice shown in Fig. 4 is represented (in two-dimensions) as where the components, computed via Eq. (44), are found to be where the components h ijkl are function of the parameter set {p 1 , p 2 , Λ 1 , Λ 2 , κ, α}. Using the same notation as in Part I, the contribution of the rods' grid and the springs are denoted as h G ijkl and h S ijkl , respectively, so that The full expression for the functions h G ijkl and h S ijkl is the following The components h S ijkl , ruling the effect of diagonal springs, can be written as An inspection of the above expressions shows that the dependence on the prestress state is quite complex for the components h G ijkl and that the limit for p → 0 yields the acoustic tensor of the unloaded grid For the special case of a square grid α = π/2 and in the absence of prestress, the acoustic tensor can be further simplified into It is recalled that the homogenized Cauchy material is defined to be strongly elliptic (SE) when the acoustic tensor is positive definite, while ellipticity is defined as the non-singularity of the acoustic tensor, As shown in Section 3.3, the acoustic tensor resulting from homogenization is symmetric, which implies that its eigenvalues are always real. This means that, letting {c 2 1 , c 2 2 } be the eigenvalues of a symmetric A (C) (n), (SE) is equivalent to the strict positiveness of the eigenvalues, c 2 1 > 0, c 2 2 > 0, and (E) is equivalent to the condition of non-vanishing eigenvalues, c 2 1 = 0, c 2 2 = 0 (for all unit vectors n). The objective is now to characterize failure of (E) by studying the eigenvalues of (47) as functions of the preload state applied to the grid-like lattice. To this end, solutions are sought for the following loss of ellipticity condition c 2 1 (n, p) c 2 2 (n, p) = 0 , where n is the usual unit vector defining the direction of propagation and p = {p 1 , p 2 } is a vector simply collecting the preload parameters. In Eq. (51) the dependence on the geometric parameters {α, Λ 1 , Λ 2 } is omitted for brevity and moreover, without loss of generality, it is assumed that c 2 1 ≤ c 2 2 . For every solution {n E , p E } of Eq. (51), the eigenvector g E = g(n E ) associated to the vanishing eigenvalue can be computed. Vectors n E and g E will be respectively referred as the direction (more precisely, the normal to) and mode of the strain localization band.  Table 1. Loss of ellipticity for different geometric configurations of the preloaded grid-like lattice (Fig. 4) in the absence of diagonal springs (κ = 0). The symmetry class is referred to the unloaded configuration. The preload p E and the localization direction n E = cos(θcr)e 1 + sin(θcr)e 2 are obtained by solving the loss of ellipticity condition (51) (assuming a radial path p = {p 1 , p 1 }) and therefore identifies a material state on the elliptic boundary.
ig. 6. Polar plots of the lowest eigenvalue of the acoustic tensor (47), for a prestress state at 50% (dashed gray line) and 100% (continuous blue line) of the limit value for ellipticity loss (p E reported in Table 1). A square lattice and a rhombic lattice are considered with Λ 1 = Λ 2 = 10 and Λ 1 = 7, Λ 2 = 15. The corresponding prestress states are represented on the the elliptic boundary of Fig. 7 as the terminal point 4 of the loading path indicated by the green dashed arrow. In the orthotropic cases, both for squared and rhombic lattices, the loss of ellipticity is characterized by the vanishing of the eigenvalue of the acoustic tensor along two directions denoted as n 1 E and n 2 E , with the associated wave amplitudes reported as g 1 E and g 2 E . The fully anisotropic cases, displays the vanishing of the eigenvalue along a single direction n 1 E . Note that for the rhombic lattice the relative orientation of n E with respect to g E shows that the mode of localization is neither a pure shear nor a pure expansion wave, but a mixing of the two. Conversely, the square lattice always reaches loss of ellipticity through the formation of pure shear bands.
It can be directly verified that in the absence of preload, p → 0, the grid-like lattice considered has c 2 1 > 0 and c 2 2 > 0, so that (SE) holds, except in the case of an extreme material, where the stiffness of the the rods becomes vanishing small, as in [38]. Due to the symmetry of A (C) (n), it follows that, starting from the unloaded state p = 0 with rods of finite stiffness and continuously varying the prestress, the material remains (SE) as long as (E) holds. Therefore, solutions of Eq. (51) are sought as pairs {n E , p E } such that p E represents the terminal point of a path starting at p = 0 and entirely contained in the (SE) domain; in other words, p E is on the boundary of the (SE) domain. The set of these points p E is referred as the elliptic boundary.
In order to explore loss of ellipticity for lattice configurations characterized by different symmetry classes, a square and a rhombic grid are considered, respectively with α = π/2 and α = π/3. For both examples, the following sets of slenderness values are selected Λ 1 = Λ 2 = 10, and Λ 1 = 7, Λ 2 = 15.
In Table 1, for each geometry considered, the first solution of Eq. (51) for equal prestress components p 1 = p 2 are reported together with the associated directions of localization, denoted as n E = cos(θ cr )e 1 + sin(θ cr )e 2 . Note that the symmetry class is referred here to the unloaded configuration, so that the symmetry of incremental response may change as an effect of loading. With the assumed values for grid angle α and slenderness, the cubic, orthotropic, and fully anisotropic cases (10 constants for planar elasticity) can be investigated.
In order to better visualize the direction n E and the associated mode g E , a polar plot of the square root of the lowest eigenvalue c 1 (n, p) is reported in Fig. 6, for the cases listed in Table 1, at two levels of preload, namely 0.5 p E (dashed gray line) and p E (continuous blue line). In Fig. 6a the square lattice with Λ 1 = Λ 2 = 10 is subject to a isotropic prestress in the two directions, p 1 = p 2 , and therefore the cubic symmetry is maintained in the prestressed state. Owing to this symmetry, ellipticity is lost along two orthogonal directions n 1 E and n 2 E . Moreover, the associated wave amplitudes g 1 E and g 2 E are perpendicular to the vector n 1 E and n 2 E respectively, hence indicating that the modes of localization are pure shear waves, the so-called shear bands.
For the orthotropic square lattice (Λ 1 = 7, Λ 2 = 15), the polar plot is given in Fig. 6b. In this case, owing to the orthotropy, waves propagating along the horizontal and vertical direction possess different velocities and therefore ellipticity is lost when the smallest of these velocities vanishes, leading to a single shear band (in this case with a normal n 1 E aligned parallel to the horizontal direction). Quite remarkably, the shear wave responsible for the ellipticity loss is the one propagating along the direction of the 'stiffer' elastic link (having the lowest slenderness), while intuitively a 'shear mechanism' would be expected in the direction of the 'soft' elastic link. This effect will be confirmed and explained further with the computation of the forced response in Section 5. Moreover, it is worth noting that the shear band directions for the square lattice, both cubic and orthotropic, are aligned parallel to the directions of the rods forming the lattice.
For the rhombic lattice with Λ 1 = Λ 2 = 10 and isotropic preload, shown in Fig. 6c, the mechanical behaviour is orthotropic and therefore two directions of localization are obtained. The associated wave amplitudes g 1 E and g 2 E both have respectively a component orthogonal and parallel to the vector n 1 E and n 2 E , so that a 'mixture' of shear and compression waves is involved.
The fully anisotropic version for the rhombic lattice (Fig. 6d) can be obtained by changing the slenderness values (Λ 1 = 7, Λ 2 = 15), so that one of the two localizations is suppressed, while the other is preserved. It is also worth noting that, in contrast to the square case, the directions of localization for the rhombic lattice are not perfectly aligned parallel to the rods' normal, instead, they result slightly inclined, as will be confirmed by the computations reported in Section 5.

The effect of the prestress directionality on the ellipticity loss and localization directions
In the previous section, the loss of ellipticity has been determined for the four lattice's configurations reported Table 1, under the assumption of isotropic preload p 1 = p 2 . This hypothesis is now removed to explore the influence of the 'directionality' of the preload state on the ellipticity loss, in terms of number and direction of strain localizations. The analysis is performed by computing the full set of solutions {n E , p E } of Eq. (51), to determine the elliptic boundary in the {p 1 , p 2 }-space, for several values of rods' inclination α and slenderness λ 1,2 . Fig. 7a and 7b show the ellipticity domains computed, respectively, for the orthotropic (Λ 1 = Λ 2 = 10) and the anisotropic grid (Λ 1 = 7, Λ 2 = 15). The following observations can be drawn.
(i) The elliptic region is unbounded for tensile (positive) preload and bounded when at least one of the preloads is compressive (negative); this is an expected feature, as the contribution of a tensile preload to the potential energy, Eq. (8), is positive definite; (ii) any deviation of the angle α from π/2 reduces the size of the elliptic region; (iii) in the orthotropic case, the domain is symmetric with respect to the bisector p 1 = p 2 (Fig. 7a); (iv) the elliptic boundary appears to be smooth everywhere except at a corner point; (v) the anisotropy induced by different values of slenderness causes the corner to move; the elliptic region is reduced in size along the direction of the smallest slenderness. . The region of positive definiteness of the elasticity tensor C is also reported (orange region). The loading path (gray dashed arrow) refers to condition of isotropic preload p 1 = p 2 on which the markers (in red for the case α = π/2 and green for α = π/3) identify the preload levels used in the computation of the forced response (Section 5). Parts (c) and (d) report the relationship between the direction in the preload space (defined by the angle ψ = arctan(p 1 , p 2 ), see part (a)) and the dyad n E ⊗ g E defined by the angles θcr and θg singling out, respectively, the direction of ellipticity loss and of the localization mode (i.e. n E = cos θcre 1 + sin θcre 2 andg E = cos θge 1 + sin θge 2 ). Notice that the case of localization along two directions (corresponding to the jumps visible in the graphs of part (c) and (d)) only occurs for a prestress state at the corner point of the elliptic boundary, which corresponds to an incrementally orthotropic material behaviour.
In the same figure, the region of positive definiteness (PD) of the elastic tensor, X · CX > 0 ∀X = 0, is also reported. Surprisingly, the (PD) region is found to be independent of the rods' angle α.
The directions of localization in terms of critical angle θ cr are reported in Figs. 7c and 7d for different values of the lattice inclination α, with the purpose of quantifying the influence of the loading path (parameterized by means of the angle ψ = arctan(p 1 , p 2 )). These plots show that only one localization occurs, except for a loading direction corresponding to a corner point on the elliptic boundary, where two simultaneous localizations occur. Even though this is expected for the orthotropic case (due to symmetry), the anisotropic grid can also exhibit a sort of 'orthotropic' response (displaying two localizations) when a proper value of preload p 1 = p 2 is applied (along the dashed black line of Fig. 7b). Furthermore, the angle between the vectors n 1,2 E and the horizontal axis ( Fig. 7c and 7d), which identifies the direction of the localization bands, depends weakly on the prestress (except at the corner point of the elliptic boundary), namely, the directions of n 1,2 E follow the orientations of the grid's ligaments, even if they result slightly misaligned.

Time-harmonic forced response near the elliptic boundary
The analysis of the homogenized continuum, equivalent to a preloaded grid of elastic rods (presented in Section 4) predicts that the incremental response can display strain localizations due to prestress-induced loss of ellipticity. In this section, the homogenization is validated through comparisons between the actual low-frequency forced response of the lattice (simulated numerically with a finite element technique) and the time-harmonic Green's function (which can be found in [25]) for the equivalent continuum, both computed at increasing levels of the preload, so that the elliptic boundary is approached. The comparison is performed for the four geometric configurations reported in Table 1 and for four prestress levels, namely {0, 0.8, 0.9, 0.99} p E (with p E being the prestress state leading to ellipticity loss), which are marked on the loading path of Fig. 7a and 7b. The lattice response is numerically analyzed using the COMSOL Multiphysics finite elements program in the frequency response mode. A square finite-size computational window with a width of 350 unit cells (of dimension 350l 1 , with l 1 denoting the cell edge) is considered, with a perfectly matched layer (PML) along the boundaries, so that here waves are not reflected, rather absorbed, and the response of an infinite body is simulated.  Fig. 8. Slowness contours and dispersion surfaces for a rhombic anisotropic lattice with Λ 1 = 7, Λ 2 = 15 (in blue) and for the effective Cauchy continuum (in red) at frequency Ω = 0.01. The evolution of the contours and dispersion surfaces induced by a compressive preload of equal components p 1 = p 2 along two inclined directions (the considered loading path is depicted in Fig. 7b) demonstrates that the nonlinear dispersion of the lattice is negligible, except when the material is very close to the elliptic boundary, namely, for a prestress above 0.9 p E . The comparison between the behaviours of the lattice and its equivalent continuum shows a great agreement.
A pulsating concentrated force, in-plane applied, is considered acting at the center of the computational domain. For a given load, the complex displacement field u(x), with horizontal and vertical components u 1 (x 1 , x 2 ) and u 2 (x 1 , x 2 ), is computed and the results are plotted in terms of the modulus of the displacement associated to its real part only, δ R (x 1 , x 2 ) = ( u 1 ) 2 + ( u 2 ) 2 (the plots of the imaginary part of the displacement is omitted for brevity).
In all the following analyses the frequency of the pulsating force is set to be Ω = ω l γ/A = 0.01, a low value providing a reasonable match, in term of acoustic properties, between the effective continuum and the lattice. In fact, it is important to note that the mismatch between the two is different from zero for any non-vanishing frequency, even though vanishing in the limit Ω → 0. In particular, as the elliptic boundary is approached, this mismatch is expected to become wider for those waves propagating along the direction of ellipticity loss. This is easily explained by the fact that, as the linear term in the dispersion relation tends to vanish (in a direction n E ), the nonlinear dispersion of the lattice becomes non-negligible for any non-vanishing frequency.
By considering for instance the rhombic anisotropic grid (Λ 1 = 7, Λ 2 = 15), the deviation between the responses of the lattice (reported in blue in Fig. 8) and its equivalent continuum (reported in red in Fig. 8) can be visualized in terms of slowness contours computed at the frequency Ω = 0.01. By comparing the contours for the four preload states, it can be appreciated that these are superimposed up 0.9 p E , so that the nonlinear dispersion of the lattice becomes non-negligible only when the material is very close to the elliptic boundary, namely, at a preload 0.99 p E , and only for waves close to the direction of ellipticity loss. It is also worth noting that when p ≈ p E , the slowness contour of the lattice (reported in blue) is always contained inside the contour of the continuum (reported in red), so that the nonlinear dispersion implies that waves speeds are slightly higher for the lattice than for the effective Cauchy medium.

Square lattice
Cubic and orthotropic square grids are considered, subject to a pulsating diagonal force (inclined at 45 • with respect to the rods' axes), with the purpose of revealing the emergence of strain localizations aligned parallel and orthogonal to the rod axes.  The case of cubic symmetry (Λ 1 = Λ 2 = 10) is analyzed in Fig. 9, where the displacement field, numerically computed for the square grid (subject to a pulsating concentrated force, upper row), is compared to the response of the homogenized continuum (analytically solved with applied the same concentrated force, lower row), for four values of prestress (increasing from left to right, p 1 = p 2 = {0, −3.347, −4.891, −5.380}). As the elliptic boundary is approached, the emergence of two strain localizations becomes evident and confirms the predicted vanishing of an eigenvalue (wave speed) reported in Fig. 6a. Fig. 10. Deformed configurations of a square lattice (Λ 1 = Λ 2 = 10) near (the zone is indicated in Fig. 9d) the point of application of a pulsating force (denoted with a red arrow) at a level of prestress close to the elliptic boundary (p = 0.99 p E ). The pattern shows a motion resulting from the superposition of two shear localizations induced by a pulsating diagonal load.
By taking different snapshots of the displacement map at different instants of time during the dynamic response of the grid and near the point of application of the pulsating force, it is possible to detect the actual localization mode activated by the applied diagonal force. In fact, deformed configurations calculated in the grid (the zone is indicated in Fig. 9d) through a finite element simulation and plotted in Fig. 10 reveal a characteristic motion resulting from the superposition of two shear localizations emanating from the loading point.   Fig. 7b). Note the emergence of a single vertical strain localization, actually a shear band, as predicted in Fig. 6b.
The square grid displays a strain localization into a single shear band when the slenderness values of the two orthogonal elastic links are set to be different, thus breaking the cubic symmetry, but preserving orthotropy. The response of the grid with Λ 1 = 7 and Λ 2 = 15 is reported in Fig. 11 for four preload states corresponding to p 1 = p 2 = {0, −1.657, −1.864, −2.050}. As already revealed by Fig. 6b, a single vertical shear band emerges, thus confirming the counter-intuitive result obtained in the previous section, namely that the shear wave responsible for the ellipticity loss is the one propagating along the direction of the 'stiffest' elastic link (which possesses the lowest slenderness). The mechanism underlying this effect is revealed by analyzing the actual deformed configuration of the grid reported in Fig. 12 (plotted at different instants of time and obtained via f.em. simulations). The figure, which refers to the zone indicated in Fig. 11d, reveals that the vertical strain localization emerges from a prevalent bending deformation of the 'soft' vertical links accompanied by an approximately rigid rotation of the 'stiff' horizontal rods, which is allowed by a large rotation at the nodes. Fig. 12. Deformed configurations of the square lattice (Λ 1 = 7, Λ 2 = 15) near (the zone is indicated in Fig. 11d) the point of application of a pulsating concentrated force (denoted with a red arrow) at a level of prestress close to the elliptic boundary (p = 0.99 p E ). The pattern shows the characteristic motion of a single vertical shear band localization induced by the pulsating load. It can be noticed that the localization emerges from a prevalent bending deformation of the 'soft' vertical links and an approximately rigid rotation of the 'stiff' horizontal rods allowed by a significant rotation of the junctions.
The comparison between the responses of the grid and of the homogenized continuum that has been presented in Figs. 9 and 11 shows an almost perfect agreement from low to high prestress levels, up to (a) Cubic square lattice (Λ 1 = Λ 2 = 10) (b) Orthotropic square lattice (Λ 1 = 7, Λ 2 = 15) Fig. 13. Fourier transform of the complex displacement fields of the cubic (a) and orthotropic (b) square lattice subject to a diagonal force in the center and a prestress state of p = 0.99 p E (i.e. close to the elliptic boundary). The slowness contours of the lattice (dashed green) and the effective continuum (dashed red) are superimposed to highlight the Bloch spectrum of waves excited by the forcing source. Note that the strong focus of the spectrum indicates that few plane waves, namely those 'slow' waves that are close to cause the ellipticity loss, prevail on the response as it expected for a material near the elliptic boundary. In fact it is worth noting that for the orthotropic grid (b) the waves propagating vertically remains mostly inactivated compared to the ones propagating horizontally.
values close to the elliptic boundary. Nonetheless, the slowness contours reported in Fig. 8 shows a mismatch between the responses of the lattice and of the equivalent continuum, which becomes clearly visible when the lattice is close to the failure of ellipticity. This can be revealed by considering the lattice's complex displacement field (reported in the last column of Figs. 9 and 11, prestressed at p = 0.99 p E ), computing its Fourier transform and superimposing this to the corresponding slowness contour. This is reported in Fig. 13, where the Fourier transform shows that, for both considered square grids, the Bloch spectrum of waves excited by the diagonal load matches the slowness contour of the lattice (reported in green) and is also highly focused around the directions of ellipticity loss where it is at the maximum distance from the contour of the continuum (reported in red). It is worth noting that the strong focus of the spectrum confirms the fact that few plane waves, namely those 'slow' waves that are close to cause the ellipticity loss, prevail on the response, as it is expected for a material near the elliptic boundary.

Rhombic lattice
In the previous section, the square lattice was shown to display only localizations in the form of pure shear bands perfectly aligned parallel to the elastic ligaments. However, on the basis of the analysis performed in Section 4, the formation of localizations is expected along different directions and with different deformation modes when a rhombic grid is considered, α = π/2.
In order to investigate the response of the orthotropic (Λ 1 = Λ 2 = 10) and anisotropic (Λ 1 = 7, Λ 2 = 15) rhombic lattices (α = π/3), both horizontal and vertical concentrated forces will be considered, so to observe a dependence of the number of strain localizations on the loading orientation. Furthermore, in contrast to what happens in the square grid, the directions of localization are expected to occur with a slight misalignment with respect to the directions of the rods, as predicted in Figs. 6c and 6d.
In Fig. 14 Figs. 14m-14p, the agreement between the lattice response and its homogenized continuum becomes evident at each state of lattice's preload. With reference to a prestress state p = 0.99 p E (last column on the right of the figure), while two localization bands are activated by the vertical force, only one is generated by the horizontal force.
Note also that the slight misalignment between the localization direction and the rod angle remains hardly visible until the material is close to elliptic boundary (compare for example the case of vanishing prestress, Fig. 14a, to the case p = 0.99 p E , Fig. 14d).
The localization modes are analyzed in Fig. 15 by inspecting the lattice deformation computed via f.e.m.   Fig. 7a). The same comparison is reported in parts (i)-(l) and (m)-(p) for a pulsating vertical force. Note that, even though ellipticity is lost along the two directions predicted in Fig. 6c, the activation of strain localization depends on the orientation of the load, so that two bands are activated by the vertical force, while only one is generated by the horizontal load. Furthermore, note that the directions of the localization bands (with angles of the normal equal to θcr = 88.2 • , 151.8 • ) is slightly misaligned with respect to the rod's inclination.
at different temporal instants through snapshots taken in the neighbourhood of the loading point [two zones are considered, which do not include the concentrated force and are shown in Fig. 14 parts (d) and (l)]. By comparing the localization band induced by the horizontal load (upper row of Fig. 15) to the one generated by the vertical force (second row in Fig. 15), it is clear that the horizontal band is characterized by an almost perfectly straight wavefront, while the inclined band displays a periodic modulation along the front. This modulation is due to the superposition of the two localization patterns that are activated by the vertical force, where the inclined band prevails on the almost horizontal one, as can be seen in Figs. 14l and 14p. The relative contribution of the two localizations can be further investigated through a Fourier transform of the lattice response, to be compared with the Bloch spectrum generated by the forcing source. Fig. 16 shows the Fourier transform of the field generated in the rhombic grid when the material is close to the elliptic boundary (p = 0.99 p E ). Fig. 16a and 16b correspond, respectively, to the Fourier transform of Fig. 14d and Fig. 14l. The two sharp peaks of Fig. 16a clearly show that the source is emanating pure plane waves t = 0  The pattern on the first (the second) row shows the motion of the localization induced by the pulsating horizontal (pulsating vertical) load. Note that the effect of the load is the generation of bands almost parallel to each other and possessing almost constant amplitude, except for a modulation of the inclined localization (shown in the lower row). In both cases the deformation mode is mostly of shear-type even though an 'expansion' component is also present, as indicated by the vectors g E , Fig. 6c.  propagating almost vertically (θ cr = 88.2 • ) (the slight tilt exactly matches the sub-horizontal wavefronts of the response). Instead, the four peaks of Fig. 16b demonstrate that two families of plane waves are activated: the prevailing ones propagate along the inclined direction (θ cr = 151.8 • ) while the vertically-propagating waves result dimmer (in agreement with the response of Fig. 14l).
For a complete anisotropic material only a single localization is expected to occur and in the case of the   Fig. 7b). Note that, even though ellipticity is lost along the direction predicted in Fig. 6d (with normal angle θcr = 151.4 • ), the activation of the strain localization depends on the orientation of the applied pulsating force, so that the band is activated by the vertical force while localization is inhibited when the grid is loaded horizontally.
anisotropic grid considered the localization direction has been predicted in Fig. 6d to occur at an inclination θ cr = 151.4 • of the band normal. However, similarly to the case of the orthotropic grid, the activation of the localization depends on the orientation of the perturbing force. This can be observed by comparing the lattice response generated by a horizontal and a vertical pulsating force, both reported in Fig. 17 and showing that strain localization is absent when a horizontal force is applied, regardless of the prestress level (see Figs. 17a-17h). On the other hand, the vertical concentrated force triggers an inclined localization when the material is brought close to ellipticity loss (see Figs. 17i-17p).
Results reported in Figs. 18 and 19, referred to the anisotropic rhombic lattice have been obtained with the same setting of Figs. 15 and 16, referring to the orthoptropic case. Fig. 18 shows that, as only one localization band is present, the deformation pattern is characterized by the generation of essentially straight wavefronts propagating outwards from the localization band. The generation of these parallel waves is perfectly captured by the sharp peaks in the Fourier transform of the lattice response, reported in Fig. 19. Fig. 19a and 19b correspond, respectively, to the Fourier transform of Fig. 17d and Fig. 17l. The two sharp peaks of Fig. 19b, sitting at the long tips of the contour aligned parallel to the direction θ cr = 151.4 • , clearly shows that the response induced by the vertical force involves pure plane t = 0 Detail of Fig. 17l t = π 2ω t = π ω t = 3π 2ω Fig. 18. As for Fig. 15, except that an anisotropic rhombic lattice (Λ 1 = 7, Λ 2 = 15) is considered, for four different levels of preload (the path is shown in Fig. 7b). The zone reported in the figure is shown in Fig. 17l. The pattern shows the motion of the localization induced by the pulsating vertical load, which generates waves emanating outwards from the localization band, parallel to each other, and possessing an almost constant amplitude. The deformation mode is mostly of shear-type even though an 'expansion' component is also present, as indicated by the vector g E reported in Fig. 6d. waves propagating with fronts inclined at θ cr = 151.4 • − 90 • = 61.4 • with respect to the horizontal axis.
It is also important to note that waves do not propagate vertically when the load is vertical while these are the only waves composing the response triggered by a horizontal force (see the peaks on the short tips of the contour in Fig. 19a). This is in agreement with the fact that localization is not generated by the horizontal force and the 'slow waves' leading the homogenized continuum to ellipticity loss remain inactive.

Conclusions
An analytic formulation has been developed for the dynamics of lattices of elastic rods (with distributed mass density and rotational inertia) subject to axial forces of arbitrary amount, incrementally loaded in the plane. This formulation leads, through an asymptotic expansion of Floquet-Bloch waves, to a low frequency approximation for an equivalent Cauchy elastic material which turns out to be coincident with that determined on the basis of an energy match (see Part I of this study). The developed technique can be employed to systematically analyze arbitrary lattice geometries and preloaded configurations, therefore predicting both local and global material instabilities, in other words micro-buckling and strain localizations. Loss of ellipticity has been analyzed for a skewed grid-like lattice, to (i) explore cubic, orthotropic and fully anisotropic homogenized material responses, (ii) compute the elliptic domain of the homogenized continuum as a function of lattice parameters, (iii) analyse the structure of the acoustic branches close to ellipticity loss, and (iv) investigate forced vibrations (both in physical and Fourier spaces) revealing low-frequency wave localizations. The homogenization is proven to yield a superb approximation to the behaviour of the rod's lattice, except when micro-buckling occurs, which remains undetected in the equivalent continuum.
[37] D. Bigoni  A Linearized equations of motion for an axially prestretched elastica The linearized dynamics of an axially stretchable Rayleigh elastic rod can be obtained through a linearization (around a stretched equilibrium configuration) of the equations governing the dynamics of large deflections and flexure of an elastic rod.
As the first-order term of (A.3a) has to vanish when the configuration {u 0 , v 0 } = {(λ 0 − 1)x 0 , 0} satisfies equilibrium, the prestretch λ 0 is the solution of the condition ψ λ (λ 0 )−P = 0, indicating that the applied load P is indeed equal to the axial preload. Moreover, it is important to note that the second-order term of (A.3a) involves the strain energy functions only in terms of second derivatives, ψ λ (λ 0 ) and ψ χ (0), evaluated on the straight stretched configuration.
It is now instrumental to update the reference configuration from the stress-free configuration to the stretched configuration, so that the second-order functionals (A.3) can be adopted to govern the incremental response of the rod. This can be performed by changing the variable of integration from where l = λ 0 l 0 is the current rod's length and the symbol has to be understood as differentiation with respect to s 8 , while ψ λ (λ 0 )λ 0 = A(λ 0 ) is the current value of axial stiffness, and ψ χ (0)λ 0 = B(λ 0 ) the current value of bending stiffness, both functions of the axial stretch λ 0 . As these second-order functionals have been derived from the large deformation beam theory, they describe the correct incremental response superimposed upon a give pre-stretched state. Therefore, the correct form of the equations of motion governing the incremental displacements can be derived employing the following functionals where now the fields {u(s, t), v(s, t)} are the current incremental fields. Note that the linear mass density appears divided by the prestretch, indicating that the current density governs the incremental inertia of the rod, in fact the conservation of mass requires γ 0 /λ 0 = γ(λ 0 ), where γ is the current linear mass density of the stretched rod. Similarly, the current rotational inertia is denoted as γ r,0 /λ 0 = γ r (λ 0 ). The governing equations (3) are directly obtained through the application of Hamilton's principle on the Lagrangian L = T − V constructed using the second-order functionals (A.5).