High-frequency homogenization in periodic media with imperfect interfaces

In this work, the concept of high-frequency homogenization is extended to the case of one-dimensional periodic media with imperfect interfaces of the spring-mass type. In other words, when considering the propagation of elastic waves in such media, displacement and stress discontinuities are allowed across the borders of the periodic cell. As is customary in high-frequency homogenization, the homogenization is carried out about the periodic and antiperiodic solutions corresponding to the edges of the Brillouin zone. Asymptotic approximations are provided for both the higher branches of the dispersion diagram (second-order) and the resulting wave field (leading-order). The special case of two branches of the dispersion diagram intersecting with a non-zero slope at an edge of the Brillouin zone (occurrence of a so-called Dirac point) is also considered in detail, resulting in an approximation of the dispersion diagram (first-order) and the wave field (zeroth-order) near these points. Finally, a uniform approximation valid for both Dirac and non-Dirac points is provided. Numerical comparisons are made with the exact solutions obtained by the Bloch–Floquet approach for the particular examples of monolayered and bilayered materials. In these two cases, convergence measurements are carried out to validate the approach, and we show that the uniform approximation remains a very good approximation even far from the edges of the Brillouin zone.

RCA, 0000-0001-9848-3482; BL, 0000-0002-5648-2737 In this work, the concept of high-frequency homogenization is extended to the case of onedimensional periodic media with imperfect interfaces of the spring-mass type. In other words, when considering the propagation of elastic waves in such media, displacement and stress discontinuities are allowed across the borders of the periodic cell. As is customary in high-frequency homogenization, the homogenization is carried out about the periodic and antiperiodic solutions corresponding to the edges of the Brillouin zone. Asymptotic approximations are provided for both the higher branches of the dispersion diagram (second-order) and the resulting wave field (leading-order). The special case of two branches of the dispersion diagram intersecting with a non-zero slope at an edge of the Brillouin zone (occurrence of a so-called Dirac point) is also considered in detail, resulting in an approximation of the dispersion diagram (first-order) and the wave field (zeroth-order) near these points. Finally, a uniform approximation valid for both Dirac and non-Dirac points is provided. Numerical comparisons are made with the exact solutions obtained by the Bloch-Floquet approach for the particular examples of monolayered and bilayered materials. In these two cases, convergence measurements are carried out to validate the approach, and we show that the uniform approximation remains a very good approximation even far from the edges of the Brillouin zone.

Introduction
Classically, dynamic homogenization is understood as a low-frequency approximation to wave propagation in heterogeneous media such as laminates, composites, or more generally any microstructured media. It consists in approximating such media by effective homogeneous media with specific properties. Homogenization is the mathematical process that allows one to find such properties. Much work on this topic has been carried out since the 1970s and it is not the aim of this introduction to be exhaustive in that regard. A particularly successful approach is the twoscale asymptotic expansion method and the notion of slow or fast variables (see e.g. [1][2][3]). Periodic media, of interest in this paper, are dealt with very efficiently by such a method.
In periodic media, waves can propagate at angular frequencies ω that are not necessarily small. The set of wavenumbers k (Bloch wavenumbers) at which waves propagate depends on the angular frequency through dispersion relations. In particular it can be shown that these dispersion relations can be entirely understood using diagrams restricting the Bloch wavenumbers to lie within the Brillouin zone. In one dimension, when the periodicity of the structure is h say, such a Brillouin zone is given by k ∈ [0, π/h]. Typically, the dispersion diagram displays band gaps, i.e. regions in the angular frequency space where waves cannot propagate. There tends to be infinitely many branches of the dispersion diagram, i.e. for a given Bloch wavenumber, one can find an infinite (countable) set of angular frequencies leading to propagating waves. Again, a lot has been written about this, and we do not aim to give an exhaustive literature review on this point, though we can refer the interested reader to [4] for example. See also [5] for a discussion of the band gaps in periodic materials from a physics point of view.
The idea of high-frequency homogenization is to approximate how the dispersion relation (and hence the media) will behave for angular frequencies ω that are close to the angular frequencies ω 0 corresponding to an edge of the Brillouin zone on the dispersion diagram. In [6], a work that has largely inspired the present paper, Craster et al. applied a two-scale asymptotic expansion method in order to achieve this for perfect interfaces. Adopting the terminology of [7], we will refer to the homogenization near the left edge of the Brillouin zone k ≈ 0 as finite frequency low wavenumber (FFLW), while the homogenization near the right edge (k ≈ π/h) will be referred to as finite frequency finite wavenumber (FFFW). Upon introducingk as (FFLW):k = k and (FFFW):k = π h − k, (1.1) the end result of the high-frequency homogenization technique is an approximation of the type where L is a macroscopic characteristic length of the material and the parameter T ∈ R can be determined explicitly. This angular frequency approximation comes together with an associated leading-order approximation to the wave field U h (X) of the form One should note that in [7], the authors generalized the technique to work in any dimension, and also pushed the asymptotic work one order further than [6]. It is also worth mentioning the more recent work [8], in which the technique was developed while including a source term in the initial evolution equation. As shown in [9], high-frequency homogenization of wave equations in the time domain can also be carried out, and the methodology has also been applied in a discrete setting to structural mechanics [10] and elastic lattices [11]. In this work, following the approach of [6], we wish to focus on extending this high-frequency homogenization technique to one-dimensional periodic media that have an imperfect interface at the edges of the periodic cell, which, to our knowledge, has not been considered before.

Problem formulation (a) The physical problem
We consider linear elastic waves propagation at a given angular frequency ω through a periodic medium of periodicity h > 0 and with a macroscopic characteristic length L > 0. We denote the physical space variable X; the density ρ h (X) and the Young's modulus E h (X) of the elastic medium are assumed to be h-periodic piecewise smooth, L ∞ and strictly positive. We may assume without loss of generalities that the edges of the periodic cell are located at X n = nh for n ∈ Z, as .  illustrated in figure 1a. We further assume that the interfaces across the edges of the periodic cells are imperfect, and of the linear spring-mass type, characterized by some mass and stiffness parameters denoted respectively M and K that are both strictly positive. This results in the following governing equation and jump conditions for the displacement field U h (X): where · X n and · X n are respectively called the jump and mean brackets at the interface X n , and are defined for any function g(X) by Due to the h-periodicity of ρ h and E h , it is possible to write them as ρ h (X) = ρ(X/h) and E h (X) = E(X/h) for some 1-periodic functions ρ and E.
Note 2.1. Throughout this work, we assume that any discontinuity of ρ and E on (0, 1) are of a perfect contact nature. This means that apart from at the interfaces X = X n , we will assume that the displacement U h (X) and the stress E h (X)(dU h /dX) are continuous.
Classically, such wave propagation problem in periodic media can be understood by using the so-called Bloch-Floquet analysis that consists in seeking solutions of the form that propagate at the Bloch wavenumber k, and where U h (X) is h-periodic, that is, U h (X + h) = U h (X). In certain simple cases (see e.g. §6a,b), this leads to an explicit dispersion relation relating ω to k, the graphical representation of which is the so-called dispersion diagram.
One should note that for k = 0, we have U h (X) = U h (X + h), i.e. the solution is periodic, the solution is antiperiodic. Such values of k correspond to the edges of the so-called Brillouin zone. The aim of the present work is to approximate the Bloch-Floquet solutions propagating at wavenumbers that are close to k = 0 (FFLW) or k = π h (FFFW) with finite angular frequency ω = O(1).

(b) Non-dimensionalization
In order to simplify the mathematical notations, we start by non-dimensionalizing the physical problem (2.1). In order to do so, we define the characteristic dimensional density ρ = ρ , Young's modulus E = 1/E −1 and wavespeed c = √ E /ρ , where the average operator · is defined for any function g as These can be used in order to define the following non-dimensional quantities The starred quantities ρ and E are chosen for convenience to be the effective properties of the medium obtained by low-frequency homogenization, implying that α = 1/β = 1, but this choice is somewhat arbitrary. Moreover, we can also show that Using these quantities, (2.1) can be rewritten as the non-dimensional governing equation subject to the jump conditions in the geometry setting of figure 1b. The equations (2.5)-(2.6) constitute our non-dimensional problem. Note that in this non-dimensional setting, the Bloch-Floquet analysis is still valid, and consists in looking for solutions of the form where κ is the non-dimensional Bloch wavenumber and u δ is δ-periodic. Note that this implies that u δ and its derivative u δ satisfy The FFLW case corresponds to κδ ≈ 0 and the FFFW case corresponds to κδ ≈ π . When κδ is exactly 0 (resp. π ), then the solution u δ (x) is δ-periodic (resp. δ-antiperiodic).
(c) Two-scale asymptotic expansion We will now make the assumption that the macroscopic characteristic length L is much bigger than the periodicity h, implying that δ 1. The material parameters α and β will hence vary on a fine scale associated with the rescaled coordinate y = x/δ (see figure 1c for the associated geometrical configuration).
Following the two-scale expansion technique, we further assume that the displacement field will have small scale features described by y, and slow continuous variations described by x. We hence pose the following ansatz for the wave field u δ and the reduced frequency μ: and treat x and y as two independent variables (scale separation), implying that d/dx ↔ ∂/∂x + (1/δ)(∂/∂y). In the FFLW case we will assume that u j is 1-periodic in y, that is u j (x, y) = u j (x, y + 1), while in the FFFW case, we will assume that u j is 1-antiperiodic, that is u j (x, y) = −u j (x, y + 1), see remark 2.2. In what follows, we will treat both cases simultaneously. The non-dimensional problem (2.5)-(2.6) can hence be rewritten as the governing equation (2.10) subject to the jump conditions at y n = n: where for any function g(x, y), the jump and mean brackets are naturally defined for n ∈ Z as g(x, y) y n = g(x, n + ) − g(x, n − ) and g y n = 1 2 (g(x, n + ) + g(x, n − )).
(2.13) Note 2.3. In order for our expansion to be compatible with the assumption regarding potential discontinuities of ρ and E within the unit cell made in note 2.1, we seek the fields such that u 0 , β(∂u 0 /∂y), u j and β(∂u j /∂y + ∂u j−1 /∂x) for j ≥ 1 are continuous functions of y on (0, 1).
Before pushing the asymptotic analysis further, we need to discuss some important properties of the jump and mean brackets. We will start with a very useful property (that can be proved directly), namely that for any two functions f (x, y) and g(x, y), and any n ∈ Z, the following relation is valid: fg y n = f y n g y n + f y n g y n . (2.14) It is equally important to note that if the function subjected to the brackets is either periodic or antiperiodic, the following result holds.
The lemma 2.4 means that as long as the function subjected to either the jump or the mean bracket is 1-periodic or 1-antiperiodic in y, then the mean and jump brackets values are independent of n, and the y n subscript can be dropped. Since β(y) is 1-periodic and u j (x, y) is either 1-periodic or 1-antiperiodic in y, this is the case for all the brackets in the conditions (2.11) and (2.12), we will hence drop the subscript from now on and just use · and · . We will make sure that whenever this notation is used, the function inside the brackets is either periodic or antiperiodic. In particular, from lemma 2.4 and (2.14), we obtain directly the following lemma that will prove very important in what follows.

Lemma 2.5.
For any two functions f (x, y) and g(x, y) that are either 1-periodic or 1-antiperiodic in y, we have fg = f g + f g .
Finally, in order to link the average operator and the jump bracket, the following result will be very useful. Lemma 2.6. For any function g per (x, y) that is 1-periodic and continuous for y ∈ (0, 1), we have ∂g per ∂y = − g per .
Notation 2.7. From now on, in order to efficiently deal with the FFLW (κδ ≈ 0) and the FFFW (κδ ≈ π ) cases simultaneously, we will assume that whenever the symbols ± or ∓ are used, the top sign corresponds to FFLW while the bottom sign corresponds to FFFW.
We are now well equipped to start developing the core theoretical part of the paper.

The case of simple eigenvalues
We will now deploy the two-scale asymptotic procedure in order to derive the approximation (1.2). In order to do so we will need to consider the contributions of (2.10)-(2.12) at the orders δ 0 , δ 1 and δ 2 .

(a) Zeroth-order field
Upon collecting the terms of order δ 0 in (2.10)-(2.12), we obtain the system Upon considering the (Sturm-Liouville) differential operator L:=(−1/α)(d/dy)(β(d/dy)), one can see that the system (3.1) constitutes an eigenvalue problem for L. It is a bit unusual since the eigenvalue also appears in the boundary conditions, but one can show, using the tailored inner product ·,· defined for some functions f (y) and g(y) (either both FFLW or both FFFW) by that this operator is symmetric (i.e. self-adjoint) and non-negative in both the FFLW and FFFW cases, the proof being deferred to appendix A. Therefore, it has a discrete set of (possibly repeated) real positive eigenvalues associated with real eigenfunctions. We denote the square root of such eigenvalues by μ 0 , and we note that eigenfunctions associated to different FFLW (resp. FFFW) eigenvalues are orthogonal for the inner product (3.2). These reduced frequencies μ 0 correspond to the intersection of the dispersion diagram with the left (FFLW) or the right (FFFW) border of the Brillouin zone. From now on we will choose one of these eigenvalues, denote it by μ 0 , and endeavour to approximate the solutions for some parameters (μ,κ) close to (μ 0 , 0), where we defineκ to be allowing us to treat the FFLW and FFFW cases simultaneously. We will also assume that μ 0 is a simple eigenvalue (multiplicity 1). The case of a double eigenvalue will be dealt with in §4. Hence, there is only one eigenfunction that we denoteû 0 (y) and that is either periodic (FFLW) or antiperiodic (FFFW). The associated solution u 0 to (3.1) can therefore be rewritten u 0 (x, y) = U 0 (x)û 0 (y), (3.4) for some function U 0 (x). It is worth mentioning at this stage that when inputting (3.4) into (3.1), we find thatû 0 (y) satisfies both the equation (3.1a) and the jump conditions (3.1b), a fact that will be used throughout the paper. The main aim of what follows is to derive a differential equation with constant coefficients satisfied by U 0 (x). Note that in the case of low-frequency homogenization, the zeroth-order field u 0 (x, y) can be shown to be independent of y; this is one of the main differences between low-and high-frequency homogenization.

(b) First-order field
Let us now set μ 0 to be one of the reduced frequencies found in the previous section, we can collect the terms of order δ 1 in (2.10)-(2.12) to obtain the following system governing the first-order field u 1 : subject to the jump conditions We will show below that there is only one possible value of μ 1 , and that it has to be zero. (i) Proving that μ 1 = 0 The terms in αμ 2 0 u 0 u 1 cancel out, and using the fact that Note that due to the form (3.4) of u 0 , we have u 0 (∂ 2 u 0 /∂x∂y) − (∂u 0 /∂x)(∂u 0 /∂y) = 0, so the first bracket in the right-hand side (RHS) of the equation above is actually zero. Now, using note 2.3 and lemma 2.6, (3.7) becomes When perfect interfaces are considered, the jump bracket term in (3.8) is automatically zero, a fact that is used in [6] to conclude that μ 1 = 0. Such reasoning cannot be used directly in our case. Instead, make use of lemma 2.5 to rewrite (3.8) as where the jump conditions (3.1a) and (3.6) have been used. Many terms cancel out and we get because m and α are strictly positive, and, in the representation (3.4), U 0 cannot be identically zero andû 0 is a real function. This result is very important and lies at the heart of the success of the high-frequency homogenization method.
(ii) An expression for u 1 (x, y) We can now simplify the equation (3.5) governing u 1 to Upon noting that, becauseû 0 is a solution to (3.1a), the field −yU 0 (x)û 0 (y) is a particular solution to (3.11), and that the differential operator applied to u 1 is exactly the same as that of (3.1a), we can conclude that u 1 can be written as for some function U 1 (x) that will be shown not to play any role in what follows, and a function v 1 (y), which is another solution to (3.1a), independent ofû 0 (y), and that is chosen to ensure that the jump conditions are satisfied. These jump conditions come from (3.6), where we have used that μ 1 = 0. Note that because of (3.12), and the periodicity properties of u 1 andû 0 , the function v 1 (y) − yû 0 (y) has to be periodic (FFLW) or antiperiodic (FFFW). Because this function will appear many times in what follows, it is worth giving it a name. Hence, we define Inputting the form (3.12) into (3.13), leads to two conditions on f 1 : One should notice, in particular, that no terms involving U 1 (x) appear in these conditions.
where Notation 2.7 has been used. Note that L 1,2 are the same operators as those applied toû 0 when determining μ 0 , though in this case the RHS was 0. Here we have these non-zero K 1,2 terms.
As in §3b(i), we can now make use of lemma 2.5 to simplify the left-hand side (LHS) of (3.22): where the jump conditions (3.1b) and (3.18) have been used. Finally, using (3.23) and dividing through by U 0 (x), (3.22) can be rewritten as which is the effective equation for U 0 .

(d) Approximation of the dispersion branches
Note that in (3.24) T = 0, but it can be either negative or positive. Since we are looking for standing waves, we seek μ 2 2 such that μ 2 2 /T ≥ 0. Remember that μ 2 is a correction term to the reduced frequency μ, such that μ 2 = μ 2 0 + δ 2 μ 2 2 + o(δ 2 ). This means that for each branch of the dispersion diagram (determined by our initial choice of eigenvalue μ 0 ), we look for a function μ 2 (κ) that will lead to an approximation of μ(κ) at the second order in δ, where κ is the reduced Bloch wavenumber. In particular, by definition of the FFLW and FFFW cases, we should have In order for our asymptotic representation (2.9) to be compatible with the fact that u δ should satisfy the Bloch-Floquet conditions ( Hence, due to the fact thatû 0 is periodic (FFLW) or antiperiodic (FFFW), we can cancel out the termsû 0 in (3.26) to get U 0 (x + δ) = ±U 0 (x) e iκδ , (3.27) where, in (3.27) and below, Notation 2.7 is being used. The second Bloch-Floquet condition in (2.8), combined with (3.27), implies that 2 /Tx for some constants A and B. The Bloch-Floquet conditions (3.27) and (3.28) lead to Since κ is restricted to (0, π/δ), i.e. to the first Brillouin zone in the dispersion diagram, and since it is assumed that μ 2 2 /T ≥ 0, (3.29) implies that, using (3.3) , we have which gives the following approximation for the reduced frequency μ(κ) The non-dimensional wave field is approximated by Note that in §6, we will find it more convenient to test the validity of (3.32) when it is written in terms of the variable y as follows u δ (δy) = U 0 (δy)û 0 (y) u 0 (δy,y) +O(κδ).
(b) In this case, we cannot conclude that μ 1 = 0 We will now apply the same methodology as in §3b(i) and consider the quantity u 1 × (3.1a (1) The exact same reasoning leads to the counterpart to (3.7): The only difference being that this time, the first bracket in the RHS of (4.2) is not zero. Instead, it can be shown directly using (4.1) that where w 0 is the Wronskian defined by The same methodology to simplify the LHS bracket in (4.2) as that used in §3b(i) can be used: first use lemma 2.6 to reduce the average bracket to a jump bracket, then use lemma 2.5 to decompose the jump bracket into four simpler jump/mean brackets that can be computed using the jump conditions (3.1b (1) ) and (3.6) to obtain which, upon regrouping the terms, dividing through by U (1) 0 (x) and using (4.1), can be rewritten where we have defined In the exact same way, we consider the quantity u 1 × (3.1a (2) ) − u (2) 0 × (3.5) to obtain 0 (x) . (4.8) Note that w 0 is the Wronskian determinant associated with (3.1a), and βw 0 can be shown to be continuous on the unit cell, hence we conclude that βw 0 is actually constant, and hence we can see that βw 0 = β(0 + )w 0 (0 + ). Sinceû (1) 0 andû (2) 0 are linearly independent, w 0 (being the associated Wronskian) is also non-zero and in this case, we cannot conclude that μ 1 = 0.
0 ) T , the two equations (4.6) and (4.8) can be recast as the first-order ODE system . (4.9) The two eigenvalues λ 1,2 of N are given by where, by the Cauchy-Schwarz inequality associated with the inner product (3.2), the quantity inside the square root is positive. The associated eigenvectors are given by and hence, upon introducing T d to be the solution to (4.9) can be written as for some constants c 1,2 . At this stage, we need to remember the Bloch-Floquet conditions (2.8), which, when applied to u 0 , imply that where Notation 2.7 has been used. Sinceû (1) 0 andû (2) 0 are linearly independent, this implies that ±U (x + δ) = e iκδ U (x). Applying this condition to (4.13), leads to the value of μ 2 1 as follows (4.14) Hence, since μ 2 = μ 2 0 + δμ 2 1 + o(δ), we obtain the linear approximations Here, the symbol ± should not be understood as per Notation 2.7, but as two different slopes, so that, near each double eigenvalue μ 0 of the dispersion diagram, we have two linear approximations with opposite slopes emerging from μ 0 . Such behaviour of the dispersion diagram, is that of so-called Dirac points.

The case of two nearby eigenvalues
In what has been done above, there is no uniform transition from the simple eigenvalue case to the double eigenvalue case. As will be seen in the examples of §6, when two simple eigenvalues are close to each other, the agreement between the dispersion diagram and the asymptotic of §3 is somewhat short-lived. In order to remedy this issue, following some ideas developed in [7,28], we will derive an asymptotic expansion for two nearby eigenvalues. Let us assume that we have two nearby eigenvalues μ (1) 0 and μ (2) 0 with their associated eigenfunctionsû (1) 0 (y) andû (2) 0 (y) that solve (3.1a)-(3.1b) and such that μ (1) 0 < μ (2) 0 . Since we are seeking an approximation that remains valid when two eigenvalues merge into one, we assume that μ (1,2) 0 (x) and u (1,2) 0 (x, y) accordingly. With such a choice of u 0 , one can show directly that we have 0 . (5. 2) The terms involving δ in the RHSs of (5.2) should be considered when collecting the terms of order δ in (2.10)-(2.12) to obtain the equation governing u 1 : and the associated jump conditions After following the exact same strategy as in §4, consider first the term u 1 × (3.1a (1) ) − u (1) 0 × (5.3) , and then the term u 1 × (3.1a (2) ) − u (2) 0 × (5.3) , to obtain two equations that can be recast in the first-order ODE system 0 ) T , the function w 0 is defined as in (4.4) and N γ is given by the parameters B 1,2 , C 1,2 and D being defined as in (4.7). Note that when taking γ → 0 in (5.5), we recover exactly (4.9), showing the consistency of our approach. Note that when deriving the second line of (5.5), we neglect the term δγ u 1 ,û (2) 0 that occurs in the process since it is of order δ.
The eigenvalues λ (γ ) 1,2 of N γ can be found to be given by while the associated eigenvectors are We can hence follow what we have done in §4c, and use the Bloch-Floquet conditions to obtain (FFLW): which are implicit relationships between μ 2 1 and κ. Fortunately, these can be inverted exactly to obtain it implies that μ 2 1 is actually complex. However, this issue is settled by realizing that upon using the inner product defined in (3.2), we can write 0 .
Hence, sinceû (1,2) 0 correspond to two different eigenvalues, they are orthogonal and we get mB 1 B 2 + D = 0. Therefore, the complex part of (5.8) disappears and it simplifies to In the limit γ → ∞ orκ → 0, μ 2 1 behaves like μ 2 1 → 0 or μ 2 1 ∼ γ , depending on the sign chosen in (5.8). This allows us to conclude that the − sign corresponds to the branch emanating from μ (1) 0 , while the + sign corresponds to the branch emanating from μ (2) 0 . There are two other interesting limits to consider.
The first is to see what happens when two eigenvalues are merging, i.e. we fixκδ and let δγ → 0. In this case (5.9) simplifies to μ 2 ≈ (μ (1) 0 ) 2 + (δγ /2) ∓ T dκ δ, which, when plotted against κδ are two straight lines with opposite slopes emanating from a point between the two nearby eigenvalues. When the two eigenvalues merge (i.e. δγ = 0), we recover exactly the double eigenvalue approximation (4.14).
The second is to see how well (5.9) approximates the dispersion diagram at the edges of the Brillouin zone for a given δγ = 0. So we fix δγ and letκδ → 0. In this case, (5.9) simplifies to Hence, for the lower branch it becomes μ 2 ≈ (μ 2 , while for the upper branch, it reads μ 2 ≈ (μ (2) 0 ) 2 + (T 2 d /δγ )(κδ) 2 . This approximation resembles (3.31), but with an incorrect quadratic coefficient (since in general ∓(T 2 d /δγ ) = T), so it is only a first-order approximation, slightly less precise than the simple eigenvalue approximation.
Hence, (5.9) is a uniform approximation, in the sense that it is valid for both simple and double eigenvalues. Moreover, we will see in the next section that using (5.9) leads to a much longer-lived fit to the exact dispersion diagram than the simple eigenvalue method.

Examples and numerical experiments
The theory developed above has the advantage to be valid for any spatially varying periodic material properties, even in cases when the dispersion diagram cannot be obtained analytically or is computationally intricate to obtain. However, in order to validate the method, we now consider two simpler examples, for which the dispersion diagram can be obtained directly by the Bloch-Floquet analysis.

(a) Monolayer
The simplest example that can be considered is the case of a monolayer material with imperfect interface. By this we mean that the density and Young's modulus are constant, so that ρ h (X) = ρ and E h (X) = E . This implies that α = β = 1. The geometry of the physical problem is represented in figure 2.
The Bloch-Floquet analysis (see appendix B) gives the following dispersion relation   The dispersion diagram classically displays band gaps as can be seen in figure 3. We will now apply the high-frequency homogenization technique to derive an analytical approximation to the higher branches of the diagram and to the associated wave fields.
In the case of a single eigenvalue, using (3.1), we find thatû 0 can be written asû 0 = A cos(μ 0 y) + B sin(μ 0 y) for some constants A and B, and, using lemma 2.4, it is subject to the jump conditions ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩û where the 2 × 2 matrix M mo = (M mo ij ) is given by The only way for non-trivial solutions to (6.3) to exist is for its determinant to be zero, which after some algebraic manipulations leads to a dispersion relation of the form D mo (μ 0 ; m, k) = 2(1 ∓ cos(μ 0 )) + m 2k μ 2 0 (1 ± cos(μ 0 )) ± μ 0 sin(μ 0 ) m + 1 k = 0, (6.4) where the fact that μ 0 = 0 has been used (we are not interested here in the low frequency limit). In practice, when calculating μ 0 and reconstructingû 0 , it can be useful to note that so that μ 0 is either a zero of M mo 11 or M mo 21 , and in the former (resp. latter) case, the top (resp. bottom) line of M mo is zero (it can be shown that sin(μ 0 ) = 0). The computed eigenvalues coincide with the edges of the Brillouin zone of the dispersion diagram of figure 3. To obtainû 0 , we set A = 1, so thatû 0 (0 + ) = 1, and compute B using the first (resp. second) line in (6.3) if M mo 11 = 0 (resp. M mo 21 = 0). Because we are ultimately interested in the value of T in (3.31), we need to calculate βw 1 in (3.24), and hence v 1 on the interval (0, 1). We do this using the fact that it satisfies the same second-order equation (3.1a) asû 0 and can hence be written v 1 (y) = C cos(μ 0 y) + D sin(μ 0 y), for some constants C and D, which, using remark 3.1, can be found by solving M mo (C, D) T Since M mo is singular, (6.6) does not have a unique solution so we set C = 1 say and use the nontrivial line of the system to determine D. This works well since it can be shown that b mo is such that b mo j = 0 whenever M j,1 = 0. Once v 1 is found, the resulting value of T is directly obtained using (3.24). Note that, in this simple case, no numerical integration is required and calculations can be performed analytically. The resulting second-order approximations μ ≈ μ 0 + (T/2μ 0 )(κδ) 2 are superposed to the dispersion diagram in figure 3a, and one can see that they approximate the branches well in the vicinity of the edges of the Brillouin zone. One should also note that, as seen in figure 4 (left), this approximation remains valid within the band gaps, where κδ is complex and is such that Re(κδ) = 0. This is expected since in these cases, (κδ) 2 remains real. Having computed all the eigenvalues and eigenfunctions, it is now straightforward to compute T d as per (4.12) and δγ as per (5.1), where the pairs of eigenvalues are chosen naturally according to the dispersion diagram (first-second) and (third-fourth) in both the FFLW and FFFW cases. Hence we can evaluate the nearby eigenvalue approximation (5.9) derived in §5. It is displayed in figure 3b, and as can clearly be seen, even if this approximation is only first-order in the vicinity of the edges of the Brillouin zone, its agreement with the dispersion diagram is much longer lived than that of the simple eigenvalue approximation. Similar observations are true within the band gaps as can be seen in figure 4b.
We will now investigate the accuracy of the zeroth-order field approximation obtained in the simple eigenvalue case. Using the Bloch-Floquet analysis (see appendix B), we can have access to the exact standing wave field u δ (x) = u δ (δy), and, to be compatible with the asymptotic expansion (2.9), we normalize it such that u δ (0 + ) = u 0 (0, 0 + ). Note that because of (2.7) and (3.33), the difference between the exact and approximated field can be written as  To illustrate the validity of our approximation, we hence compare u δ (δy) andû 0 (y) for various values of κδ in figure 5, showing, as expected, that asκδ gets smaller the zeroth-order field is a good approximation to the exact field. A similar investigation can be carried out for the nearby eigenvalue approximation of the field. As can be seen in figure 6, the approximation is good even for a value ofκδ = 0.5 that is not particularly small, and for which the agreement of the simple eigenvalue zeroth-order field is poor.
In figure 7, we plot the error between the exact field and the homogenized fields obtained using both the simple and nearby eigenvalue approximations. Due to the periodicity properties of u δ andû 0 and (6.7), the error is relatively easy to compute in the simple eigenvalue case and is defined as follows:  In the nearby eigenvalue approximation, using (4.1), (5.7) and (4.13), the errors can be written as (1) 0 (y) +û  One can see that we recover the expected behaviour u δ (δy) = u 0 (δy, y) + O(κδ), but that the nearby eigenvalue approximation performs much better for the whole range of values ofκδ used in figure 7.
We now endeavour to study how the eigenvalues μ 0 depend on m and k. In order to visualize this we display a heat map of the first and second FFLW (periodic) μ 0 (m, k) in figure 8. One can clearly see two distinct regions, on the left and on the right of the curve k = 1/m. On each side of these curves, the eigenvalues depend solely on one of the two parameters (m, k), which  In the case of a double eigenvalue, we have to follow the procedure of §4 by representing u 0 as in (4.1). We hence need to find two independent solutionsû (1,2) 0 of (3.1), which can both be writtenû (1,2) 0 (y) = A (1,2) cos(μ 0 y) + B (1,2) sin(μ 0 y). Because of the dimension of the system, any two independent vectors (A (1) , B (1) ) T and (A (2) , B (2) ) T would work, and we can hence choose (A (1) , B (1) ) T = (1, 0) T and (A (2) , B (2) ) T = (0, 1) T . As seen in §4, the two functions U (1,2) 0 (x) appearing in (4.1) satisfy the first-order ODE system (4.9), where in our case we have  (1) b (1) a (2) b (2) (a) ( b) Figure 10. Geometry of the bilayer problem (a) and non-dimensional (b) settings. (Online version in colour.) Using these, one can easily compute the associated eigenvalues λ 1,2 and eigenvectors U λ 1,2 via (4.10) and (4.11). Note that in this case, one can show that so that the eigenvalues λ 1,2 and eigenvectors U λ 1,2 of the matrix of the ODE system are simply and T d = μ 0 /(mB 2 1 + C 1 ). We can hence superpose the resulting linear approximation (4.15) onto the dispersion diagram, revealing an excellent fit, as can be seen in figure 9b. It is quite remarkable that in this case, every eigenvalue μ 0 corresponds to Dirac points. In fact this can be understood by considering a homogeneous material with only one spring-mass interface. Upon sending a wave onto this interface, one can naturally derive a coefficient of reflection Ref(μ) and a coefficient of transmission Trans(μ). It turns out that Trans(μ) = 2 1 + (mμ 2 /4k) and therefore the reflection coefficient is zero if and only if the condition k = 1/m is satisfied. Hence, in the periodic medium considered, no internal reflection can be present, no destructive/constructive interference can take place and no band gaps occur.

(b) Bilayer
We now consider the case of a bilayer material characterized by the phase fraction r ∈ (0, 1), and hence provide the imperfect interface extension to the example given in [6]. The unit cell is made up of two homogeneous materials. The first one has a length rh, density ρ 1 and Young's modulus E 1 , while the second has length (1 − r)h, density ρ 2 and Young's Modulus E 2 . The two respective wave speeds are c 1 = √ E 1 /ρ 1 and c 2 = √ E 2 /ρ 2 . The important non-dimensional functions α and β are hence defined by for y ∈ (r, 1), for y ∈ (r, 1), The interface at y = r is assumed perfect, and the geometry of this physical problem is summarized in figure 10. The classic Bloch-Floquet analysis will, in this case, give the following dispersion relation where Z = ρ c , and C i = cos(μH i ), S i = sin(μH i ), H 1 = rc /c 1 , H 2 = (1 − r)c /c 2 and, naturally, c = √ E /ρ . As per the monolayer case, the dispersion diagram displays band gaps as can be seen in figure 11. We will now apply the high-frequency homogenization technique to derive an analytical approximation to the higher branches of the diagram and to the associated wave fields.
We now need to find v 1 to obtain a second-order approximation. Since it is a solution to the same equation (3.1)a, we can write v 1 (y) = C (1) cos(Ω (1) y) + D (1) sin(Ω (1) y) on (0,r), C (2) cos(Ω (2) y) + D (2) sin(Ω (2) y) on (r, 1). (6.17) Using remark 3.1 for the conditions at the unit cell interfaces, and remembering that, according to remark 3.2, both v 1 and βv 1 should be continuous at y = r, one obtains a system of the form M bi (C (1) , D (1) , C (2) , D (2) Because the matrix M bi is singular, we compute (C (1) , D (1) , C (2) , D (2) ) T via the Moore-Penrose Pseudo-inverse [29]. Once v 1 is found, the resulting value of T is directly obtained using (3.24). The resulting approximations μ ≈ μ 0 + (T/2μ 0 )(κδ) 2 are superposed to the Bloch-Floquet diagram in figure 11a, and one can see that they approximate the branches well in the vicinity of the edges of the band gaps of the dispersion diagram. It is apparent from figure 11 that the highest antiperiodic eigenvalue displayed seems to be a double eigenvalue (in fact we will see further that it is not exactly a double eigenvalue), and that the approximation is particularly short-lived in this neighbourhood. Having computed all the eigenvalues and eigenfunctions, we can once again evaluate the nearby eigenvalue approximation (5.9). It is displayed in figure 11b, and as can clearly be seen, its agreement with the dispersion diagram is much longer-lived that of the simple eigenvalue approximation, even more so for the near-double eigenvalue.
Using the Bloch-Floquet analysis, in a very similar way to the monolayer case, we can have access to the exact standing wave field u δ (x) = u δ (δy), and we normalize it such that u δ (0 + ) = u 0 (0, 0 + ). As per the monolayer case, and because of (6.7), we illustrate the convergence of the simple eigenvalue method by comparing u δ (δy) andû 0 (y) for various values of κδ in figure 12.
A similar investigation can be carried out for the nearby eigenvalue approximation. As can be seen in figure 13, the approximation is good even for a value ofκδ = 0.5 that is not particularly small, and for which the agreement of the simple eigenvalue zeroth-order field is poor.  In figure 14, we plot the error between the exact field and the homogenized fields obtained using both the simple and nearby eigenvalue approximations. The errors can be expressed as in (6.8)-(6.11). Again, one can see that we recover the expected behaviour u δ (δy) = u 0 (δy, y) + O(κδ) , but that the nearby eigenvalue approximation performs better for the whole range of values ofκδ used in figure 14.
As for the monolayer example displayed in figure 9, it is possible to find physical parameters such that all the eigenvalues become simultaneously double eigenvalues (Dirac points). In order to do so one needs to ensure that k = 1/m, and that ρ 1 c 1 = ρ 2 c 2 . The first condition imposes that the reflection coefficient due to the imperfect interface is zero, and the second imposes that the two homogeneous materials are "impedance matched" so that no reflection occurs from their perfect interface either.  However, for the bilayer, it appears that certain parameters lead to only two of the eigenvalues merging into a double eigenvalue, as appears to be the case in figure 11. In order to visualize this phenomena, one could look at the evolution of the eigenvalues μ 0 for fixed physical parameters, but for varying r within (0, 1). The results are displayed in figure 15, and it seems that double eigenvalues or near-double eigenvalues may occur for some specific values of r, though, in this case, all the eigenvalues do not become double simultaneously. In fact, as illustrated in figure 15, if one zooms in on the areas of the graphs where eigenvalues seem to coincide, it appears that the curves do not actually touch each other. We will call these points almost-Dirac points. As highlighted above, the nearby eigenvalue approximation to the dispersion diagram is excellent for such almost-Dirac points, while the simple eigenvalue method leads to a very short-lived approximation, see figure 11.
To explore the parameter space further, we will keep the values of r, α (1,2) and β (1,2) used in figure 11 and study the variation of the fifth and sixth antiperiodic eigenvalues μ 0 that correspond to an almost-Dirac point according to figure 15. As can be seen in figure 16, in this case, we observe a similar behaviour as that of figure 8, where two distinct regions seem to be separated by a smooth curve, on which the values of m and k chosen in figure 11 (represented by a black star) seem to lie.   (1,2) and β (1,2) used in figure 11. The black star corresponds to the values of m and k used in figure 11. (Online version in colour.) All in all, it seems to be the case, that for a given integer j, the jth eigenvalue is double on some manifold given by F j (ρ 1 , ρ 2 , E 1 , E 2 , M, K, r) = 0 for some function F j , though finding an analytical expression for F j is beyond the scope of the present work.

Conclusion
In this work, we have extended the technique of high-frequency homogenization to onedimensional periodic media with linear imperfect interfaces of the spring-mass type. The extension was not direct, and many of the proofs for the classic case needed to be extended in order to deal with the extra technical difficulties arising from an imperfect interface. We have also described how the technique should be modified in the specific case of Dirac points within the dispersion diagram, and also proposed a uniform approximation taking into consideration competing nearby eigenvalues.
We have illustrated the validity of our theoretical development with the two examples of monolayered and bilayered materials. In the case of the monolayered material, we quantified the error between the exact and the homogenized fields, and we found a simple condition on the non-dimensional stiffness and mass values k and m for all the points at the edges of the Brillouin zone to become Dirac points. Similar conditions were given in the bilayered case. Moreover, in the bilayered case, almost-Dirac points have been identified while exploring the parameter space (much larger than that of the monolayered material). For both examples considered, we found that the nearby eigenvalue approximation led to a much longer lived approximation to both the dispersion diagram and the wave fields. This formulation is also convenient because it does not break down when two eigenvalues merge, unlike the simple eigenvalue approximation. As mentioned previously, one of the advantages of high-frequency homogenization is that it works even when the dispersion diagrams cannot be obtained analytically or are computationally intricate to obtain. In such cases, this high-frequency homogenization approach may, for example, provide a way of back-engineering the values of m and k describing imperfect interface, or pick specific material properties and contacts in order to ensure the presence of Dirac points, which are known to display very interesting physical properties.
In the future, we hope to be able to push the asymptotics presented in this paper to higher order, i.e. to propose first-and maybe second-order corrections to the leading-order wave fields exhibited in the present work.
Data accessibility. This article has no additional data.