Joseph Gazalet
email: joseph.gazalet@univ-valenciennes.fr

Samuel Dupont

Jean-Claude C Kastelik

Quentin Rolland

Bahram Djafari-Rouhani

A tutorial survey on waves propagating in periodic media: electronic, photonic and phononic crystals. Perception of Bloch theorem in both real and Fourier domains

Introduction

The advent of photonic and phononic crystals has enlarged the application fields of wave propagation in periodic media. Initially mostly limited to electronic wavefunctions and lattice vibrations i.e., solid state crystals [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF][START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF], periodic media include now studies concerning also optical and acoustical waves in man-made periodically nano-structured materials [START_REF] Yablonovitch | Inhibited spontaneous emission in solid-state physics and electronics[END_REF][START_REF] Djafari-Rouhani | Rayleigh waves on a superlattice stratified normal to the surface[END_REF][START_REF] Kushwaha | Theory of acoustic band structure of periodic elastic composites[END_REF][START_REF] Touratier | Floquet waves in a body with slender periodic structure[END_REF][START_REF] Lourtioz | Les cristaux photoniques ou la lumière en cage[END_REF][START_REF] Meade | Existence of a photonic band gap in two dimensions[END_REF][START_REF] Ouchani | Defect modes in one-dimensional anisotropic photonic crystal[END_REF][START_REF] Hassouani | Dual phononic and photonic band gaps in a periodic array of pillars deposited on a thin plate[END_REF][START_REF] Kushwaha | Acoustic band structure of periodic elastic composites[END_REF][START_REF] Rolland | Acousto-optic couplings in twodimensional phoxonic crystal cavities[END_REF][START_REF] Craster | Bloch dispersion and high frequency homogenization for separable doubly-periodic structures[END_REF][START_REF] Liu | Efficient omnidirectional couplers achieved by anisotropic photonic crystal waveguides[END_REF][START_REF] Vellendera | Eigenfrequency correction of Bloch-Floquet waves in a thin periodic bi-material strip with cracks lying on perfect and imperfect interfaces[END_REF][START_REF] Linton | Resonant effects in scattering by periodic arrays[END_REF][START_REF] Yariv | Equivalence of the coupled-mode and Floquet-Bloch formalisms in periodic optical waveguides[END_REF][START_REF] Magnusson | Equivalence of multiwave coupled-wave theory and modal theory for periodic-media diffraction[END_REF]. Common to all these topics is the fact that the propagating wave functions are governed by a second order partial differential equation with periodic coefficients. These periodic coefficients stand for the material parameters characterizing the crystal under consideration; for example, periodic potential for solid state crystals, electric permittivity and/or magnetic Abbreviations: RD-CFE, Real space Domain Closed-Form Expression of a periodic function; FD-CFE, Fourier Domain Closed-Form Expression of a periodic function.

1 permeability for photonic crystals, mass density, and stiffness coefficient for phononic crystals. The general form of the eigensolutions, usually called Bloch functions or Bloch modes, is dictated by the Bloch theorem. Explicitly, this theorem states that only a special type of waveforms may propagate in periodic media. This result is a direct consequence of the discrete translation symmetry of the crystal structure. It is a universal feature in the sense that it does not depend on the nature of the propagating wave: electronic, photonic or phononic, nor on the particular distribution of the periodic parameters inside a period. Traditionally, the Bloch theorem is expressed by one of the two equivalent statements: any eigenfunction may be expressed in the form of a propagating plane wave modulated by a function which exhibits the periodicity of the crystal lattice. This corresponds to the Floquet statement, restricted to purely propagating wave (i.e. not subject to attenuation): when the wave function undergoes a translation by any direct lattice translation vector R, it differs from the original function by at most a phase factor e ik•R . It is worth noting that none of these statements completely defines the wave function but only gives its general form. Indeed, the periodic modulating function depends on the particular distribution of the periodic parameters characterizing the propagating medium.

Over time, different methods have been initiated to demonstrate the Bloch theorem. These methods may be classified in two main approaches according to whether they are developed in real space domain or within the space of Fourier domain. The later, according to the crystallographic terminology, is frequently called reciprocal domain as opposed to the direct domain which refers to the real three-dimensional space. In the real space approaches, we can distinguish between the classical method of solving differential equations [START_REF] Floquet | Sur les équations différentielles à coefficients périodiques[END_REF] from the subsequent methods [START_REF] Bloch | Über die Quantenmechanik der Elektronen in Kristallgittern[END_REF], inherited from quantum mechanics which are based on group theory and/or theory of operators [START_REF] Laloë | Mécanique quantique[END_REF][START_REF] Rosencher | [END_REF][START_REF] Bir | Symetry and Strain-Induced Effects in Semiconductors[END_REF]. Similarly, in the reciprocal domain, approaches in which the periodic parameters were described by a single sinusoidal function [START_REF] Mathieu | Le mouvement vibratoire d'une membrane de forma elliptique[END_REF] or by their Fourier series expansion [START_REF] Hill | On The Lunar Perigee Which Is A Function Of The Lean Motion Of The Sun And Moon[END_REF], may be viewed as the classical precursor methods for later derivations which takes full advantage of the uninterrupted developments in Fourier analysis techniques [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF].

It is advisable not to oppose these two methods but rather to take advantage from their complementarities. Indeed, understanding the two domain correspondence is a sine qua non condition in order to have a firm grasp about the wave propagation in periodic media. The two domain representations of wave propagation in periodic media are intricately linked via Fourier transform. Brillouin had exploited this two domain correspondence associated to its solid experience on the electronic filters to develop the concept of band structure i.e. ''the Fourier domain characteristics'' of wave propagation in periodic media [START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF]. The concept of band structure includes: the establishment of the periodicity in the reciprocal domain, the definition of the Brillouin zones and the understanding of the existence of allowed and forbidden (bandgap) frequency bands as well as the slow wave principle.

In this paper following the Brillouin conception of the two domain correspondence, we give a Bloch theorem interpretation based only on Fourier transform considerations. This has been made possible since the approach is constructed directly in the distribution function framework both in the direct and Fourier domains. Also, this approach appears to be very versatile and relevant whatever is the complexity of the medium (its dimensionality or anisotropy) or that of the field under consideration (scalar or vector). The method is so straightforward that, even in the most complex cases, the reciprocal lattice wave equation can be established by a simple inspection of its expression in the direct domain.

In Section 2, we give a tutorial overview on the classical interpretations of wave propagation in periodic media. Real and Fourier domain methods are separated in two independent subsections. Wherever needs arise, basic concepts such as: eigensolutions and their k-labeling, reciprocal lattice, Brillouin zones, etc., are introduced and explained. In order to lighten this overview and at the same time preserve its completeness, all the nonessential notions in a first reading are deferred in appendixes. Section 3 is devoted to illustrate the Brillouin vision of the problem i.e. the indispensable perception of the correspondence between the two domain representations in order to draw meaningful conclusions. In Section 4, following the Brillouin illustrative example of elastic wave propagation in regularly spaced point masses, we present a versatile formalism well adapted to investigate the two domains correspondence in the field of wave propagation in periodic media. This formalism is inspired from digital electronic filter and/or discrete signal analysis techniques which use well established mathematical tools originated from the distribution functions theory. In Section 5, we illustrate the application of this approach to different situations. We begin by the simplest case of scalar waves (de Broglie waves) in solid state crystals. Then, we consider vector fields in such media that does not exhibit any other anisotropy than that introduced by the periodic non-homogeneities illustrated by electromagnetic fields in photonic crystals. Finally we address phononic crystals as an illustration of the general case of vector field in periodic media where the constitutive materials exhibit anisotropic as well as piezoelectric properties. Finally, in Section 6 we present a comprehensive discussion of the physical aspects of the presented approach.

Overview on the classical demonstrations of the Bloch theorem

The aim of this section is to present an overview on the standard methods used to introduce the Bloch nature of propagating waves in periodic media. Special attention is devoted on their basic concepts and physical interpretations rather than on mathematical rigor. Readers may wish to refer to the provided literature for more details.

To introduce these classical methods, we limit ourselves to the simplest case of scalar waves i.e. electronic waves in a solid state crystal. The wave equation governing the electronic wave is the Schrödinger equation. Under the assumptions of harmonic dependence ψ (x, t) = ψ (x) e -i E h t , the one-dimensional time independent Schrödinger equation reduces to the well-known form [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF]:

- h2 2m ∂ 2 ψ (x) ∂x 2 + U (x) ψ (x) = Eψ (x) (1)
using standard notations: U (x) for the periodic potential of period a, E = hω for the eigenenergy, h for the reduced Planck's constant, ω for the angular frequency and m for the mass of particle. The correspondent three-dimensional equation, reads:

-h2 2m ∇ 2 ψ (r) + U (r) ψ (r) = Eψ (r) .

(

) 2
The operator ∇ 2 stands for the Laplacian operator and the three-dimension periodicity of the potential U (r) is depicted using the concept of a primitive 3-D unit cell i.e., the smallest region of the space which replicate itself periodically. We say that each unit cell is assigned to a node of a virtual lattice known as the Bravais or direct lattice defined by three basis vectors a 1 , a 2 and a 3 . Our aim is to find the general form of the solutions of partial differential equations with periodic coefficients.

We separate the techniques used to derive the Bloch theorem in two distinct subsections i.e. the real space domain and the Fourier domain approaches. Each method reveals some aspect of the physical problem. So, to have a comprehensive perception of the problem, it would be interesting and insightful to have a good idea on both methods. These two subsections are constructed independently and may be read in any sequence. In order to lighten this overview, all the nonessential notions in a first reading are deferred in appendixes.

Real space methods

We begin this subsection by introducing the classical method used by Floquet for the investigation of a differential equation with periodic coefficients. This will enable one to familiarize with the main concepts and have a better understanding of the more elaborated Bloch method based on group theory.

Classical approach of Floquet

The paper of Floquet published in 1883 [START_REF] Floquet | Sur les équations différentielles à coefficients périodiques[END_REF] does not directly concern wave propagation but it falls under the classical theory of differential equations. This work remains very instructive, since he carried out a complete investigation of the one dimensional linear homogeneous differential equations of order n with periodic coefficients, explicitly this equation writes:

d n ψ (x) dx n + p 1 (x) d n-1 ψ (x) dx n-1 + p 2 (x) d n-2 ψ (x) dx n-2 + • • • p n (x)ψ (x) = 0 (3)
where, the p n stand for periodic coefficients, all of them exhibiting the same period a. i.e. p n (x + a) = p n (x).

In a first step Floquet observed that if ψ i (x) is a given solution then ψ i (x + a) will also be a solution. Indeed, the periodic coefficients are invariant under the change of variables x  → (x + a) i.e. p n (x + a) = p n (x) whereas, ψ i (x + a) ̸ = ψ i (x). He concluded that, if a complete set {ψ i (x)} of independent solutions ψ i (x); i = 1, . . . , n have been established, then the respective expressions of any of the ψ i (x + a) may be written as linear combinations of this set:

ψ i (x + a) = n  j=1 A ij ψ j (x) (4)
where, the matrix elements A ij represent the weighting factor of the ψ j (x) components in the expression ψ i (x + a) i.e. image of the ith basis function ψ i (x) after the change of variables x  → (x + a). In the second step, he seeks for the general form of the solutions χ (x) which he baptized periodic function of 2nd kind. That is, a function satisfying the condition:

χ (x + a) = εχ (x), (5)
where, ε is a complex constant. Again, he expanded this general solution in the basis {ψ i (x)}.

χ (x) = n  p=1 c p ψ p (x) . (6)
Combining these three steps, and identifying the coefficients of ψ p (x) he ended in an eigenproblem in terms of the weighting constants of the form:

n  l=1 c l A lp = εc p . (7)
Then, Floquet undergoes a complete mathematical study (which extended about forty pages) to establish the existence and the uniqueness of this general form of the solution. He ended with the theorem which states that: At least one of the solutions of the differential equation with the general form of χ (x + a) = εχ (x) exists for each distinct eigenvalue ε of Eq. [START_REF] Touratier | Floquet waves in a body with slender periodic structure[END_REF]. In other words, this solution is unique if its eigenvalue corresponds to a single-root of the eigenproblem. Whereas, it exists but it is not unique in case of multiple roots eigenvalues. Some important physical aspects are implicitly associated to the general form of the Floquet solution given in Eq. ( 5) and are introduced in the next subsection.

Physical contents of the Floquet theorem

Eigennature of the Floquet theorem.

Let us draw the attention on the ''eigennature'' of the Floquet theorem. Looking beyond the mathematical solution, it is instructive to realize the implications of the eigenproblem.

On the one-hand, we see that the expression of ψ i (x + a) in Eq. ( 4) contains not only the original function ψ i (x) but also contributions from all other independent solutions ψ j (x). This obviously leads to dissimilar overall shape for ψ i (x) and its image ψ i (x + a). On the other hand, it is well known that the choice of the set of independent solutions {ψ i (x)} is not unique. So, by introducing the condition given in Eq. ( 5) Floquet imposes to choose the specific set of independent solutions in such a manner that their overall shape is not altered as they undergo a discrete shifting. Indeed according to this condition χ (x +a) differs from χ (x) by at most a complex factor ε.

The fact that the function remains globally unchanged when it is subjected to the change of variable x  → (x + a) explains the German label eigen which means: own, inherent, appropriate . . . . For the same reason, some authors employ instead of eigen the terms normal or proper to qualify this type of modes because they represent the modes which will ''propagate'' without distortion when injected in the media.

Bloch-Floquet theorem.

Concerning the problem of wave propagation in periodic media, we are interested with the special case where |ε| = 1 (i.e. neither amplification/attenuation phenomena are considered). In other words the complex eigenvalue takes the form of a phase factor ε = e iθ ; -π < θ ≤ π Modulo (2π ). Without loss of generality, introducing the change of variable k = θ/a, we can write:

ε = e ika ; - π a < k ≤ π a Modulo  2π a  . (8)
Finally, substituting in Eq. ( 5), the solution wavefunction has to satisfy the property:

χ k (x + a) = e ika χ k (x). (9)
This is equivalent to say that the general form of the solution is:

χ k (x) = e ikx µ k (x) (10)
where, µ(x) is a periodic function. Any of these two equivalent statements characterizing the property of a propagating wave in a periodic media are known as the Bloch-Floquet Theorem.

k-labeling of the modes.

As it is well known, expressing Eq. [START_REF] Touratier | Floquet waves in a body with slender periodic structure[END_REF] in the specific set of eigenvectors ( 5) is equivalent to matrix diagonalization. Reciprocally, the solution of the eigenproblem leads to the complete set of eigensolutions {χ i (x)} and each solution corresponds to an eigenvalue ε. Finally, the later being uniquely defined by the principal value of k via Eq. ( 8) explains the k-labeling of the modes. This observation is a special case of the more general concept of mode classification according to symmetry operations as it will be mentioned in the following section.

The ''raison d'être'' of the variable k.

In this section, we have introduced the variable k as a normalization of the phase shift with respect to the period a. But, this variable known as the wavevector owns implicit physical meanings that will be revealed as we go along in this tutorial. For the moment, let us highlight that if k (equivalently θ ) is complex or imaginary, the wave is damped or decreases exponentially. It does not represent a propagating wave. We say that the corresponding angular frequency ω belongs to a ''forbidden band'' also termed ''bandgap''.

Principle of the 1st Brillouin zone.

The reason of the restriction of the k value to the interval -π /a < k ≤ π /a (i.e. to its principal value) is obvious. Indeed, according to Eq. ( 8) the eigenvalue ε is a periodic function of k with period 2π /a. So, for a given k value inside this interval, any k ′ value such that k ′ = k ± 2nπ /a will correspond to the same eigenvalue ε and subsequently the same eigenfunction χ k ′ (x) = χ k (x). Thus, all distinct solutions are completely defined on the principal value of k. This interval -π /a < k ≤ π /a is the well-known 1st Brillouin zone. The solutions replicate themselves periodically outside this interval i.e. in the higher order Brillouin zones.

Clearly, there is a degree of arbitrariness in centering the 1st Brillouin zone on the k = 0. In fact, mathematically speaking, any interval of period 2π /a will equivalently cover all the possible eigensolutions. Brillouin justified the practical interest of this convention for wave propagation. The fundamental reason of this choice is to insure the even parity of the band diagram. We will return to this point later in Section 2.2 and the fundamental reason will become evident in Section 3.2.

Viewpoint of Bloch

In his paper of 1928 concerning quantum mechanics of the electrons in crystal lattices, Bloch [START_REF] Bloch | Über die Quantenmechanik der Elektronen in Kristallgittern[END_REF] seeks the solutions for the Schrödinger equation in the presence of three-dimensional periodic potential. This problem corresponds to an extension of the Floquet differential equation to the three dimensions but restricted to a second order differential equation.

In the derivation of the theorem, he makes use of the suggestion of Wigner [START_REF] Wigner | Einige Folgerungen aus der Schrödingerschen Theorie für die Termstrukturen[END_REF], to apply well established lemma of the group theory to quantum mechanics. This constitutes an elegant way which enables shortcuts as compared to the classical differential equations theory used by Floquet [START_REF] Floquet | Sur les équations différentielles à coefficients périodiques[END_REF]. Readers non familiar with group theory may, among the abundant available literature, begin with handy introduction text books for example: [START_REF] Kireev | Semiconductor Physics[END_REF] in the field of the solid state, [START_REF] Hümmer | Group Theory in Semiconductor Optics[END_REF] in semiconductor optics and [START_REF] Joannopoulos | Photonic Crystals. Molding the Low of Light[END_REF] within the framework of photonic crystals.

What we only need to know here, among some dedicated terminology, is the definition of the ''representation'' of a group. So, rather than going into the details of the group theory, it may be worth having, in this tutorial, an overview of the method used by Bloch himself [START_REF] Bloch | Über die Quantenmechanik der Elektronen in Kristallgittern[END_REF]. Indeed, at this time the theory of groups had just been introduced in the field of quantum mechanics. So, on the one-hand, Bloch relies on physical foundations to introduce the use abstract groups and on the other-hand he gives a brief demonstration of steps which are generally eluded when using the group theory shortcuts.

2.1.3.1. Some Terminology of group theory. Bloch begins with the same observation pointed out by Floquet i.e. the invariance of the Schrödinger equation under the change of variables. He also expands the resulting functions in terms of basis functions. He makes use of the terminology of the operator theory and symmetry transformations. So, let us introduce some knowledge of the dedicated terminology taking translation symmetry as an illustrative example.

• The translation Group.

As stated above, classically, the periodic parameters of the crystal are described on a ''primitive unit cell'' and the periodicity is depicted by replicating this primitive cell to all the sites of the Bravais lattice described, according to its dimensionality, by one, two or three basis vectors (a i ; i = 1, 2, 3). Any lattice site is uniquely addressed by a lattice vector

R = n 1 a 1 + n 2 a 2 + n 3 a 3 ; n 1 , n 2 , n 3 ∈ Z .
The terminology used in group theory originates from the mathematical branch concerning symmetry. The periodicity of a crystal is described as the ''discrete translation symmetries''. That is, all the possible displacements of an unlimited crystal (as a whole) that make the crystal to coincide into itself. We say that ''the crystal is ''invariant'' after a symmetry transformation''. We also say that each translation, corresponding to a different lattice vector R = n 1 a 1 + n 2 a 2 + n 3 a 3 , i.e. to a dissimilar triplet of integers (n 1 , n 2 , n 3 ), constitutes a distinct ''symmetry element'' belonging to the ''set'' of discrete translation symmetries. We will denote a translation symmetry element as T R . Clearly, the elements of the discrete symmetry translations constitute an infinite set. This set is known as the ''translation group'' because it fulfills the four criteria defining an abstract mathematical group [START_REF] Kireev | Semiconductor Physics[END_REF][START_REF] Hümmer | Group Theory in Semiconductor Optics[END_REF][START_REF] Landau | Physique Theorique: Tome 3, Mécanique quantique[END_REF] (See Appendix A). For the moment, let us define the ''product of two ''translation elements''. It is the result obtained when we apply successively these two translations. Clearly, it is also a translation with a translation vector equals to the sum or the two translation vectors:

T R 1 T R 2 = T (R 1 +R 2 ) .
Because this product operation is commutative, we say that the translational group is ''Abelian''. This special type of groups owns important properties used in the determination of eigenvalues.

• The translation operators.

In the framework of operator theory, a space translation T R , relative to a given lattice vector R, is described by the ''translation operator'' TR . That is, an operator which when acting on a function f (r) return the value f (r + R):

TR f (r) ≡ f (r + R) . [START_REF] Hassouani | Dual phononic and photonic band gaps in a periodic array of pillars deposited on a thin plate[END_REF] Following this definition it is clear that the action of the translational operator on a function is equivalent to the change of variables r  → (r + R) pointed out in the Floquet method.

In general, f (r) differs from f (r + R), the equality holds only if the function exhibit the periodicity of the crystal or equivalently, its discrete translational symmetry. For instance, this is the case for the periodic potential U(r) = U(r + R).

We say that the periodic potential is ''invariant'' under these translational operations.

The short cuts introduced by the group theory.

Using the operator notation Schrödinger equation ( 2) is written in the form of an eigenproblem:

Ĥψ E (r) = Eψ E (r) (12)
where, Ĥ stands for the Hamiltonian operator:

Ĥ ≡  -  h2 /2m  ∇ 2 + U (x)
 and the subscript E recalls that a specific solution corresponds to a given eigenvalue. Here, the value of the eigenenergy E.

So, it is straightforward to show that, the invariance of the Hamiltonian under the change of variables r  → (r + R), mentioned by Floquet (see Section 2.1.1) is equivalent to the commutation property of the translation operators with the Hamiltonian (Appendix B). A direct corollary is that the commutation property of the translation operators with the Hamiltonian implies that any eigensolution of the translational operator is also an eigensolution of the Schrödinger equation with the same eigenvalue E (Appendix C). While in this tutorial we limit ourselves to the translation group, these two statements are not restricted to translation symmetry, but still hold for all symmetry operations including point symmetries. In group theory the space group is defined as the generalized group including translational and point symmetry groups.

The fact that translation and Hamiltonian operators have a common set of eigenvectors represents the essential benefit of the theory of operator associated to Group theory. The properties of translational groups are well known. Specifically, the general form of the eigensolutions of the translation group is known [START_REF] Kireev | Semiconductor Physics[END_REF][START_REF] Hümmer | Group Theory in Semiconductor Optics[END_REF][START_REF] Joannopoulos | Photonic Crystals. Molding the Low of Light[END_REF][START_REF] Landau | Physique Theorique: Tome 3, Mécanique quantique[END_REF]. For instance, thanks to the above mentioned Abelian nature of the translational groups it is easy to show that the eigenvalues of the translational operators are: e ik•R . This implies that the eigenfunctions of the translational operator verify: TR ψ k (r) ≡ ψ k (r + R) = e ik•R ψ k (r). This later is an expression of the Bloch theorem. Thus, the demonstration is done.

So, this is an elegant straightforward method which relegates to the well-known group theory all the mathematical determinations of the solutions as well as the associated investigations of the existence and/or uniqueness of the solution. It also allows one to clearly dissociate the commonly used double subscripts E and k to label the solution wavefunctions ψ E,k (r) of the Schrödinger equation. The former label is related to the eigenvalues of the Hamiltonian operator while the later is referred to the eigenvalues of the symmetry translation operator.

However, this method has the drawbacks that go with its qualities: anyone who is not initiated to this theory would feel hard done by, not being able to grasp the ins and outs of the physical content. Indeed, all the basic concepts are intrinsically linked to group theory and are inaccessible without a brief presentation of this theory.

Roughly speaking, the approach followed by Bloch to introduce the group theory is a generalization of the Floquet method within a well structured elegant formalism. The major advantage of introducing the group theory remains to transfer the problem to already resolved problems in the field of symmetry operations. Indeed, the crystal symmetry operations are limited to few ones and have been extensively studied and the results tabulated. The starting point resides in the fact that the eigenfunctions of the symmetry operators are also eigensolutions of the differential equation. For completeness we present in the Appendix D the details of the approach followed by Bloch to introduce the group theory. It may be viewed as an illustrative example that provides more details on the group theory techniques.

Foundation of the three-dimensional reciprocal lattice.

In the framework of the virtual box defined as the restriction of the unlimited 3-D crystal to a subsystem of parallelepiped shape sustained by (N 1 ×N 2 ×N 3 ) unit cells disposed respectively along the three directions a 1 , a 2 and a 3 (N 1 , N 2 and N 3 being very large integers), we shown in Appendix D, that a given eigensolution ϕ ′ j of the translation operator has the general form:

ϕ ′ j (x + n 1 a 1 , y + n 2 a 2 , z + n 3 a 3 ) = e i  2π m 1,j N 1 n 1 + 2π m 2,j N 2 n 2 + 2π m 3,j N 3 n 3  ϕ ′ j (x, y, z); ( 13
)
where the index j labels the triplet of integers:

 m 1,j , m 2,j , m 3,j
 which define the mode under consideration.

Eq. ( 13) is nothing but the 3-D extension of the Bloch-Floquet theorem introduced in Section 2.1.2.2. To perceive this assertion, it is convenient to express Eq. ( 13) in a form analogue to Eq. ( 9). That is, using vector notations Eq. ( 13) may be written in the form:

ϕ ′ k (r + R) = e i(k•R) ϕ ′ k (r). (14)
This is made possible using the direct lattice vector R = n 1 a 1 + n 2 a 2 + n 3 a 3 and introducing the discrete wavevector k so that the dot product k • R verifies the exponential argument in Eq. (13). This imposes the following expression for the wavevector:

k = m 1 N 1 b 1 + m 2 N 2 b 2 + m 3 N 3 b 3 (15)
where

: b 1 = 2π a 2 × a 3 a 1 • (a 2 × a 3 ) ; b 2 = 2π a 3 × a 1 a 2 • (a 3 × a 1 ) and b 3 = 2π a 1 × a 2 a 3 • (a 1 × a 2 ) (16)
are known as the reciprocal lattice basis vectors 1 The k-labeling of the eigensolutions in Eq. ( 14) takes the place of the j subscript in Eq. [START_REF] Rolland | Acousto-optic couplings in twodimensional phoxonic crystal cavities[END_REF]. Note that the discreteness of k-values is due to the artificially introduced ''virtual box''. We say that k is quasi-continuous. Indeed, since the values of the N i are very large numbers, the step between two adjacent values of k is vanishingly small.

To close this section, remember that since the translational operator commutes with the Hamiltonian, the general form of the eigensolutions ϕ ′ k (r) of the translational operator given in Eq. ( 14) stands also for the eigensolutions ψ E,k of the Schrödinger equation. The indices E and k respectively recall that each distinct solution corresponds to an eigenenergy of the Hamiltonian and an eigenvalue of the translational operator. 2.1.3.4. Brillouin zone in the three-dimensional space. As in Section 2.1.2.5, the Brillouin convention to confine the complex exponential argument into its principal value ]-π , π] leads to the definition of the three-dimensional 1st Brillouin zone. It is easy to show that the boundaries of Brillouin zone consist of the perpendicular bisector planes of the segments connecting one reciprocal lattice sites chosen as the origin (k = 0) to its nearest neighbor [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF][START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF]. These segments correspond to all the basis vectors b i of the reciprocal lattice which are equivalent with respect to the crystal point symmetry operations. Indeed, extending the principle introduced in Section 2.1.2.5 to the three dimensions case i.e. confining the principal value of the phase shift associated to any of the equivalent basis vector translation Ta i in the interval ]-π , π] gives us [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF][START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF]:

ψ E,k (r + R) = e i(k•R) ψ E,k (r) (17)
-π < k • a i ≤ π . (18)
In order to determine the 1st Brillouin zone boundaries we take the strict equality of the right-hand and/or left-hand sides of Eq. (18). For instance, for i = 1 the right-hand side equality gives us k • a 1 = π . Then, expressing k in the basis of reciprocal lattice vectors {b 1 , b 2 , b 3 }, we can write:

k • a 1 ≡ (k 1 b 1 + k 2 b 2 + k 3 b 3 ) • a 1 = π . ( 19
)
Following the definition of the reciprocal lattice basis vectors b i • a j = 2π δ ij , Eq. ( 19) reduces to:

k 1 = 1/2. ( 20
)
This means that the component, along the direction of b 1 , of any k vector having its terminal point on the 1st Brillouin Zone Boundary, would be exactly equal to |b 1 /2|. In other words, the locus of the terminal points of these k-vectors is the plane that passes perpendicularly through the midpoint of b 1 . The later statement gives the geometrical procedure to construct the 3-D boundaries as shown in Fig. 1. While, the first statement gives the mathematically expression of the boundaries:

k • (b 1 /2) |b 1 /2| = |b 1 /2| ⇔ k • b 1 2 =     b 1 2     2 . ( 21
)
Referring to the equivalent perpendicular bisector plane ′ of the reciprocal vector -b 1 (instead of b 1 ) this equation is usually written in the form [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Kireev | Semiconductor Physics[END_REF]:

2k • b 1 + |b 1 | 2 = 0. ( 22
)
Doing the same for all other symmetrically equivalent directions, we obtain the other boundaries of the 1st Brillouin zone. The mathematical expressions of higher order Brillouin zone boundaries are easily obtain using the general reciprocal lattice vector G in place of b 1 in Eq. ( 22). This writes: (2k

+ G) • G = 0 [1].

Concluding remarks

Roughly speaking, we can say that Floquet and Bloch approaches are analogous. The major difference resides in the fact that Bloch elegantly eludes the classical exhaustive mathematical analysis undergone by Floquet. This is made possible because these mathematical problems are already solved in the framework of group theory. Besides giving the general form of the solution, the Bloch-Floquet theorem involves other important physical aspects of wave propagation in periodic media.

For instance: Reciprocal space (i.e. the need of introducing a wavevector k defined in a reciprocal lattice), first Brillouin zone (i.e. the loci of all the distinct principal values of the wavevector), higher Brillouin zones (i.e. the periodicity in the reciprocal domain) etc. A clearer perception of these topics will be acquired in the following subsection since all these quantities will acquire a Fourier domain significance. The correspondence between the direct and reciprocal domains is addressed in Sections 3 and 4.

Fourier domain methods

The second standard methods are developed directly in the Fourier Domain. These approaches are intimately linked to the Plane Wave Expansion (PWE) representation of wavefunctions. Indeed in both cases we represent a propagating wave by the superposition (i.e. linear combination) of a set of monochromatic sinusoidal waves i.e. their Fourier components. Each of these elementary waves can be expressed as:

ϕ E,k (r, t) = ϕ E,k e -i(ωt-k•r) , (23)
which is completely defined by four quantities. Namely: the complex amplitude ϕ E,k , the frequency (or period) ω = 2π /T , the wavevector k (or the wavelength λ and direction of propagation ℓ; k = (2π /λ) ℓ), and their relative phase shift represented by the argument of the complex amplitude Arg  ϕ E,k 

. As indicated above, in solid state physics, the angular frequency is related to the energy by the reduced Planck constant ω = E/h. This set of monochromatic waves constitutes the spectrum of the wavefunction. This spectrum may be continuous or discrete and/or unlimited or finite.

This method conveys an improved understanding of the physical meaning for the k-vector which we have introduced in the real domain method as a more or less mathematical trick. It is worth noting that our purpose is not to oppose these two methods but rather to reveal their complementarities. Indeed, understanding the two domain correspondence is a sine qua non condition in order to have a firm grasp about the wave propagation in periodic media. This correspondence is the subject of Sections 3 and 4. For the moment let us specify this intricate linking by drawing attention to the fact that Fourier domain methods can be viewed as a special case of the translation group representation in a specific basis. Namely, ''the Fourier basis'' [START_REF] Laloë | Mécanique quantique[END_REF].

In this section we illustrate the concept of Fourier domain method using the simple case of one dimensional wavefunction. The extension to the three-dimensional is deferred in Appendix E. Indeed, it is straightforward and do not carry any new physical aspect to our survey.

One-dimension method of solving the Schrödinger equation in a periodic potential

The classical method of solving differential equations in the Fourier domain implies that the solutions sought are expressed in the form of the inverse Fourier transform of their spectral components ϕ E (k). Unfortunately, the case of wave propagating in an infinite periodic media corresponds to a situation that violates the conditions of the existence of a Fourier transform. 2 To overcome this difficulty, standard methods [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF] apply the concept of the ''virtual box'' associated to the Born-von Karman periodic boundary conditions as introduced in Appendix D. For a 1-D crystal consisting of N unit cells this boundary condition writes [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF][START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF][START_REF] Bloch | Über die Quantenmechanik der Elektronen in Kristallgittern[END_REF]:

ψ E (x + Na) = ψ E (x). (24)
Then, taking advantage of this artificially introduced periodicity one expands the solution in terms of Fourier series instead of Fourier transform. So, following the Fourier method, we assume in a first step that the Fourier components ϕ E,m of all the harmonics exist and write the Fourier series:

ψ E (x) =  m ϕ E,m e im 2π Na x ( 25
)
where the replicating period is the length of the virtual box L = Na. Note that for shortness we have omitted the time dependent exponential. In the same way, the potential U(x) being periodic with period a, it is also expressed in term of a Fourier series 3 :

U (x) =  n U n e -inbx (26)
where, b stands for the reciprocal basis vector b = 2π /a and U n for the Fourier components of the periodic potential. The later are given by:

U n = 1 a  period dxU (x) e inbx ( 27
)
2 As we will see in Section 4 a Fourier transform exists but in the field of distribution functions.

3 In all cases, except for the solution of the Schrödinger equation we will use this sign convention. Note that in Eq. ( 25) we used the opposite convention. This is because the rule of thumb in time independent Schrödinger equation is to retain the exponential e -i(ωt-k•r) for the propagating wave traveling toward the positive k direction instead of e i(ωt-k•r) (see Appendix H for more details).

Each of the Fourier components in Eq. ( 25) is of the form of Eq. ( 23). We can identify a discrete form for the wavevector

k m = mk N .
Where, k N = b/N stands for the step of this discretization. Since N tends toward infinity this step is vanishingly small and so, k m is quasi-continuous as compared to the 1st Brillouin zone dimension b. This discretization is introduced by the artificial confinement of the wavefunction in the virtual box. For convenience let us rewrite Eq. ( 25) in term of k m :

ψ E (x) =  k m ϕ E,k m e ik m x (28)
2.2.1.1. The central equation and the dispersion relation. Now, in order to determine the Fourier components ψ E,m , we substitute Eqs. ( 26) and ( 28) in the Schrödinger equation ( 1). This gives us:

- h2 2m ∂ 2 ∂x 2  k m ϕ E,k m e ik m x +  n U n e -inbx  k m ϕ E,k m e ik m x = E  k m ϕ E,k m e ik m x . ( 29
)
Rearranging the terms we obtain:

 k m  h2 k 2 m 2m -E  ϕ E,k m e ik m x +  k ′ m  n U n ϕ E,k ′ m e i(k ′ m -nb)x = 0. (30)
Now, we change of dummy variables k ′ m to k m = k ′ mnb in the second summation in Eq. [START_REF] Joannopoulos | Photonic Crystals. Molding the Low of Light[END_REF]. Note that we have anticipated this change by adjoining a prime to this dummy variable in going from ( 29) to [START_REF] Joannopoulos | Photonic Crystals. Molding the Low of Light[END_REF]. So, Eq. ( 30) writes:

 k m  ϕ E,k m  h2 k 2 m 2m -E  +  n U n ϕ E,(k m +nb)  e ik m x = 0 (31)
This sum must be zero whatever is the value of x. This is only possible if each term of the sum vanishes independently. 4 Thus, in the Fourier domain, the expression of the Schrödinger equation reduces to an infinite set of linear equations involving the set of variables ϕ E,(k m +nb) ; -∞ < n < ∞. The lines of this system are of the form 5 :

ϕ E,k m  h2 k 2 m 2m -E  +  n U n ϕ E,(k m +nb) = 0. (32)
In the field of solid state physics this set of linear equation is known as the central equation [START_REF] Kittel | Introduction to Solid State Physics[END_REF]. It will lead to the band structure E(k) i.e. the expression of the eigenenergy as a function of the wavevector. Remembering that E = hω, this is analogue to the well-known dispersion relation [START_REF] Lourtioz | Les cristaux photoniques ou la lumière en cage[END_REF][START_REF] Joannopoulos | Photonic Crystals. Molding the Low of Light[END_REF][START_REF] Born | Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light[END_REF][START_REF] Royer | Elastic Waves in Solids. Part 1: Free and Guided Propagation[END_REF] in the field of propagating waves which expresses the angular frequency in terms of the wavevector. In what follows, we will omit the index m of k since the wavevector is quasicontinuous.

Bloch theorem.

The terms of the summation in the central equation [START_REF] Born | Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light[END_REF] reveal that the real space periodicity couple all those, and only those, Fourier spectral components which differ from one another by any reciprocal lattice vector G = nb.

Thus, once we are handling a given wave vector k, we are, in fact, dealing with a subsystem of linear equations involving all components having wave vectors of the form (k + nb) for all n belonging to the relative integers Z (excluding all other k values). This justifies the common ''k-labeling'' of eigensolutions ψ E,k (x) referring to the k-subsystem. So, for a fixed k value, only the terms belonging to k values of the form (k + nb) will be retained in Eq. ( 28). This later may be rewritten as:

ψ E,k (x) =  n ϕ E,(k+nb) e i(k+nb)x . ( 33
)
Or equivalently:

ψ E,k (x) =   n ϕ E,(k+nb) e inbx  e ikx . (34)
The term between brackets has the form of a Fourier series. So, it is a periodic function of x with period a (remember that b = 2π /a). Denoting this function µ k (x) Eq. ( 34) writes:

ψ E,k (x) = µ k (x) e ikx . ( 35
)
4 A more academic demonstration exists using the orthogonal property of the Fourier components [START_REF] Laloë | Mécanique quantique[END_REF][START_REF] Kireev | Semiconductor Physics[END_REF]. 5 Note that, according to Fourier transform properties, since U (x) is a real function, the Fourier components will exhibit the property

U -n = U * n and if in addition U (-x) = U (x), then U -n = U * n .
These properties are very helpful in solving Eq. ( 32).

This is one of the expressions the one dimensional Bloch theorem. Substituting x by x + a, we obtain the second expression:

ψ E,k (x + a) = e ika ψ E,k (x) . (36)
We recall that Bloch theorem gives only the general form of the solution. To completely determine the solution i.e. find its periodic part, we have to solve the linear system. The index k recalls that the Fourier components of µ k (x) are the solutions of the k-subsystem (32) of linear equations. Fortunately, in typical cases the Fourier components U n of the periodic parameter decreases rapidly so, the number of equations retained in the linear system (32) usually reduces to a few ones.

1st

Brillouin zone as introduced in the Fourier domain methods. The ''raison d'être'' of the Brillouin zones is the periodicity of the solution in the Fourier domain. In other words the complete solutions are totally defined in any one of these periods. This is analogue to the conclusion introduced in Section 2.1.2.5. Here, we present the different Fourier domain approaches generally used to ascertain this periodicity.

In the first approach we consider the k-subsystem (32). The substitution of k by k+nb = k ′ is equivalent to a permutation of the subsystem rows. Thus the end result is an equivalent subsystem with the same eigensolutions. This reveals the periodicity of the eigenvalues E nk in the Fourier (reciprocal) space. The second subindex k refers to the subsystem being under consideration and the first subindex n to the classification of the eigenenergies by ascending order. Obviously, each eigenvalue E nk is associated to an eigensolution giving the set of the Fourier (i.e. spectral) components ϕ E,(k+nb) which completely define the periodic part µ k (x).

Other methods consist in considering any of the Eqs. ( 33) to [START_REF] Ortigueira | The comb signal and its Fourier transform[END_REF]. The simplest one is to consider Eq. ( 33) and substitute k by k + mb. This is equivalent to a change of the summation dummy variable. We directly obtain:

ψ E,k+mb (x) =  n ϕ E,[k+(m+n)b] e i[k+(m+n)b]x =  n ′ ϕ E,[k+n ′ b] e i[k+n ′ b]x = ψ E,k (x) . (37)
In other words, the Bloch functions are invariant if the k-value is modified by any reciprocal vector G = nb. Finally, as in Section 2.1.2.5, this periodicity justifies the restriction of k to its principal values that is the first Brillouin Zone:

-b/2 ≤ k ≤ b/2.

Concluding remarks

A plus point of the Fourier domain approach resides in the fact that it gives a new regard to the reciprocal space. It conveys an improved understanding of the physical meaning for the k-vector which we have introduced in the real domain method as a more or less mathematical trick. This enables one not only to take advantage of all the mathematical techniques and tools developed in this field but also to exploit the properties relating the real and Fourier domain representations to make theoretical predictions. Sections 3 and 4 are devoted to this two domains correspondence.

Brillouin vision: or the ''real and Fourier domain correspondence''

In the historical survey concerning the special case of elastic wave propagation in an array of regularly spaced point masses, Brillouin [START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF] provides an explanation of the reciprocal space periodicity of the wavefunction that is worthy of note. He underlines that the physical reason of this periodicity resides in the fact that the particle displacements are measurable only at the ''discrete'' regularly spaced positions of the point masses. In other words, the periodicity of the angular frequency ω with respect to the wave number ℵ = 1/λ (or equivalently the wave vector k = 2π /λ) appears as a direct corollary of the discretization of the wavefunction in the real space domain.

Also, throughout his work, Brillouin [START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF] has recurrently exploited the analogies between the mechanical and electrical systems. Indeed, at this time electronics was extensively used to investigate physical problems of other natures. 6 These analogies are based upon the mathematical one to one correspondence of the elements and the measurable quantities of the physical objects under consideration with the electronic basic elements (resistances, inductances, and capacitances) and the measurable electrical quantities (voltage, current intensity, electric charge. . . ). Several examples are given in Ref. [START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF], where mechanical situations are represented by analog electrical filters.

Proceeding by analogy is very helpful to grasp the subtleties and the intricacy of the two domain correspondence. Indeed, electronic systems have been extensively studied and a lot of useful mathematical tools have been developed for electronic filters as well as for transfer lines analysis. The most striking example is the Fourier-transform correspondence between the time impulse response of a filter and its transfer function i.e. its frequency response. It constitutes the cornerstone of understanding the existence of pass-band, cutoff frequencies and stop band (i.e. forbidden band or bandgap). Also, the concept of transfer function not only allows but supports the harmonic analysis of a system i.e. to focus the analysis on a single Fourier component of the function and analyze the behavior of the transfer function as a function of the Fourier domain independent variable: frequency in the case of time domain functions or wavevector in the case of real space domain functions. The underlying principle is the representation of a given function by its Fourier expansion. 7 These tools, being exclusively based on Fourier transform considerations, are easily transposable to any other physical domains.

sin(k-b)x sin(kx) sin(k+b)x
In this section, following the same illustrative example used by Brillouin [START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF] i.e. elastic wave propagation in an array of regularly spaced point masses, we summarize the fundamental points he used to reveal the physical properties of band structures i.e. its periodicity, existence of allowed and forbidden bands, and the relation between the number of bands and the media degrees of freedom.

Brillouin correspondence principle between the two domains: ''Origin of the periodicity in the Fourier domain''

Brillouin had illustrated the origin of the periodicity in the Fourier domain using the above-mentioned example of regularly spaced point mass array. He examined the propagation of a Fourier component of the displacement u k (x, t) = Ae i(ωt-kx) where A is the amplitude and plot the corresponding sinusoidal curves u k n (x, t) (at a fixed time) for different k nvalues such that: k n = k + 2nπ /a. Where: -π /a < k ≤ π /a, and a is the distance between two adjacent particles and n is an integer. As illustrated in Fig. 2 all these sinusoidal curves will look like very dissimilar. Higher absolute k-value plots will appear more shrunk along the x-axis direction. But a more careful inspection shows that all these curves intersect at the discrete abscissa x m = ma.

The reason of these intersections becomes obvious if we focus on the values of the Fourier component u k n (x m , t) at the discrete abscissas x m = ma, i.e. we consider the series:

u k n ,m (t) ≡ u k n (ma, t) = Ae i(ωt-mk n a) . ( 38
)
It is clear that the series u k n ,m (t) remain invariant for different values of k n = k + 2nπ /a, even if the corresponding functions u kn (x, t) are quite different. Indeed, the product (k n a) differs from  ka  by integer multiples of 2π . In other words, the series (38) (i.e. the discretization of the function u k n (x, t) at the abscissa x = ma) are equivalently defined by any of the Fourier components k n = k + 2nπ /a. This corroborates the fact that the 2π /a periodicity of the Fourier domain expression of the displacements is a direct result of discretization in the direct domain. In conclusion, the solution and consequently its angular frequency ω are both periodic with period (b = 2π /a) with respect to k.

One who is familiar with discrete signal analysis and/or digital electronic filters fields will identify the well-known sampling phenomenon. This shows how much Brillouin had taken advantage of its experience on the electronic filters. Even if, at this time, the discrete signal analysis had not yet achieved the great development which it has reached today with the advent of the fully digital systems (at this time, electronic filtering concerns analog rather than digital filters).

The first Brillouin zone convention, flat band location on zone boundaries and existence of slow waves

The periodicity of the wavefunction in the Fourier domain allows one to restrict the analysis on a unique period (b = 2π /a). On the other-hand, in order to comply with the fact that the behavior of a propagating wave in an unlimited media is invariant whatever it propagates along the positive or the negative x-axis, Brillouin decided to define the principal value of k to be centered on k = 0 i.e. -π /a < k ≤ π /a. According to this convention named after him ''the first Brillouin Zone (FBZ)'', the dispersion relation ω(k) turns up to be an even function.

The conjunction of the periodicity and the even parity of ω(k) gives rise to the localization of extrema (maxima or minima) every k = nπ /a; n ∈ Z including zero. These extrema points explain the existence of ''flat bands'' (i.e. dω/dk = 0) at the edges of the 1st Brillouin Zone as well as at k = 0. The condition dω/dk = 0 is equivalent to say that wavefunctions exhibit zero group velocity at the band edges. This explains the existence of slow waves near band edges which tend toward standing waves as they approach the band edges.

Bandgaps and enumerating the allowed bands

Throughout this paper we consider lossless unlimited media. So, by propagating wave, we mean a wave traveling without any attenuation. That is, its harmonic components are of the form u(x, t) = Ae i (ωt-kx) where k is a real number. In some situations, the solutions give complex or imaginary k-values. In the former case, the wave is attenuated and vanishes after some damped oscillations over relatively short distances and in the later case the media does not oscillate at all but the displacement decreases exponentially. In both cases, these waves are classified as non propagating waves. The range of angular frequencies giving complex or pure imaginary k-values constitutes what we call a bandgap. In ideal lossless media there is no absorption. This attenuation must be viewed as the inability for a local vibration to be transmitted step by step to its neighborhood in a certain frequency range and thus explains the terminology of stop band or forbidden band. In other words, to be able to propagate, a harmonic component has to fulfill a drastic condition: the ''dispersion relation'' relating its frequency to its wavelength as introduce in Section 2.2.1.1. So, the existence of bandgaps is not so difficult to admit. It stands for the range of frequencies where the dispersion relation cannot be satisfied.

Allowed band, bandgap and cutoff frequencies

Brillouin enlighten the uniqueness of the frequency band in the case of a Kelvin array of identical and regularly spaced point masses. 8 He also explains the existence of a high-cutoff frequency. For the simplest case where the interaction is limited to the nearest-neighbors of the regularly spaced point masses, the understanding of this classification in frequency bands is straightforward. Equating the force of inertia of a point mass to the internal forces exerted by its nearest-neighbor [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF][START_REF] Royer | Elastic Waves in Solids. Part 1: Free and Guided Propagation[END_REF] the wave equation (which, in this special case, is equivalent to the equation of motion of one of the particles) reduces to:

m ∂ 2 u n ∂t 2 = η (u n-1 -2u n + u n+1 ) (39)
where, m is the particle mass and η is analogue to a spring constant. Substituting u(x, t) by its expression [START_REF] Ho | Existence of a photonic gap in periodic dielectric structures[END_REF] we directly obtain the dispersion relation:

ω = 2  η m     sin  ka 2     . ( 40
)
Considering Eq. ( 40) the high-cutoff frequency ω c corresponds to the maximum of the sine function. This gives us:

ω c = 2  η/m. (41)
Indeed, if ω > ω c Eq. ( 40) imposes a complex values for k. The expression of k writes (Appendix F):

k = ± π a ± i 2 a cosh -1  ω ω c  . (42)
The complex nature of k means that the harmonic component is attenuated exponentially as it propagates. So the cutoff frequency is the ''frontier'' separating allowed and forbidden bandgap. Also, Eq. ( 40) highlights that the flatness of the dispersion relation coincides with k = ±π /a (maximum of the sine function), in agreement with the conclusion of Section 3.2.

The multiplicity of allowed bands and the number of distinct point masses in a unit cell (degrees of freedom)

Again considering the simplest case where the interaction is limited to the nearest-neighbors of the regularly spaced point masses. But this time, the point masses constitute an alternate sequence of two distinct masses m 1 and m 2 . We can define the unit cell of the system in such a manner that the point masses m 1 are located on the unit cell boundaries (i.e. at

x n = na) and m 2 at the middle of the periods (i.e. at x n+1/2 = (n + 1/2) a) (See Fig. 3).

Two harmonic components of the form of Eq. ( 38) are introduced, one for each particle: where, the first subscript refers to the unit cell while the second labels the type of particle. At a fixed angular frequency ω, these two components may pertain to the same mode only if they propagate with the same phase velocity v φ = ω/k i.e., they have the same wavevector k. So, they can only differ by their amplitudes. The overall analysis is the same as in the preceding section. But, due to the interaction between the two types of particles, the equation of motion (39) splits in a system of two coupled equations:

u n,1 (t) ≡ u 1 (na, t) = A 1 e i(ωt-kna) and u n,2 (t) ≡ u 2 (na, t) = A 2 e i(ωt-k(n+1/2)a) (43)
       m 1 ∂ 2 u n,1 ∂t 2 = η  u n,2 -2u n,1 + u n-1,2  m 2 ∂ 2 u n,2 ∂t 2 = η  u n+1,1 -2u n,2 + u n,1  . ( 44
)
Substituting u n,1 and u n,2 , the dispersion relation come out as the condition of existence of solutions for a 2 × 2 system of linear equations. This dispersion relation writes (Appendix G):

ω 4 -2η  m 1 + m 2 m 1 m 2  ω 2 + 4 η 2 m 1 m 2  sin 2 (ka/2)  = 0. (45)
The two possible solutions: ω 2 + (k) and ω 2 -(k) of Eq. ( 45) give rise to two allowed bands (Appendix G). It is easy to show that when k increases from k = 0 to k = π /a the lower band frequency ω -(k) increases from ω 0 = 0 to a maximum ω 1 = √ 2η/m 1 where, m 1 being the heavier of the two masses and at the same time, the upper band frequency ω + (k) decreases from a maximum value [3]. On the other-hand, if we invert Eq. ( 45) i.e. we take ω as the independent variable. This reads:

ω 3 = √ 2η (m 1 + m 2 ) /m 1 m 2 to a minimum value ω 2 =  2η/m 2
sin 2 (ka/2) =  m 1 + m 2 2η  -  m 1 m 2 4η 2  ω 2  ω 2 , (46)
it is straightforward to demonstrate [START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF] that the values of k becomes complex if ω belongs to the interval ω 1 < ω < ω 2 while it is purely imaginary if ω > ω 3 . These two frequency intervals correspond to bandgaps.

So, we say that doubling ''the degrees of freedom'' (represented here by the number of particles per unit cell) splits the allowed band in two distinct bands. Then, we can logically admit that increasing the number of particles per unit cell to 3 will end with a 3 × 3 system of linear equation and thus 3 bands, and so on. . . . Finally, arrive at the conclusion that for N particles per unit cell we will end with an N × N system with N allowed bands (3N if we consider the transverse displacements).

About the Bloch behavior of the solution in the Brillouin development

Brillouin had not given a lengthy dissertation on Bloch behavior of the wavefunction. But the method he used implicitly includes the physical Bloch nature of the propagating modes.

In his model describing the microscopic wave equation, he introduces by physical arguments the general form of the expected Fourier components of the solution. More especially, he specifies the required relative phase shift between the different harmonic components such that these harmonics belong to the same propagating mode. Eq. ( 43) is a good case in point.

A versatile formalism well adapted to investigate the two domain correspondence in the field of wave propagation in periodic media

We place ourselves within the two domain vision of wave propagation in periodic media initiated by Brillouin. We take advantage of the now available mathematic tools inherited from discrete signal analysis and/or digital electronic filters and exploit the correspondence between the real and reciprocal domains in order to investigate the wave behavior in periodic media. Unlike the classical methods introduced in Section 2 which are intrinsically linked within a specific procedure of finding wave equation solutions, the formalism that we develop here (like Brillouin approach) is essentially based on general considerations relating the two domain representations of a function (i.e. on Fourier analysis techniques). And so, one of the advantages of our formalism resides in the fact that it may be easily transposed to any physical domains: Solid state, photonic or phononic problems.

Specifically, we make use of properties associated to the Dirac-comb: replicating, sifting and Fourier transform properties [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF]. The most important tools along with their main properties are tabulated in Appendix H. Proceeding this way, we associate physical meaning to mathematical expressions. This enables one to have a larger look of the problem, anticipate the results and introduce shortcuts.

We restrict this section to an extension of the Brillouin illustrative example studied in Section 3.3.2. That is, the special case of the elastic waves propagation in a periodically varying mass density along a given axis. Obviously, this model is not realistic in the case of continuous media. In fact, continuous media implies a very high density of particles and so, the interaction cannot be restricted to the nearest neighbors. The extension to more realistic models is delayed to Section 5 which is devoted to solid-state, photonic and phononic crystal. In the first part of this section, we introduce a real domain closed-form expression (RD-CFE) describing the periodic parameter of the media. Then we highlight the correspondence with the physical conclusions established in the preceding section i.e. how this RD-CFE implicitly fulfills the Brillouin vision of the Fourier domain periodicity in case of periodic media (not necessarily discrete). Finally, we introduce a Fourier domain closed-form expression (FD-CFE) of the Bloch theorem as a direct result of the Brillouin two domains correspondence principle (presented in Section 3.1); the real domain expression of the Bloch theorem appears as the inverse Fourier transform of this FD-CFE of Bloch function.

Real domain closed-form expression of a periodic parameter

Consider a man-made material in which the mass m(x) varies periodically along the x-axis. We introduce the real space domain closed-form expression (RD-CFE) of m(x) using the replicating property of a Dirac-comb see Appendix H, [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF]. That is, periodically replicating its restriction to a unit cell m(x) via its convolution products with a Dirac-comb  n δ (xna). The underlying mathematical concept resides in the property of the convolution of a given function say m(x) with a Dirac-delta function 9 δ (xna) centered at x = na. This convolution is known to return the shifted image m(xna) of the function centered at this abscissa [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF]. In other words, we reconstruct the periodic function m(x) using the real space domain closedform expression (RD-CFE):

m(x) = m(x) ⊗  n δ (x -na) (47)
where, the symbol ⊗ denotes the convolution product.

So defined, the RD-CFE closely matches the physical aspect of periodic media. Indeed, a crystal is described by a lattice of nodes exhibiting the translational symmetry (i.e. the periodicity) then, a unit cell (describing the distribution over an individual period a) is assigned to each node. In the same way, the Dirac-comb  n∈Z δ (xna) in Eq. ( 47) depicts the lattice symmetry, while the convolution product mathematically assigns the restriction m(x), which represents the unit cell, to each Dirac-delta function.

Correspondence with the Brillouin approach: discretization in the real domain ⇔ periodicity in the Fourier domain

Inside a given period we can express the restriction using the convolution with a Dirac-delta function:

m(x) =  period dχ m (χ ) δ (x -χ ) . ( 48
)
Considering a discretization along the x-axis with a uniform step 10 dχ q → a/N i.e., N steps per period a, the integral in [START_REF] Wilcox | Scattering Theory for Diffraction Gratings[END_REF] will turn to be a discrete sum:

m(x) =  q a N m  q a N  δ  x -q a N  ⇔ m(x) =  q a N ρ  q a N  (49)
where the discrete abscissa χ q ≡ qa/N designates to the qth mass element and ρ(x

) = m  q a N  δ  x -q a N
 represents the mass density at this abscissa. Substituting in Eq. [START_REF] Linton | Lattice Sums for the Helmholtz Equation[END_REF] gives us:

m(x) = N-1  q=0 a N m  q a N  δ  x -q a N  ⊗  n δ (x -na) . ( 50
)
9 A generalized function in the sense of distribution functions [START_REF] Schwartz | Généralisation de la notion de fonction, de dérivation, de transformation de Fourier et applications mathématiques et physiques[END_REF].

10 These steps are considered from a macroscopic point of view. That is, their dimensions are assumed sufficiently large as compared to inter-atomic dimensions.

Using the convolution property of two Dirac-delta functions δ (x

-x 1 ) ⊗ δ (x -x 2 ) = δ (x -x 1 -x 2 )
and rearranging the terms we obtain:

m(x) =  q   n a N m  q a N  δ  x -na - qa N   . (51)
The inner summations (one for each fixed q-value) represent juxtaposed arrays of elementary masses: ∆m

 χ q  = a N m  χ q  δ 
xχ q  equally spaced by a distance a. As shown in Section 3.3.2, each of these arrays is characterized by a unique displacement expression u k (x, t) ≡ u k (x, t) = A q e i(ωt-kx) . The later having physical existence only at discrete abscissa:

x n,q = na + qa N and so, reduces to a discrete series according to (43):

u n,q (t) ≡ u q (na, t) = A q e i(ωt-nka) (52) where, we have absorbed the phase shift e -ikaq/N in the complex amplitude A q . So, we retrieve the same inherent sampling principle in the real domain emphasized by Brillouin and according to the sampling theorem, their Fourier domain representation should be periodic with a period b = 2π /a. Again, in order to constitute a propagating mode different q-value arrays (i.e. different masses), even if they may exhibit different complex amplitudes, have to oscillate with the same frequency and wavevector i.e. u n,q ′ (t) ≡ u q ′ (na, t) = A q ′ e i(ωt-kna) 0 ≤ q ′ ≤ (N -1). Consequently, the problem splits in a system of N coupled equations. Finally, since the pitch between two distinct arrays is vanishingly small in the continuum (N → +∞), the number of bands tends toward infinity.

The two domains correspondence for the periodic parameter 4.2.1. Representation of the direct and the reciprocal lattice by Dirac-combs

A fundamental point concerning the proposed formalism resides in the Fourier transform of a Dirac-comb which is known to be itself a Dirac-comb. For completeness, a brief demonstration is given in Appendix I where the Dirac-delta is approached as the limit of a sequence of functions [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF][START_REF] Ortigueira | The comb signal and its Fourier transform[END_REF]. This Fourier transform writes:

F   n δ (x -na)  = b  m δ (k -mb) ; b = 2π /a. (53)
This shows that the analogy we pointed out in the real domain, i.e. the representation of the Bravais lattice by a Dirac-comb, may be pursued in the Fourier domain. Indeed, we notice the correspondence between the spacing between two successive Dirac-delta functions in Eq. ( 53) and the reciprocal lattice basis vector b = 2π /a. Thus, this Fourier domain Dirac-comb may be viewed as the representation of the reciprocal lattice. This analogy is still valid in case of 3-Dimensional media since the following relation holds in this case [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF][START_REF] Ortigueira | The comb signal and its Fourier transform[END_REF][START_REF] Baake | Mathematical diffraction theory in Euclidean spaces: an introductory survey[END_REF]:

F   R∈R δ (r -R)  ⇔ v RL  G∈G δ (k -G) (54)
where, R and G are the direct and reciprocal lattice vectors as defined in Section 2.1.3 and describe respectively the nodes positions in the two domains and v

RL = b 1 • (b 2 × b 3 ) = (2π ) 3 /v DL is the unit cell volume in the reciprocal lattice while, v DL = a 1 • (a 2 × a 3
) is the unit cell volume in the direct lattice.

A direct inspection of the Fourier transform pair in Eq. (54) suffices to convince oneself how much the Dirac-comb tools are well adapted to the implementation of periodic media: the concept of the reciprocal lattice vectors is already included. The Dirac-delta functions of the right-hand side are located at the nodes of the reciprocal lattice as dictated by the expression G = m 1 b 1 + m 2 b 2 + m 3 b 3 and the common factor v RL specifies their strength. Hence, the determination of the reciprocal lattice basis vectors is considered as an already solved mathematical problem. This observation, associated to the self-transforming of the Dirac-comb explains the straightforwardness the proposed formalism.

Fourier Domain closed-form expression for the periodic parameter

We define the Fourier domain closed-form expression (FD-CFE) as the Fourier transform of the RD-CFE given in Eq. ( 47). Since the property of the Fourier transform of a convolution product holds in the field of the distribution functions, the convolution product in (47) converts to the simple product of the Fourier transform m(k) of the expression of the mass restriction m(x) by the Fourier domain Dirac-comb (53). Finally, the Fourier transform of the periodic parameter turns out to be a weighted Dirac-comb:

m(k) = bm(k)  m δ (k -mb) =  m bm(mb)δ (k -mb) .

Bloch theorem and the two domains correspondence principles 4.3.1. Bloch theorem in the Fourier domain as a corollary to the periodicity of the band structure established by Brillouin

As introduced in Section 3.1, Brillouin has established that the periodicity of the band structure (in the Fourier domain) is a direct result of an implicit discretization of the wavefunction in the real space domain. In Section 4.1.1, we determine the correspondence between our formalism and the Brillouin approach. This enables one to ascertain, if any proof were still necessary, that the periodicity of the band structure ω(k) established by Brillouin is still valid in the case of periodic but not discrete media. This constitutes the Brillouin correspondence principle between the two domains (real and Fourier domains). Now, consider the inverse problem. That is, given a mode labeled by its angular eigenfrequency ω, we look for the corresponding wave vectors k. According to the above mentioned periodicity, k(ω) is a multi-valued function. More specifically, for a fixed angular eigenfrequency ω, the k values are of the form:

k = k -mb; m ∈ Z (56)
where, the continuum k stands for the principal value of the wavevector i.e. the one within the first Brillouin zone and b = 2π /a is the reciprocal lattice basis vector.

In other words, the general form of a propagating mode consists of all the k-values of the form of Eq. ( 56). Thus, carrying on the problem of periodically varying mass m(x), the displacement spectrum u(k) will be composed a set of discrete components of the form:

u( k -mb).
(57) This point may be viewed as an intrinsic discretization of the wavefunction in the Fourier domain as a corollary of the band structure periodicity. But this discretization is shifted with respect to the origin by the principal value k of the wavevector.

The transcription of these results, in the framework of distribution functions, writes:

u k (k) =  m∈Z u  k -mb  δ  k -  k -mb  . (58)
In other words, the Fourier domain general expression of the wavefunction u k (k) appears as a weighted Dirac-comb. Thus, u k (k) assembles in a closed form all the discrete components u( kmb) constituting an eigenmode. Now, we emphasize that periodicity is the unique constraint we have imposed on the medium. So, this last relation necessarily involves the same physical contents as the Bloch theorem. In other words, expression (58) is nothing but the Fourier domain closed-form expression (FD-CFE) of the Bloch theorem as it will be confirmed in the following section.

The real domain Bloch theorem as a Fourier inverse transform of the FD-CFE

Using the property F -1 {δ [k ± B]} = (1/2π) e ±iBx , the inverse Fourier transform of the FD-CFE directly gives us:

u k (x) = 1 2π   m∈Z u  k -mb  e imbx  e -i kx (59)
where, we have factorized e -i kx . Now, recalling that b = 2π /a, the summation between brackets represents a Fourier series expansion and so, it stands for a periodic function µ (x) of period a µ (x + a) = µ (x) .

(60)

Finally, Eq. ( 59) exhibits the Bloch nature of the displacement wave and confirms our previous allegation. We can write:

u k (x) = 1 2π µ (x) e -i kx . ( 61
)
The minus sign in the phase factor in Eqs. ( 59)-( 61), as compared to the plus sign obtained in Eqs. ( 34)- [START_REF] Schwartz | Généralisation de la notion de fonction, de dérivation, de transformation de Fourier et applications mathématiques et physiques[END_REF], is a consequence of the different Fourier transform conventions used in solid-state and phononic crystals (see footnote 3 in Section 2.2.1 and Appendix H).

Bloch theorem as a result of the conjunction of sampling and shift Fourier transform theorems

The inherent discretization pointed out in the Section 4.3.1 enables one to establish the real domain Bloch theorem as a second principle of the two domain correspondence. Eq. ( 58) can be viewed as a periodic sampling of u n (k) in the Fourier domain with a sampling rate b = 2π /a dictated by the periodicity a of the media but, these samples are collectively shifted by a wavevector k in the Fourier Domain. So, in the real domain, the periodic part of the Bloch function results from the discretization of the Fourier spectrum and the phase factor originates from the collective shifting of these discrete Fourier components. We will return to this issue in Section 6.2.2.

Conclusion

The main advantage of the approach presented in this section resides in the fact that it is based on general considerations of Fourier analysis techniques relating the two domain representations and is not embedded within a mathematical procedure of solving wave equations. This permits to have a good perception of the phenomenon of wave propagation in periodic media and thus leads to short cuts. For instant, this highlights the universal character of the Bloch theorem since it appears as a direct consequence of general Fourier analysis theorems, namely: the combination of the well-known shift and sampling theorems.

In addition, this approach introduces other short cuts which are especially appreciable in complex cases as it will be illustrated in Section 5. The use of RD-CFE and FD-CFE, will end with a universal systematic two step procedure applicable in the treatment of periodic media. The method is so straightforward that, even in the most complex cases, the reciprocal lattice wave equation can be established by a simple inspection of its expression in the direct domain. Finally, this formalism permits to preserve the two domain correspondence at different levels: it concerns not only the wavefunction representation but also the periodic media itself with a clear dissociation between the lattice and the unit cell expressions in the two domains.

Solid state, photonic and phononic periodic media in the formalism of the closed-form expressions

In this section we illustrate the application of the method introduced in Section 4 to Solid-state, Photonic and phononic materials. In Section 5.1 we examine the simplest case of scalar electronic de Broglie waves propagation in one-dimension periodic potential and widen the method to three dimensions. The expressions of the Bloch theorem in both the real and Fourier domain are derived. The extension to vector fields is presented in Section 5.2 where, we consider electromagnetic fields in photonic crystals in the special case where the material does not exhibit any other anisotropy than that introduced by the periodic non-homogeneities. Finally, Section 5.3 addresses phononic crystals as an illustration of the general case of vector field in periodic media where the constitutive materials exhibit anisotropic as well as piezoelectric properties.

Solid state crystals

In the first subsection we apply the formalism in a one-dimension periodic solid-state media. Then we show that the extension of the method to the three-dimension crystals is straightforward and does not require special efforts.

The de Broglie waves in 1-d unlimited solid state crystals

As in Section 2, we consider the one dimensional Schrödinger equation ( 1) but, this time in the framework introduced in Section 4

- h2 2m ∂ 2 ψ (x) ∂x 2 + U (x) ψ (x) = Eψ (x) . (1)
The first step consists of writing closed form equations of the periodic parameter in the two domains i.e. the RD-CFE U (x) and its Fourier transform the FD-CFE U (k). These closed-form expressions write:

U (x) = u (x) ⊗  q∈Z δ (x -qa) ⇔ U (k) = u (k)  q∈Z bδ (k -qb) = b  q∈Z u (qb) δ (k -qb) (62)
where, we identify the Fourier transform pairs:

• U (x) ⇔ U (k):
The two-domain expressions of the periodic potential.

• u (x) ⇔ u (k): The two-domain contribution of a unit cell to the periodic potential

•  q∈Z δ (x -qa) ⇔  q∈Z bδ (k -qb):
The two-domain lattice representations, the direct and reciprocal ones with their respective basis vectors a ⇔ b = 2π /a (scalar in the 1-D case) Using Eq. ( 62) and the Fourier transform of the differential operator recalled in Appendix H, the Fourier transform of Schrödinger equation is straightforward. Indeed, the simple product ψ (x) U (x) in the real space transforms into a convolution product [U (k) ⊗ ψ (k)] /2π in the reciprocal domain then using the replicating property Appendix H (in the reciprocal domain) of the wave function with the weighted Dirac-comb directly gives us:

h2 k 2 2m ψ (k) + b 2π  q∈Z u (qb) • ψ (k -qb) = Eψ (k) (63)
where, ψ (k) stands for the Fourier transforms of ψ (x).

Obviously, Eq. ( 63) corresponds exactly to the central equation ( 32) obtained in Section 2.2 by the traditional method. But, with the fundamental differences, that in our approach we use neither virtual box concept nor the Fourier series. The straightforwardness is evident as compared to the classical method and it will be confirmed in the 3-D case which is still a two step procedure while the classical approach becomes more and more complex (see Appendix E). Note that the factor b/2π = 1/a which appears in Eq. ( 63) fulfills the correspondence between Fourier transform and the coefficients of Fourier series expansion used in the traditional method.

• Bloch Function in Fourier domain.

As introduced in Section 4.3, taking account of the intrinsic k-shifted discretization of the wavefunction kn in the Fourier domain, it may be written in the form of the weighted Dirac-comb:

kn  k ′  =  q∈Z ψ n (k -qb) δ  k ′ -(k -qb)  (64)
where, the weighting factors stand for the unknowns of the linear equations system (63). In other words, concerning the linear system of equations, kn may be viewed as the eigenvector whose coordinates correspond to the eigensolutions ψ n (kqb). As explained in Section 4.3.1, Eq. ( 64) stands for the FD-CFE of the Bloch theorem, its inverse Fourier transform gives the real space Bloch function:

kn (x) = 1 2π   q∈Z ψ (k -qb) e -iqbx  e ikx = µ (x) e ikx ; µ (x + a) = µ (x) . ( 65
)
Note the sign inversion with respect to Eqs. ( 60)-( 61) since we use the property:

F -1 {δ [k ± B]} = (1/2π ) e ∓iBx relative to
the Fourier transform convention in solid-state crystals (see Appendix H).

Three-dimensions extension

Following the same procedure, we consider the 3-D time independent Schrödinger equation:

- h2 2m ∇ 2 ψ (r) + U (r) ψ (r) = Eψ (r) . ( 66
)
We begin with the introduction of the three-dimension RD-CFE and FD-CFE for the periodic potential:

U (r) = u (r) ⊗  R∈R δ (r -R) ⇔ U (k) = u (k)  v RL  G∈G δ (k -G)  = v RL  G∈G u (G) δ (k -G) . ( 67
)
Again we identify the Fourier transform pairs:

• U (r) ⇔ U (k): The two-domain expressions of the 3-D periodic potential.

• u (r) ⇔ u (k): The two-domain contribution of the 3-D unit cell to the periodic potential. Following the same procedure as in the one-dimension case, the simple product ψ (r) U (r) in Eq. (66) converts into a convolution product [U (k) ⊗ ψ (k)] / (2π) 3 (see Appendix H); where U (k) is the weighted Dirac-comb given in Eq. (67).

•  R∈R δ (r -R) ⇔ v RL  G∈G δ (k -G)
Like in the 1-D case, we obtain the 3-D Schrödinger equation in the Fourier domain by simple inspection:

h2 k 2 2m ψ (k) + v RL (2π ) 3  G∈G u (G) ψ (k -G) = Eψ (k) . ( 68
)
As before, we can write the FD-CFE for the Bloch theorem in the 3-D space:

kn  k ′  =  G∈G ψ n (k -G) δ  k ′ -(k -G)  (69)
which straightforwardly inverse Fourier transforms to:

n,k (r) =  1 (2π) 3  G∈G ψ n (k -G) e -iG•r  e ik•r . ( 70
)
Again, identifying the brackets to a 3-D periodic function µ kn (r) we obtain the Bloch function in real space:

kn (r) = µ kn (r) e ik•r . ( 71
)
Following the formalism, the treatment of the 3-D problem is analogue to the 1-D one and does not complicate the solution process. In both cases the Bloch nature of the solution appears as an intrinsic character of wave propagation in periodic media. As indicated in Section 4.3.3, the real Bloch waveform originates from, the second principle of the two domain correspondence i.e., the conjunction of the sampling and shifting theorem in the Fourier domain.

Photonic crystals

In what follows we will point out the appreciable straightforwardness of the method when applied to more complex problems such as vector fields.

Vector field equation in the real space domain

Usually, the electromagnetic field equations may be expressed in terms of H [START_REF] Ho | Existence of a photonic gap in periodic dielectric structures[END_REF] (or E [START_REF] Sözüer | Photonic bands: Convergence problems with the plane-wave method[END_REF]) by respectively eliminating electric (or magnetic) field from Maxwell's equations. Traditionally, one use the terminology introduced by Sözüer et al. [START_REF] Sözüer | Photonic bands: Convergence problems with the plane-wave method[END_REF] and refers as ''H-method'' to the procedure consisting in solving wave equation for H, and then finds out E. The opposite method is known as ''E-method''. The H-method shows two advantages concerning numerical calculations: First, it leads to an ordinary Hermitian eigenproblem whereas; the E-method leads to a generalized Hermitian eigenproblem [START_REF] Sözüer | Photonic bands: Convergence problems with the plane-wave method[END_REF]. Second, one can take advantage of the transversality property of the magnetic vector k • H = 0 while the inhomogeneities introduced by the periodicity impose that the electric field is not transverse [START_REF] Rigneault | La nanophotonique[END_REF]. Here, we will point out the versatility and the straightforwardness of our formalism when applied to the H-method. In the case of nonmagnetic, non dispersive and sources free dielectrics unlimited materials with no other anisotropy than the 3-D periodic non-homogeneity; only the relative permittivity ε r (or equivalently the impermittivity η = ε -1 r ) is position dependent. Thus, the real space wave equation reads [START_REF] Ho | Existence of a photonic gap in periodic dielectric structures[END_REF]:

∇ × [η (r) ∇ × H (r)] + ε 0 µ 0 ∂ 2 H (r) ∂t 2 = 0.
(72)

Vector field equation in the reciprocal domain

As in the previous section, the method begins by introducing a RD-CFE for the periodic parameter η (r) using the convolution products of a Dirac-comb displaying the 3-D periodicity with the position dependent functions υ (r) pertaining to a unit cell. Thus we can write:

η (r) = υ (r) ⊗  R∈R δ (r -R) .
(73) Its Fourier transform straightforwardly gives us:

η (k) = v RL  G∈G υ (G) δ (k -G) . (74)
The simple product [η (r) ∇ × H (r)] in Eq. (72) will transform as a convolution product in the reciprocal domain.

F {[η (r) ∇ × H (r)]} = 1 (2π) 3 η (k) ⊗ [-ik × H (k)] . (75)
Then, using the replicating property with the weighted Dirac-comb (74) we directly rewrite Eq. ( 75) as:

F {[η (r) ∇ × H (r)]} = -iv RL (2π) 3  G∈G υ (G) (k -G) × H (k -G) . ( 76
)
Finally, using the compact notation k 2 0 = µ 0 ε 0 ω 2 = (ω/c) 2 , the Fourier transform of Eq. ( 72) is simply given by:

1 v DL k ×  G∈G υ (G) [(k -G) × H (k -G)] + k 2 0 H (k) = 0. ( 77
)
Eq. ( 77) stands for an infinite subset of coupled linear equations where the unknowns are the magnetic field spectral components H(k -G). Since it is more convenient to deal with scalar quantities, we express the H vectors in terms of their components in the following section.

Expression of the H vector field equation in terms of its scalar components

The application of Fourier transform implicitly assumes that the general solution consists of superposition of plane waves (i.e. plane wave packet). Taking advantage of the transverse nature of the magnetic field of each individual plane wave, the magnetic vector H(k -G) of a given plane wave in Eq. ( 77) may be expressed by its components h 1 (k -G) and h 2 (k -G) along two orthogonal unit vectors, say h 1,k-G and h 2,k-G , laying in the plane perpendicular to its wavevector (k -G). Thus, we substitute in Eq. (77) the following expressions of H(k -G), (k -G) and k:

H (k -G) = h 1 (k -G) h 1,k-G + h 2 (k -G) h 2,k-G (78) (k -G) = |k -G|  h 1,k-G × h 2,k-G  and k = |k|  h 1,k × h 2,k  . (79)
Then, performing the double cross product in Eq. ( 77) we obtain the following vector relation expressed in the {h 1,k , h

 G∈G υ (G)

v DL |k| |k -G|  +  h 2,k • h 2,k-G  -  h 2,k • h 1,k-G  -  h 1,k • h 2,k-G  +  h 1,k • h 1,k-G    h 1 (k -G) h 2 (k -G)  = k 2 0  h 1 (k) h 2 (k)  . (80)
Thus each vector equation of the subsystem (77) splits into 2 scalar equations.

Correspondence with the traditional relations

In our derivation of wave propagation in the reciprocal domain, we have implicitly expanded the field components in Fourier basis functions e -ik•r while, the conventional method (using Fourier series and/or Poisson summation formulas), retains spectral components basis of the form e -i(k+G ′′ )•r . Thus, to show the correspondence with the usual relations [START_REF] Ho | Existence of a photonic gap in periodic dielectric structures[END_REF][START_REF] Sözüer | Photonic bands: Convergence problems with the plane-wave method[END_REF][START_REF] Rigneault | La nanophotonique[END_REF],

we have to change the variables k into k +G ′′ in Eq. (80) and the summation dummy variable G into G ′ = G ′′ -G. We, obtain:

 G ′ ∈G υ  G ′′ -G ′  V DL   k + G ′′     k + G ′    +  h 2,k+G ′′ • h 2,k+G ′  -  h 2,k+G ′′ • h 1,k+G ′  -  h 1,k+G ′′ • h 2,k+G ′  +  h 1,k+G ′′ • h 1,k+G ′    h 1  k + G ′  h 2  k + G ′   = k 2 0  h 1  k + G ′′  h 2  k + G ′′   . ( 81
)
This relation is in full agreement with the results derived by Ho et al. [START_REF] Ho | Existence of a photonic gap in periodic dielectric structures[END_REF]. It is worth noting that the volume of the unit cell VDL which appears in our equations fulfills the correspondence between Fourier transform and the coefficients of Fourier series expansion used in the traditional method. Actually, Eqs. ( 80) and ( 81) are two representations of the same system of coupled linear equations since the summation runs over all the reciprocal lattice vectors G ∈ G. 20

Phononic crystals

This section illustrates the effectiveness of our method when applied to the general case of vector field in 3-D periodic media whose constituent materials being anisotropic without any restriction concerning their symmetry. Actually, our approach turns out to be a systematic ''two steps'' procedure for determining wave equations in the reciprocal space. Thus its straightforwardness is not altered whatever is the complexity of the constituent material anisotropy.

In Section 5.3.1, as an introduction we consider the simple case of non-piezoelectric materials. Then, the most universal case, where the constituent materials are piezoelectric, is addressed in Section 5.3.2.

ρ (r) ∂ 2 u i (r) ∂t 2 =  j ∂σ ij (r) ∂x j . (82)
where, ρ(r) is the mass density and σ ij (r) the stress tensor. This later, in the linear approximation, writes:

σ ij (r) =  kl c ijkl (r) ∂u k (r) ∂x l (83)
where, c ijkl (r) is the stiffness tensor of rank 4. Eliminating the stress tensor between (82) and ( 83) we obtain the acoustic wave equation expressed in terms of displacements, in the form of the following coupled partial differential equations:

ρ (r) ∂ 2 u i (r) ∂t 2 =  jkl ∂ ∂x j  c ijkl (r) ∂u k (r) ∂x l  . (84)
Again, we introduce a RD-CFE for the periodic parameters ρ(r) and c ijkl (r) using the replicating property of their respective restrictions o(r) and γ ijkl (r) to the unit cell we can write:

ρ(r) = o(r) ⊗  R∈R δ(r -R); c ijkl (r) = γ ijkl (r) ⊗  R∈R δ(r -R).
(85)

Vector field equation in the reciprocal domain.

The FD-CFE are the Fourier transforms of RD-CFE (85) and writes:

ρ(k) = V RL  G∈G o(G)δ(k -G); c ijkl (k) = V RL  G∈G γ ijkl (G)δ(k -G). (86)
Knowing these periodic parameter expressions in the two domains and assuming harmonic dependence of the form e i(ωt-k•r) , the wave equation (84) immediately transform to:

 1 (2π ) 3  V RL  G∈G o(G)δ(k -G)  ⊗  -ω 2 u i (k)   =  jkl  -ik j   1 (2π ) 3  V RL  G γ ijkl (Γ )δ(k -G)  ⊗ [-ik l u k (k)]  . (87)
By direct inspection, one can identify the correspondence between the two domains wave equations (84) and (87). That is, the transformation of the simple product of functions in the real space domain (Eq. ( 84)) into the convolution product of their respective Fourier transforms in the reciprocal domain (Eq. ( 87)) as well as the transformation of the partial derivatives ∂/∂x j to -ik j . Finally, applying the replicating property of Dirac-delta combs, gives us the phononic wave equation in the reciprocal domain:

ω 2  G∈G o(G)u i (k -G) =  jkl  G∈G γ ijkl (G)k j (k l -G l )u k (k -G) . (88)
This illustrates the straightforwardness and the universality of our method. Indeed, whatever is the subject matter under consideration: solid-state, photonic or phononic crystal, the two domain expressions of periodic parameters always display the same general forms (Eqs. ( 85), ( 86)). Thus, only two truly systematic steps are needed to obtain the wave equation ( 88) in the reciprocal domain, knowing its real space expression (84). Eq. ( 88) may be written in the more suitable matrix form for numerical calculations, just by rearranging the terms. We obtain:

 G∈G  ω 2 o (G) I 3 -Λ   u 1 (k -G) u 2 (k -G) u 3 (k -G)  = 0 (89)
where the elements of the matrix Λ are:

Λ ik =  jl γ ijkl (G)k j (k l -G l
) and I 3 is the identity matrix.

Thus, we arrive at a compact equation not only easy to demonstrate but as well easy to handle. One can convince itself, of the effortlessness as well as the effectiveness of our approach, by trying to reproduce the widely published special cases. For instance, starting from Eq. (89) one can directly find out equivalent expressions for the 27 matrices derived in Ref. [START_REF] Wu | Surface and bulk acoustic waves in two-dimensional phononic crystal consisting of materials with general anisotropy[END_REF] for the special case of two-dimensional phononic structures. By equivalent we mean that they lead to the same eigensolutions. Nonetheless, as established in Section 5.2.4, in order to find the strict correspondence, one should first substitute k by (k+G ′′ ) and the dummy variable G by (G ′ = G ′′ -G). Indeed, Ref. [START_REF] Wu | Surface and bulk acoustic waves in two-dimensional phononic crystal consisting of materials with general anisotropy[END_REF] postulates the Bloch theorem i.e. assumes space and time harmonic dependence of the form e i(ωt-(k+G ′′ )•r) whereas our approach is based on e i(ωt-k•r) harmonic dependence.

Phononic crystals with piezoelectric constituent materials

Vector field equation in the real domain.

Considering piezoelectric materials, not only expression (83) must be modified to take account of the electric field E contribution to the stress tensor but also, a subsidiary relation is needed to express the contribution of the strain to the electric displacement vector D. Choosing as independent variables the displacement u k and the electric potential φ instead of the electric field strength E = -∇φ, this system of equations writes [START_REF] Royer | Elastic Waves in Solids. Part 1: Free and Guided Propagation[END_REF]:

       σ ij (r, t) =  kl c E ijkl (r) ∂u k (r, t) ∂x l +  l e lij (r) ∂φ (r, t) ∂x l D i (r, t) =  kl e ikl (r) ∂u k (r, t) ∂x l -  l ε S il ( r) ∂φ (r, t) ∂x l (90)
where, c E ijkl is the stiffness constant at fixed electric fields, ε S il the permittivity at fixed strain and e ikl the direct piezoelectric constants [START_REF] Royer | Elastic Waves in Solids. Part 1: Free and Guided Propagation[END_REF]. Substituting σ ij in the dynamic equation ( 82) and assuming charge free media i.e. zero divergence of the electric displacement vector

 i ∂D i (r, t) /∂x i = 0, we obtain:            ρ (r) ∂ 2 u i (r, t) ∂t 2 =  j  ∂ ∂x j   kl c E ijkl (r) ∂u k (r, t) ∂x l +  l e lij (r) ∂φ (r, t) ∂x l  0 =  i ∂ ∂x i   kl e ikl (r) ∂u k (r, t) ∂x l -  l ε S il (r) ∂φ (r, t) ∂x l  . (91)
In order to determine the electromechanical wave equation in non-homogeneous periodic media, we do not eliminate the electric potential in the coupled equations system (91) as it is the usual procedure in case of homogeneous media. Instead, one [START_REF] Laude | Full band gap for surface acoustic waves in a piezoelectric phononic crystal[END_REF] defines a generalized displacement vector ũ in which u 4 represents the electrical potential φ. i.e. ũ = [u 1 , u 2 , u 3 , φ] T , and a generalized stress vectors in which: σ i4 = D i . i.e. σi = [σ i1 , σ i2 , σ i3 , D i ] T . So doing, the electromechanical wave equation (91) takes exactly of the same general form of the partial differential system (84) but this time, it couples the four variables u 1 , u 2 , u 3 and u 4 = φ:

ρ (r) ∂ 2 ũi (r, t) ∂t 2 =  j=1,2,3  ∂ ∂x j   k=1,2,3,4  l=1,2,3 cijkl (r) ∂ ũk (r, t) ∂x l  . (92)
By identification with Eq. (91), the tilde constants are: (93)

cijkl (r) = c E ijkl (r) ; . . . k = 1,

Vector field equation in the reciprocal domain.

Following our method i.e. introducing the RD-CFE of the tilde periodic material parameters and their corresponding FD-CFE and using the replicating property of the weighted Dirac-comb we obtain the required wave equation:

ω 2  G∈G õ(G)ũ i (k -G) =  jkl  G∈G γijkl (G)k j (k l -G l )ũ k (k -G) . ( 94
)
This system is analogue to Eqs. (88). Again, rearranging the terms it may be written in a matrix form:

 G ′′′ ∈G  ω 2 o (G) Ĩ4 - Λ    ũ1 (k -G) ũ2 (k -G) ũ3 (k -G) ũ4 (k -G)    = 0 (95)
where

Λik =  jl γijkl (G)k j (k l -G l ) and Ĩ4 ≡    1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0    . (96)
As usual replacing k by (k + G ′′ ) and the dummy variable G by (G ′ = G ′′ -G) one obtains the correspondence with Fourier series expansion method. Namely, here Eq. ( 94) is in total agreement with the coupled equation system ( 19) and ( 22) of Ref. [START_REF] Wilm | A full 3D plane-wave-expansion model for 1-3 piezoelectric composite structures[END_REF].

Advantages of our method, as an ending remark concerning ''periodic anisotropic materials'' treatment

This section reveals the effectiveness as well of the effortless of our approach. It is a systematic way for obtaining ''in only two steps'' the reciprocal domain expressions of the wave equation for any periodic media. This, whatever is the complexity of its anisotropy and without any additional intermediate algebra steps. This fact is not surprising since we have turned toward well-known theorems in the field of the distribution functions which have already resolved the mathematical problems. One who is familiar with the conventional method knows how tricky the implicit algebra procedure is. Our twosteps systematic method will be time saving for mathematical formulation of practical problems.

Physical contents of the formalism

Discussion

The formalism is essentially based on Fourier analysis techniques and thus, as illustrated, it suits well any physical domain.

1. Mathematically speaking, taking the Fourier transform of the wave equation implies that we operate in two steps. In the former, we express the real space field under the differential operators in terms of their corresponding inverse Fourier transforms (plane wave packets). In the later, we obtain the specific equation characterizing the subset equations system by taking the Fourier transform i.e., by projecting each term of the differential equations on a selected Fourier basis function. By essence, in 3-D spatial Fourier transform, this selected basis function is e ik•r . In other words we follow the plane wave expansion method (PWE), but in a Fourier continuum basis. In practice one does not have to go through all the steps of this mathematical procedure. But, saves time by directly applying the well known properties of Fourier transforms, namely: transformations of simple or convolution products of functions and distribution functions as well as transformation properties of differential operators. These are summarized in Appendix H.

2. Correspondence with the methods involving Poisson summation formulas.

Other approaches exist which use Fourier transform instead of Fourier series and avoid the use of virtual box concept. They are based on Poisson summation formulas [START_REF] Sözüer | Photonic bands: Convergence problems with the plane-wave method[END_REF][START_REF] Poisson | Traité de Physique Mathématique: Théorie Mathématique de la chaleur[END_REF][START_REF] Mace | Periodically stiffened fluid-loaded plates, 1: response to convected harmonic pressure and free wave propagation[END_REF][START_REF] Linton | One-and two-dimensional lattice sums for the three-dimensional Helmholtz equation[END_REF][START_REF] Linton | Lattice Sums for the Helmholtz Equation[END_REF][START_REF] Wilcox | Scattering Theory for Diffraction Gratings[END_REF]. These approaches express the periodic parameters in terms of distribution functions in the Fourier domain but, in terms of usual functions in the real space. For example, in solid-state crystals, the real domain periodic potential is written as:

U(x) =  n∈N u (x -na) . (97)
Then, the general form of the Poisson summation formulas enables one to convert a summation (of a regularly shifted usual function) on the direct lattice into a sum on the reciprocal lattice. That is, another writing of (97) also in the real domain is:

U(x) =  n∈N u (x -na) = 1 a  q∈N u (qb) e -iqbx . (98)
This constitutes a first step in obtaining the Fourier transform of the periodic parameter:

U(k) = F  1 a  q∈N u (qb) e -iqbx  = 2π a  q∈N u (qb) δ (k -qb) . (99)
This is exactly the same expression we obtained in Eq. (62). Even if the two approaches show the same straightforwardness and gives the same end results. The two frameworks are separately self-consistent even if they can be deduced from one another (as shown in Appendix J). This self-consistency is justified by the fact that, as mentioned in Section 4.2.1 and Appendix I, derivations exist for the self-transforming property of a Dirac-comb [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF] which do not rely; neither on Poisson formulas nor on Fourier series. Beside this overall mathematical equivalence between the two methods differences exist. Taking full advantage of well known properties of Dirac-comb, we introduce a RD-CFE and its corresponding FD-CFE for the periodic material parameters. So doing permits to visualize the two domain correspondence at different levels including lattice, unit cell and the periodic parameter expression as a whole. In other words, the analogy between these mathematical tools and the periodic structures permits the dissociation between the lattice and the unit cell expressions and this dissociation is preserved in both the real and Fourier domains by a one to one correspondence. Furthermore, applying the same Dirac-comb tools to the wavefunction, the Bloch theorem is found to be a direct consequence of the conjunction of the well-known shift and sampling Fourier analysis theorems.

3. Our definition of the function u(r) used in Eq. ( 62) or (67), requires some clarifications. Usually, in solid state crystals the periodic potential U(r) results from the superposition of potential overlaps between the surrounding unit cells contributions u(r). In our formalism, U(r) is defined as the convolution product of a Dirac comb with u(r). So, it is fully adapted to take into account these physically existing overlaps in the case where u(r) spreads out of the unit cell. However, in other situations, like photonic or phononic crystal, it is the restriction to the unit cell u ′ (r) of the periodic parameter which is known rather than the unit cell contribution u(r). So, a very practical implementation of the periodic potential U(r) consists in its reconstruction by a cyclic replication of its restriction u ′ (r) to the unit cell:

U (r) = u ′ (r) ⊗  R∈R δ (r -R).
It is interesting to note, that since our method is developed in the framework of distribution functions, the spectrum discreteness of U(k) appears naturally as the simple product of the of u ′ (k) by a Dirac-comb:

U (k) = u ′ (k)  V RL  G∈G δ (k -G)  = V RL  G∈G u ′ (G) δ (k -G).
Thus, the correspondence with Fourier coefficients is implicit and so, in the framework of the approach presented here, the procedure steps are strictly the same whatever the given function is u (r) or u ′ (r).

The origin of the Bloch theorem

As a consequence of the Brillouin two domain correspondence principle which establishes the periodicity of the band structure, we clarified in Section 4.3 how the spectrum of a given mode with an eigenfrequency ω consists of discrete Fourier components which are equally spaced from one another by reciprocal lattice vectors G and are collectively shifted in the reciprocal domain by a given vector k ̸ = G. These observations are transcribed by Eq. ( 69) whose inverse Fourier transform imposes the real Bloch function behavior. In the following first subsection, we establish the correspondence with the classical Bloch theorem demonstration in Fourier domain method. In the second subsection we give a more detailed view of the two domain expressions of the Bloch theorem.

Comparison with the establishment of the Bloch theorem within the classical method

The correspondence with the classical method is as follows. Starting from Eq. (68), we consider the case of homogeneous materials (no periodic potential: U(r) = 0) there is no coupling at all between individual monochromatic plane waves. We are in the presence of an infinite set of uncoupled equations. Thus, monochromatic plane waves are the eigensolutions. The expression of any equation of this subset in the reciprocal domain gives the well-known dispersion relation: ω = E/h = hk 2 /2m.

In the case of periodic media, the material parameter appears as simple products U (r) ψ (r) in the real space central equation (66) so; it will convert to convolution products in the reciprocal domain: (U (k) ⊗ ψ (k)) / (2π ) 3 . It is this convolution product of ψ (k) with the Dirac-comb, appearing in the expression (67) of U (k) which introduces the selective coupling between all spectral components ψ n (k -G) which are exactly spaced from the selected k Fourier component by regular intervals equal to reciprocal lattice vectors G. Finally, the mathematical transcription of this assertion leads to (69): ''the FD-CFE for the Bloch theorem''.

This writing of the Bloch theorem highlights the dissociation of the contributions related to the lattice structure from those related to the unit cell. The contribution of the lattice structure is depicted by the 3-D sampling intervals represented by the Dirac-comb while, the contribution of the unit cell is characterized by the weighting coefficients ψ n (k -G) which depend via Eq. (68), on the Fourier components u (G) of the unit cell potential distribution.

Bloch theorem in the direct domain

It is instructive to determine the inverse Fourier transform of Eq. ( 69) by another way. Using the convolution property of two Dirac-delta functions δ (yy 1 ) ⊗ δ (yy 2 ) = δ [y -(y 1 + y 2 )], expression (69) may be rewritten as:

kn  k ′  =   G∈G ψ n (k -G) δ  k ′ + G   ⊗ δ  k ′ -k  . ( 100
)
Here the convolution product may be interpreted as the shifting of the weighted Dirac-delta comb (term between brackets) to be centered at the fixed k-value which labels the mode. To who is familiar with Fourier analysis, this form clearly evidences the Fourier correspondence with Eq. (71) in the real space. Indeed, the convolution product in Eq. ( 100) converts, in the real domain, into a simple product of a periodic factor µ k,n (r) resulting from the reciprocal domain sampling (term between brackets) by the phase factor e ik•r introduced according to the well-known shift theorem [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF]. We conclude that, the inverse Fourier transform of Eq. (100) naturally displays the Bloch wave behavior.

F -1  δ  k ′ -k  = e ik•r

Conclusion

The study of wave propagation in periodic media is a thrilling subject that has drawn many attentions since pioneer works a century ago about electronic wave functions and quantum mechanics. We have witnessed over the last decades a renewed interest with the emergence of man-made periodic media: photonic and phononic materials together with the availability of powerful computation machines required to accompany the development of numerical methods.

Throughout this paper, we have presented how classical methods have been developed both in the real and Fourier domains. We have pointed out the two domains correspondence as reported by Brillouin and, on the basis of his perception; we have introduced mathematical tools originating from signal processing to highlight the two domains correspondence. These tools, developed in the framework of Fourier analysis and distribution function theory, permit to establish a well adapted approach to investigate wave propagation in periodic media.

Taking advantage of the two domains correspondence, the Bloch theorem appears to be a consequence of general considerations concerning the Brillouin periodicity of the band structure together with Fourier analysis theorems. The Brillouin periodicity implies an inherent ''shifted discretization'' of the eigenmode wavefunction expressions in the Fourier domain. Then, the discretization in the Fourier domain introduces the periodic part of the Bloch theorem, while the k-shifting is responsible of the existence of the phase factor.

Also, the approach permits to identify the two domain correspondence at different levels: expressions for Fourier transform pairs have been introduced for the periodic parameters as a whole or separately for the crystal lattice and unit cell expression and so, the contributions of the lattice structure are clearly dissociated from those related to the unit cell. Finally, we end with a formalism which is straightforward whatever the complexity of problem is: one or three dimensions, scalar or vector fields, in materials with or without anisotropy, as illustrated by the different analyzed examples.

Ĥψ (r + R) = Eψ (r + R) ⇔ Ĥ  TR ψ (r)  = E  TR ψ (r)  . (B.1)
Now, applying the translation operators on both sides of the Schrödinger equation gives us:

TR  Ĥψ (r)  = TR {Eψ (r)} . (B.2)
The brackets used are not mandatory they have been used in order to clarifying the different situations. The right-hand sides of the two previous equations are identical since E is a constant. So, identifying the two left-hand sides, we get:

Ĥ TR ψ (r) = TR Ĥψ (r) ⇔  Ĥ TR -TR Ĥ = 0; ψ (r) ̸ = 0. (B.3)
In the language of operators we say that the translational operator TR commutes with the Hamiltonian operator Ĥ.

The opposite statement is also true. Indeed, starting this time by the assumption that the commutation of the operators TR and Ĥ holds. We can write:

Ĥ  TR ψ (r)  = TR  Ĥψ (r)  . Since we assume that ψ(r) is an eigensolution of the Hamiltonian the second term reduces to : TR

 Ĥψ  = E  TR ψ 
. Thus the commutation property reduces to:

Ĥ  TR ψ (r)  = E  TR ψ (r)  ⇔ Ĥψ (r + R) = Eψ (r + R) . (B.4)
We conclude that the commutation property of translation and Hamiltonian operators is strictly equivalent to the property that ψ(r) and its translation TR ψ (r) ≡ ψ (r + R) are eigensolution of the Hamiltonian with the same eigenvalue E.

D.1.1. Matrix element of the operator

One defines the integral in left-hand side of Eq. (D.3)) as the matrix elementA ij :

A ij ≡  d 3 rϕ * i (x, y, z) Ta 1 ϕ j (x,
A ij =  k b k δ ik = b i . (D.5)
Eq. (D.5) tells us that the matrix element A ij is numerically equal to the weighting coefficient b i of the basis function ϕ i (x, y, z) in the expansion (D.2) for the image ϕ j (x + a, y, z) corresponding to the basis function ϕ j (x, y, z) after its translation by the operator Ta 1 . Thus Eq. (D.2) writes:

Ta 1 ϕ j (x, y, z) ≡ ϕ j (x + a 1 , y, z) =  i A ij ϕ i (x, y, z). (D.2 ′ ) D.1.

Matrix representation of the operator

Proceeding in the same way for the images of the other basic functions one obtains equation lines similar to (D.2 ′ ). In other words, the Eq. (D.2 ′′ ) stands for the matrix representation 12

A =  A ij  of the operator Ta 1 :          ϕ 1 (x + a, y, z) ϕ 2 (x + a, y, z) . . . ϕ j (x + a, y, z) . . . ϕ n (x + a, y, z)          =          A 11 A 21 • • • A j1 • • • A n1 A 21 A 22 • • • A j2 • • • A n2 . . . . . . . . . . . . . . . . . . A 1j A 2j • • • A jj • • • A nj . . . . . . . . . . . . . . . . . . A 1n A 2n • • • A jn • • • A nn                   ϕ 1 (x, y, z) ϕ 2 (x, y, z) . . . ϕ j (x, y, z) . . . ϕ n (x, y, z)          . (D.2 ′′ )
D.2. Concluding remarks 13(i) Using the terminology or group theory this square matrix A = 

A ij
 is known as a ''Representation'' of the element of the translation group corresponding to the translation vector a 1 .

(ii) The fact that following a translation the function ϕ j (r + a 1 ) acquires contributions from of all other basis functions ϕ i (r); i ̸ = j, is referred to as its ''contamination'' by these basis functions. The relative weighting factors of these contaminations are given by the constants A ij .

(iii) The transpose of the representation matrix defined in (D.2 ′′ ) expresses the weighting coefficients d j of the expansion of the image ψ(x+a 1 , y, z) =  j d j ϕ j (x, y, z) of an arbitrary function ψ(x, y, z) =  j c j ϕ j (x, y, z) in terms of the weighting coefficients c j of the latter in the {ϕ i (x)} basis. So, this expression writes:

d i =  j A ij c
 ϕ ′ i (x)
 in term of the old ones {ϕ i (x)}.

D.3. Matrix representation of the translation operator and eigensolutions

(a) case of translational symmetry element T a 1 where, the translation vector is a single lattice basis vector a 1 The step in determining the eigensolutions consisting of diagonalization of the representation matrix A is termed ''finding the irreducible representation'' A ′ . Where, the prime symbol labels the irreducible representations. Obviously, A ′ is an equivalent representation of the operator in the eigenfunction basis set

 ϕ ′ i (r)  .
Since symmetry translations constitute an Abelian group, the matrices are diagonalizable [START_REF] Laloë | Mécanique quantique[END_REF][START_REF] Bir | Symetry and Strain-Induced Effects in Semiconductors[END_REF][START_REF] Kireev | Semiconductor Physics[END_REF][START_REF] Hümmer | Group Theory in Semiconductor Optics[END_REF]. We say that the ''irreducible representation'' is one-dimensional 14 referring to the fact that its diagonal terms are formed by only one constant A ′ ii (i.e. a matrix of dimensions 1 × 1) A ′ ij ⇒ A ′ ij δ ij where, δ ij is the Kronecker-delta.

(b) case of the translational symmetry element T n 1 a 1 (n 1 successive basis vectors a 1 )

Using an iterative method, beginning with: ϕ j (x + 2a 1 , y, z) =  i A ij ϕ i (x + a 1 , y, z), it is easy to show that the representation of the translation symmetry element Tn 1 a 1 corresponds the matrix  A ij  n 1 . So, in the irreducible representation basis this matrix is also diagonal and writes:

 A ′ ij δ ij  n 1 . (D.6) (

c) Expression of the Eigenvalues and the Born-von Karman boundary conditions

In order to determine the eigensolutions, Bloch relies on thermodynamic principles which establish the equivalence between the physical properties pertaining to the entire crystal or to its restriction to a subsystem of parallelepiped shape sustained by a large number of unit cells, provided that the dimensions of this parallelepiped are sufficiently large. We will refer to this representation as the virtual box concept. This concept rely on the fact that, ''the general form of wavefunctions depends on the periodicity and has nothing to do with the artificially introduced boundaries''. Then, he justifies the use of Born-von Karman cyclic boundary conditions rather than the trivial boundary conditions imposed by the continuity of the wavefunction with the vanishing electronic wavefunction outside the crystal. The reason is because the former leads to propagating wave solutions (which we are dealing with) while the later imposes standing waves solutions.

The cyclic Born-von Karman boundary conditions impose the periodicities of N 1 a 1 , N 2 a 2 and N 3 a 3 respectively in the three-directions of a 1 , a 2 and a 3 . Where, N 1 , N 2 and N 3 stand for the number of unit cells constituting the crystal in the three directions. For instance, in the direction a 1 considered in this subsection this boundary condition reads:

ϕ ′ j (x + N 1 a 1 , y, z) = ϕ ′ j (x, y, z). (D.7)
Thus substituting n 1 by N 1 in the irreducible expression form of Eq. (D.6) leads to the identity matrix:

 A ′ ij δ ij  N 1 = I. (D.8) Now, since  A ′ ij  is a diagonal matrix, we obtain A ′ ij N 1 δ ij = A ′ jj N 1 = 1.
Thus, the diagonal terms and so the eigenvalues are the roots of unity. So, the representation of the translation operator Ta 1 is a diagonal matrix whose elements are the set of eigenvalues:

A ′ jj = e i 2π m 1 N 1 ; - N 1 2 < m 1 ≤ N 1 2 . (D.9)
Note that the interval -N 1 /2 < m 1 ≤ N 1 /2 has been chosen in order to confine the principal value of the complex exponential argument in the interval ]-π , π].

Finally for a translation operator Tn 1 a 1 relative to a translation vector n 1 a 1 for any integer n 1 , the matrix element writes:

A ′ jj = e i 2π m 1 N 1 n 1 ; - N 1 2 < m 1 ≤ N 1 2 . (D.10) (d) Bloch theorem Since  A ′ ij 
is diagonal, the diagonal terms given in Eq. (D.10) correspond to the eigenvalues of Tn 1 a 1 . Thus we can write:

ϕ ′ j (x + n 1 a 1 , y, z) = e i 2π m 1,j N 1 n 1 ϕ ′ j (x, y, z); - N 1 2 < m 1,j ≤ N 1 2
and n 1 ∈ Z .

(D.11)

The added second index j to m 1,j recalls the mutual correspondence between eigenvalues A ′ jj and eigenfunctions ϕ ′ j (x, y, z).

Doing the same for the two other directions gives us:

ϕ ′ j (x, y + n 2 a 2 , z) = e i 2π m 2,j N 2 n 2 ϕ ′ j (x, y, z); - N 1 2 < m 2,j ≤ N 1 2 and n 2 ∈ Z (D.12) ϕ ′ j (x, y, z + n 3 a 3 ) = e i 2π m 3,j N 3 n 3 ϕ ′ j (x, y, z); - N 1 2 < m 3,j ≤ N 1 2
and n 3 ∈ Z .

(D.13) Now, regrouping (D.11)-(D.13), we obtain 15 for the general case of a 3-D translation vector corresponding to any arbitrary lattice vector R = n 1 a 1 + n 2 a 2 + n 3 a 3 :

ϕ ′ j (x + n 1 a 1 , y + n 2 a 2 , z + n 3 a 3 ) = e i  2π m 1,j N 1 n 1 + 2π m 2,j N 2 n 2 + 2π m 3,j N 3 n 3  ϕ ′ j (x, y, z). (D.14)
This last relation is the general form of the solution and thus it stands for the Bloch-Floquet theorem. In Section 2.1.3.3, using coordinate vector notations we give a more explicit form of this relation. This requires the introduction of the reciprocal lattice concept.

Appendix E. 3-D Fourier domain method to determine the Bloch theorem

Consider a crystal consisting of N 1 × N 2 × N 3 unit cells along the tree direction axes as defined by its lattice basis vectors a 1 , a 1 and a 3 not necessarily orthogonal. Defining the unit vector along these axis ℓ i = a i a i we can write any position vector as: r = ℓ 1 ℓ 1 + ℓ 2 ℓ 2 + ℓ 3 ℓ 3 . Where ℓ 1 , ℓ 2 and ℓ 3 are the coordinates along the oblique axes. The 3-D Born-von Karman periodic boundary conditions write:

 ψ E (ℓ 1 + N 1 a 1 , ℓ 2 , ℓ 3 ) = ψ E (ℓ 1 , ℓ 2 , ℓ 3 ) ψ E (ℓ 1 , ℓ 2 + N 2 a 2 , ℓ 3 ) = ψ E (ℓ 1 , ℓ 2 , ℓ 3 ) ψ E (ℓ 1 , ℓ 2 , ℓ 3 + N 3 a 3 ) = ψ E (ℓ 1 , ℓ 2 , ℓ 3 )
or in vector notation ψ E (r + N i a i ) = ψ E (r).

(E.1)

Using the same sign conventions introduced in Section 2.2.1, and assuming that the Fourier harmonic components ϕ E,p 1 ,p 2 ,p 3 are known for all integers p 1 , p 2 and p 3 , the 3-D Fourier series for the artificially introduced periodicity writes: where, for shortness we have omitted the k indices p i and k stands for the set including all the wavevectors k p 1 ,p 2 ,p 3 . Note that, in the limit where the N i tend toward infinity the wavevector k recovers its continuous nature and the Fourier series (E.3) turns out to be an inverse Fourier transform.

ψ E (ℓ 1 , ℓ 2 , ℓ 3 ) =  p i ϕ E,
In the same way, the potential U (r) being periodic with period a, it is also expressed in term of a Fourier series:

U (r) =  G∈G U G e -iG•r (E.4)
where, G stands for the set of all the reciprocal lattice vectors G and U G are the Fourier components of the periodic potential. The later are given by: 15 A more rigorous derivation using the concept of the matrix ''direct product'' may be found in Refs. [START_REF] Kireev | Semiconductor Physics[END_REF][START_REF] Hümmer | Group Theory in Semiconductor Optics[END_REF].

U G = 1 V DL  V DL
Now, we change of dummy variable k ′ to k = k ′ -G in the second summation in Eq. (E.7). Note that we have anticipated this change by adjoining a prime to this dummy variable in going from (E.6) to (E.7). So, Eq. (E.7) writes:

 k∈k  ϕ E,k  h2 |k| 2 2m -E  +  G∈G U G ϕ E,k+G  e ik•r  = 0. (E.8)
This sum must be zero for any coordinate vector r. This is only possible if each term of the sum vanishes independently. Thus, in the Fourier domain, the expression of the Schrödinger equation reduces to an infinite set of linear equations involving the set of variables ϕ E,(k+G) . This system of linear equations constitutes the three dimension central equation:

ϕ E,k  h2 |k| 2 2m -E  +  G∈G U G ϕ E,k+G = 0.
(E.9)

The three-dimensional Bloch theorem.

The terms of the summation in the central equation (E.9) reveal that the real space periodicity couple all those, and only those, Fourier spectral components which differ from one another by any reciprocal lattice vector G ∈ G. Thus, once we are handling a given wave vector k, we are, in fact, dealing with a subsystem of linear equations involving all components having a wave vector of the form (k + G) for all G ∈ G (excluding all other k vectors). This justifies the common ''k-labeling'' of eigensolutions ψ E,k (r) referring to the k-subsystem. So, for a fixed k vector, only the terms of the form (k + G) will be retained in Eq. (E.3). This later may be rewritten as: ψ E,k (r + R) = e ik•R ψ E,k (r) .

(E.13)

It is instructive to notice the correspondence between, (E.13) and [START_REF] Linton | Resonant effects in scattering by periodic arrays[END_REF]. One can observe that the Fourier domain method naturally leads to the Bloch function in the form of Eq. (E.12) and then deduces the other form whereas; the group theory ends with the second form (E.13) and then the first is determined.

Appendix F. Complex expression of wavevector in the bandgap above the cutoff frequency ω > ω c

We are looking for the solutions of the equation: sin

 ka 2  = ω ω c > 1.
(F.1)

In the case ω > ω c , real solutions do not exist. So let the complex solutions be of form: (F.7)

Finally, the principal values of k write:

⇒ k = ± π a ± i 2 a cosh -1  ω ω c
 .

(F.8)

Appendix G. Splitting of the dispersion relation in two frequency bands ω + and ω - Substituting the displacement expressions from Eq. ( 43) in Eq. ( 44) and eliminating the exponential terms e i(ωt-nka) we obtain: (G.5)

 -ω 2 m 1 A 1 = η  A 2 e -ika
The solutions may be written as:

ω 2 ± = η    m 1 + m 2 m 1 m 2  ±   m 1 + m 2 m 1 m 2  2 - 4 sin 2 (ka/2) m 1 m 2   (G.6)
or:

ω 2 ± = η m 1 m 2  m 1 + m 2 ±  (m 1 + m 2 ) 2 -4m 1 m 2 sin 2 (ka/2)
 .

(G.7)

These two solutions give rise to two allowed branches.

(a) Qualitative considerations Whenever the phase ak is different from an integer q times 2π , the phase factors e jnak summation will add destructively to zero (like in a destructive interference situation). However, when ak = 2π q, the phase factors add up constructively and the infinite summation tends toward infinity. So, we see that we are in the presence of periodic Dirac-delta functions spaced by k = 2π /a i.e. it correspond to a Dirac-comb with a period equal to the reciprocal lattice basis vector b = 2π /a.  .

(I.3)

• At the pole k = 0: we can specify the value of the limit in (I.3), using l'Hopital's rule. This allows us to lift the indeterminate form 0/0 and find out the limits. We obtain: (I.6)

F  +∞  n=-∞ δ (x -na)  = Lim N→+∞  (2N + 1) a 2 cos [(2N + 1) π ]

  which define the reciprocal lattice G = m 1 b 1 +m 2 b 2 +m 3 b 3 . The definition of the reciprocal lattice basis vectors implies the condition b i • a j = 2π δ ij . Where, δ ij is the Kronecker-delta.

Fig. 1 .

 1 Fig. 1. Boundaries of the first Brillouin zones. The boundary planes corresponding to the reciprocal basis vectors b 1 and -b 1 are represented by their respective cross sections and ′ with the plane containing b 1 and arbitrary wave vectors k and k ′ having their end points on these boundaries.

Fig. 2 .

 2 Fig. 2. Sinusoidal curves standing for the Fourier components u kn (x, t) = Ae i(ωt-kn x) corresponding to three values of k n = k + 2nπ /a; n = -1, 0 and 1, for the special case a = 0.1 units of length i.e. b ≡ 2π /a = 20π and k = b/12. The points of intersection at x m = ma = 0.1m are marked with bold dots.

Fig. 3 .

 3 Fig.3. Layout of the array. The unit cells are defined in such a manner that the point masses m 1 are located on the unit cell boundaries (i.e. at x n = na) and m 2 at the middle of the periods (i.e. at x n+1/2 = (n + 1/2) a).

  clearness, the two domains correspondence is schematically illustrated in Fig.4. This correspondence may be broken down into three levels: the unit-cell expression m(x) ⇔ m(k), the direct and reciprocal lattice n δ (xna) ⇔ b  m δ (knb)as well as the two closed-form expressions m(x) ⇔ m(k) as a whole.

Fig. 4 .

 4 Fig. 4. Schematic representation of a periodic parameter using the replicating property in the direct domain and its corresponding representation in the reciprocal domain. (a) Two domain representations of a unit cell. (b) Two domain representations of the lattice. (c) Two domain representations of the periodic parameter.

  The 3-D direct and reciprocal lattice representations, with their respective lattice vectors: R = n 1 a 1 +n 2 a 2 +n 3 a 3 ⇔ G = m 1 b 1 +m 2 b 2 +m 3 b 3 constructed respectively on the direct and reciprocal lattice basis vectors a i ⇔ b i . The sets R and G include respectively all the direct and reciprocal lattice vectors i.e. all triplets of integers (n 1 , n 2 , n 3 ) and (m 1 , m 2 , m 3 ); v RL = b 1 • (b 2 × b 3 ) = (2π ) 3 /v DL stands for the volume of the reciprocal lattice unit cell and v DL is the volume of the direct lattice unit cell.

5. 3 . 1 .

 31 Phononic crystals with non-piezoelectric constituent materials 5.3.1.1. Vector field equation in the real domain. Starting from dynamic equation of motion for the displacements u i , (i = 1, 2, 3) [33]:

2 , 3

 23 and i = 1, 2, 3 = e lij (r) ; . . . k = 4 and i = 1, 2, 3 = e ikl (r) ; . . . k = 1, 2, 3 and i = 4 = ε S il (r) ; . . . k = 4 and i = 4 and ρi (r) = ρ (r) ; . . . i = 1, 2, 3 = 0; . . . i = 4.

  Substituting (D.4) in (D.3) and taking account of the orthonormal basis condition (D.1), Eq. (D.3) writes:

  j with the usual subindex notations. (iv) Different basis function sets leads to different matrix representations. The change of basis functions from the initial set {ϕ i (x)} to the new one  ϕ ′ i (x)  will introduce a new representation matrix A ′ . The determination of the later in term of the initial representation matrix A follows the usual matrix formula: A ′ = B -1 AB. Where, the matrix B expresses the new basis functions

p 1 ,p 2 ,p 3 , e i  p 1 2π N 1 a 1 ℓ 1 +p 2 2π N 2 a 2 ℓ 2 +p 3 2π N 3 a 3 ℓ 3 .(E. 2 )

 121122332 In a similar way followed in Section 2.1.3.3 in order to introduce the reciprocal basis vectors, we can identify a discrete form for the wavevectork p 1 ,p 2 ,p 3 = p 1 k N 1 +p 2 k N 2 +p 3 k N 3 . Where, k N i = b i /N i .This discretization results from the artificially limited crystal dimensions. Now, using this definition of k p 1 ,p 2 ,p 3 , Eq. (E.2) reduces to the Fourier series in vector notation:ψ E (r) =  k∈k ϕ E,k e i(k•r) (E.3)

dr 3 U

 3 (r) e iG•r . (E.5) Now, in order to determine the Fourier componentsψ E,m , we substitute Eqs. (E.3) and (E.4) in the Schrödinger equation. This gives us: ϕ E,k e ik•r +  k ′  G∈G U G ϕ E,k ′ e i(k ′ -G)•r = 0. (E.7)

ψ

  E,k (r) =  G∈G ϕ E,(k+G) e i(k+G)•r . brackets has the form of a Fourier series. So, it is a 3-D periodic function. Denoting this function µ k (r) equation (E.11) writes: ψ E,k (r) = µ k (r) e ik•r . (E.12) This is one of the expression the three dimensional Bloch theorem. Substituting r by r + R, we obtain the second expression:

ka 2 =(F. 2 ) 2  2 .(F. 5 )

 22225 A + iB.Then using the identity:sin (A + iB) = sin A cosh B + i cos A sinh B, (F.3) we can write: sin  ka = sin A cosh B + i cos A sinh B = ω ω c > 1.(F.4)Identifying the real and imaginary part in (F.4):⇒  cos A sinh B = 0 sin A cosh B = sin kaRejecting the trivial solution: sinh B = 0; sin A = sin ka/2, Eq. (F.5) gives us:

( b ) 2 )

 b2 Quantitative determination of the strength of one of the Dirac-delta functions constituting the Dirac-combWe can have a more specific determination of the Dirac-comb using limiting geometrical figures approaching the central Dirac-function of the comb. Let us begin by approximating the right-hand infinite summation in (I.1) as the limit of the following finite sum:And so, we have to determine the sum of a geometric series of 2N + 1 terms with Lim N→+∞  e -jNak  as first term, and e jak as common ratio. The sum of this series writes: e j(2N+1)ak -1 e jak -1  = Lim N→+∞  sin [(2N + 1) ak/2] sin [ak/2]

  pole at k = 0 tends toward infinity as: (2N + 1).• Determining the Zeros of the central Dirac-delta function: Since we have restricted the summation to N, the values of k are discrete. Thus, the first zeros, which border on the central pole, correspond to those of the numerator of Eq. (I.3).That is: (2N + 1) ak/2 = ±π. Finally, we obtain: identified the reciprocal lattice basis vector b = 2π /a.• Geometrical figure approaching the central Dirac-delta function: Since N → +∞ the two zeros given in Eq. (I.5) occur at k values which are vanishingly small on both sides of an infinitely pole of height (2N + 1). So we can approach the geometrical aspects of the central Dirac-delta function constituting the Dirac-comb by a triangle of height (2N + 1) and base = 2b/ (2N + 1). Its area equals b and gives the strength of the Dirac-delta functions forming the comb. Finally we can write the Fourier transform: qb) .

  /2 -2A 1 + A 2 e ika/2  -ω 2 m 2 A 2 = η  A 1 e -ika/2 -2A 2 + A 1 e ika/2  .  2ηω 2 m 1  A 1 -2η cos (ka/2) A 2 = 0 -2η cos (ka/2) A 1 +  2ηω 2 m 2  A 2 = 0. 4 m 1 m 2 -2ηω 2 (m 1 + m 2 ) + 4η 2  1 -cos 2 (ka/2)

												(G.1)
	Rearranging the terms we obtain:
												(G.2)
	Solutions exist if:
	thus:	    -2η cos (ka/2)  2η -ω 2 m 1 	 -2η cos (ka/2)     = 0  2η -ω 2 m 2			(G.3)
		or : ω 4 -	 2η m 1	+	2η m 2		ω 2 +	(2η) 2 m 1 m 2	 sin 2 (ka/2)		 = 0. = 0	(G.4)

ω

The origin of the term reciprocal lattice resides in the representation of the dot product b i • a j = δ ij by the product of two 3 × 3 matrices B and A.Where, assuming a rectangular coordinate system, the ith row of B contains the coordinates of the vector b i normalized with respect to

2π while the jth column of A contains the coordinates of the vector a j . It is obvious that, according to the dot product b i • a j = δ ij , the matrices A an B are mutually the reciprocal (inverse) of each other[START_REF] Brillouin | Propagation des ondes dans les milieux périodiques[END_REF].

This is the origin of the so-called analog electronic as opposed to the digital electronic the later concerns discrete signal analysis.

Fourier series or inverse Fourier transform depending on the discrete or continuous nature of the spectrum.

More accurately, 3 bands if we take into account the longitudinal and transverse vibrations: two for the transverse and one for the longitudinal vibrations.

Note the unusual order notation of the subscript in the matrix element which rises by way of its definition in Eq. (D.4).

Floquet, outside the framework of group theory, has exhaustively studied all these situations. Once again, this shows equivalence between the two approaches.

Note that, in general, for other symmetry space transformations this diagonalization is not always possible. In this case one reduces the matrix as far as possible and obtains at the end a block diagonal matrix. These blocks are known as the irreducible representation of order >1. Since we are interested by translation operators, we will not treat this case here.

Acknowledgment

This work has been carried out within project no. ANR-09-NANO-004 funded by ANR (P3N2009).

Appendix A

The set of symmetry translation constitutes a group. Indeed, it verifies the four following criteria:

1. The ''composition law'' associated to the set elements (here, the multiplication operation of two translations) is an internal one. i.e. the resulting product is also an element belonging to the set. (It is clear that the product of two translations is also a translation

2. The multiplication is associative:

The existence of an identity element:

The existence of an inverse for each element.

Appendix B. Invariance of the Hamiltonian under the change of variables r  → (r + R), is equivalent to the commutation property of the translation operator with the Hamiltonian

Here we present a heuristic demonstration of this statement, a more rigorous demonstration may be found in any quantum mechanics or solid-state textbooks for example [START_REF] Laloë | Mécanique quantique[END_REF][START_REF] Rosencher | [END_REF][START_REF] Kireev | Semiconductor Physics[END_REF][START_REF] Hümmer | Group Theory in Semiconductor Optics[END_REF][START_REF] Joannopoulos | Photonic Crystals. Molding the Low of Light[END_REF][START_REF] Landau | Physique Theorique: Tome 3, Mécanique quantique[END_REF].

Applying the same line of reasoning as Floquet (see Section 2.1.1) the assumption that the Hamiltonian is invariant under the change of variables r  → (r + R) imposes that: if ψ(r) is a solution of the Schrödinger equation then ψ(r + R) = TR ψ(r)

is also a solution with the same eigenvalue E. In terms of translational operator, this writes:

Appendix C. Corollary: property of two commutating operators ''they own a common set of eigenvectors with the same eigenvalue''

This may be easily verified in the case of a non degenerate state. Since in this case, the fact that the two functions ψ (r) and TR ψ (r) are solutions of Schrödinger equation belonging to the same eigenvalue E, means that they represent the same eigenfunction and thus differ from one another by at most a constant factor. In other words: TR ψ (r) = t R ψ (r) where t R is a constant which depends on the translation vector R. The opposite statement is also verified. That is, if ψ (r) is a solution of the translational operator TR ψ (r) = t R ψ (r) then TR ψ (r) is proportional to ψ (r). Now, because we assume that ψ (r) is a solution of Schrödinger equation, then TR ψ (r) will also be a solution with the same eigenvalue.

In case of degenerate states, the problem is somewhat different but, we can always find linear combinations that are simultaneously eigenfunctions of Schrödinger equation and translation operator [START_REF] Laloë | Mécanique quantique[END_REF][START_REF] Kireev | Semiconductor Physics[END_REF].

Appendix D. An illustrative example of group theory ''The translational symmetry group: representation and eigensolutions''

Our aim in this appendix is to introduce the representation of the translational symmetry operation following the way paved by Bloch to find out its eigenvalues. For sake of simplicity, and without loss of generality, Bloch restricts his study to the case of crystals with rectangular lattices with their basis vectors aligned along the coordinate axes. This assumption does not have any incidence on the general conclusions he is seeking for. So, the outcomes he obtained are still valid for any other crystallographic axes. More general demonstrations can be found in Refs. [START_REF] Laloë | Mécanique quantique[END_REF][START_REF] Kireev | Semiconductor Physics[END_REF][START_REF] Hümmer | Group Theory in Semiconductor Optics[END_REF] which confirm Bloch assumptions. This assumption enables one to use scalar coordinate notations.

D.1. Matrix representation of a space transformation operator

In group theory, the convention is to restrict the use of the term ''representation'' to the representation of the group by a set of square matrices. Each matrix represents an element of the group. 11 Here, we use translation operator as an illustrative example.

Let {ϕ i (x)} stands for a complete set of linearly independent orthonormal basis functions i.e.:

Consider the image of one of these basic functions after a space transformation. For instance the translation: Ta 1 ϕ j (x, y, z) ≡ ϕ j (x + a 1 , y, z). In an arbitrary chosen basis, this image is no more a basis function. It can be expanded as a linear combination:

Multiplying the both sides of Eq. (D.2) by ϕ * i (x, y, z) and integrating all over the space we obtain: 

11 Also, the product of matrices represents the symmetry elements product defined in Section 2.1.3.1.

Finally inverting the relation (G.5) we obtain expression [START_REF] Linton | One-and two-dimensional lattice sums for the three-dimensional Helmholtz equation[END_REF], this equation will be used to investigate the forbidden bands.

Appendix H. List of the most important Fourier transform formulas used in the framework of our two domains correspondence formalism [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF][START_REF] Laverriere | Electrodynamique classique[END_REF] Whenever needed, the Fourier transforms are separate according to the two conventions used in this paper.

The symbol ⇔ refers to Fourier transform pairs. Consider a Dirac-comb of period a:

Using the Fourier transform pair δ (x ± a) ⇔ e ∓ika The Fourier transform of the Dirac-comb writes:

Our aim is to write down the right-hand side of expression (I.1) in terms of a Dirac-comb in the reciprocal space without using the Poisson formulas.

Appendix J. Correspondence between the Dirac-comb concept and the Poisson summation formulas

In this appendix, we outline the correspondence by establishing that Poisson summation formulas may be viewed as an immediate result of FT property of Dirac-comb which is the basis of our approach. Obviously, the opposite statement is also true. Thus the two approaches are equivalent in the mathematical sense.

Poisson summation formulas as an immediate result of the property of the Fourier transform of a Dirac-comb used in our approach (three-dimensional space)

We apply the same notations used in the paper and the Fourier transform convention adopted for photonic crystal. That is:

where, we distinguish a function from its Fourier transform by its argument. If the argument contains real space (reciprocal) coordinates it denotes the real space (Fourier transform) function.

(a) Special case of the Dirac-comb/exponential pair Starting from the Fourier transform pair:

Comparing D (r) to the inverse Fourier transform of D (k), obtained by carrying out the inverse transform of individual Dirac-delta functions inside the summation (linearity of FT), gives:

Taking the Fourier Transform on both sides gives its reciprocal domain expression:

Following the same procedure but this time beginning with the Fourier transform pair:

Comparing U (r) to the inverse Fourier transform of U (k) performed by using inverse transform of individual Dirac-delta functions inside the summation gives:

Taking the Fourier Transform on both sides gives its reciprocal domain expression:

where, we applied the shift theorem in order to determine the left-hand side FT.