# A stationarity principle generating effective boundary conditions for second-order homogenization

Manon Thbaut, Basile Audoly, Claire Lestringant

<figure><img src="image_1.png" /></figure>

## To cite this version:

Manon Thbaut, Basile Audoly, Claire Lestringant. A stationarity principle generating effective boundary conditions for second-order homogenization. 2025. hal-05129635

HAL Id: hal-05129635

https://hal.science/hal-05129635v1

Preprint submitted on 3 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A stationarity principle generating effective boundary conditions for second-order homogenization

MANON THBAUT

Laboratoire de Mécanique des Solides
CNRS, Institut Polytechnique de Paris
91120 Palaiseau, France

BASILE AUDOLY

Laboratoire de Mécanique des Solides
CNRS, Institut Polytechnique de Paris
91120 Palaiseau, France

CLAIRE LESTRINGANT

Institut Jean Le Rond d'Alembert,
Sorbonne Université, CNRS,
75005 Paris, France

We derive an effective model for a periodic chain of linearly-elastic springs, achieving second-order accuracy in the scale separation parameter $\varepsilon \ll 1$. The chain has finite length and is made up of springs connecting both nearest- and next-nearest-neighbors: it serves as a one-dimensional prototype for higher-order periodic homogenization problems with boundaries. This type of problem has been approached by inserting two-scale expansions into the equations of equilibrium in the bulk and by matching them with boundary-layer solutions. We explore an alternative method operating at the energy level, bypassing the cumbersome matching procedure. We start from an ansatz of the microscopic displacement accounting for both boundary layers and for small-scale fluctuations in the bulk, and insert it into the discrete energy. This yields a continuous energy functional depending on the macroscopic displacement $u$, in the form of a series expansion in powers of $\varepsilon$. We call it a pseudo-energy $\Phi_{\varepsilon}[u]$ as it is not positive when truncated at order $\varepsilon^2$. The boundary terms in the pseudo-energy account for boundary layers in an effective way. By making the pseudo-energy stationary order by order in $\varepsilon$, we derive the homogenized equations of equilibrium along with effective boundary conditions. We provide quantitative validation showing that the effective model is correct to second order. We point out the special form of the effective higher-order tractions, which has been overlooked in strain-gradient theories proposed so far.

**Keywords** Variational asymptotic homogenization, Higher-order homogenization, Linear elasticity

# 1. INTRODUCTION

## 1.1. Motivation and context

Asymptotic homogenization [BLP78; San80; BP89] is a powerful tool to derive effective continuum models that accurately capture the mechanical behavior of complex microstructures in the limit of a small scale separation parameter, $\varepsilon \ll 1$. It is known that higher-order terms, such as the terms of order $\varepsilon^2$ involving the second gradient of the displacement, substantially improve the accuracy of the homogenized models [S18; Bou96; BS11; GK89; SC00; TB96; BP89]. A blind spot of existing approaches is the definition of effective boundary conditions preserving higher-order accuracy. The effective boundary conditions cannot be obtained using homogenization only, due to the presence of 'microscopic' boundary layers that require a dedicated analysis [BLP79; San86; Dum86; AA99; MV97; YAL24; Vin16]. By combining boundary layer analysis with asymptotic homogenization in the bulk, effective boundary conditions have been derived to order $\varepsilon^1$ [AA99; Vin16] and more recently to order $\varepsilon^2$ [TAL24].

It is unclear whether the effective boundary-value problem thus obtained has a variational structure, i.e., whether it can be obtained as the minimum condition of an effective energy. Identifying a variational structure would be beneficial in multiple ways. (i) When the number of kinematic variables increases (for example in 2D or 3D structures or when additional kinematic descriptors such as micro-rotations are involved [NCH20]), the matching step must be performed on each variable separately, which can be tedious. This step is replaced by a straightforward optimization in the variational approach. (ii) An energy approach provides a flexible and compact formulation, in which the boundary conditions have a simple interpretation: dual mechanical quantities can be identified directly in the energy. It thus lends itself easily to extensions such as nonlinear elastic constitutive laws, pre-stress or geometric non-linearities [KAL22; LA20]. (iii) When a variational principle is available, the effective equilibrium equations can be implemented using standard and robust numerical methods, such as the finite-elements method.

In prior work approaching asymptotic homogenization at the energy level [S18; AL23; BKS11; LM18; Lio81; AE09], boundary effects have been ignored. A notable exception is the work of [Pan09], whose general strategy we adopt here: he derives an effective energy by inserting relocalization formulas for the microscopic displacement into the microscopic energy. The relocalization formulas introduce the macroscopic displacement $u_{\varepsilon}$ as the primary unknown of the formulation. They include boundary-layer correctors, which makes them uniformly accurate. When made stationary, the effective energy generates an equilibrium problem whose solution is a high-order approximation of the macroscopic displacement.
The approach of [Pan09] provides a mathematically sound answer to the difficult problem of high-order homogenization in the presence of boundaries. Unfortunately, the boundary-value problems it delivers are cumbersome, even for the simple one-dimensional illustration chosen by the author [Pan09]. In the present work, we propose an alternative approach delivering boundary-problems that are significantly more compact, easier to interpret, and preserve the second-order accuracy. Our starting point is the observation that, in the traditional approach to high-order homogenization, the macroscopic displacement $u_{\epsilon}$ is derived *order by order* in $\epsilon$. The existence of this perturbative solution puts a severe constraint on the expression of the boundary load appearing in the effective model of Panasenko [Pan09]. The special form of boundary load has not been appreciated in earlier work on strain-gradient theory [Tou62; Min64; Min65] nor in the seminal work of [Pan09], and we leverage it to simplify the effective model.

We demonstrate our approach on the example of a 1D spring network featuring both nearest- and next-nearest-neighbors interactions. This is one of the simplest lattice examples exhibiting boundary layers [GW65; CT02; CZWK19; TAL24]. The task is broken up into the following sequence of steps.

* The starting point of the asymptotic expansion is an ansatz for the microscopic displacement that extends those used in the context of homogenization, with an additional contribution accounting for boundary layers: we write the nodal displacement as the sum of (i) a contribution tied to the slowly variable macroscopic displacement $u$, relevant to the interior of the domain (bulk), and (ii) a localized contribution decreasing exponentially quickly away from the endpoints, representing boundary layers.

* Following the standard procedure in homogenization, the first contribution (i) relevant to the bulk is itself written as the sum of the effective displacement $u$ and a small-amplitude, quickly oscillatory microscopic shift, which we solve in terms of the successive gradients of the effective displacement.

* The form of the boundary layer term (ii) is identified by solving the original discrete equilibrium equation at interior nodes with the applied load set to zero. The amplitude of the boundary layer remains undetermined at this stage.

* We then insert this ansatz in the expression of the discrete energy, and transform it into an integral plus boundary terms. To do so, we first perform a Taylor-expansion of the nodal quantities at the interior nodes [AG05; TB93] and then use the Euler-Maclaurin formula to replace the discrete sums over the nodes and springs by a (continuous) integral. The Euler-Maclaurin formula generates boundary terms of order $\epsilon$ and $\epsilon^2$ that add up to those produced by the boundary-layer solution.

* We identify the boundary-layer amplitude by an energy minimum condition.

The result of this procedure is a formal expansion $\Phi_{\epsilon}[u]$ of the energy approximating the microscopic energy to order $\epsilon^2$, see Equation (1.1) below. We will refer to this expansion as a *pseudo-energy*. It is a functional of the macroscopic displacement $u$, and is made up of an integral term (bulk contribution) and boundary terms. The strain energy is present in both the integral and in the boundary terms, and is expressed in terms of the macroscopic strain $u'(x)$ and its successive derivatives: the pseudo-energy is of the 'strain-gradient' type.

By making this pseudo-energy stationary, we derive a boundary-value problem that can be solved order by order. As we show, the solution $u$ is a macroscopic displacement that is asymptotically exact up to order $\epsilon^2$. The *perturbative* boundary-value problem thus obtained agrees, to order $\epsilon^2$, included, with the one we derived by matched asymptotic expansions in our previous work [TAL24]. The effective boundary conditions include higher-order corrections that capture the boundary effects in an effective way.

The *perturbative stationarity principle* for the pseudo-energy derived in this paper is weaker than the energy-minimum principle of [Pan09]. Our pseudo-energy does not meet the standard positivity requirements associated with the notion of energy. In addition, the order-$\epsilon^2$ solution $u=u_{(0)}+\epsilon u_{(1)}+\epsilon^2 u_{(2)}$ which we obtain perturbatively is not the minimum, nor even the stationary point, for fixed $\epsilon>0$, of a well-identified functional parameterized by $\epsilon$ [Bou20]. We will lift these important limitations in a follow-up paper [TAL25a] where we extend the perturbative stationarity principle to a full-fledged energy-minimum principle. The present work is thus only an intermediate step towards the identification of a proper strain-gradient formulation of the effective model. The strain-gradient theory obtained in this way will appear to be significantly simpler than that delivered by the truncation-limited-to-strain approach of [Pan09].

This paper is focused on *one-dimensional* lattices. To ease future extension to periodic continua in dimension 2 or 3, we refrain as much as possible from using any property that is specific to dimension 1, such as integrating by parts strain-gradient energy contributions of the form $\int_0^\ell u' u''' dx = -\int_0^\ell u''^2 dx + [u'u'']_0^\ell$. We limit attention to *discrete* lattices but the extension to *elastic continua* possessing periodic elastic properties, such as an elastic bar having a periodic traction modulus, raises no difficulty. We also emphasize that our goal is not to generate original results, but instead to revisit the boundary-value problems generated by matched-asymptotic expansions in earlier work, such as [TAL24], from an energy perspective. We use the method of formal asymptotic expansions, assuming that the loading is smooth enough for the macroscopic displacement $u(x)$ to have a regular dependence on $\epsilon$. Rigorous proofs of the kind that can be offered by $\Gamma$-developments and expansions [BT08] are left for future work.
## 1.2. Outline

In Section 2, we introduce some mathematical notations and tools that are used throughout the paper. In Section 3, we introduce the one-dimensional discrete lattice model. It is made up of a unit cell repeated periodically in the interior, the pattern being modified locally to accommodate the endpoints. The nodal displacements, spring elon- gations (microscopic strain measure) and discrete total potential energy are defined, and the discrete equations of equilibrium are presented.

The two-scale analysis is implemented in the subsequent sections as follows.

In Section 4, we introduce the ansatz for the nodal displacement, which is the key ingredient of the two-scale analysis. It makes use of a macroscopic displacement $u(x)$, a microscopic shift $y(x)$ accounting for the small, short-scale oscillations in the bulk, and a boundary-layer solution $Z_i$, see (4.2).

In Section 5, we use a classical homogenization procedure at second order to solve the microscopic shift $y(x)$ in terms of the successive derivatives of the macroscopic strain $u'(x)$. From this, we obtain the expression of the microscopic displacement in terms of $u(x)$ and its derivatives. This expression is valid far away from the endpoints, as the boundary layers have been neglected temporarily.

In Section 6, the boundary layer solution is derived at each of the endpoints, up to a small set of degrees of freedom ($v^0, v^\ell$), see (6.12). At this stage, the nodal displacement given by the ansatz (4.2) has been determined in terms of the macroscopic displacement $u(x)$ and of the boundary-layer amplitudes $v^0$ and $v^\ell$ only.

In Section 7, we insert this expression of the nodal displacement into the discrete energy, and use an order-$\epsilon^2$- accurate continualization procedure to rewrite the discrete sums as an integral. This yields a first expression of the pseudo-energy in terms of $u$ and ($v^0, v^\ell$), see (7.3). Next, we minimize this energy with respect to the boundary-layer amplitudes ($v^0, v^\ell$). The matching procedure used in the classical approach to two-scale expansions is bypassed, and delivers the pseudo-energy in the form of a functional $\Phi_\epsilon[u]$ that is a series expansion in powers of $\epsilon$, see (7.7),

$$ \Phi_{\epsilon}[u] = \int_{0}^{\ell} \left( \frac{1}{2} (u' \epsilon u'' - \epsilon^2 u''') \cdot \bar{\mathbf{K}} \cdot (u' \epsilon u'' - \epsilon^2 u''') - fu \right) dx \\ + \sum_{b \in \{0, \ell\}} \left( \begin{array}{l} \frac{\epsilon}{2} (u'(b) \epsilon u''(b)) \cdot \bar{\mathbf{K}}^b \cdot (u'(b) \epsilon u''(b)) \\ - \bar{\mathbf{K}}^b \cdot (p^b \epsilon f(b) \epsilon^2 f'(b)) u(b) \\ - \epsilon (\bar{\mathbf{R}}^b \cdot (p^b \epsilon f(b))) \cdot (u'(b) \epsilon u''(b)) \\ - \frac{\epsilon}{2} (p^b \epsilon f(b)) \cdot \bar{\mathbf{C}}^b \cdot (p^b \epsilon f(b)) \end{array} \right) + \mathcal{O}(\epsilon^3). \quad (1.1) $$

The effective coefficients are denoted with a bar on top and are obtained in terms of the lattice properties. The various terms are discussed in detail in the main text. Briefly, and in order of appearance in equation above, we have a bulk strain-gradient energy ($\bar{\mathbf{K}}$), a potential of distributed forces ($f$) in the integral, a boundary elasticity ($\bar{\mathbf{K}}^b$), an applied traction ($\bar{\mathbf{I}}^b$), higher-order tractions ($\bar{\mathbf{R}}^b$) and a pure-loading energy ($\bar{\mathbf{C}}^b$).

In Section 8, we derive the boundary-value problem governing equilibrium, by making the pseudo-energy $\Phi_\epsilon[u]$ stationary. The resulting equilibrium problem is solved order by order in $\epsilon$. We observe that, for this solution to exist, the effective higher-order tractions have to take a special form, see Equation (8.14).

In Section 9, we extend the above results, applicable to force-controlled boundaries, to the cased of fixed end-points.

Finally, we provide validation examples in Section 10. We compare our predictions to the solution of the microscopic problem. The comparison is done in two geometries, one that recycles extensive validation done in our previous work, and another one where we generate the microscopic solution from scratch. In both cases, the effective model is accurate to second order.

All the calculations presented in this paper were done in the symbolic language Wolfram Mathematica. The notebook which was used to generate the results is available as a Supplementary material [TA25].

## 2. MATHEMATICAL TOOLBOX

In this preliminary section, we introduce a convenient mathematical notation that makes it possible to express pseudo-energies in compact form (§2.1). The notation can also handle Taylor expansions efficiently (§2.2). We also introduce the Euler–Maclaurin formula, which will be instrumental in continualizing the pseudo-energy (§2.3), and specify the index-pairing convention for the double contraction of two tensors (§2.4).

### 2.1. Vector-of-derivatives notation

For an arbitrary function $g(x)$, we denote as $\dot{g}$ (respectively, $\ddot{g}$) the vector collecting $g$ and its first (respective first two) derivative(s), scaled by the powers of the scale separation parameter $\epsilon$, to be defined later,

$$ \dot{\mathbf{g}}(x) = \begin{pmatrix} g(x) \\ \epsilon g'(x) \end{pmatrix}, \quad \ddot{\mathbf{g}}(x) = \begin{pmatrix} g(x) \\ \epsilon g'(x) \\ \epsilon^2 g''(x) \end{pmatrix}. \quad (2.1) $$
The following identities will be used.

First, $\ddot{g}(x)$ is obtained from $\dot{g}(x)$ by appending a quantity of order $\epsilon^2$, so we have

$$ 
\begin{aligned}
\ddot{\mathbf{g}}(x) &= \begin{pmatrix} g(x) \\ \epsilon g'(x) \\ \mathcal{O}(\epsilon^2) \end{pmatrix} \\
&= \mathcal{P}^{\top} \cdot \dot{\mathbf{g}}(x) + \mathcal{O}(\epsilon^2),
\end{aligned}
\qquad (2.2) $$

where $\mathcal{P}$ is the truncation matrix, and $\mathcal{P}^{\top}$ appends a trailing 0 to a vector of dimension 2,

$$ \mathcal{P} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix}. \qquad (2.3) $$

Next, we define the shift matrix $\mathcal{R}$ discarding the first entry of a vector of dimension 3,

$$ \mathcal{R} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}. \qquad (2.4) $$

It connects the vector-of-derivatives $\dot{\mathbf{g}}'=(\dot{\mathbf{g}}')$ built from the derivative $g'$ to the original vector $\ddot{\mathbf{g}}$, as in

$$ 
\begin{aligned}
\epsilon \dot{\mathbf{g}}' &= \begin{pmatrix} \epsilon g'(x) \\ \epsilon^2 g''(x) \end{pmatrix} \\
&= \mathcal{R} \cdot \ddot{\mathbf{g}}(x).
\end{aligned}
\qquad (2.5) $$

Combining (2.2) and (2.5), we have

$$ 
\begin{aligned}
\epsilon \ddot{\mathbf{g}}'(x) &= \epsilon \mathcal{P}^{\top} \cdot \dot{\mathbf{g}}'(x) + \mathcal{O}(\epsilon^3) \\
&= \mathcal{P} \cdot \ddot{\mathbf{g}}(x) + \mathcal{O}(\epsilon^3).
\end{aligned}
\qquad (2.6) $$

where $\mathcal{S}$ is the 3x3 shift matrix,

$$ \mathcal{S} = \mathcal{P}^{\top} \cdot \mathcal{R} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (2.7) $$

## 2.2. Taylor expansions in compact form

Given a scalar $a$, and denoting the matrix exponential as 'exp', we have

$$ \exp(a\mathcal{S}) = \mathbf{I}_3 + a\mathcal{S} + \frac{a^2}{2}\mathcal{S}^2 = \begin{pmatrix} 1 & a & \frac{a^2}{2} \\ 0 & 1 & a \\ 0 & 0 & 1 \end{pmatrix}, \qquad (2.8) $$

where $\mathbf{I}_n$ denotes the identity matrix in dimension $n$. The infinite series defining the matrix exponential stops at $a^2/2\mathcal{S}^2$ as $\mathcal{S}^3=0$ for the particular matrix $\mathcal{S}$.

The matrix exponential makes it possible to write the Taylor expansion of $g$ and its derivatives about a point $x$ in compact form as

$$ 
\begin{aligned}
\ddot{\mathbf{g}}(x+a\epsilon) &= \begin{pmatrix} g(x) + a\epsilon g'(x) + \frac{a^2\epsilon^2}{2}g''(x) \\ \epsilon(g'(x) + a\epsilon g''(x)) \\ \epsilon^2 g''(x) \end{pmatrix} + \mathcal{O}((a\epsilon)^3), \\
&= \exp(a\mathcal{S}) \cdot \ddot{\mathbf{g}}(x) + \mathcal{O}((a\epsilon)^3).
\end{aligned}
\qquad (2.9) $$

## 2.3. Euler-Maclaurin formula

Given a smooth function $g(x)$ that decays for $x \to \infty$ and a scalar offset parameter $a$, the discrete sum of $g$ over the infinite grid $(\epsilon a, \epsilon(a+1), \epsilon(a+2),...)$ can be approximated by an integral,

$$ \epsilon \sum_{i=0}^{\infty} g(\epsilon(i+a)) = \int_{0}^{\infty} g dx + \epsilon \left(\frac{1}{2} - a\right) g(0) - \frac{\epsilon^2}{2} \left(a^2 - a + \frac{1}{6}\right) g'(0) + \mathcal{O}(\epsilon^3). \qquad (2.10) $$

This equality is established in Appendix A, see (A.4), based on the Euler-Maclaurin summation formula [Apo99]. The standard form of the Euler-Maclaurin formula corresponds to the particular case $a=0$ in equation above, where the first grid point $\epsilon a=0$ matches the first bound of the integral.

The boundary terms at $x=0$ in the right-hand side of (2.10) ensure $\epsilon^2$ accuracy. Different values will be assigned later to the offset parameter $a$ depending on whether we continually the strain energy or the potential energy.

A particular case of interest is when $g$ is of the special form

$$ g(x) = (\mathbf{n} \cdot \dot{\mathbf{p}}(x)) \cdot \ddot{\mathbf{q}}(x), \qquad (2.11) $$
for some 3×3 matrix $n$ and for smooth functions $p(x)$ and $q(x)$. The derivative $g'$ is then given by (2.6) as

$$ 
\begin{aligned}
g'(x) &= (n \cdot \ddot{p}'(x)) \cdot \ddot{q}(x) + (n \cdot \ddot{p}(x)) \cdot \ddot{q}'(x) \\
&= \frac{1}{\varepsilon} ((n \cdot \mathcal{S} \cdot \ddot{p}(x)) \cdot \ddot{q}(x) + (n \cdot \ddot{p}(x)) \cdot \mathcal{S} \cdot \ddot{q}(x)) + \mathcal{O}(\varepsilon^2) \\
&= \frac{1}{\varepsilon} ((n \cdot \mathcal{S} + \mathcal{S}^\top \cdot n) \cdot \ddot{p}(x)) \cdot \ddot{q}(x);
\end{aligned} 
\quad (2.12) $$

and the Euler-Maclaurin formula (2.10) becomes

$$ 
\begin{aligned}
\varepsilon \sum_{i=0}^{\infty} (n \cdot \ddot{p}(\varepsilon(i+\alpha))) \cdot \ddot{q}(\varepsilon(i+\alpha)) &= \int_{0}^{\infty} (n \cdot \ddot{p}) \cdot \ddot{q} dx \\
&\quad + \varepsilon \left( \left( (1/2-\alpha)n - \frac{\alpha^2-\alpha+1/6}{2} (n \cdot \mathcal{S} + \mathcal{S}^\top \cdot n) \right) \cdot \ddot{p}(0) \right) \cdot \ddot{q}(0) + \mathcal{O}(\varepsilon^3).
\end{aligned} 
\quad (2.13) $$

## 2.4. Double contraction

Given two tensors $b$ and $c$ with respective dimensions $n_1 \times \cdots \times n_p \times k \times k'$ and $k \times k' \times m_1 \times \cdots \times m_q$, their double contraction $b:c$ is the tensor of dimension $n_1 \times \cdots \times n_p \times m_1 \times \cdots \times m_q$ obtained by contracting the last two dimensions of $b$ with the first two dimensions of $c$, $(b:c)_{i_1 \cdots i_p a a' c a a' j_1 \cdots j_q} = \sum_{\alpha=1}^{k} \sum_{\alpha'=1}^{k'} b_{i_1 \cdots i_p a a' c a a' j_1 \cdots j_q}$.

## 3. DISCRETE LATTICE

In this Section we introduce a discrete one-dimensional lattice of springs, sketched in Figure 3.1. It is made up of a finite number $N$ of identical unit cells (dashed grey boxes in Figure 3.1) that connect to the terminal nodes using a different spring pattern (black half-springs appearing outside any dashed box). The lattice is made up of two different types of nodes (open circles and solid disks, respectively) and of four different types of springs (colored rectangles). The definition of the unit cell is to some extent arbitrary, and we use the convention that a unit cell contains one instance of each of the springs labelled $II$ to $IV$ in the figure (brown, green and purple rectangles), plus one half of each of the two instances of the springs labelled $I$ and $I'$ (blue half-rectangles in the figure)—we will see shortly what we mean by this exactly. This particular definition of the unit cell is invariant by mirror-symmetry, up to a permutation of the spring constants. This will make it possible to skip the analysis of the boundary layer forming at the second endpoint entirely using a symmetry argument.

### 3.1. Lattice properties

<figure><img src="image_3.png" /><figcaption>**Figure 3.1.** A one-dimensional lattice with nearest- and next-nearest-neighbor interactions. There are four different types of springs (colored boxes, roman numbers) and two different types of nodes labelled $\pm$ (open versus solid disks). Springs are grouped into cells labelled by an integer index $j$ in the range $0 \le j \le N-1$ (dashed rounded boxes), with average position $x^j = \varepsilon(j+1/2)$. Nodes are grouped by pairs indexed by an integer $i$ in the range $0 \le i \le N$. The average position of the pair is $\varepsilon_i$. The spatial periodicity is $\varepsilon$ and there are incomplete cells on both endpoints (half-springs of type $I$ in black). The equivalent continuum covers the domain $0 \le x \le \ell = N\varepsilon$, by convention, and is represented by the thick orange line.</figcaption></figure>

Given an integer number of cells $N$ and the length $\ell$ of the effective continuum (thick orange line in the figure), one can generate a discrete lattice by setting the cell size $\varepsilon$ as $\varepsilon = \ell/N$. In our convention, the continuous limit is obtained by keeping $\ell$ fixed, and letting $N$ go to infinity and $\varepsilon$ go to zero: in the asymptotic analysis, we will take the cell size $\varepsilon$ as the small expansion parameter, $\varepsilon \to 0$, and will view $\ell = \mathcal{O}(\varepsilon^0)$ as a parameter independent of $\varepsilon$.
The two types of nodes are labelled using a sign ±. Translating the lattice by a multiple of the spatial period ε maps nodes of the same type to one another. The initial position of the node of type ± in the pair labelled i, with 0≤i≤N, is given by

$$x_i^{\pm} = \epsilon (i \pm \frac{1}{4}). \qquad (3.1)$$

The nodes are connected to their nearest- and next-nearest neighbors by four different types of springs, each type being denoted by a Roman numeral, I, ..., IV. The spring constants ε⁻¹QI>0, ..., ε⁻¹QIV>0 of each type of spring are all assumed to scale with the scale separation parameter as ε⁻¹. This allows the effective properties of the equivalent bar to remain finite in the limit ε→0. The four scaled spring constants (QI, ..., QIV) can be chosen independently. It will be convenient to parameterize them using the elasticity parameters (Q̄, α, β, β') as

$$ \begin{align} \bar{Q} &= \frac{1}{Q_I^{-1} + Q_{II}^{-1} + Q_{III} + Q_{IV}} & Q_I &= \bar{Q}\alpha \frac{1}{(1-\beta)/2} \nonumber \\ a &= \left(1 + (Q_{III} + Q_{IV})\left(\frac{1}{Q_I} + \frac{1}{Q_{II}}\right)\right)^{-1} & Q_{II} &= \bar{Q}\alpha \frac{1}{(1+\beta)/2} \nonumber \\ \beta &= \frac{Q_I - Q_{II}}{Q_I + Q_{II}} & Q_{III} &= \bar{Q}(1-\alpha) \frac{1-\beta'}{2} \nonumber \\ \beta' &= \frac{Q_{IV} - Q_{III}}{Q_{III} + Q_{IV}} & Q_{IV} &= \bar{Q}(1-\alpha) \frac{1+\beta'}{2}. \tag{3.2} \end{align} $$

These lattice parameters are designed to yield the effective properties of the lattice using compact formulas. As we will see later, the first parameter Q̄, for instance, is nothing but the traction modulus of the equivalent bar model at leading order.

The original spring constants have to be positive, QI>0, ..., QIV>0. This constraints the elasticity parameters to be chosen in the ranges

$$\bar{Q}>0, \quad 0 \le \alpha \le 1, \quad -1 \le \beta \le +1, \quad -1 \le \beta' \le +1. \qquad (3.3)$$

We leave the elastic parameters of the lattice undetermined in most of the forthcoming analysis but from time to time, we will illustrate our findings using the particular values listed in Table 3.1. These values match those used in our previous work [TAL25b], which will make it easy to verify our results later in Section 10.1.

<table><caption>**Table 3.1.** Special values of spring constants used in the numerical illustrations (first three columns) and corresponding values of the elastic parameters in (3.2) (four last columns).</caption><thead><tr><th>Q<sub>II</sub></th><th>Q<sub>III</sub></th><th>Q<sub>IV</sub></th><th>Q&#773;</th><th>α</th><th>β</th><th>β'</th></tr></thead><tbody><tr><td>3Q<sub>I</sub></td><td>Q<sub>I</sub></td><td>2Q<sub>I</sub></td><td>15Q<sub>I</sub>/4</td><td>1/5</td><td>-1/2</td><td>+1/3</td></tr></tbody></table>

## 3.2. Strain

The lattice is one-dimensional and deforms along the x-axis when (longitudinal) forces are applied. The unknowns of the elastic problem are the set of displacements ŴI, where ± is the type of node and i is the node pair index, 0≤i≤N. We group the two displacements relative to a pair of nodes into a vector,

$$\mathbf{W}_i = \begin{pmatrix} W_i^- \\ W_i^+ \end{pmatrix}, \quad (0 \le i \le N). \qquad (3.4)$$

Next, we define the scaled change of length ÊIj of the spring of type ϕ in cell j by dividing the change of length by the cell length ε. As mentioned earlier, a cell contains two half-springs, which are of the same type I. The halves are denoted as I and I', respectively. The scaled elongations of the five types of springs appearing in a given cell j are grouped in a vector E_j, which we will refer to as the microscopic strain,

$$\mathbf{E}_j = (E_j^I E_j^{II} E_j^{III} E_j^{IV} E_j^{I'}) . \quad (0 \le j \le N-1), \qquad (3.5)$$

What we denote by ÊIj is the scaled change of length of the (full) spring of type I appearing on the left of the cell, and by ÊI'j that of the (full) spring of type I appearing on the right of the cell. The change of length of springs of type I is thus duplicated among adjacent cells, ÊI'j = ÊIj+1, which will be convenient.

The scaled changes of length collected in E_j can be expressed in terms of the displacements W_j and W_{j+1} in the adjacent pairs of nodes as

$$\mathbf{E}_j = \frac{1}{\epsilon} \mathcal{B} : (\mathbf{W}_j \mathbf{W}_{j+1}), \quad (0 \le j \le N-1) \qquad (3.6)$$

where the double dot denotes double contraction, see Section 2.4, the 1/ε coefficient is a conventional scaling factor introduced in the definition of E_j, and B is the scaled displacement-to-strain tensor, of size 5×2×2,

$$\mathcal{B} = \begin{pmatrix} (-1 & 0) \\ 0 & +1 \end{pmatrix} \begin{pmatrix} 0 & -1 \\ +1 & 0 \end{pmatrix} \begin{pmatrix} -1 & 0 \\ +1 & 0 \end{pmatrix} \begin{pmatrix} 0 & -1 \\ 0 & +1 \end{pmatrix} \begin{pmatrix} 0 & 0 \\ -1 & +1 \end{pmatrix}. \qquad (3.7)$$

In the tensor of rank 3 appearing in the right-hand side of (3.7), the successive indices refer to (i) a spring type, which can be I,II,III,IV or I', running column-wise, (ii) the relative position of an adjacent node (left for i=j versus right for i=j+1), running row-wise, and (iii) a node type (- first, + second), running column-wise. With our conventions, the first strain in (3.6), for instance, is given with the help of (3.7) as ÊIj = ε⁻¹(Β₁₁₁Ŵj + Β₁₁₂Ŵj⁺ + Β₁₂₁Ŵj+1 + Β₁₂₂Ŵj+1) = ε⁻¹(-Ŵj + Ŵj⁺).
A rigid-body displacement of the lattice with unit amplitude corresponds to identical nodal displacements $W_i=\mathcal{T}$ for all $i \ge 0$, where $\mathcal{T}$ represents the unit translation of a particular cell,

$$ \mathcal{T} = (1 \ 1). \qquad (3.8) $$

The invariance of strain by a rigid-body translation is expressed by the identity

$$ \mathcal{B}: (\mathcal{T} \ \mathcal{T}) = 0, \qquad (3.9) $$

which can be checked directly using (3.7): the sum of the entries of each of the five 2×2 sub-matrix appearing in the right-hand side is zero.

### 3.3. External loading

Two types of forces are applied to the nodes, representing a density of force $f(x)$ per unit length applied to all nodes on the one hand, and point-like forces $p^0$ and $p^\ell$ applied at the terminal nodes ($i=0,-$) and ($i=N,+$), respectively. The profile $f$ of distributed force and the terminal forces ($p^0,p^\ell$) are parameters of the problem. The nodal force relevant to the pair of nodes labelled $i$, with $0 \le i \le N$, is thus described by a vector $F_i$,

$$ F_i = \begin{pmatrix} p^0 \\ 0 \end{pmatrix} \delta_i + \frac{\epsilon}{2} \begin{pmatrix} f(x_i^-) \\ f(x_i^+) \end{pmatrix} + \begin{pmatrix} 0 \\ p^\ell \end{pmatrix} \delta_{i-N}, \qquad (3.10) $$

where $\delta_p$ is Kronecker's symbol defined by $\delta_0=1$ and $\delta_p=0$ if $p \ne 0$. The vector $F_i$ is filled using the same convention as the nodal displacements $W_i$ in (3.4), i.e., with the force on the node of type - coming first, and that on the node of type + coming second. The nodal force associated with the force density $f(x)$ could in principle be defined in different ways. In (3.10), we use the convention that the force at the node ($i, \pm$) is obtained by evaluating the continuous force profile $f(x)$ at the initial position $x_i^\pm$ of the node.

### 3.4. Energy

The energy of the lattice is the sum of the strain energy $\Phi_s$ stored in the springs and the potential energy $\Phi_e$ of the applied forces,

$$ \Phi_d(\mathbf{W}_0, \mathbf{W}_1, \dots, \mathbf{W}_N) = \Phi_s(\mathbf{E}_0 = \frac{1}{\epsilon}\mathcal{B}: (\mathbf{W}_0 \ \mathbf{W}_1), \dots, \mathbf{E}_{N-1} = \frac{1}{\epsilon}\mathcal{B}: (\mathbf{W}_{N-1} \ \mathbf{W}_N)) + \Phi_e(\mathbf{W}_0, \mathbf{W}_1, \dots, \mathbf{W}_N). \qquad (3.11) $$

The strain energy and potential energy of the applied forces are given, respectively, by

$$ \begin{align} \Phi_s(\mathbf{E}_0, \mathbf{E}_1, \dots, \mathbf{E}_{N-1}) &= \sum_{j=0}^{N-1} \frac{1}{2} (\epsilon \mathbf{E}_j) \cdot \frac{(\mathcal{Q}^{\text{bk}} + \delta_j \mathcal{Q}^0 + \delta_{N-1-j} \mathcal{Q}^\ell)}{\epsilon} \cdot (\epsilon \mathbf{E}_j) \tag{3.12} \\ \Phi_e(\mathbf{W}_0, \mathbf{W}_1, \dots, \mathbf{W}_N) &= -\sum_{i=0}^{N} \mathbf{F}_i \cdot \mathbf{W}_i, \notag \end{align} $$

where the matrices of spring constants are the three positive, symmetric, square matrices of size 5×5,

$$ \begin{align} \mathcal{Q}^{\text{bk}} &= \operatorname{diag}(1/2 \mathcal{Q}_I \ \mathcal{Q}_II \ \mathcal{Q}_III \ \mathcal{Q}_IV \ 1/2 \mathcal{Q}_I) \notag \\ \mathcal{Q}^0 &= \operatorname{diag}(1/2 \mathcal{Q}_I \ 0 \ 0 \ 0 \ 0) \notag \\ \mathcal{Q}^\ell &= \operatorname{diag}(0 \ 0 \ 0 \ 0 \ 1/2 \mathcal{Q}_I). \tag{3.13} \end{align} $$

In the strain energy (3.12)$_1$, we have restored the physical (unscaled) forms of both the change of length $\epsilon \mathbf{E}_j$ and of the spring constants $\epsilon^{-1} \mathcal{Q}^{\text{bk}} = \operatorname{diag}(\epsilon^{-1} \mathcal{Q}_I, \dots, \epsilon^{-1} \mathcal{Q}_IV)$. The right-hand side has been set up in such a way that

* springs of types II, III and IV are counted once in the unique cell $j$ to which they belong, with their actual spring constant;

* interior springs of type $I$ are attached to both two adjacent cells and are thus counted twice, but with half their actual spring constant each time, see the definition of $\mathcal{Q}^{\text{bk}}$ in (3.13)$_1$;

* the terminal springs of type $I$ are counted once, being part of the single adjacent cell; the Kronecker symbols in (3.12) and definitions of $\mathcal{Q}^{\text{bk}}$, $\mathcal{Q}^0$ and $\mathcal{Q}^\ell$ ensure that the correct spring constant $\epsilon^{-1}(1/2 \mathcal{Q}_I + 1/2 \mathcal{Q}_I) = \epsilon^{-1} \mathcal{Q}_I$ is used.

The contributions $1/2 \mathcal{Q}_I$ included in $\mathcal{Q}^0$ and $\mathcal{Q}^\ell$ thus correspond to the half-springs shown in black in Figure 3.1, that are 'missed' by the generic cell stiffness $\mathcal{Q}^{\text{bk}}$.

### 3.5. Equilibrium at interior nodes

In Appendix B, we derive the equilibrium equation of the lattice by making the energy $\Phi_d(\mathbf{W}_0, \mathbf{W}_1, \dots, \mathbf{W}_N)$ stationary. As we are interested in identifying a continuum approximation for the chain of springs by an energy method, we will not attempt to homogenize the discrete equilibrium equations, by contrast with the classical approach.
Throughout this paper, we will only need the discrete equilibrium equations applicable at an interior pair of nodes, to analyze boundary layers. They are derived in (B.4)$_2$ in the form

$$ \forall i \in [1, N-1], \quad \frac{1}{\epsilon} \mathcal{D}^{\text{bk}} \begin{pmatrix} \mathbf{W}_{i-1} \\ \mathbf{W}_i \\ \mathbf{W}_{i+1} \end{pmatrix} + \mathbf{F}_i = \mathbf{0}, \qquad (3.14) $$

where the 2×3×2 constant tensor $\mathcal{D}^{\text{bk}}$ is given in (B.3)$_2$ in terms of $\mathcal{B}$ and $\mathcal{Q}^{\text{bk}}$. By (3.14), the equilibrium of an interior pair of nodes is expressed in terms of the displacement in this pair and in the two adjacent pairs of nodes.

## 4. TWO-SCALE EXPANSION

Our goal is to derive a continuous approximation to the discrete energy (3.11–3.12) in the limit where the cell size $\epsilon$ is much smaller than the lattice length,

$$ \epsilon \ll 1. \qquad (4.1) $$

Specifically, our goal is to ensure accuracy to order $\epsilon^2$ included.

We approach the limit $\epsilon \to 0$ using a two-scale expansion, starting from the ansatz

$$ \mathbf{W}_i = \underbrace{\begin{pmatrix} u(x_i^-) \\ u(x_i^+) \end{pmatrix}}_{\mathbf{W}_i^{\text{bk}}} + \epsilon \underbrace{\begin{pmatrix} -y(x_i^-) \\ +y(x_i^+) \end{pmatrix}}_{\mathbf{W}_i^{\text{bk}}} + \epsilon \mathbf{Z}_i \qquad (4.2) $$

where $u$ is the macroscopic displacement, $y$ the microscopic shift and $Z_i$ is the boundary-layer corrector. Both $u$ and $y$ are smooth functions of $x$ in the sense of Equation (4.3) below, and $Z_i$ is a function of the index $i$ for a pair of node, decreasing exponentially quickly away from the endpoints, for both $i \gg 1$ and $N-i \gg \epsilon$.

The first part of the ansatz denoted by the under-brace $\mathbf{W}_i^{\text{bk}}$ concerns the bulk (hence the superscript ‘bk’) and is relevant to the homogenization procedure. With this ansatz, we have written the displacement of the node $(i, \pm)$ as $u(x) \pm \epsilon y(x)$ where $x = x_i^{\pm}$ is the reference position, see (3.4) and (3.1). The microscopic displacement $u(x)$ is thus the average of the nodal (microscopic) displacement over the two types of nodes, $\pm$. The microscopic shift $y(x)$ comes with a small, alternating coefficient $\pm \epsilon$, and captures the short-scale variations of the microscopic displacement. The main goal of homogenization is to express $y$ in terms of the derivatives of $u$, and to eliminate it.

The boundary layer $\epsilon Z_i$ arises as the solution predicted by homogenization, which ignores the endpoints, violates in general the equilibrium of the terminal nodes.

**Remark 4.1.** The coefficients $\epsilon$ included in factor of both the microscopic shift $y$ and of the boundary-layer corrector $Z_i$ in (4.2) ensure that $u, y$ and $Z_i$ all contribute to the microscopic strain at order $\epsilon^0$, as can be checked by using (3.6).

**Remark 4.2.** In our previous work based on matched-asymptotic expansions, we obtained a final form of the displacement which was precisely of the form (4.2), i.e., that included a boundary-layer displacement of the form $\epsilon Z_i$, see Equation [38] of [TAL24]. Matched asymptotic expansions introduce separate approximations of the displacement in the boundary layer on the one hand and in the bulk on the other hand, and we had to combine them using a composite approximation to arrive at Equation [38]. In the present work, we take the opposite route and use the composite approximation (4.2) as a starting point.

The functions $u(x)$ and $y(x)$ introduced in (4.2) do not only depend on their argument $x$ (macroscopic coordinate) but also on the scale separation parameter $\epsilon$. The dependence on $\epsilon$ is implicit in our notation. We require that this dependence can be captured by a regular series expansion,

$$ u(x) = \sum_{j=0}^{\infty} \epsilon^j u_{(j)}(x), \quad y(x) = \sum_{j=0}^{\infty} \epsilon^j y_{(j)}(x). \qquad (4.3) $$

where $u_{(j)}(x)$ and $y_{(j)}(x)$ are independent of $\epsilon$. Thus, we do not allow $u$ or $y$ to be functions of the fast variable $x/\epsilon$, for instance. By contrast, the boundary-layer corrector $Z_i$ depends on the spatial index $i$, hence varies on a short length scale $\sim \epsilon$. It is also implicitly a series in $\epsilon$ and will later be expressed in terms of the series $u'$ and $u''$, see Equation (7.6).

Following Bakhvalov [BP89], we always keep the series appearing in (4.3) hidden inside the notation $u(x)$ and $y(x)$. Doing so will lead us directly to a compact expression of the energy, avoiding the detour taken in most of the existing work on higher-order homogenization, whereby a series expansion of the energy (or of the equations of equilibrium) is first obtained, with similar terms repeated at the different orders, and these redundancies are next eliminated using resummation.

The only practical consequence of the ‘smoothness’ of the functions $u$ and $y$ embodied in the series expansions (4.3) is that the successive derivatives of $u$ and $y$ are all of order $\mathcal{O}(\epsilon^0)$, like the original functions $u$ and $y$:

$$ u = \mathcal{O}(\epsilon^0), \quad u' = \mathcal{O}(\epsilon^0), \quad u'' = \mathcal{O}(\epsilon^0), \quad \dots, \quad y = \mathcal{O}(\epsilon^0), \quad y' = \mathcal{O}(\epsilon^0), \quad y'' = \mathcal{O}(\epsilon^0), \quad \dots \qquad (4.4) $$

This follows by differentiating formally the right-hand sides of (4.3) term by term. We will repeatedly use these scaling assumptions in the forthcoming calculations. They are at the heart of the method of formal series expansions which we adopt here.
The next step of the procedure is to solve for the microscopic shift y using homogenization. This is done in Section 5. The corrector $Z_i$ will be treated in two steps: using boundary-layer analysis, we will obtain it in Section 6 up to two amplitudes $(v^0, v^\ell)$, which we will next set by energy minimization in Section 7. In the final sections, we will thus be able to express the energy in terms of the macroscopic displacement $u$ only.

# 5. HOMOGENIZATION IN THE BULK

In this section, we limit attention to nodes that are in the bulk, i.e., well inside the lattice,

$$i \gg 1 \quad \text{and} \quad N-i \gg 1. \qquad (5.1)$$

It is then possible to ignore the boundary-layer contribution $\epsilon Z_i$ in the ansatz (4.2), and we obtain the form $W_i = W_i^{bk}$ of the ansatz relevant to homogenization,

$$W_i^{bk} = \begin{pmatrix} u(x_i^\tau) \\ u(x_i^\ell) \end{pmatrix} + \epsilon \begin{pmatrix} -y(x_i^\tau) \\ +y(x_i^\ell) \end{pmatrix}. \qquad (5.2)$$

We solve the microscopic shift $y$ in terms of the derivatives of $u(x)$. This part of the paper follows the standard procedure in second-homogenization as no energy argument is required yet. We adopt a concise presentation style and invite the reader to refer to our previous work [TAL24] for details.

## 5.1. Microscopic shift

We anticipate that the microscopic shift can be found as a series expansion in terms of the successive gradients of the macroscopic strain $u'$,

$$y(x) = \mathcal{Y}_0 u'(x) + \epsilon \mathcal{Y}_1 u''(x) + \epsilon^2 \mathcal{Y}_2 u'''(x) + \dots \qquad (5.3)$$

The coefficients $\mathcal{Y}_i$ are found by inserting (5.3) into the nodal displacement $W_i^{bk}$ in (5.2), and then into the discrete equilibrium in the bulk in (3.14) relevant to the pair of nodes indexed by $i$. In the companion notebook [TA25], we do a Taylor expansion of the equilibrium equation thus obtained, in powers of $\epsilon$, about the center $x=\epsilon i$ of the pair, and solve it order by order. This yields the first coefficients as

$$\mathcal{Y}_0 = -\frac{\beta}{4}, \qquad \mathcal{Y}_1 = \frac{(1-\alpha)(1-\beta^2)\beta'}{16\alpha}. \qquad (5.4)$$

In the process, we also obtain the leading orders of the homogenized equilibrium equations such as $\overline{Qu}_0'(x) + f(x) = 0$, by averaging the discrete equilibrium (3.14) over the node types $\pm$, but we will ignore them for the moment. We will derive them again later by an energy method.

Anticipating that the higher-order constants $\mathcal{Y}_2, \mathcal{Y}_3, \dots$ enter at order $\epsilon^3$ and beyond in the energy, we make no attempt to determine them.

In the two coming subsections, we insert the microscopic solution just obtained into the expressions of the two fundamental quantities entering into the discrete energy, namely the nodal displacement and the microscopic strain.

## 5.2. Microscopic displacement

Given an index $i$ of a pair of nodes in the bulk, we insert the solution (5.3–5.4) for the microscopic shift into (5.2) and do a Taylor expansion about the center $x_i=\epsilon i$ of the pair. Using the vector-of-derivative notation introduced in (2.1), we get the displacement $W_i^{bk}$ in the form

$$W_i^{bk} = u(\epsilon i) \mathcal{T} + \epsilon \mathcal{G}(\mathcal{Y}_2) \cdot \dot{u}'(\epsilon i) + \mathcal{O}(\epsilon^4) \qquad (5.5)$$

where the leading-order term $\mathcal{T}$ is the rigid translation identified earlier in (3.8), and the $2 \times 3$ matrix $\mathcal{G}$ collects the coefficients produced by the expansion. With the help of symbolic calculations, it is calculated in the companion notebook [TA25] in terms of the lattice properties as

$$\mathcal{G}(\mathcal{Y}_2) = \begin{pmatrix} \frac{1-\beta}{4} & \cdots & \cdots \\ +\frac{1-\beta}{4} & \cdots & \cdots \end{pmatrix} + \begin{pmatrix} -1 \\ +1 \end{pmatrix} \otimes (0 \ 0 \ \mathcal{Y}_2). \qquad (5.6)$$

We decided not to calculate the constant $\mathcal{Y}_2$ as it is not available in the context of second-order homogenization. The quantity $\mathcal{G}(\mathcal{Y}_2)$ thus remains partially indeterminate. We cannot simply discard $\mathcal{Y}_2$ by setting it to zero forcefully, as this would introduce a larger error than the remainder $\mathcal{O}(\epsilon^4)$ in (5.5). We thus keep $\mathcal{Y}_2$ undetermined, treating it as a pure symbol in the companion notebook. We will check that it disappears from the final results.

Using the identity (2.6), we can degrade the approximation just obtained as $W_i^{bk} = \mathcal{L}^{bk} \cdot \dot{u}(\epsilon i) + \mathcal{O}(\epsilon^3)$, where

$$\mathcal{L}^{bk} = (\mathcal{T} \otimes (1 \ 0 \ 0) + \mathcal{G}(\mathcal{Y}_2) \cdot \mathcal{S}). \qquad (5.7)$$

Note that $\mathcal{L}^{bk}$ does not depend on the higher-order homogenization constant $\mathcal{Y}_2$ as $(0 \ 0 \ \mathcal{Y}_2) \cdot \mathcal{S} = \mathbf{0}$ by (2.7).

To sum up, we have obtained two approximations of the bulk displacement, one valid at order $\epsilon^3$ and another one valid at order $\epsilon^2$, in the form

$$W_i^{bk} = \begin{cases} u(\epsilon i) \mathcal{T} + \epsilon \mathcal{G}(\mathcal{Y}_2) \cdot \dot{u}'(\epsilon i) + \mathcal{O}(\epsilon^4) \\ \mathcal{L}^{bk} \cdot \dot{u}(\epsilon i) + \mathcal{O}(\epsilon^3). \end{cases} \qquad (5.8)$$
The approximation on the first line will be required to resolve the strain with an error $\mathcal{O}(\epsilon^3)$, as the definition of the strain in (3.6) involves a division by $\epsilon$; this is done in the next section. The lower-order approximation on the second line will be sufficient to approximate the potential of the external loading, see (3.12)$_2$; this is done in Appendix D.2.

**Remark 5.1.** We use the scaling laws (4.4) repeatedly, albeit implicitly. The expansion $W_i^{bk}=\mathcal{L}^{bk}\ddot{u}(\epsilon i)+\mathcal{O}(\epsilon^3)$ obtained in (5.8)$_2$, for instance, can be spelled out as

$$W_i^{bk} = \left( \mathcal{L}_{21}^{11} \right) u(\epsilon i) + \epsilon \left( \mathcal{L}_{22}^{12} \right) u'(\epsilon i) + \epsilon^2 \left( \mathcal{L}_{23}^{13} \right) u''(\epsilon i) + \mathcal{O}(\epsilon^3),$$

and the successive terms in the right-hand side are respectively $\mathcal{O}(\epsilon^0)$, $\mathcal{O}(\epsilon^1)$ and $\mathcal{O}(\epsilon^2)$ by (4.4). Besides, the series expansion (4.3), which we normally keep hidden, yields the explicit expression of the contributions to $W_i^{bk}$ of order $\epsilon^0$, $\epsilon^1$ and $\epsilon^2$ as

$$W_i^{bk} = \left( \mathcal{L}_{21}^{11} \right) u_{(0)}(\epsilon i) + \epsilon \left( \left( \mathcal{L}_{22}^{12} \right) u'_{(0)}(\epsilon i) + \left( \mathcal{L}_{21}^{11} \right) u_{(1)}(\epsilon i) \right) \\ + \epsilon^2 \left( \left( \mathcal{L}_{23}^{13} \right) u''_{(0)}(\epsilon i) + \left( \mathcal{L}_{22}^{12} \right) u'_{(1)}(\epsilon i) + \left( \mathcal{L}_{21}^{11} \right) u_{(2)}(\epsilon i) \right) + \mathcal{O}(\epsilon^3).$$

This cumbersome expression is only intended to show the benefit of using the compact vector-of-derivatives notation, and of treating $u(x)$ implicitly as a series in $\epsilon$. It will not be used.

## 5.3. Microscopic strains

Inserting the microscopic displacement $W_i^{bk}$ in (5.8)$_1$ into the microscopic strain $E_j$ in (3.6), we obtain the microscopic strain $E_j=E_j^{bk}$ ‘in the bulk’ as

$$E_j^{bk} = \frac{1}{\epsilon} \mathcal{B} : \left( \begin{array}{c} u(\epsilon j) \mathcal{T} \\ u(\epsilon(j+1)) \mathcal{T} \end{array} \right) + \mathcal{B} : \left( \begin{array}{c} \mathcal{G}(\mathcal{Y}_2) \cdot \ddot{u}'(\epsilon j) \\ \mathcal{G}(\mathcal{Y}_2) \cdot \ddot{u}'(\epsilon(j+1)) \end{array} \right) + \mathcal{O}(\epsilon^3). \quad (5.9)$$

By ‘in the bulk’, we mean that this expression ignores the boundary-layer corrector $\epsilon Z_i$ which we will restore later.

Using the invariance by translation in (3.9) and the Taylor expansion (2.9) around $x=\epsilon(j+1/2)$, with the shift parameter set to $\alpha=\pm 1/2$, we rewrite the bulk microscopic strain as

$$E_j^{bk} = \frac{u(\epsilon(j+1)) - u(\epsilon j)}{\epsilon} \mathcal{B} : \left( \begin{array}{c} -\mathcal{T}/2 \\ +\mathcal{T}/2 \end{array} \right) + \mathcal{B} : \left( \begin{array}{c} \mathcal{G}(\mathcal{Y}_2) \cdot \exp(-\mathcal{S}/2) \\ \mathcal{G}(\mathcal{Y}_2) \cdot \exp(+\mathcal{S}/2) \end{array} \right) \ddot{u}' \left( \epsilon \left(j + \frac{1}{2}\right) \right) + \mathcal{O}(\epsilon^3). \quad (5.10)$$

Combining with the Taylor expansion $u(\epsilon(j+1))-u(\epsilon j)=\epsilon\left(u'+\frac{\epsilon^2}{24}u''\right)\left(\epsilon\left(j+\frac{1}{2}\right)\right)=(1\ 0\ \frac{1}{24})\ddot{u}'(\epsilon\left(j+\frac{1}{2}\right))+\mathcal{O}(\epsilon^3)$, we obtain

$$E_j^{bk} = \mathcal{E}^{bk}(\mathcal{Y}_2) \cdot \ddot{u}' \left( \epsilon \left( j + \frac{1}{2} \right) \right) + \mathcal{O}(\epsilon^3), \quad (5.11)$$

where

$$\mathcal{E}^{bk}(\mathcal{Y}_2) = \mathcal{B} : \left( \left( \begin{array}{c} -\mathcal{T}/2 \\ +\mathcal{T}/2 \end{array} \right) \otimes \left( \begin{array}{cc} 1 & 0 \\ 0 & \frac{1}{24} \end{array} \right) + \left( \begin{array}{c} \mathcal{G}(\mathcal{Y}_2) \cdot \exp(-\mathcal{S}/2) \\ \mathcal{G}(\mathcal{Y}_2) \cdot \exp(+\mathcal{S}/2) \end{array} \right) \right). \quad (5.12)$$

An explicit expression of $\mathcal{E}^{bk}(\mathcal{Y}_2)$ can be obtained by inserting the expression (5.6) of $\mathcal{G}(\mathcal{Y}_2)$,

$$\mathcal{E}^{bk}(\mathcal{Y}_2) = \begin{pmatrix} \frac{1-\delta}{2} & \cdots & \cdots \\ \frac{1+\delta}{2} & \cdots & \cdots \\ 1 & \frac{1-\delta}{4} & \cdots \\ 1 & \frac{1+\delta}{4} & \cdots \\ \frac{1-\delta}{2} & \cdots & \cdots \end{pmatrix} + \begin{pmatrix} +2 \\ -2 \\ 0 \\ 0 \\ +2 \end{pmatrix} \otimes (0 \ 0 \ \mathcal{Y}_2). \quad (5.13)$$

It still depends on the unknown constant $\mathcal{Y}_2$ which will ultimately disappear from the energy, see Equation (D.4) in the Appendix.

By the notation $\ddot{u}'$ introduced in (2.1), Equation (5.11) yields the change of length of the different types of springs as an expansion in the successive derivatives of the macroscopic strain $u'$. The elongation of the spring of type II attached to cell j, for instance, is given by $E_j^{II} = \mathcal{E}_{21}^{bk} u'(x^j) + \epsilon \mathcal{E}_{22}^{bk} u''(x^j) + \epsilon^2 \mathcal{E}_{23}^{bk} u'''(x^j) + \mathcal{O}(\epsilon^3)$ where $x^j = \epsilon(j + \frac{1}{2})$ is the center of the cell (see Figure 3.1). The rows in $\mathcal{E}^{bk}$ thus correspond to the 5 different types of springs, and the columns to the macroscopic strain $u'$ and its scaled gradients $\epsilon u''$ and $\epsilon^2 u'''$.

By using the solution (5.3) for the microscopic shift, we have obtained an approximation for the displacement $W_i^{bk}$ in (5.8) and for the strain $E_j^{bk}$ in (5.11). Having fully characterized the solution in the bulk, we proceed in the next section to address the boundary layers.

# 6. BOUNDARY LAYERS

In this section, we derive a solution for the boundary layer $Z_i$ in the ansatz (4.2). This complements the results of the previous section where we analyzed a slowly varying solution at a distance from the boundaries, for $i \gg 1$ and $N-i \gg 1$, see (5.1). As the boundary layer varies significantly from one cell to the next, we go back to the discrete equations of equilibrium. The solutions we derive ultimately decay exponentially quickly for $i \gg 1$ and $N-i \gg 1$.
The boundary-layer solutions are a correction of order $\epsilon$ to the displacement, see the ansatz (4.2), but will contribute to the strain at order $\epsilon^0$, see Equation (6.14) below. Being significant in a layer of width $\sim \epsilon$, they will thus contribute to the energy at order $\epsilon$. It is not possible to achieve the desired $\epsilon^2$-accuracy without taking them into account.

## 6.1. Boundary-layer solution at interior pairs of nodes

We start by solving the discrete equations of equilibrium at interior pairs of nodes, by limiting the pair index $i$ to the range

$$1 \le i \le N-1. \tag{6.1}$$

Doing so enables us to ignore for the moment (i) the special form of the equilibrium equations arising from the exceptions on the topology of the lattice at $i=0,N$, see (B.4)$_{1,3}$, and (ii) the point-like contributions $p^b$ in the loading term $\mathbf{F}_i$, see (3.10). We will properly take them into account later—the boundary-layer solution precisely arises because of the exceptions in the lattice topology and in the applied forces.

The only relevant term in the loading vector $\mathbf{F}_i$ at interior pairs of nodes ($1 \le i \le N-1$) is thus the distributed force (f-term), see (3.10). This is a slowly-varying term, which therefore affects the slowly varying (homogenized) part of the solution $\mathbf{W}^{\text{bk}}$ in (4.2). As far as the analysis of the boundary layer is concerned, we can therefore ignore the external loading entirely, without any loss in generality.

We will therefore identify the boundary-layer corrector $\mathbf{Z}_i$, by solving the homogeneous, discrete equation of equilibrium (3.14) at interior pairs of nodes,

$$\forall i \in [1, N-1], \quad \mathcal{D}^{\text{bk}} \begin{pmatrix} \mathbf{Z}_{i-1} \\ \mathbf{Z}_i \\ \mathbf{Z}_{i+1} \end{pmatrix} = \mathbf{0}. \tag{6.2}$$

This is a linear recurrence equation with constant coefficients. We seek exponential solutions of the form

$$\mathbf{Z}_i = \mu^i \tilde{\mathcal{Z}}_\mu, \tag{6.3}$$

where $\mu$ is a complex eigenvalue and $\tilde{\mathcal{Z}}_\mu \in \mathbb{R}^2$ is an eigenvector. Inserting into (6.2), we obtain the eigenvalue problem

$$\tilde{\mathcal{D}}(\mu) \cdot \tilde{\mathcal{Z}}_{\mu} = \mathbf{0}, \quad \text{where } \tilde{\mathcal{D}}_{i,j}(\mu) = \mu^{-1} \mathcal{D}_{i1j}^{\text{bk}} + \mathcal{D}_{i2j}^{\text{bk}} + \mu \mathcal{D}_{i3j}^{\text{bk}}. \tag{6.4}$$

The eigenvalues $\mu$ are found by solving the characteristic equation $\det \tilde{\mathcal{D}}(\mu)=0$, which is rewritten with the help of symbolic calculations in the companion notebook [TA25] as

$$(1 - \mu^2)(\mu^2 - 2(1+\beta)\mu + 1) = 0. \tag{6.5}$$

Here, $\beta$ is the real constant

$$\beta = \frac{8a}{(1-a)^2(1-b^2)(1-b'^2)} \ge 0. \tag{6.6}$$

It is positive by (3.3).

The characteristic equation has real roots $\mu$ only, namely the double root $\mu=1$ and two simple roots $\mu=z$ and $\mu=1/z$, where

$$z = 1 + \beta + \sqrt{2\beta + \beta^2} > 1. \tag{6.7}$$

* The double root $\mu=1$ yields two modes that are slow varying and, thus, already captured by the homogenized solution:

* a rigid translation mode, $\mathbf{Z}_i = \mathcal{T}$ for all $i \ge 0$, which can also be obtained by setting $u(x)=1$ in (5.8)$_1$;

* the mode with uniform traction, $\mathbf{Z}_i = i \mathcal{T} + \mathcal{G}(y_2) \cdot (1 \ 0 \ 0)$ which can also be obtained by setting $u(x)=x$ in (5.8)$_1$, up to a conventional factor $\epsilon$;

* In addition, we have two simple roots

* The root $\mu=1/z \in (0,1)$ yields a boundary-layer that decreases exponentially away from the first boundary, $i=0$, as $\mu^i \propto z^{-i}$.

* The root $\mu=z \in (1,\infty)$ yields a boundary-layer that decreases exponentially away from the second boundary, $i=N$, as $\mu^i \propto z^{-(N-i)}$.

In the companion notebook, the eigenvalues associated with each one of the simple roots are obtained in closed form,

$$\begin{align}
\mathcal{Z}^0 &= \tilde{\mathcal{Z}}_{(z^{-1})} = \begin{pmatrix} 2za((1+b)+z(1-b)) \\ (4+\beta)za - ((1-b')(1-(1-a)b^2)+a b')(\frac{1-z}{2})^2 \end{pmatrix}, \tag{6.8} \\
\mathcal{Z}^\ell &= \tilde{\mathcal{Z}}_{(z^{+1})} = \begin{pmatrix} (4+\beta)za - ((1+b')(1-(1-a)b^2)-a b')(\frac{1-z}{2})^2 \\ 2za((1+b)+z(1-b)) \end{pmatrix}. \nonumber
\end{align}$$

Discarding the two slowly-varying modes ($\mu=1$), we conclude that the boundary layer solution can be written as

$$\forall i \in [1, N-1], \quad \mathbf{Z}_i = v_{\text{exp}}^0 z^{-i} \mathcal{Z}^0 + v_{\text{exp}}^\ell z^{-(N-i)} \mathcal{Z}^\ell \tag{6.9}$$
in terms of two amplitudes $v_{\text{exp}}^0$ and $v_{\text{exp}}^\ell$, which will remain undetermined until we solve the equilibrium of the terminal pairs of nodes later.

**Remark 6.1.** The eigenvector $\mathcal{Z}^\ell$ in (6.8) describing the boundary layer solution at $b=\ell$ can be obtained by swapping the entries of the eigenvector $\mathcal{Z}^0$ describing the boundary layer at $b=0$, and by changing $\ell'$ to $-\ell'$. This is a consequence of the invariance of the lattice by mirror-symmetry, which exchanges the node types $\pm$ (hence the swap in the entries of vectors) and exchanges the springs types $Q_{III} \leftrightarrow Q_{IV}$ (hence the change of sign in $\ell'$, see (3.2)). By a similar argument, $z$ in (6.7) remains unchanged when the sign of $\ell'$ is changed, as can be shown using (6.6).

## 6.2. Adding degrees of freedom for the terminal pairs of nodes

The boundary-layer solution (6.9) could be extended to the terminal pairs of nodes $i=0,N$ by solving the special form of the discrete equilibrium equations applicable there, see (B.4)$_{1.3}$. We use a different but equivalent approach: we simply parameterize the boundary-layer displacements $Z_0$ and $Z_N$ using two new unknowns each, which we denote as $(v_{-\ell}^0, v_{+\ell}^0)$ and $(v_{-\ell}^\ell, v_{+\ell}^\ell)$:

$$ \forall i \in [0, N], \quad Z_i = v_{\text{exp}}^0 z^{-i} \mathcal{Z}^0 + \begin{pmatrix} v_{-\ell}^0 \\ v_{+\ell}^0 \end{pmatrix} \delta_i + v_{\text{exp}}^\ell z^{-(N-i)} \mathcal{Z}^\ell + \begin{pmatrix} v_{-\ell}^\ell \\ v_{+\ell}^\ell \end{pmatrix} \delta_{N-i}. \quad (6.10) $$

In the right-hand side of (6.10), the Kronecker symbols $\delta_i$ introduced in (3.10) ensure that the solution (6.9) is recovered at interior pairs of nodes, for $1 \le i \le N-1$. We will determine the boundary-layer displacements $v_{\pm}^b$ of the terminal nodes later, by energy minimization. As stated above, this is equivalent to solving the equilibrium equation (B.4)$_{1.3}$. It is also simpler, given that we will have to set up energy minimization anyways to determine the global boundary-layer amplitudes $v_{\text{exp}}^b$.

Collecting the boundary-layer parameters relevant to either endpoint $b \in \{0, \ell\}$ into a vector

$$ \mathbf{v}^b = (v_{\text{exp}}^b \ v_{-\ell}^b \ v_{+\ell}^b) \in \mathbb{R}^3, \quad (6.11) $$

we rewrite (6.10) in compact form as

$$ \forall i \in [0, N], \quad Z_i = v^0 \cdot \mathcal{L}_i^0 + v^\ell \cdot \mathcal{L}_i^\ell \quad (6.12) $$

where $\mathcal{L}_i^b$ are known 3×2 matrices depending on the nodal pair index $i \in [0, N]$ and endpoint $b \in \{0, \ell\}$,

$$ \mathcal{L}_i^0 = \begin{pmatrix} z^{-i} \mathcal{Z}^0 \\ \delta_i \mathbf{I}_2 \end{pmatrix}, \quad \mathcal{L}_i^\ell = \begin{pmatrix} z^{-(N-i)} \mathcal{Z}^\ell \\ \delta_{N-i} \mathbf{I}_2 \end{pmatrix}. \quad (6.13) $$

The boundary layer amplitudes $\mathbf{v}^b$ will ultimately be solved in terms of the macroscopic displacement $u$ and its gradients, which are formal series in powers of $\varepsilon$ by (4.3). In Equation (6.12) above, the $Z_i$'s are therefore series in powers of $\varepsilon$ as well, as we have announced earlier.

## 6.3. Microscopic strain arising from the boundary layer

By Equations (3.6) and (4.2), the boundary-layer solution (6.12) contributes to the microscopic strain $E_j$ as

$$ \mathcal{E}_j^0 \cdot \mathbf{v}^0 + \mathcal{E}_j^\ell \cdot \mathbf{v}^\ell + \mathcal{O}(\varepsilon^3), \quad (6.14) $$

where the 5×3 matrices $\mathcal{E}_j^b$ are given in terms of $j \in [0, N-1]$ and $b \in \{0, \ell\}$ as

$$ \mathcal{E}_j^0 = \mathcal{B} \cdot \begin{pmatrix} (\mathcal{L}_j^0)^\top \\ (\mathcal{L}_{j+1}^0)^\top \end{pmatrix}, \quad \mathcal{E}_j^\ell = \mathcal{B} \cdot \begin{pmatrix} (\mathcal{L}_j^\ell)^\top \\ (\mathcal{L}_{j+1}^\ell)^\top \end{pmatrix}. \quad (6.15) $$

The expressions (6.12) and (6.14) of the boundary-layer displacement and strain are used in the next section to derive a continuous approximation of the energy accounting for them in an effective way.

# 7. CONTINUOUS PSEUDO-ENERGY

In this Section, we combine the two parts of the solution obtained previously in Sections 5 and 6. By the ansatz (4.2), this yields the microscopic displacement $W_i$ in terms of a macroscopic displacement $u$ and boundary-layer parameters $\mathbf{v}^b$. We insert it into the discrete energy and derive a continuous approximation of the energy functional in terms of $u$ and $\mathbf{v}^b$. We call pseudo-energy this formal expansion of the energy in powers of $\varepsilon$.

## 7.1. Expressions of displacement and strain to be inserted into discrete energy

Inserting the homogenization solution (5.8)$_2$ obtained in Section 5 and the boundary-layer solution (6.12) for $Z_i$ obtained in Section 6 into the ansatz (4.2), we obtain the microscopic displacement in the explicit form

$$ W_i = \mathcal{L}^{\text{bk}} \cdot \ddot{u} (\varepsilon i) + \varepsilon (\mathbf{v}^0 \cdot \mathcal{L}_i^0 + \mathbf{v}^\ell \cdot \mathcal{L}_i^\ell) + \mathcal{O}(\varepsilon^3). \quad (7.1) $$

By the principle of linear superposition, the strain $E_j$ in (3.6) is then the sum of the strain (5.11) predicted by homogenization and the contribution (6.14) coming from the boundary layer,

$$ E_j = \mathcal{E}^{\text{bk}}(Y_2) \cdot \ddot{u}' \left(\varepsilon \left(j + \frac{1}{2}\right)\right) + \mathcal{E}_j^0 \cdot \mathbf{v}^0 + \mathcal{E}_j^\ell \cdot \mathbf{v}^\ell + \mathcal{O}(\varepsilon^3). \quad (7.2) $$
**Remark 7.1.** It is not possible to derive (7.2) directly from (7.1). Indeed, because of the division by $\epsilon$ in (3.6), an $\mathcal{O}(\epsilon^3)$ accuracy in the strain $E_j$ requires an $\mathcal{O}(\epsilon^4)$ accuracy in the displacement $W_i$, which is more than what (7.1) can offer. We have instead derived (7.2) from the $\mathcal{O}(\epsilon^4)$-accurate approximation of the bulk displacement (5.8)$_1$ and from the boundary-layer solution (6.12), which is valid at all orders—recall that $v^b$ is implicitly a formal series in $\epsilon$.

## 7.2. Energy approximation in terms of boundary-layer amplitudes

In Appendix D, we insert the relocalization formulas (7.1) and (7.2) of the displacement and strain into the discrete energy $\Phi_d$ in (3.11)–(3.12), and derive a continuous approximation $\tilde{\Phi}_\epsilon$ in terms of the macroscopic displacement $u$ and of the boundary-layer amplitudes $v^b$, see (D.28),

$$ \tilde{\Phi}_{\epsilon}[u, v^0, v^\ell] = \int_0^\ell \left( \frac{1}{2} \ddot{u}' \cdot \bar{\mathbf{K}} \cdot \ddot{u}' - (\bar{\mathbf{H}} \cdot \ddot{\mathbf{f}}) u \right) dx \\ + \sum_{b \in \{0, \ell\}} \left( \frac{\epsilon}{2} \dot{u}'(b) \cdot \bar{\mathbf{k}}^b \cdot \dot{u}'(b) - \left( \bar{\mathbf{h}}^b \cdot \left( \begin{pmatrix} p^b \\ \epsilon \dot{f}(b) \end{pmatrix} \right) \right) \cdot \dot{u}(b) \right. \\ \left. + \epsilon \left( \frac{1}{2} v^b \cdot \bar{\mathbf{g}}^b \cdot v^b - v^b \cdot \left( \begin{pmatrix} \bar{\mathcal{J}}_u^b \\ \bar{\mathcal{J}}_f^b \end{pmatrix} \right) \cdot \left( \begin{pmatrix} \dot{u}'(b) \\ \epsilon f(b) \end{pmatrix} \right) \right) \right) + \mathcal{O}(\epsilon^3). \quad (7.3) $$

In the right-hand side, $\bar{\mathbf{K}}$ is the 3×3, symmetric matrix of homogenized elastic moduli, see (D.4), $\bar{\mathbf{H}} \in \mathbb{R}^3$ is the 'force condensation' vector $\bar{\mathbf{H}} \in \mathbb{R}^3$ that converts the distributed force and force gradients $\ddot{\mathbf{f}} = (f \ \epsilon f' \ \epsilon^2 f'')$ into an effective distributed force $\bar{\mathbf{H}} \cdot \ddot{\mathbf{f}}$, see (D.27)$_1$, $\bar{\mathbf{k}}^b$ are 2×2 symmetric matrices capturing an effective boundary elasticity at the endpoints, the 3×3 'traction condensation' matrices $\bar{\mathbf{h}}^b$ convert the boundary load into effective tractions, the symmetric 3×3 matrix $\bar{\mathbf{g}}^b \ge 0$ captures a stiffness effect associated with the boundary layer, and the 3×2 matrices $\bar{\mathbf{J}}_u^b$ and $\bar{\mathbf{J}}_f^b$ capture the coupling of the boundary layer with the macroscopic displacement and with the applied load, respectively.

The derivation of the continuous approximation (7.3) in Appendix D uses the Euler–Maclaurin formula to transform discrete sums into integrals, and involves a number of Taylor expansions. Boundary terms are produced along the way, both by the Euler–Maclaurin formula and by integrations by parts. We work them out carefully to ensure order-$\epsilon^2$ accuracy.

What we analyze in Appendix D is in fact a semi-infinite lattice ($N=\infty$). Appendix D delivers a special form of the pseudo-energy (7.3) having only one boundary term, $b=0$, see (D.28). We have extended it in (7.3) to two endpoints. By proceeding so, we neglect the exponentially small coupling between the boundary layers at the opposite ends.

With the particular definition (3.10) of the microscopic loading used in the present work, the force condensation matrix $\bar{\mathbf{H}}$ predicted by (D.27)$_1$ takes the simple value

$$ \bar{\mathbf{H}} = (1 \ 0 \ 0). \quad (7.4) $$

We will accordingly simplify the external potential in the bulk as $\int_0^\ell - (\bar{\mathbf{H}} \cdot \ddot{\mathbf{f}}) u dx = \int_0^\ell - f u dx$, starting in Equation (7.7) below.

## 7.3. Pseudo-energy

Next, we minimize $\tilde{\Phi}_\epsilon$ in (7.3) with respect to the boundary-layer amplitudes $v^b$, by solving the linear set of equations

$$ \forall b \in \{0, \ell\}, \quad \frac{\partial \tilde{\Phi}_\epsilon}{\partial v^b}[u, v^0, v^\ell] = 0. \quad (7.5) $$

This is equivalent to solve the equilibrium equations for the terminal pairs of nodes, $i=0, N$. The linear problems for $v^0$ and $v^b$ are decoupled, and make use of the symmetric, positive square matrix of coefficients $\bar{\mathbf{g}}^b \ge 0$, see (7.3) and (D.15). The minimization delivers:

* The solution for the boundary-layer amplitudes $v^b$ in terms of the macroscopic displacement $u(b)$, its gradient $u'(b)$, the boundary load $(p^b, \epsilon f(b))$, and a known 3×4 matrix $\bar{\mathcal{L}}^b$,

$$ v^b = \bar{\mathcal{L}}^b \cdot \left( \frac{\dot{u}'(b)}{p^b \epsilon f(b)} \right) + \mathcal{O}(\epsilon^2), \quad \text{where } \bar{\mathcal{L}}^b = (\bar{\mathbf{g}}^b)^{-1} \cdot \left( \bar{\mathcal{J}}_u^b \mid \bar{\mathcal{J}}_f^b \right). \quad (7.6) $$

Note that the $\mathcal{O}(\epsilon^2)$ error on the boundary-layer amplitude $v^b$ is larger than the $\mathcal{O}(\epsilon^3)$ in the energy: $v^b$, like $Z_i$, is systematically multiplied by $\epsilon$, see for instance (4.2). In Equation (7.6), the horizontal and vertical bars are designed to help interpreting the block-matrix notation: the blocks separated by the vertical bar in the definition of $\bar{\mathcal{L}}^b$ are respectively coupled to the parts of the vector separated by the horizontal bar.

* An energy approximation $\tilde{\Phi}_\epsilon[u]$ depending on the macroscopic displacement $u$ only, obtained by inserting the optimal boundary-layer amplitudes (7.6) into (7.3), see (F.3) and (F.7) in Appendix F,

$$ \tilde{\Phi}_{\epsilon}[u] = \int_{0}^{\ell} \left( \frac{1}{2} \ddot{u}' \cdot \bar{\mathbf{K}} \cdot \ddot{u}' - f u \right) dx + \sum_{b \in \{0, \ell\}} \left( \frac{\epsilon}{2} \dot{u}'(b) \cdot \bar{\mathbf{k}}^b \cdot \dot{u}'(b) \right. \\ \left. - \left( \bar{\mathbf{h}}^b \cdot \left( \begin{pmatrix} p^b \ \epsilon f(b) \end{pmatrix} \right) u(b) \right. \right. \\ \left. \left. + \epsilon \left( \bar{\mathbf{r}}^b \cdot \left( \begin{pmatrix} p^b \ \epsilon f(b) \end{pmatrix} \right) \right) \cdot \dot{u}'(b) \right) \right. \\ \left. - \frac{\epsilon}{2} \left( p^b \ \epsilon f(b) \right) \cdot \bar{\mathbf{c}}^b \cdot \left( \begin{pmatrix} p^b \ \epsilon f(b) \end{pmatrix} \right) \right) + \mathcal{O}(\epsilon^3). \quad (7.7) $$
This expression of the pseudo-energy was announced in (1.1) and is one of our main results. In Section 8, it will be used to generate the equilibrium boundary-value problem.

The expressions of the symmetric 2×2 matrix $\bar{\mathbf{k}}^b$ (boundary elasticity) appearing in (7.7), of the vector $\bar{\mathbf{\ell}}^b \in \mathbb{R}^3$ (condensation of the effective traction, dual to $u$), of the 2×2 matrix $\bar{\mathbf{r}}^b$ (condensation for higher-order tractions, dual to $u'$ and $u''$) and of the 2×2 symmetric matrix $\bar{\mathbf{c}}^b$ (pure-loading energy arising from the elimination of the boundary-layer amplitudes) are given in Appendices D and F for the first endpoint, $b=0$. The expressions relevant to the second endpoint, $b=\ell$, are obtained in the companion notebook by applying the mirror symmetry $j \leftarrow N-j$ mentioned in Remark 6.1, which maps $b=\ell$ to $b=0$. This saves us from repeating the boundary layer analysis at $b=\ell$.

The values of the effective constants appearing in the pseudo-energy (7.7) are shown in Table 7.1. These values were obtained in the companion notebook [TA25] using the formulas presented in Appendices D and F. They are used for validation purpose in Section 10.

<table><caption>**Table 7.1.** Effective properties appearing in the pseudo-energy (7.7) for the particular choice of spring constants in Table 3.1. The star symbols denote entries that can be ignored as they produce negligible energy contributions, of order $\epsilon^3$. All values are obtained in symbolic form in the companion notebook [TA25]: the numerical value -1.1947 appearing in the top-left entry of $\bar{\mathbf{k}}^0$, for instance, is a numerical approximation of $-(77+\sqrt{345})/80$.</caption><thead><tr><th>bulk</th><th>boundary b=0</th><th>boundary b=ℓ</th></tr></thead><tbody><tr><td>$\bar{\mathbf{K}} = \bar{\mathbf{Q}} \begin{pmatrix} 1 & \frac{1}{l_0} & \frac{29}{240} \\ \text{sym} & \frac{1}{l_6} & * \\ \text{sym} & * & * \end{pmatrix}$</td><td>$\bar{\mathbf{k}}^0 = \bar{\mathbf{Q}} \begin{pmatrix} -1.1947 & -0.7855 \\ \text{sym} & * \end{pmatrix}$<br>$\bar{\mathbf{r}}^0 = \begin{pmatrix} 1 & \frac{1}{l_2} & -\frac{11}{96} \\ 1.2947 & 1.4786 & * \\ 0.9063 & * & * \end{pmatrix}$<br>$\bar{\mathbf{c}}^0 = \frac{1}{\bar{\mathbf{Q}}} \begin{pmatrix} 2.6630 & 2.1136 \\ \text{sym} & * \end{pmatrix}$</td><td>$\bar{\mathbf{k}}^\ell = \bar{\mathbf{Q}} \begin{pmatrix} -1.26968 & 0.5831 \\ \text{sym} & * \end{pmatrix}$<br>$\bar{\mathbf{r}}^\ell = \begin{pmatrix} 1 & \frac{1}{l_2} & +\frac{11}{96} \\ -1.1697 & -1.2137 & * \\ 0.7039 & * & * \end{pmatrix}$<br>$\bar{\mathbf{c}}^\ell = \frac{1}{\bar{\mathbf{Q}}} \begin{pmatrix} 1.8607 & 1.4698 \\ \text{sym} & * \end{pmatrix}$</td></tr></tbody></table>

## 7.4. Final relocalization formulas

The outcome of the two-scale analysis is the pseudo-energy (7.7), along with relocalization formulas, which we discuss now. They are obtained by inserting the solution (7.6) for the boundary-layer parameters into the expression (7.1) of the nodal displacement $\mathbf{W}_i$:

$$ \mathbf{W}_i = \mathcal{L}^{\text{bk}} \cdot \ddot{\mathbf{u}}(\epsilon i) + \epsilon \sum_{b \in (0, \ell)} [ (\mathcal{L}_i^b)^\top \cdot \bar{\mathbf{\ell}}^b ] \cdot ( \dot{\mathbf{u}}'(b) |p^b \epsilon f(b) ) + \mathcal{O}(\epsilon^3). \quad (7.8) $$

The sum runs over the two boundary layers decaying exponentially away from the endpoints. The right-hand side depends only on the macroscopic displacement $u(b)$ and its gradient $u'(b)$, in addition to the known parameters of the problem, namely the boundary load $p^b$, $f(b)$ and $f'(b)$. Once the macroscopic solution $u$ has been obtained (see the following sections), Equation (7.8) makes it possible to reconstruct the nodal displacements and thus delivers a solution to the original lattice problem.

A similar formula for the microscopic strain $\mathbf{E}_j$ can be obtained by inserting (7.6) into (7.2),

$$ \mathbf{E}_j = \mathcal{E}^{\text{bk}}(\mathcal{Y}_2) \cdot \ddot{\mathbf{u}}'(x^j) + \sum_{b \in (0, \ell)} [\mathcal{E}_j^b \cdot \bar{\mathbf{\ell}}^b] \cdot (\dot{\mathbf{u}}'(b) |p^b \epsilon f(b)) + \mathcal{O}(\epsilon^3, \epsilon^2 z^{-j}, \epsilon^2 z^{-(N-j)}). \quad (7.9) $$

The error terms $\mathcal{O}(\epsilon^2 z^{-i})$ and $\mathcal{O}(\epsilon^2 z^{-(N-i)})$ are produced when the $\epsilon^2$ error in the boundary-layer parameters $\mathbf{v}^b$, see (7.6), is multiplied by the boundary-layer strain $\mathcal{E}_j^b$ in (7.2) which, by Equations (6.15) and (6.13), can be estimated as $\mathcal{E}_j^b = \mathcal{O}(\mathcal{L}_j^b) = \mathcal{O}(z^{-i}, z^{-(N-i)})$. Overall, the error term in (7.9) is of order $\epsilon^3$, except in layers of size $\epsilon$ near the endpoints, where it is of order $\epsilon^2$, and this subtlety would be missed if the right-hand side of Equation (7.8) were inserted into the definition (3.6) of the microscopic strain $\mathbf{E}_j$, as we have already noted in Remark 7.1.

**Remark 7.2.** The first term $\mathcal{E}^{\text{bk}}(\mathcal{Y}_2) \cdot \ddot{\mathbf{u}}'(x^j)$ in (7.9) makes the relocalized strain $\mathbf{E}_j$ depend on the unknown homogenization constant $\mathcal{Y}_2$ as $\epsilon^2 \mathcal{Y}_2( +2 -2 \ 0 \ 0 +2 ) u''(x^j) = \mathcal{O}(\epsilon^2)$, see Equation (5.13). In Section 10.2, we show that this dependence on $\mathcal{Y}_2$ is not significant. We could remove this dependence entirely by relocalizing $\mathbf{E}_j$ with a coarser error, $\mathcal{O}(\epsilon^2)$, than in (7.9), but prefer not to do so as the parametric analysis in terms of $\mathcal{Y}_2$ done in Section 10.2 offers interesting insights.

## 8. EQUILIBRIUM PROBLEM

In this section, we generate the equilibrium problem for the macroscopic displacement $u$ variationally from the pseudo-energy (7.7). Effective boundary conditions are included.

The equilibrium problem is obtained perturbatively, by formally making the pseudo-energy (7.7) stationary and seeking a solution $u$ in the form of a series expansion: we solve the following variational equality order by order in $\epsilon$ for the successive $u_{(i)}$'s,

$$ \forall \delta u, \quad D\Phi_{\epsilon}[u_{(0)} + \epsilon u_{(1)} + \epsilon^2 u_{(2)} + \dots; \delta u] = 0. \quad (8.1) $$

We claim that this approach produces the same series solution as our earlier, more traditional approach [TAL25b], where higher-order homogenization was carried out at the level of equilibrium and was combined with boundary-layer solutions using matched-asymptotic expansions.
## 8.1. Raw equilibrium equations in summed form

We start by putting back the expansion $u_{(0)}+\epsilon u_{(1)}+\epsilon^2 u_{(2)}+\cdots$ in its box $u(x)$, see (4.3), focussing on the equilibrium problem associated with the stationarity condition

$$ \forall \delta u, \quad D\Phi_{\epsilon}[u; \delta u] = 0, \qquad (8.2) $$

It is understood that the equilibrium problem has to be solved order by order, see Section 8.5.

In Appendix G, the first variation $D\Phi_{\epsilon}[u;\delta u]$ is worked out, including the standard integrations by parts. By the fundamental lemma of the calculus of variations, we set to zero all the coefficients of the test functions appearing in the final expression (G.12) of the first variation. This delivers the equilibrium problem in the form

$$ \begin{align} \forall x \in (0, \ell), \quad n'(x) + f(x) &= \mathcal{O}(\epsilon^3) \notag \\ \forall b \in \{0, \ell\}, \quad s^b n(b) &= \bar{\mathbf{r}}^b \cdot (p^b \epsilon \dot{\mathbf{f}}(b)) + \mathcal{O}(\epsilon^3) \tag{8.3} \\ \forall b \in \{0, \ell\}, \quad \begin{pmatrix} \epsilon \\ \epsilon^2 \end{pmatrix} \bar{\mathbf{r}}^b \cdot (p^b \epsilon f(b)) &= \begin{pmatrix} \epsilon \\ \epsilon^2 \end{pmatrix} \cdot (\bar{\mathbf{k}}^b + s^b \mathbf{l}(\bar{\mathbf{K}})) \cdot \dot{\mathbf{u}}'(b) + \mathcal{O}(\epsilon^3). \notag \end{align} $$

These equations can be interpreted respectively as the equilibrium in the bulk, the equilibrium of the boundary, and two higher-order boundary conditions (the last equation being an equality between vectors of dimension 2). The 2×2 matrix $\mathbf{l}(\bar{\mathbf{K}})$ appearing in (8.3)$_3$ is not necessarily symmetric, and is defined in terms of the homogenized elasticity matrix $\bar{\mathbf{K}}$ in (G.8). The quantity $s^b = \pm 1$ is the unit normal at the boundary,

$$ s^0 = -1, \quad s^\ell = +1. \qquad (8.4) $$

The quantity $n(x)$ appearing in the equilibrium (8.3) is the internal force. It is given by the effective constitutive law, see (G.3) and (G.11), as

$$ n(x) = \bar{\mathbf{K}}_{11} u'(x) - \epsilon^2 (\bar{\mathbf{K}}_{22} - 2\bar{\mathbf{K}}_{13}) u'''(x) + \mathcal{O}(\epsilon^3). \qquad (8.5) $$

With $\bar{\mathbf{K}}_{11}$ as the homogenized traction modulus, the leading-order term $\bar{\mathbf{K}}_{11}u'(x)$ yields the internal force in a classical bar. The gradient effect is captured by the second term in the right-hand side of (8.5), proportional to the second gradient of strain. Even though the constitutive law (8.5) depends on the even derivatives $u'$ and $u'''$ of the strain only, the asymmetry of the lattice, say for $Q_{III} \neq Q_{IV}$, is correctly taken into account via the boundary conditions, via the fact that the higher-order modulus $\bar{\mathbf{K}}_{22}-2\bar{\mathbf{K}}_{13}$ is different for the mirror-symmetric lattice, and via the relocalization formula 7.4.

In the pseudo-energy (7.7), the symmetric matrix $\bar{\mathbf{c}}^b$ was coupled to load parameters but not to the displacement $u(x)$ and thus does not appear in the equilibrium (8.3). It only plays a role in the formulation of the fixed-endpoint conditions, see Section 9.

In the higher-order boundary condition (8.3)$_3$, the diagonal matrix diag $(\epsilon, \epsilon^2)$ found on both sides has been produced by the test function $\epsilon \delta \dot{\mathbf{u}}'(b) = \text{diag}(\epsilon \delta u'(b), \epsilon^2 \delta u''(b))$. These diagonal matrices are eliminated later in (8.8), which complicates the expression of the remainder.

Before proceeding to solve the equilibrium problem (8.3) order by order in $\epsilon$, we discuss the special form of the higher-order tractions described by the matrix $\bar{\mathbf{r}}^b$, appearing in the pseudo-energy in (7.7) and in the higher-order boundary conditions (8.3)$_3$.

## 8.2. Solvability condition

The higher-order boundary condition (8.3)$_3$ depends on the unknown $u$ through the values of the strain and strain gradient, $\dot{\mathbf{u}}'(b) = (u'(b) \epsilon u''(b))$. The latter can be found using the first two equilibrium equations in (8.3)$_{1,2}$ and the constitutive law (8.5). By inserting the equilibrium value of $\dot{\mathbf{u}}'(b)$ into (8.3)$_3$, we can derive a necessary condition for equilibrium, from which the unknown $u$ has been eliminated. We proceed to derive this *solvability condition*.

When truncated to order $\epsilon$, the constitutive law (8.5) yields $n(x) = \bar{\mathbf{K}}_{11}u'(x) + \mathcal{O}(\epsilon^2)$. Combining it with the equilibrium in the bulk (8.3)$_1$ evaluated at the endpoint $x=b$ on the first hand, and with equilibrium of the boundary (8.3)$_2$ on the other hand, we get

$$ \begin{align} \dot{\mathbf{u}}'(b) &= \frac{\dot{n}(b)}{\bar{\mathbf{K}}_{11}} + \mathcal{O}(\epsilon^2) \notag \\ &= \frac{1}{\bar{\mathbf{K}}_{11}} \left( \frac{n(b)}{\epsilon n'(b)} \right) + \mathcal{O}(\epsilon^2) \tag{8.6} \\ &= \frac{1}{\bar{\mathbf{K}}_{11}} \left( \frac{1/s^b \bar{\mathbf{t}}^b \cdot (p^b \epsilon f(b))}{-\epsilon f(b)} \right) + \mathcal{O}(\epsilon^2). \notag \end{align} $$

Having downgraded the error to $\mathcal{O}(\epsilon^2)$, we can use $(p^b \epsilon \dot{\mathbf{f}}(b)) = (p^b \epsilon f(b) \epsilon^2 f'(b)) + \mathcal{O}(\epsilon^2) = \mathcal{P}^\top \cdot (p^b \epsilon f(b)) + \mathcal{O}(\epsilon^2)$, where $\mathcal{P}^\top$ appends a 0 to a vector of dimension 2, see (2.3). Using block-matrix notation, this yields

$$ \dot{\mathbf{u}}'(b) = \frac{1}{\bar{\mathbf{K}}_{11}} \begin{pmatrix} 1/s^b \bar{\mathbf{t}}^b \cdot \mathcal{P}^\top \\ -(0 \ 1) \end{pmatrix} \cdot (p^b \epsilon f(b)) + \mathcal{O}(\epsilon^2). \qquad (8.7) $$

Besides, we remove the matrix present on both sides of the higher-order boundary condition (8.3)$_3$ using left-multiplication by diag $(\epsilon^{-1} \epsilon^{-2})$,

$$ \bar{\mathbf{r}}^b \cdot (p^b \epsilon f(b)) = (\bar{\mathbf{k}}^b + s^b \mathbf{l}(\bar{\mathbf{K}})) \cdot \dot{\mathbf{u}}'(b) + \begin{pmatrix} \mathcal{O}(\epsilon^2) \\ \mathcal{O}(\epsilon^1) \end{pmatrix}. \qquad (8.8) $$
The remainder shown in the right-hand side is produced by diag ( ε⁻¹ ε⁻² ) · ( O ( ε³ ) O ( ε³ ) ) = ( O ( ε² ) O ( ε¹ ) ). Inserting the boundary strain predicted by (8.7), we get

$$
\bar{\mathbf{r}}^b \cdot (p^b \ \epsilon f(b)) = \mathbf{r}(\bar{\mathbf{K}}, s^b, \bar{\mathbf{k}}^b, \bar{\mathbf{t}}^b) \cdot (p^b \ \epsilon f(b)) + \begin{pmatrix} \mathcal{O}(\epsilon^2) \\ \mathcal{O}(\epsilon^1) \end{pmatrix}, \quad (8.9)
$$

where $\mathbf{r}(\bar{\mathbf{K}}, s^b, \bar{\mathbf{k}}^b, \bar{\mathbf{t}}^b)$ is known in terms of the homogenized elasticity $\bar{\mathbf{K}}$ in the bulk, and of the effective elasticity $\bar{\mathbf{k}}^b$, traction $\bar{\mathbf{t}}^b$ and unit normal $s^b=\pm 1$ at the boundary $b$:

$$
\mathbf{r}(\bar{\mathbf{K}}, s^b, \bar{\mathbf{k}}^b, \bar{\mathbf{t}}^b) := \frac{\bar{\mathbf{k}}^b + s^b \mathbf{l}(\bar{\mathbf{K}})}{\bar{\mathbf{K}}_{11}} \cdot \begin{pmatrix} 1/\epsilon^b \bar{\mathbf{t}}^b \cdot \mathbf{P}^\top \\ -(0 \ 1) \end{pmatrix}. \quad (8.10)
$$

Equation (8.9) has been derived by combining the equilibrium equations in (8.3). In the course of the derivation, as in many places in this paper, we have used the scaling assumptions $u' = \mathcal{O}(\epsilon^0)$, $u'' = \mathcal{O}(\epsilon^0)$, etc. applicable to solutions $u$ that can be expanded in powers of $\epsilon$. Overall, Equation (8.9) is therefore a necessary condition for the existence of a series solution to the equilibrium (8.2). As it does not depend on the unknown $u$, we call it a solvability condition.

To recapitulate, no series solution $u(x)$ exists unless the boundary traction $\bar{\mathbf{r}}^b \cdot (p^b \ \epsilon f(b))$ is of the special form (8.9).

8.3. Higher-order traction set by solvability condition

In this section, we show that there indeed exists a series solution $u(x)$ to the equilibrium (8.2); this proves that the higher-order boundary traction are of the form (8.9), by the argument above. Our starting point is the ansatz (4.2), i.e., we postulate that the microscopic problem has a solution $W_i$ of the form (4.2). We use this microscopic solution to construct a macroscopic solution $u^*$ making the pseudo-energy $\Phi_\epsilon$ stationary.

The argument is as follows:

* Rewriting relocalization formula (7.8) symbolically as $W_i=W_i^l[u]$, the fact that $\Phi_\epsilon[u]$ has been constructed as an order-$\epsilon^2$ approximation of the discrete energy can be expressed as

$$
\Phi_{\epsilon}[u] = \Phi_{d}(W_{0}^{l}[u], \dots, W_{N}^{l}[u]) + \mathcal{O}(\epsilon^{3}), \quad (8.11)
$$

for any macroscopic field $u(x)$.

* If the form of the ansatz (4.2) is correct, the microscopic solution $W_i$ derives from a macroscopic field $u^*$ in the sense that $W_i=W_i^l[u^*]$.

* The microscopic solution $W_i=W_i^l[u^*]$ satisfies the variational form of the equilibrium (principle of virtual work), see (B.2),

$$
\forall \delta \mathbf{W}_i, \quad D\Phi_d(\mathbf{W}_0^l[u^*], \dots, \mathbf{W}_N^l[u^*]; \delta \mathbf{W}_0, \dots, \delta \mathbf{W}_N) = 0. \quad (8.12)
$$

* Using formal asymptotics, we conclude, by applying the chained derivation rule to (8.11), and by using the discrete principle of virtual work (8.12) with particular virtual displacements,

$$
\begin{align}
D\Phi_{\epsilon}[u^*, \delta u] &= D\Phi_{d}(\mathbf{W}_{0}^{\mathrm{r}}[u^*], \dots, \mathbf{W}_{N}^{\mathrm{r}}[u^*]; D\mathbf{W}_{0}^{\mathrm{r}}[u^*, \delta u], \dots, D\mathbf{W}_{N}^{\mathrm{r}}[u^*, \delta u]) + \mathcal{O}(\epsilon^3) \tag{8.13} \\
&= \mathcal{O}(\epsilon^3). \nonumber
\end{align}
$$

The macroscopic displacement $u^*$ is thus a stationary point of the pseudo-energy $\Phi_\epsilon$ in the perturbative sense of (8.1).

In the following we will thus assume that the necessary condition (8.9) is satisfied, for an arbitrary choice of the boundary load ( p^b \epsilon f(b) ). The matrix r^b capturing the effective higher-order traction must therefore match

$$
\bar{\mathbf{r}}^b = \mathbf{r}(\bar{\mathbf{K}}, s^b, \bar{\mathbf{k}}^b, \bar{\mathbf{t}}^b) + \begin{pmatrix} 0 & 0 \\ 0 & * \end{pmatrix}. \quad (8.14)
$$

The '*' symbol in the right-hand side points to the fact that the component of $\bar{\mathbf{r}}^b$ in the lower-right corner is not resolved in second-order homogenization, as it contributes to the energy at order $\mathcal{O}(\epsilon^2 u''(b) \times \epsilon f(b)) = \mathcal{O}(\epsilon^3)$. The notation is meant to exclude this particular component from the equality (8.14), which applies to three out of the four components of $\bar{\mathbf{r}}^b$ only.

Equation (8.14) reveals a dependence between the higher-order traction $\bar{\mathbf{r}}^b$ and other properties of the effec-
tive model, namely the bulk homogenized elasticity $\bar{\mathbf{K}}$, the boundary properties $s^b=\pm 1$ (unit normal), $\bar{\mathbf{k}}^b$ (boundary
elasticity) and $\bar{\mathbf{t}}^b$ (effective traction). If the equality (8.14) were not satisfied, the order-by-order solution to the
variational equilibrium problem (8.2) would not exist. This would typically indicate that the ansatz (4.2) is incorrect.

We tested the validity of Equation (8.14) in the companion notebook [TA25], by comparing the effective matrix
$\bar{\mathbf{r}}^b$ generated directly by energy continualization, see Equations (F.4)$_2$ and (F.6) in the appendix, to the prediction
$\mathbf{r}(\bar{\mathbf{K}}, s^b, \bar{\mathbf{k}}^b, \bar{\mathbf{t}}^b)$ based on (8.10). This confirms that the three components of $\bar{\mathbf{r}}^b$ that we can resolve at second order
are correctly predicted by (8.14). This test could not be completed using symbolic values of the spring constants
$(Q_I, \dots, Q_{IV})$, as this produces algebraic expressions that are too complex to be simplified in a reasonable amount
of computing time. The test passed, however, when we assigned random positive values to the spring constants
$(Q_I, \dots, Q_{IV})$. We used 100 sets of random values, and verified (8.14) within an error $\sim 10^{-7}$ consistent with numerical
accuracy. We also checked that (8.14) holds exactly with the particular values in Table 3.1.
## 8.4. Returning to the equilibrium problem

Equation (8.14) has been obtained by rewriting the higher-order boundary condition (8.3)$_3$ with the help of the other equilibrium equations in (8.3)$_{1,2}$. It is a condition on the higher-order tractions which is satisfied, as we just checked. We conclude that the higher-order boundary conditions (8.3)$_3$ are redundant with the other equilibrium equations (8.3)$_{1,2}$. This was not apparent when we first derived the equilibrium in (8.3). This redundancy is a direct consequence of the special form (8.14) of the higher-order tractions produced by the two-scale analysis.

Discarding the redundant higher-order boundary conditions, we can rewrite the equilibrium problem in (8.3) and (8.5) as

$$ \begin{array}{l@{\quad}l} \text{equilibrium} & \begin{array}{l} \forall x \in [0, \ell], \quad n'(x) + f(x) = \mathcal{O}(\epsilon^3) \\ \forall b \in \{0, \ell\}, \quad s^b n(b) = \tilde{t}^b \cdot (p^b \epsilon \hat{f}(b)) + \mathcal{O}(\epsilon^3) \end{array} \\ \text{constitutive law} & \forall x \in [0, \ell], \quad n(x) = \bar{K}_{11} u'(x) - \epsilon^2 (\bar{K}_{22} - 2\bar{K}_{13}) u'''(x) + \mathcal{O}(\epsilon^3). \end{array} \qquad (8.15) $$

The values of the coefficients appearing in this equilibrium problem are listed in Table 8.1.

$$ \begin{array}{l@{\quad}l@{\quad}l} \text{bulk} & \bar{K}_{11} = \bar{Q} & \text{boundary } b=0 \\ & \bar{K}_{22} - 2\bar{K}_{13} = -3/40\bar{Q} & \tilde{t}^0 = (1 \ 1/2 \ -11/96) \quad \tilde{t}^\ell = (1 \ 1/2 \ +11/96) \end{array} $$

<caption>**Table 8.1.** Homogenized constants entering in the equilibrium problem (8.15) for the particular choice of spring constants in Table 3.1. The values have been extracted from Table 7.1. Note that the apparent gradient modulus $\bar{K}_{22}-2\bar{K}_{13}$ is negative.</caption>

## 8.5. Order by order solution

In Equation (8.15), we have obtained the equilibrium problem produced by the variational equality (8.2). In this section, we recall the perturbative nature of the equilibrium and seek a solution order by order, as in (8.1).

Inserting the expansion $u=u_{(0)}+\epsilon u_{(1)}+\epsilon^2 u_{(2)}+\dots$ into (8.15) yields a cascade of problems at the successive orders $\epsilon^i$, for $i \in \{0,1,2\}$. They all take the form of a classical elastic bar problem,

$$ \begin{array}{l} \forall x \in [0, \ell], \quad \bar{K}_{11} u''_{(i)}(x) + \tilde{f}_{(i)}(x) = 0 \\ \forall b \in \{0, \ell\}, \quad s^b \bar{K}_{11} u'_{(i)}(x) = \tilde{t}^b_{(i)}, \end{array} \qquad (8.16) $$

where the effective distributed force $\tilde{f}_{(i)}(x)$ and boundary traction $\tilde{t}^b_{(i)}$ applicable at order $\epsilon^i$ are given in Table 8.2. They are subjected to the global load balance condition $\tilde{t}^0_{(i)} + \int_0^\ell \tilde{f}_{(i)}(x) dx + \tilde{t}^\ell_{(i)} = 0$ at any order $i$.

The effective loading at order $\epsilon^2$ depends on the solution $u_{(0)}$, from leading order: the gradient effect is taken into account as a correction of order $\epsilon^2$. In addition, there are corrections $\tilde{t}^0_{3f}(b)$ and $\tilde{t}^b_{3f'}(b)$ at orders $\epsilon^1$ and $\epsilon^2$ coming from the endpoints.

$$ \begin{array}{l@{\quad}c@{\quad}c} \text{order } i & 0 & 1 \ 2 \\ \tilde{f}_{(i)}(x) & f(x) & 0 \\ \tilde{t}^b_{(i)} & \tilde{t}^b_1 p^b & \tilde{t}^b_2 f(b) \end{array} \quad \begin{array}{l} -(\bar{K}_{22}-2\bar{K}_{13})u'''_{(0)}(x) \\ \tilde{t}^b_3 f'(b) + s^b(\bar{K}_{22}-2\bar{K}_{13})u'''_{(0)}(x) \end{array} $$

<caption>**Table 8.2.** Effective loading entering in the classical elastic-bar problems (8.16) at the successive orders $\epsilon^i$: distributed force $\tilde{f}_{(i)}(x)$ and boundary traction $\tilde{t}^b_{(i)}$. The value of the constants is available from Table 8.1.</caption>

Alternatively, the following method delivers the same solution $u=u_{(0)}+\epsilon u_{(1)}+\epsilon^2 u_{(2)}+\dots$ directly in summed form. First, we integrate the equilibrium in the bulk (8.15)$_1$, and set the constant of integration using the boundary condition (8.15)$_2$ at $b=0$,

$$ n(x) = -\tilde{t}^0 \cdot (p^0 \epsilon \hat{f}(0)) - \int_0^x f dx + \mathcal{O}(\epsilon^3). \qquad (8.17) $$

The other boundary condition (8.15)$_2$ at $b=\ell$ yields the global load balance condition, $n(\ell)=\bar{K}^\ell \cdot (p^\ell \epsilon \hat{f}(\ell)) + \mathcal{O}(\epsilon^3)$, which is a necessary condition for the existence of an equilibrium. Next, the differential operator appearing in the constitutive law (8.15)$_1$ is inverted as

$$ \left[ \bar{K}_{11} \left( 1 - \epsilon^2 \frac{\bar{K}_{22} - 2\bar{K}_{13}}{\bar{K}_{11}} \frac{\mathrm{d}^2}{\mathrm{d}x^2} \right) \right]^{-1} = \frac{1}{\bar{K}_{11}} \left( 1 + \epsilon^2 \frac{\bar{K}_{22} - 2\bar{K}_{13}}{\bar{K}_{11}} \frac{\mathrm{d}^2}{\mathrm{d}x^2} \right) + \mathcal{O}(\epsilon^3), \qquad (8.18) $$

which is always possible in the framework of formal series expansions. The inverse differential operator yields the strain $u'(x)$ in terms of the internal force $n(x)$, and an additional integration yields the displacement as

$$ u(x) = u_{\epsilon}(0) + \int_{0}^{x} \left[ \frac{1}{\bar{K}_{11}} \left( 1 + \epsilon^2 \frac{\bar{K}_{22} - 2\bar{K}_{13}}{\bar{K}_{11}} \frac{\mathrm{d}^2}{\mathrm{d}x'^2} \right) n(x') \right] \mathrm{d}x'. \qquad (8.19) $$

**Remark 8.1.** A single boundary condition per endpoint appears to be sufficient to solve the equilibrium perturbatively, see (8.15)$_2$. With two endpoints, this makes 2 boundary conditions overall, which is consistent with the order (2) of the differential problems (8.16)$_1$ that are obtained at the successive orders. The equilibrium problem (8.15) is thus well posed when interpreted in a perturbative sense.
**Remark 8.2.** The equilibrium problem (8.15), however, cannot be interpreted beyond the perturbative framework. Suppose for a moment that we can ignore the remainders. By eliminating the internal force using the constitutive equation, we can rewrite the equilibrium (8.15)$_1$ as a fourth-order differential equation, $\bar{K}_{11}u''(x)-\epsilon^2(\bar{K}_{22}-2\bar{K}_{13})u'''(x)+f(x)=0$. As we still have 2 boundary conditions overall, the problem is ill posed. We have shown in previous work [TAL24] that this can be fixed by including 1 additional boundary conditions per endpoint, obtained by duplicating the leading order of the bulk equilibrium, $\bar{K}_{11}u''(b)+f(b)=0$. Even though it is well posed, the resulting boundary-value problem has no clear variational structure. Identifying a well-posed and variationally consistent equilibrium problem is the topic of our follow-up paper [TAL25a].

## 9. EXTENSION TO FIXED ENDPOINTS

The endpoints of the lattice have been so far subjected to forces. Indeed, we have not enforced any kind of kinematic constraint on the nodal positions $W_i$, and have included a contribution representing the endpoint forces $p^0$ and $p^\ell$ in the potential of the external load,

$$ \Phi_e = \dots -p^0 W_0 - p^\ell W_N^+ \quad (9.1) $$

see (3.10) and (3.12)$_2$. In this section, we cover the case of fixed-endpoint conditions, of the form

$$ W_0 = w^0, \quad W_N^+ = w^\ell, \quad (9.2) $$

where $w^b$ are the imposed displacements at the terminal nodes.

### 9.1. Strategy

Two strategies are possible for extending the analysis to fixed endpoints:

* One can insert the relocalized expressions of $W_0$ and $W_N^+$ obtained in the previous sections into the fixed-endpoint conditions (9.2), and solve them for the reaction force $p^b$. When this solution $p^b(w^b)$ is inserted into the balance of boundary forces (8.15)$_2$ derived earlier, it becomes a fixed-endpoint condition.

* Alternatively, one can Legendre-transform the pseudo-energy $\Phi_e[u, p^b]$ with respect to $p^b$. This yields a complementary pseudo-energy $\Phi_e^*[u, w^b]$ in terms of the prescribed displacement $w^b$ of the endpoints. By making it stationary with respect to $u$ order by order, one obtains the equilibrium problem for fixed endpoints.

The two approaches yield the same result. We follow the first one here as it is simpler. The second approach has the advantage of revealing the variational structure of the fixed-endpoint problem. It will be presented in Appendix [D] of our follow-up paper [TAL25a].

### 9.2. Fixed-endpoint condition in terms of macroscopic displacement

In view of the dependence of $\Phi_e$ on the reaction force $p^b$ in (9.1), we can obtain the displacements of the terminal nodes as $W_0=-\partial\Phi_d/\partial p^0$ and $W_N^-=\partial\Phi_d/\partial p^\ell$. Combining with the continuous approximation $\Phi_e=\Phi_d+\mathcal{O}(\epsilon^3)$, we rewrite the fixed-endpoint conditions in (9.2) as

$$ -\frac{\partial \Phi_e}{\partial p^b}[u, p^b] = w^b + \mathcal{O}(\epsilon^3). \quad (9.3) $$

In the left-hand side, we have made explicit the dependence of $\Phi_e[u, p^b]$ on the reaction force $p^b$, which has been implicit so far.

The left-hand side of (9.3) can be evaluated using the expression (7.7) of the pseudo-energy,

$$ \begin{aligned} \frac{\partial \Phi_e}{\partial p^b} &= \bar{t}_1^b u(b) + \epsilon(\bar{r}^b \cdot (1 \ 0)) \cdot \dot{u}'(b) + \epsilon(\bar{c}^b \cdot (1 \ 0)) \cdot (p^b \ \epsilon f(b)) + \mathcal{O}(\epsilon^3) \\ &= \epsilon \frac{p^b}{\bar{m}^b} + \bar{d}^b \cdot (\ddot{u}(b) \ \epsilon^2 f(b)) + \mathcal{O}(\epsilon^3), \end{aligned} \quad (9.4) $$

where

$$ \bar{m}^b = \frac{1}{\bar{c}_{11}^b}, \quad \bar{d}^b = \begin{pmatrix} \bar{t}_1^b \\ \bar{r}^b \cdot (1 \ 0) \\ \bar{c}_{12}^b \end{pmatrix} + \mathcal{O}(\epsilon^3). \quad (9.5) $$

Inserting (9.4) into the fixed-endpoint condition (9.3) and solving for $p^b$, we get the reaction force in terms of the prescribed endpoint displacement as

$$ p_e^b(w^b) = \frac{\bar{m}^b}{\epsilon}(w^b - \bar{d}^b \cdot (\ddot{u}(b) \ \epsilon^2 f(b))) + \mathcal{O}(\epsilon^2). \quad (9.6) $$

With this, the boundary condition (8.15)$_2$ takes the form

$$ \epsilon s^b n(b) = \bar{t}_1^b \bar{m}^b (w^b - \bar{d}^b \cdot (\ddot{u}(b) \ \epsilon^2 f(b))) + \bar{t}_2^b \epsilon^2 f(b) + \mathcal{O}(\epsilon^3). $$

Eliminating the stress variable using the constitutive law $n(b)=\bar{K}_{11}u'(b)+\mathcal{O}(\epsilon^2)$ and the equilibrium equation $f(b)=-n'(b)=-\bar{K}_{11}u''(b)+\mathcal{O}(\epsilon^2)$ and using the expression of $\bar{d}^b$ in (9.5)$_2$, this can be rewritten as

$$ \epsilon s^b \bar{K}_{11} u'(b) = \bar{t}_1^b \bar{m}^b (w^b - (\bar{t}_1^b u(b) + \epsilon \bar{r}_{11}^b u'(b) + \epsilon^2 \bar{r}_{21}^b u''(b))) - \epsilon^2 (-\bar{t}_1^b \bar{m}^b \bar{c}_{12}^b + \bar{t}_2^b \bar{K}_{11} u''(b) + \mathcal{O}(\epsilon^3)). \quad (9.7) $$
Collecting the terms that depend on $u$ and its derivative on one side, we obtain the final form of the fixed-endpoint condition as

$$ \bar{\boldsymbol{a}}^b \cdot \ddot{\boldsymbol{u}}(b) = w^b + \mathcal{O}(\varepsilon^3), \quad (9.8) $$

where $\bar{\boldsymbol{a}}^b$ is a known vector of dimension 3,

$$ \bar{\boldsymbol{a}}^b = \begin{pmatrix} \bar{t}_1^b \\ \bar{r}_{11}^b \\ \bar{r}_{21}^b \end{pmatrix} + \frac{\bar{K}_{11}}{\bar{m}^b \bar{t}_1^b} \begin{pmatrix} 0 \\ s^b \\ \bar{t}_2^b - \bar{c}_{12}^b \bar{t}_1^b \bar{m}^b \end{pmatrix}. \quad (9.9) $$

## 9.3. Equilibrium problem

The equilibrium problem with fixed endpoints is obtained by inserting in (8.15) the special form (9.8) of the boundary condition (8.15)$_2$,

$$ \begin{array}{l|l} \text{equilibrium} & \forall x \in [0, \ell], \quad n'(x) + f(x) = \mathcal{O}(\varepsilon^3) \\ \text{fixed endpoints} & \forall b \in \{0, \ell\}, \quad \bar{a}_1^b u(b) + \varepsilon \bar{a}_2^b u'(b) + \varepsilon^2 \bar{a}_3^b u''(b) = w^b + \mathcal{O}(\varepsilon^3) \\ \text{constitutive law} & \forall x \in [0, \ell], \quad n(x) = \bar{K}_{11} u'(x) - \varepsilon^2 (\bar{K}_{22} - 2\bar{K}_{13}) u'''(x) + \mathcal{O}(\varepsilon^3). \end{array} \quad (9.10) $$

The effective constants appearing in the equilibrium problem are listed in Table 9.1. With $\bar{a}_1^b=1$, the fixed-endpoint condition (9.10)$_2$ appears to be a refined version of the condition $u(b)=w^b+\mathcal{O}(\varepsilon)$ applicable at leading order for classical bars. By contrast with what happened in the case of force control at endpoints, the coefficients $\bar{c}_{ij}^b$ of the pure-loading term in the pseudo-energy (7.7)—and in particular $\bar{c}_{11}^b=1/\bar{m}^b—$, enter into the boundary problem, see (9.5)$_1$ and (9.9).

The equilibrium problem (9.10) must be solved order by order, as we did earlier in Section 8.5.

$$ \begin{array}{c} \text{bulk} \\ \bar{K}_{11} = \bar{Q} \\ \bar{K}_{22} - 2\bar{K}_{13} = -3/40\bar{Q} \end{array} \quad \begin{array}{c} \text{boundary } b=0 \\ \bar{a}^0 = \begin{pmatrix} 1 \\ -1.3683 \\ 0.1242 \end{pmatrix} \\ \bar{a}^\ell = \begin{pmatrix} 1 \\ 0.6911 \\ 0.1645 \end{pmatrix} \end{array} \quad \begin{array}{c} \text{boundary } b=\ell \\ \bar{K}_{11} = \bar{Q} \\ \bar{K}_{22} - 2\bar{K}_{13} = -3/40\bar{Q} \end{array} $$

<caption>**Table 9.1.** Effective constants entering in the equilibrium problem corresponding to fixed endpoints (9.10) for the particular choice of spring constants in Table 3.1.</caption>

# 10. VALIDATION

In this section, we verify the equilibrium problems produced by the energy method in Section 8 and 9. In Section 10.1, we compare them to the equilibrium problem obtained, for the same lattice, by the well-established method of matched asymptotic expansions in our previous work [TAL25b]. In Section 10.2, we propose a self-contained validation, in a geometry where the lattice is pinned at the first node and subjected to centrifugal forces. A quantitative analysis of convergence is included.

## 10.1. Equivalence with matched asymptotic expansions

In previous work [TAL25b], we have analyzed the same discrete lattice as in the present work, with both endpoints fixed. We applied the technique of matched asymptotic expansions to the discrete equations of equilibrium, and have obtained an equilibrium problem similar to (9.10), see Equation [3.7] in [TAL25b]. The loading parameters $(\tilde{f}, v_0, v_L^*)$ and effective constants $(K, B, \ell_0, \tilde{\ell}_0, \ell_L, \tilde{\ell}_L)$ used in the earlier paper correspond respectively to $(f, w^0, w^\ell)$ and to $(\bar{K}_{11}, \bar{K}_{22}-2\bar{K}_{13}, \bar{a}_2^0, \bar{a}_3^0, \bar{a}_2^\ell, \bar{a}_3^\ell)$ in the present paper.

Specifically, the boundary-value problems [3.7] from [TAL25b] and (9.10) in the present work are made up of the same equilibrium equations $\bar{K}_{11}u''(x)-\varepsilon^2(\bar{K}_{22}-2\bar{K}_{13})u'''(x)+f(x)=\mathcal{O}(\varepsilon^3)$ and of the same boundary conditions $\bar{a}_1^b u(b)+\varepsilon\bar{a}_2^b u'(b)+\varepsilon^2\bar{a}_3^b u''(b)=w^b+\mathcal{O}(\varepsilon^3)$ for $b \in \{0, \ell\}$. All the homogenized constants entering in these equations have identical values, up to numerical errors—in particular, the value of $B/K=-(9/32)/(15/4)=-3/40$ from Equation [3.2] in the previous paper matches the value of $(\bar{K}_{22}-2\bar{K}_{13})/\bar{Q}$ from Table 9.1, and the constants characterizing the effective fixed-boundary condition, $\ell_0=-1.37$, $\tilde{\ell}_0=0.124$, $\ell_L=0.691$, $\tilde{\ell}_L=0.165$ in [5.11] and [6.11], agree with the values of $\bar{a}^b$ in Table 9.1.

The only difference is that an additional boundary-condition was included in the equilibrium problem [3.7] of our previous work, which can be written as $\varepsilon^2(\bar{K}_{11}u''(b)+f(b))=\mathcal{O}(\varepsilon^3)$ in our current notation. The main claim of this previous paper was that this additional boundary condition makes it possible to solve the equilibrium problem as a conventional boundary-value problem, and not just in the perturbative sense. We checked in Section [7] of the previous paper that this additional boundary condition plays no role when a solution is sought order by order, as we do in the present work.

The two versions of the equilibrium problems are therefore strictly equivalent. In Section [8] of the previous paper, the equilibrium problem has been validated by comparing its solution to the discrete solution of the lattice problem: a convergence test confirmed that the effective equilibrium problem is accurate to order $\varepsilon^2$ included. This convergence test was based on the error between the discrete solution, and the relocalized solution predicted by the effective model.
This confirms that the stationarity principle proposed in the present work (i) delivers the same result as when the method of matched asymptotic expansions is applied to the discrete equilibrium equations, and (ii) that this result is asymptotically correct to order $\epsilon^2$ included.

## 10.2. Self-contained validation: rotating chain subjected to centrifugal forces

<figure><img src="image_4.png" /><figcaption>**Figure 10.1.** A one-dimensional lattice (orange) spins at angular velocity $\omega$ and stretches under the effect of the centrifugal force $f(x)$ applied at the nodes. The first endpoint $0^-$ is pinned at its initial position $x_0=-1/(4N)$.</figcaption></figure>

**Discrete model** In order to test the prediction (10.8), we set up discrete lattices having a variable number of cells, $N \ge 1$. The cell size is

$$ \epsilon = \ell/N = \mathcal{O}(N^{-1}). \quad (10.1) $$

The position of the first vertex, $0^-$, is fixed to its reference position $x_0=-1/(4N)$, see (3.1). The centrifugal force is represented by a distributed force $f(x)=\rho\omega^2 x$, and we apply no point-like force at the end $b=\ell$, i.e., $p^\ell=0$. By (3.10), the discrete force on any node labelled $i^\pm$ with $0 \le i \le N$ and $(i, \pm) \ne (0, -)$, is $\frac{1}{2N}f(x_i^\pm)=\frac{1}{2N}\rho\omega^2 x_i^\pm$, where $x_i^\pm$ is its reference position in (3.1). This problem is of the kind discussed in Section 3, but with a fixed-endpoint condition $w^0=0$ at the endpoint $b=0$,

$$ f(x) = \rho\omega^2 x, \quad w^0 = 0, \quad p^\ell = 0. \quad (10.2) $$

**Solution to effective model** The effective model is obtained by combining the fixed-endpoint condition (9.10) at $x=0$ and the force-controlled boundary condition (8.15) with the particular expression (10.2) of the load,

$$ \begin{array}{ll} \text{equilibrium} & \forall x \in [0, \ell], \quad n'(x) + \rho\omega^2 x = \mathcal{O}(\epsilon^3) \\ \text{fixed endpoint} & b=0, \quad \bar{a}_1^0 u(0) + \epsilon \bar{a}_2^0 u'(0) + \epsilon^2 \bar{a}_3^0 u'''(0) = \mathcal{O}(\epsilon^3) \\ \text{free endpoint} & b=\ell, \quad (\bar{t}_1^\ell \quad \bar{t}_2^\ell \quad \bar{t}_3^\ell) \cdot (0 \quad \epsilon \rho \omega^2 \ell \quad \epsilon^2 \rho \omega^2) = n(\ell) + \mathcal{O}(\epsilon^3) \\ \text{constitutive law} & \forall x \in [0, \ell], \quad n(x) = \bar{K}_{11} u'(x) - \epsilon^2 (\bar{K}_{22} - 2\bar{K}_{13}) u'''(x) + \mathcal{O}(\epsilon^3). \end{array} \quad (10.3) $$

We solve this boundary-value problem order by order in powers $\epsilon$, adopting the high-level presentation proposed in the second half of Section 8.5.

The solution for the internal force is given by integrating (10.3)$_1$ with the initial condition (10.3)$_3$. The result is given by a similar expression as earlier in (8.17), but starting the integration from $b=\ell$ and not from $b=0$,

$$ \begin{aligned} n(x) &= \frac{\rho\omega^2}{2}(\ell^2 - x^2) + \epsilon\bar{t}_2^\ell \rho\omega^2 \ell + \epsilon^2\bar{t}_3^\ell \rho\omega^2 + \mathcal{O}(\epsilon^3) \\ &= \rho\omega^2 \ell^2 \left(\frac{1}{2}\left(1 - \frac{x^2}{\ell^2}\right) + \bar{t}_2^\ell \frac{\epsilon}{\ell} + \bar{t}_3^\ell \frac{\epsilon^2}{\ell^2}\right) + \mathcal{O}(\epsilon^3). \end{aligned} \quad (10.4) $$

Inverting the constitutive law, we find, as earlier in (8.19),

$$ \begin{aligned} u'(x) &= \frac{1}{\bar{K}_{11}} \left( 1 + \epsilon^2 \frac{\bar{K}_{22} - 2\bar{K}_{13}}{\bar{K}_{11}} \frac{d^2}{dx^2} \right) n(x) + \mathcal{O}(\epsilon^3) \\ &= \frac{\rho\omega^2 \ell^2}{\bar{K}_{11}} \left( \frac{1}{2} \left( 1 - \frac{x^2}{\ell^2} \right) + \bar{t}_2^\ell \frac{\epsilon}{\ell} + \left( \frac{\bar{K}_{22} - 2\bar{K}_{13}}{\bar{K}_{11}} \right) \frac{\epsilon^2}{\ell^2} \right) + \mathcal{O}(\epsilon^3) \end{aligned} \quad (10.5) $$

Next, we use the fixed-endpoint condition at $x=0$ in (10.3)$_2$ to obtain the constant of integration

$$ \begin{aligned} u(0) &= -\frac{\epsilon}{\bar{a}_1^0} (\bar{a}_2^0 u'(0) + \epsilon \bar{a}_3^0 u''(0)) + \mathcal{O}(\epsilon^3) \\ &= \frac{\rho\omega^2 \ell^3 \bar{a}_2^0}{\bar{K}_{11} \bar{a}_1^0} \left( -\frac{1}{2\ell} \bar{t}_2^\ell \frac{\epsilon}{\ell^2} \right) + \mathcal{O}(\epsilon^3). \end{aligned} \quad (10.6) $$

The displacement is then given by integration as

$$ u(x) = \frac{\rho\omega^2\ell^3}{\bar{K}_{11}} \left( \frac{1}{2\ell} \frac{x}{6\ell^3} + \frac{\epsilon}{\ell} \left( -\frac{\bar{a}_2^0}{2\bar{a}_1^0} + \bar{t}_2^\ell \frac{x}{\ell} \right) + \frac{\epsilon^2}{\ell^2} \left( -\bar{t}_2^\ell \frac{\bar{a}_2^0}{\bar{a}_1^0} + \left( \bar{t}_3^\ell - \frac{\bar{K}_{22} - 2\bar{K}_{13}}{\bar{K}_{11}} \right) \frac{x}{\ell} \right) \right) + \mathcal{O}(\epsilon^3). \quad (10.7) $$
For validation purposes, we use the homogenized constants corresponding to the particular spring values in Table 3.1, which yields

$$u(x) = \frac{\rho \omega^2 \ell^3}{\bar{K}_{11}} \left( \frac{x/\ell}{2} \left( 1 + \varepsilon + \frac{91}{240} \varepsilon^2 \right) - \frac{(x/\ell)^3}{6} + \underbrace{\left( \frac{37}{112} + \frac{2}{7} \sqrt{\frac{23}{15}} \right)}_{0.6842} (\varepsilon + \varepsilon^2) \right) + \mathcal{O}(\varepsilon^3). \quad (10.8)$$

**Comparison procedure** To assess the accuracy of the effective model, we compare, for a given lattice size $N \ge 1$:

* the discrete solution ($W_i^d$)$_{0 \le i \le N}$ (displacement) and ($E_j^d$)$_{0 \le j \le N-1}$ (strain) of the lattice problem,

* with the prediction of the effective model obtained by relocating the macroscopic displacement $u(x)$ given in (10.8), ignoring the error term $\mathcal{O}(\varepsilon^3)$. Specifically, the displacement ($W_i^e$)$_{0 \le i \le N}$ is generated using the relocalization formula (7.8) and the microscopic strain ($E_j^e$)(Y$_2$)$_{0 \le i \le N-1}$ using (7.9). In both relocalization formulas, the error terms $\mathcal{O}(\dots)$ are ignored. The dependence of the microscopic strain $E_j^e$ (Y$_2$) on the higher-order homogenization constant Y$_2$ is retained and studied systematically in Figure 10.3.

In both the discrete problem and in the effective model, we work in dimensionless units by setting $\rho\omega^2=1$, $\ell=1$, and by scaling the spring constants in such a way that the leading-order equivalent traction modulus is $\bar{K}_{11}=1$.

To compare the discrete and effective models, we use an error $e_W$ based on the displacement, and an error $e_\Phi$ based on the energy,

$$e_W(N) = \left( \frac{1}{N+1} \sum_{i=0}^{N} \| \mathbf{W}_i^e - \mathbf{W}_i^d \|^2 \right)^{1/2} \quad (10.9) \\ e_\Phi(N, \mathcal{Y}_2) = \| \mathbf{E}^e(\mathcal{Y}_2) - \mathbf{E}^d \|_{\mathcal{Q}}^2.$$

Here, $\|\cdot\|_{\mathcal{Q}}$ denotes the norm constructed from the strain energy $\Phi_s$ of the discrete lattice, see (3.12),

$$\|\mathbf{E}\|_{\mathcal{Q}} = \left( \varepsilon \sum_{j=0}^{N-1} \mathbf{E}_j^-(\mathcal{Q}^{\text{bk}} + \delta_j \mathcal{Q}^0 + \delta_{N-1-j} \mathcal{Q}^\ell) \cdot \mathbf{E}_j \right)^{1/2} . \quad (10.10)$$

**Remark 10.1.** Up to a factor 1/2, the error $e_\Phi(N, \mathcal{Y}_2)$ can be interpreted as the excess of total potential energy in the lattice $\Phi_d(W_i^e) - \Phi_d(W_i^d)$ introduced by the effective model. This follows from a classical identity in linear elasticity, giving the increase in total potential energy $\Phi_d$, in (3.11), with respect to its (minimum) value $\Phi_d(W_i^d)$,

$$\begin{align}
\Phi_d(W_i^e) &= \Phi_d(W_i^d + (W_i^e - W_i^d)) \nonumber \\
&= \Phi_d(W_i^d) + \frac{\partial \Phi_d}{\partial W_i}(W_i^e - W_i^d) + \frac{1}{2} \frac{\partial^2 \Phi_d}{\partial W_i^2}(W_i^d) : ((W_i^e - W_i^d) \otimes (W_i^e - W_i^d)) \nonumber \\
&= \Phi_d(W_i^d) + \frac{1}{2} \|\mathbf{E}^e(\mathcal{Y}_2) - \mathbf{E}^d\|_{\mathcal{Q}}^2 \tag{10.11} \\
&= \Phi_d(W_i^d) + \frac{1}{2} e_\Phi(N, \mathcal{Y}_2). \nonumber
\end{align}$$

This classical identity is established as follows. In the second line, the linear term in the Taylor expansion is zero because $W_i^d$ is a minimum of $\Phi_d$, and the expansion to order 2 in ($W_i^e - W_i^d$) is exact because $\Phi_d$ has quadratic dependence on $W_i$ in linear elasticity. In the third line, we have identified the second gradient of total potential energy as the square strain-energy norm in (10.10).

**Comparison results** The error $e_W(N)$ is plotted as a function of the number of cells $N$ in the lattice Figure 10.2(a). The convergence rate $e_W(N) \sim N^{-3} \sim \varepsilon^3$ in the displacement is consistent with the error $\mathcal{O}(\varepsilon^3)$ appearing in the relocalization formula (7.8). The second-order effective model not only features a fast convergence for $N \to \infty$, it also delivers highly accurate approximations for moderate values of $N$. This is visible in Figure 10.2(b), where a highly accurate prediction is obtained with $N$ as small as $N=4$.

<figure><img src="image_5.png" /><figcaption>**Figure 10.2.** Convergence of displacement for the spinning test. (a) Log-log plot of the error $e_W$ on the discrete displacement for different lattice sizes $N$. (b) Comparison of the nodal displacements for $N=4$. This particular solution corresponds to the green dot in (a). The error of the effective model corresponds in (b) to the slight misalignment in the green versus orange disks.</figcaption></figure>
A more stringent test of convergence is given in Figure 10.3 based on the energy error $e_{\Phi}(N, \mathcal{Y}_2)$. For generic values of the coefficient $\mathcal{Y}_2$, the excess energy converges as $e_{\Phi} \sim N^{-4} = \epsilon^4$. The black curve at the bottom of plots, corresponding to $\mathcal{Y}_2=1/48$, features an exceptionally high convergence rate, $e_{\Phi} \sim N^{-5} = \epsilon^5$ ($\mathcal{Y}_2=1/48$ is the correct value of $\mathcal{Y}_2$, as predicted by third-order homogenization). This shows that the convergence rate $e_{\Phi} = \mathcal{O}(\epsilon^4)$ observed in Figure 10.3 for arbitrary values of $\mathcal{Y}_2$ is the best that can be achieved in the context of second-order homogenization.

With $e_{\Phi} \sim \epsilon^4$ as the ultimate convergence rate in second-order homogenization, this parametric study confirms that there is no drawback to setting the parameter $\mathcal{Y}_2$ arbitrarily, say to $\mathcal{Y}_2=0$ or $\mathcal{Y}_2=1$.

<figure><img src="image_6.png" /><figcaption>**Figure 10.3.** Convergence analysis for the spinning test based on the energy norm. The different curves correspond to different values of the homogenization constant $\mathcal{Y}_2$ which is not available in the context of second-order homogenization. The error decreases as $\mathcal{O}(\epsilon^4)$ in general. In the exceptional case $\mathcal{Y}_2=1/48$ (lower curve) corresponding to the 'true' value of $\mathcal{Y}_2$ predicted by third-order homogenization, the error becomes $\mathcal{O}(\epsilon^5)$. The generic fourth-order convergence rate $e_{\Phi}=\mathcal{O}(\epsilon^4)$ is thus the best that can be achieved in the context of second-order homogenization.</figcaption></figure>

**Remark 10.2.** The norm $\|\mathbf{E}\|_{\mathcal{O}}$ introduced in (10.10) is equivalent to the discrete $L_2$-norm
$\|\mathbf{E}^e(\mathcal{Y}_2) - \mathbf{E}^d\|_2 := (\frac{1}{N}\sum_{j=0}^{N-1} \|\mathbf{E}_j^e(\mathcal{Y}_2) - \mathbf{E}_j^d\|^2)^{1/2}$, in the generic case where all the spring constants are non-zero and finite. The observed error $e_{\Phi} \sim \epsilon^4$ in Figure 10.3 thus corresponds to an error $\|\mathbf{E}^e(\mathcal{Y}_2) - \mathbf{E}^d\|_2 \sim \epsilon^{4/2} \sim \epsilon^2$ with the $L_2$-norm. The latter is caused by the unresolved dependence of $\mathbf{E}^e(\mathcal{Y}_2)$ on $\mathcal{Y}_2$ at order $\epsilon^2$, as noted in Remark 7.2 below Equation (7.9).

**Remark 10.3.** If we had made an error in the homogenized solution in the bulk at order $\epsilon^1$, this would have given rise to an error of order $\epsilon^1$ in the relocalization formula (7.9), over the entire domain, causing the residual in the $L_2$-norm to be larger, $\|\mathbf{E}^e(\mathcal{Y}_2) - \mathbf{E}^d\|_2 \sim \epsilon^1$. This would have produced a slow convergence $e_{\Phi} \sim \epsilon^{2 \times 1} \sim \epsilon^2$ in the energy. Likewise, if we had made an error in the order-$\epsilon$ correction to any of the boundary layers, we would have introduced an error of order $\epsilon^1$ in (7.9) in a region of size $\sim \epsilon$, hence larger residuals $\|\mathbf{E}^e(\mathcal{Y}_2) - \mathbf{E}^d\|_2 \sim \epsilon^{(1+2 \times 1)/2} \sim \epsilon^{3/2}$ and $e_{\Phi} \sim \epsilon^3$ than actually observed. The observed convergence rate $e_{\Phi} \sim \epsilon^4$ warrants that the analysis is correct to second order.

**Remark 10.4.** In Appendix H, we characterize more finely the maximum error that can be tolerated on the macroscopic displacement $u$ and its derivatives, in the sense that the convergence rates $e_W(N) = \mathcal{O}(\epsilon^3)$ and $e_{\Phi}(N, \mathcal{Y}_2) = \mathcal{O}(\epsilon^4)$ are preserved.

# 11. CONCLUSION

To sum up, we have analyzed a family of elastic lattices in dimension 1, which we used as a prototype of periodic elastic structures with boundaries. We have derived an effective equilibrium problem valid at order $\epsilon^2$ in the scale separation parameter $\epsilon \ll 1$. The boundary conditions effectively account for the boundary layers. They are hidden in the macroscopic formulation but can be restored explicitly using the localization formulas (Section 7.4). By making the pseudo-energy $\Phi_{\epsilon}[u]$ stationary order by order up to order $\epsilon^2$, we have recovered the equilibrium problem obtained in earlier work, up to an error consistent with second order analyses. We have implemented two-scale expansions at the energy level: the pseudo-energy has been obtained by inserting the two-scale expansion (4.2) for the nodal displacement into the energy of the discrete lattice, and not into the equilibrium equations.

By homogenizing at the energy level, we were hoping to preserve the energy structure of the original lattice problem. This goal has only been partially met yet. The pseudo-energy $\Phi_{\epsilon}[u]$ delivers the effective equilibrium problem, but in a restrictive, perturbative sense. It also features a negative higher-order modulus ($\bar{K}_{22}-2\bar{K}_{13}$), see Table 8.1. In spite of these limitations, the present work is a key step towards the identification of a simple variational principle underlying second-order homogenization, as we argue now.
In this paper, we have learnt that the effective higher-order tractions, dual to the gradients $u'(b)$ and $u''(b)$ of the displacement, can be expressed in terms of the effective elastic properties of the lattice and of the applied load as

$$ \vec{r}^b \cdot (p^b \ \epsilon f(b)) = \frac{\vec{k}^b + s^b \mathbf{l}}{\vec{k}_{11}} \cdot \left( \begin{matrix} 1/s_b \vec{l}_c^b \\ -\epsilon f(b) \end{matrix} \right) + \left( \begin{matrix} \mathcal{O}(\epsilon^2) \\ \mathcal{O}(\epsilon^1) \end{matrix} \right), \quad \text{where } \mathbf{l} = \left( \begin{matrix} \vec{k}_{12} & \vec{k}_{22} - \vec{k}_{13} \\ \vec{k}_{13} & * \end{matrix} \right). \qquad (11.1) $$

This property of the effective model has not been appreciated in earlier work.

Note that the expression (11.1) of the higher-order tractions follows by inserting (8.10) into (8.14). It uses the effective boundary traction $\vec{l}_c^b = \vec{t}^b \cdot \mathcal{P}^T \cdot (p^b \epsilon f(b)) + \mathcal{O}(\epsilon^2) = \vec{t}_1^b p^b + \epsilon \vec{t}_2^b f(b) + \mathcal{O}(\epsilon^2)$ and the matrix $\mathbf{l} = \mathcal{R} \cdot (\vec{K} - \mathcal{P} \cdot \vec{K} \cdot \mathcal{P}) \cdot \mathcal{P}^T$, see (G.8).

By the argument given in Section 8.3, the form (11.1) of the higher-order tractions is general, and we conjecture that it can be extended to any micro-structure, beyond the particular family of lattices studied here, and beyond dimension 1. The special form of the higher-order tractions make the variational derivation of the equilibrium quite peculiar, as noted in Section 8.4: out of the 3 boundary conditions per endpoint that one could expect for a second-order continuum, we found that 2 were automatically satisfied. In the context of variational formulations, it is quite unusual to have a single boundary condition surviving out of three, but this is precisely what makes the perturbative approach well-posed. We conjecture that this peculiar variational behavior is common to all effective energies produced by second-order homogenization.

In a follow-up paper [TAL25a], we will show that the special form (11.1) of the boundary tractions obtained in the present work makes it possible to truncate the pseudo-energy in a way which restores positivity. To show this, we will also use the positivity property established in Equation (D.5) in the Appendix,

$$ \vec{k}_{11}>0, \quad \vec{k}_{11}\vec{k}_{22}-(\vec{k}_{12})^2\geq0. \qquad (11.2) $$

In this follow-up work, the positive truncation of the pseudo-energy will ultimately lead us to a variational principle for second-order homogenization, which was our initial goal.

We close with a couple remarks concerning possible extensions of this work.

As noted in the introduction, the extension to a continuous periodic microstructure does not raise any difficulty: the main difference with the discrete case treated here is that the dependence of the microscopic displacement on the node type $\pm$ in the ansatz (4.2) is replaced with a dependence on a (continuous) fast variable. As an example of a continuous microstructure, we will investigate in another follow-up paper the case of a laminate that is invariant in the direction perpendicular to the fibers, hence effectively one-dimensional. The derivation is remarkably similar to that proposed in the present paper. The main difference is that the boundary layers can only be solved numerically, and the effective constants are therefore not available in closed analytical form.

The extension to (genuinely) two-dimensional or three-dimensional periodic continua is highly desirable. In this presentation of the one-dimensional case, we have refrained as much as possible from using specific tricks to dimension 1 so as to make this next step as easy as possible.

**Acknowledgements** This manuscript has been typeset using TeX\_MACS, a powerful, multi-platform and freely distributed scientific editor [vdHGG+13].

BIBLIOGRAPHY

[AA99] G. Allaire and M. Amar. Boundary layer tails in periodic homogenization. ESAIM: Control, Optimisation and Calculus of Variations, 4:209–243, 1999.

[AE09] G. Allaire and K. El Ganaoui. Homogenization of a conductive and radiative heat transfer problem. Varia, 7(3):1148–1170, 2009.

[AG05] M. Arndt and M. Griebel. Derivation of higher order gradient continuum models from atomistic models for crystalline solids. Multiscale Modeling & Simulation, 4(2):531–562, 2005.

[AL23] B. Audoly and C. Lestringant. An energy approach to asymptotic, higher-order, linear homogenization. Journal of Theoretical, Computational and Applied Mechanics, page 11414, 2023.

[Apo99] T. M. Apostol. An elementary view of Euler's summation formula. The American Mathematical Monthly, 106(5):409–418, 1999.

[BKS11] A. L. Bessoud, F. Krasucki, and M. Serpilli. Asymptotic Analysis of Shell-like Inclusions with High Rigidity. Journal of Elasticity, 103(2):153–172, 2011.

[BLP78] A. Bensoussan, J. L. Lions, and G. Papanicolaou. Asymptotic Analysis for Periodic Structures, volume 5. American Mathematical Soc., 1978.

[BLP79] A. Bensoussan, J. L. Lions, and G. C. Papanicolaou. Boundary Layers and Homogenization of Transport Processes. Publ. Res. Inst. Math. Sci., 15(1):53–157, 1979.

[Bou96] C. Boutin. Microstructural effects in elastic composites. International Journal of Solids and Structures, 33(7):1023–1051, 1996.

[Bou20] C. Boutin. Homogenization Methods and Generalized Continua in Linear Elasticity. In Holm Altenbach and Andreas Öchsner, editors, Encyclopedia of Continuum Mechanics, pages 1231–1231. Springer, Berlin, Heidelberg, 2020.

[BP89] N. Bakhvalov and G. Panasenko. Homogenisation: Averaging Processes in Periodic Media, volume 36 of Mathematics and Its Applications. Springer Netherlands, Dordrecht, 1989.

[BS11] C. Boutin and J. Soubestre. Generalized inner bending continua for linear fiber reinforced materials. International Journal of Solids and Structures, 48(3):517–534, 2011.

[BT08] A. Braides and L. Truskinovsky. Asymptotic expansions by $\Gamma$-convergence. Continuum Mechanics and Thermodynamics, 20(1):21–62, 2008.

[CT02] M. Charlotte and L. Truskinovsky. Linear elastic chain with a hyper-pre-stress. Journal of the Mechanics and Physics of Solids, 50(2):217–251, 2002.
[CZWK19] N. Challamel, H. Zhang, C.M. Wang, and J. Kaplunov. Scale effect and higher-order boundary conditions for generalized lattices, with direct and indirect interactions. *Mechanics Research Communications*, 97:1–7, 2019.

[Dum86] H. Dumontet. Study of a boundary layer problem in elastic composite materials. *ESAIM: Mathematical Modelling and Numerical Analysis*, 20(2):265–286, 1986.

[GK89] B. Gambin and E. Kröner. Higher-Order Terms in the Homogenized Stress-Strain Relation of Periodic Elastic Media. *Physica status solidi* (b), 151(2):513–519, 1989.

[GW65] D. C. Gazis and R. F. Wallis. Surface tension and surface modes in semi-infinite lattices. *Surface Science*, 3(1):19–32, 1965.

[KAL22] A. Kumar, B. Audoly, and C. Lestringant. Asymptotic derivation of a higher-order one-dimensional model for tape springs. 2022.

[LA20] C. Lestringant and B. Audoly. Asymptotically exact strain-gradient models for nonlinear slender elastic structures: A systematic derivation method. *Journal of the Mechanics and Physics of Solids*, 136:103730, 2020.

[Lio81] J. L. Lions. *Some Methods in the Mathematical Analysis of Systems and Their Control*. Science Press, Beijing, Gordon And Breach, New York, 1981.

[LM18] D. T. Le and J. J. Marigo. Second order homogenization of quasi-periodic structures. *Vietnam Journal of Mechanics*, 40(4):325–348, 2018.

[Min64] R. D. Mindlin. Micro-structure in linear elasticity. *Archive for Rational Mechanics and Analysis*, 16(1):51–78, 1964.

[Min65] R. D. Mindlin. Second gradient of strain and surface-tension in linear elasticity. *International Journal of Solids and Structures*, 1(4):417–438, 1965.

[MV97] S. Moskow and M. Vogelius. First-order corrections to the homogenised eigenvalues of a periodic composite medium. A convergence proof. *Proceedings of the Royal Society of Edinburgh: Section A Mathematics*, 127(6):1263–1299, 1997.

[NCH20] H. Nassar, H. Chen, and G. Huang. Microtwist elasticity: A continuum approach to zero modes and topological polarization in Kagome lattices. *Journal of the Mechanics and Physics of Solids*, 144:104107, 2020.

[Pan09] G. Panasenko. Boundary conditions for the high order homogenized equation: laminated rods, plates and composites. *Comptes Rendus Mécanique*, 337(1):8–14, 2009.

[S18] H. Abdoul-Anziz and P. Seppecher. Strain gradient and generalized continua obtained by homogenizing frame lattices. *Mathematics and Mechanics of Complex Systems*, 6(3):213–250, 2018.

[San80] E. Sanchez-Palencia. *Non-Homogeneous Media and Vibration Theory*. Number 127 in Lecture Notes in Physics. Springer Berlin Heidelberg, Berlin, Heidelberg, 1980.

[San86] E. Sanchez-Palencia. Homogenization in mechanics. A survey of solved and open problems. 44(1):1–45, 1986.

[SC00] V. P. Smyshlyaev and K. D. Cherednichenko. On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media. *Journal of the Mechanics and Physics of Solids*, 48:1325–1357, 2000.

[TA25] M. Thbaut and B. Audoly. Supplementary material for 'A stationarity principle generating effective boundary conditions for second-order homogenization'. https://archive.softwareheritage.org/ swh:1:dir:f3eae28e4c00e6fb88b57272de952a24719c7a37;origin=https://git.renater.fr/anonscm/git/manon-thbaut-pb/manon-thbaut-pb.git;visit=swh:1:snp:32055ee43d6942bcaee84bee4ffc1c5c17389cd1;anchor=swh:1:rev:dcc8739d3732a7bfe2025-AStationarity_principle generating_e effective b 2025.

[TAL24] M. Thbaut, B. Audoly, and C. Lestringant. Effective boundary conditions for second-order homogenization. *Journal of the Mechanics and Physics of Solids*, 190:105707, 2024.

[TAL25a] M. Thbaut, B. Audoly, and C. Lestringant. Asymptotic strain-gradient theory for one-dimensional continua. https://hal.science/hal-05129643, 2025. Preprint.

[TAL25b] M. Thbaut, B. Audoly, and C. Lestringant. Fixing non-positive energies in higher-order homogenization. *Journal of the Mechanics and Physics of Solids*, 2025.

[TB93] N. Triantafyllidis and S. Bardenhagen. On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models. *Journal of Elasticity*, 1993.

[TB96] N. Triantafyllidis and S. Bardenhagen. The influence of scale size on the stability of periodic solids and the role of associated higher order gradient continuum models. *Journal of the Mechanics and Physics of Solids*, 44(11):1891–1928, 1996.

[Tou62] R. A. Toupin. Elastic materials with couple-stresses. *Archive for Rational Mechanics and Analysis*, 11(1):385–414, 1962.

[vdHGG+13] J. van der Hoeven, A. Grozin, M. Gubinelli, G. Lecerf, F. Poulain, and D. Raux. GNU TEXmacs: a scientific editing platform. *ACM Communications in Computer Algebra*, 47(1–2):59–61, 2013.

[Vin16] V. Vinoles. *Problèmes d’interface en présence de métamatériaux: modélisation, analyse et simulations*. PhD thesis, Université Paris Saclay (COMUE), 2016.

[YAL24] Y. Ye, B. Audoly, and C. Lestringant. Asymptotic, second-order homogenization of linear elastic beam networks. *Journal of the Mechanics and Physics of Solids*, 188:105637, 2024.

## APPENDIX A. SUMMATION FORMULAS

Given an integer $M=\mathcal{O}(\varepsilon^{-1})$ and a smooth function $g$, the Euler–Maclaurin summation formula [Apo99] approximates the sum of $g$ over a regular grid $(\varepsilon\alpha, \varepsilon(\alpha+1), \dots, \varepsilon(\alpha+M))$ with step size $\varepsilon$ by an integral,

$$ \varepsilon \sum_{i=0}^{M} g(\varepsilon(i+\alpha)) = \int_{\varepsilon\alpha}^{\varepsilon(M+\alpha)} g \, dx + \left( \frac{\varepsilon}{2} g - \frac{\varepsilon^2}{12} g' \right)_{x=\varepsilon\alpha} + \left( \frac{\varepsilon}{2} g' + \frac{\varepsilon^2}{12} g' \right)_{x=\varepsilon(M+\alpha)} + \mathcal{O}(\varepsilon^3). \quad (A.1) $$

In this standard form, the bounds of the integral coincide with the first and last points of the grid, $\varepsilon\alpha$ and $\varepsilon(M+\alpha)$.

We extend this formula by shifting the bounds of the integral with respect to the terminal grid points. We introduce forcefully in the right-hand side of (A.1) an integral over the alternate domain $(x_{-}, x_{+}) := (0, \varepsilon(M+\alpha+\beta))$,

$$ \varepsilon \sum_{i=0}^{M} g(\varepsilon(i+\alpha)) = \int_{x_{-}}^{x_{+}} g \, dx + \left[ \left( \frac{\varepsilon}{2} g - \frac{\varepsilon^2}{12} g' \right)_{x=\varepsilon\alpha} - \int_{x_{-}}^{\varepsilon\alpha} g \, dx \right] + \left[ \left( \frac{\varepsilon}{2} g + \frac{\varepsilon^2}{12} g' \right)_{x=\varepsilon(M+\alpha)} - \int_{\varepsilon(M+\alpha)}^{x_{+}} g \, dx \right] + \mathcal{O}(\varepsilon^3). \quad (A.2) $$
Expanding the terms in each one of the square brackets, about $x_-$ and $x_+$ respectively, we get

$$
\begin{align}
\epsilon \sum_{i=0}^{M} g(\epsilon(i+\alpha)) &= \int_{x_{-}}^{x_{+}} g dx + \left[ \left( \frac{\epsilon}{2} g - \frac{\epsilon^2}{12} g' \right) + \epsilon \alpha \frac{\epsilon}{2} g' - \left( g \int_{x_{-}}^{x_{+}+\epsilon\alpha} dx' + g' \int_{x_{-}}^{x_{+}+\epsilon\alpha} (x'-x_{-}) dx' \right) \right]_{x=x_{-}} \nonumber \\
&\quad + \left[ \left( \frac{\epsilon}{2} g + \frac{\epsilon^2}{12} g' \right) - \epsilon \beta \frac{\epsilon}{2} g' - \left( g \int_{x_{+}-\epsilon\beta}^{x_{+}} dx' + g' \int_{x_{+}-\epsilon\beta}^{x_{+}} (x'-x_{+}) dx' \right) \right]_{x=x_{+}} + \mathcal{O}(\epsilon^3) \tag{A.3} \\
&= \int_{x_{-}}^{x_{+}} g dx + \left[ \left( \frac{1}{2} - \alpha \right) \epsilon g + \frac{1}{2} \left( -\frac{1}{6} + \alpha - \alpha^2 \right) \epsilon^2 g' \right]_{x=x_{-}} \nonumber \\
&\quad + \left[ \left( \frac{1}{2} - \beta \right) \epsilon g + \frac{1}{2} \left( +\frac{1}{6} - \beta + \beta^2 \right) \epsilon^2 g' \right]_{x=x_{+}} + \mathcal{O}(\epsilon^3) \nonumber
\end{align}
$$

The following particular cases are of interest:

* For $M=\infty$ and if $g$ decays for $x \to \infty$, we get the formula announced in (2.10),

$$
\epsilon \sum_{i=0}^{\infty} g(\epsilon(i+\alpha)) = \int_{0}^{\infty} g dx + \left(\frac{1}{2} - \alpha\right) \epsilon g(0) + \frac{1}{2} \left(-\frac{1}{6} + \alpha - \alpha^2\right) \epsilon^2 g'(0) + \mathcal{O}(\epsilon^3). \quad (A.4)
$$

* For $\alpha=\beta=1/2$ and $M=N-1$, the bounds of the integral are $(x_-, x_+) = (0, \epsilon N)$, and we get

$$
\epsilon \sum_{i=0}^{N-1} g\left(\epsilon\left(i+\frac{1}{2}\right)\right) = \int_{0}^{\epsilon N} g dx + \frac{\epsilon^2}{24} g'(0) - \frac{\epsilon^2}{24} g'(\epsilon N) + \mathcal{O}(\epsilon^3). \quad (A.5)
$$

* For $\alpha=\beta=0$ and $M=N$, we have $(x_-, x_+) = (0, \epsilon N)$ as well, and the approximation formula becomes

$$
\epsilon \sum_{i=0}^{N} g(\epsilon i) = \int_{0}^{\epsilon N} g dx + \frac{\epsilon}{2} g(0) - \frac{\epsilon^2}{12} g'(0) + \frac{\epsilon}{2} g(\epsilon N) + \frac{\epsilon^2}{12} g'(\epsilon N) + \mathcal{O}(\epsilon^3). \quad (A.6)
$$

Alternatively, this equation follows directly from the classical Euler-Maclaurin formula (A.1).

The last two formulas (A.5)-(A.6) are stated for the record and will be used in a follow-up paper.

APPENDIX B. DISCRETE EQUILIBRIUM EQUATIONS

In this Appendix, we derive the discrete equations of equilibrium of the lattice, by making its energy (3.11) sta-
tionary with respect to the nodal displacements.

Inserting the expression (3.6) of the strain $\boldsymbol{E}_j$ into the strain energy (3.12)$_1$, we get

$$
\Phi_s(\boldsymbol{E}_0, \boldsymbol{E}_1, \dots, \boldsymbol{E}_{N-1}) = \sum_{j=0}^{N-1} \frac{1}{2\epsilon} \left( \frac{\boldsymbol{W}_j}{\boldsymbol{W}_{j+1}} \right) : \mathcal{B}^{\top} : (\boldsymbol{Q}^{\text{bk}} + \delta_j \boldsymbol{Q}^0 + \delta_{N-1-j} \boldsymbol{Q}^{\ell}) : \mathcal{B} : \left( \frac{\boldsymbol{W}_j}{\boldsymbol{W}_{j+1}} \right). \quad (B.1)
$$

Adding the potential energy of the external work, we calculate the first variation of the total potential energy asso-
ciated with a perturbation ($\delta \mathbf{W}_0, \dots, \delta \mathbf{W}_N$) as

$$
\begin{align}
\mathrm{D}\Phi_d(\mathbf{W}_0, \dots, \mathbf{W}_N; \delta\mathbf{W}_0, \dots, \delta\mathbf{W}_N) &= \sum_{j=0}^{N-1} \frac{1}{\epsilon} \left( \frac{\delta\mathbf{W}_j}{\delta\mathbf{W}_{j+1}} \right) : \mathcal{B}^{\top} : (\mathbf{Q}^{\text{bk}} + \delta_j \mathbf{Q}^0 + \delta_{N-1-j} \mathbf{Q}^{\ell}) : \mathcal{B} : \left( \frac{\mathbf{W}_j}{\mathbf{W}_{j+1}} \right) - \sum_{i=0}^{N} \mathbf{F}_i \cdot \delta\mathbf{W}_i \nonumber \\
&= -\left(\frac{1}{\epsilon}\mathcal{D}^0: \left(\frac{\mathbf{W}_0}{\mathbf{W}_1}\right) + \mathbf{F}_0\right) \cdot \delta\mathbf{W}_0 - \sum_{i=1}^{N-1} \left(\frac{1}{\epsilon}\mathcal{D}^{\text{bk}}: \left(\frac{\mathbf{W}_{i-1}}{\mathbf{W}_i}\right) + \mathbf{F}_i\right) \cdot \delta\mathbf{W}_i \tag{B.2} \\
&\quad -\left(\frac{1}{\epsilon}\mathcal{D}^N: \left(\frac{\mathbf{W}_{N-1}}{\mathbf{W}_N}\right) + \mathbf{F}_N\right) \cdot \delta\mathbf{W}_N \nonumber
\end{align}
$$

where the tensors $\mathcal{D}^0, \mathcal{D}^{\text{bk}}, \mathcal{D}^N$ are given by

$$
\begin{align}
\mathcal{D}^0 &= - \left( \begin{array}{c} \mathcal{B}_{\leftarrow}^{\top} \cdot (\mathbf{Q}^{\text{bk}} + \mathbf{Q}^0) \cdot \mathcal{B}_{\leftarrow} \\ \mathcal{B}_{\leftarrow}^{\top} \cdot (\mathbf{Q}^{\text{bk}} + \mathbf{Q}^0) \cdot \mathcal{B}_{\rightarrow} \end{array} \right) \nonumber \\
\mathcal{D}^{\text{bk}} &= - \left( \begin{array}{c} \mathcal{B}_{\rightarrow}^{\top} \cdot \mathbf{Q}^{\text{bk}} \cdot \mathcal{B}_{\leftarrow} \\ \mathcal{B}_{\leftarrow}^{\top} \cdot \mathbf{Q}^{\text{bk}} \cdot \mathcal{B}_{\rightarrow} \\ \mathcal{B}_{\rightarrow}^{\top} \cdot (\mathbf{Q}^{\text{bk}} + \mathbf{Q}^{\ell}) \cdot \mathcal{B}_{\leftarrow} \\ \mathcal{B}_{\leftarrow}^{\top} \cdot (\mathbf{Q}^{\text{bk}} + \mathbf{Q}^{\ell}) \cdot \mathcal{B}_{\rightarrow} \end{array} \right) \tag{B.3} \\
\mathcal{D}^{\ell} &= - \left( \begin{array}{c} \mathcal{B}_{\rightarrow}^{\top} \cdot (\mathbf{Q}^{\text{bk}} + \mathbf{Q}^{\ell}) \cdot \mathcal{B}_{\leftarrow} \\ \mathcal{B}_{\leftarrow}^{\top} \cdot (\mathbf{Q}^{\text{bk}} + \mathbf{Q}^{\ell}) \cdot \mathcal{B}_{\rightarrow} \end{array} \right). \nonumber
\end{align}
$$

Here, $\mathcal{B}_{\leftarrow} = ((-1 \ 0 \ 0) \ (0 \ -1 \ 0) \ (-1 \ 0 \ 0))$ and $\mathcal{B}_{\rightarrow} = ((0 \ 0 \ 0) \ (+1 \ 0 \ 0) \ (+1 \ 0 \ 0))$ are 5x2
matrices extracted from the tensor $\mathcal{B}$ defined in (3.7), using $(\mathcal{B}_{\leftarrow})_{ij} = B_{i1j}$ and $(\mathcal{B}_{\rightarrow})_{ij} = B_{i2j}$ (these sub-tensors of
$\mathcal{B}$ deal with the adjacent cells on the left and on the right-hand side, respectively).
This yields the discrete equilibrium equations in the form

$$ 
\begin{aligned}
\frac{1}{\epsilon} \mathcal{D}^0: \begin{pmatrix} \mathbf{W}_0 \\ \mathbf{W}_1 \end{pmatrix} + \mathbf{F}_0 &= \mathbf{0} \quad (i=0) \\
\frac{1}{\epsilon} \mathcal{D}^{\text{bk}}: \begin{pmatrix} \mathbf{W}_{i-1} \\ \mathbf{W}_i \\ \mathbf{W}_{i+1} \end{pmatrix} + \mathbf{F}_i &= \mathbf{0} \quad (1 \le i \le N-1) \\
\frac{1}{\epsilon} \mathcal{D}^\ell: \begin{pmatrix} \mathbf{W}_{N-1} \\ \mathbf{W}_N \end{pmatrix} + \mathbf{F}_N &= \mathbf{0} \quad (i=N)
\end{aligned}
\quad (B.4) $$

The equilibrium of the terminal pairs of nodes $i=0$ and $i=N$ will be not enforced by using equations (B.4)$_{1,3}$, but by making the pseudo-energy stationary with respect to the boundary-layer amplitudes. The expressions of $\mathcal{D}^0$ and $\mathcal{D}^\ell$ are therefore not required. They are given for the sake of completeness.

## APPENDIX C. IDENTITIES RELATED TO INTEGRATION BY PARTS

In this Appendix, we integrate by parts the quantity $\int (\mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx = \int m_{11} \mathrm{f} \mathrm{g} + \epsilon(m_{12} \mathrm{f}' \mathrm{g} + m_{21} \mathrm{f} \mathrm{g}') + \epsilon^2(m_{13} \mathrm{f}'' \mathrm{g} + \cdots) + \cdots dx$, where $\mathrm{f}$ and $\mathrm{g}$ are generic functions of $x$. The goal is to remove the derivatives on $g$ present inside the integral, see (C.4). The identity thus obtained will be used in Appendix D to simplify the potential of the external load, see (D.26), and in Appendix G to rewrite the first variation of the strain energy in the bulk, see (G.2). The integrations by parts are carried out using the compact vector-of derivatives notation introduced in Section 2.1. We use the notation $e_i$ for the unit vectors in dimension 3,

$$ e_1 = (1 \ 0 \ 0), \quad e_2 = (0 \ 1 \ 0), \quad e_3 = (0 \ 0 \ 1). \quad (C.1) $$

### C.1. Single integration by parts

Given a 3×3 matrix $\mathbf{m}$ and a pair of functions $(f,g)$, one can integrate $\int_0^\ell (\mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx$ by parts as follows. Using the identity $I_3=e_1 \otimes e_1 + \mathcal{S}^\top \cdot \mathcal{S}$ with $\mathcal{S}$ defined in (2.7) and Equation (2.6), one has

$$ 
\begin{aligned}
\int_0^\ell (\mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx &= \int_0^\ell (\mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot (e_1 \otimes e_1 + \mathcal{S}^\top \cdot \mathcal{S}) \cdot \ddot{\mathbf{g}} dx \\
&= \int_0^\ell ((e_1 \otimes e_1 \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} + (\mathcal{S} \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot (\mathcal{S} \cdot \ddot{\mathbf{g}})) dx \\
&= \int_0^\ell ((e_1 \otimes e_1 \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} + \epsilon (\mathcal{S} \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}}') dx + \mathcal{O}(\epsilon^3) \\
&= \int_0^\ell ((e_1 \otimes e_1 \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} - \epsilon (\mathcal{S} \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}') \cdot \ddot{\mathbf{g}}) dx + \epsilon \left[ (\mathcal{S} \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} \right]_0^\ell + \mathcal{O}(\epsilon^3)
\end{aligned}
$$

Using Equation (2.6) again, this time with $g=f$, we have

$$ \int_0^\ell (\mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx = \int_0^\ell ((e_1 \otimes e_1 \cdot \mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx + \epsilon \left[ (\mathcal{S} \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} \right]_0^\ell + \mathcal{O}(\epsilon^3). \quad (C.2) $$

The third row of the matrix appearing in the integral is zero: $e_3 \cdot (e_1 \otimes e_1 \cdot \mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) = (e_1 \mathcal{f} / e_3) e_1 \cdot \mathbf{m} - (\mathcal{S}^\top \cdot e_3) \cdot \mathbf{m} \cdot \mathcal{S} = 0$. In the integral appearing in the right-hand side of (C.2), we have therefore gotten rid of the term proportional to the second derivative $\epsilon^2 \ddot{\mathbf{g}}' = \ddot{\mathbf{g}} \cdot e_3$ using an integration by parts. The integral is thus a linear combination of the function $g=\ddot{\mathbf{g}} \cdot e_1$ in undifferentiated form, and of its first derivative $\epsilon \ddot{\mathbf{g}}' = \ddot{\mathbf{g}} \cdot e_2$, see (2.1).

### C.2. A second integration by parts

Using the identity (C.2) twice, with the original $\mathbf{m}$ first, and with $\mathbf{m} \leftarrow (e_1 \otimes e_1 \cdot \mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S})$ next, we can integrate $\int_0^\ell (\mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{u}} dx$ by parts as

$$ 
\begin{aligned}
\int_0^\ell (\mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{u}} dx &= \int_0^\ell ((e_1 \otimes e_1 \cdot \mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx + \epsilon \left[ (\mathcal{S} \cdot \mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} \right]_0^\ell + \mathcal{O}(\epsilon^3) \\
&= \int_0^\ell \left( \begin{pmatrix} e_1 \otimes e_1 \cdot (e_1 \otimes e_1 \cdot \mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) \\ -\mathcal{S} \cdot (e_1 \otimes e_1 \cdot \mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) \cdot \mathcal{S} \end{pmatrix} \cdot \ddot{\mathbf{f}} \right) \cdot \ddot{\mathbf{g}} dx \\
&\quad + \epsilon \left[ (\mathcal{S} \cdot (\mathbf{m} + e_1 \otimes e_1 \cdot \mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} \right]_0^\ell + \mathcal{O}(\epsilon^3),
\end{aligned}
$$

where the two terms that are crossed out are zero since $\mathcal{S} \cdot e_1 = 0$ (see (2.7)).

Since the image of $\mathcal{S}^2 = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}$ is contained in the subspace spanned by $e_1$, we can insert the projector $e_1 \otimes e_1$ and rewrite $\mathcal{S}^2 = e_1 \otimes e_1 \cdot \mathcal{S}^2$. The last term in the integral can therefore be rewritten as $-\mathcal{S} \cdot (-\mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) \cdot \mathcal{S} = \mathcal{S}^2 \cdot \mathbf{m} \cdot \mathcal{S}^2 = e_1 \otimes e_1 \cdot \mathcal{S}^2 \cdot \mathbf{m} \cdot \mathcal{S}^2$. Using the identity $(e_1 \otimes e_1)^2 = e_1 \otimes e_1$ for the projector $e_1 \otimes e_1$ and factoring it across all the non-zero terms in the integral, we have

$$ 
\begin{aligned}
\int_0^\ell (\mathbf{m} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx &= \int_0^\ell ((e_1 \otimes e_1 \cdot (\mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S} + \mathcal{S}^2 \cdot \mathbf{m} \cdot \mathcal{S}^2)) \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx \\
&\quad + \epsilon \left[ (\mathcal{S} \cdot (\mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} \right]_0^\ell + \mathcal{O}(\epsilon^3) \\
&= \int_0^\ell (e_1 \otimes \mathbf{M} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} dx + \epsilon \left[ (\mathcal{S} \cdot (\mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S}) \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{g}} \right]_0^\ell + \mathcal{O}(\epsilon^3)
\end{aligned}
$$

where $\mathbf{M}$ is the vector of dimension 3,

$$ \mathbf{M} = (1 \ 0 \ 0) \cdot (\mathbf{m} - \mathcal{S} \cdot \mathbf{m} \cdot \mathcal{S} + \mathcal{S}^2 \cdot \mathbf{m} \cdot \mathcal{S}^2). \quad (C.3) $$
Using the identity $(e_1 \otimes M \cdot \ddot{f}) \cdot \ddot{g} = (M \cdot \ddot{f}) (e_1 \cdot \ddot{g}) = (M \cdot \ddot{f}) g$ in the integral term, and the approximation (2.2) for $\ddot{f}$ in the boundary term, this yields

$$ \int_0^\ell (m \cdot \ddot{f}) \cdot \ddot{g} dx = \int_0^\ell (M \cdot \ddot{f}) g dx + \left[ (\mathcal{S} \cdot (m - \mathcal{S} \cdot m \cdot \mathcal{S}) \cdot \mathcal{P}^\top \cdot (\varepsilon \dot{f})) \cdot \ddot{g} \right]_0^\ell + \mathcal{O}(\varepsilon^3). \quad (C.4) $$

In the integral in the right-hand side, we have gotten rid of the derivatives $g'$ and $g''$.

## APPENDIX D. ENERGY CONTINUALIZATION FOR THE SEMI-INFINITE CHAIN

In this Appendix, we analyze a semi-infinite lattice, having $N=\infty$ and $\ell=\infty$, and derive the continuous approximation of the total potential energy $\Phi_d$ announced in (7.3). We drop the boundary-layer mode associated with the endpoint at infinity, by setting $v^\ell=0$.

### D.1. Continualized strain energy

We start deriving a continuous approximation of the strain energy $\Phi_s$. We do so by inserting the expression (7.2) of the microscopic strain with $v^\ell=0$,

$$ \begin{aligned} E_j &= \mathcal{E}^{\text{bk}}(\mathcal{Y}_2) \cdot \ddot{u}'(\varepsilon(j+1/2)) + \mathcal{E}_j^0 v^0 + \mathcal{O}(\varepsilon^3) \\ &= (\mathcal{E}_j^0 \mathcal{E}^{\text{bk}}(\mathcal{Y}_2)) \cdot \left( \ddot{u}'(\varepsilon(j+1/2)) \right) + \mathcal{O}(\varepsilon^3), \end{aligned} \quad (D.1) $$

into the discrete strain energy (3.12)$_1$,

$$ \begin{aligned} \Phi_s &= \frac{\varepsilon}{2} E_0 \cdot Q^0 \cdot E_0 + \frac{\varepsilon}{2} \sum_{j=0}^{\infty} E_j \cdot Q^{\text{bk}} \cdot E_j \\ &= \frac{\varepsilon}{2} \left( \ddot{v}^0 \cdot \left( \ddot{u}'(\varepsilon/2) \right)^{\top} \cdot Q^0 \cdot \left( \mathcal{E}_0^0 \mathcal{E}^{\text{bk}}(\mathcal{Y}_2) \right)^{\top} \cdot \left( \ddot{v}^0 \cdot \left( \ddot{u}'(\varepsilon/2) \right)^{\top} \right) \right. \\ & \qquad \left. + \frac{\varepsilon}{2} \sum_{j=0}^{\infty} \left( \ddot{v}^0 \cdot \left( \ddot{u}'(\varepsilon(j+1/2)) \right)^{\top} \cdot \left( \begin{array}{cc} \mathcal{E}_j^0 \cdot Q^{\text{bk}} \cdot \mathcal{E}_j^0 & \mathcal{E}_j^0 \cdot Q^{\text{bk}} \cdot \mathcal{E}^{\text{bk}}(\mathcal{Y}_2) \\ \text{sym} & \bar{\mathcal{K}}(\mathcal{Y}_2) \end{array} \right) \cdot \left( \ddot{u}'(\varepsilon(j+1/2)) \right) \right) + \mathcal{O}(\varepsilon^3). \end{aligned} \quad (D.2) $$

where we have identified the matrix of homogenized moduli,

$$ \bar{\mathcal{K}}(\mathcal{Y}_2) = \mathcal{E}^{\text{bk}} \mathcal{T}^\top(\mathcal{Y}_2) \cdot Q^{\text{bk}} \cdot \mathcal{E}^{\text{bk}}(\mathcal{Y}_2). \quad (D.3) $$

As $\ddot{u}'=(\mathcal{O}(\varepsilon^0), \mathcal{O}(\varepsilon^1), \mathcal{O}(\varepsilon^2))$ by (2.1)$_2$, the components $\bar{\mathcal{K}}_{23}(\mathcal{Y}_2)$ and $\bar{\mathcal{K}}_{33}(\mathcal{Y}_2)$ located below the anti-diagonal of the symmetric matrix $\bar{\mathcal{K}}(\mathcal{Y}_2)$ contribute to the energy to order $\varepsilon^{1+2}=\varepsilon^3$ and $\varepsilon^{2+2}=\varepsilon^4$, respectively, and can thus be neglected. In the following, we will therefore replace $\bar{\mathcal{K}}(\mathcal{Y}_2)$ by its truncation

$$ \bar{\mathcal{K}} = \begin{pmatrix} \bar{\mathcal{K}}_{11} & \bar{\mathcal{K}}_{12} & \bar{\mathcal{K}}_{13} \\ \bar{\mathcal{K}}_{12} & \bar{\mathcal{K}}_{22} & 0 \\ \bar{\mathcal{K}}_{13} & 0 & 0 \end{pmatrix} = \begin{pmatrix} \bar{\mathcal{K}}_{11}(\mathcal{Y}_2) & \bar{\mathcal{K}}_{12}(\mathcal{Y}_2) & \bar{\mathcal{K}}_{13}(\mathcal{Y}_2) \\ \bar{\mathcal{K}}_{12}(\mathcal{Y}_2) & \bar{\mathcal{K}}_{22}(\mathcal{Y}_2) & 0 \\ \bar{\mathcal{K}}_{13}(\mathcal{Y}_2) & 0 & 0 \end{pmatrix}, \quad (D.4) $$

where $\bar{\mathcal{K}}_{ij}=\bar{\mathcal{K}}_{ij}$ denotes a component of $\bar{\mathcal{K}}(\mathcal{Y}_2)$ on or above the anti-diagonal ($i+j \le 4$). It is critical to note that these particular components $\bar{\mathcal{K}}_{ij}=\bar{\mathcal{K}}_{ij}$, with $i+j \le 4$, are all independent of the higher-order homogenization constant $\mathcal{Y}_2$, as reflected in our notation in (D.4). The demonstration is given in Section E.1 below. The matrix $\bar{\mathcal{K}}$ is thus entirely known in the context of second-order homogenization.

Since the microscopic moduli are positive, $Q^{\text{bk}} > 0$ by (3.13)$_1$, and the symmetric matrix $\bar{\mathcal{K}}(\mathcal{Y}_2)$ in (D.3) is positive too, for any value of $\mathcal{Y}_2$. In particular, the determinants of the two square sub-blocks $(\bar{\mathcal{K}}_{11})=(\bar{\mathcal{K}}_{11})$ and $(\bar{\mathcal{K}}_{12} \bar{\mathcal{K}}_{22})=(\bar{\mathcal{K}}_{11} \bar{\mathcal{K}}_{12} \bar{\mathcal{K}}_{22})$ on the diagonal of $\bar{\mathcal{K}}$ are both positive:

$$ \bar{\mathcal{K}}_{11} > 0, \quad \bar{\mathcal{K}}_{11} \bar{\mathcal{K}}_{22} - (\bar{\mathcal{K}}_{12})^2 \ge 0. \quad (D.5) $$

The strict positivity of $\bar{\mathcal{K}}_{11}$ follows from (3.3)$_1$. The determinant of the full $3 \times 3$ matrix $\bar{\mathcal{K}}(\mathcal{Y}_2)$ is positive as well, but this information is useless as it depends on the homogenization constant $\mathcal{Y}_2$ that remains indeterminate in the context of second-order homogenization.

We return to the expansion of the energy and split the right-hand side of (D.2) into three terms, namely a contribution coming from the boundary elasticity $Q^0$ (first term in the right-hand side), a second one mixing the boundary layer and the macroscopic field (top-left diagonal block and off-diagonal blocks in the infinite sum), and a third one coming from homogenization (lower-right diagonal block, in which we replace $\bar{\mathcal{K}}(\mathcal{Y}_2)$ with $\bar{\mathcal{K}}$),

$$ \begin{aligned} \Phi_s &= \Phi_s^{\text{be}} + \Phi_s^{\text{mx}} + \Phi_s^{\text{hm}} + \mathcal{O}(\varepsilon^3) \\ \Phi_s^{\text{be}} &= \frac{\varepsilon}{2} \left( \ddot{v}^0 \cdot \left( \ddot{u}'(\varepsilon/2) \right)^{\top} \cdot Q^0 \cdot \left( \mathcal{E}_0^0 \mathcal{E}^{\text{bk}}(\mathcal{Y}_2) \right)^{\top} \cdot \left( \ddot{v}^0 \cdot \left( \ddot{u}'(\varepsilon/2) \right)^{\top} \right) \right) \\ \Phi_s^{\text{mx}} &= \frac{\varepsilon}{2} \sum_{j=0}^{\infty} \left( \ddot{v}^0 \cdot \left( \ddot{u}'(\varepsilon(j+1/2)) \right)^{\top} \cdot \left( \begin{array}{cc} \mathcal{E}_j^0 \cdot Q^{\text{bk}} \cdot \mathcal{E}_j^0 & \mathcal{E}_j^0 \cdot Q^{\text{bk}} \cdot \mathcal{E}^{\text{bk}}(\mathcal{Y}_2) \\ \text{sym} & \mathbf{0} \end{array} \right) \cdot \left( \ddot{v}^0 \cdot \left( \ddot{u}'(\varepsilon(j+1/2)) \right)^{\top} \right) \right) \\ \Phi_s^{\text{hm}} &= \frac{\varepsilon}{2} \sum_{j=0}^{\infty} \ddot{u}'(\varepsilon(j+1/2)) \cdot \bar{\mathcal{K}} \cdot \ddot{u}'(\varepsilon(j+1/2)) \end{aligned} \quad (D.6) $$
We process the three terms as follows:

* Using (2.9) and (2.2), we approximate $\ddot{u}'(\epsilon/2) = \exp(\mathcal{S}/2) \cdot \ddot{u}'(0) + \mathcal{O}(\epsilon^3) = \exp(\mathcal{S}/2) \cdot \mathcal{P}^T \cdot \dot{u}'(0) + \mathcal{O}(\epsilon^3)$. This allows the boundary elasticity term to be rewritten as,

$$ \Phi_s^{\text{be}} = \frac{\epsilon}{2} \begin{pmatrix} \mathbf{v}^0 \\ \dot{\mathbf{u}}'(0) \end{pmatrix} \cdot \tilde{\mathbf{g}}_s^0 \begin{pmatrix} \mathbf{v}^0 \\ \dot{\mathbf{u}}'(0) \end{pmatrix} + \mathcal{O}(\epsilon^3), \quad (\text{D.7}) $$

where

$$ \tilde{\mathbf{g}}_{\text{be}}^{0} = \begin{pmatrix} \mathbf{I}_3 & 0 \\ 0 & e^{\mathcal{S}/2} \cdot \mathcal{P}^{\top} \end{pmatrix}^{\top} \cdot (\mathcal{E}_0^0 \mathcal{E}_{\text{bk}}(\mathcal{Y}_2))^{\top} \cdot \mathcal{Q}^0 \cdot (\mathcal{E}_0^0 \mathcal{E}_{\text{bk}}(\mathcal{Y}_2)) \cdot \begin{pmatrix} \mathbf{I}_3 & 0 \\ 0 & e^{\mathcal{S}/2} \cdot \mathcal{P}^{\top} \end{pmatrix}. \quad (\text{D.8}) $$

The quantity $\tilde{\mathbf{g}}_s^0$ is actually independent of $\mathcal{Y}_2$: this is proven in Section E.2 of the Appendix, and checked using symbolic calculations in the companion notebook.

* We proceed with the second contribution $\Phi_s^{\text{mx}}$ in (D.6)$_3$. Using (2.9) and (2.2) again, this time with $x=0$ and $\alpha=j+1/2$, we have

$$ \begin{aligned} \ddot{u}'(\epsilon(j+1/2)) &= e^{(j+1/2)\mathcal{S}} \cdot \ddot{u}'(0) + \mathcal{O}((\epsilon j)^3) \\ &= e^{(j+1/2)\mathcal{S}} \cdot \mathcal{P}^{\top} \cdot \ddot{u}'(0) + \mathcal{O}((\epsilon j)^2). \end{aligned} \quad (\text{D.9}) $$

The remainder $\mathcal{O}((\epsilon j)^2)$ comes from the fact that the right-hand side is a series in powers of $(\epsilon j)$, see (2.9), that has been truncated beyond the linear term. Now, the matrix appearing in the sum of (D.6)$_3$ decays exponentially for $j \gg 1$, since $\mathcal{E}_j^0$ decays exponentially by (6.15) and (6.13). Each term in the sum thus brings in a remainder of order $(\epsilon j)^2$ times an exponentially small quantity. When the sum over $j$ is carried out, the final remainder is therefore $\mathcal{O}(\epsilon^2)$. It becomes $\mathcal{O}(\epsilon^3)$ when multiplied by the overall factor $\epsilon/2$, which is comparable to the allowed error.

In view of this, we can rewrite $\Phi_s^{\text{mx}}$ as

$$ \Phi_s^{\text{mx}} = \frac{\epsilon}{2} \begin{pmatrix} \mathbf{v}^0 \\ \dot{\mathbf{u}}'(0) \end{pmatrix} \cdot \tilde{\mathbf{g}}_s^0 \begin{pmatrix} \mathbf{v}^0 \\ \dot{\mathbf{u}}'(0) \end{pmatrix} + \mathcal{O}(\epsilon^3), \quad (\text{D.10}) $$

where

$$ \tilde{\mathbf{g}}_{\text{mx}}^{0} = \sum_{j=0}^{\infty} \begin{pmatrix} \mathbf{I}_3 & 0 \\ 0 & e^{(j+1/2)\mathcal{S}} \cdot \mathcal{P}^{\top} \end{pmatrix}^{\top} \cdot \begin{pmatrix} (\mathcal{E}_j^0)^{\top} \cdot \mathcal{Q}^{\text{bk}} \cdot \mathcal{E}_j^0 & \mathcal{E}_j^0 \cdot \mathcal{Q}^{\text{bk}} \cdot \mathcal{E}_{\text{bk}}(\mathcal{Y}_2) \\ \text{sym} & \mathbf{0} \end{pmatrix} \cdot \begin{pmatrix} \mathbf{I}_3 & 0 \\ 0 & e^{(j+1/2)\mathcal{S}} \cdot \mathcal{P}^{\top} \end{pmatrix} \quad (\text{D.11}) $$

As the matrix exponential $e^{(j+1/2)\mathcal{S}}$ is actually a polynomial of degree 2, see (2.8), it is possible to calculate the boundary stiffness $\tilde{\mathbf{g}}_{\text{mx}}^0$ in symbolic form. This is done in the companion notebook. The matrix $\tilde{\mathbf{g}}_{\text{mx}}^0$ does not depend on $\mathcal{Y}_2$, as proven in Section E.3 and checked in the companion notebook.

* The term $\Phi_s^{\text{hm}}$ in (D.6)$_4$ can be calculated by setting $(\mathbf{n}, \alpha, p, q) = (\frac{1}{2}\bar{\mathbf{K}}, \frac{1}{2}u', u')$ in the Euler–Maclaurin formula (2.13), and using (2.2) twice in the resulting boundary term,

$$ \Phi_s^{\text{hm}} = \int_0^\infty \frac{1}{2} \ddot{u}' \cdot \bar{\mathbf{K}} \cdot \ddot{u}' dx + \frac{\epsilon}{2} \dot{u}'(0) \cdot \bar{\mathbf{k}}_{\text{hmg}}^0 \cdot \dot{u}'(0) + \mathcal{O}(\epsilon^3), \quad (\text{D.12}) $$

where

$$ \bar{\mathbf{k}}_{\text{hmg}}^0 = \frac{1}{24} \mathcal{P} \cdot (\bar{\mathbf{K}} \cdot \mathcal{P} + \mathcal{S}^{\top} \cdot \bar{\mathbf{K}}) \cdot \mathcal{P}^{\top}. \quad (\text{D.13}) $$

We proceed to collect the contributions to $\Phi_s$ appearing in (D.7), (D.10) and (D.12). The symmetric matrices appearing in the right-hand sides are summed up, and we identify the three matrices $\bar{\mathbf{g}}^0$, $-\bar{\mathbf{j}}_u^0$ and $\bar{\mathbf{k}}^0$ appearing as sub-blocks in the sum, with respective dimensions 3×3, 3×2 and 2×2, as

$$ \tilde{\mathbf{g}}_{\text{be}}^{0} + \tilde{\mathbf{g}}_{\text{mx}}^{0} + \begin{pmatrix} \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 2} \\ \mathbf{0}_{2 \times 3} & \tilde{\mathbf{k}}_{\text{hmg}}^{0} \end{pmatrix} = \begin{pmatrix} \bar{\mathbf{g}}^{0} & -\bar{\mathbf{j}}_{u}^{0} \\ \text{sym} & \tilde{\mathbf{k}}^{0} \end{pmatrix}, \quad (\text{D.14}) $$

The $\bar{\mathbf{g}}^0$, for instance, is given by (D.8) and (D.11) and it is symmetric,

$$ \bar{\mathbf{g}}^0 = (\mathcal{E}_0^0)^{\top} \cdot \mathcal{Q}^0 \cdot \mathcal{E}_0^0 + \sum_{j=0}^{\infty} (\mathcal{E}_j^0)^{\top} \cdot \mathcal{Q}^{\text{bk}} \cdot \mathcal{E}_j^0 \ge 0. \quad (\text{D.15}) $$

The positivity of $\bar{\mathbf{g}}^0$ comes from positivity of the spring constants, see (3.2), (3.3) and (3.13).

In the notebook, we calculate all three matrices $\bar{\mathbf{g}}^0$, $-\bar{\mathbf{j}}_u^0$ and $\bar{\mathbf{k}}^0$ in symbolic form using (D.8), (D.11) and (D.13). Overall, this delivers a continuous approximation of the strain energy (D.6)$_1$, see (D.7), (D.10), (D.12) and (D.14),

$$ \begin{aligned} \Phi_s &= \int_0^\infty \frac{1}{2} \ddot{u}' \cdot \bar{\mathbf{K}} \cdot \ddot{u}' dx + \frac{\epsilon}{2} \begin{pmatrix} \mathbf{v}^0 \\ \dot{\mathbf{u}}'(0) \end{pmatrix} \cdot \begin{pmatrix} \bar{\mathbf{g}}^0 & -\bar{\mathbf{j}}_u^0 \\ \text{sym} & \tilde{\mathbf{k}}^0 \end{pmatrix} \cdot \begin{pmatrix} \mathbf{v}^0 \\ \dot{\mathbf{u}}'(0) \end{pmatrix} \\ &= \int_0^\infty \frac{1}{2} \ddot{u}' \cdot \bar{\mathbf{K}} \cdot \ddot{u}' dx + \frac{\epsilon}{2} \dot{u}'(0) \cdot \bar{\mathbf{k}}^0 \cdot \dot{u}'(0) + \epsilon \left( \frac{1}{2} \mathbf{v}^0 \cdot \bar{\mathbf{g}}^0 \cdot \mathbf{v}^0 - \mathbf{v}^0 \cdot \bar{\mathbf{j}}_u^0 \cdot \dot{u}'(0) \right) + \mathcal{O}(\epsilon^3) \end{aligned} \quad (\text{D.16}) $$

## D.2. Continualized potential of external load

We proceed to analyze the potential of the external load $\Phi_e$. Using a Taylor expansion, we start by rewriting the nodal forces in (3.10), (3.1) as

$$ F_i = \begin{pmatrix} p^0 \\ 0 \end{pmatrix} \delta_i + \epsilon \mathcal{F} \cdot \dot{\mathbf{f}}(\epsilon i). \quad (\text{D.17}) $$
where

$$
\begin{align}
\mathcal{F} &= \frac{1}{2} \begin{pmatrix} (1 \ 0 \ 0) \cdot \exp(-\mathcal{S}/4) \\ (1 \ 0 \ 0) \cdot \exp(+\mathcal{S}/4) \end{pmatrix} \nonumber \\
&= \begin{pmatrix} +\frac{1}{2} & -\frac{1}{8} & +\frac{1}{64} \\ +\frac{1}{2} & +\frac{1}{8} & +\frac{1}{64} \end{pmatrix}. \tag{D.18}
\end{align}
$$

The potential energy in (3.12)$_2$ can then be written with the help of (7.1)$_2$ and (D.17) as

$$
\begin{align}
\Phi_e &= -\sum_{i=0}^{\infty} \left( \left( \begin{pmatrix} p^0 \\ 0 \end{pmatrix} \right) \delta_i + \varepsilon \mathcal{F} \cdot \ddot{\mathbf{f}}(\varepsilon i) \right) \cdot (\mathcal{L}^{\text{bk}} \cdot \ddot{\mathbf{u}}(\varepsilon i) + \varepsilon \mathbf{v}^0 \cdot \mathcal{L}_i^0) + \mathcal{O}(\varepsilon^3) \nonumber \\
&= -\left( \begin{pmatrix} p^0 \\ 0 \end{pmatrix} \right) \cdot \mathcal{L}^{\text{bk}} \cdot \ddot{\mathbf{u}}(0) - \varepsilon \sum_{i=0}^{\infty} (\mathcal{H} \cdot \ddot{\mathbf{f}}(\varepsilon i)) \cdot \ddot{\mathbf{u}}(\varepsilon i) - \varepsilon \mathbf{v}^0 \cdot \mathcal{L}_i^0 \cdot \left( \begin{pmatrix} p^0 \\ 0 \end{pmatrix} \right) - \varepsilon^2 \mathbf{v}^0 \cdot \sum_{i=0}^{\infty} (\mathcal{L}_i^0 \cdot \mathcal{F} \cdot \ddot{\mathbf{f}}(\varepsilon i)) + \mathcal{O}(\varepsilon^3), \tag{D.19}
\end{align}
$$

where

$$
\mathcal{H} = (\mathcal{L}^{\text{bk}})^{\top} \mathcal{F}. \qquad (\text{D.20})
$$

In the last term in the right-hand side of (D.19), which is of order $\varepsilon^2$, we can use $\ddot{\mathbf{f}}(\varepsilon i) = \ddot{\mathbf{f}}(0) + \mathcal{O}(\varepsilon) = \mathbf{f}(0)(1 \ 0 \ 0) + \mathcal{O}(\varepsilon)$. This yields

$$
-\varepsilon^2 \mathbf{v}^0 \cdot \sum_{i=0}^{\infty} (\mathcal{L}_i^0 \cdot \mathcal{F} \cdot \ddot{\mathbf{f}}(\varepsilon i)) = -\varepsilon \mathbf{v}^0 \cdot \left( \sum_{i=0}^{\infty} \mathcal{L}_i^0 \right) \cdot \mathcal{F} \cdot (1 \ 0 \ 0) \varepsilon f(0) + \mathcal{O}(\varepsilon^3). \quad (D.21)
$$

In view of this, the potential energy in (D.19) can be rewritten as

$$
\Phi_e = -\varepsilon \sum_{i=0}^{\infty} (\mathcal{H} \cdot \ddot{\mathbf{f}}(\varepsilon i)) \cdot \ddot{\mathbf{u}}(\varepsilon i) - ((\mathcal{L}^{\text{bk}})^{\top} \cdot (1 \ 0 \ 0)) \cdot \ddot{\mathbf{u}}(0) p^0 - \varepsilon \mathbf{v}^0 \cdot \bar{\mathcal{J}}_f^0 \cdot \begin{pmatrix} p^0 \\ \varepsilon f(0) \end{pmatrix} + \mathcal{O}(\varepsilon^3), \quad (D.22)
$$

where $\bar{\mathcal{J}}_f^0$ is the 3×2 matrix

$$
\bar{\mathcal{J}}_f^0 = (\mathcal{L}_0^0 \cdot (1 \ 0 \ 0) \cdot (\sum_{i=0}^{\infty} \mathcal{L}_i^0) \cdot \mathcal{F} \cdot (1 \ 0 \ 0)). \quad (D.23)
$$

Rewriting the sum in (D.22) as an integral using the Euler–Maclaurin formula (2.13) with $(n, a, p, q) = (\mathcal{H}, 0, f, u)$ and using $\ddot{\mathbf{f}}(0) = \mathcal{P}^{\top} \cdot \ddot{\mathbf{f}}(0) + \mathcal{O}(\varepsilon^2)$ in the boundary term, see (2.2), we have

$$
\Phi_e = - \int_0^\infty (\mathcal{H} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{u}} dx - \left( \tilde{\mathbf{h}}^0 \cdot \begin{pmatrix} p^0 \\ \varepsilon \dot{\mathbf{f}}(0) \end{pmatrix} \right) \cdot \ddot{\mathbf{u}}(0) - \varepsilon \mathbf{v}^0 \cdot \bar{\mathcal{J}}_f^0 \cdot \begin{pmatrix} p^0 \\ \varepsilon f(0) \end{pmatrix} + \mathcal{O}(\varepsilon^3) \quad (D.24)
$$

where $\tilde{\mathbf{h}}^0$ is the 3×3 matrix defined in block-matrix notation as

$$
\tilde{\mathbf{h}}^0 = \left( (\mathcal{L}^{\text{bk}})^{\top} \cdot (1 \ 0 \ 0) \left| \left( \frac{\mathcal{H}}{2} - \frac{\mathcal{H} \cdot \mathcal{S} + \mathcal{S}^{\top} \cdot \mathcal{H}}{12} \right) \cdot \mathcal{P}^{\top} \right) \right). \quad (D.25)
$$

The last step is to integrate the integral term $-\int_0^\infty (\mathcal{H} \cdot \ddot{\mathbf{f}}) \cdot \ddot{\mathbf{u}} dx$ in (D.24) by parts twice. To do so, we use the identity (C.4) with $(\mathbf{m}, f, g, \ell) = (-\mathcal{H} \cdot f, u, \infty)$, which yields

$$
\Phi_e = - \int_0^\infty (\bar{\mathbf{H}} \cdot \ddot{\mathbf{f}}) u dx - \left( \tilde{\mathbf{h}}^0 \cdot \begin{pmatrix} p^0 \\ \varepsilon \dot{\mathbf{f}}(0) \end{pmatrix} \right) \cdot \ddot{\mathbf{u}}(0) - \varepsilon \mathbf{v}^0 \cdot \bar{\mathcal{J}}_f^0 \cdot \begin{pmatrix} p^0 \\ \varepsilon f(0) \end{pmatrix} + \mathcal{O}(\varepsilon^3) \quad (D.26)
$$

where

$$
\bar{\mathbf{H}} = (1 \ 0 \ 0) \cdot (\mathcal{H} - \mathcal{S} \cdot \mathcal{H} \cdot \mathcal{S} + \mathcal{S}^2 \cdot \mathcal{H} \cdot \mathcal{S}^2) \quad (D.27)
$$

$$
\tilde{\mathbf{h}}^0 = \tilde{\mathbf{h}}^0 + (\mathbf{0}_{3 \times 1} | -\mathcal{S} \cdot (\mathcal{H} - \mathcal{S} \cdot \mathcal{H} \cdot \mathcal{S}) \cdot \mathcal{P}^{\top}). \quad (D.27)
$$

D.3. Summary

Combining the strain energy $\Phi_s$ in (D.16) and the potential energy of the external load $\Phi_e$ in (D.26), we obtain the
continuous energy $\Phi = \Phi_s + \Phi_e$ of a semi-infinite lattice in the form

$$
\begin{align}
\tilde{\Phi}_{\varepsilon}[u, \mathbf{v}^b] ={}& \int_0^{\infty} \left( \frac{1}{2} \ddot{u}' \cdot \bar{\mathbf{K}} \cdot \ddot{u}' - (\bar{\mathbf{H}} \cdot \ddot{\mathbf{f}}) u \right) dx \nonumber \\
& + \frac{\varepsilon}{2} \ddot{u}'(0) \cdot \tilde{\mathbf{k}}^0 \cdot \ddot{u}'(0) - \left( \tilde{\mathbf{h}}^0 \cdot \begin{pmatrix} p^0 \\ \varepsilon \dot{\mathbf{f}}(0) \end{pmatrix} \right) \cdot \ddot{\mathbf{u}}(0) \nonumber \\
& + \varepsilon \left( \frac{1}{2} \mathbf{v}^0 \cdot \bar{\mathbf{y}}^0 \cdot \mathbf{v}^0 - \mathbf{v}^0 \cdot \begin{pmatrix} \bar{\mathcal{J}}_u^0 \\ \bar{\mathcal{J}}_f^0 \end{pmatrix} \cdot \begin{pmatrix} \dot{\mathbf{u}}'(0) \\ \varepsilon f(0) \end{pmatrix} \right) + \mathcal{O}(\varepsilon^3). \tag{D.28}
\end{align}
$$

After restoring the second endpoint, this yields the form of the continualized energy announced in (7.3).

APPENDIX E. QUANTITIES THAT DO NOT DEPEND ON $\mathcal{Y}_2$

E.1. Homogenized elastic moduli

In this section, we prove that the entries $\tilde{K}_{ij}(\mathcal{Y}_2)$ appearing on or above the anti-diagonal $(i+j \le 4)$ of the matrix of homogenized elastic moduli $\tilde{\mathbf{K}}(\mathcal{Y}_2)$ in (D.3) do not depend on $\mathcal{Y}_2$. This result has been announced in (D.4).
## APPENDIX F. OPTIMIZING THE BOUNDARY-LAYER AMPLITUDES

In this Appendix, we insert the optimal boundary-layer amplitudes $v^b$ given by (7.6) into the energy $\tilde{\Phi}_\varepsilon$ in (7.3). Our goal is to justify the expressions of the energy $\Phi_\varepsilon[u]$ announced in (7.7).

To start with, we calculate the optimal value of boundary term coming from the boundary layer in (7.3),

$$ \Phi_{\text{bl}} = \varepsilon \left( \frac{1}{2} \mathbf{v}^b \cdot \bar{\mathbf{g}}^b \cdot \mathbf{v}^b - \mathbf{v}^b \cdot (\bar{\mathbf{g}}_u^b | \bar{\mathbf{g}}_f^b) \cdot \left( \frac{\dot{\mathbf{u}}'(b)}{\varepsilon f(b)} \right) \right). \quad (\text{F.1}) $$

Inserting the value of $v^b$ given by (7.6), and using the fact that $\bar{\mathbf{g}}^b$ is symmetric, we have

$$ \begin{aligned} \Phi_{\text{bl}} &= \varepsilon \left( \frac{\dot{\mathbf{u}}'(b)}{\varepsilon f(b)} \right) \cdot (\bar{\mathbf{g}}_u^b | \bar{\mathbf{g}}_f^b)^\top \cdot (\bar{\mathbf{g}}^b)^{-1} \cdot \left( \frac{1}{2} \underbrace{\bar{\mathbf{g}}^b \cdot (\bar{\mathbf{g}}^b)^{-1}}_{\mathcal{Y}_3} - \mathbf{I}_3 \right) \cdot (\bar{\mathbf{g}}_u^b | \bar{\mathbf{g}}_f^b) \cdot \left( \frac{\dot{\mathbf{u}}'(b)}{\varepsilon f(b)} \right) + \mathcal{O}(\varepsilon^3) \\ &= -\frac{\varepsilon}{2} \left( \frac{\dot{\mathbf{u}}'(b)}{\varepsilon f(b)} \right) \cdot (\bar{\mathbf{g}}_u^b | \bar{\mathbf{g}}_f^b)^\top \cdot (\bar{\mathbf{g}}^b)^{-1} \cdot (\bar{\mathbf{g}}_u^b | \bar{\mathbf{g}}_f^b) \cdot \left( \frac{\dot{\mathbf{u}}'(b)}{\varepsilon f(b)} \right) + \mathcal{O}(\varepsilon^3) \quad (\text{F.2}) \\ &= \frac{\varepsilon}{2} \left( \frac{\dot{\mathbf{u}}'(b)}{\varepsilon f(b)} \right) \cdot \begin{pmatrix} -(\bar{\mathbf{g}}_u^b)^\top \cdot (\bar{\mathbf{g}}^b)^{-1} \cdot \bar{\mathbf{g}}_u^b & -(\bar{\mathbf{g}}_u^b)^\top \cdot (\bar{\mathbf{g}}^b)^{-1} \cdot \bar{\mathbf{g}}_f^b \\ \text{sym} & -(\bar{\mathbf{g}}_f^b)^\top \cdot (\bar{\mathbf{g}}^b)^{-1} \cdot \bar{\mathbf{g}}_f^b \end{pmatrix} \cdot \left( \frac{\dot{\mathbf{u}}'(b)}{\varepsilon f(b)} \right) + \mathcal{O}(\varepsilon^3). \end{aligned} $$

Expanding the block-matrix multiplication, using the identities (2.5) to rewrite $\varepsilon\dot{\mathbf{u}}'(b)=\mathcal{R}\cdot\dot{\mathbf{u}}(b)$ and (2.1) to rewrite $\varepsilon f(b)=\varepsilon(1\ 0)\cdot\dot{\mathbf{f}}(b)$ in the cross-terms, inserting the result into (7.3) and using (7.4), we obtain the relaxed expression of the pseudo-energy $\Phi_\varepsilon[u]$ as

$$ \Phi_{\varepsilon}[u] = \int_{0}^{\ell} \left( \frac{1}{2} \ddot{\mathbf{u}}' \cdot \bar{\mathbf{K}} \cdot \ddot{\mathbf{u}}' - \mathbf{fu} \right) dx \\ + \sum_{b \in \{0, \ell\}} \left( \frac{\varepsilon}{2} \dot{\mathbf{u}}'(b) \cdot \bar{\mathbf{k}}^b \cdot \dot{\mathbf{u}}'(b) - \left( \bar{\mathbf{h}}^b \cdot \left( \frac{p^b}{\varepsilon f(b)} \right) \right) \ddot{\mathbf{u}}(b) - \frac{\varepsilon}{2} \left( \frac{p^b}{\varepsilon f(b)} \right) \cdot \bar{\mathbf{c}}^b \cdot \left( \frac{p^b}{\varepsilon f(b)} \right) + \mathcal{O}(\varepsilon^3) \right), \quad (\text{F.3}) $$

where

$$ \begin{aligned} \bar{\mathbf{k}}^b &= \tilde{\mathbf{k}}^b - (\bar{\mathbf{g}}_u^b)^\top \cdot (\bar{\mathbf{g}}^b)^{-1} \cdot \bar{\mathbf{g}}_u^b, \\ \bar{\mathbf{h}}^b &= \tilde{\mathbf{h}}^b + \mathcal{R}^\top \cdot (\bar{\mathbf{g}}_u^b)^\top \cdot (\bar{\mathbf{g}}^b)^{-1} \cdot \bar{\mathbf{h}}_f^b \cdot \begin{pmatrix} 1 & (0\ 0) \\ 0 & (1\ 0) \end{pmatrix}, \\ \bar{\mathbf{c}}^b &= (\bar{\mathbf{g}}_f^b)^\top \cdot (\bar{\mathbf{g}}^b)^{-1} \cdot \bar{\mathbf{g}}_f^b. \end{aligned} \quad (\text{F.4}) $$

The matrix $\bar{\mathbf{h}}^b$ identified above is multiplied by $(p^b \ \varepsilon f \ \varepsilon^2 f')$ on the right-hand side, and by $(u \ \varepsilon u' \ \varepsilon^2 u'')$ on the left-hand side. As a result, the entries which we can resolve in second-order homogenization are

$$ \bar{\mathbf{h}}^b = \begin{pmatrix} \bar{\mathbf{h}}_{11}^b & \bar{\mathbf{h}}_{12}^b & \bar{\mathbf{h}}_{13}^b \\ \bar{\mathbf{h}}_{21}^b & \bar{\mathbf{h}}_{22}^b & \bar{\mathbf{h}}_{23}^b \\ \bar{\mathbf{h}}_{31}^b & \bar{\mathbf{h}}_{32}^b & \bar{\mathbf{h}}_{33}^b \end{pmatrix}. \quad (\text{F.5}) $$

In the main text, we have extracted the top row of $\bar{\mathbf{h}}^b$ into a vector $\bar{\mathbf{t}}^b \in \mathbb{R}^3$ and in the two remaining rows, we have extracted a 2x2 matrix $\bar{\mathbf{r}}^b$,

$$ \bar{\mathbf{h}}^b = \begin{pmatrix} \bar{\mathbf{t}}_1^b & \bar{\mathbf{t}}_2^b & \bar{\mathbf{t}}_3^b \\ \bar{\mathbf{r}}_{11}^b & \bar{\mathbf{r}}_{12}^b & \bar{\mathbf{r}}_{13}^b \\ \bar{\mathbf{r}}_{21}^b & \bar{\mathbf{r}}_{22}^b & \bar{\mathbf{r}}_{23}^b \end{pmatrix}, \quad (\text{F.6}) $$

This delivers the equivalent expression of the potential energy associated with boundary tractions which we use in (7.7),

$$ \left( \bar{\mathbf{h}}^b \cdot \left( \frac{p^b}{\varepsilon f(b)} \right) \right) \ddot{\mathbf{u}}(b) = \bar{\mathbf{t}}^b \cdot \left( \frac{p^b}{\varepsilon f(b)} \right) u(b) + \varepsilon \left( \bar{\mathbf{r}}^b \cdot \left( \frac{p^b}{\varepsilon f(b)} \right) \right) \dot{\mathbf{u}}'(b) + \mathcal{O}(\varepsilon^3). \quad (\text{F.7}) $$

## APPENDIX G. FIRST VARIATION OF THE EFFECTIVE STRAIN ENERGY

In this appendix, the first variation $D_u\Phi_\varepsilon[u, \delta u]=0$ of the pseudo-energy $\Phi_\varepsilon[u]$ in (7.7) is worked out. We carry out the integrations by parts that make it possible to identify the equilibrium problem.

### G.1. First variation of the strain energy in the bulk

In this section, we start by working out the expression of the first variation of the strain energy in the bulk,

$$ \int_{0}^{\ell} (\bar{\mathbf{K}} \cdot \ddot{\mathbf{u}}') \cdot \delta \ddot{\mathbf{u}}' dx. \quad (\text{G.1}) $$
We integrate it by part twice using the identity (C.4) with $(m,f,g)=(\bar{\mathbf{K}},u',\delta u')$:

$$ \int_0^\ell (\bar{\mathbf{K}} \cdot \ddot{\mathbf{u}}') \cdot \delta \ddot{\mathbf{u}}' dx = \int_0^\ell (\bar{\mathbf{A}} \cdot \ddot{\mathbf{u}}') \delta u' dx + \left[ (\mathcal{P} \cdot (\bar{\mathbf{K}} - \mathcal{P} \cdot \bar{\mathbf{K}}) \cdot \mathcal{P}) \cdot \mathcal{P}^\top \cdot (\varepsilon \dot{\mathbf{u}}') \cdot \delta \dot{\mathbf{u}}' \right]_0^\ell + \mathcal{O}(\varepsilon^3). \quad (G.2) $$

where $\bar{\mathbf{A}} \in \mathbb{R}^3$ is the quantity denoted as $\mathbf{M}$ in Appendix C. It is given by (C.3) as

$$ \begin{aligned} \bar{\mathbf{A}} &= (1 \ 0 \ 0) \cdot (\bar{\mathbf{K}} - \mathcal{P} \cdot \bar{\mathbf{K}} \cdot \mathcal{P} + \mathcal{P}^2 \cdot \bar{\mathbf{K}} \cdot \mathcal{P}^2) \\ &= (\bar{\mathbf{K}}_{11} \ (\bar{\mathbf{K}}_{12} - \bar{\mathbf{K}}_{21}) \ \bar{\mathbf{K}}_{13} - \bar{\mathbf{K}}_{22} + \bar{\mathbf{K}}_{31}) \\ &= (\bar{\mathbf{K}}_{11} \ 0 \ -(\bar{\mathbf{K}}_{22} - 2\bar{\mathbf{K}}_{13})). \end{aligned} \quad (G.3) $$

The boundary term in (G.2) is rewritten as

$$ [f]_0^\ell = f(\ell) - f(0) = \sum_{b \in \{0, \ell\}} s^b f(b), \quad (G.4) $$

where $s^b$ is the outward unit normal,

$$ s^b = \begin{cases} -1 & \text{if } b=0 \\ +1 & \text{if } b=\ell. \end{cases} \quad (G.5) $$

Besides, we note that the matrix $\mathcal{P}$ in (2.7), $\mathcal{P}$ in (2.3) and $\mathcal{R}$ in (2.4) satisfy

$$ \mathcal{P} \cdot \mathcal{P} = \mathcal{P} \cdot \mathcal{P}^\top \cdot \mathcal{R} = \mathcal{R}. \quad (G.6) $$

Using these relations, together with the identity (2.2) with $g=\delta u'$, we rewrite the boundary term in (G.2) as

$$ \begin{aligned} \left[ (\mathcal{P} \cdot (\bar{\mathbf{K}} - \mathcal{P} \cdot \bar{\mathbf{K}}) \cdot \mathcal{P}) \cdot \mathcal{P}^\top \cdot (\varepsilon \dot{\mathbf{u}}') \cdot \delta \dot{\mathbf{u}}' \right]_0^\ell &= \varepsilon \sum_{b \in \{0, \ell\}} s^b (\mathcal{P} \cdot (\bar{\mathbf{K}} - \mathcal{P} \cdot \bar{\mathbf{K}}) \cdot \mathcal{P}) \cdot \mathcal{P}^\top \cdot \dot{\mathbf{u}}' (b)) \cdot (\mathcal{P}^\top \cdot \delta \dot{\mathbf{u}}' (b)) \\ &= \varepsilon \sum_{b \in \{0, \ell\}} s^b (\mathcal{P} \cdot \mathcal{P} \cdot (\bar{\mathbf{K}} - \mathcal{P} \cdot \bar{\mathbf{K}}) \cdot \mathcal{P}) \cdot \mathcal{P}^\top \cdot \dot{\mathbf{u}}' (b)) \cdot \delta \dot{\mathbf{u}}' (b) \quad (G.7) \\ &= \varepsilon \sum_{b \in \{0, \ell\}} (s^b \mathbf{l}(\bar{\mathbf{K}}) \cdot \dot{\mathbf{u}}' (b)) \cdot \delta \dot{\mathbf{u}}' (b), \end{aligned} $$

where $\mathbf{l}(\bar{\mathbf{K}})$ is the 2x2 matrix

$$ \mathbf{l}(\bar{\mathbf{K}}) = \mathcal{R} \cdot (\bar{\mathbf{K}} - \mathcal{P} \cdot \bar{\mathbf{K}} \cdot \mathcal{P}) \cdot \mathcal{P}^\top. \quad (G.8) $$

Inserting this into (G.2) and using (G.3), we obtain the first variation of the strain energy in the bulk in the form

$$ \int_0^\ell (\bar{\mathbf{K}} \cdot \ddot{\mathbf{u}}') \cdot \delta \ddot{\mathbf{u}}' dx = \int_0^\ell (\bar{\mathbf{A}} \cdot \ddot{\mathbf{u}}') \delta u' dx + \varepsilon \sum_{b \in \{0, \ell\}} (s^b \mathbf{l}(\bar{\mathbf{K}}) \cdot \dot{\mathbf{u}}' (b)) \cdot \delta \dot{\mathbf{u}}' (b) + \mathcal{O}(\varepsilon^3). \quad (G.9) $$

## G.2. First variation of the complete pseudo-energy

Adding to (G.9) the first variation of the other terms in (7.7), we get

$$ D_u \Phi_{\varepsilon}[u, \delta u] = \int_0^{\ell} ((\bar{\mathbf{A}} \cdot \ddot{\mathbf{u}}') \delta u' - f \delta u) dx + \sum_{b \in \{0, \ell\}} \left( \varepsilon (\bar{\mathbf{k}}^b + s^b \mathbf{l}(\bar{\mathbf{K}}) \cdot \dot{\mathbf{u}}' (b)) \cdot \delta \dot{\mathbf{u}}' (b) - \bar{\mathbf{t}}^b \cdot (p^b \varepsilon \dot{f}(b)) \delta u(b) + \varepsilon (\bar{\mathbf{r}}^b \cdot (p^b \varepsilon f(b))) \cdot \delta \dot{\mathbf{u}}' (b) \right) + \mathcal{O}(\varepsilon^3). \quad (G.10) $$

Identifying $\bar{\mathbf{A}} \cdot \ddot{\mathbf{u}}'$ as the internal force $n$,

$$ n(x) = \bar{\mathbf{A}} \cdot \ddot{\mathbf{u}}' (x), \quad (G.11) $$

and integrating by parts the term $(\bar{\mathbf{A}} \cdot \ddot{\mathbf{u}}') \delta u' = n \delta u'$ in (G.10), we get

$$ D_u \Phi_{\varepsilon}[u, \delta u] = - \int_0^{\ell} (n' + f) \delta u dx - \sum_{b \in \{0, \ell\}} \left( \bar{\mathbf{t}}^b \cdot (p^b \varepsilon \dot{f}(b)) - s^b n(b) \delta u(b) + \varepsilon (\bar{\mathbf{r}}^b \cdot (p^b \varepsilon f(b))) - (\bar{\mathbf{k}}^b + s^b \mathbf{l}(\bar{\mathbf{K}})) \cdot \dot{\mathbf{u}}' (b) \right) \delta \dot{\mathbf{u}}' (b) + \mathcal{O}(\varepsilon^3). \quad (G.12) $$

## APPENDIX H. ACCEPTABLE ERROR ON MICROSCOPIC DISPLACEMENT

We have generated the prediction $u$ for the macroscopic displacement by making the pseudo-energy stationary order by order in $\varepsilon$, see (8.2). In future work [TAL25a], we will explore an alternative method for generating $u$ by a non-perturbative method. This alternative approach introduces additional error terms $\Delta u(x)$ in the expression of $u(x)$ that depend (mildly) on the fast variable $x/\varepsilon$. In this Appendix, we characterize the perturbations $\Delta u$ that preserve the convergence rate $e_\Phi = \mathcal{O}(\varepsilon^4)$. As discussed in Section 10.2, this convergence rate ensures the correctness of the relocalized solution to second order in $\varepsilon$.

More specifically, we establish the following sufficient condition. If the macroscopic displacement is perturbed by a quantity $\Delta u(x)$ such that

$$ \| (\ddot{\Delta \mathbf{u}})' \|_2 = \mathcal{O}(\varepsilon^2) \quad \text{and} \quad \forall b \in \{0, \ell\}, \quad \| (\Delta \mathbf{u})'(b) \| = \mathcal{O}(\varepsilon^{3/2}), \quad (H.1) $$
where $\|\cdot\|_2$ represents the $L_2$-norm, then the perturbation to the relocalized solution is negligible in the sense that the relocalized solution will still satisfy $e_\Phi=\mathcal{O}(\varepsilon^4)$.

The proof is as follows. The perturbation $\Delta u$ modifies the relocalized strain $E_j$ in (7.9) in two ways:

* the bulk contribution is perturbed by $\mathcal{E}^{\text{bk}}(\mathbb{Y}_2) \cdot (\ddot{\Delta}u)'(x^j) = \mathcal{O}(\|(\ddot{\Delta}u)'(x^j)\|)$,

* the boundary layers are modified by $\sum_{b \in \{0, \ell\}} [\mathcal{E}_j^b \cdot \bar{\mathcal{E}}^b] \cdot (\Delta u)'(b) |_{p^b, \varepsilon_f(b)} = \mathcal{O}(\|(\Delta u)'(0)\| \exp(-x^j/\varepsilon_\lambda), \|(\Delta \dot{u})'(\ell)\| \exp(-\ell^{-x^j}/\varepsilon_\lambda))$, where $\lambda = \mathcal{O}(1)$ is the scaled length of the boundary layers entering into $\mathcal{E}_j^b$.

Altogether, the perturbation $\Delta u$ modifies the relocalized strain by no more that

$$ \|\Delta \mathbf{E}_j\|_2 = \mathcal{O}(\|(\ddot{\Delta} \mathbf{u})'\|_2, \|(\Delta \dot{\mathbf{u}})'(0)\| \varepsilon^{1/2}, \|(\Delta \dot{\mathbf{u}})'(\ell)\| \varepsilon^{1/2}), \quad (H.2) $$

where the coefficients $\varepsilon^{1/2}$ arise from the size $\sim \varepsilon$ of the exponential boundary layers. In the left-hand side of (H.2), we have used the discrete $L_2$-norm of the perturbation $\Delta \mathbf{E}_j$, defined as $\|\Delta \mathbf{E}_j\|_2 = (\frac{1}{N} \sum_{j=0}^{N-1} \|\Delta \mathbf{E}_j\|^2)^{1/2}$.

In view of the assumptions (H.1), we thus have $\|\Delta \mathbf{E}_j\|_2 = \mathcal{O}(\varepsilon^2)$. Assuming in addition that none of the spring constants is non-zero, it is easy to show that the discrete $L_2$-norm is equivalent to the strain-energy norm $\|\cdot\|_{\mathcal{Q}}$ introduced in (10.9)$_2$. This shows that $\|\Delta \mathbf{E}_j\|_{\mathcal{Q}} = \mathcal{O}(\varepsilon^2)$, i.e.,

$$ \|\Delta \mathbf{E}_j\|_{\mathcal{Q}}^2 = \mathcal{O}(\varepsilon^4). \quad (H.3) $$

In view of (10.9)$_2$, a perturbation $\Delta u$ satisfying (H.1) does not affect the convergence criterion $e_\Phi = \mathcal{O}(\varepsilon^4)$.