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

HAL
open science

# An energy approach to asymptotic, higher-order, linear homogenization

Basile Audoly, Claire Lestringant

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

## To cite this version:

Basile Audoly, Claire Lestringant. An energy approach to asymptotic, higher-order, linear homogenization. Journal of Theoretical, Computational and Applied Mechanics, inPress, 10.46298/jt-cam.11414. hal-04112136v4

HAL Id: hal-04112136

https://hal.science/hal-04112136v4

Submitted on 18 Dec 2023

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.11414

HAL hal-04112136v4

**History**

Received Jun 01, 2023

Accepted Sep 26, 2023

Published Dec 18, 2023

**Associate Editor**

Laurence BRASSART

**Reviewers**

Pierre SEPPECHER

Damiano PASINI

Anonymous

Jean-François GANGHOFFER

**Open Review**

HAL hal-04232718v1

**Supplementary Material**

Software
SWHID [SoftwareHeritage]

**Licence**

CC BY 4.0

©The Authors

# An energy approach to asymptotic, higher-order, linear homogenization

Basile AUDOLY<fnref n="1" /> and Claire LESTRINGANT<fnref n="2" />

A higher-order homogenization method for linear elastic structures is proposed. While most existing approaches to homogenization start from the equations of equilibrium, this one works at the energy level. We start from an energy functional depending on microscopic degrees of freedom on the one hand and on macroscopic variables on the other hand; the homogenized energy functional is derived by relaxing the microscopic degrees of freedom and applying a formal two-scale expansion. This method delivers the energy functional of the homogenized model directly, including boundary terms that have not been discussed in previous work. Our method is formulated in a generic setting which makes it applicable to a variety of geometries in dimension 1, 2 or 3. An implementation using a symbolic calculation language is proposed and it is distributed as an open-source library. Simple illustrations to elastic trusses having pre-stress or graded elastic properties are presented. The approach is presented in the context of discrete elastic structures and the connection with previous work on the higher-order homogenization of periodic continua is discussed.

**Keywords:** asymptotic homogenization, discrete elasticity, strain-gradient models

## 1 Introduction

This work addresses the growing need for effective models capturing accurately the mechanical response of architected materials produced, e.g., by additive manufacturing techniques. In these materials, finite-size effects are not captured by standard (Cauchy) continuum models. Generalized continua theories, including higher-order (gradient) terms or micropolar fields have been introduced to overcome this limitation. Among various strategies proposed to obtain such models, asymptotic analysis offers a rigorous and fully predictive tool to derive exact higher-order or micropolar contributions without resorting to ad-hoc kinematic assumptions (Boutin 2019).

Asymptotic homogenization is a well-established technique for either discrete or continuous periodic microstructures. It aims at deriving an equivalent macroscopic set of equilibrium equations by means of a formal two-scale expansion, see (Sanchez-Palencia 1980; Bakhvalov and Panasenko 1989; Cioranescu and Donato 1999; Cioranescu and Paulin 1999) among others. The results of classical periodic homogenization are recovered at leading order. Later contributions focussed on pushing the asymptotic expansion to higher-orders, with the aim of deriving gradient contributions, either as small correctors to the leading-order prediction (Gambin and Kröner 1989; Boutin 1996; Smyshlyaev and Cherednichenko 2000; Hans and Boutin 2008; Bacigalupo 2014; Le and Marigo 2018; Abali and Barchiesi 2021) or as a non-local leading-order model when the standard equivalent medium is degenerate (Boutin and Soubestre 2011; Abdoul-Anziz and Seppecher 2018a; Abdoul-Anziz and Seppecher 2018b; Abdoul-Anziz et al. 2019; Durand et al. 2022). Asymptotic analysis has also been used to derive micropolar effective theories (Bažant and Christensen 1972; Dos Reis and Ganghoffer 2012; Nassar et al. 2020).

In spite of a large body of existing work, some aspects of higher-order homogenization remain elusive. Besides the question of applicable boundary conditions which is largely overlooked in existing work, most contributions typically start from the strong form of the equilibrium

<footnote n="1">¹ Laboratoire de Mécanique des Solides, CNRS, Institut Polytechnique de Paris, 91120 Palaiseau, France</footnote>

<footnote n="2">² Institut Jean Le Rond d'Alembert, Sorbonne Université, CNRS, 75005 Paris, France</footnote>
equations thus losing the variational structure of the initial problem (some try to reconstruct the effective energy a posteriori, but systematically ignore boundary terms). Besides, natural extensions such as pre-stress or slow spatial variations of the elastic or geometric properties of the microstructure are rarely addressed, especially in higher-order contributions.

In this paper, we adapt our earlier work on higher-order, nonlinear asymptotic dimension reduction (Lestringant and Audoly 2020) to linear homogenization of discrete elastic structures, such as elastic trusses or networks of elastic beams. Based on formal arguments, we propose a homogenization method that is directly applicable to any given discrete microstructure (the connection with the homogenization of periodic continua is discussed in Section 7. We carry out homogenization at the energy level, thereby following the pioneering work of (Le and Marigo 2018). As a benefit, important simplifications are made during the homogenization procedure, allowing the higher-order homogenized model to be obtained in compact form via a formal, systematic expansion. In addition, our method is versatile and allows for natural extensions such as the case of pre-stressed structures. We also keep track of boundary terms that have been ignored in previous work. Lastly, our method is implemented in a symbolic calculation language and distributed as an open-source library named `shoal`, for Second-order HOMogenization Automated in a Library (Audoly 2023).

Our method starts from an energy formulation derived by applying a continuum approximation to a *linear* discrete microstructure. This preliminary continualization step is not strictly part of this paper; it is briefly illustrated based on an example. We apply to this continualized energy a two-scale expansion, assuming slow variations for macroscopic quantities such as strains and lattice properties, and rapid variations for the degrees of freedom at the microstructural scale. Our method aims at condensing out these rapidly varying fields, by solving an energy stationarity problem order by order. Constraints are included from the onset in the energy formulation, which helps making the derivation compact. The procedure involves solving a set of elementary linear-algebra problems, and is implemented in a formal calculation language.

## 2 Input to the homogenization procedure

In this section, we present the elastic model serving as the starting point of the homogenization procedure. It is formulated in a continuous domain $\Omega \subset \mathbb{R}^d$ of the $d$-dimensional Euclidean space. The elastic model is specified in an abstract and generic form, making the homogenization procedure applicable to a broad range of situations including both periodic continua and discrete structures—such as one-dimensional (beam-like), two-dimensional (plate-like) and three-dimensional (bulk) elastic trusses. It is also versatile enough to handle structures possessing slowly modulated properties or pre-stress, as illustrated in Sections 4.1 and 4.2, respectively.

The task of casting a given problem into the generic form proposed in this section is not particularly difficult but has to be carried out on a case-by-case basis: this aspect is barely touched upon in Section 4 and will be illustrated in follow-up papers. When the structure under study is a discrete elastic truss, this preliminary step involves a so-called *continualization* assumption that turns the original discrete energy into the *continuous* energy functional used as a starting point in this paper, see Equation (4) and the discussion in Section 4.1.

We limit attention to *linear* homogenization problems. The extension to nonlinear homogenization can be done by adapting our previous work on nonlinear dimension reduction (Lestringant and Audoly 2020; Audoly and Lestringant 2021).

### 2.1 Energy formulation of the input model

We proceed to specify the continuous elastic model used as an input to the homogenization procedure. The presentation is intentionally abstract: illustrations will be provided in Section 4.

The model is formulated over a continuous domain $\Omega \subset \mathbb{R}^d$, and we denote by $X \in \Omega$ the space variable, see Figure 1. A deformed configuration of the elastic body is parameterized by three vector fields $y(X)$, $l(X)$ and $m(X)$ defined over $\Omega$:

* *microscopic degrees of freedom* $y(X) \in \mathbb{R}^{n_y}$ which we seek to eliminate using the homogenization procedure,
<figure><img src="image_4.png" /><figcaption>**Figure 1** Typical application of the method: homogenization of a beam lattice in a domain $\Omega \subset \mathbb{R}^d$, here with $d = 2$. In this stylized representation of the beam lattice, the gradient of material parameters $\nabla m$ includes both gradients of elastic properties (varying thickness of beams) and geometry (mesh curvature). The mathematical domain $\Omega$ where we carry out homogenization is strictly included in the physical domain of the lattice to avoid boundary layers. The elastic lattice depicted above serves as an illustration: the homogenization method is presented in an abstract setting that does not make any reference to a lattice.</figcaption></figure>

* macroscopic variables which are held fixed during homogenization, namely:

- the macroscopic strain $l(X) \in \mathbb{R}^{n_l}$

- variable material parameters $m(X) \in \mathbb{R}^{n_m}$.

The integers $n_y, n_l$ and $n_m$ are input parameters of the homogenization procedure. The goal of the procedure is to slave the microscopic degrees of freedom $y$ to the macroscopic variables $l$ and $m$, thereby delivering a homogenized model depending on $l$ and $m$ only. The difference between the macroscopic variables $l$ and $m$ is that $m$ captures the slowly variable properties of the elastic structure which are prescribed once for all, although $l$ is considered fixed during the homogenization procedure but is actually an unknown of the structural problem that the homogenized energy helps solving.

In addition, the input model makes use of microscopic strain variables, which are collected into a vector $E \in \mathbb{R}^{n_E}$. The geometric definition of the strain $E$ is assumed to be of the form

$$E = E(m(X); l(X), \nabla l(X), \nabla^2 l(X), \dots; y(X), \nabla y(X), \nabla^2 y(X), \dots) \quad (1)$$

where the dependence on the macroscopic strain $l$ and of the microscopic degrees of freedom $y$ and their gradients is linear in the context of linear elasticity,

$$E(m; l, l', l'', \dots; y, y', y'', \dots) = E_l(m) \cdot l + E'_l(m) : l' + E''_l(m) : l'' + \dots \\ + E_y(m) \cdot y + E'_y(m) : y' + E''_y(m) : y'' + \dots \quad (2)$$

Here, $l \in \mathbb{R}^{n_l}$, $l' = \nabla l \in \mathbb{T}^{(n_l, d)}$, $l'' = \nabla^2 l \in \mathbb{T}^{(n_l, d, d)}$, $y \in \mathbb{R}^{n_y}$, $y' = \nabla y \in \mathbb{T}^{(n_y, d)}$, $y'' = \nabla^2 y \in \mathbb{T}^{(n_y, d, d)}$ are dummy variables representing the local values of $l$, $y$ and their successive gradients. The homogenization procedure is implemented in a symbolic calculation language, and the tensors $E_l(m) \in \mathbb{T}^{(n_E, n_l)}$, $E'_l(m) \in \mathbb{T}^{(n_E, n_l, d)}$, $E''_l(m) \in \mathbb{T}^{(n_E, n_l, d, d)}$, ..., $E_y(m) \in \mathbb{T}^{(n_E, n_y)}$, $E'_y(m) \in \mathbb{T}^{(n_E, n_y, d)}$ and $E''_y(m) \in \mathbb{T}^{(n_E, n_y, d, d)}$ are provided as input, see Table 1, in the form of symbolic, tensor-valued functions of $m$. Expressions for these tensors that are applicable to specific examples are provided in Section 4, see Table 4. Our tensor notations such as $\mathbb{T}^{(n_l, d, d)}$ for tensor spaces are given in Appendix A.

**Remark 1** We work in the discrete case, i.e., assume a finite number ($n_y < \infty$) of microscopic degrees of freedom. This is typically relevant to elastic truss structures, see Section 4. The extension of the method to infinitely many microscopic degrees of freedom ($n_y = \infty$), relevant to periodic continua, is discussed in Section 7.

**Remark 2** We will present the homogenization method without providing definitions of the macroscopic strain $l$, variable material parameters $m$ and microscopic degrees of freedom $y$. These quantities
can be anything as long as they match the postulated forms of the strain, constraints and energy in Equations (2) to (5). This abstract presentation makes the method quite versatile, and opens up the way for extensions that we will cover in future work. Specific choices of **l**, **m** and **y** are proposed in the illustration examples (Section 4).

The model used on input may make use of constraints and we denote by $n_c \ge 0$ the number of applicable (scalar) kinematic constraints. By convention, the left-hand sides of these kinematic constraints are incorporated into the microscopic strain vector $E \in \mathbb{R}^{n_E}$, and are extracted from $E$ using an appropriate matrix $Q \in \mathbb{T}^{(n_c, n_E)}$. In view of this, the kinematic constraints are written in the form

$$Q \cdot E(m(X); l(X), \nabla l(X), \dots; y(X), \nabla y(X), \dots) = 0_{n_c} \quad \forall X. \qquad (3)$$

Section 4 provides a specific illustration on how $E$ and $Q$ can be set up to conform to Equation (3). The constant tensor $Q$ is provided as an input of the homogenization procedure, see Table 1. In the input model, the strain energy is assumed to be of the form

$$\Phi[m, l; y] = \int_{\Omega} W(m(X), E(m(X); l(X), \nabla l(X), \dots; y(X), \nabla y(X), \dots))dX, \qquad (4)$$

where the strain energy density in the bulk term is given in the context of linear elasticity by

$$W(m, E) = \frac{1}{2} E \cdot \mathcal{K}(m) \cdot E. \qquad (5)$$

Explicit expressions of the strain $E$, energy $\Phi$ and elasticity matrix $\mathcal{K}$ appearing in (4–5) will be provided in the illustration Section 4, see for example (33–34) The elasticity tensor $\mathcal{K}(m) \in \mathbb{T}^{(n_E, n_E)}$ is provided as an input to the homogenization method in the form of a tensor-valued, symbolic function of $m$, see Table 1.

The square brackets around the arguments of $\Phi[m, l; y]$ in the left-hand side of Equation (4) denote a *functional* dependence: the strain energy $\Phi$ depends on the *functions* $m, l$ and $y$ over the entire domain $\Omega$.

The list of parameters passed as an input to the homogenization procedure is recapitulated in Table 1.

<table><caption>**Table 1** List of parameters passed as an input to the homogenization procedure. The notation used in the last two columns is defined in Appendix A.</caption><thead><tr><th></th><th>nature</th><th>tensor space</th><th>symmetries</th></tr></thead><tbody><tr><td>$d \ge 1$</td><td>space dimension</td><td></td><td></td></tr><tr><td>$n_m \ge 0$</td><td>number of material parameters</td><td></td><td></td></tr><tr><td>$n_l \ge 0$</td><td>number of macroscopic degrees of freedom</td><td></td><td></td></tr><tr><td>$n_y \ge 0$</td><td>number of microscopic degrees of freedom</td><td></td><td></td></tr><tr><td>$n_E \ge n_c$</td><td>number of strain variables</td><td></td><td></td></tr><tr><td>$n_c$</td><td>number of kinematic constraints</td><td></td><td></td></tr><tr><td>$E_l(m)$</td><td>dependence of strain on $l$</td><td>$\mathbb{T}^{(n_E, n_l)}$</td><td></td></tr><tr><td>$E'_l(m)$</td><td>dependence of strain on $\nabla l$</td><td>$\mathbb{T}^{(n_E, n_l, d)}$</td><td></td></tr><tr><td>$E''_l(m)$</td><td>dependence of strain on $\nabla^2 l$</td><td>$\mathbb{T}^{(n_E, n_l, d, d)}$</td><td>S<sub>34</sub></td></tr><tr><td>$E_y(m)$</td><td>dependence of strain on $y$</td><td>$\mathbb{T}^{(n_E, n_y)}$</td><td></td></tr><tr><td>$E'_y(m)$</td><td>dependence of strain on $\nabla y$</td><td>$\mathbb{T}^{(n_E, n_y, d)}$</td><td></td></tr><tr><td>$E''_y(m)$</td><td>dependence of strain on $\nabla^2 y$</td><td>$\mathbb{T}^{(n_E, n_y, d, d)}$</td><td>S<sub>34</sub></td></tr><tr><td>$Q$</td><td>constraint extraction, see Equation (3)</td><td>$\mathbb{T}^{(n_c, n_E)}$</td><td></td></tr><tr><td>$\mathcal{K}(m)$</td><td>stiffness matrix, see Equation (5)</td><td>$\mathbb{T}^{(n_c, n_E)}$</td><td>S<sub>12</sub></td></tr></tbody></table>

**Remark 3** This formulation of the input model is designed to be versatile. For instance, the presence of pre-strain can be accommodated by adding a coefficient capturing the pre-strain intensity as an
additional entry in the vector $l$ (whose definition is up to the user), and by propagating it to $E$ by an appropriate definition of $E_l(m)$, as illustrated in Section 4.2. Similarly, the presence of pre-stress can be accommodated by including a constant entry with value 1 in $l$, propagating it to $E$, and inserting the pre-stress into the corresponding rows and column of $\mathcal{K}$.

## 2.2 Assumption of slow variations

One of the key assumptions of homogenization is that there is a separation of scales between a microscopic length $\ell$ (typically the spatial period of the underlying discrete lattice or periodic continuum) and the size $L$ of the domain, $\ell \ll L$, as sketched in Figure 1. The goal of homogenization is to deliver an effective model applicable at the macroscopic scale $L$, by hiding the ‘details’ taking place at the microscopic scale $\ell$.

Mathematically, this separation of scale is captured by the small parameter

$$ \eta = \frac{\ell}{L} \ll 1. \qquad (6) $$

The various fields $f(X)$, such as $f = m$, $f = l$ or $f = y$, are assumed to depend on the variable $X$ evolving on the slow scale $L = \ell/\eta$, implying that their successive gradients scale as

$$ \nabla^k f(X) = O(\eta^k). \qquad (7) $$

In the following, the gradient $\nabla = \partial/\partial X$ will therefore be treated implicitly as a small quantity of order $\eta$. This implicit notation has the benefit of avoiding a large number of predictable occurrences of the parameter $\eta$, as discussed in Remark 4 below.

The scaling assumption Equation (7) is not applicable in the layers that are present near the boundaries or near the point of application of point-like force. The homogenization domain $\Omega$ therefore needs to be slightly smaller than the actual physical domain of the elastic body, as shown in Figure 1 (see also Equation [1] in the work of (Abdoul-Anziz and Seppecher 2018b) for an accurate description of how $\Omega$ can be shrunk). Alternatively, boundary layers can be solved rigorously and represented in the homogenized model by means of effective boundary terms, see for example (David et al. 2012), but this is beyond the scope of the present work.

**Remark 4** The scaling assumption Equation (7) can be motivated as follows. By convention, we consider the microscopic length $\ell$ to be $\ell = O(1)$ and the macroscopic length to be $L = O(\eta^{-1})$. In our notation, any macroscopic field $f$ such as $m$, $l$ or $y$ is implicitly a function of the slow variable $\bar{X} = \eta X$, i.e., what we write as $f(X)$ should be spelled out as $f(X) = \bar{f}(\eta X)$, where $\bar{f}$ is a dimensionless function, independent of $\eta$. The gradients can then be obtained as

$$ \nabla^k f(X) = \frac{\partial^k \bar{f}}{\partial X^k}(X) = \eta^k \frac{\partial^k \bar{f}}{\partial \bar{X}^k}(\eta X) = \eta^k \bar{\nabla}^k \bar{f}(\eta X), \qquad (8) $$

where $\bar{\nabla}^k \bar{f} = O(1)$ denotes the gradient with respect to the slow variable: this implies the scaling assumption $\nabla^k f(X) = O(\eta^k)$ in Equation (7). The formal rule (7) dispenses with a notation for the slow variable.

**Remark 5** The dependence on $m(X)$ of the strain $E$ in Equations (1) to (2) and of the energy density $W$ in Equations (4) to (5) allows one to handle the case of structures whose elastic properties vary over the large scale $L = \ell/\eta$, as conveyed by the variations in thickness of the microstructure sketched in Figure 1. The definition of $m(X)$ is entirely up to the user. For a 2D elastic truss possessing rotational symmetry, for instance, one could define $m(X) = (X_1^2 + X_2^2)^{1/2}$ and $n_m = 1$ to capture the dependence of the local truss properties on the distance to the center of symmetry. In the case of variable properties without any particular symmetry, one should set $m(X) = X$ and $n_m = d$. When specifying the input model, one should ensure that any dependence of the properties of the elastic medium on the slow variable $X$ takes place through the quantity $m(X)$, as illustrated in Section 4. For structures having uniform properties over the large scale, one can ignore any dependence on $m$ and set $n_m = 0$, see Appendix D.
# 3 Summary of the main results

## 3.1 Homogenization as a partial energy relaxation

In this paper, we identify an equivalent continuum by making stationary the functional $\Phi$ in Equation (4) over the microscopic degrees of freedom $\mathbf{y}(X)$ for a fixed distribution of the macroscopic variables $\mathbf{m}(X)$ and $\mathbf{l}(X)$. The stationary point is denoted as $\mathbf{y}^*[\mathbf{m}, \mathbf{l}]$, and will be assumed to be unique—as indicated by the square brackets, the stationary point $\mathbf{y}^*$ is a functional of $\mathbf{m}$ and $\mathbf{l}$; it is also a function of $X$, whose values are denoted as $\mathbf{y}^*[\mathbf{m}, \mathbf{l}](X)$.

In view of (2-3), the kinematic constraint can be written as $Q \cdot (E_y(\mathbf{m}(X)) \cdot \mathbf{y}(X) + E'_y(\mathbf{m}(X)) : \nabla \mathbf{y}(X) + \dots + E'_l(\mathbf{m}(X)) : \nabla \mathbf{l}(X) + \dots) = 0 \quad \forall X$. It is treated using Lagrange multipliers $g(X) \in \mathbb{R}^{n_c}$. The variational problem takes the following form: for given $\mathbf{m}$ and $\mathbf{l}$, we seek the solution $(\mathbf{y}, g) = (\mathbf{y}^*[\mathbf{m}, \mathbf{l}], g^*[\mathbf{m}, \mathbf{l}])$ of

$$ \left\{ \begin{array}{l} Q \cdot \left( \begin{array}{l} E_y(\mathbf{m}(X)) \cdot \mathbf{y}(X) + E'_y(\mathbf{m}(X)) : \nabla \mathbf{y}(X) + \dots \\ + E'_l(\mathbf{m}(X)) : \nabla \mathbf{l}(X) + \dots \end{array} \right) = 0, \quad \forall X \quad (9a) \\[1.5em] D_y \Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}; \delta\mathbf{y}] + \int_{\Omega} g(X) \cdot Q \cdot \left( \begin{array}{l} E_y(\mathbf{m}(X)) \cdot \delta\mathbf{y}(X) \\ + E'_y(\mathbf{m}(X)) : \nabla \delta\mathbf{y}(X) + \dots \end{array} \right) dX = 0, \quad \forall \delta\mathbf{y}. \quad (9b) \end{array} \right. $$

where $D_y\Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}; \delta\mathbf{y}] = \lim_{\tau \to 0} (\Phi[h, \mathbf{y} + \tau\delta\mathbf{y}] - \Phi_c[h, \mathbf{y}]) / \tau$ denotes the directional derivative and $\delta\mathbf{y}(X)$ is an arbitrary perturbation. Equation (9a) warrants that the stationary point $\mathbf{y} = \mathbf{y}^*[\mathbf{m}, \mathbf{l}]$ satisfies the kinematic constraint, while Equation (9b) warrants that it is indeed a stationary point among the $\mathbf{y}$'s satisfying the kinematic constraints—the integral term takes care of these constraints by multiplying the Lagrange multipliers by the incremental form of the constraint, as usual in the calculus of variations.

Having slaved the microscopic degrees of freedom $\mathbf{y} = \mathbf{y}^*[\mathbf{m}, \mathbf{l}]$ to the macroscopic variables, we can define a homogenized energy functional $\Phi^*[\mathbf{m}, \mathbf{l}]$ by inserting the stationary point $\mathbf{y}^*[\mathbf{m}, \mathbf{l}]$ into the original $\Phi$,

$$ \Phi^*[\mathbf{m}, \mathbf{l}] = \Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}^*[\mathbf{m}, \mathbf{l}]]. \quad (10) $$

The main goal of this paper is to derive an explicit expression of the homogenized functional $\Phi^*$.
The stationary point problem in (9) can be written formally as

$$ \mathbf{y}^*[\mathbf{m}, \mathbf{l}] = \underset{\mathbf{y} \text{ such that } Q \cdot E=0 \ \forall X}{\text{stpt}} \Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}]. \quad (11) $$

Equations (10) and (11) are at the heart of our variational approach to homogenization. They can be motivated as follows, by considering the broader structural problem of interest: a discrete truss, for instance, is governed by a total potential energy $\Psi[\mathbf{m}, \mathbf{l}] + \Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}]$, where $\Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}]$ is the strain energy of the truss and $\Psi[\mathbf{m}, \mathbf{l}]$ is the potential energy of the applied loading (under standard scaling assumptions, the latter does not depend on the microscopic degrees of freedom $\mathbf{y}$). Recalling that the variable elastic properties $\mathbf{m}(X)$ are fixed by design, the full elastic problem is solved by making the total potential energy $\Psi[\mathbf{m}, \mathbf{l}] + \Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}]$ stationary with respect to both the macroscopic unknowns $\mathbf{l}(X)$ and the microscopic ones $\mathbf{y}(X)$, subjected to the kinematic conditions $Q \cdot E = 0, \forall X$. Homogenization consists simply in enforcing the stationarity conditions sequentially, with respect to the microscopic degrees of freedom $\mathbf{y}$ first and to the macroscopic strain $\mathbf{l}$ next. Indeed, the stationarity condition of $\Psi[\mathbf{m}, \mathbf{l}] + \Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}]$ with respect to $\mathbf{y}$ is nothing but that considered in (11), given that $\Psi[\mathbf{m}, \mathbf{l}]$ does not depend on $\mathbf{y}$. Next, it can be checked that the stationarity condition with respect to $\mathbf{l}$ of $\Psi[\mathbf{m}, \mathbf{l}] + \Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}]$ is equivalent to the stationarity condition of the modified functional $\Psi[\mathbf{m}, \mathbf{l}] + \Phi^*[\mathbf{m}, \mathbf{l}]$ based on the homogenized strain energy $\Phi^*$ introduced in (10). This argument not only provides a justification to Equations (10) and (11), it also explains why the homogenized energy $\Phi^*[\mathbf{m}, \mathbf{l}]$ can be used as a substitute for the original energy $\Phi[\mathbf{m}, \mathbf{l}, \mathbf{y}]$ in the analysis of the structural problem.
The homogenization works under the assumption that the energy is positive-definite in the subspace of admissible microscopic degrees of freedom, i.e.,

$$ (\forall (\mathbf{m}, \mathbf{y} \neq 0) \text{ such that } \mathbf{Q} \cdot \mathbf{E}_y(\mathbf{m}) \cdot \mathbf{y} = 0) \quad (\mathbf{E}_y(\mathbf{m}) \cdot \mathbf{y}) \cdot \mathcal{K}(\mathbf{m}) \cdot (\mathbf{E}_y(\mathbf{m}) \cdot \mathbf{y}) > 0. \quad (12) $$

As we will see, this is a necessary condition for the variational problem (9) to have a unique solution at leading order. It is also a sufficient condition for the homogenization procedure to produce a result up to second-order.

## 3.2 Homogenization results in compact form

The variational problem (9) is impossible to solve in closed form in general but thanks to the assumption of scale separation (Section 2.2), we can derive the following approximation of $\Phi^*[\mathbf{m}, \mathbf{l}]$,

$$ \Phi^*[\mathbf{m}, \mathbf{l}] = \int_{\Omega} \left( K[\mathbf{m}] : \frac{l \otimes l}{2} + A[\mathbf{m}] : (l \otimes \nabla l) + B[\mathbf{m}] :: \frac{\nabla l \otimes \nabla l}{2} \right) dX + \\ \oint_{\partial\Omega} \left[ k[\mathbf{m}] :: \left( \frac{l \otimes l}{2} \otimes \mathbf{n} \right) + a[\mathbf{m}] :: (l \otimes \nabla l \otimes \mathbf{n}) \right] da + O(L^d \eta^3). \quad (13) $$

In the boundary terms $\oint_{\partial\Omega} \dots da$ in the second line, $\partial\Omega$ denotes the boundary of the domain, $\mathbf{n}$ is the unit outward normal, and $da$ the area (if $d=3$) or the length (if $d=2$) of a boundary element, see Figure 1. The typographical variant of the integral sign $\oint$ will be used throughout for boundary integrals.

The dimension and symmetries of the homogenized tensors $K$, $A$, $B$, $k$ and $a$ are specified in Table 2. Our main result is to derive their expansion in the successive gradients $\nabla^k \mathbf{m} = O(\eta^k)$,

$$ \begin{aligned} K[\mathbf{m}] &= K_0(\mathbf{m}) + K_1(\mathbf{m}) : \nabla \mathbf{m} + K_2(\mathbf{m}) :: (\nabla \mathbf{m} \otimes \nabla \mathbf{m}) + \dots \\ A[\mathbf{m}] &= A_0(\mathbf{m}) + A_1(\mathbf{m}) : \nabla \mathbf{m} + \dots \\ B[\mathbf{m}] &= B_0(\mathbf{m}) + \dots \\ k[\mathbf{m}] &= k_1(\mathbf{m}) : \nabla \mathbf{m} + \dots \\ a[\mathbf{m}] &= a_0(\mathbf{m}) + \dots \end{aligned} \quad (14) $$

The tensors $K_i(\mathbf{m})$, $A_i(\mathbf{m})$, $B_0(\mathbf{m})$, $k_1(\mathbf{m})$ and $a_0(\mathbf{m})$ appearing in the right-hand sides of (14) are obtained in explicit form in terms of the local material parameters $\mathbf{m}$ in Appendix B, see (B.12), and in Appendix C, see Appendices C.5 and C.10. Their properties are listed in Table 3.

This completes the specification of the homogenized model (13) up to order $\eta^2$ included.

<table><caption>**Table 2** Dimensions and symmetries of the tensors appearing in the homogenized energy functional in (13).</caption><thead><tr><th></th><th>tensor space</th><th>symmetry</th></tr></thead><tbody><tr><td><b>K</b>[<b>m</b>]</td><td>T<sup>(n<sub>l</sub>,n<sub>l</sub>)</sup></td><td>S<sub>12</sub></td></tr><tr><td><b>A</b>[<b>m</b>]</td><td>T<sup>(n<sub>l</sub>,n<sub>l</sub>,d)</sup></td><td>-</td></tr><tr><td><b>B</b>[<b>m</b>]</td><td>T<sup>(n<sub>l</sub>,d,n<sub>l</sub>,d)</sup></td><td>S<sub>{12}</sub>{34}</td></tr><tr><td><b>k</b>[<b>m</b>]</td><td>T<sup>(n<sub>l</sub>,n<sub>l</sub>,d)</sup></td><td>S<sub>12</sub></td></tr><tr><td><b>a</b>[<b>m</b>]</td><td>T<sup>(n<sub>l</sub>,n<sub>l</sub>,d,d)</sup></td><td>-</td></tr><tr><td><b>Y</b>[<b>m</b>]</td><td>T<sup>(n<sub>y</sub>,n<sub>l</sub>)</sup></td><td>-</td></tr><tr><td><b>Y'</b>[<b>m</b>]</td><td>T<sup>(n<sub>y</sub>,n<sub>l</sub>,d)</sup></td><td>-</td></tr></tbody></table>

The expansion (13-14) is established in Section 5 by solving the variational problem (9) for $\mathbf{y}$ order by order in $\eta$. The solution is found in the form

$$ \mathbf{y}^*[\mathbf{m}, \mathbf{l}](X) = Y[\mathbf{m}](X) \cdot l(X) + Y'[\mathbf{m}](X) : \nabla l(X) + O(\eta^2) \quad (15) $$

where $Y[\mathbf{m}]$ and $Y'[\mathbf{m}]$ are given as expansions in the successive gradients of $\mathbf{m}$,

$$ \begin{aligned} Y[\mathbf{m}](X) &= Y_0(\mathbf{m}(X)) + Y_1(\mathbf{m}(X)) : \nabla \mathbf{m}(X) + \dots \\ Y'[\mathbf{m}](X) &= Y'_0(\mathbf{m}) + \dots \end{aligned} \quad (16) $$
<table><caption>**Table 3** Tensors delivered by the homogenization procedure, defining the homogenized model in Equations (13) and (14).</caption><thead><tr><th></th><th>tensor space</th><th>symmetry</th><th>usage</th></tr></thead><tbody><tr><td>K<sub>0</sub>(m)</td><td>T(n<sub>l</sub>,n<sub>l</sub>)</td><td>S<sub>12</sub></td><td>Φ<sup>*</sup><sub>[0]</sub> = ∫<sub>Ω</sub> K<sub>0</sub> : <sup>l⊗l</sup>&frasl;<sub>2</sub> dX</td></tr><tr><td>K<sub>1</sub>(m)</td><td>T(n<sub>l</sub>,n<sub>l</sub>,n<sub>m</sub>,d)</td><td>S<sub>12</sub></td><td>Φ<sup>*</sup><sub>[1]</sub> = ∫<sub>Ω</sub> (K<sub>1</sub> : ∇m) : <sup>l⊗l</sup>&frasl;<sub>2</sub> dX + ...</td></tr><tr><td>K<sub>2</sub>(m)</td><td>T(n<sub>l</sub>,n<sub>l</sub>,n<sub>m</sub>,d,n<sub>m</sub>,d)</td><td>S<sub>12</sub>, S<sub>{34}</sub>{56}</td><td>Φ<sup>*</sup><sub>[2]</sub> = ∫<sub>Ω</sub> (K<sub>2</sub> :: (∇m ⊗ ∇m)) : <sup>l⊗l</sup>&frasl;<sub>2</sub> dX + ...</td></tr><tr><td>A<sub>0</sub>(m)</td><td>T(n<sub>l</sub>,n<sub>l</sub>,d)</td><td>-</td><td>Φ<sup>*</sup><sub>[1]</sub> = ∫<sub>Ω</sub> A<sub>0</sub> :: (I ⊗ ∇l) dX + ...</td></tr><tr><td>A<sub>1</sub>(m)</td><td>T(n<sub>l</sub>,n<sub>l</sub>,d,n<sub>m</sub>,d)</td><td>-</td><td>Φ<sup>*</sup><sub>[2]</sub> = ∫<sub>Ω</sub> (A<sub>1</sub> : ∇m) :: (I ⊗ ∇l) dX + ...</td></tr><tr><td>B<sub>0</sub>(m)</td><td>T(n<sub>l</sub>,d,n<sub>l</sub>,d)</td><td>S<sub>{12}</sub>{34}</td><td>Φ<sup>*</sup><sub>[2]</sub> = ∫<sub>Ω</sub> B<sub>0</sub> :: <sup>∇l⊗∇l</sup>&frasl;<sub>2</sub> dX + ...</td></tr><tr><td>k<sub>1</sub>(m)</td><td>T(n<sub>l</sub>,n<sub>l</sub>,d,n<sub>m</sub>,d)</td><td>S<sub>12</sub></td><td>Φ<sup>*</sup><sub>[2]</sub> = ∮<sub>∂Ω</sub> (k<sub>1</sub>[m] : ∇m) :: (<sup>l⊗l</sup>&frasl;<sub>2</sub> ⊗ n) da + ...</td></tr><tr><td>a<sub>0</sub>(m)</td><td>T(n<sub>l</sub>,n<sub>l</sub>,d,d)</td><td>-</td><td>Φ<sup>*</sup><sub>[2]</sub> = ∮<sub>∂Ω</sub> a<sub>0</sub>[m] :: (I ⊗ ∇l ⊗ n) da + ...</td></tr><tr><td>Y<sub>0</sub>(m)</td><td>T(n<sub>y</sub>,n<sub>l</sub>)</td><td>-</td><td>y<sup>*</sup><sub>[0]</sub> = Y<sub>0</sub>(m) ⋅ l</td></tr><tr><td>Y<sub>1</sub>(m)</td><td>T(n<sub>y</sub>,n<sub>m</sub>,d,n<sub>l</sub>)</td><td>-</td><td>y<sup>*</sup><sub>[1]</sub> = (Y<sub>1</sub>(m) : ∇m) ⋅ l + ...</td></tr><tr><td>Y'<sub>0</sub>(m)</td><td>T(n<sub>y</sub>,n<sub>l</sub>,d)</td><td>-</td><td>y<sup>*</sup><sub>[1]</sub> = Y'<sub>0</sub>(m) : ∇l + ...</td></tr><tr><td>G<sub>0</sub>(m)</td><td>T(n<sub>c</sub>,n<sub>l</sub>)</td><td>-</td><td>g<sup>*</sup><sub>[0]</sub> = G<sub>0</sub>(m) ⋅ l</td></tr></tbody></table>

and the localization tensors $Y_0(m)$, $Y_1(m)$ and $Y'_0(m)$ are derived in explicit form in terms of the variable material parameters $m$ in the Appendix, see (B.8) and (C.48).

### 3.3 Homogenization results in the form of a systematic expansion

The various contributions to $\Phi^*$ in (13) can be grouped order by order as follows, by inserting (14) into (13) and using (7):

* The leading-order contribution $\Phi^*_{[0]} = O(L^d\eta^0)$ is given by

$$ \Phi^*_{[0]}[m, l] = \int_{\Omega} K_0(m) : \frac{l \otimes l}{2} dX, \qquad (17) $$

and characterizes an equivalent Cauchy medium through a homogenized stiffness tensor $K_0(m)$ depending only on the local material parameters $m$: this homogenized stiffness $K_0(m)$ matches that predicted by classical homogenization.

* The first correction $\Phi^*_{[1]} = O(L^d\eta^1)$ is given by

$$ \Phi^*_{[1]}[m, l] = \int_{\Omega} \left( (K_1(m) : \nabla m) : \frac{l \otimes l}{2} + A_0(m) : (l \otimes \nabla l) \right) dX. \qquad (18) $$

* The second correction $\Phi^*_{[2]} = O(L^d\eta^2)$ is given by

$$ \Phi^*_{[2]}[m, l] = \int_{\Omega} \left( \begin{aligned} & (K_2(m) :: (\nabla m \otimes \nabla m)) : \frac{l \otimes l}{2} \\ & + (A_1(m) : \nabla m) : (l \otimes \nabla l) + B_0(m) :: \frac{\nabla l \otimes \nabla l}{2} \end{aligned} \right) dX \\ + \oint_{\partial\Omega} \left[ (k_1(m) : \nabla m) : \left( \frac{l \otimes l}{2} \otimes n \right) + a_0(m) :: (l \otimes \nabla l \otimes n) \right] da. \qquad (19) $$

The homogenized energy $\Phi^*[m, l]$ in (13) is nothing but the sum

$$ \Phi^*[m, l] = \Phi^*_{[0]} + \Phi^*_{[1]} + \Phi^*_{[2]} + O(L^d\eta^3), \qquad (20) $$

and it is asymptotically exact up to a higher-order contribution $\Phi^*_{[3]} = O(L^d\eta^3)$ which we do not attempt to resolve.

**Remark 6** The actual derivation of the homogenized model proceeds in the reverse order than the high-level presentation above: the order-by-order expansion (17–20) is derived first, and the compact form (13–14) is obtained next by rearranging the terms.
The solution for $\mathbf{y}$ in (15–16) is derived based on the assumption that the microscopic variables $\mathbf{y} = \mathbf{y}^*[\mathbf{m}, \mathbf{l}]$ can be expanded in powers of $\eta$,

$$ \mathbf{y}(X) = \mathbf{y}_{[0]}(X) + \mathbf{y}_{[1]}(X) + \mathbf{y}_{[2]}(X) + \dots, \qquad (21) $$

where $\mathbf{y}_{[k]}(X) = O(\eta^k)$ denotes the contribution of order $\eta^k$ to $\mathbf{y}$. Specifically, the microscopic solution $\mathbf{y}^* = \mathbf{y}^*[\mathbf{m}, \mathbf{l}] = \mathbf{y}_{[0]}^*(X) + \mathbf{y}_{[1]}^*(X) + \dots$ is derived order by order as

$$ \mathbf{y}_{[0]}^*(X) = Y_0(\mathbf{m}(X)) \cdot \mathbf{l}(X) \qquad (22a) $$

$$ \mathbf{y}_{[1]}^*(X) = (Y_1(\mathbf{m}(X)) : \nabla \mathbf{m}(X)) \cdot \mathbf{l}(X) + Y'_0(\mathbf{m}(X)) : \nabla \mathbf{l}(X), \qquad (22b) $$

which yields (15–16) by rearranging the terms.

**Remark 7** As discussed in Section 2.2, there are implicit scaling factors $\eta^k$ in all our formulas. Their consistency can be checked as follows. Take Equation (22b), for instance: the subscript ‘[1]’ in the left-hand side indicates that this is a quantity of order $\eta$; this is consistent with the fact that the right-hand side is homogeneous of degree 1 with respect to the symbol $\nabla = O(\eta)$. When checking homogeneity, the boundary terms must be treated with special care: in Equation (19), for instance, the integrand of the bulk integral is as quantity of order $\eta^2$, in line with the subscript ‘[2]’ appearing in the left-hand. The integrand of the boundary integral is however a quantity of order $\eta^1$; the paradox is resolved by noting that the measure of the domain $O(L^d)$ for the bulk integral, but $O(L^{d-1}) = O(L^d\eta)$ for the boundary integral—recall that $\ell = O(1)$ and $L = O(\eta^{-1})$. Ultimately, both integrals are of order $O(L^d\eta^2)$.

# 4 Illustrations

In this section, we provide simple illustrations of the homogenization method. Equivalent high-order beam models are derived for various truss lattices, in the same line as a number of earlier works on periodic 1D structures, including (Hans and Boutin 2008; Abdoul-Anziz and Seppecher 2018b). This is not a fundamentally new contribution, our main goal being to illustrate how the abstract method can be applied to specific problems. The first two examples demonstrate that the homogenization can naturally handle elastic structures possessing graded properties (Section 4.1) and pre-strain (Section 4.2), two features that are not commonly addressed in the literature. A variant of this truss, this time including *rigid* bars arranged in a way that the macroscopic strain is constrained, is proposed in the Appendix (Appendix E.9). We also demonstrate how the method can be extended to a frame made up of beams (rather than springs), and show that the homogenization method can be adapted to deliver a Timoshenko beam model (Section 4.3).

## 4.1 A truss lattice having slowly variable elastic properties

We consider a truss lattice comprising elastic bars connected by perfect hinges, arranged in rectangular geometry made up of square cells with side length $a$, see Figure 2.

<figure><img src="image_5.png" /><figcaption>**Figure 2** Truss lattice example. (a) General view, (b) a specific bar.</figcaption></figure>
**Discrete model.** Let $(e_1, e_2)$ be a frame of orthonormal directors with $e_1$ aligned with the longitudinal direction of the truss. The nodes are labelled with indices $\alpha_{\pm} = (\alpha, \pm)$ where $\alpha$ in an integer spanning the longitudinal direction, and the symbol $\pm$ is a discrete transverse coordinate labelling the lower versus upper layer of nodes. We denote by $X_{\alpha_{\pm}} = \bar{S}_{\alpha}e_1 \pm \frac{a}{2}e_2$ with $\bar{S}_{\alpha} = a\alpha$ the (undeformed) position of node $\alpha_{\pm}$, and by $\delta_{\alpha_{\pm}}$ its infinitesimal displacement.

Given a bar labelled by $\varphi$, we assign to it an orientation and denote by $\beta_{\varphi}^1$ and $\beta_{\varphi}^2$ its ordered end-nodes in reference configuration, by $l_{\varphi} = \|X_{\beta_{\varphi}^2} - X_{\beta_{\varphi}^1}\|$ its undeformed length and by $t_{\varphi} = (X_{\beta_{\varphi}^2} - X_{\beta_{\varphi}^1})/l_{\varphi}$ its undeformed unit tangent. The (dimensionless) axial strain in the bar is given by

$$ \varepsilon_{\varphi} = \frac{(\delta_{\beta_{\varphi}^2} - \delta_{\beta_{\varphi}^1}) \cdot t_{\varphi}}{l_{\varphi}}. \qquad (23) $$

The discrete elastic energy in the lattice is written in the form

$$ \Phi_d = \sum_{\varphi} w_{\varphi}(\varepsilon_{\varphi}) \quad \text{with} \quad w_{\varphi}(\varepsilon_{\varphi}) = \frac{1}{2} K (\bar{S}_{\varphi}^c) \varepsilon_{\varphi}^2. \qquad (24) $$

The lattice has graded properties as the elastic constant $K$ of the bars depends on the midpoint coordinate $\bar{S}_{\varphi}^c = \frac{1}{2}(\bar{S}_{\beta_{\varphi}^1} + \bar{S}_{\beta_{\varphi}^2})$. In order to keep the homogenization results as simple as possible, we make the simplifying assumption that the different types of bars have identical elastic constants $K$: the more natural assumption that all bars have identical cross-sections (and thus that their elastic constants $K$ is proportional to the length $l_{\varphi}$) could be addressed by making $K$ a function of not only $\bar{S}_{\varphi}^c$ but also of the type of the bar $\varphi$, which does not raise any particular difficulty.

**Continualization, scaling assumptions** Our continuous model is one-dimensional ($d=1$) and involves macroscopic fields that are functions of the longitudinal coordinate $\bar{S}$. As part of the continualization step, the nodal displacement $\delta_{\alpha_{\pm}}$ is sought in terms of continuous fields in the longitudinal direction as

$$ \delta_{\alpha_{\pm}} = \left( U(\bar{S}_{\alpha}) \mp \frac{a}{2} V'(\bar{S}_{\alpha}) + Y_{1}^{\pm}(\bar{S}_{\alpha}) \right) e_1 + \left( V(\bar{S}_{\alpha}) + Y_{2}^{\pm}(\bar{S}_{\alpha}) \right) e_2, \qquad (25) $$

where $U(\bar{S})$ and $V(\bar{S})$ denote the macroscopic longitudinal and transverse displacement of the equivalent rod, respectively, and $Y_i^{\pm}(\bar{S})$ are the components of the microscopic displacement for either row of nodes ($\pm$). The special case $Y_1^{\pm} = Y_2^{\pm} = 0$ corresponds to the (asymptotically incorrect) assumption of an unshearable model having rigid cross-sections—note that the term $\mp aV'/2$ represents the rigid rotation of the cross-section imposed by the centerline. We do not impose $Y_1^{\pm}$ and $Y_2^{\pm}$ to be zero.

We impose the kinematic constraint

$$ \langle Y_i^{\pm}(\bar{S}) \rangle = 0, \quad \forall \bar{S} \quad \forall i \in \{1, 2\}, \qquad (26) $$

where $\langle \cdot \rangle$ denotes the average over the top and bottom rows, i.e., $\langle Y_i^{\pm} \rangle = \frac{1}{2}(Y_i^{-} + Y_i^{+})$. This warrants that $(U(\bar{S}), V(\bar{S}))$ capture the average nodal displacement at coordinate $S$: indeed, by combining (25–26), we have $\langle \delta_{\alpha_{\pm}} \rangle = U(\bar{S}_{\alpha})e_1 + V(\bar{S}_{\alpha})e_2$. As a result, the equivalent rod passes through the midpoints of the sections $[X_{\alpha_{-}}X_{\alpha_{+}}]$. This choice is somewhat arbitrary: by using different weights in the average (26), we could introduce a lateral offset in the definition of the centerline.

We define the macroscopic length to be $a/\eta$, with $\eta \ll 1$ the scale separation parameter; see Equation (6). A standard scaling analysis yields the macroscopic stretching strain as $\varepsilon \sim U'(\bar{S})$, the macroscopic rotation as $\theta(\bar{S}) \sim V'(\bar{S})$, the relative rotation between successive transverse links as $a\theta'(\bar{S})$, and thus the differential stretching strain of bars located on the inner and outer sides (curvature effect) as $\varepsilon \sim a\theta'(\bar{S})$. Natural scales are found by balancing the two sources of stretching strain $\varepsilon$, yielding $\varepsilon \sim U' \sim a\theta' \sim aV''$, which suggests introducing dimensionless
quantities in the form

$$U(\bar{S}) = \frac{a}{\eta} u \left( \frac{\bar{S}}{a/\eta} \right) \qquad (27a)$$

$$V(\bar{S}) = \frac{a}{\eta^2} v \left( \frac{\bar{S}}{a/\eta} \right) \qquad (27b)$$

$$Y_i^{\pm}(\bar{S}) = a y_i^{\pm} \left( \frac{\bar{S}}{a/\eta} \right), \qquad (27c)$$

where $u$, $v$ and $y_i^{\pm}$ are dimensionless unknowns and $S = \bar{S}/(a/\eta) = \bar{S}\eta/a$ is a slow variable, i.e., the arclength scaled by the macroscopic length $a/\eta$.

In addition, we assume that the variations of elastic properties take place on the macroscopic scale, i.e.,

$$K(\bar{S}) = k \left( \frac{\bar{S}}{a/\eta} \right). \qquad (28)$$

In what follows, we eliminate the original quantities in favor of the dimensionless ones, $u(S)$, $v(S)$, $y_i^{\pm}(S)$ and $k(S)$ everywhere. By contrast with the rest of the paper, we keep explicit track of the small coefficient $\eta$ in the present section.

**Setting up the input of the homogenization procedure** In the classical theory of linear, planar beams, the two relevant strain measures are the stretching strain $e(S) := u'(S)$ and the bending strain $c(S) := v''(S)$. We thus anticipate that the homogenized energy will depend on the macroscopic strain $l$ ($n_l = 2$) which we define as

$$l(S) = (e(S), c(S)) := (u'(S), v''(S)). \qquad (29)$$

In terms of the unscaled displacement $(U, V)$, the dimensionless strain measures are given by $e(S) = U'(aS/\eta)$ and $c(S) = aV''(aS/\eta)$.

In addition, we define the vector of microscopic degrees of freedom ($n_y = 4$) as

$$\mathbf{y}(S) = (y_1^-(S), y_2^-(S), y_1^+(S), y_2^+(S)). \qquad (30)$$

These quantities $l$ and $y$ were chosen in such a way that the strain of a bar $\varphi$ can be expressed in terms of $l$, $y$ and their successive gradients at the midpoint coordinate $S_c$, see Equation (32) below.

In our discrete model (24), the elastic constants $k(S_c)$ depend on the midpoint coordinate $S_c$. This makes the energy density $W$ in (34) depend explicitly on $S$. Since $W$ is required by design to depend on $m$ and $E$ only, see (5), we pack up the coordinate $S$ into the list of material parameters $m(S)$,

$$m(S) = (S), \quad n_m = 1. \qquad (31)$$

An alternative (and ultimately equivalent) approach would be to define $m(S)$ as $k(S)$.

Next, we define the strain $E$ as the concatenation of (i) the discrete strains $E_i = \varepsilon_{\varphi}$ given in (23), in each of the 5 types $i$ of bars that make up the lattice, see Figure 2, together with (ii) the left-hand sides of the two kinematic constraints appearing in (26),

$$E = \left( E_1 \cdots E_5 \quad E_6 = \langle y_1^{\pm} \rangle \quad E_7 = \langle y_2^{\pm} \rangle \right). \qquad (32)$$

We therefore have $n_E = 7$ strain variables.

Having included the left-hand sides of the kinematic constraint (26) at positions 6 and 7 in $E$, we can easily express the $n_c = 2$ constraints in the form $Q \cdot E = 0$ expected in Equation (3), by defining the constraint extraction matrix as $Q = (\mathbf{0}_{2 \times 5} \ \mathbf{1}_{2 \times 2})$ using block-matrix notation.
Inserting the expression (25) of the nodal displacement in the expression of E, using the scaled quantities introduced in (27), and performing Taylor expansions about the midpoints of the bars, we get

$$
\begin{align*}
E_1 &= y_2^+ - y_2^- \\
E_2 &= l_1 + \frac{l_2}{2} + \eta y_1' + \frac{\eta^2}{48} (2l_1'' + l_2'') + O(\eta^3) \\
E_3 &= \frac{l_1 - y_1^- - y_2^- + y_1^+ + y_2^+}{2} - \frac{\eta}{24} (l_2' - 6(y_1'' + y_2'' + y_1^{++} + y_2^{++})) \\
    &\quad + \frac{\eta^2}{48} (l_1'' - 3(y_1'' + y_2'' - y_1^{++} - y_2^{++})) + O(\eta^3) \\
E_4 &= l_1 - \frac{l_2}{2} + \eta y_1' + \frac{\eta^2}{48} (2l_1'' - l_2'') + O(\eta^3) \\
E_5 &= \frac{l_1 + y_1^- - y_2^- - y_1^+ + y_2^+}{2} + \frac{\eta}{24} (l_2' + 6(y_1'' - y_2'' + y_1^{++} - y_2^{++})) \\
    &\quad + \frac{\eta^2}{48} (l_1'' + 3(y_1'' - y_2'' - y_1^{++} + y_2^{++})) + O(\eta^3) \\
E_6 &= (y_1^- + y_1^+)/2 \\
E_7 &= (y_2^- + y_2^+)/2,
\end{align*}
\tag{33}
$$

which is of the form $E = E(m; l, l', l'', \dots; y, y', y'', \dots)$ expected in Equation (2). In (33), the argument of the functions $l, l', l'', y, y', y''$ is implicitly assumed to be the midpoint $S_{\varphi}^c$ of each bar. The details of the calculations can be found in the companion Mathematica notebook.

The discrete energy (24) is finally continualized in the canonical form (4-5) as

$$ \Phi_d \approx \int_{-\infty}^{+\infty} W(m, E) dS \quad \text{where } W(m, E) = \frac{k(m_1)}{2\eta} \sum_{i=1}^{5} E_i^2. \tag{34} $$

The notation $k(m_1)$ conveys the fact that the argument $m_1 = S$ to be passed to the stiffness distribution $k(S)$ is the first (and only) component $m_1$ of $m$, see (31).

In (34), the discrete energy (24) has been continualized by using the formal rule

$$ \sum_{\varphi} w_{\varphi} = \sum_{i=1}^{5} \sum_{\varphi \in i} w_{\varphi} \approx \sum_{i=1}^{5} \int_{-\infty}^{+\infty} w_i \frac{dS}{\eta} = \sum_{i=1}^{5} \int_{-\infty}^{+\infty} \frac{1}{2} k E_i^2 \frac{dS}{\eta} = \int_{-\infty}^{+\infty} W dS, \tag{35} $$

where $i$ is an index running over the 5 different types of links, $\varphi \in i$ is included in the partial sum of all links $\varphi$ belonging to a particular family $i$, and $w_i$ denotes the expression of the energy $w_{\varphi}$ relevant to the links $\varphi$ belonging to a particular family $i$. The coefficient $1/\eta$ appearing in the definition of $W$ in (34) is nothing but the lineic density of links of each type per unit dimensionless length $S$. Thanks to the assumed periodicity of the lattice, we have been able to rewrite the discrete sum in (24) into an integral in (34).

**Remark 8** In the argument sketched above, the continualization is based on the formal approximation rule for a discrete sum in terms of an integral, $\sum_{\varphi \in i} w_{\varphi} \approx \int w_i \frac{dS}{\eta}$, which has been used in a number of earlier work. A rigorous justification of this approximation is known as the Euler-MacLaurin formula: it is correct to order $\eta^2$ in an infinite domain, but needs to be corrected by boundary terms in a finite domain.

At this point, we can identify the quantities that are required on input of the homogenization method, as specified in Table 1: they are listed in Table 4. The integer constants appearing in the left column of Table 4 and the constraint extraction matrix $Q$ have been collected from the above discussion. The tensors $E_l, \dots, E_l''$ collect the numerical coefficients appearing in (33) and are found by identification with Equation (2). The elastic stiffness tensor $K(S)$ is found by identifying (34) with (5). Expressions in Table 4 make use of the notation $\delta_i^n = (0_1 \dots 0_{i-1} 1_i 0_{i+1} \dots 0_n) \in \mathbb{R}^n$ for the $i$-th Kronecker vector with length $n$; for example, $\delta_2^5 = (0 1 0 0 0)$.
<table><caption>**Table 4** Input parameters used for the homogenization of the truss lattice. The spring constant $k(S)$ is set up as a symbolic function of $S$. For the sake of brevity, the terms denoted by ellipses are omitted: full expressions are available in the input file inhomogeneous-truss.nb included in the library (Audoly 2023).</caption><tr><td></td><td></td></tr><tr><td>d = 1</td><td>E<sub>l</sub>(S) = (δ<sub>2</sub><sup>7</sup> + δ<sub>4</sub><sup>7</sup> + δ<sub>3</sub><sup>7</sup> + δ<sub>5</sub><sup>7</sup>) ⊗ δ<sub>1</sub><sup>2</sup> + δ<sub>2</sub><sup>7</sup> - δ<sub>4</sub><sup>7</sup> ⊗ δ<sub>2</sub><sup>2</sup></td></tr><tr><td>n<sub>m</sub> = 1</td><td>E'<sub>l</sub>(S) = <sup>1</sup>&frasl;<sub>24</sub>(-δ<sub>3</sub><sup>7</sup> + δ<sub>5</sub><sup>7</sup>) ⊗ δ<sub>2</sub><sup>2</sup> ⊗ δ<sub>1</sub><sup>1</sup></td></tr><tr><td>m = (S)</td><td>E''<sub>l</sub>(S) = <sup>1</sup>&frasl;<sub>24</sub>E<sub>l</sub>(S) ⊗ δ<sub>1</sub><sup>1</sup> ⊗ δ<sub>1</sub><sup>1</sup></td></tr><tr><td>n<sub>l</sub> = 2</td><td>E<sub>y</sub>(S) = δ<sub>1</sub><sup>7</sup> ⊗ ... + δ<sub>3</sub><sup>7</sup> ⊗ ... + δ<sub>5</sub><sup>7</sup> ⊗ ... + δ<sub>6</sub><sup>7</sup> ⊗ <sup>δ<sub>4</sub><sup>4</sup> + δ<sub>3</sub><sup>4</sup></sup>&frasl;<sub>2</sub> + δ<sub>7</sub><sup>7</sup> ⊗ <sup>δ<sub>2</sub><sup>4</sup> + δ<sub>4</sub><sup>4</sup></sup>&frasl;<sub>2</sub></td></tr><tr><td>n<sub>y</sub> = 4</td><td>E'<sub>y</sub>(S) = ...</td></tr><tr><td>n<sub>E</sub> = 7</td><td>E''<sub>y</sub>(S) = <sup>1</sup>&frasl;<sub>16</sub>(δ<sub>3</sub><sup>7</sup> ⊗ (-δ<sub>1</sub><sup>4</sup> - δ<sub>2</sub><sup>4</sup> + δ<sub>3</sub><sup>4</sup> + δ<sub>4</sub><sup>4</sup>) + δ<sub>5</sub><sup>7</sup> ⊗ (δ<sub>1</sub><sup>4</sup> - δ<sub>2</sub><sup>4</sup> - δ<sub>3</sub><sup>4</sup> + δ<sub>4</sub><sup>4</sup>)) ⊗ δ<sub>1</sub><sup>1</sup> ⊗ δ<sub>1</sub><sup>1</sup></td></tr><tr><td>n<sub>c</sub> = 2</td><td>K(S) = <sup>k(S)</sup>&frasl;<sub>η</sub> Σ<sub>i=1</sub><sup>5</sup> δ<sub>i</sub><sup>7</sup> ⊗ δ<sub>i</sub><sup>7</sup></td></tr><tr><td></td><td>Q(S) = Σ<sub>i=1</sub><sup>2</sup> δ<sub>i</sub><sup>2</sup> ⊗ δ<sub>5+i</sub><sup>7</sup></td></tr></table>

**Homogenization results** We propose an implementation of the general homogenization method described in this paper (Sections 5 and 6) in the form of an open-source library named shoal (Audoly 2023). The quantities listed in Table 4 are passed to the library in symbolic form, see the input file inhomogeneous-truss.nb<fnref n="1" />. The code executes without any further input from the user. It delivers the tensors listed in Table 3, which we now proceed to interpret using (13–16) and (29–30). We note that $\nabla m = \delta_1^1 \otimes \delta_1^1$ in view of the definition (31) of $m$. The homogenization method returns

$$ 
\begin{aligned}
Y_0 &= \frac{1}{6}(\delta_2^4 - \delta_4^4) \otimes \delta_1^2, & Y'_0 &= \frac{11}{24}(\delta_1^4 - \delta_3^4) \otimes \delta_2^2 \otimes \delta_1^1, \\
Y_1 &= \frac{k'(S)}{2k(S)}(\delta_1^4 - \delta_3^4) \otimes \delta_2^2 \otimes \delta_1^1 \otimes \delta_1^1.
\end{aligned} 
\qquad (36) $$

In view of (15–16), the microscopic displacement $y = Y_0 \cdot l + Y'_0 : \nabla l + (Y_1 : \nabla m) \cdot l + \dots$ is found with the help of (29–30) as

$$ 
\begin{aligned}
y_1^- &= + \frac{11}{24} \eta c'(S) + \frac{\eta k'(S)}{2k(S)} c(S) + O(\eta^2) \\
y_2^- &= + \frac{e(S)}{6} + O(\eta^2) \\
y_1^+ &= - \frac{11}{24} \eta c'(S) - \frac{\eta k'(S)}{2k(S)} c(S) + O(\eta^2) \\
y_2^+ &= - \frac{e(S)}{6} + O(\eta^2)
\end{aligned} 
\qquad (37) $$

The quantities $\pm e(S)/6$ of order $\eta^0$ in the right-hand side represent the predictions of classical (leading-order) homogenization. Note that the correction of order $\eta^1$ includes not only a contribution proportional to the strain gradient $c'$ but also one proportional to the gradient $k'$ of elastic properties—this corresponds, respectively, to the $\nabla l$ and $\nabla m$ contributions appearing in the right-hand side of (22b).

The homogenized energy functional is obtained by interpreting the output of the code similarly using Table 3, see also (13–14),

$$ \Phi^* = \int_{-\infty}^{+\infty} \left[ \frac{1}{2} \left( \frac{7k}{3} e^2 + \left( \frac{k}{2} - \frac{\eta^2 k'^2}{2k} \right) c^2 \right) - \eta^2 k' \left( \frac{ee'}{12} + \frac{23}{48} cc' \right) - \frac{\eta^2 k}{2} \left( \frac{1}{6} e'^2 + \frac{11}{24} c'^2 \right) + O(\eta^3) \right] \frac{dS}{\eta} \qquad (38) $$

The detailed expressions of the tensors $K_0, K_1, \dots$ underlying the above expression of $\Phi^*$ can be found in the Mathematica notebook inhomogeneous-truss.nb included in the library. The energy $\Phi^*$ in (38) depends on the stretching strain $e(S)$, on the dimensionless curvature strain

<footnote n="1">¹ See shoal-library-v1.0/discrete_engine/illustrations/JTCAM2023-paper/inhomogeneous-truss.nb</footnote>
$c(S)$, as well as on their gradients $e'(S)$ and $c'(S)$, and on the gradients of elastic properties $k'(S)$ (the dependence of $e$, $c$ and $k$ on $S$ implicit in the above integral). This result is valid for any given slowly varying distribution of spring stiffness $k(S)$ since $k(S)$ is treated as a symbolic function. As we consider an infinite structure, we have ignored in (38) the boundary terms that are returned by the homogenization library.

In view of the negative coefficients $-1/6$ and $-11/24$ appearing in the last term in the integral, the energy $\Phi^*$ can be made arbitrarily large and negative by incorporating oscillations into the unknowns $e(S)$ and $c(S)$, having both small amplitude and small wavelength. This holds even in the simple case where the elastic properties are uniform ($k' = 0$). This points to the lack of lower semi-continuity of the homogenized energy $\Phi^*$, an undesirable property of higher-order gradient models that has been documented by several authors for various elastic structures, see for instance (Le and Marigo 2018). It calls for a regularization of the functional $\Phi^*$, a point which we address in future work.

**Remark 9** Alternatively, the homogenized functional $\Phi^*$ can be expressed in terms of the original unscaled quantities, namely the elastic constant $K(\bar{S})$, the axial strain $\bar{e}(\bar{S}) = U'(\bar{S}) = e(\bar{S}\eta/a)$ and the curvature $\bar{c}(\bar{S}) = V''(\bar{S}) = \frac{1}{a}c(\bar{S}\eta/a)$. The result is free of the parameter $\eta$,

$$ \Phi^* = \int_{-\infty}^{+\infty} \left[ \frac{1}{2} \left( \frac{7K}{3}\bar{e}^2 + \left( \frac{K}{2} - \frac{a^2 K'^2}{2k} \right) c^2 \right) - a^2 K' \left( \frac{\bar{e}\bar{e}'}{12} + \frac{23}{48}\bar{c}\bar{c}' \right) - \frac{a^2 K}{2} \left( \frac{1}{6}\bar{e}'^2 + \frac{11}{24}\bar{c}'^2 \right) + O(a^3) \right] \frac{d\bar{S}}{a}. \quad (39) $$

## 4.2 A truss lattice with pre-strain

In order to illustrate the ability of our method to handle pre-strain, we now include an additional tensile pre-strain +$p(S)/2$ in the lower layer of bars, and a contractile pre-strain -$p(S)/2$ in the upper layer of bars. For the sake of simplicity, we focus on the case of uniform elastic properties across the length, taking $k$ to be independent of $S$.

Accordingly, we change the strain definition for the two diagonal bars to $E_2 = E_2^0 - p(S)/2$ and $E_4 = E_4^0 - (-p(S)/2)$, where the quantities $E_i^0$ refer to the expressions of $E^i$ in the absence of pre-strain given in the right-hand sides of (33) (Section 4.1).

Since $E$ is required to be a function of $l$ and $y$ and their derivatives, we add a third component $l_3 = p$ to the macroscopic strain $l$ and rewrite $E_2 = E_2^0 - l_3(S)/2$ and $E_4 = E_4^0 + l_3(S)/2$. Note that the pre-strain is subtracted from the original element strain, so that the energy in the bar labelled 2 in Figure 2, for instance, is given by $\frac{k}{2}(E_2^0 - p(S)/2)^2$ and is minimum when $E_2^0 = p(S)/2$.

To sum up, we make the following changes in the specification of the input problem:

$$ \begin{gathered} m = (\cdot), \quad n_m = 0, \quad n_l = 3, \quad l = (e = u', c = v'', p), \\ E_2 = E_2^0 - \frac{l_3(S)}{2}, \quad E_4 = E_4^0 + \frac{l_3(S)}{2}, \quad \mathcal{K} = \frac{k}{\eta} \sum_{i=1}^{5} \delta_i^7 \otimes \delta_i^7. \end{gathered} \quad (40) $$

The tensors $E_y$, $E'_y$ and $E''_y$ are unchanged. The dimension of the tensors $E_l$, $E'_l$ and $E''_l$ is increased to reflect the new dimension $n_l = 3$ of the macroscopic strain vector $l$. This change of dimension takes place simply by adding zero entries to $E'_l$ and $E''_l$, and by including a new contribution $E_l = \cdots + \frac{1}{2}(-\delta_2^7 + \delta_4^7) \otimes \delta_3^3$ to $E_l$, to reflect the new $l_3$-terms in $E_2$ and $E_4$.

The homogenization code is run again with the modified set of input parameters, see the notebook pretrained-truss nb included in the library. The homogenized energy is obtained as

$$ \Phi^* = \int_{-\infty}^{+\infty} \left[ \frac{k}{2} \left( \frac{7}{3} e^2 + \frac{(c-p)^2}{2} \right) + \frac{\eta^2 k}{2} \left( -\frac{e'^2}{6} - \frac{11}{24} c'^2 + \frac{23}{24} c'p' - \frac{p'^2}{2} \right) + O(\eta^3) \right] \frac{dS}{\eta}. \quad (41) $$

In the leading-order term, the pre-strain sets the natural curvature of the equivalent beam as $c_0(S) = p(S)$, as could be anticipated. Note that the gradient of pre-strain $p'(S)$ contributes to the energy at order $\eta^2$.
<figure><img src="image_6.png" /><figcaption>**Figure 3** An elastic frame made up of identical beams having length a, stretching modulus EA and bending modulus EI. (a) General view, (b) the three modes of deformation of a linear beam, illustrated here for a beam β in the family φ = 1: stretching ε<sub>φ</sub>, bending κ<sub>φ</sub> and shearing τ<sub>φ</sub>.</figcaption></figure>

4.3 Homogenizing an elastic frame into a Timoshenko model

We now analyze the frame made up of linear elastic beams, shown in Figure 3. Specifically, we show that the homogenization procedure can deliver a Timoshenko beam model, when the rotation of the cross-sections is added to the list of macroscopic parameter *l*, implying that this parameter is kept fixed during the energy relaxation. This example has been treated in (Abdoul-Anziz and Seppecher 2018b) and we use their results to verify our method.

**Changes to the homogenization procedure** To address the elastic frame, we start over from the model in Section 4.1, with the following changes.

• There are 3 families of beams (and not 5) as indicated by the labels φ ∈ {1, 2, 3} in Figure 3(a).

• We revert to the case of uniform properties by discarding the parameter *m*.

• We change the scaling assumption on the transverse displacement *V* in (27b) to

$$
V(\bar{S}) = \frac{a}{\eta} v \left( \frac{\bar{S}}{a/\eta} \right), \qquad (42)
$$

in such a way that the rotation $V'$ (and not the curvature $V''$) is of order $\eta^0$: this delivers the
Timoshenko model is a more convenient form.

* To model the stiff junctions, we introduce the nodal rotation $\theta_{\alpha\pm}$ at a node $\alpha_{\pm}$ (in addition to the nodal displacement $\boldsymbol{\delta}_{\alpha\pm}$ from (25)): $\theta_{\alpha\pm}$ is given in terms of two additional unknown continuous functions $Y_{\theta}^{\pm}(\bar{S})$ and $Y_{\theta}^{\pm}(\bar{S})$ as

$$
\theta_{\alpha_{\pm}} = V'(\bar{S}_{\alpha}) + Y_{\theta}^{\pm}(\bar{S}_{\alpha}). \quad (43)
$$

Including the centerline rotation $V'$ as the first term in the right-hand side makes the microscopic displacements $Y_{\theta}^{\pm}$ invariant by rigid-body rotations, so that they can be expressed in terms of the strain as assumed in our procedure.

* For a beam of type φ, the classical theory of linear (Euler-Bernoulli) beams can be summarized as follows. In addition to the stretching strain ε<sub>φ</sub> in (23), one has to consider a curvature strain κ<sub>φ</sub> and a shearing strain τ<sub>φ</sub>, as sketched in Figure 3(b),

$$
\kappa_{\varphi} = \theta_{\beta_{\varphi}^2} - \theta_{\beta_{\varphi}^1}, \quad \tau_{\varphi} = \frac{\theta_{\beta_{\varphi}^1} + \theta_{\beta_{\varphi}^2}}{2} - \frac{(\delta_{\beta_{\varphi}^2} - \delta_{\beta_{\varphi}^1}) \cdot (\mathbf{e}_3 \times \mathbf{t}_{\varphi})}{l_{\varphi}}. \quad (44)
$$

Here, $\tau_{\varphi}$ is an *apparent* shearing strain at the scale of unit cells, which actually resolves into *bending* microscopically, see Figure 3(b): our Euler-Bernoulli ‘microscopic’ beam model is shearless. The elastic energy of the beam is then of the form

$$
w_{\varphi} = \frac{EAa}{2} \varepsilon_{\varphi}^2 + \frac{EI}{2a} (\kappa_{\varphi}^2 + 12\tau_{\varphi}^2). \qquad (45)
$$

* We define the (apparent) rotation of the cross-sections $\bar{\gamma}$ as the rotation of the line passing through the nodes $\alpha_-$ and $\alpha_+$ facing each other,

$$
\bar{\gamma}(\bar{S}) = -\frac{1}{a}(\delta_{\alpha_+} - \delta_{\alpha_-}) \cdot \mathbf{e}_1 = V'(\bar{S}) - \frac{Y_1^+(\bar{S}) - Y_1^-(\bar{S})}{a}. \quad (46)
$$
* Next, we consider the following scaled macroscopic parameters: centerline extensional strain $e(S) = u'(S) = U'(aS/\eta)$, centerline curvature $c(S) = v''(S) = \frac{a}{\eta}V''(aS/\eta)$, rotation of cross-sections $\gamma(S) = v'(S) - (y_{\theta}^{+}(S) - y_{\theta}^{-}(S)) = \bar{\gamma}(aS/\eta)$ and shear angle

$$g(S) = \gamma(S) - v'(S) = -(y_{\theta}^{+}(S) - y_{\theta}^{-}(S)). \qquad (47)$$

* The vector of microscopic degrees of freedom $\mathbf{y}$ is extended to reflect the presence of the two additional fields capturing nodal rotation: Equation (30) is changed to

$$\mathbf{y}(S) = (y_{1}^{-}(S), y_{2}^{-}(S), y_{\theta}^{-}(S), y_{1}^{+}(S), y_{2}^{+}(S), y_{\theta}^{+}(S)). \qquad (48)$$

* A key modification to the homogenization procedure is that we include the shear angle in the list of macroscopic parameters,

$$l(S) = (e(S), c(S), g(S)), \qquad (49)$$

by contrast with (29). Equation (47) can then be rewritten as

$$l_3(S) + y_{\theta}^{+}(S) - y_{\theta}^{-}(S) = 0 \quad \forall S. \qquad (50)$$

This kinematic constraint is taken care of by including the left-hand side of the above equation as a component of strain (specifically, $E_{12}$) and by extending the constraint matrix $Q$ from Table 4 to include a third row filled with zeros, except for an entry equal 1 in column 12, see Equation (3).

* The strain $E \in \mathbb{R}^{12}$ is now a vector of length 12, made up of the three strain measures ($\varepsilon_{\varphi}, \kappa_{\varphi}, \tau_{\varphi}$) for each of the three types of beams, $\varphi \in \{1, 2, 3\}$, the two left-hand sides in the zero-average-displacement constraints (26), and the left-hand side of (50).

* The discrete energy is finally approximated by a continuous integral $\Phi = \int_{-\infty}^{+\infty} W dS$ capturing the elastic potentials of the beams (45)

$$W(E) = \frac{EAa}{2\eta} (E_1^2 + E_4^2 + E_7^2) + \frac{EI}{2a\eta} ((E_2^2 + E_5^2 + E_8^2) + 12(E_3^2 + E_6^2 + E_9^2)). \qquad (51)$$

The coefficient $\eta$ in the denominators is produced when the discrete sum over beams is approximated by an integral, as earlier in (34).

The scale separation parameter $\eta \ll 1$ is unspecified so far, and we set it to be the small aspect-ratio of the beams, i.e., we set

$$\eta := \sqrt{\frac{EI}{a^2 EA}}. \qquad (52)$$

Doing so, we are anticipating that, in the forthcoming Timoshenko model, the shearing mode will relax over a typical macroscopic length $a/\eta$ (this macroscopic length $a/\eta$ has been introduced earlier in the definition of the scaled arclength $\bar{S} = S\eta/a$). The validity of this assumption will be checked later, by verifying that the various terms in the Timoshenko model are produced with consistent powers of $\eta$, see Equation (55) below.

Equation (52) is used to eliminate $EI$ in favor of $EA$ and $\eta$ in the elastic potential (51): in the homogenization procedure, the energy $\Phi$ provided on input may depend explicitly on $\eta$.

**Results of the homogenization procedure** The homogenization procedure is carried out in the Mathematica notebook frame-timoshenko nb included in the library. With the boundary terms omitted (case of an infinitely long truss), the result of the homogenization procedure are interpreted as

$$\Phi^{\star} = \int_{-\infty}^{+\infty} (W_{\text{ext}}^{\star}(e(S), e'(S)) + W_{\text{Tm}}^{\star}(g(S), c(S), g'(S)) + O(\eta^3)) \frac{dS}{\eta}, \qquad (53)$$

where $W_{\text{ext}}^{\star}(e, e')$ is a higher-order gradient bar model applicable to the longitudinal displacement $u$,

$$W_{\text{ext}}^{\star}(e, e') = \frac{aEA}{2} \left( 2e^2 - \frac{\eta^2}{6} e'^2 \right), \qquad (54)$$
and $W_{Tm}^*(g, c, g')$ is a Timoshenko beam model applicable to the transverse displacement $v$,

$$W_{Tm}^*(g, c, g') = \frac{\eta^2 a EA}{2} \left( \frac{1}{2} (c + g')^2 + 8g^2 \right). \qquad (55)$$

The details can be found in the Mathematica notebook.

**Remark 10** The square term $\eta^2(c + g')^2$ appearing in $W_{Tm}^*$ is first obtained in expanded form, with the different terms $\eta^2 c^2$, $2\eta^2 cg'$ and $\eta^2 g'^2$ appearing respectively in the successive contributions $\Phi_{[0]}^*, \Phi_{[1]}^*$ and $\Phi_{[2]}^*$ to the energy. The factored form in the equation above emerges when $\Phi^* = \Phi_{[0]}^* + \Phi_{[1]}^* + \Phi_{[2]}^*$ is simplified.

With $W_{ext}^*(e, e') = O(\eta^0)$ while $W_{Tm}^*(g, c, g') = O(\eta^2)$, the stretching energy is dominant: this agrees with the fact that Abdoul-Anziz and Seppecher (Abdoul-Anziz and Seppecher 2018b) reported an inextensible model ($e \equiv 0$).

Turning therefore attention to the transverse displacement, we eliminate the shear angle $g$ in favor of the apparent rotation $\gamma = v' + g$ of cross-sections. Recalling the definition of the scaled curvature $c = v''$, and rearranging using (52), we obtain the Timoshenko model in standard form,

$$W_{Tm}^* = \frac{1}{2} \left( \lambda_T (\eta \gamma'(S))^2 + \zeta_T (v'(S) - \gamma(S))^2 \right), \qquad (56)$$

and the elastic moduli are identified as

$$\lambda_T = \frac{aEA}{2}, \qquad \zeta_T = \frac{8EI}{a}. \qquad (57)$$

The same elastic frame has been homogenized in (Abdoul-Anziz and Seppecher 2018b) and the authors reported homogenized moduli $\lambda_T = 1/2$ and $\zeta_T = 2$. These values are a special case of (57), implying that our homogenization results are consistent. Indeed, a careful analysis of the beam model used by these authors reveals that they limited attention to the special case $EA = 1/a$ and $EI = a/4$ (when setting both their elastic moduli $a_{AS}$ and $f_{AS}$ to 1, see the note frame-timoshenko-AS18-beam-model1.pdf in the library).

**Remark 11** In terms of the original (un-scaled) deflection $V(\bar{S})$ and cross-section rotation $\bar{\gamma}(\bar{S}) = \gamma(\bar{S}\eta/a)$, the Timoshenko model appearing in (53-56) can be rewritten as

$$\Phi_{\text{transv}}^{*} = \int_{-\infty}^{+\infty} \frac{1}{2} \left( \frac{EAa^2}{2} \bar{\gamma}'^2(\bar{S}) + \frac{8EI}{a^2} (V'(\bar{S}) - \bar{\gamma}(\bar{S}))^2 \right) d\bar{S}. \qquad (58)$$

**Remark 12** We obtained a Timoshenko model as a consequence of the fact that the shear angle $g = l_3$ is treated as a macroscopic parameter, and not just because the lattice is made up of beams. Indeed, lattices made up of bars ($EI = 0$) can yield Timoshenko models as well when homogenized. The motivation for using beams in this example was to compare with the earlier work of (Abdoul-Anziz and Seppecher 2018b).

# 5 Derivation of the homogenized model

## 5.1 Leading order (classical homogenization)

At order $\eta^0$, the microscopic displacement (21) is given by $y(X) = y_{[0]}(X) + O(\eta)$. The gradients terms $\nabla^k m$, $\nabla^k l$ and $\nabla^k y$ are of order $\eta^k$ by (7) and can be ignored for $k \ge 1$. With the gradients neglected, we denote the microscopic strain in (2) as

$$E^{(0)}(m, l, y) = E(m; l, 0, 0, \dots; y, 0, 0, \dots) = E_l(m) \cdot l + E_y(m) \cdot y, \qquad (59)$$

and the bulk energy density in (5) as

$$W^{(0)}(m, l, y) = W(m, E^{(0)}(m, l, y)). \qquad (60)$$
The order $\eta^0$ approximation of the strain energy (4) can then be obtained as $\Phi[m, l, y] = \Phi_{[0]}[m, l, y_{[0]}] + O(L^d\eta)$, where

$$ \Phi_{[0]}[m, l, y_{[0]}] = \int_{\Omega} W^{(0)}(m(X), l(X), y_{[0]}(X))dX. \qquad (61) $$

At leading order $\eta^0$, the variational problem (9) can be written as

$$ \left\{ \begin{array}{l} Q \cdot (E_y(m(X)) \cdot y_{[0]}(X) + E_l(m(X)) \cdot l(X)) = 0 \quad \forall X \quad (62a) \\ D_y \Phi_{[0]}[m, l, y_{[0]}; \delta y] + \int_{\Omega} ((Q \cdot E_y(m(X)))^T \cdot g_{[0]}(X)) \cdot \delta y(X) dX = 0 \quad \forall \delta y. \quad (62b) \end{array} \right. $$

$$ \left\{ \begin{array}{l} Q \cdot (E_y(m(X)) \cdot y_{[0]}(X) + E_l(m(X)) \cdot l(X)) = 0 \quad \forall X \\ D_y \Phi_{[0]}[m, l, y_{[0]}; \delta y] + \int_{\Omega} ((Q \cdot E_y(m(X)))^T \cdot g_{[0]}(X)) \cdot \delta y(X) dX = 0 \quad \forall \delta y. \end{array} \right. \quad (62b) $$

Its solution is denoted as $(y_{[0]}, g_{[0]}) = (y_{[0]}^*, g_{[0]}^*)$.

No gradient of $\nabla y_{[0]}$ is present in the expression of $\Phi_{[0]}$ in (61) nor in the integral in (62b), implying that this variational problem is *local*: at any point $X$, we must solve the following problem for the unknowns $y_{[0]}^*(X)$ and $g_{[0]}^*(X)$,

$$ Q \cdot E_y(m(X)) \cdot y_{[0]}^*(X) + Q \cdot E_l(m(X)) \cdot l(X) = 0 \qquad (63) $$

$$ \frac{\partial W^{(0)}}{\partial y}(m(X), l(X), y_{[0]}^*(X)) + (Q \cdot E_y(m(X)))^T \cdot g_{[0]}^*(X) = 0, $$

where we have used the expression of $\Phi_{[0]}$ in (61). The solution $(y_{[0]}^*(X), g_{[0]}^*(X))$ at any particular point $X$ depends on the local values of $m(X)$ and $l(X)$ only. In Appendix B, the solution $y_{[0]}^*(X)$ is obtained in the form announced earlier in (22a),

$$ y_{[0]}^*(X) = Y_0(m(X)) \cdot l(X), \qquad (64) $$

and an explicit expression for the localization tensor $Y_0(m)$ is given in (B.7-B.8).

Inserting (64) into (61), we derive in Appendix B the dominant contribution to the energy $\Phi_{[0]}^*[m, l] = \Phi_{[0]}[m, l, y_{[0]}^*]$ that was announced in (17), namely

$$ \Phi_{[0]}^*[m, l] = \int_{\Omega} \frac{1}{2} l(X) \cdot K_0(m(X)) \cdot l(X) dX. \qquad (65) $$

The expression of the leading-order stiffness tensor $K_0(m)$ is given in Equation (B.12) in the Appendix.

## 5.2 Analysis of the gradient effect

We now proceed to solve the next orders in the microscopic displacement, see (21) and (64),

$$ y(X) = Y_0(m(X)) \cdot l(X) + y_{[1]}(X) + y_{[2]}(X) + \dots \qquad (66) $$

Inserting this into (4), we derive a Taylor expansion of the energy as

$$ \Phi = \Phi_{[0]}^*[m, l] + \Phi_{[1]}^*[m, l] + \Phi_{[2]}[m, l, y_{[1]}] \qquad (67) $$

where $\Phi_{[k]}$ denotes the term of order $\eta^k$, and the star is used to mark energy contributions that do not depend on the yet-unknown corrector $y_{[1]}$. The dominant term $\Phi_{[0]}^*$ is the functional found earlier in (65), while the next-order terms $\Phi_{[1]}$ and $\Phi_{[2]}$ are obtained in Equations (C.20) and (C.24) in Appendix C as

$$ \Phi_{[1]}^*[m, l] = \int_{\Omega} (\mathcal{A}(m) : h \otimes h) : \nabla h(X) dX \qquad (68a) $$

$$ \Phi_{[2]}[m, l, y_{[1]}] = \int_{\Omega} \left( \begin{aligned} & (\hat{\mathcal{B}}^{(0)}(m) : (h \otimes h)) \therefore \frac{\nabla h \otimes \nabla h}{2} \\ & + (\hat{\mathcal{B}}^{(1)}(m) \cdot h) \therefore (\nabla h \otimes y_{[1]}) + W_{yy}(m) : \frac{y_{[1]} \otimes y_{[1]}}{2} \end{aligned} \right) dX \qquad (68b) \\ + \int_{\Omega} \left( (C^{(0)}(m) : (h \otimes h)) \therefore \nabla^2 h + (C^{(1)}(m) \cdot h) : \nabla y_{[1]} \right) dX, $$
Closed-form expressions for the operators appearing in the right-hand side are derived in Appendix C. In the right-hand sides above, we have introduced the vector **h** = (l, m, (1)) obtained by concatenating the macroscopic variables **l** and **m**, and adding a trailing entry 1, see Appendix C.1: this notation trick simplifies the calculations significantly.

It is remarkable that $\Phi^*_{[1]}[m, l]$ does not depend on the corrector $y_{[1]}$, even though both are of order $\eta$. As a result, the first correction $\Phi^*_{[1]}$ depends on the macroscopic variables $(m, l)$ only, as conveyed by the star notation which we reserve for the output of the homogenization procedure.

For a similar reason explained in the Appendix, $\Phi_{[2]}[m, l, y_{[1]}]$ does not depend on $y_{[2]}$ even though both are of order $\eta^2$. It does depend on the unknown correction $y_{[1]}$, however. The gradient term $\nabla y_{[1]}$ appearing in the integrand of $\Phi_{[2]}$ in (68b) can be removed by integrating by parts the $C^{(1)}$ term—the benefit is that the problem of optimizing $\Phi_{[2]}$ with respect to the function $y_{[1]}(X)$ then becomes a local problem. We choose to integrate the $C^{(0)}$ by parts as well, as the result can be merged with the $\hat{B}^{(0)}(m)$ term. These two integration by parts are carried out in the Appendix and the result is

$$ \Phi_{[2]}[m, l, y_{[1]}] = \Phi_{[2]}^{\text{bt}}[m, l, y_{[1]}] + \Phi_{[2]}^{\text{it}}[m, l, y_{[1]}], \quad (69) $$

where the boundary terms $\Phi_{[2]}^{\text{bt}}$ and integral terms $\Phi_{[2]}^{\text{it}}$ are given in (C.29), respectively, as

$$ \Phi_{[2]}^{\text{bt}}[m, l, y_{[1]}] = \oint_{\partial\Omega} \left( (C^{(0)}(m) : (h \otimes h)) \circ (\nabla h \otimes n) + (C^{(1)}(m) \cdot h) : (y_{[1]} \otimes n) \right) da \quad (70) $$

and

$$ \Phi_{[2]}^{\text{it}}[h, y_{[1]}] = \int_{\Omega} \left( \begin{aligned} & (\mathcal{B}^{(0)}(m) : (h \otimes h)) \circ \frac{\nabla h \otimes \nabla h}{2} \\ & + (\mathcal{B}^{(1)}(m) \cdot h) \circ (\nabla h \otimes y_{[1]}) + \mathcal{W}_{yy}(m) : \frac{y_{[1]} \otimes y_{[1]}}{2} \end{aligned} \right) dX. \quad (71) $$

As anticipated, the gradient terms $\nabla y_{[1]}$ have all disappeared.

Having worked out the expansion of the energy, we proceed to solve the variational problem (9) at order $\eta$: inserting the expansion $y(X) = y^*_{[0]}(X) + y_{[1]}(X) + \cdots$ into the energy (67) we get a variational problem for the corrector $y_{[1]}$ and a Lagrange multiplier $g_{[1]}$,

$$ \left\{ \begin{array}{l} Q \cdot (E_y(m(X)) \cdot y_{[1]}(X) + (\mathcal{J}^1(m(X)) \cdot h(X)) : \nabla h(X)) = 0 \quad \forall X \quad (72a) \\ D_y \Phi_{[2]}[m, l, y_{[1]}; \delta y] + \int_{\Omega} g_{[1]}(X) \cdot Q \cdot E_y(m(X)) \cdot \delta y(X) dX = 0 \quad (\forall \delta y). \end{array} \right. $$

$$ \int_{\Omega} g_{[1]}(X) \cdot Q \cdot E_y(m(X)) \cdot \delta y(X) dX = 0 \quad (\forall \delta y). \quad (72b) $$

The incremental form of the kinematic constraint appearing in (72a) above is established in the Appendix in terms of an operator $\mathcal{J}^1(m)$, see (C.13), (C.10a) and (C.12a).

In the absence of any gradient of $y_{[1]}$, see (71), the variational problem for $y_{[1]}$ in (72) is local. This variational problem is solved in the Appendix, §C.8:

* The boundary integral $\Phi_{[2]}^{\text{bt}}$ from Equation (70) yields a stationarity condition applicable on the boundary $\partial\Omega$ of the domain, see Equation (C.31) in the Appendix. This condition does not depend on $y_{[1]}$ and it warrants variational consistency of the input model $\Phi$. We will analyze this condition further in future work.

* The bulk integral $\Phi_{[2]}^{\text{it}}$ from Equation (71) yields a stationarity condition applicable in the interior $\Omega^\circ$ of the domain, that yields the corrector $y^*_{[1]}(X)$ and Lagrange multiplier $g^*_{[1]}(X)$. The solution is of the form announced in (22b),

$$ y^*_{[1]}(X) = (Y_1(m(X)) : \nabla m(X)) \cdot l(X) + Y'_0(m(X)) : \nabla l(X), \quad (73) $$

where the localization tensors $Y_1(m)$ and $Y'_0(m(X))$ are derived in Equation (C.48) in the Appendix.

Inserting this solution into (69–71) yields the functional $\Phi^*_{[2]}[m, l] = \Phi_{[2]}^{\text{it}}[m, l, y^*_{[1]}] + \Phi_{[2]}^{\text{bt}}[m, l, y^*_{[1]}]$ that has been announced in (19).
# 6 Symbolic implementation: The shoal library

We have implemented the method in the symbolic calculation language Wolfram Mathematica (Wolfram Research, Inc. 2021), and we distribute it as an open-source library named shoal (for Second-order HOMogenization Automated using a Library) (Audoly 2023).

The data listed in Table 1 describing the problem at hand is passed to this library, which then returns the tensors listed in Table 3. The extension to rank-deficient problems presented in Appendix E is implemented. The homogenization proceeds by computing the following quantities, in the following order:

* $W_{yl}(m), W_{yy}(m)$ using (B.3b-B.3c),

* $P(m)$ using (B.5), its null vectors $N_P(m)$ in (E.3), its Moore-Penrose inverse $P^†(m)$, as well as $\tilde{I}$ and $\hat{I}$ in (E.8),

* $R(m)$ using (E.19), $Y_0(m), G_0(m)$ using (B.8-B.9), and then $F_0(m), K_0(m)$ and $S_0(m)$ using (B.11), (B.12) and (B.14), see also Table 5,

* $V^l, V^m, V^1$ using (C.2),

* $L(m)$ using (C.5), $L^1(m), L^{11}(m), L^2(m)$ using (C.7), $J^1(m), J^{11}(m), J^2(m)$ using (C.10),

* $A(m)$ using (C.21), and then $A_0(m), K_1(m)$ using (C.23),

* $\hat{B}^{(0)}(m), \hat{B}^{(1)}(m), C^{(0)}(m), C^{(1)}(m)$ using (E.24),

* $\Delta B^{(0)}(m), \Delta B^{(1)}(m)$ using (C.28), and then $B^{(0)}(m), B^{(1)}(m)$ using (C.30),

* $R'(m), Y'(m), G'(m)$ using (C.35-C.37), and then $Y_1(m), Y'_0(m)$ using (C.48a-C.48b),

* $B(m), C(m)$ using (C.45) and (C.39), and then $K_2(m), A_1(m), B_0(m), k_1(m), a_0(m)$ using (C.48c-C.48g),

* if $n_d > 0$, the solvability conditions appearing in (E.14), (E.25) and (E.30). The implementation makes use of standard tensor algebra operations on symbolic tensors, including general transpositions and multiple contractions, as well as symbolic differentiation with respect to $m$, see (C.28). Note that the vector $l$ never appears explicitly in the implementation.

At leading order and in the non-deficient case, the procedure is implemented by the equations listed in Table 5, corresponding to the first three bullet points above. The homogenization at the two following orders makes use of the subsequent bullet points. The special case of uniform properties, when the parameter $m$ is absent, is worked out in Appendix D, see Table D.2 in particular.

<caption>**Table 5** Implementation of the leading-order procedure in the non-deficient case ($n_d = 0$), based on the formulas referenced in the first three items in the bullet list from Section 6.</caption>

$$
\begin{align*}
\mathcal{W}_{yy}(m) &= E_y^T(m) \cdot \mathcal{K}(m) \cdot E_y(m) \\
\mathcal{W}_{yl}(m) &= E_y^T(m) \cdot \mathcal{K}(m) \cdot E_l(m) \\
P(m) &= \begin{pmatrix} \mathcal{W}_{yy}(m) & (Q \cdot E_y(m))^T \\ Q \cdot E_y(m) & \mathbf{0}_{n_c \times n_c} \end{pmatrix} \\
\mathcal{R}(m) &= -P^{-1}(m) \cdot \begin{pmatrix} \mathcal{W}_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix} \\
Y_0(m) &= \begin{pmatrix} I_{n_y} & \mathbf{0}_{n_y \times n_c} \end{pmatrix} \cdot \mathcal{R}(m) \\
G_0(m) &= \begin{pmatrix} \mathbf{0}_{n_c \times n_c} & I_{n_c} \end{pmatrix} \cdot \mathcal{R}(m) \\
F_0(m) &= E_l(m) + E_y(m) \cdot Y_0(m) \\
K_0(m) &= F_0^T(m) \cdot \mathcal{K}(m) \cdot F_0(m) \\
S_0(m) &= \mathcal{K}(m) \cdot F_0(m) + Q^T \cdot G_0(m)
\end{align*}
$$

# 7 Connection with the second-order homogenization of periodic continua

In this section, we show that the extension of the proposed homogenization method to the case of periodic elastic continuum gives similar results to the classical approach (Smyshlyaev and Cherednichenko 2000; Boutin 2019; Durand et al. 2022). Similarly to the discrete example discussed earlier, we consider an infinite medium, so that boundary terms do not matter. We
consider a periodic microstructure: the extension to slowly varying material properties or geometry could be considered as well, and would yield a derivation similar to the one proposed in (Le and Marigo 2018).

<figure><img src="image_7.png" /><figcaption>**Figure 4** A continuum with a periodic microstructure.</figcaption></figure>

It is possible to derive the homogenization for periodic continua from the homogenization results for discrete systems established in the previous sections but this is quite technical. We thus prefer to carry out homogenization from scratch, following the same sequence of steps again. In addition, we do not attempt to link the various quantities relevant to the periodic continua to those relevant to discrete systems in a detailed way: we will content ourselves with formal analogies between them.

## 7.1 Canonical form

We consider a linear elastic continuum of dimension $d = 2$ or $d = 3$ with a periodic microstructure of unit cell $\mathcal{U}$, as illustrated in Figure 4 for $d = 2$. The displacement field in the composite depends on a slow coordinate $X$ and a fast coordinate $x$, with $X = \eta x$ where $\eta$ is a small parameter, see (6), as

$$u(X, x) = U(X) + y(X, x) \qquad (74)$$

where $y(X, x) \in \mathbb{R}^d$ denotes a rapidly fluctuating field, $\mathcal{U}$-periodic with respect to its second argument and subject to a zero-average constraint over one unit-cell,

$$\langle y(X, \cdot) \rangle = \int_{\mathcal{U}} y(X, x) dx = 0, \quad \forall X. \qquad (75)$$

This constraint ensures that $U(X)$ captures the average displacement, $\langle u(X, \cdot) \rangle = U(X)$. We define the macroscopic strain $l$ as the symmetrized gradient of $U(X)$

$$l(X) = \nabla^s U \in \mathbb{T}^{(d,d)}, \qquad (76)$$

where the symbol $\nabla^s$ denotes the symmetrized gradient with respect to the slow coordinate $X$. The microscopic strain is the pair $E = (\varepsilon(X, x), \langle y(X, \cdot) \rangle)$ with $\varepsilon \in \mathbb{T}^{(d,d)}$,

$$\varepsilon(X, x) = l(X) + \nabla^s y(X, x) + \partial^s y(X, x), \qquad (77)$$

where $\partial^s$ denotes the symmetrized gradient with respect to the fast coordinate $x$. This definition is formally similar to (2): we can identify $E_l \cdot l \sim (l(X), 0)$, $E_y \cdot y \sim (\partial^s y(X, x), 0)$ and $E_{y'} \cdot \nabla y \sim (\nabla^s y(X, x), 0)$ and cast the constraint (75) in a form similar to (3): $Q \cdot E(X, x) = \langle y(X, \cdot) \rangle = 0$ $\forall X$, with $Q = (0, 1)$.

The energy is postulated in the form

$$\Phi = \int_{\Omega} W(X) dX \quad \text{where } W(X) = \frac{1}{2|\mathcal{U}|} \int_{\mathcal{U}} \varepsilon(X, x) : \mathbb{C}(x) : \varepsilon(X, x) dx, \qquad (78)$$

and where $\varepsilon(X, x) = (1, 0) \cdot E(X, x)$. In (78), $\mathbb{C}(x) \in \mathbb{T}^{(d,d,d,d)}$ is a $\mathcal{U}$-periodic elasticity tensor (slowly varying properties could be considered here by including an additional dependency on $X$, i.e., by considering $\mathbb{C}(X, x)$, see (Le and Marigo 2018)). Note that no boundary integral is present in the above expression of $\Phi$. Boundary integrals do emerge, however, when one
attempts to justify rigorously the hierarchy in the integrations postulated in (78)—over **x** first and over **X** next—, starting from the elastic energy of a periodic continuum in the limit of scale separation $\eta \to 0$. In addition, boundary terms make it possible to account for the incomplete unit cells located at the boundary of the structure. The analysis of these boundary terms is however beyond the scope of this paper.

## 7.2 Leading order homogenization

At leading order, we neglect slow gradients by setting $\nabla^s \mathbf{y}(\mathbf{X}, \mathbf{x}) = \mathbf{0}$, the energy writes in the form (61) with

$$W^{(0)}(\mathbf{l}, \mathbf{y}) = \frac{1}{2|\mathcal{U}|} \int_{\mathcal{U}} (\mathbf{l} + \partial^s \mathbf{y}(\mathbf{x})) : \mathbb{C}(\mathbf{x}) : (\mathbf{l} + \partial^s \mathbf{y}(\mathbf{x})) \, d\mathbf{x}. \quad (79)$$

For a given symmetric strain tensor $\mathbf{l} \in \mathbb{T}^{(d,d)}$, we seek the stationary point $\mathbf{y}(\mathbf{x}) = \mathbf{y}_{[0]}(\mathbf{x})$ subject to the constraint $\langle \mathbf{y} \rangle = \mathbf{0}$. This writes, using Lagrange multiplier $\mathbf{g} = \mathbf{g}_{[0]}$,

$$\left\{ \begin{array}{l} \langle \mathbf{y} \rangle = \mathbf{0} \\ \int_{\mathcal{U}} (\mathbf{l} + \partial^s \mathbf{y}(\mathbf{x})) : \mathbb{C}(\mathbf{x}) : \partial^s \delta \mathbf{y}(\mathbf{x}) d\mathbf{x} + \mathbf{g} \cdot \langle \delta \mathbf{y} \rangle = \mathbf{0} \quad \forall \delta \mathbf{y} \end{array} \right. \quad (80)$$

where $\delta \mathbf{y}(\mathbf{x})$ is a test function. This problem is similar to (63).

Replacing by a constant $\delta \mathbf{y}$ shows that $\mathbf{g} = \mathbf{0}$. Integration by parts provides the strong form of equilibrium.

$$\begin{array}{l} \mathrm{div}_{\mathbf{x}} \sigma^{(0)}(\mathbf{l}, \mathbf{y}, \mathbf{x}) = \mathbf{0} \quad \text{on } \mathcal{U}, \\ \\ \langle \mathbf{y} \rangle = \mathbf{0} \end{array} \quad (81)$$

where we have defined the microscopic stress

$$\sigma^{(0)}(\mathbf{l}, \mathbf{y}, \mathbf{x}) = \mathbb{C}(\mathbf{x}) : (\mathbf{l} + \partial^s \mathbf{y}(\mathbf{x})). \quad (82)$$

Note that $\sigma^{(0)}$ is automatically periodic in $\mathbf{x}$, which warrants the equilibrium at the boundary of a cell $\mathcal{U}$. In (81), we use the notation $\mathrm{div}_{\mathbf{x}}$ for the divergence with respect to the microscopic (fast) coordinate $\mathbf{x}$.

We denote the solution of (81) as $\mathbf{y} = \mathbf{y}_{[0]}(\mathbf{l}, \mathbf{x})$ and build catalog such that $\mathbf{y}_{[0]}(\mathbf{l}, \mathbf{x}) = Y_0(\mathbf{x}) : \mathbf{l}$. We further define the strain in homogeneous solution as $\varepsilon_{[0]} = F_0(\mathbf{l}, \mathbf{x})$ where

$$F_0(\mathbf{l}, \mathbf{x}) = (\mathbf{l} + \partial^s (Y_0(\mathbf{x}) : \mathbf{l})), \quad (83)$$

as well as the stress $S_0(\mathbf{l}, \mathbf{x}) = \sigma^{(0)}(\mathbf{l}, \mathbf{y}_{[0]}^*, \mathbf{x}) = \mathbb{C}(\mathbf{x}) : F_0(\mathbf{l}, \mathbf{x})$.

The principle of virtual work at dominant order (81) then takes the form

$$\int_{\mathcal{U}} S_0(\mathbf{l}, \mathbf{x}) : \partial^s \delta \mathbf{y}(\mathbf{x}) d\mathbf{x} = \mathbf{0} \quad \forall \delta \mathbf{y}. \quad (84)$$

We eventually obtain the leading order energy of the form (17), (65) as

$$\Phi_{[0]}^*[\mathbf{l}] = \int W_{[0]}^*(\mathbf{l}(\mathbf{X})) d\mathbf{X} \quad \text{where} \quad W_{[0]}^*(\mathbf{l}) = \frac{1}{2|\mathcal{U}|} \int_{\mathcal{U}} F_0(\mathbf{l}, \mathbf{x}) : \mathbb{C}(\mathbf{x}) : F_0(\mathbf{l}, \mathbf{x}) d\mathbf{x}. \quad (85)$$

The expression (85) is formally identical to the leading order energy (35) in (Durand et al. 2022), which has been derived based on the approach of (Smyshlyaev and Cherednichenko 2000). This is akin to standard periodic homogenization (see also Equation (22) in (Le and Marigo 2018)).
## 7.3 Gradient correction

As done earlier in (C.12), we further expand the microscopic strain E order by order:

$$E^{(1)} = (\boldsymbol{\varepsilon}^{(1)}(X, \boldsymbol{x}), \langle \boldsymbol{y}_{[1]}(X, \cdot) \rangle), \quad E^{(2)} = (\boldsymbol{\varepsilon}^{(2)}(X, \boldsymbol{x}), \langle \boldsymbol{y}_{[2]}(X, \cdot) \rangle), \quad (86)$$

with

$$\begin{aligned}
\boldsymbol{\varepsilon}^{(1)}(X, \boldsymbol{x}) &= \nabla^s (Y_0(\boldsymbol{x}) : \boldsymbol{l}(X)) + \partial^s \boldsymbol{y}_{[1]}(X, \boldsymbol{x}) \\
\boldsymbol{\varepsilon}^{(2)}(X, \boldsymbol{x}) &= \nabla^s \boldsymbol{y}_{[1]}(X, \boldsymbol{x}) + \partial^s \boldsymbol{y}_{[2]}(X, \boldsymbol{x}).
\end{aligned}
\quad (87)$$

Following (C.17), (C.18), we eliminate $\boldsymbol{y}_{[1]}$ from the energy at order 1 using (84) and obtain

$$\Phi_{[1]} = \int_{\Omega} S_0(l(X), \boldsymbol{x}) : \nabla^s (Y_0(\boldsymbol{x}) : l(X)) dX. \quad (88)$$

In (Durand et al. 2022), odd-order tensors are considered to be zero due to centro-symmetry assumption: here we work in a more general setting. A similar order-1 energy contribution is derived in (Le and Marigo 2018), see their equation (36), with an additional contribution coming from the macroscopic gradient of elastic properties.

The energy at order 2 further writes, see (C.17),

$$\Phi_{[2]} = \int_{\Omega} \frac{1}{|\mathcal{U}|} \int_{\mathcal{U}} \left( \frac{1}{2} \boldsymbol{\varepsilon}^{(1)}(X, \boldsymbol{x}) : \mathbb{C}(\boldsymbol{x}) : \boldsymbol{\varepsilon}^{(1)}(X, \boldsymbol{x}) + S_0(l(X), \boldsymbol{x}) : \boldsymbol{\varepsilon}^{(2)}(X, \boldsymbol{x}) \right) d\boldsymbol{x} dX. \quad (89)$$

Using relation (84), we can eliminate $\boldsymbol{y}_{[2]}$ from the second term

$$\int_{\mathcal{U}} S_0(l(X), \boldsymbol{x}) : \boldsymbol{\varepsilon}^{(2)}(X, \boldsymbol{x}) d\boldsymbol{x} = \int_{\mathcal{U}} S_0(l(X), \boldsymbol{x}) : \nabla^s \boldsymbol{y}_{[1]}(X, \boldsymbol{x}) d\boldsymbol{x}. \quad (90)$$

Next, we identify the terms (C.26) and (C.27) that require an integration by parts

$$\begin{aligned}
\int_{\Omega} (\mathcal{C}^{(0)} : (\boldsymbol{h}(X) \otimes \boldsymbol{h}(X))) : \nabla^2 \boldsymbol{h}(X) dX &= 0, \\
\int_{\Omega} (\mathcal{C}^{(1)} \cdot \boldsymbol{h}(X)) : \nabla \boldsymbol{y}_{[1]}(X) dX &= \int_{\Omega} \frac{1}{|\mathcal{U}|} \int_{\mathcal{U}} S_0(l(X), \boldsymbol{x}) : \nabla^s \boldsymbol{y}_{[1]}(X, \boldsymbol{x}) d\boldsymbol{x} dX.
\end{aligned}
\quad (91)$$

After integration by parts, the energy (89) writes

$$\begin{aligned}
\Phi_{[2]} = \int_{\Omega} \frac{1}{|\mathcal{U}|} \int_{\mathcal{U}} \left( \frac{1}{2} \boldsymbol{\varepsilon}^{(1)}(X, \boldsymbol{x}) : \mathbb{C}(\boldsymbol{x}) : \boldsymbol{\varepsilon}^{(1)}(X, \boldsymbol{x}) \right) d\boldsymbol{x} dX & \\
+ \oint_{\partial\Omega} \left( \frac{1}{|\mathcal{U}|} \int_{\mathcal{U}} S_0(l(X), \boldsymbol{x}) \cdot \boldsymbol{y}_{[1]}(X, \boldsymbol{x}) d\boldsymbol{x} \right) \cdot \boldsymbol{n} d\boldsymbol{a}, &
\end{aligned}
\quad (92)$$

where $\boldsymbol{n}$ is the unit normal to the domain boundary. In (92), we use the notation $\operatorname{div}_X$ for the divergence with respect to the slow coordinate $X$. The domain being infinite, we can ignore the boundary term and the associated stationarity condition on the boundary.

The corrector $\boldsymbol{y} = \boldsymbol{y}_{[1]}(\nabla l, \boldsymbol{x})$ is eventually found as a solution to the local variational problem, similar to (72),

$$\left\{
\begin{aligned}
\langle \boldsymbol{y} \rangle &= 0 \\
\int_{\mathcal{U}} (\boldsymbol{\sigma}^{(1)} \cdot \partial^s \delta \boldsymbol{y} - \operatorname{div}_X S_0(l, \boldsymbol{x}) \cdot \delta \boldsymbol{y}) d\boldsymbol{x} + \boldsymbol{g}_{[1]} \cdot \langle \delta \boldsymbol{y} \rangle &= 0 \quad \forall \delta \boldsymbol{y},
\end{aligned}
\right.
\quad (93)$$

where we have identified the first order correction to the microscopic stress $\boldsymbol{\sigma}^{(1)}(\nabla l, l, \boldsymbol{y}, \boldsymbol{x}) = \mathbb{C}(\boldsymbol{x}) : \nabla^s (Y_0(\boldsymbol{x}) : l) + \partial^s \boldsymbol{y}(\boldsymbol{x})$. After integrating by parts, we obtain

$$\begin{aligned}
\operatorname{div}_X \boldsymbol{\sigma}^{(1)}(l, \boldsymbol{y}, \boldsymbol{x}) - \boldsymbol{g}_{[1]} &= - \operatorname{div}_X S_0(l, \boldsymbol{x}) \quad \text{on } \mathcal{U}, \\
\langle \boldsymbol{y} \rangle &= 0
\end{aligned}
\quad (94)$$
The microscopic stress $\sigma^{(1)}$ is automatically periodic in $x$, which warrants the equilibrium at the boundary of a cell $U$. This microscopic problem (94) is similar to equation (16) in (Durand et al. 2022).

The correction $\mathbf{y}^*_{[1]}(\nabla l, x)$ to the microscopic degrees of freedom writes, following (73), $\mathbf{y}^*_{[1]}(\nabla l, x) = Y'_0(x) : \nabla l$. We can thus write the first order correction to strain as $\varepsilon_{[1]} = F_1(l, \nabla l, x)$ where

$$F_1(l, \nabla l, x) = \nabla^s (Y_0(x) : l) + \partial^s (Y'_0(x) : \nabla l). \quad (95)$$

Inserting into (92) allows us to identify the second order contribution to the homogenized energy in the form (19) as

$$\Phi^*_{[2]}[l, \nabla l] = \int_{\Omega} B_0 :: \frac{\nabla l \otimes \nabla l}{2} dX, \quad (96)$$

with

$$B_0 :: \frac{\nabla l \otimes \nabla l}{2} = \frac{1}{|U|} \int_{U} \left( \begin{array}{l} \frac{1}{2} F_1(l, \nabla l, x) : \mathbb{C}(x) : F_1(l, \nabla l, x) \\ - \operatorname{div}_X S_0(l, x) \cdot (Y'_0(x) : \nabla l) \end{array} \right) dx. \quad (97)$$

Boundary integrals have been removed in (96) due to the fact that we consider an infinite domain. This result is similar to equation (35) in (Durand et al. 2022) and equation (50) in (Le and Marigo 2018) (this latter work reports additional terms that capture the effect of a gradient of elastic properties).

Note that it is possible to identify the tensors $K_0$, $A_0$ and $B_0$ appearing in (13–14), by factoring out $l$, $\nabla l$ and $n$ in equations (85), (88), (96) and (97). However, this introduces cumbersome expressions and we prefer to ignore this step.

# 8 Discussion and conclusion

We have proposed an asymptotically exact, second-order homogenization procedure for linear, discrete elastic structures, such as elastic trusses or networks of elastic beams. Our homogenization method works at the energy level, see Equations (10) and (11), and is similar in this respect to the approach to dimension reduction developed by (Berdichevskii 1981; Hodges 2006; Lestringant and Audoly 2020). It uses an abstract energy formulation as a starting point, see Section 2, and as a result can equally cover two-dimensional or three-dimensional lattices made up of beams or springs. It is designed to be generic, and addresses the case of pre-stress or pre-strain as well as slowly modulated elastic or geometric properties—in their work on the homogenization of periodic continua, (Le and Marigo 2018) address the case of slowly modulated elastic properties but assume that the geometry of the unit cell is invariant. The method can also account for kinematic constraints, and can be readily applied to, e.g., a lattice made up of inextensible beams, see Appendix E. Besides, the connection with existing approaches on the homogenization of periodic continua has been pointed out in Section 7.

The homogenization procedure involves a series of linear algebra calculations that have been implemented once for all in the form of an open-source library, named shoal and based on the symbolic calculation language Wolfram Mathematica. This work will serve as a foundation for a series of applications which we will cover in follow-up papers. In one particular follow-up paper, we will analyze more advanced beam lattices than those treated in the illustration section (Section 4).

The uniqueness of the solutions $\mathbf{y}^*_{[0]}$ and $\mathbf{y}^*_{[1]}$ to the variational problems (B.4) and (C.33) follows from the assumption (12) that $\mathbf{W}_{yy}(\mathbf{m})$ is positive definite on the subspace of kinematically admissible microscopic degrees of freedom $\mathbf{y}$ defined by $Q \cdot E_y(\mathbf{m}(X)) \cdot \mathbf{y} = 0$. In these circumstances, $\mathbf{y}^*_{[0]}$ and $\mathbf{y}^*_{[1]}$ do not only make the strain energy stationary, as mentioned in the paper, but also minimum.

Classical work on homogenization postulates the displacement in the form of an expansion $\mathbf{u}(X, x) = \sum_i \eta^i v_i(X, x)$, where $(X, x)$ are the slow and fast variables, respectively, and identifies
a potential $\tilde{\Phi}^*[\langle v_0 \rangle, \langle v_1 \rangle, \dots]$ depending separately on each one of the macroscopic averages $\langle v_i \rangle = \langle v_i(X, \cdot) \rangle$, see (Le and Marigo 2018) for instance. As a result, the stationarity condition for $\tilde{\Phi}^*$ takes the form of separate problems for $\langle v_0 \rangle, \langle v_1 \rangle$, etc. It has apparently not been appreciated that the homogenized energy has to be a functional $\Phi^*[\langle v_0 \rangle + \eta \langle v_1 \rangle + \eta^2 \langle v_2 \rangle + \dots]$ of the total macroscopic displacement $\langle u \rangle = \langle v_0 \rangle + \eta \langle v_1 \rangle + \eta^2 \langle v_2 \rangle + \dots$: this warrants that its value is unaffected by a re-parameterization of $\eta$ leaving the physical displacement $u$ unchanged. Our approach works differently: with the notation of the continuous case from Section 7, we postulate $u(X, x) = U(X) + \sum_i \eta^i y_i(X, x)$, subject to the constraint $\langle y_i(X, \cdot) \rangle = 0$ (for any $X$ and any $i$), see (74-75), compute the strain $l$ in terms of the (total) macroscopic displacement $U = \langle u \rangle$, and obtain a homogenized energy $\Phi^*[l]$ depending on the total macroscopic strain $l$ (and not on each of the various contributions $l_i$ to the strain separately). This is simpler than the traditional approach, and ultimately equivalent.

Among the perspectives opened up by the present work, we can mention the extension to non-linear elastic structures (which can be addressed by adapting our previous work on non-linear dimension reduction (Lestringant and Audoly 2020)) and a careful treatment of the boundary layers (which we could incorporate into the homogenized model by means of effective boundary terms, along the lines of what has been done by (David et al. 2012)).

# Appendix

## A Tensor algebra

The dimension of the Euclidean space is denoted as $d$. The Euclidean space is endowed with an orthonormal Cartesian basis ($e_1, \dots, e_d$). A generic point in the Euclidean space is denoted as $X \in \mathbb{R}^d$.

We denote as $\mathbb{T}^{(n_1, n_2, \dots, n_p)}$ the tensor space $\mathbb{T}^{(n_1, n_2, \dots, n_p)} = \mathbb{R}^{n_1} \otimes \mathbb{R}^{n_2} \otimes \dots \otimes \mathbb{R}^{n_p}$ made of tensors $R$ of rank $p$ and dimensions $n_1 \times n_2 \times \dots \times n_p$. In particular, $I_k \in \mathbb{T}^{(k,k)}$ denotes the identity matrix in dimension $k$ and $0_{n_1 \times \dots \times n_p} \in \mathbb{T}^{(n_1, n_2, \dots, n_p)}$ the null tensor with dimensions $n_1 \times \dots \times n_p$.

Tensors and vectors are denoted using bold symbols, while scalars (including tensor components) are denoted using non-bold symbols.

Given two tensors $R \in \mathbb{T}^{(n_1, n_2, \dots, n_p)}$ and $R' \in \mathbb{T}^{(n'_1, n'_2, \dots, n'_p)}$, we denote as

* $R \cdot R' \in \mathbb{T}^{(n_1, n_2, \dots, n_{p-1}, n'_1, \dots, n'_p)}$ their simple contraction (whose existence requires $n_p = n'_1$),

* $R : R' \in \mathbb{T}^{(n_1, n_2, \dots, n_{p-2}, n'_3, \dots, n'_p)}$ their double contraction (whose existence requires $n_{p-1} = n'_1$ and $n_p = n'_2$),

* $R \cdot R' \in \mathbb{T}^{(n_1, n_2, \dots, n_{p-3}, n'_4, \dots, n'_p)}$ their triple contraction, (whose existence requires $n_{p-2} = n'_1$, $n_{p-1} = n'_2$ and $n_p = n'_3$),

* etc.

If they exist, the contracted tensors are given by

$$ (R \cdot R')_{i_1 \dots i_{p-1} i'_2 \dots i'_p} = R_{i_1 \dots i_{p-1} j} R'_{j i'_2 \dots i'_p} \\ (R : R')_{i_1 \dots i_{p-2} i'_3 \dots i'_p} = R_{i_1 \dots i_{p-1} j k} R'_{j k i'_2 \dots i'_p} \qquad (A.1) $$

$$ (R : R')_{i_1 \dots i_{p-3} i'_4 \dots i'_p} = R_{i_1 \dots i_{p-1} j k l} R'_{j k l i'_2 \dots i'_p} $$

Note the ordering of the contracted indices $j, k, l$, etc. in the right-hand sides—in particular, the double contraction of two rank-2 tensors $A$ and $B$ is given in our notation by $A : B = \text{tr}(A \cdot B^T)$. Here and elsewhere in the paper, we use Einstein summation whereby any index that is repeated on one side of the equal sign is implicitly summed.

The action of a matrix $R$ on a vector $v$ is viewed as a special case of the contraction of a tensor of rank 2 with a tensor of rank 1, and is denoted as $R \cdot v$, with a dot.

The outer product of two tensors $T$ and $T'$ is denoted as $T \otimes T'$. In particular, the outer product of two vectors is denoted as $v \otimes v'$. Vector transposition is not a meaningful operation in our notation.
Given a tensor $R \in \mathbb{T}^{(n_1, n_2, ..., n_p)}$ and a permutation $(\sigma_1, ..., \sigma_p)$ of the levels $(1, ..., p)$ of the tensor, we denote as $R^{T_{\sigma_1}...\sigma_p}$ the generalized transpose of $R$, such that the level $i$ in the original tensor becomes level $\sigma_i$ in the transpose:

$$ (R^{T_{\sigma_1}...\sigma_p})_{i_1...i_p} = R_{i_1...i_p} \qquad (A.2) $$

For a tensor of rank $p = 4$ and the permutation $(\sigma_1, \sigma_2, \sigma_3, \sigma_4) = (1, 3, 4, 2)$, for instance, we have $(R^{T_{1342}})_{ijkl} = R_{ijkl}$.

Transposing will allow us to reorder the indices of a tensor in any desired order. Suppose for instance that we wish to rewrite an expression $R_{ijkl}$ as the component $R'_{ijkl}$ of another tensor whose indices must appear in alphabetical order: $R'$ is clearly a transpose of $R$, and the permutation is found by noting that the levels (1, 2, 3, 4) in the original tensor $R$, corresponding to the indices $(i, k, l, j)$, become respectively the levels (1, 3, 4, 2) = $(\sigma_1, \sigma_2, \sigma_3, \sigma_4)$ in $R'$. This yields

$$ R_{ijkl} = (R^{T_{1342}})_{ijkl}. \qquad (A.3) $$

Index reordering using transposition will be routinely used in combination with contractions to remove indices in tensor algebra, as in $R_{ijkl}R'_{ijkl} = R^{T_{1324}} :: R'$.

The transpose $R^T$ of a matrix $R \in \mathbb{T}^{(n_1, n_2)}$ is a particular case of the generalized transpose, $R^T = R^{T_{21}}$.

The symmetrization of a tensor $R$ with respect to a pair of indices $(i, j)$ is denoted as $R^{S_{ij}}$. For a tensor $R$ of rank $p = 4$, for instance, the symmetrization with respect to the first and third indices is given by

$$ R^{S_{13}} = \frac{1}{2}(R + R^{T_{3214}}), \qquad (A.4) $$

where $R^{T_{3214}}$ is obtained from $R$ by permuting the first and third levels.

A tensor invariant by a permutation of its levels $i$ and $j$ will be said to satisfy the $S_{ij}$ symmetry; for instance, $R^{S_{13}}$ is $S_{13}$ symmetric, by construction. More generally, a tensor will be said to satisfy the $S_{\{ij\}\{kl\}}$ symmetry if it is symmetric by the combined permutation of indices $i \leftrightarrow k$ and $j \leftrightarrow l$. For instance,

$$ R \text{ is } S_{\{23\}\{45\}} \text{ symmetric} \Leftrightarrow R_{ijklm} = R_{ilmjk}. \qquad (A.5) $$

The symmetrization with respect to pairs of indices works similarly,

$$ R^{S_{\{23\}\{45\}}} = \frac{1}{2}(R + R^{T_{14523}}). \qquad (A.6) $$

Clearly, $R$ is $S_{\{23\}\{45\}}$ symmetric if and only if $R^{S_{\{23\}\{45\}}} = R$.

The composition of symmetrizations is represented by the symbol $\circ$. In Equation (C.7), for instance, it stands for

$$ R^{S_{23} \circ S_{\{45\}\{67\}}} = (R^{S_{\{45\}\{67\}}})^{S_{23}}. \qquad (A.7) $$

Given a tensor $R(q) \in \mathbb{T}^{(n_1, n_2, ..., n_p)}$ taking an argument $q \in \mathbb{R}^{n_q}$, we denote as $dR/dq \in \mathbb{T}^{(n_1, n_2, ..., n_p, n_q)}$ its gradient,

$$ \left(\frac{dR}{dq}\right)_{i_1...i_p, j} = \frac{\partial R_{i_1...i_p}}{\partial q_j}. \qquad (A.8) $$

By a standard convention, the index $j$ corresponding to differentiation appears last in the gradient. When the parameter $q$ coincides with the spatial variable $X \in \mathbb{R}^d$, we use the nabla notation,

$$ \nabla R = \frac{dR}{dX}. \qquad (A.9) $$

The alternate notation $R\nabla$ has the advantage of respecting the order of indices but is also less standard.
## B Detailed analysis of leading order (classical homogenization)

Inserting the microscopic strain $E^{(0)}(m, l, y)$ given in (59) into the expression (60) of the strain energy density $W^{(0)}(h, y)$, we have

$$W^{(0)}(m, l, y) = W(m, E_l(m) \cdot l + E_y(m) \cdot y) \quad (B.1)$$

Using the expression of $W$ in (5) and expanding, we rewrite this in block-matrix notation as

$$W^{(0)}(m, l, y) = \frac{1}{2} \begin{pmatrix} l \\ y \end{pmatrix} \cdot W(m) \cdot \begin{pmatrix} l \\ y \end{pmatrix}, \quad \text{where} \quad W(m) = \begin{pmatrix} W_{ll}(m) & W_{yl}^T(m) \\ W_{yl}(m) & W_{yy}(m) \end{pmatrix}, \quad (B.2)$$

and the tensors $W(m) \in \mathbb{T}^{(n_l+n_y, n_l+n_y)}$, $W_{ll}(m) \in \mathbb{T}^{(n_l, n_l)}$, $W_{yl}(m) \in \mathbb{T}^{(n_y, n_l)}$ and $W_{yy}(m) \in \mathbb{T}^{(n_y, n_y)}$ are given by

$$W_{ll}(m) = E_l^T(m) \cdot K(m) \cdot E_l(m) \quad (B.3a)$$

$$W_{yl}(m) = E_y^T(m) \cdot K(m) \cdot E_l(m) \quad (B.3b)$$

$$W_{yy}(m) = E_y^T(m) \cdot K(m) \cdot E_y(m). \quad (B.3c)$$

Using (B.2), the optimality condition (63) for $y_{[0]}^*$ and $g_{[0]}^*$ takes the form

$$P(m(X)) \cdot \begin{pmatrix} y_{[0]}^*(X) \\ g_{[0]}^*(X) \end{pmatrix} + \begin{pmatrix} W_{yl}(m(X)) \\ Q \cdot E_l(m(X)) \end{pmatrix} \cdot l(X) = 0. \quad (B.4)$$

where

$$P(m) = \begin{pmatrix} W_{yy}(m) & (Q \cdot E_y(m))^T \\ Q \cdot E_y(m) & \mathbf{0}_{n_c \times n_c} \end{pmatrix} \in \mathbb{T}^{(n_y+n_c, n_y+n_c)}. \quad (B.5)$$

We focus on the case where $P(m)$ is invertible: the non-invertible (rank-deficient) case is treated in Appendix E. The solution $(y_{[0]}^*(X), g_{[0]}^*(X))$ is then found by inverting this linear system as

$$\begin{pmatrix} y_{[0]}^*(X) \\ g_{[0]}^*(X) \end{pmatrix} = \mathcal{R}(m(X)) \cdot l(X) \quad (B.6)$$

where

$$\mathcal{R}(m) = -P^{-1}(m) \cdot \begin{pmatrix} W_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix}. \quad (B.7)$$

In the code, we implemented a more general expression of $\mathcal{R}$ that applies to rank-deficient matrices, see Appendix E and Equation (E.19) in particular.

Equation (B.6) matches the form of the solution $y_{[0]}^*(X) = Y_0(m) \cdot l$ announced in (64) and the localization tensor is identified as

$$Y_0(m) = \begin{pmatrix} I_{n_y} & \mathbf{0}_{n_y \times n_c} \end{pmatrix} \cdot \mathcal{R}(m). \quad (B.8)$$

The solution for the Lagrange multipliers is $g_{[0]}^*(X) = G_0(m) \cdot l$ where

$$G_0(m) = \begin{pmatrix} \mathbf{0}_{n_c \times n_y} & I_{n_c} \end{pmatrix} \cdot \mathcal{R}(m). \quad (B.9)$$
We proceed to introduce important additional quantities that characterize the leading-order solution.

The strain $E^{[0]}(m, l) = E^{(0)}(m, l, y = Y_0(m) \cdot l)$ is found using (59) as

$$E^{[0]}(m, l) = F_0(m) \cdot l, \qquad (B.10)$$

where the strain localization tensor $F_0(m)$ is given by

$$F_0(m) = E_l(m) + E_y(m) \cdot Y_0(m). \qquad (B.11)$$

The leading-order strain energy $W^{\star}_{[0]}(m, l) = W^{(0)}(m, l, y = Y_0(m) \cdot l)$ can then be written with the help of (B.2-B.3) and (B.11) in a form that matches that announced in (65), namely $W^{\star}_{[0]}(m, l) = \frac{1}{2} l \cdot K_0(m) \cdot l$, where the elasticity tensor $K_0(m)$ characterizing the equivalent Cauchy-type elastic continuum at order $\eta^0$ is identified as

$$K_0(m) = F_0^T(m) \cdot K(m) \cdot F_0(m). \qquad (B.12)$$

To complete the analysis of solutions at order $\eta^0$, we derive a useful identity that will help simplify the higher orders in the energy expansion. Inserting the solution $y^{\star}_{[0]}(X) = Y_0(m) \cdot l$ and $g^{\star}_{[0]}(X) = G_0(m) \cdot l$ into the stationarity condition (B.4) and using (B.3) and identifying $F_0$ from (B.11), we get

$$E_y^T(m) \cdot (K(m) \cdot F_0(m) + Q^T \cdot G_0(m)) \cdot l = 0 \qquad (B.13a)$$

$$Q \cdot E^{[0]}(m, l) = 0. \qquad (B.13b)$$

The quantity $(K(m) \cdot F_0(m) + Q^T \cdot G_0(m)) \cdot l$ appearing in the first equation can be identified as the leading-order stress, consisting of the elastic stress $K(m) \cdot F_0(m) \cdot l = K(m) \cdot E^{[0]}$ plus the stress $Q^T \cdot G_0(m) \cdot l$ enforcing the constraint. We therefore introduce the stress localization tensor as

$$S_0(m) = K(m) \cdot F_0(m) + Q^T \cdot G_0(m), \qquad (B.14)$$

and rewrite Equation (B.13), after simplification by the arbitrary factor $l$, as

$$E_y^T(m) \cdot S_0(m) = 0 \qquad (B.15a)$$

$$Q \cdot F_0(m) = 0. \qquad (B.15b)$$

Equation (B.15a) is the principle of virtual work at leading order: multiplying by a virtual displacement $\delta y$ on the left-hand side and by $l$ on the right-hand side, and rearranging, it takes the usual form $(S_0(m) \cdot l) \cdot (E_y(m) \cdot \delta y) = 0$, where the left-hand side is the stress contracted with the virtual increment of strain.

# C Detailed analysis of the gradient effect

## C.1 Packed macroscopic variables

For the analysis of the gradient effect, it is convenient to introduce the quantity $h(X)$ obtained by concatenating the microscopic variables $l$ and $m$ together with a trailing 1: in block-vector notation,

$$h(X) = (l(X) \quad m(X) \quad 1) \in \mathbb{R}^{n_h} \quad \text{where } n_h = n_l + n_m + 1. \qquad (C.1)$$

With the help of the matrices $V^l \in T^{(n_l, n_h)}$, $V^m \in T^{(n_m, n_h)}$ and the vector $V^1 \in \mathbb{R}^{n_h}$ defined in block-matrix notation by

$$V^l = (I_{n_l} \quad 0_{n_l \times (n_m+1)}), \quad V^m = (0_{n_m \times n_l} \quad I_{n_m} \quad 0_{n_m \times 1}), \quad V^1 = (0_{(n_l+n_m)} \quad 1), \qquad (C.2)$$
<table><caption>**Table C.1** Summary of the tensors used internally by the homogenization procedure in Appendix C. All these tensors make use of the compact *h* notation, see (C.1).</caption><thead><tr><th></th><th>tensor space</th><th>symmetry</th><th>usage</th><th>content</th></tr></thead><tbody><tr><td>L<sup>(m)</sup></td><td>T(n<sub>y</sub>,n<sub>h</sub>)</td><td>-</td><td>y<sup>*</sup><sub>[0]</sub> = L · h</td><td>l</td></tr><tr><td>L<sup>1</sup>(m)</td><td>T(n<sub>y</sub>,d,n<sub>h</sub>,d,n<sub>h</sub>)</td><td>-</td><td>∇y<sup>*</sup><sub>[0]</sub> = (L<sup>1</sup> · h) : ∇h</td><td>(l ⊗ ∇m, ∇l)</td></tr><tr><td>L<sup>11</sup>(m)</td><td>T(n<sub>y</sub>,d,d,n<sub>h</sub>,d,n<sub>h</sub>,d,n<sub>h</sub>)</td><td>S<sub>23</sub>, S<sub>{45}</sub>{67}</td><td>∇<sup>2</sup>y<sup>*</sup><sub>[0]</sub> = (L<sup>11</sup> · h) :: (∇h ⊗ ∇h) + ...</td><td>(l ⊗ ∇m, ∇l) ⊗ ∇m</td></tr><tr><td>L<sup>2</sup>(m)</td><td>T(n<sub>y</sub>,d,d,n<sub>h</sub>,d,d,n<sub>h</sub>)</td><td>S<sub>23</sub>, S<sub>56</sub></td><td>∇<sup>2</sup>y<sup>*</sup><sub>[0]</sub> = (L<sup>2</sup> · h) :: ∇<sup>2</sup>h + ...</td><td>(∇<sup>2</sup>l, l ⊗ ∇<sup>2</sup>m)</td></tr><tr><td>J<sup>1</sup>(m)</td><td>T(n<sub>E</sub>,n<sub>h</sub>,d,d,n<sub>h</sub>)</td><td>-</td><td>E<sup>1</sup> = (J<sup>1</sup> · h) : ∇h + ...</td><td>(l ⊗ ∇m, ∇l)</td></tr><tr><td>J<sup>11</sup>(m)</td><td>T(n<sub>E</sub>,n<sub>h</sub>,d,d,n<sub>h</sub>,d,n<sub>h</sub>)</td><td>S<sub>{23}</sub>{45}</td><td>E<sup>2</sup> = (J<sup>11</sup> · h) :: (∇h ⊗ ∇h) + ...</td><td>(l ⊗ ∇m, ∇l) ⊗ ∇m</td></tr><tr><td>J<sup>2</sup>(m)</td><td>T(n<sub>E</sub>,n<sub>h</sub>,d,d,n<sub>h</sub>)</td><td>S<sub>34</sub></td><td>E<sup>2</sup> = (J<sup>2</sup> · h) :: ∇<sup>2</sup>h + ...</td><td>(∇<sup>2</sup>l, l ⊗ ∇<sup>2</sup>m)</td></tr><tr><td>A<sup>(m)</sup></td><td>T(n<sub>h</sub>,d,n<sub>h</sub>,n<sub>h</sub>)</td><td>S<sub>34</sub></td><td>Φ<sub>[1]</sub> = ∫<sub>Ω</sub> (A : (h ⊗ h)) : ∇hdX</td><td>(l ⊗ l ⊗ ∇m, l ⊗ ∇l)</td></tr><tr><td>B̂<sup>(0)</sup>(m)</td><td>T(n<sub>h</sub>,d,n<sub>h</sub>,d,n<sub>h</sub>,n<sub>h</sub>)</td><td>(delayed)</td><td>Φ<sub>[2]</sub> = ∫<sub>Ω</sub> (B<sup>(0)</sup> : (h ⊗ h)) :: <sup>∇h⊗∇h</sup>⁄<sub>2</sub> dX + ...</td><td>(l ⊗ ∇m, ∇l) ⊗<sup>2</sup></td></tr><tr><td>ΔB<sup>(0)</sup>(m)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>B<sup>(0)</sup>(m)</td><td>beginequation2opt[ΔB<sup>(0)</sup>(m)]<br>B<sup>(1)</sup>(m)<br>B̂<sup>(1)</sup>(m)</td><td>-</td><td>Φ<sub>[2]</sub> = ∫<sub>Ω</sub> (B<sup>(1)</sup> · h) :: (∇h ⊗ y<sub>1</sub>)dX + ...</td><td>(l ⊗ ∇m, ∇l) ⊗ y<sub>1</sub></td></tr><tr><td>C<sup>(0)</sup>(m)</td><td>T(n<sub>h</sub>,d,d,n<sub>h</sub>,n<sub>h</sub>)</td><td>S<sub>23</sub>, S<sub>45</sub></td><td>Φ<sub>[2]</sub> = ∫<sub>Ω</sub> (C<sup>(0)</sup> : (h ⊗ h)) :: ∇<sup>2</sup>hdX + ...</td><td>l ⊗ (l ⊗ ∇<sup>2</sup>m, ∇<sup>2</sup>l)</td></tr><tr><td>C<sup>(1)</sup>(m)</td><td>T(n<sub>y</sub>,d,d,n<sub>h</sub>)</td><td>-</td><td>Φ<sub>[2]</sub> = ∫<sub>Ω</sub> (C<sup>(1)</sup> · h) : ∇y<sub>1</sub>dX + ...</td><td>l ⊗ ∇y<sub>1</sub></td></tr><tr><td>B<sup>(m)</sup></td><td>T(n<sub>h</sub>,d,n<sub>h</sub>,d,n<sub>h</sub>,n<sub>h</sub>)</td><td>S<sub>{12}</sub>{34}, S<sub>56</sub></td><td>Φ<sup>*</sup><sub>[2]</sub> = ∫<sub>Ω</sub> (B<sup>(0)</sup> : (h ⊗ h)) :: <sup>∇h⊗∇h</sup>⁄<sub>2</sub> dX + ...</td><td>(l ⊗ ∇m, ∇l) ⊗<sup>2</sup></td></tr><tr><td>C<sup>(m)</sup></td><td>T(n<sub>h</sub>,d,d,n<sub>h</sub>,n<sub>h</sub>)</td><td>S<sub>45</sub></td><td>Φ<sup>*</sup><sub>[2]</sub> = ... + ∮<sub>∂Ω</sub> (C : (h ⊗ h)) :: (∇h ⊗ n)da</td><td>l ⊗ (l ⊗ ∇m, ∇l) ⊗ n</td></tr><tr><td>Y<sup>(m)</sup></td><td>T(n<sub>y</sub>,d,d,n<sub>h</sub>)</td><td>-</td><td>y<sup>*</sup><sub>[1]</sub> = (Y' · h) : ∇h</td><td>(l ⊗ ∇m, ∇l)</td></tr><tr><td>G<sup>'</sup>(m)</td><td>T(n<sub>c</sub>,d,d,n<sub>h</sub>)</td><td>-</td><td>g<sup>*</sup><sub>[1]</sub> = (G' · h) : ∇h</td><td>(l ⊗ ∇m, ∇l)</td></tr></tbody></table>
on can rewrite the definition of h in (C.1) as

$$h(X) = (\mathcal{V}^l)^T \cdot l(X) + (\mathcal{V}^m)^T \cdot m(X) + \mathcal{V}^1. \quad (C.3)$$

The converse (unpacking) operation is implemented as

$$l(X) = \mathcal{V}^l \cdot h(X), \quad m(X) = \mathcal{V}^m \cdot h(X). \quad (C.4)$$

## C.2 Structure coefficients

The leading-order prediction (22a) for the microscopic degrees of freedom, $\mathbf{y}^*_{[0]} = Y_0(\mathbf{m}(X)) \cdot l(X)$ can be expressed in terms of $h(X)$ with the help of (C.4) as

$$\mathbf{y}^*_{[0]} = \mathcal{L}(\mathbf{m}(X)) \cdot h(X) \quad \text{where} \quad \mathcal{L}(\mathbf{m}) = Y_0(\mathbf{m}) \cdot \mathcal{V}^l. \quad (C.5)$$

The successive gradients of $\mathbf{y}(X)$ in (66) can then be calculated as

$$\begin{align}
\mathbf{y}(X) &= \mathcal{L}(\mathbf{m}(X)) \cdot h(X) + \mathbf{y}_{[1]}(X) + \mathbf{y}_{[2]}(X) + \mathcal{O}(\eta^3) \nonumber \\
\nabla \mathbf{y}(X) &= (\mathcal{L}^1(\mathbf{m}(X)) \cdot h(X)) : \nabla h(X) + \nabla \mathbf{y}_{[1]}(X) + \mathcal{O}(\eta^3) \nonumber \\
\nabla^2 \mathbf{y}(X) &= (\mathcal{L}^{11}(\mathbf{m}(X)) \cdot h(X)) :: (\nabla h(X) \otimes \nabla h(X)) \nonumber \\
&\quad + (\mathcal{L}^2(\mathbf{m}(X)) \cdot h(X)) :: \nabla^2 h(X) + \mathcal{O}(\eta^3) \tag{C.6}
\end{align}$$

where the symbol $\mathcal{O}(\eta^3)$ stands for terms of order $\eta^3$ and higher, such as $\nabla^2 \mathbf{y}_{[1]} = \mathcal{O}(\eta^{2+1})$ and $\nabla \mathbf{y}_{[2]} = \mathcal{O}(\eta^{1+2})$.

By design, the tensor $\mathcal{L}^1(\mathbf{m})$, $\mathcal{L}^{11}(\mathbf{m})$ and $\mathcal{L}^2(\mathbf{m})$ capture the successive gradients of $\mathbf{y}^*_{[0]}$ as $\nabla \mathbf{y}^*_{[0]} = (\mathcal{L}^1(\mathbf{m}) \cdot h) : \nabla h$ and $\nabla^2 \mathbf{y}^*_{[0]} = (\mathcal{L}^{11}(\mathbf{m}) \cdot h) :: (\nabla h \otimes \nabla h) + (\mathcal{L}^2(\mathbf{m}) \cdot h) :: (\nabla^2 h)$. They are identified by differentiating (C.5) with respect to $X$, which yields

$$\mathcal{L}^1(\mathbf{m}) = \left( \frac{d\mathcal{L}}{dm} \cdot \mathcal{V}^m \otimes I_d \right)^{T_{15324}} + (\mathcal{L}(\mathbf{m}) \otimes \mathcal{V}^1 \otimes I_d)^{T_{13524}} \quad (C.7a)$$

$$\mathcal{L}^{11}(\mathbf{m}) = \left[ \left( \frac{d\mathcal{L}^1}{dm} \cdot \mathcal{V}^m \otimes I_d \right)^{T_{12678435}} + (\mathcal{L}^1(\mathbf{m}) \otimes \mathcal{V}^1 \otimes I_d)^{T_{12456837}} \right]^{S_{23} \circ S_{\{45\} \{67\}}} \quad (C.7b)$$

$$\mathcal{L}^2(\mathbf{m}) = [(\mathcal{L}^1(\mathbf{m}) \otimes I_d)^{T_{1245736}}]^{S_{23} \circ S_{56}} \quad (C.7c)$$

Table C.1 lists the properties of all the tensors used in this appendix, starting with the tensors $\mathcal{L}$, $\mathcal{L}^1$, $\mathcal{L}^{11}$ and $\mathcal{L}^2$ just defined.

The symmetrization operations outside the square brackets in (C.7b-C.7c) are a matter of convention. They reflect the symmetries of the tensors with which the operators $\mathcal{L}$ are contracted.

The ‘content’ column in Table C.1 can be explained as follows. By design, $\mathbf{y}^*_{[0]} = \mathcal{L}(\mathbf{m}) \cdot h$ can be ‘unpacked’ (i.e., expressed in terms of $l$ and $m$) as $\mathbf{y}^*_{[0]} = Y_0(\mathbf{m}) \cdot l$, which is a function of $m$ contracted with $l$: the dependence on $m$ will be treated implicitly, and we express this by writing that the content of $\mathbf{y}^*_{[0]} = \mathcal{L}(\mathbf{m}) \cdot h$ is $l$, hence the symbol $l$ appearing in the ‘content’ column for the row $\mathcal{L}$. Similarly, $\nabla \mathbf{y}^*_{[0]} = (\mathcal{L}^1(\mathbf{m}) \cdot h) : \nabla h$ can be unpacked as

$$\begin{align}
\nabla \mathbf{y}^*_{[0]} &= \nabla(Y_0(\mathbf{m}(X)) \cdot l(X)) = (dY_0/d\mathbf{m})_{yl,m} l_l \nabla \mathbf{m}_{m,a} + (Y_0)_{yl} \nabla l_{l,a} \nonumber \\
&\quad (dY_0/d\mathbf{m}) (\mathbf{m}) : (l \otimes \nabla \mathbf{m}) + Y_0(\mathbf{m}) \cdot \nabla l \tag{C.8}
\end{align}$$

this is the sum of two terms, one being a function of $m$ contracted with $l \otimes \nabla m$, the other one being a function of $m$ contracted with $\nabla l$: this is conveyed by the content column in the table, which shows $(l \otimes \nabla m, \nabla l)$ for the row labelled $\mathcal{L}^1(\mathbf{m})$ used for reconstructing $\nabla \mathbf{y}^*_{[0]} = (\mathcal{L}^1(\mathbf{m}) \cdot h) : \nabla h$.

The point of the $h$ notation is to deal in a simple way with the multiplicity of terms appearing in the last column of Table 1. The remainder of the appendix will make use of this higher-level $h$ notation. On the other hand, we use Table 1 to keep track of the actual content of the various tensors.
Now, we proceed to represent the strain in terms of h and its successive gradients. Inserting (C.6) into the strain expression E in (2), identifying $E^{[0]}(m, l)$ using (B.10), and rearranging the other terms, we get

$$ 
\begin{aligned}
E = E^{[0]}(m, l) &+ (\mathcal{J}^1(m)(X) \cdot h(X)) : \nabla h(X) \\
&+ (\mathcal{J}^{11}(m)(X) \cdot h(X)) :: (\nabla h(X) \otimes \nabla h(X)) \\
&+ (\mathcal{J}^2(m)(X) \cdot h(X)) :: \nabla^2 h(X) \\
&+ E_y(m) \cdot (\mathbf{y}_{[1]}(X) + \mathbf{y}_{[2]}(X)) + E'_y(m) : \nabla \mathbf{y}_{[1]}(X) + O(\eta^3)
\end{aligned} 
\quad (C.9) $$

where the so-called structure coefficients are identified by

$$ \mathcal{J}^1(m) = (E'_l^{T_{132}}(m) \cdot \mathbf{V}^l \otimes \mathbf{V}^1)^{T_{1324}} + E'_y(m) : \mathcal{L}^1(m) \quad (C.10a) $$

$$ \mathcal{J}^{11}(m) = E''_y(m) : \mathcal{L}^{11}(m) \quad (C.10b) $$

$$ \mathcal{J}^2(m) = (E''_{l}^{T_{1423}}(m) \cdot \mathbf{V}^l \otimes \mathbf{V}^1)^{T_{1425}} + E''_y(m) : \mathcal{L}^2(m) \quad (C.10c) $$

As indicated in Table 1, $\mathcal{J}^{11}$ is $S_{\{23\}\{45\}}$-symmetric: this is a consequence of the fact that $\mathcal{L}^{11}$ is $S_{\{45\}\{67\}}$-symmetric. The $S_{34}$-symmetry of $\mathcal{J}^2$ can be justified by a similar argument. Grouping the terms order by order, we can rewrite the strain in (C.9) as

$$ E = E^{[0]}(m, l) + E^{[1]} + E^{[2]} + O(\eta^3), \quad (C.11) $$

where the contributions $E^{[1]} = O(\eta)$ and $E^{[2]} = O(\eta^2)$ are given, respectively, by

$$ E^{[1]} = (\mathcal{J}^1(m) \cdot h) : \nabla h + E_y(m) \cdot \mathbf{y}_{[1]} \quad (C.12a) $$

$$ 
\begin{aligned}
E^{[2]} = (\mathcal{J}^{11}(m) \cdot h) :: & (\nabla h \otimes \nabla h) \\
& + (\mathcal{J}^2(m) \cdot h) :: \nabla^2 h + E'_y(m) : \nabla \mathbf{y}_{[1]} + E_y(m) \cdot \mathbf{y}_{[2]}.
\end{aligned} 
\quad (C.12b) $$

Inserting (C.11-C.12) into the kinematic constraint $Q \cdot E = 0$, see (3), we obtain the expression of the constraint order by order in $\eta$ as

$$ Q \cdot E^{[i]} = 0 \quad \text{for } i = 0, 1, 2, \dots \quad (C.13) $$

The constraint $Q \cdot E^{[0]}(m, l) = 0$ has been enforced during the solution of the leading order, see (B.13b).

### C.3 Strain energy expansion in terms of corrective displacement

Inserting (C.11) into the energy (4), and using the quadratic expression (5) of the energy density, we obtain a Taylor expansion of the energy as

$$ \Phi = \int_{\Omega} \left( \frac{1}{2} E^{[0]}(m, l) \cdot \mathcal{K}(m) \cdot E^{[0]}(m, l) + (\mathcal{K}(m) \cdot E^{[0]}(m, l)) \cdot (E^{[1]} + E^{[2]}) \right) dX. \quad (C.14) $$

Grouping the terms in the integrand order by order, identifying the term of order $\eta^0$ as $\frac{1}{2}l \cdot K_0 \cdot l$ using (B.12), and the quantity $\mathcal{K}(m) \cdot E^{[0]}(m, l) = \mathcal{K}(m) \cdot F^{[0]}(m) \cdot l = (S_0(m) - Q^T \cdot G_0(m)) \cdot l$ by (B.10) and (B.14), we obtain the energy expansion as

$$ \Phi = \Phi^*_{[0]}[m, l] + \Phi_{[1]} + \Phi_{[2]} + \dots \quad (C.15) $$

where $\Phi_{[i]} = O(\eta^i)$ are the successive terms in the expansion, the leading term is the quantity $\Phi^*_{[0]}[m, l]$ identified earlier in (65) and the higher-order terms are given by

$$ 
\begin{aligned}
\Phi_{[1]} &= \int_{\Omega} ((S_0(m) - Q^T \cdot G_0(m)) \cdot l) \cdot E^{[1]} dX \\
\Phi_{[2]} &= \int_{\Omega} \left( \frac{1}{2} E^{[1]} \cdot \mathcal{K}(m) \cdot E^{[1]} + ((S_0(m) - Q^T \cdot G_0(m)) \cdot l) \cdot E^{[2]} \right) dX.
\end{aligned} 
\quad (C.16) $$
The $Q^T$ terms appearing in both bulk integrals can be removed, as $-(Q^T \cdot G_0(m) \cdot l) \cdot E^{[i]} = -(G_0(m) \cdot l) \cdot (Q \cdot E^{[i]}) = 0$ by (C.13). This yields

$$ \Phi_{[1]} = \int_{\Omega} (S_0(m) \cdot l) \cdot E^{[1]} dX \qquad (C.17a) $$

$$ \Phi_{[2]} = \int_{\Omega} \left( \frac{1}{2} E^{[1]} \cdot \mathcal{K}(m) \cdot E^{[1]} + (S_0(m) \cdot l) \cdot E^{[2]} \right) dX. \qquad (C.17b) $$

Inserting the expression of $E^{[1]}$ from (C.12a) in $\Phi_{[1]}$ and using the principle of virtual work at dominant order in (B.15a), we have

$$ \begin{aligned} \Phi_{[1]} &= \int_{\Omega} (S_0(m) \cdot l) \cdot (\mathcal{J}^1(m) \cdot h) : \nabla h + E_y(m) \cdot y_{[1]} \, dX \\ &= \int_{\Omega} (S_0(m) \cdot l) \cdot (\mathcal{J}^1(m) \cdot h) : \nabla h + ([E_y^T(m) \cdot S_0(m)] \cdot l) \cdot y_{[1]} \, dX \quad (C.18) \\ &= \int_{\Omega} (S_0(m) \cdot l) \cdot (\mathcal{J}^1(m) \cdot h) : \nabla h dX \end{aligned} $$

By a similar argument, the expression of $E^{[2]}$ in (C.12b) can be inserted in the bulk integral appearing in $\Phi_{[2]}$, which shows that the term $(S_0(m) \cdot l) \cdot E_y(m) \cdot y_{[2]}$ is zero. This yields

$$ \Phi_{[1]} = \int_{\Omega} (S_0(m) \cdot l) \cdot (\mathcal{J}^1(m) \cdot h) : \nabla h dX \qquad (C.19a) $$

$$ \Phi_{[2]} = \int_{\Omega} \left( \frac{1}{2} E^{[1]} \cdot \mathcal{K}(m) \cdot E^{[1]} + (S_0(m) \cdot l) \cdot \left( (\mathcal{J}^{11}(m) \cdot h) :: (\nabla h \otimes \nabla h) + (\mathcal{J}^2(m) \cdot h) :: \nabla^2 h + E'_{y}(m) : \nabla y_{[1]} \right) \right) dX. \quad (C.19b) $$

### C.4 Correction at order $\eta$

We can rewrite (C.19a) in terms of the packed macroscopic variable $h$ as

$$ \Phi_{[1]}^{\star} = \int_{\Omega} (\mathcal{A}(m) : h \otimes h) : \nabla h(X) dX, \qquad (C.20) $$

where

$$ \mathcal{A}(m) = [(\mathcal{J}^1(m))^{T_{4123}} \cdot S_0(m) \cdot \mathcal{V}^l]^{S_{34}}. \qquad (C.21) $$

In (C.20), we use the star notation $\Phi_{[1]}^{\star}$ to emphasize that the right-hand side no longer depends on the unknown corrector $y_{[1]}$ and only depends on the macroscopic fields $m$ and $h$, thanks to the elimination of the unknown corrector $y_{[1]}$ done earlier in (C.18).

### C.5 Extraction of $A_0$ and $K_1$

As announced in Table C.1, the content of $\mathcal{A}(m)$ is $l \otimes l \otimes \nabla m$ and $l \otimes \nabla l$, which means that the right-hand side of Equation (C.20) can be unpacked using (C.3) as

$$ \Phi_{[1]}^{\star} = \int_{\Omega} \left( A_0(m) : (l \otimes \nabla l) + (K_1(m) : \nabla m) : \frac{l \otimes l}{2} \right) dX, \qquad (C.22) $$

where the tensors $A_0(m)$ and $K_1(m)$ are extracted from $\mathcal{A}(m)$ as

$$ \begin{aligned} A_0(m) &= 2 \left[ \mathcal{A}^{T_{2134}}(m) : \left( (\mathcal{V}^l)^T \otimes (\mathcal{V}^l)^T \otimes \mathcal{V}^1 \right)^{T_{14253}} \right]^{T_{321}} \\ K_1(m) &= 2 \left[ \mathcal{A}^{T_{2134}}(m) : \left( (\mathcal{V}^m)^T \otimes (\mathcal{V}^l)^T \otimes (\mathcal{V}^l)^T \right)^{T_{142536}} \right]^{T_{4312}} \end{aligned} \qquad (C.23) $$

No other term can be present in the right-hand side of (C.22): a term such as $D_0(m) :: (l \otimes l \otimes \nabla l)$, for instance, would be inconsistent with the fact that the energy is homogeneous with degree 2 in $l$ and $\nabla l$, see (2) and (5).

The properties of the tensors $A_0(m)$ and $K_1(m)$ are listed in Table 3.

The expression of the first correction $\Phi_{[1]}^{\star}$ to the energy in (C.22) was announced in (18).
### C.6 Energy correction at order $\eta^2$

We now turn attention to the correction $\Phi_{[2]}$ in (C.19b). Inserting the expression of $E^{[1]}$ in (C.12a), we get after rearranging the terms

$$ \Phi_{[2]} = \int_{\Omega} \left( \begin{aligned} & (\hat{\mathcal{B}}^{(0)}(\mathbf{m}) : (\mathbf{h} \otimes \mathbf{h})) :: \frac{\nabla \mathbf{h} \otimes \nabla \mathbf{h}}{2} \\ & + (\hat{\mathcal{B}}^{(1)}(\mathbf{m}) \cdot \mathbf{h}) :: (\nabla \mathbf{h} \otimes \mathbf{y}_{[1]}) + \mathcal{W}_{yy}(\mathbf{m}) : \frac{\mathbf{y}_{[1]} \otimes \mathbf{y}_{[1]}}{2} \end{aligned} \right) dX + \int_{\Omega} \left( (\mathcal{C}^{(0)}(\mathbf{m}) : (\mathbf{h} \otimes \mathbf{h})) :: \nabla^2 \mathbf{h} + (\mathcal{C}^{(1)}(\mathbf{m}) \cdot \mathbf{h}) : \nabla \mathbf{y}_{[1]} \right) dX, \quad (C.24) $$

where $\mathcal{W}_{yy}(\mathbf{m})$ is the tensor introduced in the analysis of the leading order, see (B.3c), and

$$ \begin{aligned} \hat{\mathcal{B}}^{(0)}(\mathbf{m}) &= ((\mathcal{J}^1(\mathbf{m}))^{T_{4123}} \cdot \mathcal{K}(\mathbf{m}) \cdot \mathcal{J}^1(\mathbf{m}))^{T_{125346}} + 2 ((\mathcal{J}^{11}(\mathbf{m}))^{T_{612345}} \cdot \mathcal{S}_0(\mathbf{m}) \cdot \mathcal{V}^l) \\ \hat{\mathcal{B}}^{(1)}(\mathbf{m}) &= ((\mathcal{J}^1(\mathbf{m}))^{T_{4123}} \cdot \mathcal{K}(\mathbf{m}) \cdot E_y(\mathbf{m}))^{T_{1243}} \quad (C.25) \\ \mathcal{C}^{(0)}(\mathbf{m}) &= [(\mathcal{J}^2(\mathbf{m}))^{T_{51234}} \cdot \mathcal{S}_0(\mathbf{m}) \cdot \mathcal{V}^l]^{S_{45}} \\ \mathcal{C}^{(1)}(\mathbf{m}) &= (E'_y(\mathbf{m}))^{T_{312}} \cdot \mathcal{S}_0(\mathbf{m}) \cdot \mathcal{V}^l. \end{aligned} $$

As indicated by the 'delayed' keyword in Table C.1, we do not yet enforce the natural symmetries of $\hat{\mathcal{B}}^{(0)}$, which reflect the symmetries of the tensor $\frac{1}{2}(\nabla\mathbf{h} \otimes \nabla\mathbf{h}) \otimes (\mathbf{h} \otimes \mathbf{h})$ with which it gets contracted: they will be enforced later on the children of $\hat{\mathcal{B}}^{(0)}$.

In the code, we implemented an extension of (C.25) that covers the rank-deficient case as well, see Equation (E.24) in Appendix E.

Note that $\Phi_{[2]}$ no longer depends on $\mathbf{y}_{[2]}$ thanks to the work done in Appendix C.3. It still depends on $\mathbf{y}_{[1]}(\mathbf{X})$ and its gradient, however. We proceed to remove the dependence on the gradient $\nabla\mathbf{y}_{[1]}(\mathbf{X})$ by integrating by parts.

### C.7 Integration by parts

The $\mathcal{C}^{(i)}$ terms appearing in (C.24) can be integrated by parts as

$$ \begin{aligned} \int_{\Omega} (\mathcal{C}^{(0)}(\mathbf{m}(\mathbf{X})) : (\mathbf{h}(\mathbf{X}) \otimes \mathbf{h}(\mathbf{X}))) :: & \nabla^2 \mathbf{h}(\mathbf{X}) dX = \\ \oint_{\partial\Omega} (\mathcal{C}^{(0)}(\mathbf{m}(\mathbf{X})) : (\mathbf{h}(\mathbf{X}) \otimes \mathbf{h}(\mathbf{X}))) :: & (\nabla \mathbf{h}(\mathbf{X}) \otimes \mathbf{n}(\mathbf{X})) d\mathbf{a} \\ + \int_{\Omega} (\Delta \mathcal{B}^{(0)}(\mathbf{m}(\mathbf{X})) : (\mathbf{h}(\mathbf{X}) \otimes \mathbf{h}(\mathbf{X}))) :: & \frac{\nabla \mathbf{h}(\mathbf{X}) \otimes \nabla \mathbf{h}(\mathbf{X})}{2} dX \end{aligned} \quad (C.26) $$

and

$$ \begin{aligned} \int_{\Omega} (\mathcal{C}^{(1)}(\mathbf{m}(\mathbf{X})) \cdot \mathbf{h}(\mathbf{X})) : \nabla \mathbf{y}_{[1]}(\mathbf{X}) dX = \\ \oint_{\partial\Omega} (\mathcal{C}^{(1)}(\mathbf{m}(\mathbf{X})) \cdot \mathbf{h}(\mathbf{X})) : (\mathbf{y}_{[1]}(\mathbf{X}) \otimes \mathbf{n}(\mathbf{X})) d\mathbf{a} \\ + \int_{\Omega} (\Delta \mathcal{B}^{(1)}(\mathbf{m}(\mathbf{X})) \cdot \mathbf{h}(\mathbf{X})) : (\nabla \mathbf{h}(\mathbf{X}) \otimes \mathbf{y}_{[1]}(\mathbf{X})) dX \end{aligned} \quad (C.27) $$

where

$$ \begin{aligned} \Delta \mathcal{B}^{(0)}(\mathbf{m}) &= -2 \left( \frac{d\mathcal{C}^{(0)}}{dm}(\mathbf{m}) \cdot \mathcal{V}^m \right)^{T_{342561}} - 4 (\mathcal{C}^{(0)}(\mathbf{m}))^{T_{12435}} \otimes \mathcal{V}^1 \\ \Delta \mathcal{B}^{(1)}(\mathbf{m}) &= - \left( \frac{d\mathcal{C}^{(1)}}{dm}(\mathbf{m}) \cdot \mathcal{V}^m \right)^{T_{3241}} - (\mathcal{C}^{(1)}(\mathbf{m}))^{T_{321}} \otimes \mathcal{V}^1, \end{aligned} \quad (C.28) $$
Note that the operator acting on $(\nabla h \otimes \nabla h)$ in the first equation above has been symmetrized with respect to an exchange of the $\nabla h$'s.

Inserting (C.26-C.27) into (C.24), we have

$$ \Phi_{[2]} = \oint_{\partial\Omega} \left( (C^{(0)}(m) : (h \otimes h)) : (\nabla h \otimes n) + (C^{(1)}(m) \cdot h) : (y_{[1]} \otimes n) \right) da \\ + \int_{\Omega} \left( \begin{array}{l} (\mathcal{B}^{(0)}(m) : (h \otimes h)) : \frac{\nabla h \otimes \nabla h}{2} \\ + (\mathcal{B}^{(1)}(m) \cdot h) : (\nabla h \otimes y_{[1]}) + W_{yy}(m) : \frac{y_{[1]} \otimes y_{[1]}}{2} \end{array} \right) dX \quad (C.29) $$

where

$$ \mathcal{B}^{(0)}(m) = \hat{\mathcal{B}}^{(0)}(m) + \Delta\mathcal{B}^{(0)}(m) \quad (C.30) $$

$$ \mathcal{B}^{(1)}(m) = \hat{\mathcal{B}}^{(1)}(m) + \Delta\mathcal{B}^{(1)}(m) $$

Equation (C.29) has been announced in (69-71).

## C.8 Optimal corrective displacement

The integrand of $\Phi_{[2]}$ in (C.29) depends on $y_{[1]}$ but not on its gradient, thanks to the integration by parts, Appendix C.7. It does not depend on $y_{[2]}$ either.

The variational problem (72) for $y_{[1]}$ is therefore local and we proceed to solve it. The functional $\Phi_{[2]}$ is given in (C.29) as the sum of a bulk integral and a boundary integral. The variational problem (72) therefore yields two sets of local conditions:

* at any point $X \in \partial\Omega$ on the boundary, the increment $(C^{(1)}(m) \cdot h) : (\delta y \otimes n)$ coming from the boundary integral should vanish for any perturbation $\delta y$ satisfying the incremental constraint $Q \cdot E_y(m) \cdot \delta y = 0$; using a Lagrange multiplier $\tilde{g}$, we must solve $(C^{(1)}(m) \cdot h) : (\delta \tilde{y} \otimes n) + \tilde{g} \cdot Q \cdot E_y(m) \cdot \delta \tilde{y} = 0$ for any $\delta \tilde{y} \in \mathbb{R}^{n_y}$. Eliminating the virtual quantity $\delta \tilde{y}$, we can rewrite this as $C^{(1)}(m) : (n \otimes h) + (Q \cdot E_y(m))^T \cdot \tilde{g} = 0$. A solution $\tilde{g}$ exists if and only if the vector $C^{(1)}(m) : (n \otimes h)$ is contained in the image of the operator $(Q \cdot E_y(m))^T$. Inserting the expression of $h$ in (C.3), we therefore obtain the stationarity condition in the form

$$ C^{(1)}(m) : (n \otimes ((V^l)^T \cdot l + (V^m)^T \cdot m + V^1)) \in \operatorname{Im}(Q \cdot E_y(m))^T. \quad (C.31) $$

* at any point $X \in \Omega^\circ$ in the interior of the domain $\Omega$, the problem for $y_{[1]}^*(X)$ is a quadratic optimization problem with a linear constraint. Using a Lagrange multiplier $g_{[1]}^*(X)$, its solution $(\tilde{y}, \tilde{g}) = (y_{[1]}^*(X), g_{[1]}^*(X))$ must satisfy

$$ \left\{ \begin{array}{l} Q \cdot (E_y(m(X)) \cdot \tilde{y} + (J^1(m(X)) \cdot h(X)) : \nabla h(X)) = 0 \\ W_{yy}(m(X)) : (\tilde{y} \otimes \delta \tilde{y}) + (\mathcal{B}^{(1)}(m) \cdot h) : (\nabla h(X) \otimes \delta \tilde{y}) \\ \qquad + \tilde{g} \cdot Q \cdot E_y(m(X)) \cdot \delta \tilde{y} = 0, \quad (\forall \delta \tilde{y} \in \mathbb{R}^{n_y}). \end{array} \right. \quad (C.32) $$

We take note of the stationarity condition (C.31) on the boundary, which we will address in future work, and proceed to solve the problem (C.32) in the interior. The latter can be rewritten in matrix form as

$$ P(m(X)) \cdot \begin{pmatrix} y_{[1]}^*(X) \\ g_{[1]}^*(X) \end{pmatrix} + \begin{pmatrix} (\mathcal{B}^{(1)}(m(X)))^{T_{2314}} \\ Q \cdot J^1(m(X)) \end{pmatrix} : (\nabla h(X) \otimes h(X)) = 0. \quad (C.33) $$

The matrix $P$ in the left-hand side above is identical to the one which appeared in the analysis of the leading order, see (B.4).

The solution $(y_{[1]}^*(X), g_{[1]}^*(X))$ is obtained by inverting $P$ (the case where $P$ is non-invertible is treated in Appendix E). This yields the correction $y_{[1]}^*(X)$ to the microscopic degrees of freedom in form

$$ y_{[1]}^*(X) = (\mathcal{Y}'(m(X)) \cdot h(X)) : \nabla h(X), \quad (C.34) $$
with the localization tensor given by

$$ \mathbf{Y}'(\mathbf{m}) = \left( \mathbf{I}_{n_y} \quad \mathbf{0}_{n_y \times n_c} \right) \cdot \mathbf{R}'(\mathbf{m}) \qquad (\text{C.35}) $$

where

$$ \mathbf{R}'(\mathbf{m}) = -\mathbf{P}^{-1}(\mathbf{m}) \cdot \begin{pmatrix} (\mathcal{B}^{(1)}(\mathbf{m}))^{T_{2314}} \\ \mathcal{Q} \cdot \mathcal{J}^{1}(\mathbf{m}) \end{pmatrix}. \qquad (\text{C.36}) $$

Compared to the leading-order in (B.7), only the second factor in the right-hand side has changed. To handle the rank-deficient case, we implement an extension of (C.36) that uses the pseudo-inverse rather than the inverse, see Equation (E.27) in Appendix E.

The Lagrange multiplier is given as a byproduct as

$$ \mathbf{g}^*_{[1]}(\mathbf{X}) = (\mathcal{G}'(\mathbf{m}(\mathbf{X})) \cdot \mathbf{h}(\mathbf{X})) : \nabla \mathbf{h}(\mathbf{X}) \text{ where } \mathcal{G}'(\mathbf{m}) = (\mathbf{0}_{n_c \times n_y} \quad \mathbf{I}_{n_c}) \cdot \mathbf{R}'(\mathbf{m}). \qquad (\text{C.37}) $$

### C.9 Relaxed energy correction at order $\eta^2$

Inserting the solution $\mathbf{y}_{[1]}(\mathbf{X}) = \mathbf{y}^*_{[1]}(\mathbf{X})$ given in (C.34), we can rewrite the boundary terms in (C.29) as

$$ (\Phi_{[2]}^{\text{bt}})^* [\mathbf{h}] = \oint_{\partial\Omega} (\mathcal{C}(\mathbf{m}) : (\mathbf{h} \otimes \mathbf{h})) \circ (\nabla \mathbf{h} \otimes \mathbf{n}) d\mathbf{a}. \qquad (\text{C.38}) $$

where

$$ \mathcal{C}(\mathbf{m}) = \mathcal{C}^{(0)}(\mathbf{m}) + [((\mathcal{C}^{(1)}(\mathbf{m}))^{T_{312}} \cdot \mathbf{Y}'(\mathbf{m}))^{T_{34125}}]^{S_{45}}. \qquad (\text{C.39}) $$

The optimality condition (C.33) can be split in two sets of equations, namely

$$ (\mathcal{B}^{(1)}(\mathbf{m}))^{T_{2314}} \circ (\nabla \mathbf{h} \otimes \mathbf{h}) + \mathcal{W}_{yy}(\mathbf{m}) \cdot \mathbf{y}^*_{[1]}(\mathbf{X}) + \mathbf{g}^*_{[1]}(\mathbf{X}) \cdot \mathcal{Q} \cdot \mathcal{E}_y(\mathbf{m}) = \mathbf{0} \qquad (\text{C.40}) $$

and

$$ \mathbf{Q} \cdot \mathcal{E}_y(\mathbf{m}) \cdot \mathbf{y}^*_{[1]}(\mathbf{X}) + \mathbf{Q} \cdot \mathcal{J}^{1}(\mathbf{m}) \circ (\nabla \mathbf{h} \otimes \mathbf{h}) = \mathbf{0}. \qquad (\text{C.41}) $$

Taking the dot products of the first equation by $\mathbf{y}^*_{[1]}(\mathbf{X})$ and of the second equation by $\mathbf{g}^*_{[1]}(\mathbf{X})$ and subtracting, we obtain

$$ (\mathcal{B}^{(1)}(\mathbf{m}) \cdot \mathbf{h}) \circ (\nabla \mathbf{h}(\mathbf{X}) \otimes \mathbf{y}^*_{[1]}(\mathbf{X})) + \mathcal{W}_{yy}(\mathbf{m}) \circ (\mathbf{y}^*_{[1]}(\mathbf{X}) \otimes \mathbf{y}^*_{[1]}(\mathbf{X})) - \mathbf{g}^*_{[1]}(\mathbf{X}) \cdot \mathcal{Q} \circ ((\mathcal{J}^{1}(\mathbf{m}) \cdot \mathbf{h}) : \nabla \mathbf{h}) = \mathbf{0}. \qquad (\text{C.42}) $$

This identity can be used to eliminate the $\mathcal{B}^{(1)}$ term appearing in the bulk integral $\Phi_{[2]}^{\text{it}}$ over $\Omega$ in (C.29) as

$$ \Phi_{[2]}^{\text{it}} = \int_{\Omega} \left( \begin{aligned} & (\mathcal{B}^{(0)}(\mathbf{m}) : (\mathbf{h} \otimes \mathbf{h})) \circ \frac{\nabla \mathbf{h} \otimes \nabla \mathbf{h}}{2} \\ & + \mathbf{g}^*_{[1]} \cdot \mathbf{Q} \circ ((\mathcal{J}^{1}(\mathbf{m}) \cdot \mathbf{h}) : \nabla \mathbf{h}) - \mathcal{W}_{yy}(\mathbf{m}) : \frac{\mathbf{y}_{[1]} \otimes \mathbf{y}_{[1]}}{2} \end{aligned} \right) d\mathbf{X}. \qquad (\text{C.43}) $$

Inserting the expressions of $\mathbf{y}^*_{[1]}$ and $\mathbf{g}^*_{[1]}$ obtained earlier in (C.34) and (C.37), we can rewrite this as

$$ (\Phi_{[2]}^{\text{it}})^* [\mathbf{h}] = \int_{\Omega} (\mathcal{B}(\mathbf{m}) : (\mathbf{h} \otimes \mathbf{h})) \circ \frac{\nabla \mathbf{h} \otimes \nabla \mathbf{h}}{2} d\mathbf{X}, \qquad (\text{C.44}) $$
where

$$ \mathcal{B}(\mathbf{m}) = \left[ \mathcal{B}^{(0)}(\mathbf{m}) + \begin{pmatrix} 2\mathcal{G}'^{T_{4123}}(\mathbf{m}) \cdot \mathbf{Q} \cdot \mathcal{J}^1(\mathbf{m}) \\ -\mathbf{y}'^{T_{4123}}(\mathbf{m}) \cdot \mathcal{W}_{yy}(\mathbf{m}) \cdot \mathbf{y}'(\mathbf{m}) \end{pmatrix} \right]^{S_{\{12\}}(34) \circ S_{56}} \quad (C.45) $$

By Equations (C.38) and (C.44), the energy contribution $\Phi_{[2]}^*[h] = (\Phi_{[2]}^{\text{it}})^*[h] + (\Phi_{[2]}^{\text{bt}})^*[h]$ is given by

$$ \Phi_{[2]}^*[h] = \int_{\Omega} (\mathcal{B}(\mathbf{m}):(\mathbf{h} \otimes \mathbf{h})) \dfrac{\nabla \mathbf{h} \otimes \nabla \mathbf{h}}{2} dX + \oint_{\partial \Omega} (\mathcal{C}(\mathbf{m}):(\mathbf{h} \otimes \mathbf{h})) \cdot (\nabla \mathbf{h} \otimes \mathbf{n}) da. \quad (C.46) $$

## C.10 Final extraction

The last step in the homogenization procedure is to unpack the tensors $\mathbf{y}', \mathcal{B}$ and $\mathcal{C}$ appearing in (C.34) and (C.46), following a similar procedure as earlier in Appendix C.5. This allows to remove any reference to the quantity $\mathbf{h}$ used internally in this Appendix in favor of $\mathbf{l}$ and $\mathbf{m}$.

With the help of the ‘content tracking’ done in Table C.1, we obtain the unpacked form of these tensors as

$$ (\mathbf{y}'(\mathbf{m}) \cdot \mathbf{h}) : \nabla \mathbf{h} = (Y_1(\mathbf{m}) : \nabla \mathbf{m}) \cdot \mathbf{l} + Y'_0(\mathbf{m}) : \nabla \mathbf{l} \quad (C.47a) $$

$$ \begin{aligned} (\mathcal{B}(\mathbf{m}): (\mathbf{h} \otimes \mathbf{h})) :: & \frac{\nabla \mathbf{h} \otimes \nabla \mathbf{h}}{2} = (K_2(\mathbf{m}) :: (\nabla \mathbf{m} \otimes \nabla \mathbf{m})) : \frac{\mathbf{l} \otimes \mathbf{l}}{2} \\ & + (A_1(\mathbf{m}) : \nabla \mathbf{m}) :: (\mathbf{l} \otimes \nabla \mathbf{l}) + B_0(\mathbf{m}) :: \frac{\nabla \mathbf{l} \otimes \nabla \mathbf{l}}{2} \end{aligned} \quad (C.47b) $$

$$ (\mathcal{C}(\mathbf{m}): (\mathbf{h} \otimes \mathbf{h})) \cdot (\nabla \mathbf{h} \otimes \mathbf{n}) = (k_1(\mathbf{m}) : \nabla \mathbf{m}) \cdot \left( \frac{\mathbf{l} \otimes \mathbf{l}}{2} \otimes \mathbf{n} \right) + a_0(\mathbf{m}) \cdot (\mathbf{l} \otimes \nabla \mathbf{l} \otimes \mathbf{n}) \quad (C.47c) $$

where the sub-tensors $Y_1(\mathbf{m})$, $Y'_0(\mathbf{m})$, $K_2(\mathbf{m})$, $A_1(\mathbf{m})$, $B_0(\mathbf{m})$, $k_1(\mathbf{m})$ and $a_0(\mathbf{m})$ are identified as

$$ Y_1(\mathbf{m}) = \left( \mathbf{y}'^{T_{1324}} : \left( (\mathbf{V}^m)^T \otimes (\mathbf{V}^l)^T \right)^{T_{1324}} \right)^{T_{1432}} \quad (C.48a) $$

$$ Y'_0(\mathbf{m}) = \left( \mathbf{y}'^{T_{1324}} : \left( (\mathbf{V}^l)^T \otimes \mathbf{V}^1 \right)^{T_{132}} \right)^{T_{132}} \quad (C.48b) $$

$$ K_2(\mathbf{m}) = \left( \mathcal{B}^{T_{314256}}(\mathbf{m}) :: \left( (\mathbf{V}^m)^T \otimes (\mathbf{V}^m)^T \otimes (\mathbf{V}^l)^T \otimes (\mathbf{V}^l)^T \right)^{T_{15263748}} \right)^{T_{463512}} \quad (C.48c) $$

$$ A_1(\mathbf{m}) = 2 \left( \mathcal{B}^{T_{314256}}(\mathbf{m}) :: \left( (\mathbf{V}^l)^T \otimes (\mathbf{V}^m)^T \otimes (\mathbf{V}^l)^T \otimes \mathbf{V}^1 \right)^{T_{1526374}} \right)^{T_{35241}} \quad (C.48d) $$

$$ B_0(\mathbf{m}) = \left( \mathcal{B}^{T_{314256}}(\mathbf{m}) :: \left( (\mathbf{V}^l)^T \otimes (\mathbf{V}^l)^T \otimes \mathbf{V}^1 \otimes \mathbf{V}^1 \right)^{T_{152634}} \right)^{T_{2413}} \quad (C.48e) $$

$$ k_1(\mathbf{m}) = 2 \left( \mathcal{C}^{T_{31245}}(\mathbf{m}) :: \left( (\mathbf{V}^m)^T \otimes (\mathbf{V}^l)^T \otimes (\mathbf{V}^l)^T \right)^{T_{142536}} \right)^{T_{53412}} \quad (C.48f) $$

$$ a_0(\mathbf{m}) = 2 \left( \mathcal{C}^{T_{31245}}(\mathbf{m}) :: \left( (\mathbf{V}^l)^T \otimes (\mathbf{V}^l)^T \otimes \mathbf{V}^1 \right)^{T_{14253}} \right)^{T_{3421}} \quad (C.48g) $$

The expression of $\mathbf{y}_{[1]}^*$ announced in (73) then follows from (C.34) and (C.47a), whereas the expression of $\Phi_{[2]}^*$ announced in (19) follows from (C.46) and (C.47b-C.47c).

## D Special case of homogeneous properties

The special case of homogeneous properties (applicable to a perfectly periodic elastic truss for instance) is considered here. In this special case, the parameter $\mathbf{m}(\mathbf{X})$ goes away ($n_m = 0$).
The analysis of leading order is unchanged: it delivers $Y_0$ and $S_0$, which are no longer functions of $m$, as illustrated in Table 5 for the non-deficient case. There are many simplifications at the next orders and the corresponding, specialized formulas are provided in Table D.2. Note that the tensor dimensions are changed compared to the general case, see also the ‘usage’ column in the table. As a consequence, the indices used in the transpose operations are affected. Although it is straightforward in principle, the specialization of the general formulas to this special case is cumbersome—to a point that we found it easier to re-derive them from scratch. Their consistency with the general formulas is checked in a dedicated Mathematica notebook<fnref n="2" />.

As can be expected, all tensors that get contracted with gradients of $m$ are zero, $K_1 = 0$, $K_2 = 0$, $A_1 = 0$, $k_1 = 0$, $Y_1 = 0$.

# E Extension to a rank-deficient matrix

## E.1 Special form of null vectors

Assuming that they exist, let us first characterize the null vectors of the symmetric matrix $P(m)$ introduced in (B.5), entering in both the leading order problem (B.4) and in the determination of the corrective displacement (C.33).

For any $z = (z_y, z_c) \in \mathbb{R}^{n_y+n_c}$ such that $P(m) \cdot z = 0$, we have

$$ \begin{aligned} \mathcal{W}_{yy}(m) \cdot z_y + z_c \cdot \mathbf{Q} \cdot E_y(m) &= 0 \\ \mathbf{Q} \cdot E_y(m) \cdot z_y &= 0 \end{aligned} \qquad (E.1) $$

Multiplying the first equation by $z_y$, and using the second equation, we get $z_y \cdot \mathcal{W}_{yy}(m) \cdot z_y = 0$. We observe that the assumption (12) (positive-definiteness of the energy on the subspace of admissible microscopic degrees of freedom) can be rewritten as: $\mathbf{Q} \cdot E_y(m) \cdot \mathbf{y} = 0$ and $\mathbf{y} \cdot \mathcal{W}_{yy}(m) \cdot \mathbf{y} > 0$ implies $\mathbf{y} = 0$. Therefore, we have $z_y = 0$, which then yields $(\mathbf{Q} \cdot E_y(m))^T \cdot z_c = 0$, i.e., the $z_c$ block is a null vector of $(\mathbf{Q} \cdot E_y(m))^T$. We have just shown

$$ P(m) \cdot z = 0 \quad \Rightarrow \quad z = (0_{n_y}, z_c) \quad \text{with} \quad (\mathbf{Q} \cdot E_y(m))^T \cdot z_c = 0. \qquad (E.2) $$

The only way that the matrix $P(m)$ can be singular is because of the $\mathbf{Q} \cdot E_y(m)$ block.

With $n_d$ denoting the rank deficiency of the matrix $P$ or $(\mathbf{Q} \cdot E_y(m))^T$ (both are the same by the argument above), we denote as $\mathbf{N}(m) \in \mathbb{T}^{(n_d, n_c)}$ a list of null vectors of $(\mathbf{Q} \cdot E_y(m))^T$, arranged in rows. Equation (E.2) then shows that the null vectors of $P(m)$ are the rows of

$$ N_P(m) = \mathbf{N}(m) \cdot \begin{pmatrix} 0_{n_c \times n_y} & I_{n_c} \end{pmatrix} \in \mathbb{T}^{(n_d, n_y+n_c)}. \qquad (E.3) $$

## E.2 Solutions of the linear equation

We consider the linear equation for a vector $X \in \mathbb{R}^{n_y+n_c}$,

$$ P(m) \cdot X = Y. \qquad (E.4) $$

Multiplying by any null vector $z$ of $P$ and using the symmetry $P^T = P$, we obtain $z \cdot Y = 0$. Repeating this argument with all the null vectors that have been arranged into $N_P(m)$, we obtain $n_d$ solvability conditions

$$ N_P(m) \cdot Y = 0. \qquad (E.5) $$

When (E.5) is satisfied, the solutions $X$ of (E.4) can be expressed with the help of the Moore-Penrose inverse $P^\dagger(m)$ of $P(m)$ as

$$ X = P^\dagger(m) \cdot Y + N_P^T(m) \cdot \hat{l}. \qquad (E.6) $$

<footnote n="2">² See shoal-library-v1.0/discrete_engine/tests/verifySpecialFormulasHomogeneousCase.nb</footnote>
<table><caption>**Table D.2** Higher-order homogenization in the special case of homogeneous properties. A complete implementation of the method in this special case is possible based on Table 5 (leading order, ignoring any dependence on *m*) and on the definitions appearing in the first column of the table above. The quantities appearing in the grey rows are the main results (localization tensors for corrective displacement $Y'_0$, Lagrange multipliers $G'_0$, bulk energy contributions $A_0$ and $B_0$ and boundary contribution $A_0$).</caption><thead><tr><th>definition</th><th>dimension</th><th>sym.</th><th>usage</th></tr></thead><tbody><tr><td>L<sup>1</sup> = (Y<sub>0</sub> ⊗ I<sub>d</sub>)<sup>T<sub>324</sub></sup></td><td>T<sub>0</sub>(n<sub>y</sub>,d,n<sub>1</sub>,d)</td><td></td><td>∇y<sub>l</sub><sup>*</sup><sub>0</sub> = L<sup>1</sup> : ∇l</td></tr><tr><td>L<sup>2</sup> = [(Y<sub>0</sub> ⊗ I<sub>d</sub> ⊗ I<sub>d</sub>)<sup>T<sub>4236</sub></sup>]<sup>S<sub>23</sub>○S<sub>56</sub></sup></td><td>T<sub>0</sub>(n<sub>y</sub>,d,d,n<sub>1</sub>,d,d)</td><td>S<sub>23</sub>, S<sub>56</sub></td><td>∇<sup>2</sup>y<sub>l</sub><sup>*</sup><sub>0</sub> = L<sup>2</sup> : ∇<sup>2</sup>l</td></tr><tr><td>S<sup>1</sup> = E<sub>l</sub>' + E'<sub>y</sub> : L<sup>1</sup></td><td>T<sub>0</sub>(n<sub>E</sub>,n<sub>1</sub>,d)</td><td></td><td>E<sup>[1]</sup> = S<sup>1</sup> : ∇l + ...</td></tr><tr><td>S<sup>2</sup> = E<sub>l</sub>'' + E'<sub>y</sub> : L<sup>2</sup></td><td>T<sub>0</sub>(n<sub>E</sub>,n<sub>1</sub>,d,d)</td><td>S<sub>34</sub></td><td>E<sup>[2]</sup> = S<sup>2</sup> : ∇<sup>2</sup>l + ...</td></tr><tr><td>A<sub>0</sub> = [(S<sup>1</sup>)<sup>T<sub>312</sub></sup> · S<sub>0</sub>]<sup>T<sub>231</sub></sup></td><td>T<sub>0</sub>(n<sub>1</sub>,n<sub>1</sub>,d)</td><td></td><td>Φ<sup>[1]</sup> = ∫<sub>0</sub>A<sub>0</sub> : (l ⊗ ∇l)dX</td></tr><tr><td>B̂<sup>(0)</sup> = (S<sup>1</sup>)<sup>T<sub>312</sub></sup> · K̂ · S<sup>1</sup></td><td>T<sub>0</sub>(n<sub>1</sub>,d,n<sub>1</sub>,d)</td><td>delayed</td><td>Φ<sup>[2]</sup> = ∫<sub>0</sub>B̂<sup>(0)</sup> : ∇l/2 ∇l dX + ...</td></tr><tr><td>B̂<sup>(1)</sup> = (S<sup>1</sup>)<sup>T<sub>312</sub></sup> · K̂ · E<sub>y</sub></td><td>T<sub>0</sub>(n<sub>1</sub>,d,n<sub>y</sub>)</td><td>–</td><td>Φ<sup>[2]</sup> = ∫<sub>0</sub>B̂<sup>(0)</sup> : (∇l ⊗ y<sub>l</sub>)dX + ...</td></tr><tr><td>C<sup>(0)</sup> = (S<sup>2</sup>)<sup>T<sub>423</sub></sup> · S<sub>0</sub></td><td>T<sub>0</sub>(n<sub>1</sub>,d,d,n<sub>1</sub>)</td><td>S<sub>23</sub></td><td>Φ<sup>[2]</sup> = ∫<sub>0</sub>(C<sup>(0)</sup> · l) : ∇<sup>2</sup>l dX + ...</td></tr><tr><td>C<sup>(1)</sup> = (E'<sub>y</sub>)<sup>T<sub>312</sub></sup> · S<sub>0</sub></td><td>T<sub>0</sub>(n<sub>y</sub>,d,n<sub>1</sub>)</td><td>–</td><td>Φ<sup>[2]</sup> = ∫<sub>0</sub>(C<sup>(1)</sup> · l) : ∇y<sub>1</sub> dX + ...</td></tr><tr><td>B̂<sup>(0)</sup> = B̂<sup>(0)</sup> - 2(C<sup>(0)</sup>)<sup>T<sub>423</sub></sup></td><td>T<sub>0</sub>(n<sub>1</sub>,d,n<sub>1</sub>,d)</td><td>delayed</td><td>same as B̂<sup>(0)</sup></td></tr><tr><td>B̂<sup>(1)</sup> = B̂<sup>(1)</sup> - (C<sup>(1)</sup>)<sup>T<sub>321</sub></sup></td><td>T<sub>0</sub>(n<sub>1</sub>,d,n<sub>y</sub>)</td><td>–</td><td>same as B̂<sup>(1)</sup></td></tr><tr><td>R̃ = -P<sup>†</sup> · (B̂<sup>(1)</sup>)<sup>T<sub>231</sub></sup> Q̂ · F̂<sup>1</sup></td><td>T<sub>0</sub>(n<sub>y</sub>+n<sub>c</sub>,n<sub>1</sub>,d)</td><td>–</td><td>(y<sub>l</sub><sup>*</sup><sub>1</sub> g<sub>l</sub><sup>*</sup><sub>1</sub>) = R̃' : ∇l</td></tr><tr><td>Y<sub>0</sub> = (I<sub>n<sub>y</sub></sub> 0<sub>n<sub>y</sub>×n<sub>c</sub></sub>) · R̃'</td><td>T<sub>0</sub>(n<sub>y</sub>,n<sub>1</sub>,d)</td><td>–</td><td>y<sub>l</sub><sup>*</sup><sub>1</sub> = Y'<sub>0</sub> : ∇l</td></tr><tr><td>G<sub>0</sub> = (0<sub>n<sub>c</sub>×n<sub>y</sub></sub> I<sub>n<sub>c</sub></sub>) · R̃'</td><td>T<sub>0</sub>(n<sub>c</sub>,n<sub>1</sub>,d)</td><td>–</td><td>g<sub>l</sub><sup>*</sup><sub>1</sub> = G'<sub>0</sub> : ∇l</td></tr><tr><td>B<sub>0</sub> = [B̂<sup>(0)</sup>+2G<sub>0</sub><sup>T<sub>312</sub></sup> · Q̂ · S<sup>1</sup>]<sup>S<sub>(12)</sub>|<sup>(34)</sup></sup> - Y'<sub>0</sub><sup>T<sub>312</sub></sup> · W<sub>y</sub>y · Y'<sub>0</sub></td><td>T<sub>0</sub>(n<sub>1</sub>,d,n<sub>1</sub>,d)</td><td>S<sub>{12}</sub>{34}</td><td>Φ<sup>*</sup><sub>l</sub><sup>[2]</sup> = ∫<sub>0</sub>B<sub>0</sub> : ∇l/2 ∇l dX + ...</td></tr><tr><td>a<sub>0</sub> = (C<sup>(0)</sup>+Y'<sub>0</sub><sup>T<sub>312</sub></sup> · C<sup>(1)</sup>)<sup>T<sub>2341</sub></sup></td><td>T<sub>0</sub>(n<sub>1</sub>,d,d)</td><td>–</td><td>Φ<sup>*</sup><sub>l</sub><sup>[2]</sup> = ... + ∫<sub>0</sub>a<sub>0</sub> : (l ⊗ ∇l ⊗ n)da</td></tr></tbody></table>
for an arbitrary choice of the coefficients $\hat{l} \in \mathbb{R}^{n_d}$. In the right-hand side, the first term is a particular solution furnished by the pseudo-inverse $P^\dagger(m)$, and the second term is a linear combination of the column-vectors in $N_P^T(m)$ forming a basis of $\ker P(m)$, with arbitrary coefficients $(\hat{l}_i)_{1\le i \le n_d}$.

### E.3 Extended macroscopic strain vector

When the matrix $P$ is rank deficient, we append the $n_d$ coefficients $\hat{l}_i$ appearing in (E.6) to the macroscopic strain vector $l$, and write

$$l = \begin{pmatrix} \check{l} & \hat{l} \end{pmatrix} \in \mathbb{R}^{(\check{n}_l+n_d)}, \qquad (E.7)$$

where $\check{l}$ are the usual macroscopic strain vector that defines the microscopic strain $E$, see (2), referred to as $l$ in the main body of the paper, while $\hat{l}$ are the additional parameters parametrizing the solution $X$ of the rank-deficient linear problem. The dimension of $l$ is now $n_l = \check{n}_l + n_d$. We denote the injection matrices $\check{I}$ and $\hat{I}$ of $\check{l}$ and $\hat{l}$ into $l$, respectively,

$$\check{I} = \begin{pmatrix} I_{\check{n}_l} \\ 0_{n_d \times \check{n}_l} \end{pmatrix} \in \mathbb{T}^{(n_l \times \check{n}_l)} \qquad (\text{E.8a}) \qquad \hat{I} = \begin{pmatrix} 0_{\check{n}_l \times n_d} \\ I_{n_d} \end{pmatrix} \in \mathbb{T}^{(n_l \times n_d)} \qquad (\text{E.8b})$$

which enable us to rewrite $l$ as $l = \check{I} \cdot \check{l} + \hat{I} \cdot \hat{l}$. Since $\check{I} \cdot \check{I}^T$ and $\hat{I} \cdot \hat{I}^T$ are orthogonal projections from the space $\mathbb{R}^{n_l}$ in which $l$ lives onto the subspaces with equations $\hat{l} = 0$ and $\check{l} = 0$, respectively, the following identity holds,

$$\check{I} \cdot \check{I}^T + \hat{I} \cdot \hat{I}^T = I_{n_l}. \qquad (E.9)$$

To capture the fact that the macroscopic strain $E$ in (2) is a function of the original set of macroscopic degrees of freedom $\check{l}$, but not of the added $\hat{l}$ part, we require that any sub-block in $E_l(m)$, $E'_l(m)$ or $E''_l(m)$ corresponding to a range of indices $\hat{l}$ vanishes, i.e.,

$$E_l(m) \cdot \hat{I} = 0, \quad E_l^{T132}(m) \cdot \hat{I} = 0 \quad E_l^{T1432}(m) \cdot \hat{I} = 0, \qquad (\text{E.10})$$

Indeed, Equation (E.10) warrants that the strain in (2) can be rewritten in terms of $\check{l}$ and its gradients as

$$E = \check{E}_l(m) \cdot \check{l} + \check{E}'_l(m) : \nabla \check{l} + \check{E}''_l(m) : \nabla^2 \check{l} + \dots \\ + E_y(m) \cdot y + E'_y(m) : y' + E''_y(m) : y'' + \dots \qquad (\text{E.11})$$

where

$$\check{E}_l(m) = E_l(m) \cdot \check{I}, \quad \check{E}'_l(m) = (E_l^{T132}(m) \cdot \check{I})^{T132}, \quad \check{E}''_l(m) = (E_l^{T1432}(m) \cdot \check{I})^{T1432}. \qquad (\text{E.12})$$

The proof of (E.11) is left to the reader.

In view of Equation (E.10) and Table 1, $\hat{l}$ does not appear anywhere in the specification of the problem: it is a set of free parameters that are reserved for parameterizing the solution (E.6) of the rank-deficient linear problem.

### E.4 Changes to leading-order analysis

The leading-order problem (B.4) is of the form (E.4) with

$$Y = - \begin{pmatrix} W_{yl}(m(X)) \\ Q \cdot E_l(m(X)) \end{pmatrix} \cdot l(X). \qquad (\text{E.13})$$
The solvability condition (E.5) yields

$$ \left[ N_P(m) \cdot \begin{pmatrix} W_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix} \right] \cdot l(m) = 0. \qquad (E.14) $$

The vector in square brackets is an output of the homogenization procedure, representing $n_d$ conditions that are linear in the macroscopic strain $l$.

When (E.14) is satisfied, the solution is given by (E.6) as

$$ \begin{pmatrix} y^*_{[0]}(X) \\ g^*_{[0]}(X) \end{pmatrix} = -P^\dagger(m) \cdot \begin{pmatrix} W_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix} \cdot l + N_P^T(m) \cdot \hat{l}. \qquad (E.15) $$

Thanks to the definition of the extended macroscopic strain vector $l$ in (E.7), this solutions matches the form (B.6) used in the non-deficient case, provided we replace the inverse by the Moore-Penrose inverse and include a new term in the definition of $\mathcal{R}$ in (B.7),

$$ \mathcal{R}(m) = -P^{\dagger}(m) \cdot \begin{pmatrix} W_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix} + N_P^T(m) \cdot \hat{I}^T. \qquad (E.16) $$

It is convenient to rewrite this equation in a slightly different form, for a reason that will be discussed later. Using (E.10), one can show that

$$ \begin{pmatrix} W_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix} \cdot \hat{I} = 0, \qquad (E.17) $$

i.e., the operators $W_{yl}$ and $E_l$ do not sense the added degrees of freedom $\hat{l}$. Combining with (E.9), this shows that

$$ \begin{pmatrix} W_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix} = \begin{pmatrix} W_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix} \cdot \check{I} \cdot \check{I}^T. \qquad (E.18) $$

Inserting into the expression of $\mathcal{R}(m)$, we obtain

$$ \mathcal{R}(m) = -P^{\dagger}(m) \cdot \begin{pmatrix} W_{yl}(m) \\ Q \cdot E_l(m) \end{pmatrix} \cdot \check{I} \cdot \check{I}^T + N_P^T(m) \cdot \hat{I}^T. \qquad (E.19) $$

We use this expression of $\mathcal{R}(m)$ in the code and not that proposed earlier in (B.7). Indeed, the latter can be recovered as a particular case: when the matrix $P$ is invertible, $n_d = 0$, implying that $P^{\dagger}(m) = P^{-1}(m)$ and that $N_P$ and $\hat{I}$ are zero-dimension array, and the last term in (E.19) vanishes.

The definitions (B.8-B.15) of the other leading-order quantities such as $Y_0$, $G_0$, etc. are unchanged.

The following identity can be established using (B.8), (E.19), (E.3) and the orthogonality of the
projectors $\check{\mathcal{I}}^T \cdot \hat{\mathcal{I}} = 0$ which follows from (E.8),

$$ 
\begin{aligned}
Y_0 \cdot \hat{\mathcal{I}} &= \begin{pmatrix} I_{n_y} & \mathbf{0}_{n_y \times n_c} \end{pmatrix} \cdot \mathcal{R}(\mathbf{m}) \cdot \hat{\mathcal{I}} \\
&= \begin{pmatrix} I_{n_y} & \mathbf{0}_{n_y \times n_c} \end{pmatrix} \cdot \mathcal{N}_P^T(\mathbf{m}) \\
&= \begin{pmatrix} \mathcal{N}_P(\mathbf{m}) \cdot \begin{pmatrix} I_{n_y} \\ \mathbf{0}_{n_y \times n_c} \end{pmatrix} \end{pmatrix}^T \quad (E.20) \\
&= \begin{pmatrix} \mathcal{N}(\mathbf{m}) \cdot \begin{pmatrix} \mathbf{0}_{n_c \times n_y} & I_{n_c} \end{pmatrix} \cdot \begin{pmatrix} I_{n_y} \\ \mathbf{0}_{n_y \times n_c} \end{pmatrix} \end{pmatrix}^T \\
&= \mathbf{0}.
\end{aligned}
$$

As a result, the microscopic displacement $\mathbf{y}_{[0]} = Y_0(\mathbf{m}) \cdot \mathbf{l}$ can be expressed as $\mathbf{y}_{[0]} = Y_0(\mathbf{m}) \cdot (\check{\mathcal{I}} \cdot \check{\mathcal{I}}^T + \hat{\mathcal{I}} \cdot \hat{\mathcal{I}}^T) \cdot \mathbf{l} = \check{Y}_0(\mathbf{m}) \cdot \check{\mathcal{I}}$ where $\check{Y}_0(\mathbf{m}) = Y_0(\mathbf{m}) \cdot \check{\mathcal{I}}$: it depends on the original set of degrees of freedom $\check{\mathcal{I}}$ only and the matrix $Y_0$ has a zero block in the range of indices associated with $\hat{\mathcal{I}}$, $Y_0(\mathbf{m}) = (\check{Y}_0(\mathbf{m}) \mathbf{0}_{n_y \times n_c})$.

The stress $G_0 \cdot \mathbf{l}$, however, can depend on the $\hat{\mathcal{I}}$-block of $\mathbf{l}$ as well: by adapting the calculation in (E.20), one can show that $G_0 \cdot \hat{\mathcal{I}} = \mathcal{N}^T(\mathbf{m})$ is non-zero in the rank-deficient case. The components of $\hat{\mathcal{I}}$ can therefore be interpreted as the stress associated with the macroscopic kinematic constraint (E.14); this stress is akin to a Lagrange multiplier, i.e., is not set by any constitutive law.

Combining (E.20) with (B.11) and (B.12), one can show that the strain localization tensor $F_0$ and the equivalent stiffness $K_0$ are also uncoupled to $\hat{\mathcal{I}}$, implying zero $\hat{\mathcal{I}}$-sub-blocks,

$$ F_0(\mathbf{m}) = \begin{pmatrix} \check{F}_0(\mathbf{m}) & \mathbf{0}_{n_y \times n_d} \end{pmatrix}, \quad (E.21a) $$

$$ K_0(\mathbf{m}) = \begin{pmatrix} \check{K}_0(\mathbf{m}) & \mathbf{0}_{\check{n}_l \times n_d} \\ \mathbf{0}_{n_d \times \check{n}_l} & \mathbf{0}_{n_d \times n_d} \end{pmatrix}. \quad (E.21b) $$

In addition, we obtain the following identity by combining (B.14), (E.20) and (E.21a),

$$ 
\begin{aligned}
S_0 \cdot \hat{\mathcal{I}} &= \mathcal{K}(\mathbf{m}) \cdot [F_0(\mathbf{m}) \cdot \hat{\mathcal{I}}] + \mathcal{Q}^T \cdot G_0(\mathbf{m}) \cdot \hat{\mathcal{I}} \quad (E.22) \\
&= \mathcal{Q}^T \cdot \mathcal{N}^T(\mathbf{m}).
\end{aligned}
$$

The reason we prefer the expression of $\mathcal{R}$ in (E.19) to that derived first in (E.16) is that it makes it much more evident $Y_0$, $F_0$ and $K_0$ are insensitive to the added degrees of freedom $\hat{\mathcal{I}}$, see the identities (E.20-E.21).
## E.5 Changes to the energy expansion

With the help of (E.9), (E.22), (C.13), (C.12b) and (B.15a), the second term in the integrand in (C.17b) can be written as

$$ 
\begin{aligned}
(S_0(m) \cdot l) \cdot E^{[2]} &= (S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot l) \cdot E^{[2]} + (S_0(m) \cdot \hat{I} \cdot \hat{I}^T \cdot l) \cdot E^{[2]} \\
&= (S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot l) \cdot E^{[2]} + (Q^T \cdot N^T(m) \cdot \hat{I}^T \cdot l) \cdot E^{[2]} \\
&= (S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot l) \cdot E^{[2]} + (N^T(m) \cdot \hat{I}^T \cdot l) \cdot Q \cdot E^{[2]} \\
&= (S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot l) \cdot ((\mathcal{J}^{11}(m) \cdot h) :: (\nabla h \otimes \nabla h) \\
&\quad + (\mathcal{J}^2(m) \cdot h) \cdot \nabla^2 h + E'_y(m) \cdot \nabla y_{[1]}) + (E_y^T(m) \cdot S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot l) \cdot y_{[2]} \\
&= (S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot l) \cdot ((\mathcal{J}^{11}(m) \cdot h) :: (\nabla h \otimes \nabla h) + (\mathcal{J}^2(m) \cdot h) \cdot \nabla^2 h + E'_y(m) \cdot \nabla y_{[1]}),
\end{aligned}
\quad (E.23) $$

where the terms that have been crossed out are zero. In view of this, the definition of the operators $\hat{\mathcal{B}}^{(0)}$, $\hat{\mathcal{B}}^{(1)}$, $\mathcal{C}^{(0)}$, $\mathcal{C}^{(1)}$ can be modified by including the projector $\check{I} \cdot \check{I}^T$ onto the subspace $\hat{l} = 0$, to the right of $S_0(m)$,

$$ 
\begin{aligned}
\hat{\mathcal{B}}^{(0)}(m) &= ((\mathcal{J}^1(m))^{T_{4123}} \cdot \mathcal{K}(m) \cdot \mathcal{J}^1(m))^{T_{125346}} \\
&\quad + 2((\mathcal{J}^{11}(m))^{T_{612345}} \cdot S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot \mathcal{V}^l) \\
\hat{\mathcal{B}}^{(1)}(m) &= ((\mathcal{J}^1(m))^{T_{4123}} \cdot \mathcal{K}(m) \cdot E_y(m))^{T_{1243}} \\
C^{(0)}(m) &= [(\mathcal{J}^2(m))^{T_{51234}} \cdot S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot \mathcal{V}^l]^{S_{45}} \\
C^{(1)}(m) &= (E'_y(m))^{T_{312}} \cdot S_0(m) \cdot \check{I} \cdot \check{I}^T \cdot \mathcal{V}^l.
\end{aligned}
\quad (E.24) $$

The original and amended definitions in (C.25) and (E.24) are equally valid, but the latter has the advantage that it yields final tensors $a_1$ and $B_0$ having zero $\hat{l}$-sub-blocks (the sub-blocks obtained with the former set of definitions do evaluate to zero when the solvability constraints are considered but this is much less evident, and potentially confusing).

Our implementation makes use of (E.24) and not (C.25), both definitions being identical for non-singular matrices: when $n_d = 0$, we have $\check{I} \cdot \check{I}^T = I_{n_d}$ by (E.9).

## E.6 Changes to the corrective displacement

The linear problem (C.33) makes use of the same matrix $P$ as the leading-order problem. The solvability condition for $y_{[1]}$ is furnished by (E.5) as

$$ 
\left[ N_P(m) \cdot \begin{pmatrix} (\mathcal{B}^{(1)}(m(X)))^{T_{2314}} \\ Q \cdot \mathcal{J}^1(m(X)) \end{pmatrix} \right] \therefore (\nabla h(X) \otimes h(X)) = 0. \quad (E.25) $$

The quantity in square brackets is an output of the homogenization procedure that encodes $n_d$ conditions depending linearly on $l$ and $\nabla l$.

As we did earlier at the leading order, the solution $(y_{[1]}^*(X), g_{[1]}^*(X))$ of (C.33) is the sum of a particular solution furnished by the Moore-Penrose inverse and a linear combination of the null vectors of $P$ with new coefficients $\hat{l}$,

$$ 
\begin{pmatrix} y_{[1]}^*(X) \\ g_{[1]}^*(X) \end{pmatrix} = -P^\dagger(m) \cdot \begin{pmatrix} (\mathcal{B}^{(1)}(m(X)))^{T_{2314}} \\ Q \cdot \mathcal{J}^1(m(X)) \end{pmatrix} \therefore (\nabla h(X) \otimes h(X)) + N_P^T(m) \cdot \hat{l}. \quad (E.26) $$

The first term in the right-hand side is taken care of by replacing the inverse of $P(m)$ by its Moore-Penrose inverse in the definition of $\mathcal{R}'$ in (C.36), as we did earlier with $\mathcal{R}$. The second
term $N_P^T(m) \cdot \hat{l}$ in $y_{[1]}^*(X)$ and $g_{[1]}^*(X)$ adds up in a straightforward way to the term $N_P^T(m) \cdot \hat{l}$ from leading order: the full microscopic displacement $y^*(X) = y_{[0]}^*(X) + y_{[1]}^*(X) + \dots$ and Lagrange multipliers $g^*(X) = g_{[0]}^*(X) + g_{[1]}^*(X) + \dots$ obtained by summing up the $\eta^0$ and $\eta^1$ contributions are now linear combinations of the null vectors in $N_P(m)$ with coefficients $\hat{l} + \hat{\hat{l}}$ having both a leading order contribution ($\hat{l}$) and an order $\eta$ contribution ($\hat{\hat{l}}$). A simple way to deal with this complication is (i) to discard the $N_P^T(m) \cdot \hat{l}$ contribution in (E.26), and (ii) agree that $\hat{l}$ is a series in $\eta$. This avoids extending the vector $l$ with $n_d$ new entries for every homogenization order.

Concretely, we simply replace the inverse appearing in (C.36) by the Moore-Penrose inverse: our implementation uses

$$ \mathcal{R}'(m) = -P^{\dagger}(m) \cdot \begin{pmatrix} (\mathcal{B}^{(1)}(m))^{T_{2314}} \\ \mathcal{Q} \cdot \mathcal{J}^{1}(m) \end{pmatrix}. \qquad (E.27) $$

### E.7 Solvability condition for $y_{[2]}$

The quantity $y_{[2]}$ entering in

$$ E^{[2]} = [(\mathcal{J}^{11} \cdot h) :: (\nabla h \otimes \nabla h) + (\mathcal{J}^2 \cdot h) :: \nabla^2 h + E'_y : \nabla y_{[1]}] + E_y^T \cdot y_{[2]} \qquad (E.28) $$

has been eliminated from (E.23) using the constraint $Q \cdot E^{[2]} = 0$. For $y_{[2]}$ to exist, one must have

$$ 0 = Q \cdot E^{[2]} = Q \cdot [(\mathcal{J}^{11} \cdot h) :: (\nabla h \otimes \nabla h) + (\mathcal{J}^2 \cdot h) :: \nabla^2 h + E'_y : \nabla y_{[1]}] + Q \cdot E_y^T \cdot y_{[2]}. \qquad (E.29) $$

This leads to the compatibility condition

$$ Q \cdot [(\mathcal{J}^{11} \cdot h) :: (\nabla h \otimes \nabla h) + (\mathcal{J}^2 \cdot h) :: \nabla^2 h + E'_y : \nabla y_{[1]}] \in \operatorname{Im}(Q \cdot E_y^T). \qquad (E.30) $$

In the code, a basis of vectors perpendicular to $\operatorname{Im}(Q \cdot E_y^T)$ is produced using a row-reduction algorithm, and we print out the conditions that each of these vectors is perpendicular to the vector $Q \cdot [...]$ appearing in the left-hand side above: this yields conditions depending linearly on $\nabla h \otimes \nabla h \otimes h$ and $\nabla^2 h \otimes h$.

### E.8 Summary: extension to rank-deficient problems

The following extension of the code enables us to deal with a rank-deficient matrix $P(m)$:

* provide integers $\tilde{n}_l$ and $n_d$ and the injection matrix $\tilde{I}$ as an optional arguments to the homogenization procedure and check the condition (E.10) on the tensors $E_l$, $E'_l$ and $E''_l$ passed in argument;

* compute a set of null vectors of the symmetric matrix $P(m)$, check that there are $n_d$ such vectors and that they are of the form (E.2), compute the Moore-Penrose inverse $P^{\dagger}(m)$ if $n_d > 0$;

* return the solvability conditions (E.14), (E.25) and (E.30) whenever $n_d > 0$;

* replace Equations (B.7), (C.25) and (C.36) yielding $\mathcal{R}, \hat{\mathcal{B}}^{(0)}, \hat{\mathcal{B}}^{(1)}, \mathcal{C}^{(0)}, \mathcal{C}^{(1)}$ and $\mathcal{R}'$ with their extensions (E.19), (E.24) and (E.27)

### E.9 Illustration: a truss lattice with inextensible beams

We consider again the truss lattice shown in Figure 2, but assume this time that the bars on the upper side (+) are inextensible: the corresponding stretching strain is constrained to be zero, $E_4 = 0$. In view of (33) and (29), the inextensibility condition for the upper beams can be rewritten as

$$ e(S) - \frac{1}{2}c(S) = 0 + \mathcal{O}(\eta). \qquad (E.31) $$
This is a constraint that applies to the macroscopic strain $l = (e, c)$: this contrasts with the illustrations presented in Section 4, where the macroscopic strain was unconstrained. The linear problems encountered in the homogenization become rank-deficient, as we show now by running the homogenization procedure.

We first run the homogenization procedure with similar input parameters as those listed in Table 4, except that

* we limit attention to uniform spring constants, setting $m = 0$ and $n_m = 0$ and $k(S) = k$;

* we add the inextensibility constraint by incrementing $n_c = 2 + 1 = 3$ and defining $Q(S) = \delta_1^3 \otimes \delta_6^7 + \delta_2^3 \otimes \delta_7^7 + \delta_3^3 \otimes \delta_4^7$, where the first two terms represent the zero-average conditions $\langle y_1^\pm \rangle = 0$ and $\langle y_2^\pm \rangle = 0$ already used earlier, and the last term encodes the new constraint $E_4 = 0$ (we using the notation introduced in Table 4 for unit vectors);

* optionally, we drop the term $i = 4$ in the assembly of the stiffness matrix $K$ in Table 4—this does not matter as the corresponding strain is zero anyway.

When run with these parameters, the homogenization returns an error message indicating that the linear problem has a rank deficiency of 1, and that a single coefficient $\hat{l} = (\hat{l}_1)$ must be appended to the list of macroscopic parameters $l = (e, c)$, see (E.7).

We therefore modify the input further as follows:

* with $l = (e, c, \hat{l}_1)$, we now set $n_l = 3$, and we extend the dimensions of the tensors $E_l, E'_l, E''_l$ accordingly, by filling them up with zeros;

* we include the optional argument $\hat{I} = \delta_3^3 \otimes \delta_1^1$ describing the injection $\hat{l}$ into $l$, see (E.8b). When run with these parameters, the homogenization procedure returns

* the solvability condition for $y_{[0]}$, see (E.14), in the form that we had announced in (E.31)

$$e(S) - \frac{1}{2}c(S) = 0, \qquad (E.32)$$

* the solvability condition for $y_{[1]}$ in (E.25) in the form $0 = 0$ (it is automatically satisfied),

* the homogenized energy is obtained by interpreting the output of the code as earlier in Section 4, except that we can now eliminate $e(S) = \frac{1}{2}c(S)$ in favor of $c(S)$ using the solvability condition: the result is

$$\Phi^*[c] = \int_{-\infty}^{+\infty} \left[ \frac{13k}{24}c^2 - \frac{k\eta^2}{4}c'^2 + O(\eta^3) \right] \frac{dS}{\eta}. \qquad (E.33)$$

## References

Abali, B. E. and E. Barchiesi (2021). Additive manufacturing introduced substructure and computational determination of metamaterials parameters by means of the asymptotic homogenization. *Continuum Mechanics and Thermodynamics* 33(4):993–1009. [DOI], [OA].

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

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

Abdoul-Anziz, H., P. Seppecher, and C. Bellis (2019). Homogenization of frame lattices leading to second gradient models coupling classical strain and strain-gradient terms. *Mathematics and Mechanics of Solids* 24(12):3976–3999. [DOI], [HAL].

Audoly, B. (2023). *The shoal library*. Software Heritage persistent IDentifiers (SWHIDs) permlink [SWHID]. Version 1.0. Git repository [GIT].

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

Bacigalupo, A. (2014). Second-order homogenization of periodic materials based on asymptotic approximation of the strain energy: formulation and validity limits. *Meccanica* 49(6):1407–1425. [DOI], [ARXIV].

Bakhvalov, N. and G. Panasenko (1989). *Homogenisation: averaging processes in periodic media: mathematical problems in the mechanics of composite materials*. Springer. [DOI], [OA].
Bažant, Z. and M. Christensen (1972). Analogy between micropolar continuum and grid frame-
works under initial stress. *International Journal of Solids and Structures* 8(3):327–346. [DOI],
[OA].

Berdichevskii, V. (1981). On the energy of an elastic rod. *Journal of Applied Mathematics and Mechanics* 45(4):518–529. [DOI].

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

Boutin, C. (2019). Homogenization Methods and Generalized Continua in Linear Elasticity. *Encyclopedia of Continuum Mechanics*. Ed. by H. Altenbach and A. Öchsner. Springer, pp 1–35. [DOI].

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

Cioranescu, D. and P. Donato (1999). *An Introduction to Homogenization*. Vol. 17. Oxford University Press. ISBN: 9780198565543.

Cioranescu, D. and J. S. J. Paulin (1999). *Homogenization of Reticulated Structures*. Vol. 136. Springer. ISBN: 9781461274377. [DOI].

David, M., J.-J. Marigo, and C. Pideri (2012). Homogenized interface model describing inhomogeneities located on a surface. *Journal of Elasticity* 109(2):153–187. [DOI], [HAL].

Dos Reis, F. and J. Ganghoffer (2012). Construction of micropolar continua from the asymptotic homogenization of beam lattices. *Computers & Structures* 112-113:354–363. [DOI].

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

Gambin, B. and E. Kröner (1989). Higher-order terms in the homogenized stress-strain relation of periodic elastic media. *physica status solidi (b)* 151(2):513–519. [DOI].

Hans, S. and C. Boutin (2008). Dynamics of discrete framed structures: a unified homogenized description. *Journal of Mechanics of Materials and Structures* 3(9):1709–1739. [DOI], [OA].

Hodges, D. H. (2006). *Nonlinear composite beam theory*. Vol. 213. Progress in astronautics and aeronautics. American Institute of Aeronautics and Astronautics. [DOI].

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

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

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

Sanchez-Palencia, E. (1980). Non-homogeneous media and vibration theory. *Lecture Note in Physics* 320:57–65. [DOI], [OA].

Smyshlyaev, V. and K. Cherednichenko (2000). On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media. *Journal of the Mechanics and Physics of Solids* 48(6-7):1325–1357. [DOI].

Wolfram Research, Inc. (2021). *Mathematica, Version 13.0.0*. Champaign, IL, 2021.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long

as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** BA implemented the library. Both authors contributed equally to all other aspects of this work.

**Supplementary Material** The homogenization method presented in this paper is distributed in the form of a library, named 'shoal', available at SoftwareHeritage https://archive.softwareheritage.org/swh:1:dir:788abad82a752c44def322f4bb6c64bb63367f7e;origin=https://git.renater.fr/anonscm/git/shoal/shoal.git;visit=swh:1:snp:5013a4223e0901d9c6aa52da15e7doff9b5606e4;anchor=swh:1:rev:a48502c10dd2a54b3fa83acb245632fb320d4ee9. The git repository of the library can be accessed through GIT https://git.renater.fr/anonscm/git/shoal/shoal.git. It is implemented in the symbolic calculation language Wolfram Mathematica, and distributed under CeCILL v2.1 license. The code used to produce the examples presented in this paper is included in the library.

**Acknowledgements** We are grateful to Yang Ye for his useful feedback on the technical document that we used internally for this project.

**Funding** Not applicable.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.