On the Failure of Classic Elasticity in Predicting Elastic Wave Propagation in Gyroid Lattices for Very Long Wavelengths

In this work we investigate the properties of elastic waves propagating in gyroid lattices. First, we rigorously characterize the lattice from the point of view of crystallography. Second, we use Bloch–Floquet analysis to compute the dispersion relations for elastic waves. The results for very long wavelengths are then compared to those given by classic elasticity for a cubic material. A discrepancy is found in terms of the polarization of waves and it is related to the noncentrosymmetry of the gyroid. The gyroid lattice results to be acoustically active, meaning that transverse waves exhibit a circular polarization when they propagate along an axis of rotational symmetry. This phenomenon is present even for very long wavelengths and is not captured by classic elasticity.


Introduction
Architectured materials are those that possess an inner geometry [1]. This multi-scale spatial arrangement of the constitutive materials allows for achieving mechanical properties that are not present in bulk material itself [2]. Although this appears to be an engineering-based approach to materials design, it should be noted that this strategy is, in fact, central in nature where biomaterials must perform many functions from a small and limited set of elementary chemical elements [3,4]. Therefore, to enhance some target properties, regular patterns often emerge. The best-known example is the honeycomb, where bees need to maximize the volume of the cells while minimizing the quantity of matter (wax) used [5]. Another example is the iridescent color of the wings of some butterflies. This phenomenon is due to the non-centrosymmetric mesostructure of the material constituting the wings, which acts a photonic crystal [6,7].
The study of elastic waves propagating in architectured materials is of particular interest since unconventional effects due to the local organization of the matter can emerge on a macroscale. In order to study these phenomena adequately, two points of view can be adopted. Either to describe all the details of the architecture or to consider an effective continuum as replacement. The first option is very general since no particular modeling assumptions are involved. However, since the inner geometry of the material has to be explicitly described and meshed, the numerical cost is often prohibitive for actual applications. Moreover, the computed solution often contains many unnecessary details for practical use. The second option, which is based on elastodynamic homogenization [8][9][10] amounts to substitute the initial heterogeneous material by an equivalent homogeneous continuum. This equivalence is only valid under specific assumptions on the range of variation of some intrinsic parameters, and hence more restrictive than the first approach. However, within the validity domain Several works can be found in the literature that investigate 2D chiral elasticity and focus on their unusual mechanical properties such as negative Poisson ratio [24,25]. Concerning wave propagation, these architectures have been extensively studied [26][27][28] and the need for a generalized continuum theory in order to capture the onset of dispersion and anisotropy at higher frequencies has been pointed out [15,16]. In all these cases, the unit cells under investigation are chiral and centrosymmetric. Again, it is worth noting that such a combination is possible pnly in 2D. When moving to 3D, the picture becomes more complex and due to habits coming from the 2D situation, an ambiguity between the two definitions can be usually found in the literature. For instance, the well-known and studied in-plane hexachiral and tetrachiral patterns [29,30], are no more chiral once extruded in 3D.
Concerning wave propagation in non-centrosymmetric materials, if the phenomenon can be studied in 2D [15], the effects become even more interesting in 3D as the polarization of waves are then richer. As a consequence, interest in 3D non-centrosymmetric metamaterials have recently emerged for their strech-twist coupling [31] or their acoustical activity [32]. The effects related to size-dependent properties and characteristic lengths are also exploited and investigated in [33] and a micropolar generalized model is used to investigate acoustical activity in [32].
The present work is focused on the features of elastic wave propagating in non-centrosymmetric architectured materials. Among them, those based on gyroid architectures are probably the most commonly used. In electromagnetics they are widely studied as metamaterials [6], in acoustics as phononic crystals [34], and in biomechanics as bone substitutes [35,36]. In this paper, we will highlight a particular situation for which the solution predicted by classic continuum mechanics is wrong even for very long wavelengths. It is important to note that the sensitivity of the mechanical behavior to the lack of centrosymmetry can also manifest in statics [37][38][39].
The paper is organized as follows: In Section 2 the gyroid lattice is described. In Section 3 the Bloch-Floquet analysis is introduced along with some necessary definitions for polarization studies. Dispersion analysis is performed and discussed in Section 4. Section 5 compares the results from Section 4 to those obtained in the LW-LF approximation. Finally, some conclusions are drawn in Section 6.

Notations
Throughout this paper, the Euclidean space E 3 is equipped with a rectangular Cartesian coordinate system with origin O and an orthonormal basis B = {e 1 , e 2 , e 3 }. Upon the choice of a reference point O in E 3 , the Euclidean space and its underlying vector space E 3 can be considered as coincident. As a consequence, points will be designated by their vector positions with respect to O. For the sake of simplicity, E 3 will, from now on, simply be denoted E . In the following, r will designate the position vector of a point P, and, with respect to B, When needed, Einstein summation convention is used, i.e., when an index appears twice in an expression, it implies summation of that term over all the values of the index. The dot operator (·) stands for the scalar product, the ∧ for the cross product, and δ ij is the Kronecker delta. Moreover, the following convention is retained: • Blackboard fonts will denote tensor spaces: T; • Tensors of order > 1 will be denoted using uppercase Roman Bold fonts: T; • Vectors will be denoted by lowercase Roman Bold fonts: t.
The orthogonal group in R 3 is defined as O(3) = {Q ∈ GL(3)|Q T = Q −1 }, in which GL(3) denotes the set of invertible transformations acting on R 3 .

The Gyroid Lattice
The gyroid is a triply periodic minimal surface introduced for the first time in [40,41]. Since it is a minimal surface, it has a zero mean curvature, meaning that every point on the surface is a saddle point with equal and opposite principal curvatures [42]. It is periodic with respect to three orthogonal space vectors and is chiral, meaning that the surface only possesses rotation symmetry elements or, equivalently, that it does not possess any symmetry plane nor symmetry center [43].

Parametrization of the Gyroid Lattice
The gyroid's morphology is usually described using a level surface given by the following equation: with a and b real parameters and x, y, z coordinates of the position vector r.
In this paper, we will focus on a particular gyroid defined by : φ(x, y, z; a, b) := sin 2πax cos 2πay + sin 2πay cos 2πaz + sin 2πaz cos 2πax − b, that exists only in the range |b| < √ 2 [43]. Indeed, beyond this peculiar value of b, it is possible to show that the surface described by φ presents discontinuities located at the borders of the fundamental (or asymmetrical) unit cell. A proof of this geometric constraint is provided in Appendix C. Due to its chiral nature, the gyroid surface exists in two enantiomorphic forms: dextrogyre and levogyre. The surface described by the implicit Equation (1) will be, arbitrarily, chosen to be the dextrogyre form. The levogyre form of the implicit equation is easily obtained by applying, for instance, the transformation x → −x in Equation (1) . From the definition of the surface, one can then obtain a volume by defining the presence of matter for points satisfying the following inequality: The parameter a controls the spatial period while the parameter b controls the porosity p, defined as the unity minus the ratio between the volume of the gyroid lattice and that of the unit cell. Examples of unit cells of such solids obtained with different values of b are plotted in Figure 1.
Despite what the angle of view may suggest, all these structures are simply connected.
The relationship between the porosity p and the parameter b is not analytical but can be estimated numerically. As can be observed in Figure 2 this relationship is almost linear and for porosities between 0.2 and 0.8, the following linearized formula can be used to estimate the porosity:

Symmetry Properties
Thanks to its triply periodic nature, the gyroid structure can be considered as a crystal [34] with Body Centred Cubic (BCC) Bravais lattice and point group O, using group theoretic notation [44], or 432 using the Hermann-Maugin notation. It should be noted that this cubic group only contains rotations geometry of the gyroid is hence chiral. To be more specific there are 3 different cubic point groups: O, O − , and O ⊕ Z c 2 . The first one just contains rotations, the group is hence chiral and non-centrosymetric. The second, O − , possesses symmetry planes but not the inversion, the group is achiral and non-centrosymetric. The last group O ⊕ Z c 2 is centrosymmetric hence achiral. Some details are provided in the Appendix A, and more information can be found in [23]. The symmetry of the spatial structure is described by the space group, which details how transformations from the Bravais lattice and the point group are combined in the actual crystal. The space group (S G) of the gyroid crystal is, using Hermann-Maugin notation, I4 1 32 (space group #214 in the International Tables of Crystallography [45]), where the I stands for Body-Centered (BC), meaning that the conventional unit cell defined in crystallography is not primitive, but body-centered (more details provided in Section 2.3). This space group contains screw axes and, as such, is not symmorphic. A space group is called symmorphic if, apart from the lattice translations, all generating symmetry operations leave one common point fixed. Permitted as generators are thus only the point-group operations: Rotations, reflections, inversions, and rotoinversions. The symmorphic space groups may be easily identified because their Hermann-Mauguin symbol does not indicate a glide or screw operation. If C stands for the "crystal" structure, and for the group action as defined in Appendix B: ∀r ∈ C, ∀g ∈ S G, r = g r ∈ C.
From the generating transformations defined in Appendix B and using the equation of the gyroid surface (c.f. Equation (1)) it is straightforward to verify that: ∀g ∈ S G, φ(g r) = φ(r).

Unit Cell
Due to the periodicity of the geometry, the study of the gyroid structure can be restricted to a unit cell: ∀r ∈ E, ∃(r 0 , t) ∈ (T , R) , r = r 0 + t (4) where T designates a unit cell and R a periodicity lattice. Note that for a given lattice R, the choice of T is not unique. It can be convenient to chose a reference unit cell from which other unit cells will be defined using lattice vectors a i and a triplet (n 1 , n 2 , n 3 ) combined as follow: t = n 1 a 1 + n 2 a 2 + n 3 a 3 , n i ∈ Z. The triplet (0, 0, 0) is then associated to the reference unit cell. Once a lattice basis chosen, the considered unit cell is defined as follows, and the associated periodicity lattice is given by: R C = {t ∈ E|t = n 1 a 1 + n 2 a 2 + n 3 a 3 , n i ∈ Z}. These geometrical sets can be described using lattice vectors a i , gathered into a basis as B = {a 1 , a 2 , a 3 }. Note that they can as well be described with respect to the basis B = {e 1 , e 2 , e 3 } of E.
Among all the possible unit cells, some are special and have been given a standard name in the crystallography community: The Conventional Unit Cell (CUC) and Primitive Unit Cell (PUC).
The CUC of the a BCC lattice is depicted in Figure 3. It is defined as the smallest cell having its edges along the symmetry directions of the Bravais lattice. Notice that, for BC lattice, this unit cell is not minimal and a so-called PUC can be considered instead. For a continuous structure tilling of the space, the primitive unit cell is defined as the smallest tile that generates the whole tiling using only translations. As such the primitive unit cell is a fundamental domain with respect to translational symmetries only.

BCC Conventional Unit Cell
For a BCC lattice, the conventional unit cell is defined as depicted in Figure 3. As its faces are perpendicular to Bravais lattice directions, despite its non minimality, this unit cell is easy to use for numerical computations. In this case, the conventional lattice vectors a i , are chosen such that a i ∧ e i = 0.

BCC Primitive Unit Cell
For a BCC lattice, two possible PUC are represented in Figure 4.  The primitive lattice vectors a i are not unique and the ones for the PUC depicted in Figure 4a are defined as: while the ones presented in Figure 4b are defined as: They form P = {b i } 1≤i≤3 and P = {b i } 1≤i≤3 two other bases of E 3 , which metric tensors are given by: Being defined by a more symmetrical set of vectors P, only the first primitive unit cell will be considered here after.

Reciprocal Basis and Brillouin Zone
The vector space dual to E is symbolized by E , and is formally defined as the space of linear forms on E, ∀l ∈ E , ∀u ∈ E, l(u) = α ∈ R.
Upon the choice of a scalar product the two spaces can be identified: and from a basis of E a basis of E can be constructed. In the field of physics, E corresponds to the space of wavevectors, and a generic element of E is denoted by k. For our applications, it is fundamental to introduce the reciprocal lattice R of R: The vectors (a 1 , a 2 , a 3 ) constitute the lattice basis B of R and verify: B can be computed from any lattice basis (a 1 , a 2 , a 3 ) of R according to: Due to the following property, vectors of the reciprocal lattice are the supports of R-periodic functions on E since, In addition, use will be made of the First Brillouin Zone (FBZ) T of the reciprocal lattice R defined as: Using the reciprocal lattice R and the FBZ T , any wavevector k can be expressed as: We can geometrically interpret T as the set of wavevectors k which are closer to the null wavevector than to any other wavevector ξ of the reciprocal lattice R . It is the Wigner-Seitz cell of the reciprocal lattice, this cell is uniquely defined and independent of the choice of T . Similarly to the primitive unit cell in the direct space, the FBZ is a fundamental domain with respect to translations. Physically, the wavelength λ is defined as the inverse of the wavenumber, which is the norm of the wavevector : λ = 1/ k . Then, wavevectors belonging to T have wavelengths that are greater than the wavelength of the periodicity lattice. When k → 0 the wavelength becomes infinite, solicitations varying with almost null wavenumber are said to be scale separated with respect to the periodicity lattice. This is usually the regime in which the LW approximation of elastodynamics homogenization holds.
The FBZ can be further reduced if we consider also the symmetry operations of the point group. The result is an Irreducible Brillouin Zone (IBZ), that is delimited by points of high symmetry, summarized for the considered gyroid lattice in Table 1. In this table, the high symmetry points are given in the non-orthogonal reciprocal lattice basis P , dual to the primitive lattice basis P as well as in the orthonormal lattice basis B which coincides with its dual in the reciprocal space. The path obtained connecting these high symmetry points along the edges of the IBZ is often used to characterize the photonic and phononic properties of the lattice [6,26]. However, it has been pointed out that this choice is not always reliable as some relevant information, e.g., about band gaps, could be missing [46].
In this paper we consider the basis P, that corresponds to the one depicted in Figure 4a. The reciprocal lattice is itself a Bravais lattice, and in the case of BCC lattice, it is a Face Centered Cubic (FCC) lattice. Using Equation (6), the reciprocal basis P is equal to: The metric tensor of P is the inverse of the one of P:

Analysis Tools
In the previous section we characterized the lattice from the point of view of crystallography. In this next section we will use these results to compute the elastodynamic response of the lattice. The objective of the present section is to provide the analysis tools to be used to perform the computation and to interpret the results.

Bloch-Floquet Analysis
Since the material is periodic, the dispersion diagram will be computed using Bloch-Floquet analysis [47]. The elastodynamics equation for the periodic continuum reads: where ρ(r) is the R-periodic mass density and C(r) is the R-periodic fourth-order elasticity tensor. As we saw in Section 2.3, each cell of the assembly can be identified by the triplet (n 1 , n 2 , n 3 ), where the triplet (0, 0, 0) is conventionally assigned to the reference unit cell. The position of a point r of the (n 1 , n 2 , n 3 )-cell is obtained from the position of a point in the reference unit cell r 0 by Equation (4) where t = n p a p . Being R-periodic, ρ(r)and C(r) verify: Thanks to the Bloch-Floquet theorem [47], elementary solutions to the Equation (9) over C can be searched for in the form of Bloch-waves : where U k is the complex polarisation vector which is R-periodic in space and constant in time, f is the frequency of the Bloch-wave, and k its wavevector U k describes the movement of matter as the wave propagates. It is important to remark that wavevectors k follow the so-called "crystallographer's definition" which consists in dropping the often seen 2π coefficient. This implies, for instance, that the norm k , i.e., the wavenumber, is directly the inverse of the wavelength λ, which is more convenient for physical interpretation of the results. In the case of an homogeneous material the polarization vector becomes constant in space and the classical plane wave solution is retrieved. Since the displacement field in Equation (10) is complex valued, its real part should be computed in order to retrieve the physical solution. From its definition as a Bloch-wave, the displacement at a point r image of the r 0 ∈ T by a translation t ∈ R has the following expression: The physical meaning is that the displacement vector at two homologous points (two points r 1 and r 2 are said homologous if r 1 − r 2 ∈ R) only differs by a phase factor. Additionally, the Bloch-wave expression in Equation (10) has the interesting property to be also R -periodic. Indeed, due to its R-periodicity, U k (r 0 ) can be decomposed as a Fourier series, leading to the equivalent expression for u k (r 0 ) : whereŨ k+ξ stands for the Fourier coefficients of the series expansion. Using this particular form, the R -periodicity of the Bloch-waves is easily proven using the change of variableξ = ξ + ξ : The main consequence of this property is that the characterization of the elastodynamics behavior of a periodic material does not require to investigate the mechanical response to all the k ∈ E but can be restricted to the study of k ∈ T via R -periodicity of the wavevector, where T corresponds to the First Brillouin Zone (FBZ), as introduced in Section 2.3.

Polarization of Waves in Homogeneous Materials
Before presenting the results, it is useful to recall some definitions concerning the polarization of elastic waves in homogeneous materials. Let us go back to the Bloch-wave ansatz introduced in Equation (10), since the material is now considered homogeneous, the polarization vector U k is constant both in space and time. In the most general case, the complex polarization vector U k , that will be denoted U from now on for the sake of simplicity, can be decomposed in its real and imaginary parts as follows: An interpretation of this decomposition, and of its consequences on wave propagation, can be obtained by considering the real part of Equation (10) : Since the vectors U R and U C are independent, the polarization of the displacement can be very rich. Its precise nature is directly related to conditions on U R and U C , as summarized in Table 2. It is important to remark that different conventions are used to define the handedness of circularly polarized waves. In this paper, we will consider that a wave is right handed if it follows the curl of the fingers of a right hand whose thumb is directed towards the wave propagation, away from the source. In the table, the unit normal vector defining the direction of propagation is defined by n = k k . Table 2. The polarizations of plane waves and conditions on the complex amplitude.

Polarization Condition
Longitudinal polarization U = αn with α ∈ C Transverse polarization U R and U C belong to the plane orthogonal to n Linear polarization A complex polarization vector can lead to a phase shift between the components of the displacement vectorû, and thus to a polarization that is other than linear, see Figure 5 for illustration of a linear and two circular polarizations with opposite handedness.

Dispersion Analysis Using Finite Elements Analysis (FEA)
In order to investigate the ultrasonic properties of the gyroid lattice, and given the periodicity of the architecture as described in the previous sections, an approach based on Bloch-Floquet analysis will be followed. For the sake of simplicity, the conventional unit cell depicted in Figure 6 is used to define the numerical model. For the investigation of the elastodynamic properties of the gyroid crystal, the wavevector will be restricted to the boundaries of the Irreducible Brillouin Zone (IBZ), as depicted in red in Figure 7a. The high symmetry points of this IBZ are defined in Table 1. The model has been implemented using the commercial software Comsol Multiphysics and considering titanium as constitutive material, the parameters of which are displayed in Table 3. The mesh of the unit cell is presented in Figure 6 and consists of 66,232 tetrahedral elements. Lagrange quadratic elements are used, for a total of 329,277 degrees of freedom. Periodic Bloch-Floquet conditions are implemented by imposing them as displacement conditions at the boundaries, following Equation (11). Then, the wavenumber in k is imposed and the corresponding frequencies are retrieved by solving the eigenvalue problem. The computation of each wavenumber took an average of 109 seconds on a workstation equipped with an Intel(R) Xeon(R) CPU E5-1650 v2 at 3.50 GHz using six cores.  The results of the dispersion analysis are depicted in Figure 7a. It can be observed, qualitatively, that these results are similar to those obtained for electromagnetic waves in [6] (see Figure 8). In particular, the behavior of the acoustic branches, i.e., those branches starting from the origin Γ, corresponding to transverse waves (gray lines in Figure 7a) is remarkably similar.  Since the objective of the paper is to investigate the behavior of the lattice within the LW-LF approximation given by classic continuum mechanics, the phase velocities and polarization of waves have been computed for large values of the wavelength with respect to the size of the unit cell. In this case, k = 16.7 m −1 , which corresponds to a wavelength close to 60 times the size of the unit cell. For each mode, the polarization vector has been estimated by computing the average of the complex displacement of the eigen-mode over the unit cell, and the results are listed in Table 4. We will now analyze the results along the following directions of propagation, also depicted in Figure 7b:

•
[001]: This direction is going from the center of the fundamental cell to the middle of a face. It corresponds to an axis of rotation of order 4 (rotations of π/2 rad); • [011]: This direction is going from the center of the fundamental cell to the middle of an edge. It corresponds to an axis of rotation of order 2 (rotations of π rad); • [111]: This direction is going from the center of the fundamental cell to a vertex. It corresponds to an axis of rotation of order 3 (rotations of 2π/3 rad).
Using the conditions listed in Table 2, we have characterized the polarization for each of the above propagation directions. The results are summarized in Table 4.
Along the direction [001] a longitudinal wave propagating at 3018.5 m/s can be observed. The two transverse waves have eigenmodes with complex amplitude and propagate with close phase velocity that tends to a common value of of 1932.2 m/s for infinite wavelengths, as it can be also observed in Figure 9. These complex amplitudes correspond to two circularly polarized transverse waves with opposite handedness. In the direction In summary, circularly polarized waves exist only if the direction of propagation is along a rotation axis of symmetry of order greater than 2. Moreover, as can be seen in Figure 9, for both [001] and [111] directions, the circularly polarized waves with opposite handedness propagate with the same phase velocity only in the infinite wavelength limit, and they start to diverge as the wavenumber increases. In particular, for direction [001] the right handed wave becomes slower than the left handed one, while the opposite phenomenon can be observed for direction [111]. This is due to the chirality of the unit cell. As one can notice, phase velocities are different even for a very large wavelength compared to the size of the unit cell, i.e., λ/a ∼ 10. Since only the phase velocity is affected, and not the amplitude, this effect can be interpreted as the elastic equivalent of circular birefringence in optics. This means that if a linearly polarized wave passes through a gyroid lattice, the polarization plane of the incident wave will be rotated. This is due to the phase difference (retardance) between the two circular components, which produces a rotation of the polarization plane. The concept is illustrated in Figure 10. Moreover, since phase velocity is involved in reflection of waves at boundaries via the Snell-Descartes law, and in particular in the definition of the Brewster angle of total reflection, gyroid lattices show the potential for being used as elastic circular polarizing filters.
Furthermore, the overall dispersion is normal for direction [001], i.e., phase velocity decreases when increasing the frequency, and anomalous for direction [001] (see [48] for the definition). The anisotropy of the material and the dispersive properties could also have consequences on surface and guided waves propagating in presence of boundaries [49,50], as well as in reflection/transmission problems [51]. These effects will be investigated in further works.

Long-Wavelength and Low-Frequency Approximation and Classic Elasticity
In this last section we will introduce and identify the equivalent homogenized model in the framework of classic linear elasticity. This equivalent homogenized model is characterized by a couple of effective tensors ρ and C H in such a way that the displacement field v is solution of the following problem: div C H : (v(r) ⊗ ∇) = ρv(r) (12) where v verifies < u >= v, < . > denotes the spatial average operator over T and u is the displacement field solution to the heterogeneous problem Equation (9), as done for instance in [10].
Since the effective continuum is homogeneous, we consider a plane wave solution with k = f /c n where c is the phase velocity of the wave and n a unitary vector. The substitution of this wave solution and of a linear elastic constitutive law into Equation (12) leads to following equation: where Γ = n · C H · n is the Christoffel, or acoustic, tensor. The solution of the eigenvalue problem stated in Equation (13) for a given direction n gives the phase velocities and polarizations of waves in the effective continuum.
In classical elasticity, a material with cubic symmetry is defined by three independent material constants. Using Mandel notation [52], the corresponding elastic tensor for a material having its symmetry axis parallel to B reads: It is worth noting that classical elasticity is insensitive to the lack of centrosymmetry [53]. The symmetry class of the elasticity tensor in the cubic system is O ⊕ Z c 2 meaning that even if the material symmetry group of the unit cell does not contain the inversion, the symmetry group of elasticity tensor will inherit it.
In the case of cubic materials, the solutions of Equation (13) listed in Table 5 are directly obtained. Using the phase velocities computed from the Bloch-Floquet analysis, parameters c 11 , c 12 , c 44 , and ρ can be identified and the homogeneous equivalent properties listed in Table 6 are then deduced. It can be noticed that, as presented in Table 5, since we considered the propagation along the rotational axes of symmetry, for each direction we observe a purely longitudinal wave and two purely transverse waves. We now move on to comparing phase velocities and polarizations obtained from the Bloch-Floquet analysis with the ones forecast by the Long-Wavelength and Low-Frequency approximation using classical elasticity. We start with direction [001]. As already mentioned, this direction corresponds to a rotational axis of symmetry of order 4. In this case, as the elasticity tensor is non sensitive to chirality, the symmetry group of the physical phenomenon (the symmetry group of the physical phenomenon is the intersection of the symmetry group of the constitutive equations and the symmetry group of the mechanical solicitation) is conjugate to D 4 ⊕ Z c 2 . Indeed, as the acoustic tensor defined in Equation (13) is a second-order tensor, the Hermann theorem of Crystal physics [11] predicts its behavior as transversely isotropic, i.e., O(2) ⊕ Z c 2 . As a consequence Γ must have an eigenspace of dimension 2. All the directions of wave propagation for which this is verified are called acoustic axis of the material. Moving to the results presented in Table 5, we can see that the classic theory indeed predicts one faster longitudinal wave and two slower transverse waves propagating with the same phase velocity.
In the previous section we saw that Bloch-Floquet analysis identifies these waves as circularly polarized with opposite handedness, which is of course equivalent. Indeed, even for very large wavelengths, the numeric evaluation of the polarization provided by Table 4 corresponds to the eigensystem (eigenvalues plus eigenvectors) of a Hermitian acoustic tensor. In such case, the space corresponding to the double eigenvalue (which is real due to Hermitian nature of the matrix) is two dimensional over the complex field C. However, in the case of a classic Cauchy continuum the acoustic tensor is symmetric real, and the eigen-space corresponding to the double eigenvalue is two dimensional over the real field R. Since the 2D vector space over C can be considered as a four dimensional vector space over R the span is not equivalent. Moreover, as presented in Secton 4, when the ratio between the wavelength and the size of the unit cell becomes finite, this 2D eigenspace splits into two 1D eigenspaces with different phase velocities. It is important to notice again that this effect occurs in the LF-LW regime, where the classical elastodynamic homogenization is supposed to hold, or give at least approximated results while preserving the physics of the problem. Similar results are obtained for propagation along [111], which corresponds to the rotational axis of symmetry of order 3. In this case the symmetry group of the physical phenomenon is conjugate to D 3 ⊕ Z c 2 , and thus again transverse isotropy is imposed to the acoustic tensor. Finally, the direction of propagation [011] is along to a rotation axis of symmetry of order 2, the physical point group is thus conjugate to D 2 ⊕ Z c 2 . Here, the symmetry class of the acoustic tensor is D 2 ⊕ Z c 2 , and all the eigenspaces are unidimensional. In this last case, as this kind of symmetry can be seen by second order tensors, the results from FEA on the heterogeneous material and classic elasticity are in agreement in the LF-LW regime.
In this section we have shown that some discrepancies can be observed when using an overall homogeneous continuum of Cauchy type. Classical elasticity (as opposed to generalized elasticity) is not rich enough to capture certain specific physical phenomena related to the symmetries of the material. In particular, if phase velocities are correctly estimated the polarizations are incorrectly predicted. Moreover, as it is well known, the onset of dispersion when frequency or wavenumber increase cannot be described in the classic Cauchy model.

Conclusions
In this work it was shown that a classical continuum model could not capture the correct behavior of elastic waves propagating in gyroid lattices. This is due to the fact that the classic continuum mechanics was insensitive to the lack of centrosymmetry of the architectured material. However, it is a well-established belief that the effects of non-centrosymmetry are only related to waves having a wavelength which compares to the size of the microstructure. Here we demonstrate that the solution given by the classical theory failed to predict the correct response, even in the Long Wavelength-Low Frequency domain.
In order to capture the onset of this unconventional behavior, called acoustical activity, the elastic continuum model needs to be enriched. Different strategies of enrichment are possible. In particular, the use of strain-gradient elasticity will be investigate in a forthcoming study.
The main practical consequence of the results presented in this work was in the evidence that circularly polarized waves could be observed in gyroid lattices, and that classic models failed to describe such effect. This could have a practical interest, since devices based on manipulation of circular polarization are frequently used in optics and electromagnetism. However, if one wants to exploit the same effects in mechanics, it appears important not to rely on classical theories of elasticity. Finally, it should be noted that, in this work, we addressed bulk wave propagation in an infinite medium. The interaction of these waves with boundaries, in the case of reflection/transmission problems or in the case of guided propagation, also deserves to be investigated.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Dictionary
To obtain the normal forms for the different classes the generators provided in the following table have been used :  Where the following elements of O(3) will be used in this study: (3) the reflection through the plane normal to n (P n = 1 − 2n ⊗ n). Table A2. Dictionary between different group notations for Type I subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.  Table A3. Dictionary between different group notations for Type II subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.

Hermann-Maugin Schonflies
Group Nature Type III Subgroups Table A4. Dictionary between different group notations for Type III subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.

Appendix B. Generators of Space Group #214
Consider the affine space E 3 , the vector space R 3 acts on E 3 by translations. The affine group Aff(E 3 ) of E 3 , which is the set of all affine invertible transformations is constructed as the semidirect product of R 3 by GL (3), the general linear group of R 3 : as such, an affine transformation is given by a pair (Q, v) ∈ GL(R 3 ) × R 3 . Composition of transformations follows from the construction of Aff(E 3 ) as a semi-direct product, to be explicit:

Elements of Aff(E 3 ) can be nicely represented by(4x4) block matrices:
Q v 0 1 the internal law in Aff(E 3 ) following the matrix product in M 4,4 .
For our needs, we are interested not in the full affine group but in the group of isometries of E 3 , this group Euc(E 3 ) is a subgroup of Aff(E 3 ) and defined as the semi direct product of the orthogonal group and the spatial translation of R 3 : Space groups can be considered as discrete subgroups of Euc(E 3 ). The generators of the space group I4 1 32 (No. 214) are given in the following table in various notations [54]:

Appendix C. Proof
The gyroid lattice is defined from an implicit equation (Equation (1)) that creates a periodic surface. For a given value of parameter b (= √ 2), this surface is found to become singular, thus creating an unrealistic discontinuous solid. This section presents an explanation for the admissible variation range of gyroid parameter :|b| < √ 2. Let us first restrict the variation range of variables x, y, z in Equation (2) to [0, 1 /2] in order to work in the fundamental domain of function φ. The gyroid lattice restricted to this domain is presented in Figure A1a). The fundamental domain of the gyroid is invariant with respect to the following symmetry operations: -Rotation of angle 2π /3 along the axis defined by equations y = z = x, plotted in plain line in Figure A1 and corresponding to the transformation (x, y, z) → (y, z, x). The directing vector of this axis is (1, It is trivial to see that these transformations leave Equation (2) unchanged thus defining symmetry operations. As a consequence, the point symmetry group of the fundamental domain of the gyroid is conjugated to D 3 . For the sake of simplicity, we will only consider generating operations of 2π /3 rotation about (x, x, x) axis and π rotation about (x, 1/4 − x, 1/8) axis, denoted C 3 and C 2 in the following, respectively; the two other π rotations being generated by combination of these two generators.
If the gyroid surface intersects one of the rotation axes non-orthogonally, then the surface automatically becomes degenerate. The expression of the normal director to the gyroid surface at its intersection point with generating symmetry axes and the equation defining this intersection point are summarized in the following Table A6: Table A6. The expression of the normal director to the gyroid surface at the intersection point with its symmetry axes and expression of this intersection point.

Sym. axis & Directing Vector
Normal to the Gyroid Surface Intersection Point Note that the normal director to the gyroid surface depends on variable x which, itself, is determined by parameter b through the non-linear equation defining the intersection point at which the normal director is computed.
One can easily check that the normal directors are generically colinear with directing vectors of C 3 and C 2 operations. However, for given values of variable x (or equivalently of parameter b), the normal to the gyroid surface is null, thus leading to singularity of the gyroid surface. These values are x = 0-and thus b = √ 2-leading to singularity of the gyroid surface at its intersection with axis C 2 and x = 1 /8-and thus b = 3 /2-leading to singularity of the gyroid surface at its intersection with both axes.
Finally, the equation defining the intersection point between the gyroid surface and the C 2 axis (see Figure A1b) depends on x and parameter b : √ 2 cos x + sin x 2 − b = 0. By plotting this equation considering b as a function of x, we can see that there are two intersection points between the gyroid surface and the C 2 axis for values of b over √ 2 thus showing that the gyroid surface forms a closed domain in these directions leading to an unrealistic discontinuous solid.