# Fixing non-positive energies in higher-order homogenization

Manon Thbaut, Basile Audoly, Claire Lestringant

► To cite this version:

Manon Thbaut, Basile Audoly, Claire Lestringant. Fixing non-positive energies in higher-order homogenization. Journal of the Mechanics and Physics of Solids, 2025, pp.106168. 10.1016/j.jmps.2025.106168. hal-04760885v2

HAL Id: hal-04760885

https://hal.science/hal-04760885v2

Submitted on 5 May 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fixing non-positive energies
in higher-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

## Abstract

Energy functionals produced by second-order homogenization of periodic elastic structures commonly feature negative gradient moduli. We show that this undesirable property is caused by the truncation of the energy expansion in powers of the small scale separation parameter. By revisiting Cholesky's LDLT decomposition, we propose an alternative truncation method that restores positivity while preserving the order of accuracy. We illustrate this method on a variety of periodic structures, both continuous and discrete, and derive compact analytical expressions of the homogenized energy that are positive and accurate to second order. The method can also cure the energy functionals produced by second-order *dimension reduction*, which suffer similar non-positivity issues. It naturally extends beyond second order.

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

## 1 Introduction

Asymptotic homogenization is a rigorous tool to derive the effective properties of a given periodic elastic microstructure [AS77, BP89, BLP11, San80]. It delivers the microscopic strain as a series expansion involving the successive gradients of a macroscopic displacement. The terms in this expansion depend on localization tensors that are obtained as solutions of a cascade of cell problems. At leading order, a standard effective Cauchy continuum is obtained. Pushed to higher-order, the method yields an equivalent continuum whose energy depends on the successive gradients of the macroscopic strain [Bou96, BS11, GK89, LM18]. This correction allows to capture size-effects that typically occur when macroscopic fields vary on the scale of a few cells [DLSS22, YAL24].

Our long-term goal is to approach asymptotic second-order homogenization using an energy perspective, i.e., to come up with a recipe producing an *effective energy* whose minimizer is a *macroscopic displacement*. Equally important is to identify relocalization formulas for the microscopic displacement, strain and stress, that connect the macroscopic displacement with the microscopic features of the original problem. The goal is to produce approximations of the initial problem within some prescribed accuracy, say at second order in the scale separation parameter. A natural procedure to generate the effective energy is to insert the relocalization formulas into the initial energy of the microscopic problem, and to eliminate the microscopic features by integrating over the period of the microstructure [LM18]. This program runs into two major difficulties.

First, boundary layers are ignored by the homogenization procedure, which is therefore only applicable in the bulk [Dum86, MV97]. An energy treatment of boundary layers is left for future work. We give up on identifying effective boundary conditions for the moment and limit attention to an infinite medium in the present paper.

The second stumbling block is that the effective second-order elastic stiffness is often non-positive [ABV16, AL23, DLSS22, LM18, SC00]. The corresponding second-order effective energy is therefore non-convex, which makes the equilibrium boundary-value problem mathematically ill-posed. Typically, spurious solutions having short-scale oscillations optimize the energy by making it large and negative [TAL24a]. The negative stiffness associated with short-scale perturbations leads to severe issues when the homogenized energy functional is minimized numerically [KA24]. The goal of this paper is to address these non-positivity issues.

Some remedies have been proposed in the literature. The so-called 'Boussinesq trick' is often applied to the strong form of the equilibrium, effectively turning the gradient stiffness into a positive one [ABV16]. The trick, however, introduces higher-order derivatives of the applied load, degrading its smoothness. Besides, in the special case of structures possessing 'floppy' modes, one can project the second-order stiffness tensor on the subspace of those floppy modes, where the classical elasticity predicts a zero stiffness and the gradient elasticity is therefore dominant [DLSS22]. This makes the second-order stiffness tensor positive, but only in the direction of floppy modes. Another option, proposed by [SC00], is to truncate the microscopic strain *only*, and derive the effective energy without any further truncation. The energy obtained in this way is always positive but a *second-order approximation* of the strain is required, which forces homogenization to be pushed to *third order*. Ultimately, the positivity is restored at the price of calculating a number of non-asymptotic higher-order energy terms whose expressions tend to be cumbersome.

In the present work, we propose an original procedure for truncating the energy, applicable to homogenization at any order. It preserves the positivity of the original 'microscopic' energy that gets homogenized. In the context of second-order homogenization, the effective energy obtained in this way is the sum of two terms: (i) a first term in which the leading-order homogenized elastic moduli are combined with a *generalized strain*, which we derive in the form of the usual macroscopic strain corrected by strain gradients, and (ii) a second term in which higher-order elastic moduli are combined with the strain-gradient. The procedure warrants that both the classical (leading-order) elastic moduli and the higher-order ones are non-negative. The proposed truncation method preserves asymptotic
correctness: the difference between the positive energy functional delivered by the new method and the non-positive one produced by the standard approach is not larger than the homogenization error.

The proposed truncation method is inspired from Cholesky decomposition, a popular numerical method used for solving linear systems of equations involving a positive matrix $t$, see [Str07]. The Cholesky method produces a diagonal matrix $D \ge 0$ having positive entries and a lower triangular matrix $L$ having 1's on the diagonal, such that the original matrix can be decomposed as $t = L \cdot D \cdot L^T$. Whereas the Cholesky method is primarily used as a numerical method, we revisit it here in the context of homogenization (see Section 5). Specifically, we apply the first few iterations of the Cholesky method to a block-matrix $t$ representing the asymptotic expansion of the energy, and observe that this produces a truncated energy possessing the desired properties of positivity and asymptotic correctness. The positive, block-diagonal entries in the matrix $D \ge 0$ are interpreted as homogenized moduli and the entries in the lower triangular matrix $L$ as coefficients defining the generalized strain.

The approach is illustrated with a variety of periodic structures, including continuous microstructures (Section 2) and discrete truss and beam lattices (Section 3 and 2.4). We also cover the case of non-centro-symmetric structures, whose homogenized energy couples strain and strain gradient (Section 6). In passing, we explain the empirical observation that strain-gradient moduli are never negative in *soft* structures, i.e., when the leading-order elasticity is degenerate [DLSS22, S18a, Dur22]. We further show that our method naturally extends to higher-order *dimension reduction* and fixes the negative higher-order effective moduli that are found in, e.g., the bending of elastic beams (Section 7). Finally, force profiles varying on a short scale are used in Section 8 as a benchmark for comparing our method to existing ones from the literature.

In all our examples, we follow the *energy-based* higher-order homogenization method developed in [Le15, AL23]. The macroscopic displacement is prescribed and the small-scale, oscillatory component of the displacement (microscopic shift) is found by optimizing the elastic energy. The potential of the external loading is ignored in this process as the external loading is coupled to the macroscopic displacement but not to the microscopic shift under standard scaling assumptions on the load magnitude—the external loading can be restored later, when the homogenized model is used to solve a structural problem, see the discussion in Section 3.1 in [AL23] and Section 8 of the present paper. Next, the solution for the microscopic shift is inserted into the elastic energy, delivering a homogenized energy functional $\Phi^*$. The present work is concerned with the final step, which is to truncate $\Phi^*$ in a way that preserves positivity.

The positive energy we construct is not meant to improve the predictions of second-order homogenization schemes. Instead, our goal is to identify a variational structure underpinning these schemes, which is not known to date. Such a variational structure would considerably facilitate their numerical implementation.

## 2 A laminate subjected to antiplane deformation

In this section we analyze a laminate, i.e., an elastic medium made up of layers having dissimilar elastic properties. In the interest of simplicity, we limit attention to antiplane deformation: the layers are periodic in the direction denoted as $i_1$ in Figure 2.1, and the displacement takes place along an orthogonal direction $i_3$ and is invariant along $i_3$. There is no major difficulty in lifting the assumption of antiplane elasticity but the calculations are more involved.

The beginning of our analysis closely follows the work of [Le15], where an asymptotic homogenization procedure has been used to derive the equivalent elastic properties up to second order in the scale separation parameter $\eta$. In Section 2.4, we recall the deficiencies of the homogenized energy functionals obtained in this prior work. In the following subsections, we discuss alternate truncation procedures, including one in Section 2.6 that preserves positivity and is simple to set up.

### 2.1 Geometry and elastic energy

<figure><img src="image_2.png" /><figcaption>**Figure 2.1.** Geometry of a laminate having periodicity $\eta$ in the direction $i_1$. The laminate is invariant in the direction $i_3$ perpendicular to the plane $(i_1, i_2)$ shown. The shear modulus of the phases (shown in blue and grey) are $\mu_I$ and $\mu_{II}$, respectively. We limit attention to the case of linear, antiplane elasticity: the displacement is along the transverse direction $i_3$, and invariant in this direction.</figcaption></figure>
We consider the infinite laminate sketched in Figure 2.1. It is periodic along $i_1$ with spatial period $\eta > 0$. The phases are made up of a linear isotropic material, with respective shear moduli $\mu_I > 0$ and $\mu_{II} > 0$. We denote as $a$ the fraction of the first phase, with $0 < a < 1$. The domain of the first phase is the set of points $\mathbf{x} = (x_1, x_2) \in \mathbb{R}^2$ such that

$$ \frac{x_1}{\eta} \in \chi \quad (\text{phase I}) \qquad (2.1) $$

where $\chi$ is the set of scaled coordinates $y = \eta^{-1}x_1$ covering phase I, as shown at the bottom of Figure 2.1,

$$ \chi = \bigcup_{i \in \mathbb{Z}} [i, i+a]. \qquad (2.2) $$

The set $\chi$ contains the real numbers $y$ of the form $i+b$ where $i$ is a signed integer and $b$ is the fractional part of $y$, with $0 \le b \le a < 1$.

We limit attention to out-of-plane loading, i.e., assume that the displacement $v$ is along the transverse direction $i_3$, and invariant in the transverse direction,

$$ v(x_1, x_2) = v(x_1, x_2) i_3. \qquad (2.3) $$

The in-plane gradient of the transverse displacement $v(x_1, x_2)$ is denoted as

$$ \mathbf{E}(x_1, x_2) = \nabla v(x_1, x_2) = \left( \begin{array}{c} \frac{\partial v}{\partial x_1}(x_1, x_2) \\ \frac{\partial v}{\partial x_2}(x_1, x_2) \end{array} \right) \in \mathbb{R}^2. \qquad (2.4) $$

The classical, linearized strain $\tilde{\mathbf{E}}$ from 3D continuum mechanics can be expressed in terms of the reduced strain vector $\mathbf{E}$ as

$$ \tilde{\mathbf{E}}(x_1, x_2) = \frac{1}{2} (\mathbf{E}(x_1, x_2) \otimes i_3 + i_3 \otimes \mathbf{E}(x_1, x_2)), \qquad (2.5) $$

where $\otimes$ denotes the outer product. With $\mathrm{tr}\,\tilde{\mathbf{E}} = 0$, the strain energy $\Phi$ of the laminate, per unit thickness in the transverse direction $i_3$, can be expressed as

$$ \begin{aligned} \Phi[\mathbf{E}] &= \frac{1}{2} \int_{\mathbb{R}^2} 2\mu(x_1) \tilde{\mathbf{E}} : \tilde{\mathbf{E}} \, d\mathbf{x} \\ &= \frac{1}{2} \int_{\mathbb{R}^2} \mu(x_1) \|\mathbf{E}(x_1, x_2)\|^2 d\mathbf{x} \\ &= \frac{1}{2} \int_{\mathbb{R}^2} \mathbf{E}(x_1, x_2) \cdot \mathbf{k} \left( \frac{x_1}{\eta} \right) \cdot \mathbf{E}(x_1, x_2) \, d\mathbf{x} \end{aligned} \qquad (2.6) $$

where $\mu(x_1)$ denotes the local shear modulus, $\mu(x_1) = \mu_I$ if $\eta^{-1}x_1 \in \chi$ and $\mu(x_1) = \mu_{II}$ if $\eta^{-1}x_1 \notin \chi$, and $\|\cdot\|$ is the Euclidean norm of the vector $\mathbf{E}$. The elastic stiffness $\mathbf{k}(y)$ is the periodic function with period 1 of the stretched variable $y = \eta^{-1}x_1$

$$ \mathbf{k}(y) = \mathbf{I}_2 \times \begin{cases} \mu_I & \text{if } y \in \chi, \\ \mu_{II} & \text{if } y \notin \chi. \end{cases} \qquad (2.7) $$

Here, $\mathbf{I}_2$ is the identity matrix of dimension 2. Since both $\mu_I$ and $\mu_{II}$ are strictly positive, the symmetric matrix $\mathbf{k}(y)$ is also strictly positive in the sense of quadratic forms,

$$ \forall y, \quad \mathbf{k}(y) > 0. \qquad (2.8) $$

## 2.2 Series solution for the microscopic strain

We now consider the limit of well-separated scales,

$$ \eta \ll 1, \qquad (2.9) $$

in which homogenization is applicable. We use formal two-scale expansions and seek the displacement $v(x_1, x_2)$ in (2.3) in the form of a series expansion having dependence on both the slow variable $\mathbf{x} = (x_1, x_2)$ and on the fast variable $y = \eta^{-1}x_1$,

$$ v(\mathbf{x}) = \sum_{i=0}^{+\infty} \eta^i v_i(\mathbf{x}, \eta^{-1}x_1). \qquad (2.10) $$

The 'coefficients' $v_i(\mathbf{x}, y)$ are assumed to be smooth functions of their arguments $\mathbf{x}$ and $y$, and can be sought as periodic functions of $y$ without loss of generality. In (2.10), we have introduced a single stretched variable $y = \eta^{-1}x_1$ as the laminate is periodic in the direction $i_1$.

The leading-order displacement $v_0(\mathbf{x}, y)$ can be identified with the macroscopic displacement $u(\mathbf{x})$ and does not depend on $y$. The higher-order contributions $v_i(\mathbf{x}, y)$ for $i > 0$ are found [BP89] as the solutions of successive problems formulated on the unit cell $y \in [0, 1]$, using asymptotic homogenization. They are proportional to the successive gradients of the macroscopic strain $\varepsilon(\mathbf{x})$: the solution derived in [Le15] can be written as

$$ v(\mathbf{x}) = u(\mathbf{x}) + \eta \Psi_0(\eta^{-1}x_1) \cdot \varepsilon(\mathbf{x}) + \eta^2 \Psi_1(\eta^{-1}x_1) \cdot \nabla \varepsilon(\mathbf{x}) + \eta^3 \Psi_2(\eta^{-1}x_1) \cdot \nabla^2 \varepsilon(\mathbf{x}) + \mathcal{O}(\eta^4). \qquad (2.11) $$
In our notation, the macroscopic (effective) strain $\epsilon$ and its gradients are presented as vectors that are defined in block notation as

$$ \epsilon(\mathbf{x}) = \nabla u(\mathbf{x}) = \left( \frac{\partial u}{\partial x_1}(\mathbf{x}) \right) \in \mathbb{R}^2, \quad \nabla \epsilon(\mathbf{x}) = \left( \frac{\partial \epsilon}{\partial x_1}(\mathbf{x}) \right) \in \mathbb{R}^4, \quad \nabla^2 \epsilon(\mathbf{x}) = \left( \frac{\partial \nabla \epsilon}{\partial x_1}(\mathbf{x}) \right) \in \mathbb{R}^8, \quad \nabla^2 \nabla \epsilon(\mathbf{x}) = \left( \frac{\partial^2 \epsilon}{\partial x_1 \partial x_2}(\mathbf{x}) \right) \in \mathbb{R}^8. \quad (2.12) $$

The macroscopic strain vector $\epsilon(\mathbf{x}) \in \mathbb{R}^2$ is a condensed representation of the classical strain tensor $\tilde{\mathbf{e}}(\mathbf{x})$ from continuum mechanics, and the latter can be expressed in terms of the former by a formula similar to (2.5). Note that some of the entries in the strain gradients are redundant, such as $(\nabla^2\epsilon)_3 = \partial[(\nabla\epsilon)_3]/\partial x_1 = \partial^2\epsilon_1/(\partial x_1 \partial x_2) = \partial[(\nabla\epsilon)_1]/\partial x_2 = (\nabla^2\epsilon)_5$.

The first relocalization vectors $\Psi_0(y) \in \mathbb{R}^2$ and $\Psi_1(y) \in \mathbb{R}^4$ appearing in the homogenization solution (2.11) are obtained in subsection 2.4.2 of [Le15] in explicit form as

$$ \begin{aligned} \Psi_0(y) &= \begin{pmatrix} \frac{(\mu_I - \mu_{II})}{(1-a)\mu_I + a\mu_{II}} & \times \begin{cases} (1-a)\left(\frac{a}{2} - y\right) & \text{if } y \in \chi \\ a\left(y - \frac{1+a}{2}\right) & \text{if } y \notin \chi \end{cases} \\ &-\frac{(\mu_I - \mu_{II})}{(1-a)\mu_I + a\mu_{II}} \times \begin{cases} (1-a)\left(\frac{ay}{2} - \frac{y^2}{2} + \frac{a(1-2a)}{12}\right) & \text{if } y \in \chi \\ a\left(\frac{a^2}{2} - \frac{1+a}{2}y + \frac{2n^2+3a+1}{12}\right) & \text{if } y \notin \chi \end{cases} \end{pmatrix} \\ \Psi_1(y) &= \begin{pmatrix} 0 \\ 0 \\ -(\mu_I - \mu_{II}) \left( \frac{a(1-a)}{12} \left( \frac{a^2}{\mu_I} - \frac{(1-a)^2}{\mu_{II}} \right) + \begin{cases} \frac{(1-a)}{\mu_I} \left( \frac{y^2}{2} - \frac{ay}{2} \right) & \text{if } y \in \chi \\ \frac{a}{\mu_{II}} \left( \frac{1+a}{2}y - \frac{y^2}{2} - \frac{a}{2} \right) & \text{if } y \notin \chi \end{cases} \right) \end{pmatrix}. \end{aligned} \quad (2.13) $$

By design, the operators $\Psi_i(y)$ are all $y$-periodic and their average is zero,

$$ \langle \Psi_0 \rangle_y = \mathbf{0}_2, \quad \langle \Psi_1 \rangle_y = \mathbf{0}_4, \quad \dots \quad (2.14) $$

where $\langle \cdot \rangle_y$ represents the average operator on a unit-cell

$$ \langle \cdot \rangle_y = \int_0^1 dy. \quad (2.15) $$

The microscopic strain $\mathbf{E}$ can be obtained by inserting the solution (2.11) for the microscopic displacement into (2.4). The result is of the form $\mathbf{E}(\mathbf{x}) = \mathbf{E}^*(\mathbf{x}, \eta^{-1}x_1)$ where $\mathbf{E}^*(\mathbf{x}, y)$ is the smooth function, defined in block-matrix notation as

$$ \begin{aligned} \mathbf{E}^*(\mathbf{x}, y) &= \mathbf{F}_0(y) \cdot \boldsymbol{\epsilon}(\mathbf{x}) + \eta \mathbf{F}_1(y) \cdot \nabla \boldsymbol{\epsilon}(\mathbf{x}) + \eta^2 \mathbf{F}_2(y) \cdot \nabla^2 \boldsymbol{\epsilon}(\mathbf{x}) + \mathcal{O}(\eta^3) \\ &= (\mathbf{F}_0(y) \quad \mathbf{F}_1(y) \quad \mathbf{F}_2(y)) \cdot \begin{pmatrix} \boldsymbol{\epsilon}(\mathbf{x}) \\ \eta \nabla \boldsymbol{\epsilon}(\mathbf{x}) \\ \eta^2 \nabla^2 \boldsymbol{\epsilon}(\mathbf{x}) \end{pmatrix} + \mathcal{O}(\eta^3). \end{aligned} \quad (2.16) $$

The strain relocalization matrices $\mathbf{F}_0, \mathbf{F}_1, \mathbf{F}_2, \dots$ have respective dimensions 2×2, 2×4, 2×8, ... and are given by

$$ \begin{aligned} \mathbf{F}_0(y) &= \mathbf{I}_2 + \mathbf{i}_1 \otimes \mathbf{\Psi}_0'(y) \\ \mathbf{F}_1(y) &= \mathbf{i}_1 \otimes (\mathbf{\Psi}_0(y) \mathbf{0}_2) + \mathbf{i}_2 \otimes (\mathbf{0}_2 \mathbf{\Psi}_0(y)) + \mathbf{i}_1 \otimes \mathbf{\Psi}_1'(y) \\ \mathbf{F}_2(y) &= \mathbf{i}_1 \otimes (\mathbf{\Psi}_1(y) \mathbf{0}_4) + \mathbf{i}_2 \otimes (\mathbf{0}_4 \mathbf{\Psi}_1(y)) + \mathbf{i}_1 \otimes \mathbf{\Psi}_2'(y), \end{aligned} \quad (2.17) $$

where $'$ denotes derivation with respect to $y$ and $(s, t)$ denotes the vector obtained by concatenating the vectors $s$ and $t$ in block-vector notation. In the right-hand side of (2.17), the terms in the first and second column come from differentiating (2.11) with respect to $\mathbf{x}$ and $\eta^{-1}x_1$ respectively. As indicated by the brace, we denote as $\mathbf{F}_2^*$ the contribution to $\mathbf{F}_2$ that does not depend on the higher-order displacement localizer $\mathbf{\Psi}_2$.

In view of (2.14), the averages of the relocalization matrices are $\langle \mathbf{F}_0 \rangle_y = \mathbf{I}_2$, $\langle \mathbf{F}_1 \rangle_y = \mathbf{0}$, $\langle \mathbf{F}_2 \rangle_y = \mathbf{0}$, etc. By (2.16), this ensures that the macroscopic strain coincides with the average of the microscopic strain with respect to the fast variable, $\boldsymbol{\epsilon}(\mathbf{x}) = \langle \mathbf{E}^*(\mathbf{x}, y) \rangle_y$.

## 2.3 Series expansion of the energy

In the limit of well-separated scales $\eta \to 0$ and up to boundary terms which we can ignore since the laminate covers the entire plane $\mathbb{R}^2$, the strain energy $\Phi[\mathbf{E}]$ in (2.6) can be integrated sequentially, over the fast variable $y$ first and over the slow variable $\mathbf{x}$ next,

$$ \begin{aligned} \Phi[\mathbf{E}] &= \frac{1}{2} \int_{\mathbb{R}^2} \mathbf{E}(\mathbf{x}) \cdot \mathbf{k} \left(\frac{x_1}{\eta}\right) \cdot \mathbf{E}(\mathbf{x}) \, d\mathbf{x} \\ &= \frac{1}{2} \int_{\mathbb{R}^2} \langle \mathbf{E}^*(\mathbf{x}, y) \cdot \mathbf{k}(y) \cdot \mathbf{E}^*(\mathbf{x}, y) \rangle_y d\mathbf{x}. \end{aligned} \quad (2.18) $$
This identity is at the heart of energy-based homogenization and is used, e.g., in Equation [2.5] of [Le15].

Inserting the series solution (2.16) into the right-hand side in equation above, we obtain the strain energy $\Phi^*[\varepsilon] = \Phi[\mathbf{E}]$ in the form of a functional of the macroscopic strain $\varepsilon$,

$$ \Phi^*[\varepsilon] = \int_{\mathbb{R}^2} \frac{1}{2} \begin{pmatrix} \varepsilon(\mathbf{x}) \\ \eta \nabla \varepsilon(\mathbf{x}) \\ \eta^2 \nabla^2 \varepsilon(\mathbf{x}) \end{pmatrix} \cdot \mathbf{t} \cdot \begin{pmatrix} \varepsilon(\mathbf{x}) \\ \eta \nabla \varepsilon(\mathbf{x}) \\ \eta^2 \nabla^2 \varepsilon(\mathbf{x}) \end{pmatrix} \mathrm{d}\mathbf{x} + \mathcal{O}(\eta^3) \quad (2.19) $$

where the matrix of homogenized elastic properties $\mathbf{t}$ is the symmetric matrix of dimensions 14 × 14 defined by

$$ \mathbf{t} = \langle (\mathbf{F}_0(y) \ \mathbf{F}_1(y) \ \mathbf{F}_2(y))^T \cdot \mathbf{k}(y) \cdot (\mathbf{F}_0(y) \ \mathbf{F}_1(y) \ \mathbf{F}_2(y)) \rangle_y \quad (2.20) $$

Note that the block matrix $(\mathbf{F}_0 \ \mathbf{F}_1 \ \mathbf{F}_2)$ is of size $2 \times (2+4+8) = 2 \times 14$, and that both contractions take place over the 2 rows of this matrix. The quantity $\mathbf{t}$ can be expressed in terms of $\mu_I$, $\mu_{II}$ and $a$ using (2.13) and (2.17) but the resulting expressions are cumbersome and are not included here.

The symmetric matrix $\mathbf{t}$ is clearly positive in view of the positivity of $\mathbf{k}(y)$ in (2.8),

$$ \mathbf{t} \ge 0, \quad (2.21) $$

and our goal will be to preserve this positivity in the final homogenized energy.

The laminate being mirror-symmetric in both the horizontal and vertical directions, the homogenized energy $\Phi^*$ can only couple gradients of the strain whose orders are either both even, or both odd: these are the only energy terms whose sign is preserved by the mirror symmetries. This is a classical result known as centro-symmetry in the context of periodic homogenization, see Appendix B of [Bou96] and [ABB09, TB96]. An explicit calculation of $\mathbf{t}$ using a symbolic calculation language confirms that the blocks in the matrix $\mathbf{t}$ that couple even and odd orders in the strain gradients are zero, so that $\mathbf{t}$ can be rewritten in block-matrix notation as

$$ \mathbf{t} = \begin{pmatrix} \mathbf{K} & \mathbf{0}_{2 \times 4} & \mathbf{C} \\ \mathbf{0}_{4 \times 2} & \mathbf{B} & \mathbf{0}_{4 \times 8} \\ \mathbf{C}^T & \mathbf{0}_{8 \times 4} & \mathbf{Z} \end{pmatrix}. \quad (2.22) $$

In view of (2.19), the diagonal terms $\mathbf{K}, \mathbf{B}, \mathbf{Z}$ are square matrices with respective dimensions $2 \times 2, 4 \times 4$ and $8 \times 8$ that give rise to energy contributions proportional to $\varepsilon \otimes \varepsilon$, $\eta^2 \nabla \varepsilon \otimes \nabla \varepsilon$ and $\eta^4 \nabla^2 \varepsilon \otimes \nabla^2 \varepsilon$, respectively, whereas the off-diagonal block $\mathbf{C}$ has dimensions $2 \times 8$ and gives rise to an energy contribution proportional to $\eta^2 \varepsilon \otimes \nabla^2 \varepsilon$.

Equation (2.19) can be rewritten in expanded form as

$$ \Phi^*[\varepsilon] = \int_{\mathbb{R}^2} \frac{1}{2} (\varepsilon \cdot \mathbf{K} \cdot \varepsilon + \eta^2 (\nabla \varepsilon \cdot \mathbf{B} \cdot \nabla \varepsilon + 2 \varepsilon \cdot \mathbf{C} \cdot \nabla^2 \varepsilon) + \eta^4 \nabla^2 \varepsilon \cdot \mathbf{Z} \cdot \nabla^2 \varepsilon) \mathrm{d}\mathbf{x} + \mathcal{O}(\eta^3). \quad (2.23) $$

The sub-blocks of $\mathbf{t}$ associated with the successive orders in $\eta$ are therefore arranged along lines parallel to the SW-NE diagonal, stacked form NW (order $\eta^0$) to SE (order $\eta^4$), using standard abbreviations for intercardinal directions:

* the contribution of order $\eta^0$ comes from the top-left block $\mathbf{K}$ which contains the results of leading-order (classical) homogenization,

* the diagonal $\mathbf{C}^T \mathbf{B} \mathbf{C}$ gives rise to energy contributions proportional to $\eta^2$,

* the lower-right block $\mathbf{Z}$ gives rise to energy contributions proportional to $\eta^4$—these contributions are negligible in view of the remainder $\mathcal{O}(\eta^3)$ in (2.23) but we keep it for a moment.

With help of (2.20), we obtain the following expressions of the blocks making up the matrix $\mathbf{t}$ in (2.20)

$$ \begin{array}{c|c|c} \text{order} & \text{diagonal blocks} & \text{off-diagonal blocks} \\ \hline \eta^0 & \mathbf{K} = \langle \mathbf{F}_0^T(y) \cdot \mathbf{k}(y) \cdot \mathbf{F}_0(y) \rangle_y & \\ \eta^2 & \mathbf{B} = \langle \mathbf{F}_1^T(y) \cdot \mathbf{k}(y) \cdot \mathbf{F}_1(y) \rangle_y & \mathbf{C} = \langle \mathbf{F}_0^T(y) \cdot \mathbf{k}(y) \cdot \mathbf{F}_2(y) \rangle_y \\ \eta^4 & \mathbf{Z} = \langle \mathbf{F}_2^T(y) \cdot \mathbf{k}(y) \cdot \mathbf{F}_2(y) \rangle_y & \end{array} \quad (2.24) $$

By (2.8) or (2.21), the symmetric blocks $\mathbf{K}$ and $\mathbf{B}$ appearing on the diagonal of $\mathbf{t}$ are positive in the sense of quadratic forms,

$$ \begin{aligned} \mathbf{K} &> 0 \\ \mathbf{B} &\ge 0. \end{aligned} \quad (2.25) $$

The negligible diagonal block $\mathbf{Z}$ is positive as well but this property will not be used. The strict positivity of $\mathbf{K}$ in (2.25) comes from the fact that $\mathbf{k}(y) > 0$ and $\det \mathbf{F}_0(y) = [(1-a) \mu_I + a \mu_{II}]^{-1} \times \begin{cases} \mu_{II} & \text{if } y \in \chi > 0, \\ \mu_I & \text{if } y \notin \chi \end{cases}$, implying that $\mathbf{F}_0$ is full-rank.

**Remark 2.1.** The blocks $\mathbf{C}$ and $\mathbf{Z}$ are expressed in (2.24) in terms of $\mathbf{F}_2$, which is itself given in (2.17)$_3$ in terms of a high-order relocalization vector $\mathbf{\Psi}_2(y)$. The expression of $\mathbf{\Psi}_2(y)$ is not included in Equation (2.13) and can only be obtained by pushing homogenization to a higher order which leads to cumbersome calculations. The dependence of $\mathbf{C}$ and $\mathbf{Z}$ on $\mathbf{\Psi}_2$ is therefore bad news, but the following arguments show that it can be ignored entirely:

* (C does not actually depend on $\mathbf{\Psi}_2$) By the same argument as used at the end of Section C.3 in the Appendix of [AL23], one can establish the identity $\langle \mathbf{F}_0^T(y) \cdot \mathbf{k}(y) \cdot (\mathbf{i}_1 \otimes \mathbf{\Psi}_2'(y)) \rangle_y = \mathbf{0}$ which is a property of the leading-order solution $\mathbf{F}_0$, and holds irrespective of the expression of $\mathbf{\Psi}_2(y)$. With this identity, one can eliminate $\mathbf{\Psi}_2$ entirely from the definition of $\mathbf{C}$, using $\mathbf{C} = \langle \mathbf{F}_0^T(y) \cdot \mathbf{k}(y) \cdot \mathbf{F}_2^T(y) \rangle_y$ where $\mathbf{F}_2^*$ is identified in (2.17).
* (The dependence of Z on Ψ₂ can be neglected) As noted earlier, (i) Z gives rise to negligible energy contributions, of order η⁴, and (ii) the only benefit of Z is to preserve the positivity of t, see Section 2.5 below. This positivity however follows from (2.20), no matter what expression of Ψ₂(y) is used. Ignoring Ψ₂(y) in the definition of Z, that is setting Z = (F̄₂T(y) · k(y) · F̄₂(y))y is therefore acceptable.

## 2.4 Naïve truncation of the energy

In this section, we present the standard truncation procedure used in the literature on high-order homogenization, including in [Le15] among a number of other works, and point out its deficiencies.

The most obvious way of truncating the energy in (2.23) is to drop the negligible term of order η⁴, which yields an approximation Φ†[ε] = Φ*[ε] + O(η³) consistent with second-order homogenization in the form

$$ \Phi^{\dagger}[\varepsilon] = \int_{\mathbb{R}^2} \frac{1}{2} (\varepsilon \cdot \mathbf{K} \cdot \varepsilon + \eta^2 (\nabla \varepsilon \cdot \mathbf{B} \cdot \nabla \varepsilon + 2 \varepsilon \cdot \mathbf{C} \cdot \nabla^2 \varepsilon)) \, dx. \quad (2.26) $$

This truncation amounts to trim off the blocks appearing below the SW-NE diagonal in the expression of t in (2.22), that is to replace t with

$$ t^{\dagger} = \begin{pmatrix} K & 0 & C \\ 0 & B & 0 \\ C^T & 0 & 0 \end{pmatrix} \quad (2.27) $$

Next the C-terms, proportional to ∇²ε, are integrated by parts and merged with the B-terms. This yields another approximation Φ†[ε] of the energy consistent with the previous one up to boundary terms (which we ignore since the laminate is infinite in all directions),

$$ \begin{aligned} \Phi^{\dagger}[\varepsilon] &= \int_{\mathbb{R}^2} \left( \frac{1}{2} \varepsilon(\mathbf{x}) \cdot \mathbf{K} \cdot \varepsilon(\mathbf{x}) + \frac{\eta^2}{2} \nabla \varepsilon(\mathbf{x}) \cdot \mathbf{B}^{\dagger} \cdot \nabla \varepsilon(\mathbf{x}) \right) \, dx \\ &= \int_{\mathbb{R}^2} \frac{1}{2} \begin{pmatrix} \varepsilon(\mathbf{x}) \\ \eta \nabla \varepsilon(\mathbf{x}) \\ \eta^2 \nabla^2 \varepsilon(\mathbf{x}) \end{pmatrix} \cdot \begin{pmatrix} \mathbf{K} & 0 & 0 \\ 0 & \mathbf{B}^{\dagger} & 0 \\ 0 & 0 & 0 \end{pmatrix} \cdot \begin{pmatrix} \varepsilon(\mathbf{x}) \\ \eta \nabla \varepsilon(\mathbf{x}) \\ \eta^2 \nabla^2 \varepsilon(\mathbf{x}) \end{pmatrix} \, dx. \end{aligned} \quad (2.28) $$

where the gradient properties are now captured by the symmetric matrix

$$ \mathbf{B}^{\dagger} = \mathbf{B} - (\tilde{\mathbf{C}} + \tilde{\mathbf{C}}^{\mathrm{T}}). \quad (2.29) $$

Here, $\tilde{\mathbf{C}}$ is the 4 × 4 matrix obtained by splitting the matrix $\mathbf{C} = (c_1 \ c_2)$ of dimensions 2 × 8 into two blocks $\mathbf{C}_1$ and $\mathbf{C}_2$ having dimensions 2 × 4 each, and rearranging them as

$$ \tilde{\mathbf{C}} = \begin{pmatrix} C_1 \\ C_2 \end{pmatrix}. $$

The new tensor $\mathbf{B}^\dagger$ capturing the gradient effect is not positive. Indeed, its diagonal component $B_{44}^\dagger$, calculated with the help of Equation [2.89] in [Le15], is strictly negative,

$$ B_{44}^{\dagger} = -\frac{a^2 (1-a)^2 (\mu_I - \mu_{II})^2 ((1-a) \mu_I + a \mu_{II})}{12 \mu_I \mu_{II}} < 0. \quad (2.30) $$

With a non-positive matrix of high-order homogenized moduli $\mathbf{B}^\dagger$, the energy Φ†[ε] is entirely useless. Indeed, observe that the mode $u_-(x_1, x_2) = \sin(\eta^{-1}x_2)$ has strain $\varepsilon_-= (0 \ \eta^{-1} \cos(\eta^{-1}x_2))$, strain gradient $\nabla\varepsilon_-= -(0 \ 0 \ 0 \ \eta^{-2} \sin(\eta^{-1}x_2))$, and contributes to the energy density through the gradient term as $\frac{\eta^2}{2} B_{44}^\dagger (\nabla\varepsilon_-)^2 = \frac{\eta^{-2}}{2} B_{44}^\dagger \sin^2(\eta^{-1}x_2)$, which on average over $x_2$ and for small η is a very large, negative contribution. An arbitrary equilibrium solution u that makes the energy Φ†[ε] stationary can be combined with a perturbation along the mode $u_-$ with a small amplitude $O(\eta^{1/2}) \ll 1$ to make the energy as large and negative as one wishes. This shows that the energy Φ†[ε] is not bounded from below in any neighborhood of any particular solution. In particular, no equilibrium is stable.

In the remainder of Section 2, we explore an alternative truncation method that preserves the positivity of t and requires no integration by parts.

## 2.5 Truncation limited to strain

In [SC00], it was observed that a positive energy is obtained when all terms in the right-hand side of (2.19) are kept and the error term O(η³) is discarded,

$$ \int_{\mathbb{R}^2} \frac{1}{2} (\varepsilon \cdot \mathbf{K} \cdot \varepsilon + \eta^2 (\nabla \varepsilon \cdot \mathbf{B} \cdot \nabla \varepsilon + 2 \varepsilon \cdot \mathbf{C} \cdot \nabla^2 \varepsilon) + \eta^4 \nabla^2 \varepsilon \cdot \mathbf{Z} \cdot \nabla^2 \varepsilon) \, dx. \quad (2.31) $$

The positivity of the right-hand side follows from (2.21).

In this approach, the strain (2.16) is truncated beyond order η² and inserted in the microscopic energy. The energy is not truncated further. As a result, the terms η⁴∇²ε·Z·∇²ε = O(η⁴) survive in the energy, warranting its positivity. In the original method of [SC00], the tensor Z is calculated exactly from the relocizer Ψ₂, which requires homogenization beyond second order, see Remark 2.1. This complicates the derivation of the model significantly. Recently, a variant of this method that neglects the Ψ₂ terms in Z has been proposed [KA24]. Either way, the ‘negligible’ Z block needs to be calculated, which is tedious, as shown by the (incomplete) set of expressions given in the Appendix of [KA24].
## 2.6 Minimal truncation preserving positivity

We proceed to present the new truncation method which preserves positivity without computing the Z-block. It only uses the asymptotically meaningful blocks K, B and C as input, i.e., it starts from the truncated matrix of elastic moduli,

$$t = \begin{pmatrix} K & 0_{2 \times 4} & C \\ 0_{4 \times 2} & B & \dots \\ C^T & \dots & \dots \end{pmatrix}. \qquad (2.32)$$

where the dots denote blocks associated with higher-order effects.

Instead of integrating the $\epsilon \otimes \nabla^2 \epsilon$ terms by parts, we incorporate them into the leading-order terms: up to terms of order $\eta^3$ denoted by the ellipses, we have

$$ \begin{aligned} \Phi^*[\epsilon] &= \int_{\mathbb{R}^2} \frac{1}{2} (\epsilon \cdot K \cdot \epsilon + \eta^2 (\nabla \epsilon \cdot B \cdot \nabla \epsilon + 2 \epsilon \cdot C \cdot \nabla^2 \epsilon)) \, d\mathbf{x} + \dots \\ &= \int_{\mathbb{R}^2} \frac{1}{2} (\{\epsilon + \eta^2 K^{-1} \cdot C \cdot \nabla^2 \epsilon\} \cdot K \cdot \{\epsilon + \eta^2 K^{-1} \cdot C \cdot \nabla^2 \epsilon\} + \eta^2 \nabla \epsilon \cdot B \cdot \nabla \epsilon) \, d\mathbf{x} + \dots \end{aligned} \qquad (2.33) $$

The right-hand sides on the first and second line differ only by a negligible term, $\eta^4 (C \cdot \nabla^2 \epsilon) \cdot K^{-1} \cdot (C \cdot \nabla^2 \epsilon)$. We have just identified another asymptotically correct approximation of the energy at this order, $\Phi_+[\epsilon] = \Phi^*[\epsilon] + \mathcal{O}(\eta^3)$, in the form

$$ \Phi_+[\epsilon] = \int_{\mathbb{R}^2} \left( \frac{1}{2} \{ \epsilon + \eta^2 \mathbf{K}^{-1} \cdot \mathbf{C} \cdot \nabla^2 \epsilon \} \cdot \mathbf{K} \cdot \{ \epsilon + \eta^2 \mathbf{K}^{-1} \cdot \mathbf{C} \cdot \nabla^2 \epsilon \} \right) \, d\mathbf{x} + \frac{\eta^2}{2} \nabla \epsilon \cdot \mathbf{B} \cdot \nabla \epsilon \qquad (2.34) $$

In addition to being asymptotically correct, this energy is positive thanks to the positivity of the symmetric matrices $K > 0$ and $B \ge 0$, see (2.25). The strict positivity of $K$ also ensures that it can be inverted.

The energy density in (2.34) is made up of two terms:

* On the first line, the leading-order homogenized elastic coefficients are combined with a generalized strain appearing in curly braces, which is the classical strain $\epsilon$ at leading order corrected by a second-gradient term $\eta^2 \mathbf{K}^{-1} \cdot \mathbf{C} \cdot \nabla^2 \epsilon$;

* On the second line, the positive matrix of higher-order moduli $B$ is combined with the standard strain gradient.

The key property of $\Phi_+$ in (2.34) is that the $C$ block, which has no special positivity property—it is not even a square matrix—, enters in the definition of the generalized strain (inside the arguments of the quadratic forms) but not in the definition of the elastic moduli (in the quadratic forms themselves). By contrast, the writing in (2.26) suggests that the entries in $C$ can be interpreted as homogenized elastic moduli, which is misleading.

The homogenized energy $\Phi_+[\epsilon]$ is therefore both asymptotically correct and simple to derive.

## 2.7 Notation recap

In Table 2.1, we recapitulate the symbols used for the energy at the successive steps of homogenization. The microscopic energy we started from is denoted as $\Phi[E]$, where $E$ is the microscopic strain. Using homogenization, we have been able to eliminate the macroscopic variables in favor of the macroscopic strain $\epsilon$ and derive a functional $\Phi^*[\epsilon]$. Next, we have discussed two strategies for truncating $\Phi^*$. In existing work, all terms that are negligible are removed, yielding a functional $\Phi^\dagger[\epsilon]$ which is subsequently integrated by parts into $\Phi^\dagger[\epsilon]$. The latter, however, is non-positive. We proposed an alternative truncation yielding a non-negative energy functional $\Phi_+ \ge 0$. All these energy functional agree up to an error of order $\eta^{p+1}$, where $\eta \ll 1$ is the scale separation parameter and $p$ is the homogenization order—we use $p=2$ throughout in this paper.

<table><caption>Table 2.1. Main steps in homogenization procedure and symbols used for the different types of energy functionals.</caption><tr><td>microscopic model</td><td>$\Phi[E]$</td></tr><tr><td>output of homogenization</td><td>$\Phi^*[\epsilon]$</td></tr><tr><td>naïve truncation</td><td>$\Phi^\dagger[\epsilon], \Phi^\ddagger[\epsilon]$</td></tr><tr><td>positive truncation</td><td>$\Phi_+[\epsilon]$</td></tr></table>

## 3 A simple lattice: the triangular truss

In this section, we show that the truncation procedure proposed in Section 2.6 can be applied to discrete elastic structures as well. We start with the triangular truss in dimension 2, shown in Figure 3.1. In a simple lattice such as the triangular truss, the nodal displacement is expressed directly in terms of the macroscopic displacement and there is no need to introduce a microscopic shift. This makes homogenization straightforward. Non-positivity issues do arise in the homogenized energy of simple lattices, however, and they can be fixed by the same strategy as that presented in Section 2.6 for the laminate, as we show. The case of a complex lattice is addressed in the next section.

This section and the next one are based on the high-order homogenization method for elastic lattices presented in [YAL24]. For a detailed presentation, the reader is referred to this paper.
<figure><img src="image_3.png" /><figcaption>**Figure 3.1.** A 2D infinite periodic triangular lattice, composed of identical springs with spring constant *k*. A unit cell is shown in red, and contains three springs of types $\varphi \in \{I,II,III\}$. The center of a particular spring analyzed in the text is denoted as $\mathbf{x}^c$ and its endpoints as $\pm$.</figcaption></figure>

## 3.1 Geometry and elastic energy

The triangular lattice is shown in Figure 3.1. It is obtained by repeating periodically the unit-cell, shown in red in the figure, along the generating vectors $\eta \mathbf{i}_1$ and $\eta \mathbf{i}'_2 = \eta (\mathbf{i}_1 - \sqrt{3} \mathbf{i}_2)/2$, where $\eta \ll 1$ is the length of the truss elements. There are three families of springs with different orientations in space, having all the same rest length $\eta$ and spring constant $k$. The families are labelled by an index $\varphi \in \{I, II, III\}$, with $\varphi = I$ for the springs aligned with the $\mathbf{i}_1$ direction.

In this simple lattice, homogenization is carried out by writing the nodal displacement $\mathbf{v}_\beta$ of a node labelled $\beta$ having undeformed position $\mathbf{x}_\beta$ in terms of a macroscopic displacement $\mathbf{u}: \mathbb{R}^2 \to \mathbb{R}^2$ as

$$ \mathbf{v}_\beta = \mathbf{u}(\mathbf{x}_\beta). \qquad (3.1) $$

Next, the linearized change of length of the different springs is calculated. Consider the particular spring of type $I$ having midpoint (center) $\mathbf{x}^c$ shown in Figure 3.1, for instance. We denote its endpoints as $\beta = -$ and $\beta = +$, having undeformed position $\mathbf{X}_\pm = \mathbf{x}^c \pm \frac{\eta}{2} \mathbf{i}_1$, respectively. The linearized change of length of the spring can be written as

$$ \begin{aligned} (\mathbf{v}_+ - \mathbf{v}_-) \cdot \mathbf{i}_1 &= \left( \mathbf{u}\left(\mathbf{x}^c + \frac{\eta}{2} \mathbf{i}_1\right) - \mathbf{u}\left(\mathbf{x}^c - \frac{\eta}{2} \mathbf{i}_1\right) \right) \cdot \mathbf{i}_1 \\ &= \eta \left( \frac{\partial u_1}{\partial x_1}(\mathbf{x}^c) + \frac{\eta^2}{3! 2^2} \frac{\partial^3 u_1}{\partial x_1^3}(\mathbf{x}^c) + \mathcal{O}(\eta^4) \right) \\ &= \eta \left( \varepsilon_{11}(\mathbf{x}^c) + \frac{\eta^2}{3! 2^2} \frac{\partial^2 \varepsilon_{11}}{\partial x_1^2}(\mathbf{x}^c) + \mathcal{O}(\eta^4) \right) \\ &= \eta \left( [\mathbf{i}_1 \otimes \mathbf{i}_1] : \varepsilon(\mathbf{x}^c) + \eta^2 \left[ \frac{1}{3! 2^2} \mathbf{i}_1 \otimes \mathbf{i}_1 \otimes \mathbf{i}_1 \otimes \mathbf{i}_1 \right] : \nabla^2 \varepsilon(\mathbf{x}^c) + \mathcal{O}(\eta^4) \right), \end{aligned} \qquad (3.2) $$

where the symbols : and :: represent multiple contraction, see Appendix A for a proper definition, and $\varepsilon(\mathbf{x})$ denotes the macroscopic strain,

$$ \varepsilon(\mathbf{x}) = \frac{1}{2} (\nabla \mathbf{u}(\mathbf{x}) + \nabla \mathbf{u}^T(\mathbf{x})). \qquad (3.3) $$

The strain $\varepsilon(\mathbf{x})$ is now a symmetric, 2 × 2 matrix as we no longer use the compact, vector form of the strain used earlier in (2.12), that was relevant to antiplane elasticity. The successive strain gradients $\nabla^k \varepsilon$ are now tensors of order $k+2$ and dimensions $2 \times \dots \times 2$ and are no longer represented in compact, vector form either.

In view of (3.2), we introduce the continuous function

$$ E_I(\mathbf{x}) = [\mathbf{i}_1 \otimes \mathbf{i}_1] : \varepsilon(\mathbf{x}) + \eta^2 \left[ \frac{1}{3! 2^2} \mathbf{i}_1 \otimes \mathbf{i}_1 \otimes \mathbf{i}_1 \otimes \mathbf{i}_1 \right] : \nabla^2 \varepsilon(\mathbf{x}) + \mathcal{O}(\eta^4), \qquad (3.4) $$

which yields the change of length in a spring of type $I$ as $\eta E_I(\mathbf{x}^c)$, where $\mathbf{x}^c$ is the spring midpoint position in undeformed configuration. The dependence of $E_I(\mathbf{x})$ on $\eta$ is implicit in our notation.

The argument is repeated with the two other families of springs and the quantities $E_I(\mathbf{x})$, $E_{II}(\mathbf{x})$ and $E_{III}(\mathbf{x})$ are collected in a vector called microscopic strain,

$$ \mathbf{E}(\mathbf{x}) = \begin{pmatrix} E_I(\mathbf{x}) \\ E_{II}(\mathbf{x}) \\ E_{III}(\mathbf{x}) \end{pmatrix} \in \mathbb{R}^3. \qquad (3.5) $$

Using (3.2) and similar formulas for the other types of springs, we can express $\mathbf{E}(\mathbf{x})$ in terms of the macroscopic strain $\varepsilon$ as

$$ \mathbf{E}(\mathbf{x}) = \mathbf{F}_0 : \varepsilon(\mathbf{x}) + \eta^2 \mathbf{F}_2 : \nabla^2 \varepsilon(\mathbf{x}) + \mathcal{O}(\eta^3). \qquad (3.6) $$

where the constant tensors $\mathbf{F}_0$ and $\mathbf{F}_2$, with respective dimensions 3 × 2 × 2 and 3 × 2 × 2 × 2 × 2 are obtained by assembling the coefficients appearing in square brackets in Equation (3.2) for type-I springs and in the corresponding equations for the other types of springs.
Denoting as $i \in \varphi$ the index of a generic spring of type $\varphi$ and as $\mathbf{x}_i^c$ its midpoint, we can write the energy of the discrete lattice as a sum of all springs, grouped by family,

$$ \Phi = \sum_{\varphi \in \{I,II,III\}} \sum_{i \in \varphi} \frac{1}{2} k (\eta E_{\varphi}(\mathbf{x}_i^c))^2. \qquad (3.7) $$

Using the Euler–macLaurin formula [KC02], and up to boundary terms which we ignore as the domain is infinite, we can rewrite each sum (involving a particular type of links $\varphi$) as an integral

$$ \begin{aligned} \Phi &= \sum_{\varphi \in \{I,II,III\}} \int_{\mathbb{R}^2} \frac{1}{2} \rho \eta^2 k E_{\varphi}^2(\mathbf{x}) \, d\mathbf{x} \\ &= \int_{\mathbb{R}^2} \frac{1}{2} \rho \eta^2 k \|\mathbf{E}(\mathbf{x})\|^2 d\mathbf{x}. \end{aligned} \qquad (3.8) $$

Here, $\rho$ represents the density per unit area of any particular type of springs,

$$ \rho = \frac{2}{\sqrt{3} \eta^2}. \qquad (3.9) $$

The strain energy (3.8) can be rewritten in a way similar to (2.6–2.7) as

$$ \Phi[\mathbf{E}] = \int_{\mathbb{R}^2} \frac{1}{2} \mathbf{E}(\mathbf{x}) \cdot \mathbf{k} \cdot \mathbf{E}(\mathbf{x}) \, d\mathbf{x}, \qquad (3.10) $$

where $\mathbf{k}$ is the stiffness matrix

$$ \mathbf{k} = \frac{2k}{\sqrt{3}} \mathbf{I}_3. \qquad (3.11) $$

Equation (3.10) follows from Equations [3.38,3.39,5.1] in [YAL24].

## 3.2 Series solution for the microscopic strain

As the lattice is simple, there are neither a microscopic shift nor a fast variable $y$, and the expansion (3.6) already yields the microscopic strain $\mathbf{E}$ in terms of the macroscopic strain $\varepsilon$ and its gradients, like Equation (2.16) for the laminate.

In the companion notebook homogenize-triangular-truss nb provided as supplementary material [TAL24b] to the present paper, we use the shoal library [Aud23], that carries out automatic second-order homogenization of elastic lattices based on the symbolic calculation language Wolfram Mathematica [Wol21]. The library is presented in detail in [AL23]. The geometrical and elastic features of the lattice are given to the library, which calculates the various quantities related to its homogenization, including the strain localization tensors $\mathbf{F}_0$ and $\mathbf{F}_2^*$: their detailed expressions are available in the companion Mathematica notebook.

**Remark 3.1.** The shoal library carries out homogenization to second order, and therefore calculates the displacement localizers $\mathbf{\Psi}_0$ and $\mathbf{\Psi}_1$ but not $\mathbf{\Psi}_2$. It returns the localizer $\mathbf{F}_2^*$ that misses the $\mathbf{\Psi}_2^*$ contribution, see (2.17)$_3$. This has no consequence by Remark 2.1.

## 3.3 Series expansion of the energy

By inserting the expansion of the strain $\mathbf{E}(\mathbf{x})$ in (3.6) into the continuous energy $\Phi[\mathbf{E}]$ in (3.10), we can define an effective energy $\Phi^*[\varepsilon]$ for the triangular lattice as

$$ \Phi^*[\varepsilon] = \int_{\mathbb{R}^2} \frac{1}{2} \varepsilon(\mathbf{x}) : \mathbf{K} : \varepsilon(\mathbf{x}) + \eta^2 \varepsilon(\mathbf{x}) : \mathbf{C} : \nabla^2 \varepsilon(\mathbf{x}) \, d\mathbf{x} + \mathcal{O}(\eta^3) \qquad (3.12) $$

where $\mathbf{K}$ and $\mathbf{C}$ are the fourth and sixth-tensors, respectively, given by

$$ \begin{aligned} K_{ijkl} &= (F_0)_{pij} k_{pq} (F_0)_{qkl} \\ C_{ijklmn} &= (F_0)_{pij} k_{pq} (F_2^*)_{qklmn}. \end{aligned} \qquad (3.13) $$

Here, we follow Einstein's summation convention, whereby a repeated index implies an implicit sum on this index. These expressions of $\mathbf{K}$ and $\mathbf{C}$ are similar to those derived for the laminate in (2.24). The $\mathbf{C}$ tensor uses only the contribution $\mathbf{F}_2^*$ to the localization tensor $\mathbf{F}_2$, which by Remark 2.1 does not change the result.

The tensors $\mathbf{K}$ and $\mathbf{C}$ are directly output by the shoal library in symbolic form as the library already implements the formulas (3.13). With the help of the symbolic calculation language Wolfram Mathematica, we have rewritten their action on $\varepsilon$ and $\nabla^2\varepsilon$ in terms of invariants that reflect the symmetries of the triangular lattice. The details of the reduction are available from the companion notebook and the result is

$$ \mathbf{K} : \varepsilon(\mathbf{x}) = (\lambda + \mu) \mathbf{I}_2 \text{tr} \varepsilon + 2 \mu \varepsilon^D, \qquad (3.14) $$

where $\varepsilon^D = \varepsilon - \mathbf{I}_2 \text{tr} \varepsilon / 2$ denotes the deviatoric part of $\varepsilon$ in dimension 2 and $\lambda$ and $\mu$ are the Lamé parameters predicted by leading-order homogenization,

$$ \lambda = \mu = \frac{\sqrt{3} k}{4} > 0. \qquad (3.15) $$
The expression of $\mathbf{K}$ in (3.14) is consistent with classical homogenization results from previous work, see [VP12, S18c] among others.

Equation (3.14) shows that the quadratic form $\mathbf{K}$ acting on symmetric tensors $\boldsymbol{\varepsilon}$ is strictly positive,

$$ \mathbf{K} > 0. \qquad (3.16) $$

The tensor $\mathbf{C}$ in (3.13) is reduced similarly using a symbolic calculation language and we obtain

$$ \mathbf{C} :: \nabla^2 \boldsymbol{\varepsilon}(\mathbf{x}) = \frac{k}{32\sqrt{3}} \left( \nabla^2 \mathrm{tr} \boldsymbol{\varepsilon} + \Delta \boldsymbol{\varepsilon} + \frac{I_2}{2} \mathcal{D} :: \nabla^2 \boldsymbol{\varepsilon} + \left[ \frac{\mathcal{F} : \nabla^2 \boldsymbol{\varepsilon} : \mathcal{F}^{T[3,1,2]}}{2} \right]^{\mathrm{sym}} \right), \qquad (3.17) $$

where $\Delta \boldsymbol{\varepsilon} = \nabla^2 \boldsymbol{\varepsilon} : I_2$ denotes the Laplacian of the strain,

$$ \mathcal{F} = i_1 \otimes (-i_1 \otimes i_1 + i_2 \otimes i_2) + i_2 \otimes (i_1 \otimes i_2 + i_2 \otimes i_1) \qquad (3.18) $$

is a rank-3 tensor that is invariant by the $D_6$ symmetry of the triangular lattice, $\mathcal{D}_{ijkl} = \mathcal{F}_{aij} \mathcal{F}_{akl} / 2$ is the projector on the subspace of deviatoric strain, $\mathbf{s}^{T[3,1,2]}$ denotes the generalized transpose of a rank-3 tensor $\mathbf{s}$ such that $(\mathbf{s}^{T[3,1,2]})_{jki} = s_{jki}$, and $[\mathbf{i}]^{\mathrm{sym}} = (\mathbf{i} + \mathbf{i}^T) / 2$ denotes the symmetric part of a $2 \times 2$ tensor.

**Remark 3.2.** Comparing (2.23) and (3.13) reveals that $\mathbf{B}=\mathbf{0}$ for the triangular lattice. This is a direct consequence of the fact that $\mathbf{F}_1=\mathbf{0}$, as shown by comparing (2.16) and (3.6), see also the expression of $\mathbf{B}$ in (2.24).

**Remark 3.3.** Since for this application (as well as for the forthcoming ones) we will be applying the new truncation method presented in Section 2.6 and will be bypassing the truncation limited to strain from [SC00] (see Section 2.5), we have not even derived the $\mathbf{Z}$ terms of order $\eta^4$ in the energy (3.12), unlike earlier in (2.23).

## 3.4 Naïve truncation of the energy

If we truncate and integrate by parts the term $\boldsymbol{\varepsilon}(\mathbf{x}) : \mathbf{C} :: \nabla^2 \boldsymbol{\varepsilon}(\mathbf{x})$ as we did earlier in (2.28), and as is commonly done in the literature, we get an energy $\Phi^{\dagger}[\boldsymbol{\varepsilon}]$ in the form

$$ \Phi^{\dagger}[\boldsymbol{\varepsilon}] = \int_{\mathbb{R}^2} \frac{1}{2} \boldsymbol{\varepsilon}(\mathbf{x}) : \mathbf{K} : \boldsymbol{\varepsilon}(\mathbf{x}) - \frac{\eta^2}{2} \nabla \boldsymbol{\varepsilon}(\mathbf{x}) : (2\mathbf{C}) : \nabla \boldsymbol{\varepsilon}(\mathbf{x}) \, d\mathbf{x}. \qquad (3.19) $$

This expression is similar to that in (2.28–2.29) for the laminate, but with $\mathbf{B}^{\dagger}=-2\mathbf{C}$ since $\mathbf{B}=\mathbf{0}$. Note that the tensor $\mathbf{C}$ natively satisfies the major symmetries, $C_{ijklmn} = C_{lmnjk}$.

The quadratic form $\mathbf{B}^{\dagger}=-2\mathbf{C}$ is calculated in symbolic form in the companion notebook, and all its eigenvalues are found to be negative, $\mathbf{B}^{\dagger} \le 0$, i.e., $2\mathbf{C} \ge 0$. The energy $\Phi^{\dagger}[\boldsymbol{\varepsilon}]$ is therefore non-positive and the naïve truncation fails again.

## 3.5 Truncation preserving positivity

As argued earlier in Section 2.6, incorporating the second-order term $\boldsymbol{\varepsilon}(\mathbf{x}) : \mathbf{C} :: \nabla^2 \boldsymbol{\varepsilon}(\mathbf{x})$ in (3.12) into the leading-order term $\boldsymbol{\varepsilon}(\mathbf{x}) : \mathbf{K} : \boldsymbol{\varepsilon}(\mathbf{x})$ preserves positivity. This yields an approximation of the energy that is exact to order $\eta^2$ included, $\Phi_+[\boldsymbol{\varepsilon}] = \Phi^*[\boldsymbol{\varepsilon}] + \mathcal{O}(\eta^3)$, in the form

$$ \Phi_+[\boldsymbol{\varepsilon}] = \int_{\mathbb{R}^2} \frac{1}{2} \tilde{\boldsymbol{\varepsilon}} : \mathbf{K} : \tilde{\boldsymbol{\varepsilon}} \, d\mathbf{x} \quad \text{where } \tilde{\boldsymbol{\varepsilon}} = \{\boldsymbol{\varepsilon}(\mathbf{x}) + \eta^2 \mathbf{K}^{-1} : \mathbf{C} :: \nabla^2 \boldsymbol{\varepsilon}(\mathbf{x})\}. \qquad (3.20) $$

Here, $\mathbf{K}^{-1} = \frac{I_2 \otimes I_2}{4(\lambda+\mu)} + \frac{\mathcal{D}}{2\mu}$ is the inverse of the strictly positive quadratic form $\mathbf{K}$, satisfying $\mathbf{K}^{-1} : \mathbf{K} : \boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}$. The proposed energy $\Phi_+$ matches that in (3.12) up to a negligible term, of order $\eta^4$. It also closely resembles that of the laminate in (2.34), except for the fact that the triangular lattice has $\mathbf{B}=\mathbf{0}$. The positivity of leading-order stiffness $\mathbf{K}$ in (3.16) warrants that $\Phi_+[\boldsymbol{\varepsilon}]$ is positive.

The generalized strain $\tilde{\boldsymbol{\varepsilon}} = \{\boldsymbol{\varepsilon} + \eta^2 \mathbf{K}^{-1} : \mathbf{C} :: \nabla^2 \boldsymbol{\varepsilon}\}$ is obtained in the companion notebook in explicit form as

$$ \tilde{\boldsymbol{\varepsilon}} = \left\{ \boldsymbol{\varepsilon} + \frac{\eta^2}{48} \left( \Delta \boldsymbol{\varepsilon} + (\nabla^2 \mathrm{tr} \boldsymbol{\varepsilon})^D + \left[ \frac{\mathcal{F} : \nabla^2 \boldsymbol{\varepsilon} : \mathcal{F}^{T[3,1,2]}}{2} \right]^{\mathrm{sym}} \right) \right\}, \qquad (3.21) $$

where $(\nabla^2 \mathrm{tr} \boldsymbol{\varepsilon})^D = \mathcal{D} :: (\nabla^2 \mathrm{tr} \boldsymbol{\varepsilon}) = \nabla^2 \mathrm{tr} \boldsymbol{\varepsilon} - I_2 \Delta \mathrm{tr} \boldsymbol{\varepsilon} / 2$ denotes the deviatoric part of the second gradient of $(\mathrm{tr} \boldsymbol{\varepsilon})$. This, together with the expression of $\mathbf{K}$ in (3.14), makes the expression (3.20) of $\Phi_+$ fully explicit.

**Note 3.4.** In its first version, the library shoul was implementing the naïve truncation method from Section 3.4, which broke positivity. Along with this paper, we are releasing an update to the library that outputs the homogenized properties $\mathbf{K}$, $\mathbf{B}$ and $\mathbf{C}$ without integrating by parts. This update to the library was used to produce all the homogenization results for the lattices analyzed in this paper (triangular, honeycomb and pantograph).

**Remark 3.5. (geometric compatibility)** A symmetric tensor field $\tilde{\boldsymbol{\varepsilon}}(\mathbf{x})$ is said to be geometrically compatible if there exists a vector field $\tilde{\mathbf{u}}(\mathbf{x})$ such that $\tilde{\boldsymbol{\varepsilon}}(\mathbf{x}) = (\nabla \tilde{\mathbf{u}}(\mathbf{x}) + \nabla \tilde{\mathbf{u}}^T(\mathbf{x})) / 2$. In dimension 2, a necessary condition for geometrically compatibility is $\partial \tilde{\varepsilon}_{11} / \partial x_2^2 + \partial \tilde{\varepsilon}_{22} / \partial x_1^2 - 2 \partial \tilde{\varepsilon}_{12} / (\partial x_1 \partial x_2) = 0$ for all $\mathbf{x}$. In the companion notebook, we test this equality using the generalized strain $\tilde{\boldsymbol{\varepsilon}}$ given in (3.21) and find that it is violated at order $\eta^2$. This proves that there does not exist a generalized displacement $\tilde{\mathbf{u}}(\mathbf{x})$ in terms of which the generalized strain could be defined as $\tilde{\boldsymbol{\varepsilon}}(\mathbf{x}) = (\nabla \tilde{\mathbf{u}}(\mathbf{x}) + \nabla \tilde{\mathbf{u}}^T(\mathbf{x})) / 2$.
## 4 A complex lattice: the honeycomb

In this section, we consider the honeycomb lattice in dimension 2, as sketched in Figure 4.1. It is a complex lattice as it contains two different types of nodes (shown in black and in white in the figure) that are not congruent to one another via the generating translations of the lattice. In complex lattices, one needs to introduce microscopic degrees of freedom, denoted as $\xi$ and $\psi_{\pm}$ here, that capture the relative displacements and rotations of the two types of nodes. Homogenization links the microscopic degrees of freedom to the macroscopic strain $\varepsilon$ and its gradients.

Second-order homogenization of the honeycomb lattice has been carried out in [YAL24], and we follow closely this presentation in Sections 4.1 to 4.4. In Section 4.5, we apply the novel truncation scheme preserving positivity to the non-positive energy obtained in this previous work.

<figure><img src="image_4.png" /><figcaption>**Figure 4.1.** A 2D infinite periodic honeycomb lattice, composed of identical beams. The red cell, indexed $i$, is the minimal unit-cell and contains three beams of type $\varphi \in \{I, II, III\}$. The center of a particular beam that is analyzed in the text is denoted as $\mathbf{x}^c$.</figcaption></figure>

### 4.1 Geometry and elastic energy

We consider an infinite honeycomb lattice, obtained by repeating the unit-cell (shown in red in the figure) periodically along the generating vectors $\mathbf{i}_1' = (\sqrt{3}\mathbf{i}_1 + \mathbf{i}_2)/2$ and $\mathbf{i}_2$. The unit cell contains three families of beams of length $\eta$, each family being labelled $\varphi \in \{I, II, III\}$. The beams are considered inextensible.

When this complex lattice is homogenized, the Cauchy-Born rule states that the displacement $\mathbf{v}_{\beta}$ and rotation $\theta_{\beta}$ of a node having undeformed position $\mathbf{x}_{\beta}$ take different forms depending on the node family, see also Equation [3.4] in [YAL24],

$$ \begin{array}{c|cc} & \text{displacement } \mathbf{v}_{\beta} & \text{rotation } \theta_{\beta} \\ \text{white nodes} & \mathbf{u}(\mathbf{x}_{\beta}) - \eta \xi(\mathbf{x}_{\beta}) & \gamma(\mathbf{x}_{\beta}) + \psi_{-}(\mathbf{x}_{\beta}) \\ \text{black nodes} & \mathbf{u}(\mathbf{x}_{\beta}) + \eta \xi(\mathbf{x}_{\beta}) & \gamma(\mathbf{x}_{\beta}) + \psi_{+}(\mathbf{x}_{\beta}). \end{array} \qquad (4.1) $$

Here $\mathbf{u}(\mathbf{x})$, $\xi(\mathbf{x})$, $\gamma(\mathbf{x})$ and $\psi_{\pm}(\mathbf{x})$ are smooth functions of the macroscopic variable $\mathbf{x}$ that can be expressed as regular series in $\eta$, as earlier in (2.10), $\mathbf{u}(\mathbf{x})$ being the macroscopic displacement, $\gamma(\mathbf{x}) = (\partial u_2 / \partial x_1 - \partial u_1 / \partial x_2) / 2$ the associated rotation, $\xi(\mathbf{x})$ the microscopic shift and $\psi_{\pm}(\mathbf{x})$ the microscopic rotation of the black nodes ($\psi_{-}$) and white notes ($\psi_{+}$). Even though these functions are smooth, the displacement $\mathbf{v}_{\beta}$ and rotation $\theta_{\beta}$ in (4.1) vary quickly at the scale of a unit cell due to the alternating signs $\pm$ in front of the microscopic displacement ($\pm\eta\xi$ term) and in the microrotation ($\psi_{\pm}$ terms).

The main task of homogenization is to solve for the microscopic displacement and rotation in terms of the macroscopic strain $\varepsilon(\mathbf{x}) = (\nabla\mathbf{u} + \nabla\mathbf{u}^T)/2$, see (3.3). A perturbative solution is available from Equation [5.10] in [YAL24],

$$ \begin{aligned} \xi(\mathbf{x}) &= -\frac{1}{4} \mathcal{J} : \varepsilon(\mathbf{x}) + \mathbf{0} \times \eta + \dots \\ \psi_{\pm}(\mathbf{x}) &= \mathbf{0} \pm \frac{3\eta}{4} \mathrm{curl}(\mathcal{J} : \varepsilon(\mathbf{x})) + \dots \end{aligned} \qquad (4.2) $$

where $\mathcal{J}$ is the rank-3 tensor characteristic of $D_6$ symmetry introduced in (3.18) and $\mathrm{curl}\,\mathbf{a}(\mathbf{x})$ denotes the (scalar) curl of a vector field in dimension 2.

The inextensibility of the beams requires the macroscopic deformation to be incompressible,

$$ \mathrm{tr}\,\varepsilon = 0, \qquad (4.3) $$

as stated for instance in Equation [5.8] from [YAL24].

The energy of the beam labelled 'I' in Figure 4.1 is given by linear beam theory as $\frac{1}{2}\frac{EI}{\eta}(\kappa_I^2 + 12\tau_I^2)$ where $\eta$ is the beam length, $EI$ is the bending modulus, and $\kappa_I = \theta_+ - \theta_-$ and $\tau_I = \frac{\theta_+ + \theta_-}{2} - \frac{(v_+ - v_-) \cdot \mathbf{i}_2}{\eta}$ are two strain measures relevant to the linear theory of inextensible beams, see Equation [2.12] in [YAL24]. The elastic energy for the entire lattice can then be written as

$$ \Phi = \sum_{\varphi \in \{I, II, III\}} \sum_{i \in \varphi} \frac{1}{2} \frac{EI}{\eta} (\kappa_{\varphi}^2 + 12 \tau_{\varphi}^2), \qquad (4.4) $$
where the first sum runs over beam families $\varphi$, and the second sum over beams belonging to a particular family, see Equations [2.16] and [2.23] in [YAL24]. The dependence of the strain measures $\kappa_\varphi$ and $\tau_\varphi$ on the particular member $i$ of the family $\varphi$ is implicit in the sum above.

## 4.2 Series solution for the microscopic strain

The strain measure $\kappa_I = \theta_+ - \theta_-$ quantifies the mean curvature imposed by the rotation $\theta_\pm$ of the endpoints. By a similar calculation as earlier in (3.2), but paying attention to the fact that the (-) endpoint is now black whereas the (+) endpoint is white, one can express the strain $\kappa_I$ in the horizontal beam labelled 'I' in Figure 4.1, having midpoint $\mathbf{x}^c$, as

$$ 
\begin{aligned}
\kappa_I &= \theta_+ - \theta_- \\
&= \left( \gamma(\mathbf{x}^c + \frac{\eta}{2}\mathbf{i}_1) + \psi_-(\mathbf{x}^c + \frac{\eta}{2}\mathbf{i}_1) \right) - \left( \gamma(\mathbf{x}^c - \frac{\eta}{2}\mathbf{i}_1) + \psi_+(\mathbf{x}^c - \frac{\eta}{2}\mathbf{i}_1) \right) \\
&= \eta \nabla\gamma(\mathbf{x}^c) + \dots + \psi_-(\mathbf{x}^c) + \frac{\eta}{2} \nabla\psi_-(\mathbf{x}^c) \cdot \mathbf{i}_1 + \dots - \psi_+(\mathbf{x}^c) + \frac{\eta}{2} \nabla\psi_+(\mathbf{x}^c) \cdot \mathbf{i}_1 + \dots
\end{aligned} 
\quad (4.5) $$

The gradient of the macroscopic rotation $\nabla\gamma$ can be worked out as $\nabla\gamma = (\partial\varepsilon_{12}/\partial x_1 - \partial\varepsilon_{11}/\partial x_2)\mathbf{i}_1 + (-\partial\varepsilon_{12}/\partial x_2 + \partial\varepsilon_{22}/\partial x_1)\mathbf{i}_2 = \mathcal{G}$. $\mathcal{G} \coloneqq \nabla\varepsilon$ in terms of the strain gradient $\nabla\varepsilon$, where $\mathcal{G} = \mathbf{i}_1 \otimes \frac{\mathbf{i}_1 \otimes \mathbf{i}_2 + \mathbf{i}_2 \otimes \mathbf{i}_1}{2} \otimes \mathbf{i}_1 - \mathbf{i}_1 \otimes \mathbf{i}_1 \otimes \mathbf{i}_1 \otimes \mathbf{i}_2 - \mathbf{i}_2 \otimes \frac{\mathbf{i}_1 \otimes \mathbf{i}_2 + \mathbf{i}_2 \otimes \mathbf{i}_1}{2} \otimes \mathbf{i}_2 + \mathbf{i}_2 \otimes \mathbf{i}_2 \otimes \mathbf{i}_2 \otimes \mathbf{i}_1$ is a constant rank-4 tensor, see Appendix A.3 of [YAL24]. On using the solution $(4.2)_2$ for the microscopic shift and rewriting curl ($\mathcal{G}' : \varepsilon = \mathcal{G}' \cdot \nabla\varepsilon$ with the help of another constant rank-4 tensor $\mathcal{G}'$ (whose expression is omitted for the sake of brevity), one obtains the Taylor expansion of $\kappa_I$ as

$$ \kappa_I = [\mathbf{0}] : \varepsilon(\mathbf{x}^c) + \eta \left[ \mathcal{G} - \frac{3}{2} \mathcal{G}' \right] : \nabla\varepsilon(\mathbf{x}^c) + \eta^2 [\mathbf{0}] : \nabla^2\varepsilon(\mathbf{x}^c) + O(\eta^3). \quad (4.6) $$

Next, we introduce the full set of deformation measures relevant to all 3 families of beams in the form the microscopic strain vector,

$$ \mathbf{E} = (\kappa_I \quad \tau_I \quad \kappa_{II} \quad \tau_{II} \quad \kappa_{III} \quad \tau_{III}) \in \mathbb{R}^6. \quad (4.7) $$

A similar calculation as that done above in (4.6) yields the Taylor expansions of the strain measures $\kappa_\varphi$ and $\tau_\varphi$ for each type $\varphi$ of beams. Collecting them, one obtains a Taylor expansion for $\mathbf{E} = \mathbf{E}(\mathbf{x})$ in the same form as earlier in (3.6), where $\mathbf{E}(\mathbf{x})$ is the smooth function

$$ \mathbf{E}(\mathbf{x}) = \mathbf{F}_0 : \varepsilon(\mathbf{x}) + \eta \mathbf{F}_1 : \nabla\varepsilon(\mathbf{x}) + \eta^2 \mathbf{F}_2 : \nabla^2\varepsilon(\mathbf{x}) + \mathcal{O}(\eta^3). \quad (4.8) $$

The tensors $\mathbf{F}_0$, $\mathbf{F}_1$ and $\mathbf{F}_2^*$ are calculated in symbolic form using the homogenization library `shoal` in the companion notebook `homogenize-inextensible-honeycomb nb`, provided as supplementary material [TAL24b]. The tensor $\mathbf{F}_2^*$ matches $\mathbf{F}_2$ up to terms depending on $\Psi_2^*$ that can be ignored, see Remarks 2.1 and 3.1. The tensors $\mathbf{F}_0$, $\mathbf{F}_1$ and $\mathbf{F}_2^*$ are constant that have no dependence on any parameter.

## 4.3 Series expansion of the energy

By a similar argument as earlier in Equation (3.8), the strain energy $\Phi$ in (4.4) can be continualized as

$$ \Phi[\mathbf{E}] = \frac{1}{2} \int_{\mathbb{R}^2} \mathbf{E}(\mathbf{x}) \cdot \mathbf{k} \cdot \mathbf{E}(\mathbf{x}) \, d\mathbf{x} \quad (4.9) $$

where $\mathbf{k}$ is the diagonal matrix

$$ \mathbf{k} = \frac{\rho EI}{\eta} \begin{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \\ \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \\ \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \end{pmatrix}, \quad (4.10) $$

with missing components being zero. The 3 diagonal blocks are each associated with a family of spring, to reflect the ordering conventions in (4.7). The coefficients 1 and 12 making up each block come from Equation (4.4) and $\rho = 2/(3\sqrt{3}\eta^2)$ is the area density of each particular type of beam. The continualized form of the energy (4.9) appears in Equations [3.38,3.39,5.1] in the original paper [YAL24].

For the lattice energy to remain finite in the limit $\eta \to 0$, we require that the bending moduli $EI$ scale as $\eta^3$, i.e., we rewrite $EI = \bar{E}\bar{I}\eta^3$ where the scaled bending modulus $\bar{E}\bar{I}$ is finite. The coefficient in the right-hand side of (4.10) can then be rewritten as $\rho EI/\eta = 2\bar{E}\bar{I}/(3\sqrt{3})$.

By inserting the approximated strain $\mathbf{E}(\mathbf{x})$ (4.8) into the continuous energy $\Phi[\mathbf{E}]$ (4.9), we identify an effective energy for the honeycomb lattice in the form

$$ \Phi^*[\varepsilon^D] = 2\mu \int_{\mathbb{R}^2} \frac{1}{2} \varepsilon^D(\mathbf{x}) : \varepsilon^D(\mathbf{x}) + \eta^2 \left( \frac{1}{2} \nabla\varepsilon^D(\mathbf{x}) : \mathbf{B} : \nabla\varepsilon^D(\mathbf{x}) + \varepsilon^D(\mathbf{x}) : \mathbf{C} : \nabla^2\varepsilon^D(\mathbf{x}) \right) d\mathbf{x} + \mathcal{O}(\eta^3) \quad (4.11) $$

where $\mu$ is the homogenized shear modulus predicted by classical homogenization

$$ \mu = 4\sqrt{3}\bar{E}\bar{I}, \quad (4.12) $$

which is obviously positive, $\mu > 0$. By the incompressibility condition (4.3), the macroscopic strain is purely deviatoric, hence the notation $\varepsilon = \varepsilon^D$ in the homogenized energy. There is no term coupling $\varepsilon^D$ and $\nabla\varepsilon^D$ in the expression (4.11) of $\Phi^*$, as a consequence of the lattice being centro-symmetric.
The tensors **B** and **C** are calculated by the library *shoal* in terms of the $F_i$'s using formulas similar to (2.24) and (3.13). They are given in symbolic form in the companion notebook.

With the help of symbolic calculations, we have further reduced the operators **B** and **C** in the companion notebook: their action on the first and second gradient of strain has been rewritten in terms of invariants as

$$ 
\begin{aligned}
\nabla \varepsilon^D(\mathbf{x}) \circled{=} \mathbf{B} \circled{=} \nabla \varepsilon^D(\mathbf{x}) &= \frac{3}{16} \mathrm{curl}^2(\mathcal{F}:\varepsilon^D) + \frac{1}{48} \|\mathcal{D}:\nabla(\mathcal{F}:\varepsilon^D)\|^2 \\
\mathbf{C} \circled{=} \nabla^2 \varepsilon^D(\mathbf{x}) &= \frac{1}{16} \left( 3 \Delta \varepsilon^D - 7 \mathcal{D}: \left[ \frac{\mathcal{F}:\nabla^2 \varepsilon^D: \mathcal{F}^T [3,1,2]}{2} \right] \right),
\end{aligned}
\quad (4.13) $$

where $\mathcal{F}$ is given in (3.18), $\mathcal{D}$ and $T[3,1,2]$ are the deviatoric projector and generalized transpose already appearing in (3.17), respectively, $\Delta$ is the Laplacian and the square norm $\|\mathbf{c}\|^2$ of a (deviatoric) tensor $\mathbf{c}$ is defined as $\|\mathbf{c}\|^2 = \mathbf{c}:\mathbf{c}$. The positivity of the quadratic form $\mathbf{B}$ is evident from (4.13)$_1$,

$$ \nabla \varepsilon^D(\mathbf{x}) \circled{=} \mathbf{B} \circled{=} \nabla \varepsilon^D(\mathbf{x}) \geq 0. \quad (4.14) $$

## 4.4 Naïve truncation of the energy

The original analysis of honeycomb lattices [YAL24] was based on an early version of the library *shoal* which calculated the expression (4.11) internally but carried out an integration by parts and the naïve truncation without exposing the form $\Phi^*$ nor the tensors **B** and **C** to the user: the homogenized energy reported in Equation [5.9] of this paper [YAL24] took the form (analogous to (3.19))

$$ \Phi^{\dagger}[\varepsilon^D] = 2\mu \int_{-\infty}^{+\infty} \left( \frac{1}{2} \varepsilon^D : \varepsilon^D + \frac{\eta^2}{2 \times 8} (-7 \nabla \varepsilon^D \circled{=} \nabla \varepsilon^D + \frac{22}{3} \|\nabla \varepsilon^D : \mathbf{I}_2\|^2 + 2 \mathrm{div}^2(\mathcal{F}:\varepsilon^D)) \right) \mathrm{d}\mathbf{x}. \quad (4.15) $$

The strain-gradient contribution is non-positive due to the negative coefficient (–7), as noted by the authors in their Remark 17.

## 4.5 Truncation preserving positivity

Applying instead the positivity-preserving truncation scheme from Section 2.6 to (4.11), we obtain an alternative expression of the homogenized energy that is consistent to order $\eta^2$ included,

$$ \Phi_+[\varepsilon^D] = 2\mu \int_{\mathbb{R}^2} \left( \frac{1}{2} \{\varepsilon^D(\mathbf{x}) + \eta^2 \mathbf{C} \circled{=} \nabla^2 \varepsilon^D(\mathbf{x})\} : \{\varepsilon^D(\mathbf{x}) + \eta^2 \mathbf{C} \circled{=} \nabla^2 \varepsilon^D(\mathbf{x})\} \right. \\ \left. + \frac{\eta^2}{2} \nabla \varepsilon^D(\mathbf{x}) \circled{=} \mathbf{B} \circled{=} \nabla \varepsilon^D(\mathbf{x}) \right) \mathrm{d}\mathbf{x}. \quad (4.16) $$

The quantity $\{\dots\}$ inside the curly braces is a generalized strain measure. Equation (4.13)$_2$ shows that the strain correction $\eta^2 \mathbf{C} \circled{=} \nabla^2 \varepsilon^D(\mathbf{x})$ is a deviatoric rank-2 tensor (i.e., it is symmetric and traceless) like the leading-order strain $\varepsilon^D$.

A fully explicit expression of $\Phi_+$ is obtained by combining the expression above with (4.13), which yields

$$ \Phi_+[\varepsilon^D] = 2\mu \int_{\mathbb{R}^2} \left( \frac{1}{2} \|\varepsilon^D + \frac{\eta^2}{16} (3 \Delta \varepsilon^D - 7 \mathcal{D}: \left[ \frac{\mathcal{F}:\nabla^2 \varepsilon^D: \mathcal{F}^T [3,1,2]}{2} \right])\|^2 \right. \\ \left. + \frac{\eta^2}{2} \left( \frac{3}{16} \mathrm{curl}^2(\mathcal{F}:\varepsilon^D) + \frac{1}{48} \|\mathcal{D}:\nabla(\mathcal{F}:\varepsilon^D)\|^2 \right) \right) \mathrm{d}\mathbf{x}, \quad (4.17) $$

This functional is clearly positive.

**Remark 4.1. (geometric compatibility)** By the same method as for the triangular lattice, see Remark 3.5 at the end of section 3, it has been checked in the companion notebook that the generalized strain $\{\varepsilon^D(\mathbf{x}) + \eta^2 \mathbf{C} \circled{=} \nabla^2 \varepsilon^D(\mathbf{x})\}$ is not geometrically compatible.

# 5 Connecting with Cholesky decomposition

In this section, we connect with the LDLT Cholesky decomposition, a classical method used in numerical analysis [Str07]. The Cholesky decomposition takes a symmetric positive matrix $\mathbf{t}$ on input and decomposes it iteratively into $\mathbf{t} = \mathbf{L} \cdot \mathbf{D} \cdot \mathbf{L}^T$, where $\mathbf{L}$ is a lower-triangular matrix filled with 1's on the diagonal, and $\mathbf{D}$ is a diagonal matrix containing only positive entries, $\mathbf{D} \ge 0$.

We introduce the following variant of the Cholesky decomposition:

* We start from a symbolic matrix $\mathbf{t}$ whereas the original Cholesky method works on numerical matrices.

* We carry out only the first steps of the Cholesky decomposition: when homogenizing to order $\eta^p$ included ($p=2$ in our examples), we stop after $p+1$ iterations.

* In our variant, the input matrix $\mathbf{t}$ is defined by *blocks* and the Cholesky decomposition is carried out at the block level and not at the level of the individual entries of $\mathbf{t}$. This enables us to handle the common situation where there are multiple scalar measures of macroscopic strain—the laminate example from Section 2, for instance, uses a strain $\mathbf{E}$ having 2 components, see (2.4).
The Cholesky decomposition helps extending the truncation method proposed so far in two important ways:

* to non-centro-symmetric structures (the examples analyzed so far were all centro-symmetric, a particular case in which our truncation method simply requires moving the cross terms of the form $C\epsilon\epsilon''$ into the leading order $K\epsilon^2/2$, see for instance Section 2.6).

* to arbitrary homogenization order (even though emphasis is on second-order homogenization here).

## 5.1 Cholesky decomposition by blocks

The input to our variant of the Cholesky decomposition is symmetric matrix $t$ written in block-matrix notation as

$$t = \begin{pmatrix} K & \underline{A} & \underline{C} & \cdots \\ \underline{A}^T & \underline{B} & \cdots \\ \underline{C}^T & \underline{\underline{\underline{C}}} & \cdots \\ \vdots & \vdots & \ddots & \ddots \end{pmatrix}, \qquad (5.1)$$

where the ellipses represents blocks that we do not need to resolve as they represent higher-order effects. Obviously, the diagonal blocks $K, \underline{B}$, etc. have to be symmetric.

Anticipating on what comes next, we annotate the various blocks with multiple underlines to denote the iteration at which the various blocks enter into the procedure: $K$ will be used in the first iteration (numbered 0), $\underline{A}$ in the second iteration (numbered 1) and $\underline{B}$ and $\underline{C}$ in the third iteration (numbered 2). A key remark is that the iteration numbers match the orders of the energy terms associated with the different blocks when $t$ is interpreted as the matrix of elastic coefficients—see Equations (2.19) and (2.22) for example. Indeed, the $K$ block gives rise to an energy contribution $\epsilon \cdot K \cdot \epsilon / 2$ of order $\eta^0$, the $\underline{A}$ block gives rise to an energy contribution $\epsilon \cdot \underline{A} \cdot (\eta \epsilon')$ of order $\eta^1$, the $\underline{B}$ and $\underline{C}$ blocks give rise to energy contributions $(\eta \epsilon') \cdot \underline{B} \cdot (\eta \epsilon') / 2 + \epsilon \cdot \underline{C} \cdot (\eta^2 \epsilon'')$ of order $\eta^2$, etc. This match in the orders of magnitude of the blocks in the matrix $t$ produced by homogenization, and in the iteration number of the Cholesky decomposition makes the latter perfectly suited to the former.

We assume that the top-left (symmetric) $K$ block is positive definite, and that the matrix $t$ is positive but not necessarily definite. The only positivity properties which we will use in second-order homogenization are

$$K > 0, \quad \begin{pmatrix} K & \underline{A} \\ \underline{A}^T & \underline{B} \end{pmatrix} \ge 0. \qquad (5.2)$$

For homogenization at order $p$, we would assume that the first $p$ square super-blocks starting from the top-left corner of $t$ are all positive.

Following Cholesky, we seek a decomposition of $t$ in the form

$$t = \begin{pmatrix} \underline{I} & \underline{0} & \underline{0} & \cdots \\ \underline{\underline{L}}_1 & \underline{I} & \underline{0} & \cdots \\ \underline{\underline{\underline{L}}}_2 & \cdots & \cdots & \cdots \\ \vdots & \vdots & \ddots & \ddots \end{pmatrix} \cdot \begin{pmatrix} \underline{D}_0 & \underline{0} & \underline{0} & \cdots \\ \underline{0} & \underline{\underline{D}}_1 & \cdots & \cdots \\ \underline{0} & \cdots & \cdots & \cdots \\ \vdots & \vdots & \ddots & \ddots \end{pmatrix} \cdot \begin{pmatrix} \underline{I} & \underline{0} & \underline{0} & \cdots \\ \underline{\underline{L}}_1 & \underline{I} & \underline{0} & \cdots \\ \underline{\underline{\underline{L}}}_2 & \cdots & \cdots & \cdots \\ \vdots & \vdots & \ddots & \ddots \end{pmatrix}^T, \qquad (5.3)$$

where the first matrix is lower-triangular by blocks and has only identity blocks on the diagonal, the third matrix is the transpose of the first, and the second matrix is block-diagonal.

Calculating the matrix products on the right-hand side of (5.3) and identifying with (5.1), one obtains the following relations: at iteration 0, $D_0=K$; at iteration 1, $D_0 \cdot \underline{L}_1=\underline{A}$; at iteration 2, $D_0 \cdot \underline{\underline{L}}_2=\underline{C}$ and $\underline{L}_1 \cdot D_0 \cdot \underline{\underline{L}}_1+\underline{D}_1=\underline{B}$, etc. The pieces of the LDLT decomposition can be solved from these relations uniquely and the solution is shown in Table 5.1. Note that we have correctly anticipated the iteration number, shown by the number of underlines, at which the quantities $L_i$ and $D_i$ are solved (left-hand sides in the table) and at which the blocks of $t$ are first used (in the right-hand sides in the table).

<table><caption>**Table 5.1.** Solution for the Cholesky decomposition.</caption><thead><tr><th>iteration</th><th>0</th><th>1</th><th>2</th></tr></thead><tbody><tr><td>moduli</td><td>$D_0=K$</td><td></td><td>$\underline{D}_1=\underline{B}-\underline{A}^T \cdot K^{-T} \cdot \underline{A}$</td></tr><tr><td>coefficients</td><td></td><td>$\underline{L}_1=(K^{-1} \cdot \underline{A})^T$</td><td>$\underline{\underline{\underline{L}}}_2=(K^{-1} \cdot \underline{C})^T$</td></tr></tbody></table>

We now apply the quadratic form associated to $t$ in (5.1) to a vector written in block-vector notation as $\begin{pmatrix} \epsilon^{(0)} \\ \eta \epsilon^{(1)} \\ \vdots \end{pmatrix}$. Each of the subvectors $\epsilon^{(k)}$ is a flattened (rank-1) version of the rank-(k + 2) tensor $\nabla^k \epsilon$ representing the k-th gradient of the macroscopic strain, as earlier in Equation (2.12). Using the Cholesky decomposition (5.3), we have

$$\frac{1}{2} \begin{pmatrix} \epsilon^{(0)} \\ \eta \epsilon^{(1)} \\ \eta^2 \epsilon^{(2)} \\ \vdots \end{pmatrix} \cdot t \cdot \begin{pmatrix} \epsilon^{(0)} \\ \eta \epsilon^{(1)} \\ \eta^2 \epsilon^{(2)} \\ \vdots \end{pmatrix} = \frac{1}{2} (\hat{\epsilon}^{(0)} \cdot D_0 \cdot \hat{\epsilon}^{(0)} + \eta^2 \hat{\epsilon}^{(1)} \cdot \underline{D}_1 \cdot \hat{\epsilon}^{(1)} + \cdots) \qquad (5.4)$$
where the generalized strains $\hat{\varepsilon}^{(0)}$, $\hat{\varepsilon}^{(1)}$, ... are defined as

$$ 
\begin{pmatrix} \hat{\varepsilon}^{(0)} \\ \eta \hat{\varepsilon}^{(1)} \\ \eta^2 \hat{\varepsilon}^{(2)} \\ \dots \end{pmatrix} = \begin{pmatrix} \mathbf{I} & \mathbf{0} & \mathbf{0} & \dots \\ \underline{\mathbf{L}}_1 & \mathbf{I} & \mathbf{0} & \dots \\ \underline{\mathbf{L}}_2 & \dots & \mathbf{0} & \dots \\ \dots & \dots & \dots & \dots \end{pmatrix}^T \begin{pmatrix} \varepsilon^{(0)} \\ \eta \varepsilon^{(1)} \\ \eta^2 \varepsilon^{(2)} \\ \dots \end{pmatrix} \\
= \begin{pmatrix} \varepsilon^{(0)} + \eta \underline{\mathbf{L}}_1^T \cdot \varepsilon^{(1)} + \eta^2 \underline{\mathbf{L}}_2^T \cdot \varepsilon^{(2)} + \dots \\ \eta \varepsilon^{(1)} + \dots \\ \dots \end{pmatrix}. \tag{5.5}
$$

In the right-hand side of (5.4), the matrices $D_i$ are interpreted as homogenized moduli (satisfying the positivity properties established in the forthcoming Section 5.2) and the matrices $L_i$ appearing in (5.5) contain the coefficients entering in the generalized strains.

## 5.2 Positivity of the homogenized moduli

By the same argument as in the original Cholesky decomposition, the positivity of $t$ in (5.2) carries over to the blocks $D_i$:

$$ D_0 > 0, \quad \underline{D}_1 \ge 0, \quad \dots \tag{5.6} $$

This positivity property holds in the sense of quadratic forms, $D_0$ being positive definite while the other $D_i$'s are positive but not necessarily definite.

The proof goes as follows. Inserting (5.1) into (5.3) and extracting the top-left $2 \times 2$ super-block (or $p \times p$ super-block at higher order), we get

$$ 
\begin{pmatrix} \underline{\mathbf{K}} & \underline{\mathbf{A}} \\ \underline{\mathbf{A}}^T & \underline{\mathbf{B}} \end{pmatrix} = \begin{pmatrix} \mathbf{I} & \mathbf{0} & \mathbf{0} & \dots \\ \underline{\mathbf{L}}_1 & \mathbf{I} & \mathbf{0} & \dots \\ \underline{\mathbf{L}}_2 & \dots & \mathbf{0} & \dots \\ \dots & \dots & \dots & \dots \end{pmatrix} \cdot \begin{pmatrix} \underline{\mathbf{D}}_0 & \mathbf{0} & \mathbf{0} & \dots \\ \mathbf{0} & \underline{\mathbf{D}}_1 & \mathbf{0} & \dots \\ \mathbf{0} & \dots & \mathbf{0} & \dots \\ \dots & \dots & \dots & \dots \end{pmatrix} \cdot \begin{pmatrix} \mathbf{I} & \mathbf{0} & \mathbf{0} & \dots \\ \underline{\mathbf{L}}_1 & \mathbf{I} & \mathbf{0} & \dots \\ \underline{\mathbf{L}}_2 & \mathbf{I} & \mathbf{0} & \dots \\ \dots & \dots & \mathbf{I} & \dots \end{pmatrix}^T \tag{5.7} \\
= \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ \underline{\mathbf{L}}_1 & \mathbf{I} \end{pmatrix} \cdot \begin{pmatrix} \underline{\mathbf{D}}_0 & \mathbf{0} \\ \mathbf{0} & \underline{\mathbf{D}}_1 \end{pmatrix} \cdot \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ \underline{\mathbf{L}}_1 & \mathbf{I} \end{pmatrix}^T.
$$

Multiplying by the inverse of the triangular matrices appearing on both the left and right, we obtain

$$ \begin{pmatrix} \underline{\mathbf{D}}_0 & \mathbf{0} \\ \mathbf{0} & \underline{\mathbf{D}}_1 \end{pmatrix} = \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ -\underline{\mathbf{L}}_1 & \mathbf{I} \end{pmatrix} \cdot \begin{pmatrix} \underline{\mathbf{K}} & \underline{\mathbf{A}} \\ \underline{\mathbf{A}}^T & \underline{\mathbf{B}} \end{pmatrix} \cdot \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ -\underline{\mathbf{L}}_1 & \mathbf{I} \end{pmatrix}^T \tag{5.8} $$

which proves (5.6) using (5.2).

## 5.3 Centro-symmetry as a particular case

In the special case where the unit cell is centro-symmetric, all blocks having an odd number of underlines vanish. In the context of second-order homogenization, the $\underline{\mathbf{A}}$ block is the only one vanishing,

$$ \underline{\mathbf{A}} = \mathbf{0}. \tag{5.9} $$

The solution to the Cholesky method in Table 5.1 then boils down to $D_0 = \underline{\mathbf{K}}$, $L_1 = \mathbf{0}$, $D_1 = \underline{\mathbf{B}}$ and $L_2 = (\underline{\mathbf{K}}^{-1} \cdot \underline{\mathbf{C}})^T$, and the decomposition (5.4–5.5) takes the special form

$$ 
\begin{gathered} 
\frac{1}{2} \begin{pmatrix} \varepsilon^{(0)} \\ \eta \varepsilon^{(1)} \\ \eta^2 \varepsilon^{(2)} \\ \dots \end{pmatrix} \cdot \begin{pmatrix} \underline{\mathbf{K}} & \underline{\mathbf{0}} & \underline{\mathbf{C}} & \dots \\ \mathbf{0} & \underline{\mathbf{B}} & \dots & \dots \\ \underline{\mathbf{C}}^T & \dots & \mathbf{0} & \dots \\ \dots & \dots & \dots & \dots \end{pmatrix} \cdot \begin{pmatrix} \varepsilon^{(0)} \\ \eta \varepsilon^{(1)} \\ \eta^2 \varepsilon^{(2)} \\ \dots \end{pmatrix} \\
= \frac{1}{2} \begin{pmatrix} \{\varepsilon^{(0)} + \eta^2 \underline{\mathbf{K}}^{-1} \cdot \underline{\mathbf{C}} \cdot \varepsilon^{(2)} + \dots\} \cdot \underline{\mathbf{K}} \cdot \{\varepsilon^{(0)} + \eta^2 \underline{\mathbf{K}}^{-1} \cdot \underline{\mathbf{C}} \cdot \varepsilon^{(2)} + \dots\} \\ + \{\eta \varepsilon^{(1)} + \dots\} \cdot \underline{\mathbf{B}} \cdot \{\eta \varepsilon^{(1)} + \dots\} \\ + \dots \end{pmatrix} \tag{5.10}
\end{gathered} 
$$

This matches the form of the positive energy functionals which we have obtained for the centro-symmetric examples covered so far, namely the laminate in Equation (2.34), the triangular lattice in Equation (3.20), having $\underline{\mathbf{B}} = \mathbf{0}$, and the honeycomb lattice in Equation (4.17). The positivity of the moduli $D_0 = \underline{\mathbf{K}} > 0$ and $\underline{\mathbf{D}}_1 = \underline{\mathbf{B}} \ge 0$ follows directly from the assumptions in the centro-symmetric case, see (5.2).

In the centro-symmetric case, the Cholesky decomposition therefore reduces to the positivity-preserving truncation method used so far.

# 6 A non-symmetric pantograph

In this section, we analyze a variant of the classical pantograph truss, shown in Figure 6.1. In this variant, the springs shown in black in the figure are stiffer than the springs shown in blue, which breaks the left-right symmetry and makes the lattice non-centro-symmetric. As a result, the homogenized one-dimensional model couples the strain $\varepsilon(x)$ and the strain gradient $\varepsilon'(x)$, as we will see.
The homogenized one-dimensional energy is derived using standard high-order homogenization method first. Next, we apply the Cholesky decomposition from Section 5 to the resulting energy to truncate it in a way that preserves positivity.

Spring constants $2k/\eta$ and $k/\eta$ are assigned to the black and blue springs, respectively. The vertical springs shown by dashed lines in Figure 6.1 are assigned a spring constants $k\chi/\eta$, where $\chi$ is the stiffness contrast parameter. Pantographs are known to feature a floppy mode in the limit $\chi \to 0$ and this floppy mode is coupled to macroscopic stretching [S18b, BGP17, DLSS22, SAI11]. When we homogenize the pantograph, the parameter $\chi$ is treated as a symbol whose value remains unspecified. We discuss the finite-contrast case $\chi = \mathcal{O}(1)$ first, and the high-contrast limit $\chi \to 0$ next.

<figure><img src="image_5.png" /><figcaption>**Figure 6.1.** A non-symmetric pantograph. The black springs have been made stiffer than the blue springs so as to break the mirror-symmetry with respect to a vertical plane. The spring families $\phi$ are shown using Roman numerals. The spring constant $k\chi\eta^{-1}$ of the vertical springs (dashed line) is parameterized by $\chi$. In the limit of large contrast, $\chi \to 0$, the pantograph displays a soft mode whereby all right triangles move rigidly and the pantograph stretches or contracts horizontally.</figcaption></figure>

## 6.1 Summary of homogenization

The pantograph shown has 13 different types of springs, labelled using Roman numerals, and 6 types of nodes (Bravais sub-lattices), shown using different colors. Except for the choice of spring constants, this pantograph is identical to that studied in section 7.5 of [YAL24]. We homogenize it automatically again using the shoal library in the companion notebook homogenize-pantograph-non-symmetric nb, provided as supplementary material [TAL24b]. The homogenization follows the same steps as that of the honeycomb (Section 4) and is almost identical to that presented in Section 7.5 of [YAL24] and we limit ourselves to a summary of the main results here.

The homogenization limit $\eta \to 0$ is approached as follows. Working with a symbolic value of the contrast parameter $\chi$, the library calculates the expressions of the microscopic shifts $\xi_\phi$ of the different Bravais sub-lattices (indexed by $\phi = I, \dots, XIII$) as an expansion in terms of the successive derivatives of the one-dimensional macroscopic strain $\varepsilon(x) = u'(x)$, where $u(x)$ is the (scalar) macroscopic displacement. Having determined the microscopic shifts, it proceeds to derive the expressions of the rescaled change of length $\varepsilon_\phi$ of the different types of spring. They are collected in a vector $\mathbf{E} = (\varepsilon_1, \dots, \varepsilon_{XIII}) \in \mathbb{R}^{13}$. Using Taylor expansions, it is written in terms of the macroscopic strain $\varepsilon(x)$ in a similar form as earlier in (4.8),

$$ \mathbf{E}(x, \chi) = \varepsilon(x) \mathbf{F}_0(\chi) + \eta \varepsilon'(x) \mathbf{F}_1(\chi) + \eta^2 \varepsilon''(x) \mathbf{F}_2(\chi) + \mathcal{O}(\eta^3). \quad (6.1) $$

The strain relocalization vectors $\mathbf{F}_0(\chi) = (\frac{33}{69} + \frac{\sqrt{3}}{181}\chi, \dots)$, $\mathbf{F}_1(\chi) = (\frac{3}{2} + \frac{9}{181}\chi + \frac{8}{69}, \dots)$ and $\mathbf{F}_2^*(\chi) = (\frac{60-121\chi}{432\sqrt{3}(181\chi+69)}, \dots)$ are output by the library, see the companion notebook. A $\Psi_2^*$ contribution to $\mathbf{F}_2$ is omitted in $\mathbf{F}_2^*$, which has no consequence, see Remark 3.1.

The (discrete) strain energy of the lattice can be approximated by the integral

$$ \Phi[\mathbf{E}] = \frac{1}{2} \int_{-\infty}^{+\infty} \mathbf{E}(x, \chi) \cdot \mathbf{k}(\chi) \cdot \mathbf{E}(x, \chi) \, dx, \quad (6.2) $$

where

$$ \mathbf{k}(\chi) = k \operatorname{diag}(2, \dots, 2, 1, \dots, 1, \chi) \quad (6.3) $$

is the matrix filled with the spring constants $2k/\eta$, $k/\eta$ and $\chi k/\eta$ on the diagonal, times the density $\rho = \eta^{-1}$ of each particular type of spring per unit length, times a factor $\eta^2$ reflecting the convention that the change of the length of springs is $\eta \mathbf{E}$ and not just $\mathbf{E}$.

Inserting (6.1) into (6.2), one obtains an expansion of the energy of the form

$$ \Phi^*[\varepsilon] = \int_{-\infty}^{+\infty} \frac{1}{2} \begin{pmatrix} \varepsilon(x) \\ \eta \varepsilon'(x) \\ \eta^2 \varepsilon''(x) \end{pmatrix} \cdot \mathbf{t}(\chi) \cdot \begin{pmatrix} \varepsilon(x) \\ \eta \varepsilon'(x) \\ \eta^2 \varepsilon''(x) \end{pmatrix} \, dx + \mathcal{O}(\eta^3) \quad (6.4) $$

where $\mathbf{t}$ is a symmetric $3 \times 3$ matrix

$$ \mathbf{t}(\chi) = \begin{pmatrix} K(\chi) & A(\chi) & C(\chi) \\ A(\chi) & B(\chi) & \dots \\ C(\chi) & \dots & \dots \end{pmatrix}, \quad (6.5) $$
whose (non-asymptotic) higher-order contributions, denoted as ellipses, are ignored here. The coefficients of $t$ are given by similar formulas as earlier in (2.24),

$$ 
\begin{array}{c|ccc}
\text{order} & \text{diagonal blocks} & \text{off-diagonal blocks} \\
\eta^0 & K(\chi) = \mathbf{F}_0(\chi) \cdot \mathbf{k}(\chi) \cdot \mathbf{F}_0(\chi) & \\
\eta^1 & A(\chi) = \mathbf{F}_0(\chi) \cdot \mathbf{k}(\chi) \cdot \mathbf{F}_1(\chi) & \\
\eta^2 & B(\chi) = \mathbf{F}_1(\chi) \cdot \mathbf{k}(\chi) \cdot \mathbf{F}_1(\chi) & C(\chi) = \mathbf{F}_0(\chi) \cdot \mathbf{k}(\chi) \cdot \mathbf{F}_2(\chi)
\end{array}
\qquad (6.6) $$

Their expressions and their particular values for $\chi=1$ and $\chi=0$ are listed in Table 6.1.

Note the presence of the $A(\chi)$ term coupling the effective strain $\epsilon$ and its gradient $\epsilon'$, which is a consequence of the broken symmetry.

$$ 
\begin{array}{c|cccc}
 & K(\chi)/k & A(\chi)/k & B(\chi)/k & C(\chi)/k \\
\hline
\text{general case} & \frac{276\chi}{69+181\chi} & \frac{2\chi(69+85\chi)}{(69+181\chi)^2} & \frac{552-19\chi+638\chi^2}{(69+181\chi)^2} & \frac{\chi(1149+17242\chi)}{18(69+181\chi)^2} \\
\chi=1 & 1.104 & 4.928 \times 10^{-3} & 1.874 \times 10^{-2} & 1.635 \times 10^{-2} \\
\chi=0 & 0 & 0 & 1.159 \times 10^{-1} & 0
\end{array}
$$

<caption>**Table 6.1.** Scaled values of the homogenized elastic constants of the pantograph: general case, case of low elastic contrast ($\chi=1$) and limit of infinite elastic contrast discussed in Section 6.4 ($\chi=0$).</caption>

Since $\mathbf{k}(\chi) \ge 0$ is a positive quadratic form and $\mathbf{t}(\chi)$ in (6.5–6.6) can be rewritten as $\mathbf{t}(\chi) = (\mathbf{F}_0(\chi) \ \mathbf{F}_1(\chi) \ \mathbf{F}_2(\chi))^T$, $\mathbf{k}(\chi) \cdot (\mathbf{F}_0(\chi) \ \mathbf{F}_1(\chi) \ \mathbf{F}_2(\chi))$, we conclude that $\mathbf{t}(\chi) \ge 0$ is a positive quadratic form. In particular, its top-left $2 \times 2$ sub-block is a positive quadratic form: both the first diagonal entry $K$ and the determinant of the latter are positive,

$$ \forall \chi > 0, \quad K(\chi) > 0 \quad \text{and} \quad K(\chi) B(\chi) \ge A^2(\chi). \qquad (6.7) $$

The strict positivity of $K(\chi) = \mathbf{F}_0(\chi) \cdot \mathbf{k}(\chi) \cdot \mathbf{F}_0(\chi)$ comes from the fact that $\mathbf{k}(\chi) > 0$ and $\mathbf{F}_0(\chi) \neq 0$ for $\chi > 0$. These inequalities can also be checked directly using the expressions in Table 6.1.

In the limit $\chi \to 0$, $K(\chi) \to 0$ and the pantograph becomes soft, as discussed in Section 6.4 below.

## 6.2 Naïve truncation of the energy

When the energy $\Phi^*[\epsilon]$ (6.4) is truncated beyond order $\eta^2$, we get

$$ \Phi^{\dagger}[\epsilon] = \int_{-\infty}^{+\infty} \frac{1}{2} (K(\chi) \epsilon^2 + 2 \eta A(\chi) \epsilon \epsilon' + \eta^2 (B(\chi) \epsilon'^2 + 2C(\chi) \epsilon \epsilon'')) \, dx. \qquad (6.8) $$

Both the terms $A \epsilon \epsilon' = \frac{d}{dx} \left( \frac{A \epsilon^2}{2} \right)$ and $C \epsilon \epsilon'' = \frac{d}{dx} (C \epsilon \epsilon') - C \epsilon'^2$ can be integrated by parts. Upon neglecting the boundary terms, we obtain the approximation

$$ \Phi^{\dagger}[\epsilon] = \int_{-\infty}^{+\infty} \frac{1}{2} (K(\chi) \epsilon^2 + \eta^2 B^{\dagger}(\chi) \epsilon'^2) \, dx \qquad (6.9) $$

where

$$ B^{\dagger}(\chi) = B(\chi) - 2C(\chi). \qquad (6.10) $$

The quantity $B^{\dagger}(\chi)$ is not necessarily positive, unlike the original gradient modulus $B(\chi)$ which is a diagonal entry of $\mathbf{t}(\chi) \ge 0$, see also (6.7). For $\chi=1$, for instance, $B^{\dagger}(1) = -1.396 \times 10^{-2} k < 0$, making the energy $\Phi^{\dagger}[\epsilon]$ in (6.9) non-positive. The set of positive values of $\chi$ that make $B^{\dagger}(\chi)$ positive can be worked out as $\chi \in (0, 0.602)$: we will explain in Section 6.4 why $B^{\dagger}(\chi)$ has to be positive in the limit $\chi \to 0$ where the lattice becomes floppy.

## 6.3 Cholesky-based truncation preserving positivity

We restart from the homogenized energy given in (6.4), and apply Cholesky decomposition from Section 5 to the matrix $\mathbf{t}$ appearing in (6.5). Since there is only one scalar macroscopic strain measure $\epsilon$ in the one-dimensional case, we can ignore the block structure of $\mathbf{t}$ which we assumed in Section 5 for the sake of generality.

Setting $\epsilon^{(0)} = (\epsilon)$, $\epsilon^{(1)} = (\epsilon')$, $\epsilon^{(2)} = (\epsilon'')$, Equations (5.4–5.5) yield an alternative expression of the energy which is accurate to order $\eta^2$ included,

$$ \Phi_+[\epsilon] = \int_{-\infty}^{+\infty} \left( \frac{1}{2} D_0(\chi) \left\{ \epsilon(x) + \eta L_1(\chi) \epsilon'(x) + \eta^2 L_2(\chi) \epsilon''(x) \right\}^2 + \frac{\eta^2}{2} D_1(\chi) \epsilon'^2(x) \right) dx. \qquad (6.11) $$

The homogenized moduli $D_i(\chi)$ and coefficients $L_i(\chi)$ are read off from Table 5.1 as

$$ 
\begin{aligned}
D_0(\chi) &= K(\chi) & L_1(\chi) &= \frac{A(\chi)}{K(\chi)} \\
D_1(\chi) &= B(\chi) - \frac{A^2(\chi)}{K(\chi)} & L_2(\chi) &= \frac{C(\chi)}{K(\chi)}.
\end{aligned}
\qquad (6.12) $$
Combining (6.11) and (6.12), we obtain the positive, order-$\eta^2$-accurate expression of the energy for the pantograph as

$$ \Phi_+[\varepsilon] = \int_{-\infty}^{+\infty} \left( \frac{1}{2} K(\chi) \left\{ \varepsilon(x) + \eta \frac{A(\chi)}{K(\chi)} \varepsilon'(x) + \eta^2 \frac{C(\chi)}{K(\chi)} \varepsilon''(x) \right\}^2 + \frac{\eta^2}{2} \left( B(\chi) - \frac{A^2(\chi)}{K(\chi)} \right) \varepsilon'^2(x) \right) dx. \quad (6.13) $$

In the non-symmetric case, the generalized strain $\{\dots\}$ includes a correction of order $\eta$ proportional to $\varepsilon'$.

In the finite-contrast case, $\chi > 0$, the positivity of the energy $\Phi_+$ follows from (6.7), see also the general proof in Section 5.2. The limit of large elastic contrast, $\chi \to 0$, is studied next.

## 6.4 Floppy limit, $\chi \to 0$

In this section, we discuss the form of the energy $\Phi_+$ in (6.13) in the limit of large elastic contrast, $\chi \to 0$, where a floppy mode coupled to macroscopic extension appears. Specifically, we explain why the naïve procedure from Section 6.2 produces a positive energy in this particular case, $B^{\ddagger}(0) \ge 0$.

From the expressions of $K(\chi)$, $A(\chi)$ and $C(\chi)$ appearing in the second row of Table 6.1, one can obtain the following equivalents in the limit $\chi \to 0$

$$ K(\chi) \sim 4k\chi, \quad A(\chi) \sim \frac{2k}{69}\chi, \quad C(\chi) = \frac{383k}{28566}\chi, \quad (6.14) $$

whereas $B(\chi)$ converges to a positive constant

$$ B(\chi) \to \frac{8k}{69}. \quad (6.15) $$

The scaling behavior observed in (6.14–6.15) can be explained as follows. For $\chi=0$, the diagonal matrix of microscopic moduli $k(0)$ in (6.3) has a zero entry in the last, 13th, position on the diagonal and the first 12 components of the vector $F_0(0)$ vanish, implying $k(0) \cdot F_0(0) = 0$. Now, from (6.6), the homogenized constants $K(\chi)$, $A(\chi)$ and $C(\chi)$ appear to be the dot product of $k(\chi) \cdot F_0(\chi) \to 0$ by one of the $F_i(\chi)$'s, but $B(\chi)$ is not of this form. This is what makes $K(\chi)$, $A(\chi)$ and $C(\chi)$ proportional to $\chi$ in the limit $\chi \to 0$, whereas $B(\chi)$ has a non-zero limit. This behavior appears to be generic in structures possessing a soft mode, which is such that $k \cdot F_0$ is not full-rank, see [DLSS22].

In the limit $\chi \to 0$, the ratios $A/K$ and $C/K$ in (6.13) are therefore finite, whereas both $K$ and $A^2/K$ go to zero. In view of (6.13), the homogenized energy relevant to the soft limit is therefore

$$ \Phi_+[\varepsilon] = \int_{-\infty}^{+\infty} \frac{\eta^2}{2} B(0) \varepsilon'^2(x) dx. \quad (6.16) $$

On the other hand, the energy produced by the naïve truncation method in (6.9) becomes, in the limit,

$$ \Phi^{\ddagger}[\varepsilon] = \int_{-\infty}^{+\infty} \frac{\eta^2}{2} B^{\ddagger}(0) \varepsilon'^2 dx. \quad (6.17) $$

Equation (6.10) shows that $B^{\ddagger}(0) = B(0)$ as $C(0) = 0$ by (6.14)3. The energies (6.16) and (6.17) produced by the naïve and Cholesky methods therefore coincide in the presence of a soft mode, when the regular, leading-order elasticity is degenerate. This explains why the naïve truncation method predicts a correct (positive) strain-gradient elasticity in this soft limit.

# 7 Extension to dimension reduction: bending of an elastic slab

In this section we show that our approach extends naturally to dimension reduction. We analyze the bending of a linearly elastic slab in dimension 2. Specifically, we revisit the work of [LA20] and derive a second-order energy for the bending of an infinitely long rectangular 2D beam, fixing the non-positivity of the gradient-of-curvature modulus reported in this previous work.

## 7.1 Summary of dimension reduction

<figure><img src="image_6.png" /><figcaption>**Figure 7.1.** A linearly elastic slab in dimension 2 having Young modulus $Y$ and Poisson's ratio $\nu$. A second-order equivalent beam bending model is derived by dimension reduction, see Equation (7.13).</figcaption></figure>

The rectangular domain of the slab is shown in Figure 7.1. It is infinite in the longitudinal $x_1$ direction and its width is denoted as $\eta$. We introduce the stretched transverse coordinate

$$ y = \frac{x_2}{\eta}. \quad (7.1) $$
The stretched coordinate $y$ varies in the range $-\frac{1}{2} \le y \le +\frac{1}{2}$.

Dimension reduction addresses the limit $\eta \ll 1$, see [LA20], and can be summarized as follows. The displacement field relevant to bending is of the form (see Section 4.2 of [LA20])

$$ \mathbf{u}(x_1, y) = [-\eta y V'(x_1) + \eta^2 \xi_1(x_1, y)] \mathbf{i}_1 + [V(x_1) + \eta^2 \xi_2(x_1, y)] \mathbf{i}_2, \quad (7.2) $$

where $V(x_1)$ denotes the macroscopic displacement in the transverse direction, $\xi_i(x_1, y)$ are the components of the microscopic displacement and $(\mathbf{i}_1, \mathbf{i}_2)$ is the Cartesian basis. The latter is subject to the constraint

$$ \forall x_1, \quad \int_{-1/2}^{+1/2} \xi_i(x_1, y) \, dy = 0 \qquad (7.3) $$

warranting that the macroscopic displacement $V(x_1) \mathbf{i}_2$ is the average of the displacement $\mathbf{u}(x_1, y)$ in any particular cross-section with coordinate $x_1$.

Using asymptotic analysis, the microscopic displacement that minimizes the elastic energy for a prescribed deflection profile $V$ has been obtained in [LA20] in the form

$$ \begin{aligned} \xi_1(x_1, y) &= 0 & -y \left( \frac{6+5\nu}{24} - \frac{2+\nu}{6}y^2 \right) \eta \kappa'(x_1) + \dots \\ \xi_2(x_1, y) &= \nu \kappa(x_1) \left( \frac{y^2}{2} - \frac{1}{24} \right) + 0 \times \eta & + \dots \end{aligned} \quad (7.4) $$

where $\kappa(x_1)$ denotes the macroscopic bending strain,

$$ \kappa(x_1) = V''(x_1), \quad (7.5) $$

and $\nu$ is the Poisson's ratio of the material appearing in Equation (7.8) below. The solution (7.4) satisfies the constraint (7.3).

The microscopic 2D strain $(\nabla_x \mathbf{u} + \nabla_x \mathbf{u}^T)/2$ is obtained from (7.2) in the form $\eta \mathbf{E}(x_1, y)$, where

$$ \mathbf{E}(x_1, y) = \left( -y \kappa(x_1) + \eta \frac{\partial \xi_1}{\partial x_1} \right) \mathbf{i}_1 \otimes \mathbf{i}_1 + \left( \frac{\partial \xi_1}{\partial y} + \eta \frac{\partial \xi_2}{\partial x_1} \right) \frac{\mathbf{i}_1 \otimes \mathbf{i}_2 + \mathbf{i}_2 \otimes \mathbf{i}_1}{2} + \frac{\partial \xi_2}{\partial y} \mathbf{i}_2 \otimes \mathbf{i}_2. \quad (7.6) $$

Inserting (7.4), we obtain

$$ \begin{aligned} \mathbf{E}(x_1, y) = & \left[ -y (\mathbf{i}_1 \otimes \mathbf{i}_1 - \nu \mathbf{i}_2 \otimes \mathbf{i}_2) \right] \kappa(x_1) \\ & + \left[ -(1+\nu) \left( \frac{1}{4} - y^2 \right) \frac{\mathbf{i}_1 \otimes \mathbf{i}_2 + \mathbf{i}_2 \otimes \mathbf{i}_1}{2} \right] \eta \kappa'(x_1) \\ & + \left[ -y \left( \frac{6+5\nu}{24} - \frac{2+\nu}{6}y^2 \right) \mathbf{i}_1 \otimes \mathbf{i}_1 + \dots \right] \eta^2 \kappa''(x_1) \\ & + \dots \end{aligned} \quad (7.7) $$

The quantities appearing in square brackets in (7.7) have been denoted as $\mathbf{F}_i$ earlier, see for instance Equation (2.16). In the fourth line of (7.7), the ellipsis stands for a contribution coming from the order-$\eta^2$ term in the displacement, that has not been resolved in (7.4). This contribution does not contribute to the energy at order-$\eta^2$, however, by the argument already given in Remark 2.1.

The linear theory of elasticity in 2D yields the strain energy of the slab in the form

$$ \Phi = \frac{\eta^3 Y}{2} \int_{-\infty}^{+\infty} \int_{-1/2}^{+1/2} \left( \frac{\nu}{1-\nu^2} \mathrm{tr}^2 \mathbf{E}(x_1, y) + \frac{1}{1+\nu} \mathbf{E}(x_1, y) : \mathbf{E}(x_1, y) \right) dx_1 dy, \quad (7.8) $$

where $Y$ and $\nu$ are the Young modulus and Poisson's ratio, respectively, in the sense of two-dimensional elasticity. By using the identity $\mathbf{E}:\mathbf{E}=\mathbf{E}^D:\mathbf{E}^D+\frac{1}{2}\mathrm{tr}^2\mathbf{E}$, where $\mathbf{E}^D=\mathbf{E}-\frac{\mathbf{I}_2}{2}\mathrm{tr}\mathbf{E}$ is the deviatoric part, one can show that the range of elastic parameters corresponding to a stable material are

$$ Y > 0 \text{ and } -1 < \nu < 1. \quad (7.9) $$

Note that the incompressible limit corresponds to $\nu=1$ in 2D, and not to $\nu_{3D}=1/2$ as in 3D.

The prefactor $\eta^3$ in (7.8) comes from the fact that the actual 2D strain is $\eta \mathbf{E}$ (and not just $\mathbf{E}$) and that we are using the stretched coordinate $y$ in the transverse integration. For the energy to be finite in the limit $\eta \to 0$, we can conveniently assume that the Young modulus scales as $Y = \eta^{-3} \bar{Y}$, where $\bar{Y} = \mathcal{O}(1)$.

Inserting the expression (7.7) of the strain into (7.8) and integrating with respect to the transverse coordinate $y$, we obtain

$$ \Phi^*[\kappa] = \frac{\bar{Y}/12}{2} \int_{-\infty}^{+\infty} \left( \kappa^2(x_1) + \frac{2}{5} \left(1 + \frac{11\nu}{12}\right) \eta^2 \kappa(x_1) \kappa''(x_1) + \frac{1+\nu}{5} \eta^2 \kappa'^2(x_1) \right) dx_1 \quad (7.10) $$

where $\bar{Y}/12$ is the scaled bending modulus predicted by classical dimension reduction. Due to the left-right symmetry, there is no energy contribution coupling $\epsilon$ and $\epsilon'$ (A terms).

The form $\Phi^*=\frac{1}{2}\int(K\kappa^2+\eta^2(2C\kappa\kappa''+B\kappa'^2))dx_1$ of the energy matches that relevant to homogenization, see for instance (6.8) for the pantograph. The coefficients can be identified from (7.10) as

$$ \begin{aligned} K &= \frac{\bar{Y}}{12} \\ C &= \frac{K}{5} \left(1 + \frac{11\nu}{12}\right) \\ B &= \frac{K}{5} (1 + \nu). \end{aligned} \quad (7.11) $$
By the same arguments as earlier, the homogenized moduli $K > 0$ and $B \ge 0$ are positive over the range (7.9) of permissible elastic parameters. This is a consequence of the positivity of the energy which we started from in (7.8). The positivity of $C$ is coincidental.

## 7.2 Naïve truncation of the energy

Upon integration by parts of the $\kappa \kappa''$ term in (7.10) and after discarding the boundary terms, we recover the non-positive functional

$$ \Phi^{\ddagger}[\kappa] = \frac{\bar{Y}/12}{2} \int_{-\infty}^{+\infty} \left( \kappa^2(x_1) - \frac{1}{5} \left(1 + \frac{5\nu}{6}\right) \kappa'^2(x_1) \right) dx_1, \quad (7.12) $$

which was the final result of [LA20], see their Equation [4.10]. The expression $B^{\ddagger} = B - 2C = -\frac{K}{5}(1 + \frac{5}{6}\nu)$ of the apparent gradient modulus is similar to that derived earlier in the context of homogenization, see (2.29) and (6.10). The apparent modulus is negative, $B^{\ddagger} < 0$, which makes $\Phi^{\ddagger}$ non-convex.

## 7.3 Truncation preserving positivity

Applying instead the method from Section 2.6 to the one-dimensional energy $\Phi^*$ in (7.10), we get

$$ \Phi_{+}[\kappa] = \frac{\bar{Y}/12}{2} \int_{-\infty}^{+\infty} \left( \left\{ \kappa(x_1) + \frac{1}{5} \left(1 + \frac{11\nu}{12}\right) \eta^2 \kappa''(x_1) \right\}^2 + \frac{1+\nu}{5} \eta^2 \kappa'^2(x_1) \right) dx_1. \quad (7.13) $$

This is a convex, second-order-accurate approximation of $\Phi^*$. It is akin to the other energy functionals $\Phi_+$ derived earlier, see for instance in (2.34) for the laminate.

**Remark 7.1. (geometric compatibility)** The generalized strain appearing in (7.13) can be rewritten as $\{\kappa(x_1) + \frac{1}{5}(\cdots) \eta^2 \kappa''(x_1)\} = \hat{V}''(x_1)$ in terms of the alternate macroscopic displacement $\hat{V}$,

$$ \hat{V}(x_1) = V(x_1) + \frac{1}{5} \left(1 + \frac{11\nu}{12}\right) \eta^2 V''(x_1). \quad (7.14) $$

In this sense, the generalized bending strain is geometrically compatible, as is any strain measure in dimension 1. This yields a compact and positive expression of the homogenized energy in the form

$$ \frac{\bar{Y}/12}{2} \int_{-\infty}^{+\infty} \left( \hat{V}'''(x_1) + \frac{1+\nu}{5} \eta^2 \hat{V}'''(x_1) (x_1) \right) dx_1, \quad (7.15) $$

which differs from the energy $\Phi_+[\kappa]$ in (7.13) by negligible terms, of order $\eta^4$, coming from subdominant contributions to the strain gradient.

# 8 Comparison to existing energy formulations

In this Section, we compare our method with two second-order energy formulations from the literature. Our goal is to assess how they handle distributed forces that vary on a significantly shorter lengthscale $\delta$ than the macroscopic scale. Such force profiles defeat the scale separation assumption and the accuracy of the homogenized models is not guaranteed. Quickly-varying force profiles are still interesting as they reveal that some models tend to regularize the applied force profile (as one could hope thanks to the gradient effect) while others do not.

## 8.1 Comparison method

The comparison is made on a low-contrast version of the 1D pantograph from Section 6, with the contrast parameter set to $\chi = 1$. We drop the argument $\chi$ throughout, i.e., we denote $K(1) = K$, etc. We use a profile $f_{\delta}(x)$ of distributed force

$$ f_{\delta}(x) = \frac{1}{4} \left( 2 + \operatorname{erf} \frac{1-x}{\delta} + \operatorname{erf} \frac{1+x}{\delta} \right), \quad (8.1) $$

which is a box profile regularized with a parameter $\delta$, see the black curve in Figure 8.1(b). In the numerical examples given later, we use a relatively sharp profile by choosing $\delta = 0.03 \ll 1$.

Using prior work, we generate different positive energy functionals $\Phi$:

* the Boussinesq trick [ABV16] produces an energy

$$ \Phi = \int_{-\infty}^{+\infty} \left( \frac{1}{2} K \varepsilon^2 - \left( f_{\delta} + \eta^2 \frac{B^{\ddagger}}{K} f_{\delta}'' \right) u \right) dx, \quad (8.2) $$

* the truncation limited to strain proposed by [SC00] yields

$$ \Phi = \int_{-\infty}^{+\infty} \left[ \int_{-\infty}^{+\infty} \frac{1}{2} \begin{pmatrix} \varepsilon(x) \\ \eta \varepsilon'(x) \\ \eta^2 \varepsilon''(x) \end{pmatrix} \cdot \begin{pmatrix} K & A & C \\ A & B & t_{23} \\ C & t_{23} & t_{33} \end{pmatrix} \cdot \begin{pmatrix} \varepsilon(x) \\ \eta \varepsilon'(x) \\ \eta^2 \varepsilon''(x) \end{pmatrix} dx - f_{\delta} u \right] dx, \quad (8.3) $$
where the coefficients $t_{23}$ and $t_{33}$ are given in the first column of Table 8.1.

Our method (minimal truncation) produces an energy $\Phi$ in (6.13) of the form (8.3) but with the coefficients $t_{23}$ and $t_{33}$ given in the second column of Table 8.1.

Our comparison reveals that the Boussinesq trick produces marked oscillations at $x = \pm 1$, where the distributed force varies quickly, by a Gibbs phenomenon. These oscillations are due to the presence of the second derivative of the load $f_{\delta}''$ in (8.2). The truncation limited to strain and the minimal truncation both take advantage of the regularizing effect of the strain gradients $\epsilon'$ and $\epsilon''$ in the elastic energy, and produce both a much gentler response.

<table><caption>**Table 8.1.** Higher-order elastic constants for the pantograph with contrast $\chi = 1$, used in Equation (8.3). The constants are given by analogous formulas to those in (2.24) for the laminate. As the truncation limited to strain requires the full expression of $\mathbf{F}_2$ and not just the part $\mathbf{F}_2^2$ that does not depend on $\Psi_2$, we have pushed homogenization to order $\eta^4$ in energy to calculate $\Psi_2$ and the full $\mathbf{F}_2$ for the purpose of comparison.</caption><thead><tr><th></th><th>truncation limited to strain</th><th>minimal truncation</th></tr></thead><tbody><tr><td>$t_{23}$</td><td>$\mathbf{F}_1 \cdot \mathbf{k} \cdot \mathbf{F}_2 = 1.874 \times 10^{-2} k$</td><td>$\frac{AC}{K} = 7.297 \times 10^{-5} k$</td></tr><tr><td>$t_{33}$</td><td>$Z = \mathbf{F}_2 \cdot \mathbf{k} \cdot \mathbf{F}_2 = 8.110 \times 10^{-4} k$</td><td>$\frac{C^2}{K} = 2.421 \times 10^{-4} k$</td></tr></tbody></table>

## 8.2 Justification of the energies

The derivation of the energy (8.2) starts from the naïve truncation, which produces an energy of the form (8.3) with $t_{23} = t_{33} = 0$. This energy is not positive and we make it positive using the Boussinesq trick as follows. We write its stationary condition as $K \epsilon'(x) - \eta^2 B^\dagger \epsilon'''(x) + f_\delta(x) = 0$, and solve this equation order by order in $\eta$ as $\epsilon' = -K^{-1}(f_\delta(x) + \eta^2 B^\dagger f_\delta''(x) / K) + \mathcal{O}(\eta^3)$. This is nothing but the minimizer of (8.2) up to terms of order $\mathcal{O}(\eta^3)$.

The energy (8.3) and the coefficients in the first column of Table 8.1 (truncation limited to strain) have been obtained by applying the method from Section 2.5 to the energy expansion (6.4-6.5) of the pantograph. The details are omitted.

The coefficients in the second column of Table 8.1 (minimal truncation, our method) are obtained from Section 6.3, see Equation (6.13).

Note that we have restored the potential of the external force in (8.2-8.3), which we had omitted so far.

## 8.3 Equilibrium solutions

By minimizing the energy (8.2), we get the equilibrium equation

$$K \epsilon'(x) + f_{\delta}(x) + \eta^2 \frac{B^{\dagger}}{K} f_{\delta}''(x) = 0, \quad (\text{Boussinesq/naïve truncation}), \qquad (8.4)$$

where $B^{\dagger} = B - 2C < 0$ from (6.10). By minimizing the energy (8.3) instead, we get

$$K \epsilon'(x) - \eta^2 B^{\dagger} \epsilon^{(3)}(x) + \eta^4 t_{33} \epsilon^{(5)}(x) + f_{\delta}(x) = 0, \quad (\text{truncation limited to strain, minimal truncation}) \qquad (8.5)$$

where $t_{33}$ depends on the truncation method, see Table (6.1).

The above equations are solved using Fourier analysis as follows. The transfer function $H(\eta \omega)$ characterizing the response ($-K \epsilon'$) to a harmonic force profile $f \propto \cos \omega x$ with wavenumber $\omega$ is identified from (8.4) and (8.5) as

$$H(\eta \omega) \quad \begin{array}{c|ccc} \text{Boussinesq/naïve tr.} & \text{tr. limited to strain} & \text{minimal truncation} \\ \hline & \frac{1}{1 + \frac{\eta^2 B^{\dagger}}{K} \omega^2 + \frac{\eta^4 Z}{K} \omega^4} & \frac{1}{1 + \frac{\eta^2 B^{\dagger}}{K} \omega^2 + \frac{\eta^4 C^2}{K^2} \omega^4} \end{array} \qquad (8.6)$$

The transfer functions are shown in Figure 8.1(a). They yield the Fourier transform $\hat{K} \epsilon'(\omega)$ of the response as

$$-\hat{K} \epsilon'(\omega) = H(\eta \omega) \hat{f}_{\delta}(\omega), \qquad (8.7)$$

where $\hat{f}_{\delta}$ is the Fourier transform of the applied distributed force (8.1). The responses $-K \epsilon'(x)$ of the different models shown in Figure 8.1(b') have been produced by calculating numerically the inverse Fourier transform of the left-hand side of (8.7).

In the regions where $f_\delta(x)$ is constant, that is away from the points $x = \pm 1$, all models agree with the prediction $-K \epsilon' = f_\delta + \mathcal{O}(\eta)$ of the equivalent bar model at leading order. Around the points $x = \pm 1$ where $f_\delta(x)$ varies quickly, however, the three models behave differently. The Boussinesq trick produces a marked Gibbs phenomenon as the $\eta^2 f_\delta'' \sim (\eta / \delta)^2$ term blows up when $\delta \ll \eta$. In Fourier space, the function $H(\eta \omega)$ for this model is unbounded and the high-frequency content of the applied force is enhanced. By contrast, both the truncation limited to strain and the minimal truncation have a bounded transfer function $H(\eta \omega)$ which produces gentler oscillations.

Despite this essential difference in their response to short-wavelengths, all models agree at long wavelengths as they are all asymptotically exact, $H(\eta \omega) = 1 - \frac{\eta^2 B^{\dagger}}{K} \omega^2 + \mathcal{O}((\eta \omega)^3)$.
<figure><img src="image_7.png" /><figcaption>**Figure 8.1.** Comparison of the predictions of different second-order models: naïve truncation fixed by the Boussinesq trick (red), truncation limited to strain (blue) and the minimal truncation proposed in this paper (green). The plots were generated with a scale separation parameter $\eta = 0.3$ in the 1D pantograph, and a regularization parameter $\delta = 0.03$ in the force profile. (a) Transfer functions mapping $f_{\delta}$ to $-K\epsilon'$ in Fourier space. (b) Regularized-box profile of distributed force $f_{\delta=0.03}(x)$, and (b') response predicted by the various models.</figcaption></figure>

Although the truncation limited to strain and the minimal truncation handle quickly-varying forces in a qualitatively similar way, the latter is simpler to implement as it does not require the calculation of the second-order localizer $\Psi_2$ for the displacement, see Remark 2.1 and the discussion in subsection 2.5.

# 9 Discussion and conclusion

We propose a truncation method for energy functionals produced by higher-order homogenization. It solves a long-standing issue concerning the presence of non-positive strain-gradient stiffnesses. This issue made simulations of homogenized models ill-posed and may explain the limited popularity of high-order homogenization. In our approach, the non-positive, higher-order homogenized constants are revisited as coefficients entering in the definition of generalized strain measures that are used as arguments of a positive quadratic form. The positivity of the energy $\Phi_+$ obtained in this way is guaranteed by the positivity of the initial microscopic energy. The truncation preserves the order of accuracy of the homogenization procedure.

We reiterate that we cannot generate effective boundary conditions until we properly address boundary layers, which we will cover in follow-up work. This limits temporarily the applicability of the proposed approach to infinite continua, as we did in Section 8.

The positive functional $\Phi_+$ produced by truncation depends solely and uniquely on the homogenized constants that are obtained by homogenization at a given order $p$. This is an improvement over the earlier truncation method of [SC00] which shares the key positivity property with ours and handles irregular forces similarly, but requires calculating a number of homogenized elastic moduli beyond order $p$.

In fact, our truncation method can be implemented as post-processing step that takes as an input the order-$p$ Taylor expansion $\Phi^\dagger$ of the energy produced by conventional homogenization, such as that in Equations (2.26) and (6.8). To derive the positive energy $\Phi_+$ for the laminate which we obtained in (2.34), for example, we could simply extract the non-positive energy functional $\Phi^\dagger$ in (2.26) from previous work [Le15] (carefully avoiding the alternate form (2.28) that involves an additional integration by parts) and apply the decomposition method directly.

Thanks to the iterative nature of the underlying Cholesky decomposition, the positive functional $\Phi_+$ obtained at order $p$ refines that obtained at a lower order $p' < p$, by adding higher-order terms without modifying the existing terms. By extending the Cholesky factorization done in Section 5 at order $p' = 2$ to order $p = 4$, for instance, one obtains an energy of the form

$$ \Phi = \frac{1}{2} \int \left( \hat{\epsilon}^{(0)} \cdot \mathbf{D}_0 \cdot \hat{\epsilon}^{(0)} + \eta^2 \hat{\epsilon}^{(1)} \cdot \underline{\underline{\mathbf{L}}}_1 \cdot \hat{\epsilon}^{(1)} + \eta^4 \hat{\epsilon}^{(2)} \cdot \underline{\underline{\underline{\mathbf{L}}}}_2 \cdot \hat{\epsilon}^{(2)} + \dots \right) d\mathbf{x} \quad (9.1) $$

where the generalized strains $\hat{\epsilon}^{(k)}$ are given by

$$ \begin{align} \hat{\epsilon}^{(0)} &= \epsilon^{(0)} + \eta \underline{\underline{\mathbf{L}}}_1^T \cdot \epsilon^{(1)} + \eta^2 \underline{\underline{\underline{\mathbf{L}}}}_2^T \cdot \epsilon^{(2)} + \eta^3 \underline{\underline{\underline{\underline{\mathbf{L}}}}}_3^T \cdot \epsilon^{(3)} + \eta^4 \underline{\underline{\underline{\underline{\underline{\mathbf{L}}}}}}_4^T \cdot \epsilon^{(4)} + \dots \nonumber \\ \hat{\epsilon}^{(1)} &= \epsilon^{(1)} + \eta \underline{\underline{\underline{\tilde{\mathbf{L}}}}}_2^T \cdot \epsilon^{(2)} + \eta^2 \underline{\underline{\underline{\tilde{\mathbf{L}}}}}^T_3 \cdot \epsilon^{(3)} + \dots \tag{9.2} \\ \hat{\epsilon}^{(2)} &= \epsilon^{(2)} + \dots \nonumber \end{align} $$

and $\epsilon^{(k)}$ represents the components of the tensor $\nabla^k \epsilon$ flattened into a vector. The homogenized constants having 2 underlines or less are identical to those obtained at order $p' = 2$, see Equations (5.4–5.5). Increasing the homogenization order to $p=4$ brings in the quantities bearing 3 or 4 underlines, namely a new, positive block of homogenized moduli $\underline{\underline{\mathbf{L}}}_2$ acting on the second gradient of strain $\hat{\epsilon}^{(2)}$, and higher-order coefficients $\underline{\underline{\mathbf{L}}}_3, \underline{\underline{\tilde{\mathbf{L}}}}_2, \underline{\underline{\tilde{\mathbf{L}}}}_4$ and $\underline{\underline{\tilde{\mathbf{L}}}}_3$ refining the generalized strain measures $\hat{\epsilon}^{(0)}$ and $\hat{\epsilon}^{(1)}$ that were already present at order 2.
With $\underline{L}_1 \neq \underline{\tilde{L}}_2$ and $\underline{L}_2 \neq \underline{\tilde{L}}_3$ in general, the second generalized strain $\hat{\varepsilon}^{(1)}$ appearing in (9.2) is not simply the gradient of the first generalized strain $\hat{\varepsilon}^{(0)}$.

The Cholesky-based truncation naturally handles non-centro-symmetric structures, see Section 6, continuous as well as discrete structures, arbitrary space dimension, and arbitrary homogenization orders $p$, including odd ones—the case $p=3$, for instance, corresponds to discarding the quantities $\underline{D}_2$, $\underline{L}_4$ and $\underline{\tilde{L}}_3$ in (9.1–9.2).

The truncation method restores positiveness by introducing a higher-order non-asymptotic term, which in the case of the laminate and at homogenization order $p=2$ takes the form $\frac{\eta^4}{2} \nabla^2 \varepsilon \cdot [C^T \cdot K^{-1} \cdot C] \cdot \nabla^2 \varepsilon$, see (2.34). As a result, the order of the boundary-value problem for the displacement is 6, which is higher than the order (4) of the boundary-value problem produced by the naïve truncation method. In general, when homogenizing at an even order $p$, the boundary-value problem associated with the positive energy $\Phi_+$ is of order 2 ($p+1$), which is 2 more than that produced by the naïve truncation method, namely 2 $p$. With this relatively high differentiation order, the numerical solution of the boundary-value problem produced by the positive functional $\Phi_+$ may require special attention.

In hindsight, at least two instances of positive energy functionals of the type $\Phi_+$ have appeared in earlier work on second-order dimension reduction:

* The one-dimensional energy governing the twisting of linearly elastic beams has been obtained in Equation [5.17] of [AL21] in a form which we can rewrite as

$$ \Phi_+ = \int_0^L \frac{1}{2} \mu J \left( \{\tau(x) + \eta^2 L \tau''(x)\}^2 + \eta^2 b \tau'^2(x) \right) dx, \quad (9.3) $$

where $\mu J > 0$ is the classical twisting modulus, $\eta$ is the cross-section radius, $\{\tau + \dots\}$ is a generalized twisting strain that corrects the classical twisting strain $\tau(x) = \theta'(x)$ by a higher-order term, $\theta(x)$ is the rotation of the cross-section about the axis, $L$ is a constant depending on the cross-section properties, and $b$ is a higher-order modulus. The reasoning done in Section 5.2 of the present work proves that the higher-order modulus $b > 0$ is positive, as long as the underlying elastic material is stable. The positivity of $b$ had only been noted in the special case of elliptical cross-section in the original paper [AL21].

* The one-dimensional energy governing the stretching of an elastic cylinder subjected to finite deformations has been initially derived in [AH16] and later reformulated in Appendix B.1 of [AH19] in the form

$$ \Phi_+ = \int_0^L W_0 \left( \{\lambda(X) + \eta^2 L(\lambda(X)) \chi^2(X)\} \right) dX \quad (9.4) $$

where $\lambda(X) = dx/dX$ is the effective longitudinal stretch, $X$ and $x$ are the longitudinal coordinates in reference and deformed configurations, respectively, $\eta$ is the initial radius of the cylinder, $W_0(\lambda)$ is the energy per unit length when the cylinder is subjected to uniform simple traction, the coefficient $L$ was found as $L(\lambda) = 1/(16\lambda^4)$. This energy is similar to that obtained in (3.20) for the triangular truss, when the quadratic potential $\frac{1}{2}\varepsilon:K:\varepsilon$ applicable to linear elasticity is identified with the non-quadratic potential $W_0(\lambda)$ relevant to finite-strain elasticity. By contrast with the linear case, however, the generalized strain $\hat{\lambda} = \{\lambda + \eta^2 \chi^2/(16\lambda^4)\}$ used as an argument of $W_0$ in (9.4) includes a correction proportional to $(\eta\chi')^2$ and not $\eta^2\chi'$. The functional (9.4) will typically yield stable solution whenever the underlying one-dimensional potential $W_0(\lambda)$ is convex.

The energy functional (9.4) relevant to finite elasticity suggests that the proposed reduction method can be extended to nonlinear problems in homogenization and dimension reduction.

In the nonlinear analysis of the stretched cylinder [AL21], it has been noted that the generalized one-dimensional stretch $\hat{\lambda}(X) = \{\lambda(X) + \eta^2 \chi^2(X)/(16\lambda^4(X))\}$ has a simple physical interpretation, unlike $\lambda(X)$: it is the average over the cross-section of the microscopic nonlinear measure of longitudinal strain. This suggests that the generalized strain measures produced by the Cholesky procedure, such as $\{\varepsilon + \eta + \dots + \eta^2\dots\}$, are more meaningful than the standard strain measures, such as $\varepsilon = (\nabla u + \nabla u^T)/2$, used in previous work on higher-order homogenization. We hope that additional useful mathematical properties of the homogenized energy concerning, e.g., the regularity of the solution, can be identified by working with this alternate strain measure.

Having restored the crucial property of positivity, our method could help increasing the popularity of higher-order homogenization. There remains a last obstacle towards this goal, which will hopefully be clarified in future work: we have limited attention to infinite systems in this paper and it is yet unclear how one can account for boundary conditions in the energy approach.

Acknowledgements: we would like to thank Hussein Nassar for pointing out the connection of our work with the Cholesky method, and Antoine Gloria for insightful comments.

Bibliography

[ABB09] N. Auffray, R. Bouchet, and Y. Bréchet. Class-jump phenomenon for physical symmetries in bi-dimensional space. In *Mechanics of Microstructured Solids: Cellular Materials, Fibre Reinforced Solids and Soft Tissues*, volume 46, page 1. Springer, 2009.

[ABV16] G. Allaire, M. Briane, and M. Vanninathan. A comparison between two-scale asymptotic expansions and Bloch wave expansions for the homogenization of periodic structures. *SeMA Journal*, 73(3):237–259, sep 2016.
[AH16] B. Audoly and J. W. Hutchinson. Analysis of necking based on a one-dimensional model. *Journal of the Mechanics and Physics of Solids*, 97:68–91, dec 2016.

[AH19] B. Audoly and J. W. Hutchinson. One-dimensional modeling of necking in rate-dependent materials. *Journal of the Mechanics and Physics of Solids*, 123:149–171, 2019.

[AL21] B. Audoly and C. Lestringant. Asymptotic derivation of high-order rod models from non-linear 3D elasticity. *Journal of the Mechanics and Physics of Solids*, 148:104264, mar 2021.

[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.

[AS77] J. L. Auriault and E. Sanchez-Palencia. Etude du comportement macroscopique d'un milieu poreux saturé déformable. *J Mec*, 16:575–603, 1977.

[Aud23] B. Audoly. The shoal library. https://archive.softwareheritage.org/browse/origin/https://git.renater.fr/anonscm/git/shoal/shoal.git. 2023.

[BGP17] C. Boutin, F. dell'Isola, I. Giorgio, and L. Placidi. Linear pantographic sheets: Asymptotic micro-macro models identification. *Mathematics and Mechanics of Complex Systems*, 5(2):127–162, 2017.

[BLP11] A. Bensoussan, J. L. Lions, and G. Papanicolaou. Asymptotic Analysis for Periodic Structures, volume 374. American Mathematical Soc., 2011.

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

[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.

[DLSS22] B. Durand, A. Lebée, P. Seppecher, and K. Sab. Predictive strain-gradient homogenization of a pantographic material with compliant junctions. *Journal of the Mechanics and Physics of Solids*, 160:104773, 2022.

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

[Dur22] B. Durand. Modélisation Numérique et Caractérisation Expérimentale de Matériaux Architecturés Au Comportement Élastique de Second Gradient. These de doctorat, Marne-la-vallée, ENPC, 2022.

[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.

[KA24] A. Kumar and B. Audoly. Wavelength selection in the twist buckling of pre-strained elastic ribbons. 2024.

[KC02] V. G. Kac and P. Cheung. *Quantum calculus*, volume 113. Springer, 2002.

[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.

[Le15] D. T. Le. Modèle d'endommagement à Gradient : Approche Par Homogénéisation. These de doctorat, Paris 6, 2015.

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

[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.

[S18a] H. Abdoul-Anziz and P. Seppecher. Homogenization of periodic graph-based elastic structures. *Journal de l'École polytechnique Mathématiques*, 5:259–288, 2018.

[S18b] 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.

[S18c] 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.

[SAI11] P. Seppecher, J. J. Alibert, and F. Isola. Linear elastic trusses leading to continua with exotic mechanical interactions. *Journal of Physics: Conference Series*, 319(1):12018, 2011.

[San80] E. Sanchez-Palencia. Non-homogeneous media and vibration theory. Volume 127 of Lecture Notes in Physics. Springer, Berlin, Heidelberg, 1980.

[SC00] V. P. Smyshlyayev and K. D. Cherednichenko. On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media p. Page 33, 2000.

[Str07] G. Strang. *Computational Science and Engineering*. Wellesley-Cambridge press, 2007.

[TAL24a] 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.

[TAL24b] M. Thbaut, B. Audoly, and C. Lestringant. Supplementary material to 'Fixing non-positive energies in higher-order homogenization': mathematica notebooks implementing homogenization of the lattice examples. https://archive.softwareheritage.org/swh:1:dir:198a2af51d9e4d0215192f801680fb68fee2da67;origin=https://git.renater.fr/anonscm/git/shoal/shoal.git;visit=swh:1:snp:a6206a25ac2debd558d139910aa4344eb03289be;anchor=swh:1:rev:7b07086583aefa04223db3bf79dc09bbb32c9c74;shoal-library-v1.0/elastic_lattices/illustrations/2025-JMPS-paper-Thbaut-et-al/, 2024.

[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, nov 1996.

[VP12] A. Vigliotti and D. Pasini. Linear multiscale analysis and finite element validation of stretching and bending dominated lattice materials. *Journal of Materials*, 46:57–68, 2012.

[Wol21] Wolfram. Mathematica, Version 13.2. https://www.wolfram.com/mathematica, 2021. Champaign, IL, 2021.

[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 Notation for multiple contractions

Given two tensors $a$ and $b$ of respective order $p$ and $p'$, the contraction operations are defined as

$$ \begin{align} (a \cdot b)_{i_1 \dots i_p} & = a_{i_1 \dots i_p} b_{kli_5 \dots i'_p} , \nonumber \\ (a \circ b)_{i_1 \dots i_p} & = a_{i_1 \dots i_p} b_{klmi_4 \dots i'_p} , \tag{A.1} \\ (a \circ \circ b)_{i_1 \dots i_p} & = a_{i_1 \dots i_p} b_{klmn} b_{klmni_5 \dots i'_p} . \nonumber \end{align} $$