Near the physical edges of the lattice, boundary layers form and the homogenization breaks down [Dum86; LM18]. Until we properly address these layers, they prevent us from making accurate global predictions of the equilibrium solution based on the homogenized theory (4.5)—the analysis of the boundary layers is a significant piece of work and we leave it to a follow-up paper. In the meantime, we give up on making predictions based on the stand-alone homogenized theory (4.5) and keep it tethered to numerical simulations of the full lattice. Our verification method proceeds by extracting the microscopic fluctuations from full numerical simulations of the lattice, and by comparing them to the predictions of the homogenization procedure. This approach is purely local, warranting that the systematic errors in the boundary layers do not interfere with the verification in the bulk. Concretely, our verification method produces color maps of the lattice showing where the homogenization method is accurate and where it is not: it is considered to be valid if the accuracy is good everywhere but in a neighborhood of the boundaries, in the sense that the residual error goes sufficiently rapidly to zero with the cell size (see the convergence plot in Figure 6.4).

The verification procedure builds up on the explicit construction of the continuous fields proposed in Remark 3.2. We start from the nodal displacement and rotation ($\bar{v}_{\beta}, \theta_{\beta}) \in \mathbb{R}^{n_a}$ obtained numerically by running discrete lattice simulations based on the energy (2.23). The procedure outlined in Remark 3.2 is implemented through a series of interpolation and averaging steps detailed in Section 6.1 below, delivering interpolations for the macroscopic displacement and the microscopic degrees of freedom $\bar{y}(\mathbf{X}) = (\xi_1(\mathbf{X}) / \eta, \psi_1(\mathbf{X}), \dots)$ in the form of continuous functions of $\mathbf{X}$. Here, the overbar marks quantities extracted from the numerical simulation. By differentiating $\bar{\dot{u}}(\mathbf{X})$ symbolically, one can obtain the macroscopic strain $\bar{l}(\mathbf{X})$ and its gradient $\nabla \bar{l}(\mathbf{X})$. The accuracy of the homogenization procedure is assessed by comparing locally the microscopic degrees of freedom $\bar{y}(\mathbf{X})$ that are directly extracted from the simulations, to the prediction $\mathbf{y} = \mathbf{Y} \cdot \bar{l} + \mathbf{Y}' \cdot \nabla \bar{l}$ in (4.3) based on the localization tensors $\mathbf{Y}$ and $\mathbf{Y}'$ obtained via homogenization, combined with the interpolated macroscopic strain $\bar{l}(\mathbf{X})$ and its gradient $\nabla \bar{l}(\mathbf{X})$.

The verification is therefore focussed on the prediction (4.3) of the microscopic displacement $\mathbf{y}$. This is indeed the key prediction of the homogenization method as it slaves the microscopic displacement $\mathbf{y}$ to the macroscopic degrees of freedom. The expression of the homogenized energy (4.5) follows straightforwardly by inserting this expression of $\mathbf{y}$ into the original energy.

## 6.1. Verification procedure

As a first illustration of the verification procedure, we consider a 2D rectangular strip made of a honeycomb lattice subject to compression, as shown in Figure 6.1. The length of the strip is $L=1$ and its width is $\delta = 0.36$. The displacements of the nodes at the left vertical boundary are blocked, $v_{\beta \in lb} = 0$. Those at the right boundary are displaced by an amount $\Delta$, $v_{\beta \in rb} = -\Delta e_1$, where $e_1$ is a unit vector along the long axis of the strip. The rotations are blocked along both short edges, $\theta_{\beta \in lb} = \theta_{\beta \in rb} = 0$. We check convergence for finer and finer lattices, when the beam length $\ell = O(\eta)$ goes to zero ($\ell \to 0$).

The lattice is made up of extensible 2D beams, modeled with the discrete energy (2.19–2.20). All of them are assigned a bending modulus $EI = 1$ and a stretching modulus $EA = 10^{12}$; this warrants that the aspect-ratio parameter $\chi = EI/EA\ell^2$ remains at most of order $10^{-10}$ in all our simulations, so that all beams are effectively in the inextensible regime. The simulations results are compared to the inextensible homogenization theory (§5.3–5.4).

<figure><img src="image_7.png" /><figcaption>**Figure 6.1.** Honeycomb strip in compression: the lattice is loaded by bringing the short edges closer to one another by a distance $\Delta$. Sliding and rotation is blocked on both short edges. We check convergence of the discrete numerical solution when the beam length goes to zero, $\ell \to 0$.</figcaption></figure>

For a given beam length $\ell$, the discrete simulation yields a set of nodal displacements $\bar{v}_{\beta}$ and rotations $\bar{\theta}_{\beta}$. For each beam $\alpha$ in the lattice and for every Bravais sub-lattice $b \in \{1, 2\}$, we pick $N_b^a = 6$ nodes from the sub-lattice $b$ in the neighborhood of the beam $\alpha$, and use the corresponding nodal values $\bar{v}_{\beta}$ and $\bar{\theta}_{\beta}$ to build second-order interpolations $\bar{v}_b^a(\mathbf{X})$ and $\bar{\theta}_b^a(\mathbf{X})$ (in general, we need $N_b^a = 1 + d + d(d+1)/2 = (d+1)(d+2)/2$ neighbors to build the second-order interpolations in dimension $d$). These interpolations applied to the sub-lattice $b$ are valid in the neighborhood of the element $\alpha$. In Figure 6.2b, we plot the interpolation of the longitudinal displacement $\bar{u}_b^a(X_1) = \bar{v}_b^a(X_1) \cdot e_1$ for the particular case of simple compression (uniaxial stress): the interpolations $\bar{u}_1^a(X_1)$ and
$\tilde{u}_2^a(X_1)$ on the two Bravais sub-lattices are different (dashed lines), which points to the existence of a microscopic displacement.

<figure><img src="image_6.png" /><figcaption>**Figure 6.2.** Local interpolations $\tilde{u}_b^a(\mathbf{X}) = \tilde{v}_b^a(\mathbf{X}) \cdot e_1$ of the longitudinal displacement, illustrated in the particular case of uniaxial compression of a honeycomb lattice. The state of uniaxial compression is obtained by modifying the boundary conditions on the short edges of the lattice shown in Figure 6.1 to let the nodes to slide freely in the transverse direction. (a) Schematic view of the deformed lattice. The link $a$ about which the interpolations are constructed is highlighted in blue. Only the $N_\gamma^a = 6$ neighboring nodes from each Bravais sub-lattice $b \in \{1,2\}$ that enter into the interpolation are shown. (b) Longitudinal displacement from discrete simulation ($\tilde{u}_a$, open and filled disks), interpolated in each Bravais sub-lattice ($\tilde{u}_b^a(X_1)$, dashed lines) and macroscopic displacement ($\tilde{u}^a(X_1)$, solid line). In this simple illustration, the $X_2$ direction is an invariance direction and all displacements are functions of $X_1$ only.</figcaption></figure>

Following Remark 3.2, we obtain a local interpolation $\tilde{u}^a(\mathbf{X})$ of the macroscopic displacement by averaging over the Bravais sub-lattices,

$$ \tilde{u}^a(\mathbf{X}) = \frac{1}{n_b} \sum_{b=1}^{n_b} \tilde{v}_b^a(\mathbf{X}). \qquad (6.1) $$

Differentiating symbolically with respect to $\mathbf{X}$, we obtain the displacement gradient $\nabla \tilde{u}^a(\mathbf{X})$, from which we extract the local interpolations for the macroscopic rotation $\tilde{\gamma}^a(\mathbf{X}) = -\frac{1}{2} \nabla \tilde{u}^a(\mathbf{X}) : \mathcal{N}$ and the macroscopic strain $\tilde{\epsilon}^a(\mathbf{X}) = \nabla \tilde{u}^a(\mathbf{X}) : \mathcal{M}_b$ using (3.8). The macroscopic state vector is then defined as $\tilde{l}^a(\mathbf{X}) = (\tilde{u}^a(\mathbf{X}) \ \tilde{\gamma}^a(\mathbf{X}) \ \tilde{\epsilon}^a(\mathbf{X}))$ or $\tilde{l}^a(\mathbf{X}) = (\tilde{\epsilon}^a(\mathbf{X}))$ depending on whether we use the standard definition (3.10) or the alternate one in (3.24).

Next, the local microscopic displacement and rotation associated to each Bravais sub-lattice $b$ ($1 \le b \le n_b$) are extracted using (3.4) as

$$ \tilde{\xi}_b^a(\mathbf{X}) = \tilde{v}_b^a(\mathbf{X}) - \tilde{u}^a(\mathbf{X}) \qquad \tilde{\psi}_b^a = \tilde{\theta}_b^a(\mathbf{X}) - \tilde{\gamma}^a(\mathbf{X}). \qquad (6.2) $$

We finally collect all the microscopic degrees of freedom following (3.15), and evaluate the interpolation at the center $\mathbf{X}_a^c$ of the element being considered,

$$ \tilde{y}^a = \left( \frac{\tilde{\xi}_1^a(\mathbf{X}_a^c)}{\eta}, \tilde{\psi}_1^a(\mathbf{X}_a^c), \dots, \frac{\tilde{\xi}_{n_b}^a(\mathbf{X}_a^c)}{\eta}, \tilde{\psi}_{n_b}^a(\mathbf{X}_a^c) \right). \qquad (6.3) $$

We compare this to the leading order $(\mathbf{y}^*)^a_{[0]}$ and first correction $(\mathbf{y}^*)^a_{[1]}$ predicted by the homogenized results (5.10) for an inextensible lattice,

$$ \begin{aligned} (\mathbf{y}^*)^a_{[0]} &= \frac{a}{4} ( + \mathcal{J} : \tilde{\epsilon}^a \ 0 - \mathcal{J} : \tilde{\epsilon}^a \ 0 ), \\ (\mathbf{y}^*)^a_{[1]} &= \frac{3a\eta}{4} ( \mathbf{0}_2 \ \mathrm{curl}(\mathcal{J} : \tilde{\epsilon}^a) \ \mathbf{0}_2 - \mathrm{curl}(\mathcal{J} : \tilde{\epsilon}^a) ). \end{aligned} \qquad (6.4) $$

In the right-hand sides, we use the estimate $\tilde{\epsilon}^a$ of the strain, and $\nabla \tilde{\epsilon}^a$ of the strain gradient at the element center $\mathbf{X}_a^c$, as obtained from the discrete numerical simulation: these values are extracted from the interpolations $\tilde{l}^a(\mathbf{X}_a^c)$ and $\nabla \tilde{l}^a(\mathbf{X}_a^c)$, based on the definition (3.10) or (3.24) of $\mathbf{l}$ in terms of the strain $\tilde{\epsilon}$ and the Mandel convention (3.7) for representing the strain $\epsilon$ as a flat vector $\tilde{\epsilon}$.

Specifically, we compute at every element $a$ the error $e^a_{[0]}$ on the microscopic displacement based on leading-order homogenization, the error $e^a_{[1]}$ based on second-order homogenization, as well as its magnitude $e^a_{[\times]}$ in the raw output of numerical simulations as

$$ \begin{aligned} e^a_{[\times]} &= \left\| \frac{\tilde{y}^a}{a \Delta} \right\|, \\ e^a_{[0]} &= \left\| \frac{\tilde{y}^a - (\mathbf{y}^*)^a_{[0]}}{a \Delta} \right\|, \\ e^a_{[1]} &= \left\| \frac{\tilde{y}^a - ((\mathbf{y}^*)^a_{[0]} + (\mathbf{y}^*)^a_{[1]})}{a \Delta} \right\|. \end{aligned} \qquad (6.5) $$

where $a = \mathcal{O}(\eta^0)$ is the conventional length of the beams in the topological lattice and $\Delta$ is the conventional magni-
$\tilde{u}_2^a(X_1)$ on the two Bravais sub-lattices are different (dashed lines), which points to the existence of a microscopic displacement.

<figure><img src="image_8.png" /><figcaption>**Figure 6.2.** Local interpolations $\tilde{u}_b^a(\mathbf{X}) = \tilde{v}_b^a(\mathbf{X}) \cdot e_1$ of the longitudinal displacement, illustrated in the particular case of uniaxial compression of a honeycomb lattice. The state of uniaxial compression is obtained by modifying the boundary conditions on the short edges of the lattice shown in Figure 6.1 to let the nodes to slide freely in the transverse direction. (a) Schematic view of the deformed lattice. The link $a$ about which the interpolations are constructed is highlighted in blue. Only the $N_\gamma^a = 6$ neighboring nodes from each Bravais sub-lattice $b \in \{1,2\}$ that enter into the interpolation are shown. (b) Longitudinal displacement from discrete simulation ($\tilde{u}_a$, open and filled disks), interpolated in each Bravais sub-lattice ($\tilde{u}_b^a(X_1)$, dashed lines) and macroscopic displacement ($\tilde{u}^a(X_1)$, solid line). In this simple illustration, the $X_2$ direction is an invariance direction and all displacements are functions of $X_1$ only.</figcaption></figure>

Following Remark 3.2, we obtain a local interpolation $\tilde{u}^a(\mathbf{X})$ of the macroscopic displacement by averaging over the Bravais sub-lattices,

$$ \tilde{u}^a(\mathbf{X}) = \frac{1}{n_b} \sum_{b=1}^{n_b} \tilde{v}_b^a(\mathbf{X}). \qquad (6.1) $$

Differentiating symbolically with respect to $\mathbf{X}$, we obtain the displacement gradient $\nabla \tilde{u}^a(\mathbf{X})$, from which we extract the local interpolations for the macroscopic rotation $\tilde{\gamma}^a(\mathbf{X}) = -\frac{1}{2} \nabla \tilde{u}^a(\mathbf{X}) : \mathcal{N}$ and the macroscopic strain $\tilde{\epsilon}^a(\mathbf{X}) = \nabla \tilde{u}^a(\mathbf{X}) : \mathcal{M}_b$ using (3.8). The macroscopic state vector is then defined as $\tilde{l}^a(\mathbf{X}) = (\tilde{u}^a(\mathbf{X}) \ \tilde{\gamma}^a(\mathbf{X}) \ \tilde{\epsilon}^a(\mathbf{X}))$ or $\tilde{l}^a(\mathbf{X}) = (\tilde{\epsilon}^a(\mathbf{X}))$ depending on whether we use the standard definition (3.10) or the alternate one in (3.24).

Next, the local microscopic displacement and rotation associated to each Bravais sub-lattice $b$ ($1 \le b \le n_b$) are extracted using (3.4) as

$$ \tilde{\xi}_b^a(\mathbf{X}) = \tilde{v}_b^a(\mathbf{X}) - \tilde{u}^a(\mathbf{X}) \qquad \tilde{\psi}_b^a = \tilde{\theta}_b^a(\mathbf{X}) - \tilde{\gamma}^a(\mathbf{X}). \qquad (6.2) $$

We finally collect all the microscopic degrees of freedom following (3.15), and evaluate the interpolation at the center $\mathbf{X}_a^c$ of the element being considered,

$$ \tilde{y}^a = \left( \frac{\tilde{\xi}_1^a(\mathbf{X}_a^c)}{\eta}, \tilde{\psi}_1^a(\mathbf{X}_a^c), \dots, \frac{\tilde{\xi}_{n_b}^a(\mathbf{X}_a^c)}{\eta}, \tilde{\psi}_{n_b}^a(\mathbf{X}_a^c) \right). \qquad (6.3) $$

We compare this to the leading order $(\mathbf{y}^*)^a_{[0]}$ and first correction $(\mathbf{y}^*)^a_{[1]}$ predicted by the homogenized results (5.10) for an inextensible lattice,

$$ \begin{aligned} (\mathbf{y}^*)^a_{[0]} &= \frac{a}{4} ( + \mathcal{J} : \tilde{\epsilon}^a \ 0 - \mathcal{J} : \tilde{\epsilon}^a \ 0 ), \\ (\mathbf{y}^*)^a_{[1]} &= \frac{3a\eta}{4} ( \mathbf{0}_2 \ \mathrm{curl}(\mathcal{J} : \tilde{\epsilon}^a) \ \mathbf{0}_2 - \mathrm{curl}(\mathcal{J} : \tilde{\epsilon}^a) ). \end{aligned} \qquad (6.4) $$

In the right-hand sides, we use the estimate $\tilde{\epsilon}^a$ of the strain, and $\nabla \tilde{\epsilon}^a$ of the strain gradient at the element center $\mathbf{X}_a^c$, as obtained from the discrete numerical simulation: these values are extracted from the interpolations $\tilde{l}^a(\mathbf{X}_a^c)$ and $\nabla \tilde{l}^a(\mathbf{X}_a^c)$, based on the definition (3.10) or (3.24) of $\mathbf{l}$ in terms of the strain $\tilde{\epsilon}$ and the Mandel convention (3.7) for representing the strain $\epsilon$ as a flat vector $\tilde{\epsilon}$.

Specifically, we compute at every element $a$ the error $e^a_{[0]}$ on the microscopic displacement based on leading-order homogenization, the error $e^a_{[1]}$ based on second-order homogenization, as well as its magnitude $e^a_{[\times]}$ in the raw output of numerical simulations as

$$ \begin{aligned} e^a_{[\times]} &= \left\| \frac{\tilde{y}^a}{a \Delta} \right\|, \\ e^a_{[0]} &= \left\| \frac{\tilde{y}^a - (\mathbf{y}^*)^a_{[0]}}{a \Delta} \right\|, \\ e^a_{[1]} &= \left\| \frac{\tilde{y}^a - ((\mathbf{y}^*)^a_{[0]} + (\mathbf{y}^*)^a_{[1]})}{a \Delta} \right\|. \end{aligned} \qquad (6.5) $$

where $a = \mathcal{O}(\eta^0)$ is the conventional length of the beams in the topological lattice and $\Delta$ is the conventional magni-
tude of the applied load in the numerical simulations. The denominators in right-hand sides in (6.5) make the errors $e_{[\cdot]}^\alpha$'s independent of both these conventions. An interpolation present in the right-hand sides is evaluated at the center $X_\alpha^c$ of the element $\alpha$.

In Figure 6.3, these errors are painted on the deformed lattice, for two different beam lengths $\ell$, using a logarithmic colormap. We observe a significant and consistent reduction on the magnitude of the error when we go from no homogenization ($e_{[\times]}^\alpha$) to leading-order homogenization ($e_{[0]}^\alpha$), and then to second-order homogenization ($e_{[1]}^\alpha$). As anticipated earlier, this improvement takes place everywhere in the lattice, except in a layer that is a few cells thick at the boundaries, where it does not decrease noticeably.

<figure><img src="image_7.png" /><figcaption>**Figure 6.3.** Honeycomb strip in compression: error $e_{[i]}^\alpha$ on the microscopic displacement predicted by the leading-order ($i=0$) or second-order ($i=1$) homogenization, compared with its magnitude ($i=\times$) in the raw numerical solution, for two different lattices (coarser one in the left column, finer one in the right column). The error is painted over the deformed lattice as a base-10-logarithmic colormap. An arbitrary magnification factor is applied on the displacement for plotting the deformed lattice. The numerical simulations were produced with an imposed displacement $\Delta = 0.1$ but the results are rescaled in such a way that this value of $\Delta$ does not affect the colors. The consistent decrease of the error in logarithmic scale with the order of homogenization is a sign that the homogenization is correct.</figcaption></figure>

This first visual indication that the homogenization is correct can be confirmed by a convergence test with respect to the beam length ($\ell \to 0$). We select a point in the interior of the rectangular domain, shown by the blue cross in the inset in Figure 6.4, and we compute the local errors $e_{[0]}^\alpha$ and $e_{[1]}^\alpha$ at the nearest element center $X_\alpha^c$ for finer and finer
tude of the applied load in the numerical simulations. The denominators in right-hand sides in (6.5) make the errors $e_{[\cdot]}^\alpha$'s independent of both these conventions. An interpolation present in the right-hand sides is evaluated at the center $\mathbf{X}_\alpha^c$ of the element $\alpha$.

In Figure 6.3, these errors are painted on the deformed lattice, for two different beam lengths $\ell$, using a logarithmic colormap. We observe a significant and consistent reduction on the magnitude of the error when we go from no homogenization ($e_{[\times]}^\alpha$) to leading-order homogenization ($e_{[0]}^\alpha$), and then to second-order homogenization ($e_{[1]}^\alpha$). As anticipated earlier, this improvement takes place everywhere in the lattice, except in a layer that is a few cells thick at the boundaries, where it does not decrease noticeably.

<figure><img src="image_9.png" /><figcaption>**Figure 6.3.** Honeycomb strip in compression: error $e_{[i]}^\alpha$ on the microscopic displacement predicted by the leading-order ($i=0$) or second-order ($i=1$) homogenization, compared with its magnitude ($i=\times$) in the raw numerical solution, for two different lattices (coarser one in the left column, finer one in the right column). The error is painted over the deformed lattice as a base-10-logarithmic colormap. An arbitrary magnification factor is applied on the displacement for plotting the deformed lattice. The numerical simulations were produced with an imposed displacement $\Delta = 0.1$ but the results are rescaled in such a way that this value of $\Delta$ does not affect the colors. The consistent decrease of the error in logarithmic scale with the order of homogenization is a sign that the homogenization is correct.</figcaption></figure>

This first visual indication that the homogenization is correct can be confirmed by a convergence test with respect to the beam length ($\ell \to 0$). We select a point in the interior of the rectangular domain, shown by the blue cross in the inset in Figure 6.4, and we compute the local errors $e_{[0]}^\alpha$ and $e_{[1]}^\alpha$ at the nearest element center $\mathbf{X}_\alpha^c$ for finer and finer
lattices. The numerical results, shown in Figure 6.4, reveal the following scaling laws: $e_{[0]}^{\alpha} = O(\ell^1)$ and $e_{[1]}^{\alpha} = O(\ell^2)$ for $\ell \to 0$. The fact that the order of the error increases from 1 to 2 when we switch from leading-order to second-order homogenization is a numerical evidence that the latter is asymptotically correct. A similar behavior as in Figure 6.4 is obtained when a different target point is chosen (data not shown). For points close to the boundary, the scaling regime seen in Figure 6.4 does not start until $\ell$ becomes significantly smaller than the distance to boundary (data not shown).

<figure><img src="image_8.png" /><figcaption>**Figure 6.4.** Honeycomb strip in compression: log-log plot of the homogenization errors $e_{[0]}^{\alpha}$ and $e_{[1]}^{\alpha}$ computed at the nearest element $\alpha$ to a fixed target (blue cross in the inset) as a function of the 'microscopic' element length $\ell$. The blue dashed lines are linear fits with slope 1 and 2, respectively.</figcaption></figure>

## 6.2. Sheared strip, with or without a hole

Next, we test the rectangular strip when subjected to a shearing displacement $\Delta$, as shown in Figure 6.5. The displacement is still blocked along the short boundary on the left-hand side ($v_{\beta \in lb} = 0$) and the nodes along the short boundary on the right-hand side are moved vertically ($v_{\beta \in rb} = +\Delta e_2$). In addition, the nodes along both short boundaries can now rotate freely.

<figure><img src="image_9.png" /><figcaption>**Figure 6.5.** Honeycomb strip subjected shear: loading geometry. The nodes on both short edges can now rotate freely.</figcaption></figure>

The same testing procedure is applied. The colormaps of the homogenization errors are shown in logarithmic scale in Figure 6.6 for two different beam lengths (a base-10 logarithm is used throughout in the paper). Here again, the higher-order correction significantly increases the accuracy of the prediction, except in a small layer along the boundaries.
lattices. The numerical results, shown in Figure 6.4, reveal the following scaling laws: $e_{[0]}^{\alpha} = O(\ell^1)$ and $e_{[1]}^{\alpha} = O(\ell^2)$ for $\ell \to 0$. The fact that the order of the error increases from 1 to 2 when we switch from leading-order to second-order homogenization is a numerical evidence that the latter is asymptotically correct. A similar behavior as in Figure 6.4 is obtained when a different target point is chosen (data not shown). For points close to the boundary, the scaling regime seen in Figure 6.4 does not start until $\ell$ becomes significantly smaller than the distance to boundary (data not shown).

<figure><img src="image_10.png" /><figcaption>**Figure 6.4.** Honeycomb strip in compression: log-log plot of the homogenization errors $e_{[0]}^{\alpha}$ and $e_{[1]}^{\alpha}$ computed at the nearest element $\alpha$ to a fixed target (blue cross in the inset) as a function of the 'microscopic' element length $\ell$. The blue dashed lines are linear fits with slope 1 and 2, respectively.</figcaption></figure>

## 6.2. Sheared strip, with or without a hole

Next, we test the rectangular strip when subjected to a shearing displacement $\Delta$, as shown in Figure 6.5. The displacement is still blocked along the short boundary on the left-hand side ($v_{\beta \in lb} = 0$) and the nodes along the short boundary on the right-hand side are moved vertically ($v_{\beta \in rb} = +\Delta e_2$). In addition, the nodes along both short boundaries can now rotate freely.

<figure><img src="image_11.png" /><figcaption>**Figure 6.5.** Honeycomb strip subjected shear: loading geometry. The nodes on both short edges can now rotate freely.</figcaption></figure>

The same testing procedure is applied. The colormaps of the homogenization errors are shown in logarithmic scale in Figure 6.6 for two different beam lengths (a base-10 logarithm is used throughout in the paper). Here again, the higher-order correction significantly increases the accuracy of the prediction, except in a small layer along the boundaries.
$\ell = 0.0130$

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

$\ell = 0.00658$

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

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

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

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

<figure><img src="image_15.png" /><figcaption>**Figure 6.6.** Honeycomb strip subjected to shear: homogenization error for two different lattices (columns), using the same conventions as in Figure 6.3. The boundary conditions are represented in a simplified way, the reader is referred to Figure 6.5 for details.</figcaption></figure>

In the next test, we punch the lattice by removing a disk-like region with radius $r = 0.04$ in its center. The boundary conditions are identical to those in Figure 6.5. This relatively small hole causes stress concentration, visible on the raw output of the numerical simulations in the first row in Figure 6.7. The typical flower-like stress concentration pattern survives to leading-order homogenization, see the second row. When higher-order homogenization is used (third row), however, it almost disappears from the error map, indicating that second-order homogenization accurately captures the fast variations of the stress in the vicinity of the hole: the error goes to zero quickly away from the hole. The existence of a ‘safe zone’ surrounding the hole (extended blue region in this third row, where the relative error is as small as ~$10^{-3}$) suggests that the punched strip could be very accurately represented by combining effective boundary conditions along the edge of the hole with a higher-order continuum model for the strip. This is probably not possible with leading-order homogenization (second row) as the region of influence of the hole stretches all the way to the lateral boundaries.
$\ell = 0.0130$

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

$\ell = 0.00658$

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

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

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

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

<figure><img src="image_17.png" /><figcaption>**Figure 6.6.** Honeycomb strip subjected to shear: homogenization error for two different lattices (columns), using the same conventions as in Figure 6.3. The boundary conditions are represented in a simplified way, the reader is referred to Figure 6.5 for details.</figcaption></figure>

In the next test, we punch the lattice by removing a disk-like region with radius $r = 0.04$ in its center. The boundary conditions are identical to those in Figure 6.5. This relatively small hole causes stress concentration, visible on the raw output of the numerical simulations in the first row in Figure 6.7. The typical flower-like stress concentration pattern survives to leading-order homogenization, see the second row. When higher-order homogenization is used (third row), however, it almost disappears from the error map, indicating that second-order homogenization accurately captures the fast variations of the stress in the vicinity of the hole: the error goes to zero quickly away from the hole. The existence of a ‘safe zone’ surrounding the hole (extended blue region in this third row, where the relative error is as small as ~$10^{-3}$) suggests that the punched strip could be very accurately represented by combining effective boundary conditions along the edge of the hole with a higher-order continuum model for the strip. This is probably not possible with leading-order homogenization (second row) as the region of influence of the hole stretches all the way to the lateral boundaries.
<figure><img src="image_16.png" /><figcaption>**Figure 6.7.** Honeycomb strip with a hole, subjected to shear: homogenization error for two different lattices (columns), using the same conventions as in Figure 6.3 and 6.5. Same boundary conditions as in Figure 6.5.</figcaption></figure>

**6.3. Cracked strip in tension**

Before moving on to more complex examples, we consider the cracked version of the rectangular honeycomb strip shown in Figure 6.8, having an aspect-ratio $\delta/L = 0.78$. The crack is obtained by removing the beams that intersect the segment with length $a$, making an angle $\pi/6$ with the direction $e_1$, and whose midpoint is at the center of the lattice. Crack propagation is not considered. The same boundary conditions are used as for the shearing test, but a longitudinal displacement $v_{\beta \in rb} = \Delta e_1$ is imposed on the right-hand-side boundary.

<figure><img src="image_17.png" /><figcaption>**Figure 6.8.** Cracked honeycomb strip in tension.</figcaption></figure>

Colormaps of the homogenization error are shown in Figure 6.9 for $\ell = 0.0125$. Similar to what has been obtained for a punched lattice, some stress concentration is observed at both crack tips. When second-order homogenization is used, not only does the homogenization error drop significantly, it is also very small except in the intermediate
<figure><img src="image_18.png" /><figcaption>**Figure 6.7.** Honeycomb strip with a hole, subjected to shear: homogenization error for two different lattices (columns), using the same conventions as in Figure 6.3 and 6.5. Same boundary conditions as in Figure 6.5.</figcaption></figure>

**6.3. Cracked strip in tension**

Before moving on to more complex examples, we consider the cracked version of the rectangular honeycomb strip shown in Figure 6.8, having an aspect-ratio $\delta/L = 0.78$. The crack is obtained by removing the beams that intersect the segment with length $a$, making an angle $\pi/6$ with the direction $e_1$, and whose midpoint is at the center of the lattice. Crack propagation is not considered. The same boundary conditions are used as for the shearing test, but a longitudinal displacement $v_{\beta \in rb} = \Delta e_1$ is imposed on the right-hand-side boundary.

<figure><img src="image_19.png" /><figcaption>**Figure 6.8.** Cracked honeycomb strip in tension.</figcaption></figure>

Colormaps of the homogenization error are shown in Figure 6.9 for $l = 0.0125$. Similar to what has been obtained for a punched lattice, some stress concentration is observed at both crack tips. When second-order homogenization is used, not only does the homogenization error drop significantly, it is also very small except in the intermediate
neighborhood of the crack tips and lips and of the edges of the rectangular region. The fact that the perturbation at the crack tips are disconnected from each other suggests that this geometry can be accurately represented by a higher-order continuum model endowed with effective boundary conditions for the crack tips, crack lips and free boundaries.

<figure><img src="image_18.png" /><figcaption>Figure 6.9. Cracked honeycomb strip in tension: homogenization errors for $\ell = 0.0125$, using the same conventions as in Figure 6.3. The representation of the boundary conditions is simplified, the reader is referred to Figure 6.8 for details.</figcaption></figure>

## 7. EXTENSIONS AND FURTHER ILLUSTRATIONS

This section illustrates (i) more advanced features of the homogenization method, namely non-uniform elastic properties (§7.1), non-uniform geometric properties (§7.3) and pre-strain (§7.2), as well as (ii) its versatility, by exploring different lattice topologies (2D Kagome lattice in §7.2), different types of elements (springs §7.3–7.5), different space dimensions (3D truss in §7.4, 1D truss in §7.5) as well the case of a large elastic contrast (§7.5).

### 7.1. Inhomogeneous elastic properties

We consider an extension of the honeycomb compression test done in Section 6.1, now using a lattice whose bending rigidity is graded in space: a beam $\alpha$ with midpoint position $X_{\alpha}^c$ is assigned a stiffness $EI(X_{\alpha}^c)$, where

$$EI(X_1, X_2) = \frac{1}{6} \left( 1 + 5 \tanh^2 \left( 6 \frac{2X_2^3 - X_1}{1 + 4X_2^2} \right) \right), \quad (7.1)$$

for $-1/2 \le X_1 \le +1/2$ and $-\delta/2 \le X_2 \le +\delta/2$, with the origin of the coordinate axes at the center of the lattice and $\delta = 0.84$. The distribution (7.1) of bending stiffness is shown in Figure 7.1. The lattice is six times softer along an S-shaped region that extends to opposite corners. As earlier, we work in the inextensible regime by using a large stretching modulus $EA = 10^{12}$ for all beams. The detailed boundary conditions are shown in Figure 7.1.

<figure><img src="image_19.png" /><figcaption>Figure 7.1. Honeycomb lattice with a graded bending modulus $EI(X)$ given by (7.1). The bending modulus $EI$ is constant in every beam and is found by evaluating the function $EI(X)$ at the beam midpoint $X_{\alpha}^c$. The boundary conditions are depicted using the same conventions as earlier.</figcaption></figure>

The homogenization is based on the inextensible beam model, as earlier. Compared to the homogeneous case (§5.4), the only change is that the symbolic parameter $EI$ must now be included in the list of variable material parameters,

$$m = \{EI\}. \quad (7.2)$$
neighborhood of the crack tips and lips and of the edges of the rectangular region. The fact that the perturbation at the crack tips are disconnected from each other suggests that this geometry can be accurately represented by a higher-order continuum model endowed with effective boundary conditions for the crack tips, crack lips and free boundaries.

<figure><img src="image_20.png" /><figcaption>Figure 6.9. Cracked honeycomb strip in tension: homogenization errors for $\ell = 0.0125$, using the same conventions as in Figure 6.3. The representation of the boundary conditions is simplified, the reader is referred to Figure 6.8 for details.</figcaption></figure>

# 7. EXTENSIONS AND FURTHER ILLUSTRATIONS

This section illustrates (i) more advanced features of the homogenization method, namely non-uniform elastic properties (§7.1), non-uniform geometric properties (§7.3) and pre-strain (§7.2), as well as (ii) its versatility, by exploring different lattice topologies (2D Kagome lattice in §7.2), different types of elements (springs §7.3–7.5), different space dimensions (3D truss in §7.4, 1D truss in §7.5) as well the case of a large elastic contrast (§7.5).

## 7.1. Inhomogeneous elastic properties

We consider an extension of the honeycomb compression test done in Section 6.1, now using a lattice whose bending rigidity is graded in space: a beam $\alpha$ with midpoint position $X_{\alpha}^c$ is assigned a stiffness $EI(X_{\alpha}^c)$, where

$$EI(X_1, X_2) = \frac{1}{6} \left( 1 + 5 \tanh^2 \left( 6 \frac{2X_2^3 - X_1}{1 + 4X_2^2} \right) \right), \quad (7.1)$$

for $-1/2 \le X_1 \le +1/2$ and $-\delta/2 \le X_2 \le +\delta/2$, with the origin of the coordinate axes at the center of the lattice and $\delta = 0.84$. The distribution (7.1) of bending stiffness is shown in Figure 7.1. The lattice is six times softer along an S-shaped region that extends to opposite corners. As earlier, we work in the inextensible regime by using a large stretching modulus $EA = 10^{12}$ for all beams. The detailed boundary conditions are shown in Figure 7.1.

<figure><img src="image_21.png" /><figcaption>Figure 7.1. Honeycomb lattice with a graded bending modulus $EI(X)$ given by (7.1). The bending modulus $EI$ is constant in every beam and is found by evaluating the function $EI(X)$ at the beam midpoint $X_{\alpha}^c$. The boundary conditions are depicted using the same conventions as earlier.</figcaption></figure>

The homogenization is based on the inextensible beam model, as earlier. Compared to the homogeneous case (§5.4), the only change is that the symbolic parameter $EI$ must now be included in the list of variable material parameters,

$$m = \{EI\}. \quad (7.2)$$
This is done by changing a single line of code: m={} is replaced with m={EI}, see the companion notebook homogenize-honeycomb-inextensible-variable nb. The library delivers the microscopic displacement in the form

$$
\begin{align}
\xi_b^* &= \ell \left( \frac{(-1)^{b+1}}{4} \mathcal{F} : \varepsilon \ell^0 + 0 \ell^1 + \mathcal{O}(\ell^2) \right), \tag{7.3} \\
\psi_b^* &= 0 \ell^0 + \frac{3(-1)^{b+1}}{4} \left( -\operatorname{curl} (\mathcal{F} : \varepsilon) + \frac{(\mathcal{F} : \varepsilon) \wedge \nabla (EI)}{EI} \right) \ell^1 + \mathcal{O}(\ell^2),
\end{align}
$$

where the constant tensor $\mathcal{F}$ has been defined in (5.5), the curl operator has been introduced in (5.6), and $\wedge$ denotes the 2D wedge product, $\mathbf{a} \wedge \mathbf{b} = a_1 b_2 - a_2 b_1$. The only difference with the homogeneous case (5.10) is the presence of the gradient term $\nabla(EI)$ in the microscopic rotation $\psi_b$. The gradient effect has now two contributions, one associated with the strain gradient $\nabla\varepsilon$ and the other one with the gradient of elastic properties $\nabla(EI)$.

**Remark 7.1.** The statement (7.2) informs the library that the parameter EI entering in the stiffness $\mathcal{H}$ varies in space, implying that the gradient must be calculated by the chain rule as $\nabla \mathcal{H} = \frac{d\mathcal{H}}{d(EI)} \cdot \nabla EI$. The library does not know about the specific profile $EI(X)$ in (7.1) and treats the gradient $\nabla(EI)$ symbolically: the explicit expression of $\nabla(EI)$ must be inserted by the user into the output (7.3) of the library.

The validity of the predictions (7.3) are tested using the same procedure as earlier. The results are shown in Figure 7.2. In the discrete simulations, the deformation is concentrated in the soft part of the lattice; by a Poisson’s effect, the lattice forms a neck in the center. The rapid variations in bending stiffness creates strong gradients that make the predictions of the leading-order homogenization relatively inaccurate, see Figure 7.2(b). Taking this gradient effect into account improves the accuracy by approximately one order of magnitude, see Figure 7.2(c). Both terms in $(7.3)_2$ are equally important: dropping either one deteriorates the accuracy significantly (plots not shown).

<figure><img src="image_20.png" /><figcaption>**Figure 7.2.** Honeycomb lattice with graded bending modulus: accuracy of the prediction (7.3) for the microscopic displacement. All beams have length $\ell = 0.0125$.</figcaption></figure>

## 7.2. Kagome lattice with pre-strain

In this section, we change the lattice topology to *Kagome* and analyze the effect of a non-uniform pre-strain, as produced typically by thermal expansion in the presence of a localized heat source. The rectangular lattice is shown in Figure 7.3a: it is pinned at two adjacent corners and there are no other loading than the pre-strain. It is made up of beams with length $\ell = 0.006579$, aspect-ratio parameter $\chi = 0.002$, stretching modulus $EA = 1$ and bending modulus $EI = EA \ell^2 \chi$.

In the discrete simulations, the extensional pre-strain $p(X)$ is accounted for by changing the strain energy (2.20) of the beams to

$$ w_{\alpha}(\mathbf{E}_{\alpha}) = \frac{1}{2}(\mathbf{E}_{\alpha} - \mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c}))) \cdot \mathcal{H}_{\phi_{\alpha}}(\mathbf{E}_{\alpha} - \mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c}))) \quad \text{where} \quad \mathbf{E}_{0}(p) = (p, 0, 0). \qquad (7.4) $$

The pre-strain vector $\mathbf{E}_0$ reflects the ordering convention in (2.12): the pre-strain $p$ being extensional, it offsets the first slot of $\mathbf{E}_\alpha$ representing the stretching strain $\varepsilon_\alpha$. In (7.4), the value of the pre-strain $p$ in any particular beam $\alpha$ is picked based on a pre-strain map $p(X)$ evaluated at the midpoint $\mathbf{X}_\alpha^c$: in the simulations, we use the map

$$ p(\mathbf{X}) = \Delta \exp \left[ - \left( \frac{\|\mathbf{X}\|}{0.2} \right)^2 \right] \qquad (7.5) $$

plotted in Figure 7.3a, which has a maximum at the center $\mathbf{X}=\mathbf{0}$ of the upper boundary, and decreases exponentially away from this point. The pre-strain magnitude is set to $\Delta = 0.3$ in the simulations but this value is a matter of convention as all our results are rescaled by $\Delta$. The results of the discrete simulations are shown in Figure 2.20 (deformed lattice) along with the homogenization results (colors).
This is done by changing a single line of code: m={} is replaced with m={EI}, see the companion notebook homogenize-honeycomb-inextensible-variable nb. The library delivers the microscopic displacement in the form

$$
\begin{align}
\xi_b^* &= \ell \left( \frac{(-1)^{b+1}}{4} \mathcal{F} : \varepsilon \ell^0 + 0 \ell^1 + \mathcal{O}(\ell^2) \right), \tag{7.3} \\
\psi_b^* &= 0 \ell^0 + \frac{3(-1)^{b+1}}{4} \left( -\operatorname{curl} (\mathcal{F} : \varepsilon) + \frac{(\mathcal{F} : \varepsilon) \wedge \nabla (EI)}{EI} \right) \ell^1 + \mathcal{O}(\ell^2),
\end{align}
$$

where the constant tensor $\mathcal{F}$ has been defined in (5.5), the curl operator has been introduced in (5.6), and $\wedge$ denotes the 2D wedge product, $\mathbf{a} \wedge \mathbf{b} = a_1 b_2 - a_2 b_1$. The only difference with the homogeneous case (5.10) is the presence of the gradient term $\nabla(EI)$ in the microscopic rotation $\psi_b$. The gradient effect has now two contributions, one associated with the strain gradient $\nabla\varepsilon$ and the other one with the gradient of elastic properties $\nabla(EI)$.

**Remark 7.1.** The statement (7.2) informs the library that the parameter EI entering in the stiffness $\mathcal{H}$ varies in space, implying that the gradient must be calculated by the chain rule as $\nabla \mathcal{H} = \frac{d\mathcal{H}}{d(EI)} \cdot \nabla EI$. The library does not know about the specific profile $EI(X)$ in (7.1) and treats the gradient $\nabla(EI)$ symbolically: the explicit expression of $\nabla(EI)$ must be inserted by the user into the output (7.3) of the library.

The validity of the predictions (7.3) are tested using the same procedure as earlier. The results are shown in Figure 7.2. In the discrete simulations, the deformation is concentrated in the soft part of the lattice; by a Poisson’s effect, the lattice forms a neck in the center. The rapid variations in bending stiffness creates strong gradients that make the predictions of the leading-order homogenization relatively inaccurate, see Figure 7.2(b). Taking this gradient effect into account improves the accuracy by approximately one order of magnitude, see Figure 7.2(c). Both terms in $(7.3)_2$ are equally important: dropping either one deteriorates the accuracy significantly (plots not shown).

<figure><img src="image_22.png" /><figcaption>**Figure 7.2.** Honeycomb lattice with graded bending modulus: accuracy of the prediction (7.3) for the microscopic displacement. All beams have length $\ell = 0.0125$.</figcaption></figure>

## 7.2. Kagome lattice with pre-strain

In this section, we change the lattice topology to *Kagome* and analyze the effect of a non-uniform pre-strain, as produced typically by thermal expansion in the presence of a localized heat source. The rectangular lattice is shown in Figure 7.3a: it is pinned at two adjacent corners and there are no other loading than the pre-strain. It is made up of beams with length $\ell = 0.006579$, aspect-ratio parameter $\chi = 0.002$, stretching modulus $EA = 1$ and bending modulus $EI = EA \ell^2 \chi$.

In the discrete simulations, the extensional pre-strain $p(X)$ is accounted for by changing the strain energy (2.20) of the beams to

$$ w_{\alpha}(\mathbf{E}_{\alpha}) = \frac{1}{2}(\mathbf{E}_{\alpha} - \mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c}))) \cdot \mathcal{H}_{\varphi_{\alpha}}(\mathbf{E}_{\alpha} - \mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c}))) \quad \text{where} \quad \mathbf{E}_{0}(p) = (p, 0, 0). \qquad (7.4) $$

The pre-strain vector $\mathbf{E}_0$ reflects the ordering convention in (2.12): the pre-strain $p$ being extensional, it offsets the first slot of $\mathbf{E}_\alpha$ representing the stretching strain $\varepsilon_\alpha$. In (7.4), the value of the pre-strain $p$ in any particular beam $\alpha$ is picked based on a pre-strain map $p(X)$ evaluated at the midpoint $\mathbf{X}_\alpha^c$: in the simulations, we use the map

$$ p(\mathbf{X}) = \Delta \exp \left[ - \left( \frac{\|\mathbf{X}\|}{0.2} \right)^2 \right] \qquad (7.5) $$

plotted in Figure 7.3a, which has a maximum at the center $\mathbf{X}=\mathbf{0}$ of the upper boundary, and decreases exponentially away from this point. The pre-strain magnitude is set to $\Delta = 0.3$ in the simulations but this value is a matter of convention as all our results are rescaled by $\Delta$. The results of the discrete simulations are shown in Figure 2.20 (deformed lattice) along with the homogenization results (colors).
<figure><img src="image_21.png" /><figcaption>**Figure 7.3.** Kagome lattice with pre-strain. (a) Imposed pre-strain map $p(\mathbf{X})$ and boundary conditions. (b) Underlying topological lattice, made up of $n_\varphi = 6$ beam families and $n_b = 3$ Bravais sub-lattices. In each family, a particular beam is shown using a thick line, with the arrow denoting the conventional orientation. All beams have equal length $\ell = 0.006579$, stretching modulus $EA = 1$ and aspect-ratio parameter $\chi = 0.002$.</figcaption></figure>

We proceed to homogenize the lattice. The underlying topological Kagome lattice makes use of the conventions shown in Figure 7.3b. It is specified in the homogenization code based on the properties listed in Table 7.1.

<table><caption>**Table 7.1.** Properties of the topological Kagome lattice shown in Figure 7.3b, used as an input to the homogenization code: space dimension *d*, number of Bravais sub-lattices $n_b$, number of element families $n_\varphi$, density $\bar{\rho}_\varphi$ per unit area of element belonging to any particular family $\varphi$, connectivity matrix $b_\varphi^\pm$ and end-to-end vector $\delta_\varphi = \Delta_\varphi^+ - \Delta_\varphi^-$ for a beam of type $\varphi$ in reference configuration.</caption><thead><tr><th></th><th><i>d</i></th><th><i>n</i><sub>b</sub></th><th><i>n</i><sub>φ</sub></th><th><i>ρ</i><sub>φ</sub></th><th><i>b</i><sub>φ</sub><sup>-</sup></th><th><i>b</i><sub>φ</sub><sup>+</sup></th><th><i>δ</i><sub>φ</sub></th></tr></thead><tbody><tr><td></td><td>2</td><td>3</td><td>6</td><td><sup>1</sup>&frasl;<sub>2&radic;3 <i>ℓ</i><sup>2</sup></sub></td><td>1</td><td>1</td><td><i>ℓ</i><b>e</b><sub>1</sub> <i>ℓ</i><span style="font-size: 1.5em; vertical-align: middle;">(</span><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;"><i>e</i><sub>1</sub></div><div style="padding: 0 5px;">2</div></div><span style="font-size: 1.5em; vertical-align: middle;">-</span><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">&radic;3</div><div style="padding: 0 5px;">2</div></div><b>e</b><sub>2</sub><span style="font-size: 1.5em; vertical-align: middle;">)</span></td><td>2</td><td>3</td><td>2</td><td>1</td><td>3</td><td>1</td></tr></tbody></table>

In the homogenization library, the pre-strain is taken into account by expanding (7.4) as

$$w_{\alpha}(\mathbf{E}_{\alpha}) = \frac{1}{2} \mathbf{E}_{\alpha} \cdot [\mathcal{H}_{\varphi_{\alpha}} \cdot \mathbf{E}_{\alpha}] + 1 [-\mathcal{H}_{\varphi_{\alpha}} \cdot \mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c}))] \cdot \mathbf{E}_{\alpha} + \frac{1}{2} [\mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c})) \cdot \mathcal{H}_{\varphi_{\alpha}} \cdot \mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c}))] \mathbf{1}. \quad (7.6)$$

Contrary to what was assumed earlier in (3.29), the element energy (7.6) is no longer homogeneous of degree 2 with respect to the unknown strain $\mathbf{E}_\alpha$. To work around this, we incorporate the 1's appearing in (7.6) into the macroscopic strain $\mathbf{E}$, amending its definition (3.25) as—mind the trailing 1,

$$\mathbf{E} = \left( \varepsilon_I, \kappa_I, \tau_I, \dots, \varepsilon_{n_\varphi}, \kappa_{n_\varphi}, \tau_{n_\varphi}, \eta^{-1} \sum_{b=1}^{n_b} (\xi_b)_1, \dots, \eta^{-1} \sum_{b=1}^{n_b} (\xi_b)_d, 1 \right). \quad (7.7)$$

Equation (7.6) can then be implemented by incorporating the square brackets appearing in (7.6) into the element contribution $\mathcal{H}_\varphi^g$ (3.31), which we now define as

$$\mathcal{H}_{\varphi}^{g}(\eta, \mathbf{X}) = \begin{pmatrix} \mathbf{0}_{p_{\varphi} \times p_{\varphi}} & \mathbf{0}_{p_{\varphi} \times n_{E_{\varphi}}} & \mathbf{0}_{p_{\varphi} \times (p_{\varphi}' - 1)} & \mathbf{0}_{p_{\varphi} \times 1} \\ \mathbf{0}_{n_{E_{\varphi}} \times p_{\varphi}} & [\mathcal{H}_{\varphi}] & \mathbf{0}_{n_{E_{\varphi}} \times (p_{\varphi}' - 1)} & [-\mathcal{H}_{\varphi} \cdot \mathbf{E}_0(p(\mathbf{X}))] \\ \mathbf{0}_{(p_{\varphi}' - 1) \times p_{\varphi}} & \mathbf{0}_{(p_{\varphi}' - 1) \times n_{E_{\varphi}}} & \mathbf{0}_{(p_{\varphi}' - 1) \times (p_{\varphi}' - 1)} & \mathbf{0}_{(p_{\varphi}' - 1) \times 1} \\ \mathbf{0}_{1 \times p_{\varphi}} & [-\mathcal{H}_{\varphi} \cdot \mathbf{E}_0(p(\mathbf{X}))]^T & \mathbf{0}_{1 \times (p_{\varphi}' - 1)} & [\mathbf{E}_0(p(\mathbf{X})) \cdot \mathcal{H}_{\varphi} \cdot \mathbf{E}_0(p(\mathbf{X}))] \end{pmatrix}. \quad (7.8)$$

Now, since $\mathbf{E}$ is required to depend linearly on $l$, $y$ and their gradients by (3.34), the updated definition of $\mathbf{E}$ in (7.7) requires that we append a trailing entry 1 in the macroscopic strain $l$ as well, modifying (3.24) into

$$l = (\tilde{\varepsilon}_1 = \varepsilon_{11}, \tilde{\varepsilon}_2 = \varepsilon_{22}, \tilde{\varepsilon}_3 = \sqrt{2} \varepsilon_{12}, 1), \quad (7.9)$$

and that we propagate it to $\mathbf{E}$: to do so, we extend the dimension of the tensors $\mathbf{E}_l, \dots, \mathbf{E}_y''$ appropriately and fill the new entries with 0's, except that in the lower-right corner of $\mathbf{E}_l$ which is set to 1. In the library, the content of the vectors $l, y, \mathbf{E}$, etc. is not fixed once for all but can instead be chosen in a flexible way, as documented in the Implementation notes that are distributed along with the library: this makes the implementation of the above changes straightforward.

From the user's perspective, the pre-strain is dealt with by including an option `includeUnitStrain=True` that sets up the trailing 1's in both $\mathbf{E}$ and $l$, and by providing an optional argument `prestrain={"$\epsilon$":p}` to the constructor for beam elements. The pre-strain is represented by a pure symbol $p$ that is treated in a similar way as the
<figure><img src="image_23.png" /><figcaption>**Figure 7.3.** Kagome lattice with pre-strain. (a) Imposed pre-strain map $p(\mathbf{X})$ and boundary conditions. (b) Underlying topological lattice, made up of $n_\varphi = 6$ beam families and $n_b = 3$ Bravais sub-lattices. In each family, a particular beam is shown using a thick line, with the arrow denoting the conventional orientation. All beams have equal length $\ell = 0.006579$, stretching modulus $EA = 1$ and aspect-ratio parameter $\chi = 0.002$.</figcaption></figure>

We proceed to homogenize the lattice. The underlying topological Kagome lattice makes use of the conventions shown in Figure 7.3b. It is specified in the homogenization code based on the properties listed in Table 7.1.

<table><caption>**Table 7.1.** Properties of the topological Kagome lattice shown in Figure 7.3b, used as an input to the homogenization code: space dimension $d$, number of Bravais sub-lattices $n_b$, number of element families $n_\varphi$, density $\tilde{\rho}_\varphi$ per unit area of element belonging to any particular family $\varphi$, connectivity matrix $b_\varphi^\pm$ and end-to-end vector $\delta_\varphi = \Delta_\varphi^+ - \Delta_\varphi^-$ for a beam of type $\varphi$ in reference configuration.</caption><thead><tr><th></th><th>d</th><th>n<sub>b</sub></th><th>n<sub>φ</sub></th><th>ρ<sub>φ</sub></th><th>φ</th><th>I</th><th>II</th><th>III</th><th>IV</th><th>V</th><th>VI</th></tr></thead><tbody><tr><td></td><td>2</td><td>3</td><td>6</td><td>1 / (2√3 ℓ<sup>2</sup>)</td><td>b<sub>φ</sub><sup>-</sup></td><td>1</td><td>1</td><td>3</td><td>2</td><td>2</td><td>3</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>b<sub>φ</sub><sup>+</sup></td><td>2</td><td>3</td><td>2</td><td>1</td><td>3</td><td>1</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>δ<sub>φ</sub></td><td>ℓe<sub>1</sub></td><td>ℓ(e<sub>1</sub>/2 - √3/2e<sub>2</sub>)</td><td>ℓ(e<sub>1</sub>/2 + √3/2e<sub>2</sub>)</td><td>ℓe<sub>1</sub></td><td>ℓ(e<sub>1</sub>/2 + √3/2e<sub>2</sub>)</td><td>ℓ(e<sub>1</sub>/2 - √3/2e<sub>2</sub>)</td></tr></tbody></table>

In the homogenization library, the pre-strain is taken into account by expanding (7.4) as

$$w_{\alpha}(\mathbf{E}_{\alpha}) = \frac{1}{2} \mathbf{E}_{\alpha} \cdot [\mathcal{H}_{\varphi_{\alpha}} \cdot \mathbf{E}_{\alpha}] + 1 [-\mathcal{H}_{\varphi_{\alpha}} \cdot \mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c}))] \cdot \mathbf{E}_{\alpha} + \frac{1}{2} [\mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c})) \cdot \mathcal{H}_{\varphi_{\alpha}} \cdot \mathbf{E}_{0}(p(\mathbf{X}_{\alpha}^{c}))] \cdot \mathbf{E}_{\alpha} \quad (7.6)$$

Contrary to what was assumed earlier in (3.29), the element energy (7.6) is no longer homogeneous of degree 2 with respect to the unknown strain $\mathbf{E}_\alpha$. To work around this, we incorporate the 1's appearing in (7.6) into the macroscopic strain $\mathbf{E}$, amending its definition (3.25) as—mind the trailing 1,

$$\mathbf{E} = \left( \varepsilon_I, \kappa_I, \tau_I, \dots, \varepsilon_{n_\varphi}, \kappa_{n_\varphi}, \tau_{n_\varphi}, \eta^{-1} \sum_{b=1}^{n_b} (\xi_b)_1, \dots, \eta^{-1} \sum_{b=1}^{n_b} (\xi_b)_d, 1 \right). \quad (7.7)$$

Equation (7.6) can then be implemented by incorporating the square brackets appearing in (7.6) into the element contribution $\mathcal{H}_\varphi^g$ (3.31), which we now define as

$$\mathcal{H}_{\varphi}^{g}(\eta, \mathbf{X}) = \begin{pmatrix} \mathbf{0}_{p_{\varphi} \times p_{\varphi}} & \mathbf{0}_{p_{\varphi} \times n_{E_{\varphi}}} & \mathbf{0}_{p_{\varphi} \times (p'_{\varphi}-1)} & \mathbf{0}_{p_{\varphi} \times 1} \\ \mathbf{0}_{n_{E_{\varphi}} \times p_{\varphi}} & [\mathcal{H}_{\varphi}] & \mathbf{0}_{n_{E_{\varphi}} \times (p'_{\varphi}-1)} & [-\mathcal{H}_{\varphi} \cdot \mathbf{E}_0(p(\mathbf{X}))] \\ \mathbf{0}_{(p'_{\varphi}-1) \times p_{\varphi}} & \mathbf{0}_{(p'_{\varphi}-1) \times n_{E_{\varphi}}} & \mathbf{0}_{(p'_{\varphi}-1) \times (p'_{\varphi}-1)} & \mathbf{0}_{(p'_{\varphi}-1) \times 1} \\ \mathbf{0}_{1 \times p_{\varphi}} & [-\mathcal{H}_{\varphi} \cdot \mathbf{E}_0(p(\mathbf{X}))]^T & \mathbf{0}_{1 \times (p'_{\varphi}-1)} & [\mathbf{E}_0(p(\mathbf{X})) \cdot \mathcal{H}_{\varphi} \cdot \mathbf{E}_0(p(\mathbf{X}))] \end{pmatrix}. \quad (7.8)$$

Now, since $\mathbf{E}$ is required to depend linearly on $l$, $y$ and their gradients by (3.34), the updated definition of $\mathbf{E}$ in (7.7) requires that we append a trailing entry 1 in the macroscopic strain $l$ as well, modifying (3.24) into

$$l = (\tilde{\varepsilon}_1 = \varepsilon_{11}, \tilde{\varepsilon}_2 = \varepsilon_{22}, \tilde{\varepsilon}_3 = \sqrt{2} \varepsilon_{12}, 1), \quad (7.9)$$

and that we propagate it to $\mathbf{E}$: to do so, we extend the dimension of the tensors $\mathbf{E}_l, \dots, \mathbf{E}_y''$ appropriately and fill the new entries with 0's, except that in the lower-right corner of $\mathbf{E}_l$ which is set to 1. In the library, the content of the vectors $l, y, \mathbf{E}$, etc. is not fixed once for all but can instead be chosen in a flexible way, as documented in the Implementation notes that are distributed along with the library: this makes the implementation of the above changes straightforward.

From the user's perspective, the pre-strain is dealt with by including an option `includeUnitStrain=True` that sets up the trailing 1's in both $\mathbf{E}$ and $l$, and by providing an optional argument `prestrain={"$\epsilon$":p}` to the constructor for beam elements. The pre-strain is represented by a pure symbol $p$ that is treated in a similar way as the
variable stiffness in section 7.1: it is declared as a *variable material parameter* by the assignment $m=\{p\}$, and the homogenization results are given symbolically in terms of $\mathbf{m}=(p)$ and its gradient $\nabla \mathbf{m}=(\nabla p)$, as well as on $\mathbf{l}$ and its gradients $\nabla \mathbf{l}$, etc. To interpret these results, the user must take into account the fact that $(\nabla l)_{4i} = \partial(l_4)/\partial X_i = \partial 1/\partial X_i = 0$ in view of (7.9).

In the limit $EI \to 0$ of perfectly flexible beams, the Kagome lattice possesses a number of zero-energy modes [HF06]. One of these zero modes has infinite wavelength and therefore survives homogenization: it involves a rigid rotation of adjacent triangles in opposite directions and translates the Bravais sub-lattice having index $b$ by a vector $\mathcal{W}_b$,

$$ \mathcal{W}_1 = \frac{-e_1 - \sqrt{3}e_2}{2}, \quad \mathcal{W}_2 = \frac{-e_1 + \sqrt{3}e_2}{2}, \quad \mathcal{W}_3 = e_1. \qquad (7.10) $$

Indeed, the corresponding microscopic displacement $\mathbf{y}=(\mathcal{W}_1 0 \mathcal{W}_2 0 \mathcal{W}_3 0)$ makes all the extensional degrees of freedom in $E_l \cdot \mathbf{y}$ vanish.

The Kagome lattice is homogenized in the Mathematica notebook `homogenize-kagome-prestress nb` distributed along with the library. The lattice made up of extensible beams, corresponding to the element stiffness matrix in (2.19). The microscopic displacement computed by the library takes the form

$$ \begin{gathered} \xi_b^* = \ell \begin{pmatrix} 0 & \ell^0 + \frac{c_1(\chi)}{\chi} \operatorname{div}(\mathcal{F}:\boldsymbol{\varepsilon}) \mathcal{W}_b \\ + [\mathcal{W}_b \cdot \mathcal{F}] \cdot (c_2(\chi)(\nabla \operatorname{tr} \boldsymbol{\varepsilon} - 2 \nabla p) + c_3(\chi) \mathcal{F}:\nabla(\mathcal{F}:\boldsymbol{\varepsilon})) \\ + c_4(\chi) \operatorname{curl}(\mathcal{F}:\boldsymbol{\varepsilon}) [\mathcal{A}^T \cdot \mathcal{W}_b] \end{pmatrix} \ell + \mathcal{O}(\ell^2) \\ \psi_b^* = -\frac{1}{2} \mathcal{W}_b \wedge (\mathcal{F}:\boldsymbol{\varepsilon}) \ell^0 + 0 \ell^1 + \mathcal{O}(\ell^2), \end{gathered} \qquad (7.11) $$

where $\mathcal{A}$ is the Levi-Civita (purely anti-symmetric) symbol in dimension 2, $\mathcal{A}=e_1 \otimes e_2 - e_2 \otimes e_1$ and the coefficients $c_i$ are rational functions of $\chi$ that are bounded for $\chi \to 0$,

$$ \begin{align} c_1(\chi) &= \frac{1+18\chi}{288} & c_2(\chi) &= \frac{1}{12(1+12\chi)} \tag{7.12} \\ c_3(\chi) &= \frac{1-24\chi}{24(1+12\chi)} & c_4(\chi) &= -\chi/4. \notag \end{align} $$

The compact tensorial expressions in (7.11) has been obtained by applying to the raw output of the homogenization code a procedure similar to that yielding the irreducible form of elasticity tensors in the presence of symmetries [AKO17].

In Figure 7.4, the expression (7.11) of the microscopic displacement is verified using the usual procedure. The leading-order homogenization fails to give accurate predictions in regions where the gradient of pre-strain is important and the improvement brought about by second-order homogenization is very significant there. If the correction proportional to $\nabla p$ is dropped in the second-order homogenization results (7.11)$_2$, the agreement with the numerical solution degrades noticeably (data not shown).

<figure><img src="image_22.png" /><figcaption>**Figure 7.4.** Kagome lattice with spatially graded pre-strain: verification of the microscopic displacement (7.11) based on the discrete numerical solution. Same lattice as in Figure 7.3.</figcaption></figure>
variable stiffness in section 7.1: it is declared as a *variable* material parameter by the assignment $m=\{p\}$, and the homogenization results are given symbolically in terms of $\mathbf{m}=(p)$ and its gradient $\nabla \mathbf{m}=(\nabla p)$, as well as on $\mathbf{l}$ and its gradients $\nabla \mathbf{l}$, etc. To interpret these results, the user must take into account the fact that $(\nabla l)_{4i} = \partial(l_4)/\partial X_i = \partial 1/\partial X_i = 0$ in view of (7.9).

In the limit $EI \to 0$ of perfectly flexible beams, the Kagome lattice possesses a number of zero-energy modes [HF06]. One of these zero modes has infinite wavelength and therefore survives homogenization: it involves a rigid rotation of adjacent triangles in opposite directions and translates the Bravais sub-lattice having index $b$ by a vector $\mathcal{W}_b$,

$$ \mathcal{W}_1 = \frac{-e_1 - \sqrt{3}e_2}{2}, \quad \mathcal{W}_2 = \frac{-e_1 + \sqrt{3}e_2}{2}, \quad \mathcal{W}_3 = e_1. \qquad (7.10) $$

Indeed, the corresponding microscopic displacement $\mathbf{y} = (\mathcal{W}_1 0 \mathcal{W}_2 0 \mathcal{W}_3 0)$ makes all the extensional degrees of freedom in $E_l \cdot \mathbf{y}$ vanish.

The Kagome lattice is homogenized in the Mathematica notebook homogenize-kagome-prestress nb distributed along with the library. The lattice made up of extensible beams, corresponding to the element stiffness matrix in (2.19). The microscopic displacement computed by the library takes the form

$$ \begin{gathered} \xi_b^* = \ell \begin{pmatrix} 0 & \ell^0 + \frac{c_1(\chi)}{\chi} \operatorname{div}(\mathcal{F}:\boldsymbol{\varepsilon}) \mathcal{W}_b \\ + [\mathcal{W}_b \cdot \mathcal{F}] \cdot (c_2(\chi)(\nabla \operatorname{tr} \boldsymbol{\varepsilon} - 2 \nabla p) + c_3(\chi) \mathcal{F}:\nabla(\mathcal{F}:\boldsymbol{\varepsilon})) \\ + c_4(\chi) \operatorname{curl}(\mathcal{F}:\boldsymbol{\varepsilon}) [\mathcal{A}^T \cdot \mathcal{W}_b] \end{pmatrix} \ell + \mathcal{O}(\ell^2) \\ \psi_b^* = -\frac{1}{2} \mathcal{W}_b \wedge (\mathcal{F}:\boldsymbol{\varepsilon}) \ell^0 + 0 \ell^1 + \mathcal{O}(\ell^2), \end{gathered} \qquad (7.11) $$

where $\mathcal{A}$ is the Levi-Civita (purely anti-symmetric) symbol in dimension 2, $\mathcal{A} = e_1 \otimes e_2 - e_2 \otimes e_1$ and the coefficients $c_i$ are rational functions of $\chi$ that are bounded for $\chi \to 0$,

$$ \begin{align} c_1(\chi) &= \frac{1+18\chi}{288} & c_2(\chi) &= \frac{1}{12(1+12\chi)} \tag{7.12} \\ c_3(\chi) &= \frac{1-24\chi}{24(1+12\chi)} & c_4(\chi) &= -\chi/4. \notag \end{align} $$

The compact tensorial expressions in (7.11) has been obtained by applying to the raw output of the homogenization code a procedure similar to that yielding the irreducible form of elasticity tensors in the presence of symmetries [AKO17].

In Figure 7.4, the expression (7.11) of the microscopic displacement is verified using the usual procedure. The leading-order homogenization fails to give accurate predictions in regions where the gradient of pre-strain is important and the improvement brought about by second-order homogenization is very significant there. If the correction proportional to $\nabla p$ is dropped in the second-order homogenization results (7.11)$_2$, the agreement with the numerical solution degrades noticeably (data not shown).

<figure><img src="image_24.png" /><figcaption>**Figure 7.4.** Kagome lattice with spatially graded pre-strain: verification of the microscopic displacement (7.11) based on the discrete numerical solution. Same lattice as in Figure 7.3.</figcaption></figure>
**Remark 7.2.** We have considered a Kagome lattice made up of *extensible beams*. Indeed, it is neither possible to use *springs EI = 0* as the lattice then possess a zero-energy mode—the library aborts, printing out the zero mode (7.10)—, nor *inextensible beams* ($\chi = 0$) as this makes the lattice fully rigid—the library aborts, reporting a compatibility condition in the form $\epsilon(\mathbf{X}) = 0$ for all $\mathbf{X}$.

**Remark 7.3.** The $\chi^{-1}$ coefficient in factor of $c_1$ in the right-hand side of (7.11)$_1$ implies that the amplitude of the soft mode $W_b$ becomes larger and larger in the limit of inextensible beams, $\chi \to 0$. This is not surprising as, by definition, this mode has a vanishing elastic energy in the limit. A detailed analysis of the inextensible limit of the Kagome lattice is beyond the scope of this paper—in principle, a limit energy can be obtained by analyzing the energy of the discrete lattice computed by the library in the limit $\chi \to 0$, similar to what we did earlier in Section 5.3.

## 7.3. Circular arch

Our next example is a half-annular lattice made up of elastic springs (EI = 0), sketched in Figure 7.5. This lattice features *spatial inhomogeneity*: both the orientation and the natural length of the springs vary across the annular domain.

The discrete arch is generated as follows. We pick two integers $n_1$ and $n_2$ and define

$$ \eta = \frac{\pi}{n_2} \quad \text{and} \quad \delta = \pi \frac{n_1}{n_2}. \qquad (7.13) $$

In the simulations shown in Figure 7.7, we used $(n_1, n_2) = (13, 90)$. The nodes are indexed by a pair of integers $(i, j)$ with $-n_1 \le i \le +n_1$ and $0 \le j \le n_2$ and their position $\mathbf{X}_{(i,j)}$ in reference configuration is set to

$$ \mathbf{X}_{(i,j)} = \exp(\eta i) \mathbf{e}_r(\eta j) \quad \text{where} \quad \mathbf{e}_r(\theta) = \cos \theta \mathbf{e}_1 + \sin \theta \mathbf{e}_2. \qquad (7.14) $$

The inner and outer radii of the arch are therefore given as $\exp(\pm \eta n_1) = \exp(\pm \delta)$. Each node is assigned a Bravais index $b$ with $b=1$ if $i+j$ is even (black nodes in the figure) and $b=2$ if it is odd (white nodes). As shown in the figure, we set up elastic springs radially between all pairs of adjacent nodes, azimuthally between all pairs of adjacent nodes, and diagonally between pairs of adjacent black nodes ($b=1$).

<figure><img src="image_23.png" /><figcaption>**Figure 7.5.** Circular arch. This annular lattice is made up of spring elements connecting nodes $\mathbf{X}_{(i,j)}$ indexed by two integers, a radial coordinate $-n_1 \le i \le +n_1$ and an azimuthal coordinate $0 \le j \le n_2$. The short edge $j=n_2$ is blocked and the other edges are free. A weight-like vertical force $\mathbf{g}_{(i,j)}$ is applied on the nodes, with uniform density $\Delta$ per unit area.</figcaption></figure>

The elasticity of the spring (elements) is set up as follows. Considering a spring $\alpha$ having adjacent nodes $(i\bar{a}, j\bar{a})$ and $(i\bar{a}^+, j\bar{a}^+)$, we define the undeformed end-to-end vector as $\delta_\alpha = \mathbf{X}_{(i\bar{a}, j\bar{a})} - \mathbf{X}_{(i\bar{a}^+, j\bar{a}^+)}$, the natural spring length as $\ell_\alpha = \|\delta_\alpha\|$, the undeformed unit tangent $\mathbf{t}_\alpha = \delta_\alpha / \ell_\alpha$ and the linearized spring elongation $z_\alpha$ as

$$ z_\alpha = \mathbf{t}_\alpha \cdot (\mathbf{v}_{(i\bar{a}^+, j\bar{a}^+)} - \mathbf{v}_{(i\bar{a}, j\bar{a})}). \qquad (7.15) $$

The spring is assigned an elastic energy $w_\alpha = \frac{1}{2} k_\alpha z_\alpha^2$, where $k_\alpha$ is the spring constant which we set up as

$$ k_\alpha = \frac{EA}{\ell_\alpha} \qquad (7.16) $$

to represent elastic bars having all the same traction modulus EA. The short edge $j=n_2$ is blocked and a force

$$ \mathbf{g}_{(i,j)} = -\Delta (\eta r_i)^2 \mathbf{e}_2 \qquad (7.17) $$

is applied on the node $(i, j)$, representing a weight-like force having a constant density $(-\Delta \mathbf{e}_2)$ per unit area, where $r_i = \exp(\eta i)$ is the radial coordinate of the node. The nodal displacements $v_\beta$ at equilibrium are proportional to $\Delta$ and independent of the global elasticity constant EA. The invariable forces $\mathbf{g}_{(i,j)}$ are taken into account through their potential energy $-\sum_{i=-n_1}^{+n_1} \sum_{j=0}^{n_2} \mathbf{g}_{(i,j)} \cdot \mathbf{v}_{(i,j)}$.
**Remark 7.2.** We have considered a Kagome lattice made up of *extensible beams*. Indeed, it is neither possible to use *springs EI = 0* as the lattice then possess a zero-energy mode—the library aborts, printing out the zero mode (7.10)—, nor *inextensible beams* ($\chi = 0$) as this makes the lattice fully rigid—the library aborts, reporting a compatibility condition in the form $\epsilon(\mathbf{X}) = 0$ for all $\mathbf{X}$.

**Remark 7.3.** The $\chi^{-1}$ coefficient in factor of $c_1$ in the right-hand side of (7.11)$_1$ implies that the amplitude of the soft mode $W_b$ becomes larger and larger in the limit of inextensible beams, $\chi \to 0$. This is not surprising as, by definition, this mode has a vanishing elastic energy in the limit. A detailed analysis of the inextensible limit of the Kagome lattice is beyond the scope of this paper—in principle, a limit energy can be obtained by analyzing the energy of the discrete lattice computed by the library in the limit $\chi \to 0$, similar to what we did earlier in Section 5.3.

## 7.3. Circular arch

Our next example is a half-annular lattice made up of elastic springs (EI = 0), sketched in Figure 7.5. This lattice features *spatial inhomogeneity*: both the orientation and the natural length of the springs vary across the annular domain.

The discrete arch is generated as follows. We pick two integers $n_1$ and $n_2$ and define

$$ \eta = \frac{\pi}{n_2} \quad \text{and} \quad \delta = \pi \frac{n_1}{n_2}. \qquad (7.13) $$

In the simulations shown in Figure 7.7, we used $(n_1, n_2) = (13, 90)$. The nodes are indexed by a pair of integers $(i, j)$ with $-n_1 \le i \le +n_1$ and $0 \le j \le n_2$ and their position $\mathbf{X}_{(i,j)}$ in reference configuration is set to

$$ \mathbf{X}_{(i,j)} = \exp(\eta i) \mathbf{e}_r(\eta j) \quad \text{where} \quad \mathbf{e}_r(\theta) = \cos \theta \mathbf{e}_1 + \sin \theta \mathbf{e}_2. \qquad (7.14) $$

The inner and outer radii of the arch are therefore given as $\exp(\pm \eta n_1) = \exp(\pm \delta)$. Each node is assigned a Bravais index $b$ with $b=1$ if $i+j$ is even (black nodes in the figure) and $b=2$ if it is odd (white nodes). As shown in the figure, we set up elastic springs radially between all pairs of adjacent nodes, azimuthally between all pairs of adjacent nodes, and diagonally between pairs of adjacent black nodes ($b=1$).

<figure><img src="image_25.png" /><figcaption>**Figure 7.5.** Circular arch. This annular lattice is made up of spring elements connecting nodes $\mathbf{X}_{(i,j)}$ indexed by two integers, a radial coordinate $-n_1 \le i \le +n_1$ and an azimuthal coordinate $0 \le j \le n_2$. The short edge $j=n_2$ is blocked and the other edges are free. A weight-like vertical force $\mathbf{g}_{(i,j)}$ is applied on the nodes, with uniform density $\Delta$ per unit area.</figcaption></figure>

The elasticity of the spring (elements) is set up as follows. Considering a spring $\alpha$ having adjacent nodes $(i\bar{a}, j\bar{a})$ and $(i\bar{a}^+, j\bar{a}^+)$, we define the undeformed end-to-end vector as $\delta_\alpha = \mathbf{X}_{(i\bar{a}^+, j\bar{a}^+)} - \mathbf{X}_{(i\bar{a}, j\bar{a})}$, the natural spring length as $\ell_\alpha = \|\delta_\alpha\|$, the undeformed unit tangent $\mathbf{t}_\alpha = \delta_\alpha / \ell_\alpha$ and the linearized spring elongation $z_\alpha$ as

$$ z_\alpha = \mathbf{t}_\alpha \cdot (\mathbf{v}_{(i\bar{a}^+, j\bar{a}^+)} - \mathbf{v}_{(i\bar{a}, j\bar{a})}). \qquad (7.15) $$

The spring is assigned an elastic energy $w_\alpha = \frac{1}{2} k_\alpha z_\alpha^2$, where $k_\alpha$ is the spring constant which we set up as

$$ k_\alpha = \frac{EA}{\ell_\alpha} \qquad (7.16) $$

to represent elastic bars having all the same traction modulus EA. The short edge $j=n_2$ is blocked and a force

$$ \mathbf{g}_{(i,j)} = -\Delta (\eta r_i)^2 \mathbf{e}_2 \qquad (7.17) $$

is applied on the node $(i, j)$, representing a weight-like force having a constant density $(-\Delta \mathbf{e}_2)$ per unit area, where $r_i = \exp(\eta i)$ is the radial coordinate of the node. The nodal displacements $\mathbf{v}_\beta$ at equilibrium are proportional to $\Delta$ and independent of the global elasticity constant EA. The invariable forces $\mathbf{g}_{(i,j)}$ are taken into account through their potential energy $-\sum_{i=-n_1}^{+n_1} \sum_{j=0}^{n_2} \mathbf{g}_{(i,j)} \cdot \mathbf{v}_{(i,j)}$.
Typical results of the discrete simulations are shown in Figure 7.7 (deformed lattice), along with a comparison to the predictions of homogenization (colors) which we proceed to derive now.

We use homogenization to characterize the limit where both discretization parameters $n_1$ and $n_2$ are large: the scale separation parameter $\eta = \mathcal{O}(1/n_2)$ is then small, $\eta \ll 1$. We limit attention to the case of a stubby arch, and not a slender one, by assuming that the aspect-ratio parameter $\tilde{\delta} = \mathcal{O}(n_1/n_2)$ remains finite.

The periodic topological lattice underlying the arch is shown in Figure 7.6a. Its nodes are positioned on a square lattice spanning a rectangular domain, $\tilde{X}_{(i,j)} = i e_1 + j e_2$, with $-n_1 \le i \le +n_1$ and $0 \le j \le n_2$. There are $n_\varphi = 6$ element families (colors in Figure 7.6a) whose properties are listed in Table 7.2. In this imaginary configuration, the area density of the elements belonging to any particular family $\varphi$ is $\tilde{\rho}_\varphi = 1/2$.

<table><caption>**Table 7.2.** Properties of the topological square lattice underlying the arch shown Figure 7.6a: space dimension $d$, number $n_b$ of Bravais sub-lattices, number $n_\varphi$ of element (spring) families, element density $\tilde{\rho}_\varphi$, connectivity $b_\varphi^\pm$, and end-to-end vectors $\tilde{\delta}_\varphi$.</caption><thead><tr><th></th><th><i>d</i></th><th><i>n</i><sub>b</sub></th><th><i>n</i><sub>ϕ</sub></th><th><i>ρ̃</i><sub>ϕ</sub></th><th><i>I</i></th><th><i>II</i></th><th><i>III</i></th><th><i>IV</i></th><th><i>V</i></th><th><i>VI</i></th></tr></thead><tbody><tr><th><i>b</i><sub>ϕ</sub><sup>-</sup></th><td></td><td></td><td></td><td></td><td>1</td><td>2</td><td>2</td><td>1</td><td>1</td><td>1</td></tr><tr><th><i>b</i><sub>ϕ</sub><sup>+</sup></th><td></td><td></td><td></td><td></td><td>2</td><td>1</td><td>1</td><td>2</td><td>1</td><td>1</td></tr><tr><th><i>δ̃</i><sub>ϕ</sub></th><td></td><td></td><td></td><td></td><td><b>e</b><sub>2</sub></td><td><b>e</b><sub>2</sub></td><td><b>e</b><sub>1</sub></td><td><b>e</b><sub>1</sub></td><td><b>e</b><sub>1</sub> + <b>e</b><sub>2</sub></td><td><b>e</b><sub>1</sub> - <b>e</b><sub>2</sub></td></tr></tbody></table>

<figure><img src="image_24.png" /><figcaption>**Figure 7.6.** Circular arch: homogenization conventions. (a) Underlying periodic topological lattice, with the different element families represented by colors. (b) Close-up view of particular element of type $\varphi=I$ in reference configuration, drawn with a large value of the cell size $\sim \eta$ for the sake of legibility. Note that the element center $\mathbf{X}_I^c$ (target-like symbol, defined as the image of the element midpoint $\tilde{\mathbf{X}}_I^f$ by the diffeomorphism $f$) is distinct from the spring midpoint.</figcaption></figure>

The node positions $\mathbf{X}_{(i,j)}$ in (7.14) are produced by applying the map $\mathbf{X}_{(i,j)} = f(\eta \tilde{\mathbf{X}}_{(i,j)})$ in (2.3) to the topological lattice $\tilde{\mathbf{X}}_{(i,j)}$ using an exponential-map diffeomorphism $f$ similar to (2.4),

$$f(\tilde{\mathbf{x}}) = \exp(\tilde{x}_1) e_r(\tilde{x}_2). \qquad (7.18)$$

The arch therefore complies with the lattice generation framework presented in Section 2. The Jacobian can be calculated as $J_f(\tilde{\mathbf{x}}) = \exp(2\tilde{x}_1)$, and Equation (3.40) yields the density of springs in a given family $\varphi$ as

$$\rho_{\varphi}(r, \theta) = \frac{1}{2(\eta r)^2}, \qquad (7.19)$$

where $(r, \theta)$ are the polar coordinates in reference configuration, $\mathbf{X} = r e_r(\theta)$.

As shown by the target-like symbols in Figures 2.1b' and 7.6b, the link centers $\mathbf{X}_\alpha^c$, conventionally defined as the image by the diffeomorphism of the edge midpoints, see (2.5), do not coincide with the midpoints. The vectors $\Delta_\varphi^\pm$ yielding the endpoints $\pm$ relative to the center $\mathbf{X}_\alpha^c$ (thin black arrows in Figure 7.6b) are given by Equation (2.7) in terms of the family index $\varphi$, the scale separation parameter $\eta$ and of the polar coordinates $(r, \theta)$ of the center $\mathbf{X}_\alpha^c$ as

$$\Delta_{\varphi}^{\pm}(\eta, r, \theta) = r \exp\left(\pm \frac{\eta}{2} \tilde{\delta}_{\varphi} \cdot e_1\right) e_r\left(\theta \pm \frac{\eta}{2} \tilde{\delta}_{\varphi} \cdot e_2\right) - r e_r(\theta), \qquad (7.20)$$

where the edge vectors $\tilde{\delta}_\varphi$ of the topological lattice are listed in Table 7.2.

In the homogenization code, the lattice is initialized with the information listed in Table 7.2, along with the expressions of $\rho_\varphi$ in (7.19) and $\Delta_\varphi^\pm$ in (7.20).

The following changes are made to the code to handle spring elements, as opposed to the beam elements that have been used so far. The rotational degrees of freedom $\psi_b$ are dropped. There is a single strain ($n_{E_\varphi} = 1$) attached to a spring $a$, namely the elongation $z_a$: the element strain in (2.12) becomes $\mathbf{E}_a = (z_a)$. In view of the expression (7.15) of the elongation, we change the displacement-strain matrix in (2.14–2.15) to

$$\mathcal{D}_{\varphi}(\eta, r, \theta) = \begin{pmatrix} -\dot{\mathbf{i}}_{\varphi}(\eta, r, \theta) \\ +\dot{\mathbf{i}}_{\varphi}(\eta, r, \theta) \end{pmatrix} \in \mathbf{T}^{(1,2,2)} \quad (2D \text{ springs}), \qquad (7.21)$$
Typical results of the discrete simulations are shown in Figure 7.7 (deformed lattice), along with a comparison to the predictions of homogenization (colors) which we proceed to derive now.

We use homogenization to characterize the limit where both discretization parameters $n_1$ and $n_2$ are large: the scale separation parameter $\eta = \mathcal{O}(1/n_2)$ is then small, $\eta \ll 1$. We limit attention to the case of a stubby arch, and not a slender one, by assuming that the aspect-ratio parameter $\tilde{\delta} = \mathcal{O}(n_1/n_2)$ remains finite.

The periodic topological lattice underlying the arch is shown in Figure 7.6a. Its nodes are positioned on a square lattice spanning a rectangular domain, $\tilde{\mathbf{X}}_{(i,j)} = i\mathbf{e}_1 + j\mathbf{e}_2$, with $-n_1 \le i \le +n_1$ and $0 \le j \le n_2$. There are $n_\varphi = 6$ element families (colors in Figure 7.6a) whose properties are listed in Table 7.2. In this imaginary configuration, the area density of the elements belonging to any particular family $\varphi$ is $\tilde{\rho}_\varphi = 1/2$.

<table><caption>**Table 7.2.** Properties of the topological square lattice underlying the arch shown Figure 7.6a: space dimension $d$, number $n_b$ of Bravais sub-lattices, number $n_\varphi$ of element (spring) families, element density $\tilde{\rho}_\varphi$, connectivity $b_\varphi^\pm$, and end-to-end vectors $\tilde{\delta}_\varphi$.</caption><thead><tr><th></th><th><i>d</i></th><th><i>n</i><sub>b</sub></th><th><i>n</i><sub>φ</sub></th><th><i>ρ̃</i><sub>φ</sub></th><th><i>I</i></th><th><i>II</i></th><th><i>III</i></th><th><i>IV</i></th><th><i>V</i></th><th><i>VI</i></th></tr></thead><tbody><tr><th><i>b</i><sub>φ</sub><sup>-</sup></th><td></td><td></td><td></td><td></td><td>1</td><td>2</td><td>2</td><td>1</td><td>1</td><td>1</td></tr><tr><th><i>b</i><sub>φ</sub><sup>+</sup></th><td></td><td></td><td></td><td></td><td>2</td><td>1</td><td>1</td><td>2</td><td>1</td><td>1</td></tr><tr><th><i>δ̃</i><sub>φ</sub></th><td></td><td></td><td></td><td></td><td><b>e</b><sub>2</sub></td><td><b>e</b><sub>2</sub></td><td><b>e</b><sub>1</sub></td><td><b>e</b><sub>1</sub></td><td><b>e</b><sub>1</sub> + <b>e</b><sub>2</sub></td><td><b>e</b><sub>1</sub> - <b>e</b><sub>2</sub></td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>**Figure 7.6.** Circular arch: homogenization conventions. (a) Underlying periodic topological lattice, with the different element families represented by colors. (b) Close-up view of particular element of type $\varphi=I$ in reference configuration, drawn with a large value of the cell size $\sim \eta$ for the sake of legibility. Note that the element center $\mathbf{X}_I^c$ (target-like symbol, defined as the image of the element midpoint $\tilde{\mathbf{X}}_I^j$ by the diffeomorphism $f$) is distinct from the spring midpoint.</figcaption></figure>

The node positions $\mathbf{X}_{(i,j)}$ in (7.14) are produced by applying the map $\mathbf{X}_{(i,j)} = f(\eta \tilde{\mathbf{X}}_{(i,j)})$ in (2.3) to the topological lattice $\tilde{\mathbf{X}}_{(i,j)}$ using an exponential-map diffeomorphism $f$ similar to (2.4),

$$f(\tilde{\mathbf{x}}) = \exp(\tilde{x}_1) e_r(\tilde{x}_2). \qquad (7.18)$$

The arch therefore complies with the lattice generation framework presented in Section 2. The Jacobian can be calculated as $J_f(\tilde{\mathbf{x}}) = \exp(2\tilde{x}_1)$, and Equation (3.40) yields the density of springs in a given family $\varphi$ as

$$\rho_{\varphi}(r, \theta) = \frac{1}{2(\eta r)^2}, \qquad (7.19)$$

where $(r, \theta)$ are the polar coordinates in reference configuration, $\mathbf{X} = r e_r(\theta)$.

As shown by the target-like symbols in Figures 2.1b' and 7.6b, the link centers $\mathbf{X}_\alpha^c$, conventionally defined as the image by the diffeomorphism of the edge midpoints, see (2.5), do not coincide with the midpoints. The vectors $\Delta_\varphi^\pm$ yielding the endpoints $\pm$ relative to the center $\mathbf{X}_\alpha^c$ (thin black arrows in Figure 7.6b) are given by Equation (2.7) in terms of the family index $\varphi$, the scale separation parameter $\eta$ and of the polar coordinates $(r, \theta)$ of the center $\mathbf{X}_\alpha^c$ as

$$\Delta_{\varphi}^{\pm}(\eta, r, \theta) = r \exp\left(\pm\frac{\eta}{2} \tilde{\delta}_{\varphi} \cdot \mathbf{e}_1\right) e_r\left(\theta \pm \frac{\eta}{2} \tilde{\delta}_{\varphi} \cdot \mathbf{e}_2\right) - r e_r(\theta), \qquad (7.20)$$

where the edge vectors $\tilde{\delta}_\varphi$ of the topological lattice are listed in Table 7.2.

In the homogenization code, the lattice is initialized with the information listed in Table 7.2, along with the expressions of $\rho_\varphi$ in (7.19) and $\Delta_\varphi^\pm$ in (7.20).

The following changes are made to the code to handle spring elements, as opposed to the beam elements that have been used so far. The rotational degrees of freedom $\psi_b$ are dropped. There is a single strain ($n_{E_\varphi} = 1$) attached to a spring $a$, namely the elongation $z_a$: the element strain in (2.12) becomes $\mathbf{E}_a = (z_a)$. In view of the expression (7.15) of the elongation, we change the displacement-strain matrix in (2.14–2.15) to

$$\mathcal{D}_{\varphi}(\eta, r, \theta) = \left( \begin{pmatrix} -\dot{\mathbf{i}}_{\varphi}(\eta, r, \theta) \\ +\dot{\mathbf{i}}_{\varphi}(\eta, r, \theta) \end{pmatrix} \right) \in \mathbf{T}^{(1,2,2)} \quad (2D \text{ springs}), \qquad (7.21)$$
where $\hat{\boldsymbol{t}}_\varphi$ is still defined by (2.13). The choice of the spring constants in (7.16) is captured by setting up the element stiffness matrix as

$$ \mathcal{H}_{\varphi}(\eta, r, \theta) = \left( \frac{EA}{\hat{\ell}_{\varphi}(\eta, r, \theta)} \right) \in \mathbb{T}^{(1,1)} \quad (\text{springs}), \qquad (7.22) $$

so that (2.20) yields the elastic energy of the springs.

As anticipated in Section 2, the properties (such as $\mathcal{H}_\varphi$, $\rho_\varphi$, etc.) of this curved lattice depend on position $\mathbf{X}$ through the polar coordinates $(r, \theta)$. To inform the library that the symbols $r$ and $\theta$ vary in space, we set up the list of variable material parameters as

$$ \mathbf{m} = (r, \theta). \qquad (7.23) $$

The library furnishes the homogenized tensors in terms $\mathbf{m}$ and $\nabla \mathbf{m}$, treated as symbolic quantities: in the companion Mathematica notebook homogenize-circulararch nb, we insert the explicit expression of the gradient $\nabla \mathbf{m} = (\mathbf{e}_r(\theta), \mathbf{e}_\theta(\theta)/r)$ in this symbolic output, where $\mathbf{e}_\theta(\theta) = -\sin\theta\mathbf{e}_1 + \cos\theta\mathbf{e}_2$ is the azimuthal unit vector.

The homogenization proceeds automatically as usual, except that the symbolic expressions of $\mathbf{E}_l(\eta; r, \theta)$, ..., $\mathbf{E}_l''(\eta; r, \theta)$ obtained by (3.23) are so long that they make the procedure hang: to work around this, we replace $\mathbf{E}_l(\eta; r, \theta)$, ..., $\mathbf{E}_l''(\eta; r, \theta)$ by their Taylor expansions with respect to $\eta$ to order $\eta^3$, which are considerably simpler. These quantities being all of order $\eta$, expanding to order $\eta^3$ warrants a relative accuracy of order $\eta^2$ is preserved, which is consistent with our second-order homogenization scheme.

The second-order homogenization procedure delivers the microscopic displacement in the form

$$ \xi_b^* = (-1)^{b+1} r \eta \left( 0 \eta^0 + \frac{\eta^1}{8\sqrt{2}} \left[ \left( -\operatorname{tr}\boldsymbol{\epsilon} + r \frac{\partial \operatorname{tr}\boldsymbol{\epsilon}}{\partial r} + 2 \frac{\partial \varepsilon_{r\theta}}{\partial \theta} \right) \mathbf{e}_r + \left( 2 \varepsilon_{r\theta} + \frac{\partial \operatorname{tr}\boldsymbol{\epsilon}}{\partial \theta} + 2r \frac{\partial \varepsilon_{r\theta}}{\partial r} \right) \mathbf{e}_\theta \right] + \mathcal{O}(\eta^2) \right). \quad (7.24) $$

where $\operatorname{tr}\boldsymbol{\epsilon} = \varepsilon_{rr} + \varepsilon_{\theta\theta}$, and $\varepsilon_{rr}$, $\varepsilon_{\theta\theta}$ and $\varepsilon_{r\theta}$ are the components of the macroscopic strain in the polar basis $(\mathbf{e}_r, \mathbf{e}_\theta)$. This prediction is verified by comparing to discrete simulations following the usual procedure in Figure 7.7. The vanishing contribution $0\eta^0$ in the parenthesis in (7.24) points to the fact that the leading-order microscopic displacement is zero for this particular truss: it does not help improving the error, $e_{[0]}^a = e_{[\times]}^a$ (Figure 7.7, top). The higher-order term in (7.24) improves the error substantially and uniformly in the interior (Figure 7.7, bottom).

<figure><img src="image_25.png" /><figcaption>**Figure 7.7.** Circular arch subjected to a weight-like force $-\Delta \mathbf{e}_2$ per unit area, with discretization parameters $(n_1, n_2) = (13, 90)$ and loading intensity $\Delta = 0.1$. The logarithmic color maps represents the error on the microscopic displacement: the first non-zero prediction from homogenization in (7.24) comes from second-order (bottom).</figcaption></figure>

**Remark 7.4.** The microscopic displacement in (7.24) can be rewritten as

$$ \xi_b^* = (-1)^{b+1} r \eta \left( 0 \eta^0 + \frac{\eta^1}{8\sqrt{2}} \left[ (\mathcal{F}_p + \mathbf{e}_r \otimes (\mathcal{F}_r - \mathbf{I}_2)) : \boldsymbol{\epsilon} + r (\nabla \operatorname{tr}\boldsymbol{\epsilon} + \mathcal{F}_\theta \cdot \nabla (\mathcal{F}_\theta : \boldsymbol{\epsilon})) \right] + \mathcal{O}(\eta^2) \right), \quad (7.25) $$
where $\hat{\boldsymbol{t}}_\varphi$ is still defined by (2.13). The choice of the spring constants in (7.16) is captured by setting up the element stiffness matrix as

$$ \mathcal{H}_{\varphi}(\eta, r, \theta) = \left( \frac{EA}{\hat{\ell}_{\varphi}(\eta, r, \theta)} \right) \in \mathbb{T}^{(1,1)} \quad (\text{springs}), \qquad (7.22) $$

so that (2.20) yields the elastic energy of the springs.

As anticipated in Section 2, the properties (such as $\mathcal{H}_\varphi$, $\rho_\varphi$, etc.) of this curved lattice depend on position $\mathbf{X}$ through the polar coordinates $(r, \theta)$. To inform the library that the symbols $r$ and $\theta$ vary in space, we set up the list of variable material parameters as

$$ \boldsymbol{m} = (r, \theta). \qquad (7.23) $$

The library furnishes the homogenized tensors in terms $\boldsymbol{m}$ and $\nabla \boldsymbol{m}$, treated as symbolic quantities: in the companion Mathematica notebook homogenize-circulararch nb, we insert the explicit expression of the gradient $\nabla \boldsymbol{m} = (e_r(\theta), e_\theta(\theta)/r)$ in this symbolic output, where $e_\theta(\theta) = -\sin\theta e_1 + \cos\theta e_2$ is the azimuthal unit vector.

The homogenization proceeds automatically as usual, except that the symbolic expressions of $E_l(\eta; r, \theta)$, ..., $E_l''(\eta; r, \theta)$ obtained by (3.23) are so long that they make the procedure hang: to work around this, we replace $E_l(\eta; r, \theta)$, ..., $E_l''(\eta; r, \theta)$ by their Taylor expansions with respect to $\eta$ to order $\eta^3$, which are considerably simpler. These quantities being all of order $\eta$, expanding to order $\eta^3$ warrants a relative accuracy of order $\eta^2$ is preserved, which is consistent with our second-order homogenization scheme.

The second-order homogenization procedure delivers the microscopic displacement in the form

$$ \xi_b^* = (-1)^{b+1} r \eta \left( 0 \eta^0 + \frac{\eta^1}{8\sqrt{2}} \left[ \left( -\operatorname{tr}\boldsymbol{\varepsilon} + r \frac{\partial \operatorname{tr}\boldsymbol{\varepsilon}}{\partial r} + 2 \frac{\partial \varepsilon_{r\theta}}{\partial \theta} \right) e_r + \left( 2 \varepsilon_{r\theta} + \frac{\partial \operatorname{tr}\boldsymbol{\varepsilon}}{\partial \theta} + 2r \frac{\partial \varepsilon_{r\theta}}{\partial r} \right) e_\theta \right] + \mathcal{O}(\eta^2) \right). \quad (7.24) $$

where $\operatorname{tr}\boldsymbol{\varepsilon} = \varepsilon_{rr} + \varepsilon_{\theta\theta}$, and $\varepsilon_{rr}$, $\varepsilon_{\theta\theta}$ and $\varepsilon_{r\theta}$ are the components of the macroscopic strain in the polar basis $(e_r, e_\theta)$. This prediction is verified by comparing to discrete simulations following the usual procedure in Figure 7.7. The vanishing contribution $0\eta^0$ in the parenthesis in (7.24) points to the fact that the leading-order microscopic displacement is zero for this particular truss: it does not help improving the error, $e_{[0]}^a = e_{[1]}^a$ (Figure 7.7, top). The higher-order term in (7.24) improves the error substantially and uniformly in the interior (Figure 7.7, bottom).

<figure><img src="image_27.png" /><figcaption>**Figure 7.7.** Circular arch subjected to a weight-like force $-\Delta e_2$ per unit area, with discretization parameters $(n_1, n_2) = (13, 90)$ and loading intensity $\Delta = 0.1$. The logarithmic color maps represents the error on the microscopic displacement: the first non-zero prediction from homogenization in (7.24) comes from second-order (bottom).</figcaption></figure>

**Remark 7.4.** The microscopic displacement in (7.24) can be rewritten as

$$ \xi_b^* = (-1)^{b+1} r \eta \left( 0 \eta^0 + \frac{\eta^1}{8\sqrt{2}} \left[ (\mathcal{F}_p + e_r \otimes (\mathcal{F}_r - \mathbf{I}_2)) : \boldsymbol{\varepsilon} + r (\nabla \operatorname{tr}\boldsymbol{\varepsilon} + \mathcal{F}_\theta \cdot \nabla (\mathcal{F}_\theta : \boldsymbol{\varepsilon})) \right] + \mathcal{O}(\eta^2) \right), \quad (7.25) $$
where $\mathcal{F}_p$ is the polar version of the tensor $\mathcal{F}$ introduced in (5.5), $\mathcal{F}_p = e_r \otimes \mathcal{F}_r + e_\theta \otimes \mathcal{F}_\theta$ with $\mathcal{F}_r = -e_r \otimes e_r + e_\theta \otimes e_\theta$ and $\mathcal{F}_\theta = e_r \otimes e_\theta + e_\theta \otimes e_r$. In the presence of variable properties, the corrective displacement not only includes a term proportional to $\nabla \epsilon$ but also one proportional to $\epsilon$ (see the first term in the square bracket above). The latter is produced by inserting $\nabla m = (e_r(\theta), e_\theta(\theta)/r)$ into the corrective $\nabla m \otimes \epsilon$ term returned by the library in symbolic form.

## 7.4. Shearing and twisting of a 3D elastic truss

We now address an example in dimension $d=3$, namely an elastic lattice made up of springs whose unit cell is shown in Figure 7.8. The unit cell dimensions is $\ell \times \ell \times \sqrt{3} \ell$. There are $n_\varphi = 11$ different types of bars: their length is $\ell_\alpha = \ell$ for types $\varphi \in \{I, VI\}$, $\ell_\alpha = \sqrt{3} \ell$ for type $\varphi = XI$, and $\ell_\alpha = \sqrt{5} \ell/2$ for the other types $\varphi$. The spring constants $k_\alpha$ in Equation (7.16) are used, representing elastic bars have all the same stretching modulus $EA$. As earlier in 2D, the spring elongation $z_\alpha$ is given by (7.15) and their energy is of the form $w_\alpha = k_\alpha z_\alpha^2/2$.

<figure><img src="image_26.png" /><figcaption>**Figure 7.8.** Unit cell of a 3D lattice made up of springs, with box dimensions $\ell \times \ell \times \sqrt{3} \ell$. The lattice has $n_b = 2$ Bravais sub-lattices (solid vs. open disks) and $n_\varphi = 11$ element families (colors). The springs are treated as elastic bars having identical stretching modulus $EA$.</figcaption></figure>

In Figure 7.9, a shear test is conducted on a block comprising $20 \times 10 \times 10$ unit cells: the displacement of the nodes belonging to the face shown in the left-side is blocked, whereas that of those belonging to the opposite face is set to $\Delta e_3$, with $\Delta = 0.1 \ell$. In Figure 7.10, a twisting test is conducted on a block comprising $30 \times 15 \times 11$ unit cells: a rotation by an angle $\Delta$ is imposed on the second face, with $\Delta = 0.1$ (only the components of the displacement in the perpendicular plane $(e_2, e_3)$ are prescribed, the longitudinal component along $e_1$ remaining unconstrained).

$$ \frac{d}{3} \qquad n_b \qquad n_\varphi \qquad \frac{\rho_\varphi}{(\sqrt{3} \ell^3)^{-1}} $$

<caption>**Table 7.3.** Parameters used to set up the 3D elastic truss in the homogenization code, see Figures 7.8–7.10. The connectivity matrix $b_\varphi^\pm$ and the edge vectors $\delta_\varphi$ are omitted as they can be read off easily from Figure 7.8.</caption>

The homogenization of the 3D truss is carried out in the companion notebook homogenize-3Dtruss nb. The
where $\mathcal{F}_p$ is the polar version of the tensor $\mathcal{F}$ introduced in (5.5), $\mathcal{F}_p = e_r \otimes \mathcal{F}_r + e_\theta \otimes \mathcal{F}_\theta$ with $\mathcal{F}_r = -e_r \otimes e_r + e_\theta \otimes e_\theta$ and $\mathcal{F}_\theta = e_r \otimes e_\theta + e_\theta \otimes e_r$. In the presence of variable properties, the corrective displacement not only includes a term proportional to $\nabla \epsilon$ but also one proportional to $\epsilon$ (see the first term in the square bracket above). The latter is produced by inserting $\nabla m = (e_r(\theta), e_\theta(\theta)/r)$ into the corrective $\nabla m \otimes \epsilon$ term returned by the library in symbolic form.

## 7.4. Shearing and twisting of a 3D elastic truss

We now address an example in dimension $d=3$, namely an elastic lattice made up of springs whose unit cell is shown in Figure 7.8. The unit cell dimensions is $\ell \times \ell \times \sqrt{3} \ell$. There are $n_\varphi = 11$ different types of bars: their length is $\ell_\alpha = \ell$ for types $\varphi \in \{I, VI\}$, $\ell_\alpha = \sqrt{3} \ell$ for type $\varphi = XI$, and $\ell_\alpha = \sqrt{5} \ell/2$ for the other types $\varphi$. The spring constants $k_\alpha$ in Equation (7.16) are used, representing elastic bars have all the same stretching modulus $EA$. As earlier in 2D, the spring elongation $z_\alpha$ is given by (7.15) and their energy is of the form $w_\alpha = k_\alpha z_\alpha^2/2$.

<figure><img src="image_28.png" /><figcaption>**Figure 7.8.** Unit cell of a 3D lattice made up of springs, with box dimensions $\ell \times \ell \times \sqrt{3} \ell$. The lattice has $n_b = 2$ Bravais sub-lattices (solid vs. open disks) and $n_\varphi = 11$ element families (colors). The springs are treated as elastic bars having identical stretching modulus $EA$.</figcaption></figure>

In Figure 7.9, a shear test is conducted on a block comprising $20 \times 10 \times 10$ unit cells: the displacement of the nodes belonging to the face shown in the left-side is blocked, whereas that of those belonging to the opposite face is set to $\Delta e_3$, with $\Delta = 0.1 \ell$. In Figure 7.10, a twisting test is conducted on a block comprising $30 \times 15 \times 11$ unit cells: a rotation by an angle $\Delta$ is imposed on the second face, with $\Delta = 0.1$ (only the components of the displacement in the perpendicular plane $(e_2, e_3)$ are prescribed, the longitudinal component along $e_1$ remaining unconstrained).

$$ \frac{d}{3} \qquad n_b \qquad n_\varphi \qquad \rho_\varphi $$

<caption>**Table 7.3.** Parameters used to set up the 3D elastic truss in the homogenization code, see Figures 7.8–7.10. The connectivity matrix $b_\varphi^\pm$ and the edge vectors $\delta_\varphi$ are omitted as they can be read off easily from Figure 7.8.</caption>

The homogenization of the 3D truss is carried out in the companion notebook homogenize-3Dtruss nb. The
lattice is set up using the parameters listed in Table 7.3 (the connectivity matrix $b_{\phi}^{\pm}$ and the edge vectors $\delta_{\phi}$ are provided to the library but are not included in the table). To represent the springs, the 3D version of the displacement-to-strain matrix $\mathcal{D}_{\phi}$ in (7.21) is used, along with the element stiffness $\mathcal{K}_{\phi}$ in (7.22). For this particular truss, the leading-order microscopic displacement vanishes, $\mathbf{Y}_0=\mathbf{0}$, and the first-order one is returned by the library in the form

$$ \xi_b^* = (-1)^{b+1} \ell \left( 0 \ell^0 + \frac{5\sqrt{5}}{64} \ell \left( -\epsilon_{11,1} \mathbf{e}_1 + \epsilon_{22,2} \mathbf{e}_2 + \frac{1}{\sqrt{3}} \epsilon_{33,3} \mathbf{e}_3 \right) + \mathcal{O}(\ell^2) \right). \quad (7.26) $$

This prediction is verified against numerical simulations in Figures 7.9 and 3.4, by the usual procedure. As in all the previous examples, a significant improvement is brought about by second-order homogenization, except in the vicinity of the edges.

<figure><img src="image_27.png" /><figcaption>Figure 7.9. 3D truss, shearing test. The springs are colored according to the magnitude of the microscopic displacement (no homogenization) or of the error the displacement predicted by first-order or second-order homogenization. To aid visualization, a single layer of springs is shown, the full lattice domain being shown by the black box.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 7.10. 3D truss, twisting test. Same conventions as in Figure 7.9. On the face shown in foreground, the rotation with angle $\Delta$ is imposed on the $\mathbf{e}_2$ and $\mathbf{e}_3$ (normal) components of the displacement, the $\mathbf{e}_1$ (longitudinal) component remains unconstrained.</figcaption></figure>

## 7.5. A pantograph

Our final example is the 1D elastic truss shown in Figure 7.11, known as a pantograph. It is made up of springs having identical elastic constants $k$, except for the vertical springs whose elastic constant is $\chi k \ll k$; we address the case of a small elastic contrast, $\chi \ll 1$. The homogenization of pantograph-like lattices has been studied in a number of work, see for instance [SAd11; AS18b]. More broadly, microstructures featuring large elastic contrast have been extensively studied, see for instance the work of [PS97; Bel17; JS20]. By recovering known results about the pantograph, we demonstrate how our method can handle the special case of a large elastic contrast.

In all the previous illustrations given in this paper, the gradient effect was a correction to the leading order effect—our homogenization scheme has indeed been derived with this situation in mind. The pantograph is different as the strain-gradient term $d\epsilon_{11}/dX_1$ ends up contributing to its energy at the same order as the strain $\epsilon_{11}$, see Equation (7.29) below. Our homogenization scheme can be revisited to deal with this situation, as we show now.
lattice is set up using the parameters listed in Table 7.3 (the connectivity matrix $b_{\phi}^{\pm}$ and the edge vectors $\delta_{\phi}$ are provided to the library but are not included in the table). To represent the springs, the 3D version of the displacement-to-strain matrix $\mathcal{D}_{\phi}$ in (7.21) is used, along with the element stiffness $\mathcal{K}_{\phi}$ in (7.22). For this particular truss, the leading-order microscopic displacement vanishes, $\mathbf{Y}_0=\mathbf{0}$, and the first-order one is returned by the library in the form

$$ \xi_b^* = (-1)^{b+1} \ell \left( 0 \ell^0 + \frac{5\sqrt{5}}{64} \ell \left( -\epsilon_{11,1} \mathbf{e}_1 + \epsilon_{22,2} \mathbf{e}_2 + \frac{1}{\sqrt{3}} \epsilon_{33,3} \mathbf{e}_3 \right) + \mathcal{O}(\ell^2) \right). \quad (7.26) $$

This prediction is verified against numerical simulations in Figures 7.9 and 3.4, by the usual procedure. As in all the previous examples, a significant improvement is brought about by second-order homogenization, except in the vicinity of the edges.

<figure><img src="image_29.png" /><figcaption>**Figure 7.9.** 3D truss, shearing test. The springs are colored according to the magnitude of the microscopic displacement (no homogenization) or of the error the displacement predicted by first-order or second-order homogenization. To aid visualization, a single layer of springs is shown, the full lattice domain being shown by the black box.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 7.10.** 3D truss, twisting test. Same conventions as in Figure 7.9. On the face shown in foreground, the rotation with angle $\Delta$ is imposed on the $\mathbf{e}_2$ and $\mathbf{e}_3$ (normal) components of the displacement, the $\mathbf{e}_1$ (longitudinal) component remains unconstrained.</figcaption></figure>

## 7.5. A pantograph

Our final example is the 1D elastic truss shown in Figure 7.11, known as a pantograph. It is made up of springs having identical elastic constants $k$, except for the vertical springs whose elastic constant is $\chi k \ll k$; we address the case of a small elastic contrast, $\chi \ll 1$. The homogenization of pantograph-like lattices has been studied in a number of work, see for instance [SA11; AS18b]. More broadly, microstructures featuring large elastic contrast have been extensively studied, see for instance the work of [PS97; Bel17; JS20]. By recovering known results about the pantograph, we demonstrate how our method can handle the special case of a large elastic contrast.

In all the previous illustrations given in this paper, the gradient effect was a correction to the leading order effect—our homogenization scheme has indeed been derived with this situation in mind. The pantograph is different as the strain-gradient term $d\epsilon_{11}/dX_1$ ends up contributing to its energy at the same order as the strain $\epsilon_{11}$, see Equation (7.29) below. Our homogenization scheme can be revisited to deal with this situation, as we show now.
<figure><img src="image_29.png" /><figcaption>**Figure 7.11.** A pantograph is a 1D truss made up of springs. All springs have the same constant $k$, except for the vertical dashed springs whose spring constant is $\chi k$. We consider the limit of large elastic contrast $\chi \ll 1$. There are $n_b=6$ Bravais sub-lattices (symbols and integers) and $n_\phi = 13$ different types of springs: a unit cell is shown using colored segments. In the limit $\chi=0$ where the dashed vertical springs are absent, the truss possesses a zero-energy mechanism. This is the same lattice as studied in [AS18b], except that beams have been replaced by springs and the weak (dashed vertical) springs have been added.</figcaption></figure>

The pantograph is homogenized in the companion Mathematica notebook *homogenize-pantograph-springs-highcontrast nb*. In is drawn in the Euclidean plane ($d=2$), see Figure 7.11, but the limit model is a one-dimensional bar. We could write specialized code handling 2D lattices that are 1D in the limit but a straightforward work-around is to work in the plane and imagine that we are homogenizing an infinite number of replicas of the pantograph obtained by copying the original one, translating it along $e_2$ and pasting it, repeatedly. Concretely, we set up the dimension as $d=2$ in the library and identify the density per area $\rho_\phi$ of each family of elements with its density per unit length $\rho_\phi = 1/\ell$, where $\ell = O(\eta)$ is the size of the unit cell. We also check that the homogenized energy computed by the library does not depend on $\epsilon_{12}$ and $\epsilon_{22}$ since the replicas are disconnected.

The homogenized energy returned by the library takes the form

$$ \Phi^* = \frac{k \ell}{2} \int_{-\infty}^{+\infty} \left( 4 \chi c_5(\chi) (\epsilon_{11})^2 + \frac{2}{23} c_6(\chi) \left( \ell \frac{d\epsilon_{11}}{dX_1} \right)^2 + \mathcal{O}(\eta^3) \right) dX_1, \quad (7.27) $$

where $c_5(\chi) = (1 + \frac{7}{2}\chi)^{-1}$ and $c_6(\chi) = (1 - \frac{19}{54}\chi - \frac{109}{27}\chi^2) / (1 + \frac{7}{2}\chi)^2$ are $\mathcal{C}^\infty$-smooth functions of $\chi$ at $\chi=0$, with $c_5(0) = c_6(0) = 1$. The underlying microscopic displacement $\xi_b^*$ is available from the companion notebook. We have not yet specified how $\chi$ scales with the cell size $\eta$, and the above result is correct as long as $\chi$ remains bounded, $\chi = \mathcal{O}(\eta^0)$, see Remark 2.2.

The limit of large elastic contrast is now addressed by taking the formal limit $\chi \to 0$ in (7.27). The distinguished limit where both $\chi \to 0$ and $\eta \to 0$ with

$$ \chi = \mathcal{O}(\ell^2) = \mathcal{O}(\eta^2) \quad (7.28) $$

is particularly interesting, as it makes the two terms balanced in the integral: at leading order, we then have

$$ \Phi^* = \frac{4k \ell \chi}{2} \int_{-\infty}^{+\infty} \left( (\epsilon_{11})^2 + \frac{1}{46} \left( \frac{\ell}{\chi^{1/2}} \frac{d\epsilon_{11}}{dX_1} \right)^2 + \mathcal{O}(\eta) \right) dX_1, \quad (7.29) $$

where $\ell/\chi^{1/2} = \mathcal{O}(1)$ is the macroscopic persistence length of the gradient effect. This energy can be interpreted based on the existence of a soft mode having zero energy in the limit $\chi \to 0$ (vertical dashed springs are removed). This soft mode involves combined rigid rotations of the triangles, as shown by the light gray circular arrows in Figure 7.11, with an amplitude proportional to the macroscopic strain $\epsilon_{11}$. The first term in the integral in (7.29) can be rewritten as $\frac{1}{2} \int (k\chi)(\ell\epsilon_{11})^2 dX_1/\ell$, which describes the small energy penalty brought about by the weak springs when the soft mode is activated. Classical homogenization captures only this term, yielding a non-definite energy functional $\frac{1}{2} \int 0 \times (\ell\epsilon_{11})^2 dX_1/\ell$ in the limit of large elastic contrast $\chi=0$. The gradient effect solves this issue.

The pantograph example shows that the special behavior of high-contrast microstructures, which has been documented in the literature, can be readily recovered by taking the formal limit $\chi \to 0$ in the output of the homogenization library. The *symbolic* and *second-order* capabilities of the library are instrumental here.

**Remark 7.5.** For consistency with the scaling assumption that the energy is $\Phi^* = \mathcal{O}(\eta^0)$, one should view the spring constant $k$ as quantity of order $1/\ell^3 = \mathcal{O}(\eta^{-3})$. However, the magnitude of this overall multiplicative constant in the energy is a matter of convention; nothing changes if different scaling assumptions are used.

# 8. DISCUSSION AND CONCLUSION

We have proposed a symbolic, second-order homogenization method for elastic lattices. It is distributed as a library, named *shoal*, implemented in the symbolic calculation language Wolfram Mathematica: given a description of the lattice, the library runs automatically and returns both the localization tensors in (4.3) and the homogenized energy (4.5) in symbolic form.
<figure><img src="image_31.png" /><figcaption>**Figure 7.11.** A pantograph is a 1D truss made up of springs. All springs have the same constant $k$, except for the vertical dashed springs whose spring constant is $\chi k$. We consider the limit of large elastic contrast $\chi \ll 1$. There are $n_b=6$ Bravais sub-lattices (symbols and integers) and $n_\varphi = 13$ different types of springs: a unit cell is shown using colored segments. In the limit $\chi=0$ where the dashed vertical springs are absent, the truss possesses a zero-energy mechanism. This is the same lattice as studied in [AS18b], except that beams have been replaced by springs and the weak (dashed vertical) springs have been added.</figcaption></figure>

The pantograph is homogenized in the companion Mathematica notebook *homogenize-pantograph-springs-highcontrast nb*. In is drawn in the Euclidean plane ($d=2$), see Figure 7.11, but the limit model is a one-dimensional bar. We could write specialized code handling 2D lattices that are 1D in the limit but a straightforward work-around is to work in the plane and imagine that we are homogenizing an infinite number of replicas of the pantograph obtained by copying the original one, translating it along $e_2$ and pasting it, repeatedly. Concretely, we set up the dimension as $d=2$ in the library and identify the density per area $\rho_\varphi$ of each family of elements with its density per unit length $\rho_\varphi = 1/\ell$, where $\ell = O(\eta)$ is the size of the unit cell. We also check that the homogenized energy computed by the library does not depend on $\varepsilon_{12}$ and $\varepsilon_{22}$ since the replicas are disconnected.

The homogenized energy returned by the library takes the form

$$ \Phi^* = \frac{k \ell}{2} \int_{-\infty}^{+\infty} \left( 4 \chi c_5(\chi) (\varepsilon_{11})^2 + \frac{2}{23} c_6(\chi) \left( \ell \frac{d\varepsilon_{11}}{dX_1} \right)^2 + \mathcal{O}(\eta^3) \right) dX_1, \quad (7.27) $$

where $c_5(\chi) = (1 + \frac{7}{2}\chi)^{-1}$ and $c_6(\chi) = (1 - \frac{19}{54}\chi - \frac{109}{27}\chi^2) / (1 + \frac{7}{2}\chi)^2$ are $\mathcal{C}^\infty$-smooth functions of $\chi$ at $\chi=0$, with $c_5(0) = c_6(0) = 1$. The underlying microscopic displacement $\xi_b^*$ is available from the companion notebook. We have not yet specified how $\chi$ scales with the cell size $\eta$, and the above result is correct as long as $\chi$ remains bounded, $\chi = \mathcal{O}(\eta^0)$, see Remark 2.2.

The limit of large elastic contrast is now addressed by taking the *formal* limit $\chi \to 0$ in (7.27). The distinguished limit where both $\chi \to 0$ and $\eta \to 0$ with

$$ \chi = \mathcal{O}(\ell^2) = \mathcal{O}(\eta^2) \quad (7.28) $$

is particularly interesting, as it makes the two terms balanced in the integral: at leading order, we then have

$$ \Phi^* = \frac{4k \ell \chi}{2} \int_{-\infty}^{+\infty} \left( (\varepsilon_{11})^2 + \frac{1}{46} \left( \frac{\ell}{\chi^{1/2}} \frac{d\varepsilon_{11}}{dX_1} \right)^2 + \mathcal{O}(\eta) \right) dX_1, \quad (7.29) $$

where $\ell/\chi^{1/2} = \mathcal{O}(1)$ is the macroscopic persistence length of the gradient effect. This energy can be interpreted based on the existence of a soft mode having zero energy in the limit $\chi \to 0$ (vertical dashed springs are removed). This soft mode involves combined rigid rotations of the triangles, as shown by the light gray circular arrows in Figure 7.11, with an amplitude proportional to the macroscopic strain $\varepsilon_{11}$. The first term in the integral in (7.29) can be rewritten as $\frac{1}{2} \int (k\chi)(\ell\varepsilon_{11})^2 dX_1/\ell$, which describes the small energy penalty brought about by the weak springs when the soft mode is activated. Classical homogenization captures only this term, yielding a non-definite energy functional $\frac{1}{2} \int 0 \times (\ell\varepsilon_{11})^2 dX_1/\ell$ in the limit of large elastic contrast $\chi=0$. The gradient effect solves this issue.

The pantograph example shows that the special behavior of high-contrast microstructures, which has been documented in the literature, can be readily recovered by taking the formal limit $\chi \to 0$ in the output of the homogenization library. The *symbolic* and *second-order* capabilities of the library are instrumental here.

**Remark 7.5.** For consistency with the scaling assumption that the energy is $\Phi^* = \mathcal{O}(\eta^0)$, one should view the spring constant $k$ as quantity of order $1/\ell^3 = \mathcal{O}(\eta^{-3})$. However, the magnitude of this overall multiplicative constant in the energy is a matter of convention; nothing changes if different scaling assumptions are used.

# 8. DISCUSSION AND CONCLUSION

We have proposed a symbolic, second-order homogenization method for elastic lattices. It is distributed as a library, named *shoal*, implemented in the symbolic calculation language Wolfram Mathematica: given a description of the lattice, the library runs automatically and returns both the localization tensors in (4.3) and the homogenized energy (4.5) in symbolic form.
In the design of the library, emphasis has been placed on versatility: it works in arbitrary dimension $d$ ($1 \le d \le 3$), can handle lattices featuring graded geometric and/or elastic properties, pre-stress or a large elastic contrast. Elastic lattices made up of springs or beams have been demonstrated and new types of elements can be created easily. We are not aware of previous work addressing both graded properties and pre-stress in the context of higher-order, discrete homogenization.

For a variety of lattices, we have derived in closed analytical form both the nodal displacements and rotations—see Equations (5.10), (7.3), (7.11), (7.25), (7.26)—and the homogenized energy—see Equations (5.9), (7.29)—up to second order. Expressions of this kind are novel to the best of our knowledge. Being symbolic, they fully convey the influence of the various lattice parameters. We hope that this kind of results will be used in the future to gain analytical insights into the effective properties of various lattices.

We have proposed a procedure for verifying the correctness (in an asymptotic sense) of the microscopic displacement predicted by the library, based on comparisons with discrete simulations, see Figure 6.4. We have presented a number of validation examples for which we have plotted physical maps of the homogenization error. The improvement in accuracy brought about by second-order homogenization is significant everywhere except in layers forming at free boundaries, where the assumption of scale separation breaks down.

Until these boundary layers have been analyzed and accounted for in the homogenized energy $\Phi^*$, it is not possible to solve global equilibrium problems for the lattice by making the functional $\Phi^*$ stationary. This is one main limitation of the present work, which is just a step towards this goal—a rather ambitious one in the context of higher-order homogenization. The derivation of the *other* boundary terms in (4.5)—produced by homogenization and not by boundary layers—is original to the best of our knowledge and will need to be included in the homogenized energy functional predicting equilibrium.

There is a second reason that prevents the homogenized energy $\Phi^*$ from being used in self-contained simulations, at least in its current form. As noted in Remark 5.3 and in our previous work [AL23], it contains negative strain-gradient moduli, such as the coefficient -7 appearing in the second line of (5.9), and is therefore non-positive. This is a common finding in higher-order homogenization [LM18; DLSS22]. This strong limitation explains the limited popularity of higher-order homogenization so far. In future work, we hope to introduce a regularization strategy that preserves the accuracy of the solution while restoring the desired positivity.

The library is *symbolic* and does *second-order* homogenization. When combined together, these two original features are quite powerful. Consider for instance the challenging class of lattice homogenization problems involving a small dimensionless parameter $\chi \ll 1$ (such as the slenderness ratio $EI / [EA \ell^2]$ in beam lattices or the elastic contrast, see §7.5) in addition to the scale separation parameter $\eta \ll 1$. Non-conventional limit behaviors can be obtained in such lattices, which have been approached using specialized methods in the literature. These limit behaviors are often associated with the presence of a large-amplitude, slowly decaying soft mode (Remark 7.3 in §7.2, §7.5). Our approach provides both a unifying perspective and a straightforward approach to these problems: by homogenizing $(\eta \to 0)$ with fixed $\chi$ first, and taking a distinguished limit in $\chi$ and $\eta$ in the homogenized energy next, one can identify the limit energy based on simple dimensional analysis (§7.5). This approach involves exchanging the limits on $\chi$ and $\eta$, which is not mathematically rigorous but yields correct results in all the cases we have tried.

Along the same lines, albeit at a more basic (and less original) level, the homogenization procedure handles bending-dominated and stretching-dominated lattices in a unified way: the class of a particular lattice can be deducted from the symbolic expression of the leading-order homogenized energy, see the discussion immediately after Equation (5.3).

The proposed method is subject to the following practical limitations. (i) The library solves linear algebra problems in symbolic form, which can require a large amount of computing time especially if there are many microscopic degrees of freedom. When this happens, it is possible to revert to *numeric* homogenization; this is done easily, by assigning numerical values to the lattice parameters. The symbolic computations required to generate the results included in this paper all took from a few seconds to less than a minute on a personal computer. (ii) The homogenized energy and microscopic displacement computed by the library have been cast by hand into compact tensor form, shown for instance in Equation (5.9–5.10): this process could in principle be automated by exploiting the symmetries of the lattice. (iii) We have homogenized discrete lattices made up of 1D beams of springs, without questioning the validity of these 1D models. To address this question, one needs to revert to the equations of *bulk* elasticity in a periodic domain comprising slender junctions [AB21; DLSS22].

In this work, we have limited attention to *linear elastic* lattices and thus have not considered material or geometrical nonlinearity. When treated incrementally, homogenization problems for lattices having a more general constitutive behavior (such as non-linearly elastic, visco-elastic or elasto-plastic) yield linearized problems featuring both pre-stress and spatially graded properties (such as the tangent elastic stiffness). This is precisely the kind of problems that our approach can solve. In future work, we plan to apply our homogenization scheme incrementally. By doing so, we hope to gain insights into the homogenizing of lattices featuring a wide range of constitutive behaviors.

**Acknowledgements** Yang Ye thanks the China Scholarship Council (CSC) for his doctoral funding. This manuscript has been typeset using *TEX*, a powerful, multi-platform and freely distributed scientific editor [vdHGG+13].

BIBLIOGRAPHY
In the design of the library, emphasis has been placed on versatility: it works in arbitrary dimension $d$ ($1 \le d \le 3$), can handle lattices featuring graded geometric and/or elastic properties, pre-stress or a large elastic contrast. Elastic lattices made up of springs or beams have been demonstrated and new types of elements can be created easily. We are not aware of previous work addressing both graded properties and pre-stress in the context of higher-order, discrete homogenization.

For a variety of lattices, we have derived in closed analytical form both the nodal displacements and rotations—see Equations (5.10), (7.3), (7.11), (7.25), (7.26)—and the homogenized energy—see Equations (5.9), (7.29)—up to second order. Expressions of this kind are novel to the best of our knowledge. Being symbolic, they fully convey the influence of the various lattice parameters. We hope that this kind of results will be used in the future to gain analytical insights into the effective properties of various lattices.

We have proposed a procedure for verifying the correctness (in an asymptotic sense) of the microscopic displacement predicted by the library, based on comparisons with discrete simulations, see Figure 6.4. We have presented a number of validation examples for which we have plotted physical maps of the homogenization error. The improvement in accuracy brought about by second-order homogenization is significant everywhere except in layers forming at free boundaries, where the assumption of scale separation breaks down.

Until these boundary layers have been analyzed and accounted for in the homogenized energy $\Phi^*$, it is not possible to solve global equilibrium problems for the lattice by making the functional $\Phi^*$ stationary. This is one main limitation of the present work, which is just a step towards this goal—a rather ambitious one in the context of higher-order homogenization. The derivation of the *other* boundary terms in (4.5)—produced by homogenization and not by boundary layers—is original to the best of our knowledge and will need to be included in the homogenized energy functional predicting equilibrium.

There is a second reason that prevents the homogenized energy $\Phi^*$ from being used in self-contained simulations, at least in its current form. As noted in Remark 5.3 and in our previous work [AL23], it contains negative strain-gradient moduli, such as the coefficient -7 appearing in the second line of (5.9), and is therefore non-positive. This is a common finding in higher-order homogenization [LM18; DLSS22]. This strong limitation explains the limited popularity of higher-order homogenization so far. In future work, we hope to introduce a regularization strategy that preserves the accuracy of the solution while restoring the desired positivity.

The library is *symbolic* and does *second-order* homogenization. When combined together, these two original features are quite powerful. Consider for instance the challenging class of lattice homogenization problems involving a small dimensionless parameter $\chi \ll 1$ (such as the slenderness ratio $EI / [EA \ell^2]$ in beam lattices or the elastic contrast, see §7.5) in addition to the scale separation parameter $\eta \ll 1$. Non-conventional limit behaviors can be obtained in such lattices, which have been approached using specialized methods in the literature. These limit behaviors are often associated with the presence of a large-amplitude, slowly decaying soft mode (Remark 7.3 in §7.2, §7.5). Our approach provides both a unifying perspective and a straightforward approach to these problems: by homogenizing $(\eta \to 0)$ with fixed $\chi$ first, and taking a distinguished limit in $\chi$ and $\eta$ in the homogenized energy next, one can identify the limit energy based on simple dimensional analysis (§7.5). This approach involves exchanging the limits on $\chi$ and $\eta$, which is not mathematically rigorous but yields correct results in all the cases we have tried.

Along the same lines, albeit at a more basic (and less original) level, the homogenization procedure handles bending-dominated and stretching-dominated lattices in a unified way: the class of a particular lattice can be deducted from the symbolic expression of the leading-order homogenized energy, see the discussion immediately after Equation (5.3).

The proposed method is subject to the following practical limitations. (i) The library solves linear algebra problems in symbolic form, which can require a large amount of computing time especially if there are many microscopic degrees of freedom. When this happens, it is possible to revert to *numeric* homogenization; this is done easily, by assigning numerical values to the lattice parameters. The symbolic computations required to generate the results included in this paper all took from a few seconds to less than a minute on a personal computer. (ii) The homogenized energy and microscopic displacement computed by the library have been cast by hand into compact tensor form, shown for instance in Equation (5.9–5.10): this process could in principle be automated by exploiting the symmetries of the lattice. (iii) We have homogenized discrete lattices made up of 1D beams of springs, without questioning the validity of these 1D models. To address this question, one needs to revert to the equations of *bulk* elasticity in a periodic domain comprising slender junctions [AB21; DLSS22].

In this work, we have limited attention to *linear elastic* lattices and thus have not considered material or geometrical nonlinearity. When treated incrementally, homogenization problems for lattices having a more general constitutive behavior (such as non-linearly elastic, visco-elastic or elasto-plastic) yield linearized problems featuring both pre-stress and spatially graded properties (such as the tangent elastic stiffness). This is precisely the kind of problems that our approach can solve. In future work, we plan to apply our homogenization scheme incrementally. By doing so, we hope to gain insights into the homogenizing of lattices featuring a wide range of constitutive behaviors.

**Acknowledgements** Yang Ye thanks the China Scholarship Council (CSC) for his doctoral funding. This manuscript has been typeset using *TEX*, a powerful, multi-platform and freely distributed scientific editor [vdHGG+13].

BIBLIOGRAPHY

[21] Wolfram Research, Inc. Mathematica, Version 13.0.0. 2021. Champaign, IL, 2021.
[AB21] Bilen Emek Abali and Emilio Barchiesi. Additive manufacturing introduced substructure and computational determination of metamaterials parameters by means of the asymptotic homogenization. *Continuum Mechanics and Thermodynamics*, 33(4):993–1009, 2021.

[ABB10] Nicolas Auffray, Regis Bouchet, and Y Brechet. Strain gradient elastic homogenization of bidimensional cellular media. *International Journal of Solids and Structures*, 47(13):1698–1710, 2010.

[AG97] Michael F Ashby and Lorna J Gibson. Cellular solids: structure and properties. *Press Syndicate of the University of Cambridge, Cambridge, UK*, pages 175–231, 1997.

[AK017] N. Auffray, B. Kolev, and M. Olive. Handbook of bi-dimensional tensors: part I: harmonic decomposition and symmetry classes. *Mathematics and Mechanics of Solids*, 9(22):1847–1865, 2017.

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

[AS18a] Houssam Abdoul-Anziz and Pierre Seppecher. Homogenization of periodic graph-based elastic structures. *Journal de l'École polytechnique, Mathématiques*, 5(259–288, 2018).

[AS18b] Houssam Abdoul-Anziz and Pierre Seppecher. Strain gradient and generalized continua obtained by homogenizing frame lattices. *Mathematics and mechanics of complex systems*, 6(3):213–250, 2018.

[ASB19] Houssam Abdoul-Anziz, Pierre Seppecher, and Cédric Bellis. 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, 2019.

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

[Bel17] M. Bellieud. Homogenization of stratined elastic composites with high contrast. *SIAM Journal on Mathematical Analysis*, 49(4):2615–2665, 2017.

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

[Bou19] C. Boutin. Homogenization methods and generalized continua in linear elasticity. In H. Altenbach and A. Öchsner, editors, *Encyclopedia of Continuum Mechanics*. Springer, Berlin, Heidelberg, 2019.

[BT94] Scott Bardenhagen and Nicolas Triantafyllidis. Derivation of higher order gradient continuum theories in 2, 3-d non-linear elasticity from periodic lattice models. *Journal of the Mechanics and Physics of Solids*, 42(1):111–139, 1994.

[CMR06] Denis Caillerie, Ayman Mourad, and Annie Raoult. Discrete homogenization in graphene sheet modeling. *Journal of Elasticity*, 84(1):33–68, 2006.

[CP12] Doina Cioranescu and Jeannine Saint Jean Paulin. *Homogenization of reticulated structures*, volume 136. Springer Science & Business Media, 2012.

[Dav13] Cesare Davini. Homogenization of linearly elastic honeycombs. *Mathematics and Mechanics of Solids*, 18(1):3–23, 2013.

[DG12] F Dos Reis and JF Ganghoffer. Construction of micropolar continua from the asymptotic homogenization of beam lattices. *Computers & Structures*, 112:354–363, 2012.

[DLSS22] Baptiste Durand, Arthur Lebée, Pierre Seppecher, and Karam Sab. Predictive strain-gradient homogenization of a pantographic material with compliant junctions. *Journal of the Mechanics and Physics of Solids*, page 104773, 2022.

[DO11] Cesare Davini and Federica Ongaro. A homogenized model for honeycomb cellular materials. *Journal of Elasticity*, 104(1):205–226, 2011.

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

[DYF+20] Bolei Deng, Siqin Yu, Antonio E Forte, Vincent Tournat, and Katia Bertoldi. Characterization, stability, and application of domain walls in flexible mechanical metamaterials. *Proceedings of the National Academy of Sciences*, 117(49):31002–31009, 2020.

[FDA10] Norman A Fleck, Vikram S Deshpande, and Michael F Ashby. Micro-architectured materials: past, present and future. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 466(2121):2495–2516, 2010.

[GLTK19] Raphaël N Glaesener, Claire Lestringant, Bastian Telgen, and Dennis M Kochmann. Continuum models for stretching- and bending-dominated periodic trusses undergoing finite deformations. *International Journal of Solids and Structures*, 171:117–134, 2019.

[HF06] R. G. Hutchinson and N. A. Fleck. The structural performance of the periodic truss. *Journal of the Mechanics and Physics of Solids*, 54:756–782, 2006.

[JS20] L. Jakabčin and P. Seppecher. On periodic homogenization of highly contrasted elastic structures. *Journal of the Mechanics and Physics of Solids*, 144:104104, 2020.

[KM04] Rajesh S Kumar and David L McDowell. Generalized continuum modeling of 2-d periodic cellular solids. *International Journal of solids and structures*, 41(26):7399–7422, 2004.

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

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

[LP16] Y. Le Floch and À Pelayo. Euler–MacLaurin formulas via differential operators. *Advances in Applied Mathematics*, 73:99–124, 2016.

[LR13] Hervé Le Dret and Annie Raoult. Homogenization of hexagonal lattices. *Networks & Heterogeneous Media*, 8(2):541, 2013.

[MA02] Andrei V Metrikine and Harm Askes. One-dimensional dynamically consistent gradient elasticity models derived from a discrete microstructure: part 1: generic formulation. *European Journal of Mechanics-A/Solids*, 21(4):555–572, 2002.

[MA06] AV Metrikine and H Askes. An isotropic dynamically consistent gradient elasticity model derived from a 2d lattice. *Philosophical Magazine*, 86(21-22):3259–3286, 2006.

[ME68] R. D. Mindlin and N. N. Eshel. On first strain-gradient theories in linear elasticity. *International Journal of Solids and Structures*, 4:109–124, 1968.

[NCH20] Hussein Nassar, Hui Chen, and Guoliang Huang. Microtwist elasticity: a continuum approach to zero modes and topological polarization in kagome lattices. *Journal of the Mechanics and Physics of Solids*, 144:104107, 2020.

[PS97] C. Pideri and P. Seppecher. A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium. *Continuum Mechanics and Thermodynamics*, 9:241–257, 1997.

[PTC22] Paul Plucinsky, Ian Tobasco, and Paolo Celli. Continuum field theory for the deformations of planar kirigami. *Bulletin of the American Physical Society*, 2022.

[RA16] Giuseppe Rosi and Nicolas Auffray. Anisotropic and dispersive wave propagation within strain-gradient framework. *Wave Motion*, 63:120–134, 2016.

[RA19] Giuseppe Rosi and Nicolas Auffray. Continuum modelling of frequency dependent acoustic beam focussing and steering in hexagonal lattices. *European Journal of Mechanics-A/Solids*, 77:103803, 2019.
[AB21] Bilen Emek Abali and Emilio Barchiesi. Additive manufacturing introduced substructure and computational determination of metamaterials parameters by means of the asymptotic homogenization. *Continuum Mechanics and Thermodynamics*, 33(4):993–1009, 2021.

[ABB10] Nicolas Auffray, Regis Bouchet, and Y Brechet. Strain gradient elastic homogenization of bidimensional cellular media. *International Journal of Solids and Structures*, 47(13):1698–1710, 2010.

[AG97] Michael F Ashby and Lorna J Gibson. Cellular solids: structure and properties. *Press Syndicate of the University of Cambridge*, Cambridge, UK, pages 175–231, 1997.

[AK017] N. Auffray, B. Kolev, and M. Olive. Handbook of bi-dimensional tensors: part I: harmonic decomposition and symmetry classes. *Mathematics and Mechanics of Solids*, 9(22):1847–1865, 2017.

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

[AS18a] Houssam Abdoul-Anziz and Pierre Seppecher. Homogenization of periodic graph-based elastic structures. *Journal de l'École polytechnique, Mathématiques*, 5(259–288, 2018).

[AS18b] Houssam Abdoul-Anziz and Pierre Seppecher. Strain gradient and generalized continua obtained by homogenizing frame lattices. *Mathematics and mechanics of complex systems*, 6(3):213–250, 2018.

[ASB19] Houssam Abdoul-Anziz, Pierre Seppecher, and Cédric Bellis. 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, 2019.

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

[Bel17] M. Bellieud. Homogenization of stratined elastic composites with high contrast. *SIAM Journal on Mathematical Analysis*, 49(4):2615–2665, 2017.

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

[Bou19] C. Boutin. Homogenization methods and generalized continua in linear elasticity. In H. Altenbach and A. Öchsner, editors, *Encyclopedia of Continuum Mechanics*. Springer, Berlin, Heidelberg, 2019.

[BT94] Scott Bardenhagen and Nicolas Triantafyllidis. Derivation of higher order gradient continuum theories in 2, 3-d non-linear elasticity from periodic lattice models. *Journal of the Mechanics and Physics of Solids*, 42(1):111–139, 1994.

[CMR06] Denis Caillerie, Ayman Mourad, and Annie Raoult. Discrete homogenization in graphene sheet modeling. *Journal of Elasticity*, 84(1):33–68, 2006.

[CP12] Doina Cioranescu and Jeannine Saint Jean Paulin. *Homogenization of reticulated structures*, volume 136. Springer Science & Business Media, 2012.

[Dav13] Cesare Davini. Homogenization of linearly elastic honeycombs. *Mathematics and Mechanics of Solids*, 18(1):3–23, 2013.

[DG12] F Dos Reis and JF Ganghoffer. Construction of micropolar continua from the asymptotic homogenization of beam lattices. *Computers & Structures*, 112:354–363, 2012.

[DLSS22] Baptiste Durand, Arthur Lebée, Pierre Seppecher, and Karam Sab. Predictive strain-gradient homogenization of a pantographic material with compliant junctions. *Journal of the Mechanics and Physics of Solids*, page 104773, 2022.

[DO11] Cesare Davini and Federica Ongaro. A homogenized model for honeycomb cellular materials. *Journal of Elasticity*, 104(1):205–226, 2011.

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

[DYF+20] Bolei Deng, Siqin Yu, Antonio E Forte, Vincent Tournat, and Katia Bertoldi. Characterization, stability, and application of domain walls in flexible mechanical metamaterials. *Proceedings of the National Academy of Sciences*, 117(49):31002–31009, 2020.

[FDA10] Norman A Fleck, Vikram S Deshpande, and Michael F Ashby. Micro-architectured materials: past, present and future. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 466(2121):2495–2516, 2010.

[GLTK19] Raphaël N Glaesener, Claire Lestringant, Bastian Telgen, and Dennis M Kochmann. Continuum models for stretching- and bending-dominated periodic trusses undergoing finite deformations. *International Journal of Solids and Structures*, 171:117–134, 2019.

[HF06] R. G. Hutchinson and N. A. Fleck. The structural performance of the periodic truss. *Journal of the Mechanics and Physics of Solids*, 54:756–782, 2006.

[JS20] L. Jakabčin and P. Seppecher. On periodic homogenization of highly contrasted elastic structures. *Journal of the Mechanics and Physics of Solids*, 144:104104, 2020.

[KM04] Rajesh S Kumar and David L McDowell. Generalized continuum modeling of 2-d periodic cellular solids. *International Journal of solids and structures*, 41(26):7399–7422, 2004.

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

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

[LP16] Y. Le Floch and À Pelayo. Euler–MacLaurin formulas via differential operators. *Advances in Applied Mathematics*, 73:99–124, 2016.

[LR13] Hervé Le Dret and Annie Raoult. Homogenization of hexagonal lattices. *Networks & Heterogeneous Media*, 8(2):541, 2013.

[MA02] Andrei V Metrikine and Harm Askes. One-dimensional dynamically consistent gradient elasticity models derived from a discrete microstructure: part 1: generic formulation. *European Journal of Mechanics-A/Solids*, 21(4):555–572, 2002.

[MA06] AV Metrikine and H Askes. An isotropic dynamically consistent gradient elasticity model derived from a 2d lattice. *Philosophical Magazine*, 86(21-22):3259–3286, 2006.

[ME68] R. D. Mindlin and N. N. Eshel. On first strain-gradient theories in linear elasticity. *International Journal of Solids and Structures*, 4:109–124, 1968.

[NCH20] Hussein Nassar, Hui Chen, and Guoliang Huang. Microtwist elasticity: a continuum approach to zero modes and topological polarization in kagome lattices. *Journal of the Mechanics and Physics of Solids*, 144:104107, 2020.

[PS97] C. Pideri and P. Seppecher. A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium. *Continuum Mechanics and Thermodynamics*, 9:241–257, 1997.

[PTC22] Paul Plucinsky, Ian Tobasco, and Paolo Celli. Continuum field theory for the deformations of planar kirigami. *Bulletin of the American Physical Society*, 2022.

[RA16] Giuseppe Rosi and Nicolas Auffray. Anisotropic and dispersive wave propagation within strain-gradient framework. *Wave Motion*, 63:120–134, 2016.

[RA19] Giuseppe Rosi and Nicolas Auffray. Continuum modelling of frequency dependent acoustic beam focussing and steering in hexagonal lattices. *European Journal of Mechanics-A/Solids*, 77:103803, 2019.
[RAPG19] Ondrej Rokoš, Maqsood M Ameen, Ron HJ Peerlings, and Mark GD Geers. Micromorphic computational homogenization for mechanical metamaterials with patterning fluctuation fields. *Journal of the Mechanics and Physics of Solids*, 123:119–137, 2019.

[RDK17] Julien Réthoré, Thi Bach Tuyet Dang, and Christine Kaltenbrunner. Anisotropic failure and size effects in periodic honeycomb materials: a gradient-elasticity approach. *Journal of the Mechanics and Physics of Solids*, 99:35–49, 2017.

[RDVB19a] G Rizzi, F Dal Corso, D Veber, and D Bigoni. Identification of second-gradient elastic materials from planar hexagonal lattices. part i: analytical derivation of equivalent constitutive tensors. *International Journal of Solids and Structures*, 176:1–18, 2019.

[RDVB19b] G Rizzi, F Dal Corso, D Veber, and D Bigoni. Identification of second-gradient elastic materials from planar hexagonal lattices. part ii: mechanical characteristics and model validation. *International Journal of Solids and Structures*, 176:19–35, 2019.

[RKD+15] Julien Réthoré, Christine Kaltenbrunner, Thi Bach Tuyet Dang, Philippe Chaudet, and Manuel Kuhn. Gradient-elasticity for honeycomb materials: validation and identification from full-field measurements. *International Journal of Solids and Structures*, 72:108–117, 2015.

[SA11] P. Seppecher, J.-J. Alibert, and F. dell'Isola. Linear elastic trusses leading to continua with exotic mechanical interactions. *Journal of Physics: Conference Series*, 319:21–25, 2011.

[San80] Évariste Sanchez-Palencia. Non-homogeneous media and vibration theory. *Lecture Note in Physics*, Springer-Verlag, 320:57–65, 1980.

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

[SCO+22] Angkur Jyoti Dipanka Shaikeea, Huachen Cui, Mark O'Masta, Xiaoyu Rayne Zheng, and Vikram Sudhir Deshpande. The toughness of mechanical metamaterials. *Nature Materials*, 21(3):297–304, 2022.

[TB93] Nicolas Triantafyllidis and S Bardenhagen. On higher order gradient continuum theories in 1-d nonlinear elasticity. derivation from and comparison to the corresponding discrete models. *Journal of Elasticity*, 33(3):259–293, 1993.

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

[VDP14] Andrea Vigliotti, Vikram S Deshpande, and Damiano Pasini. Non linear constitutive models for lattice materials. *Journal of the Mechanics and Physics of Solids*, 64:44–60, 2014.

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

# APPENDIX A. NOTATIONS AND TOOLBOX

## A.1. Tensors

The dimension of the Euclidean space is denoted as $d$. We address both elastic lattices in the plane ($d=2$, see the illustration examples in §6.1–7.3 and §7.5) and in the three-dimensional space ($d=3$, §7.4). 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 $\mathbf{X} \in \mathbb{R}^d$, see Figure 2.3.

We denote as $\mathbf{T}^{(n_1, n_2, \dots, n_p)}$ the tensor space $\mathbf{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 $\mathbf{G}$ of rank $p$ and dimensions $n_1 \times n_2 \times \dots \times n_p$. In particular, $I_r \in \mathbf{T}^{(r, r)}$ denotes the identity matrix in dimension $r$ and $\mathbf{0}_{n_1 \times \dots \times n_p} \in \mathbf{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 $\mathbf{G} \in \mathbf{T}^{(n_1, n_2, \dots, n_p)}$ and $\mathbf{G}' \in \mathbf{T}^{(n'_1, n'_2, \dots, n'_p)}$, we denote as

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

* $\mathbf{G} : \mathbf{G}' \in \mathbf{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$),

* $\mathbf{G} : \mathbf{G}' \in \mathbf{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.

The contracted tensors are given by

$$ \begin{align} (G \cdot G')_{i_1 \dots i_{p-1} i'_1 \dots i'_p} &= G_{i_1 \dots i_{p-1} j} G'_{j i'_1 \dots i'_p} \tag{A.1} \\ (G : G')_{i_1 \dots i_{p-2} i'_1 \dots i'_p} &= G_{i_1 \dots i_{p-1} j} G'_{j i'_1 \dots i'_p} \notag \\ (G : G')_{i_1 \dots i_{p-3} i'_1 \dots i'_p} &= G_{i_1 \dots i_{p-1} j k} G'_{j k i'_1 \dots i'_p} \notag \end{align} $$

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. Note the ordering of the contracted indices $j, k, l$, etc. in the right-hand sides.

The action of a matrix $\mathbf{G}$ on a vector $\mathbf{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 $\mathbf{G} \cdot \mathbf{v}$, with a dot.

The outer product of two tensors $\mathbf{T}$ and $\mathbf{T}'$ is denoted as $\mathbf{T} \otimes \mathbf{T}'$. In particular, the outer product of two vectors is denoted as $\mathbf{v} \otimes \mathbf{v}'$. Vector transposition is never used nor is even meaningful.

Given a tensor $\mathbf{G} \in \mathbf{T}^{(n_1, n_2, \dots, n_p)}$, and a permutation $(\sigma_1, \dots, \sigma_p)$ of the levels $(1, \dots, p)$ of the tensor, we denote as $\mathbf{G}^{\mathbf{T}^{(\sigma_1, \dots, \sigma_p)}}$ the generalized transpose of $\mathbf{G}$, such that the level $i$ in the original tensor becomes level $\sigma_i$ in the transpose:

$$ (G^{\mathbf{T}^{(\sigma_1, \dots, \sigma_p)}})_{i_1 \dots i_p} = G_{i_{\sigma_1} \dots i_{\sigma_p}} \tag{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 $(G^{\mathbf{T}^{(1342)}})_{ijkl} = G_{iklj}$.
[RAPG19] Ondrej Rokoš, Maqsood M Ameen, Ron HJ Peerlings, and Mark GD Geers. Micromorphic computational homogenization for mechanical metamaterials with patterning fluctuation fields. *Journal of the Mechanics and Physics of Solids*, 123:119–137, 2019.

[RDK17] Julien Réthoré, Thi Bach Tuyet Dang, and Christine Kaltenbrunner. Anisotropic failure and size effects in periodic honeycomb materials: a gradient-elasticity approach. *Journal of the Mechanics and Physics of Solids*, 99:35–49, 2017.

[RDVB19a] G Rizzi, F Dal Corso, D Veber, and D Bigoni. Identification of second-gradient elastic materials from planar hexagonal lattices. part i: analytical derivation of equivalent constitutive tensors. *International Journal of Solids and Structures*, 176:1–18, 2019.

[RDVB19b] G Rizzi, F Dal Corso, D Veber, and D Bigoni. Identification of second-gradient elastic materials from planar hexagonal lattices. part ii: mechanical characteristics and model validation. *International Journal of Solids and Structures*, 176:19–35, 2019.

[RKD+15] Julien Réthoré, Christine Kaltenbrunner, Thi Bach Tuyet Dang, Philippe Chaudet, and Manuel Kuhn. Gradient-elasticity for honeycomb materials: validation and identification from full-field measurements. *International Journal of Solids and Structures*, 72:108–117, 2015.

[Ad11] P. Seppecher, J.-J. Alibert, and F. dell'Isola. Linear elastic trusses leading to continua with exotic mechanical interactions. *Journal of Physics: Conference Series*, 319:21–25, 2011.

[San80] Évariste Sanchez-Palencia. Non-homogeneous media and vibration theory. *Lecture Note in Physics, Springer-Verlag*, 320:57–65, 1980.

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

[SCO+22] Angkur Jyoti Dipanka Shaikeea, Huachen Cui, Mark O'Masta, Xiaoyu Rayne Zheng, and Vikram Sudhir Deshpande. The toughness of mechanical metamaterials. *Nature Materials*, 21(3):297–304, 2022.

[TB93] Nicolas Triantafyllidis and S Bardenhagen. On higher order gradient continuum theories in 1-d nonlinear elasticity. derivation from and comparison to the corresponding discrete models. *Journal of Elasticity*, 33(3):259–293, 1993.

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

[VDP14] Andrea Vigliotti, Vikram S Deshpande, and Damiano Pasini. Non linear constitutive models for lattice materials. *Journal of the Mechanics and Physics of Solids*, 64:44–60, 2014.

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

# APPENDIX A. NOTATIONS AND TOOLBOX

## A.1. Tensors

The dimension of the Euclidean space is denoted as $d$. We address both elastic lattices in the plane ($d=2$, see the illustration examples in §6.1–7.3 and §7.5) and in the three-dimensional space ($d=3$, §7.4). 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 $\mathbf{X} \in \mathbb{R}^d$, see Figure 2.3.

We denote as $\mathbf{T}^{(n_1, n_2, \dots, n_p)}$ the tensor space $\mathbf{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 $\mathbf{G}$ of rank $p$ and dimensions $n_1 \times n_2 \times \dots \times n_p$. In particular, $I_r \in \mathbf{T}^{(r, r)}$ denotes the identity matrix in dimension $r$ and $\mathbf{0}_{n_1 \times \dots \times n_p} \in \mathbf{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 $\mathbf{G} \in \mathbf{T}^{(n_1, n_2, \dots, n_p)}$ and $\mathbf{G}' \in \mathbf{T}^{(n'_1, n'_2, \dots, n'_p)}$, we denote as

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

* $\mathbf{G} : \mathbf{G}' \in \mathbf{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$),

* $\mathbf{G} : \mathbf{G}' \in \mathbf{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.

The contracted tensors are given by

$$ 
\begin{aligned}
(G \cdot G')_{i_1 \dots i_{p-1} i'_2 \dots i'_p} &= G_{i_1 \dots i_{p-1} j} G'_{j i'_2 \dots i'_p} \\
(G : G')_{i_1 \dots i_{p-2} i'_3 \dots i'_p} &= G_{i_1 \dots i_{p-1} j} G'_{j k i'_2 \dots i'_p} \\
(G : G')_{i_1 \dots i_{p-3} i'_4 \dots i'_p} &= G_{i_1 \dots i_{p-1} j k l} G'_{j k i'_2 \dots i'_p}
\end{aligned} 
\quad (A.1) 
$$

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. Note the ordering of the contracted indices $j, k, l$, etc. in the right-hand sides.

The action of a matrix $\mathbf{G}$ on a vector $\mathbf{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 $\mathbf{G} \cdot \mathbf{v}$, with a dot.

The outer product of two tensors $\mathbf{T}$ and $\mathbf{T}'$ is denoted as $\mathbf{T} \otimes \mathbf{T}'$. In particular, the outer product of two vectors is denoted as $\mathbf{v} \otimes \mathbf{v}'$. Vector transposition is never used nor is even meaningful.

Given a tensor $\mathbf{G} \in \mathbf{T}^{(n_1, n_2, \dots, n_p)}$, and a permutation $(\sigma_1, \dots, \sigma_p)$ of the levels $(1, \dots, p)$ of the tensor, we denote as $\mathbf{G}^{\mathbf{T}^{(\sigma_1, \dots, \sigma_p)}}$ the generalized transpose of $\mathbf{G}$, such that the level $i$ in the original tensor becomes level $\sigma_i$ in the transpose:

$$ (G^{\mathbf{T}^{(\sigma_1, \dots, \sigma_p)}})_{i_1 \dots i_p} = G_{i_{\sigma_1} \dots i_{\sigma_p}} \quad (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 $(G^{T^{1342}})_{ijkl} = G_{iklj}$.
Transposing enables us to reorder the indices of a tensor in any desired order. Suppose that we wish to rewrite an expression $G_{ijkl}$ as the component $G'_{ijkl}$ of another tensor whose indices are in alphabetical order: $G'$ is a transpose of $G$, the permutation being found by noting that the levels (1, 2, 3, 4) in the original tensor $G$, 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 $G'$. This yields

$$G_{ijkl} = (G^{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 $G_{ijkl} G'_{ijkl} = G^{T_{1324l}} G'$. The usual transpose $G^T$ of a matrix $G \in T^{(n_1,n_2)}$ is a particular case of the generalized transpose, $G^T = G^{T_{21}}$.

Given a tensor $G(Y) \in T^{(n_1,n_2,\dots,n_p)}$ taking an argument $Y \in \mathbb{R}^s$, we denote as $\nabla G(Y) \in T^{(n_1,n_2,\dots,n_p,s)}$ its gradient,

$$\nabla G_{i_1 \dots i_p} = \frac{\partial G_{i_1 \dots i_p}}{\partial Y_j}.$$

By a standard convention, the index $j$ corresponding to differentiation is therefore the last one even though the symbol $\nabla$ appears first in $\nabla G$.

## A.2. Infinitesimal rotations

In terms of the dimension $d$ of the Euclidean space, we define $n_r = \frac{d(d-1)}{2}$ as the dimension of infinitesimal rotations: $n_r = 1$ in dimension $d=2$ and $n_r = 3$ in dimension $d=3$. An infinitesimal rotation is represented by a vector $\theta \in \mathbb{R}^{n_r}$, also known as a pseudo-vector. By convention, its indices are written as

$$\theta = \begin{cases} (\theta_3) & (d=2) \\ (\theta_1, \theta_2, \theta_3) & (d=3), \end{cases} \qquad (A.4)$$

i.e., the unique index $k$ of $\theta_k$ is taken to be $k=3$ and not $k=1$ when $d=2$, as $\theta$ then lives on the line perpendicular to the Euclidean plane $\mathbb{R}^2$.

The pseudo-vector $\theta$ is uniquely associated with an antisymmetric matrix $\hat{\theta} \in T^{(d,d)}$ by

$$\hat{\theta} = -\mathcal{N} \cdot \theta \Leftrightarrow \theta = -\frac{1}{2} \hat{\theta} : \mathcal{N}, \qquad (A.5)$$

where $\mathcal{N} \in T^{(d,d,n_r)}$ is a constant rank-3 tensor that is antisymmetric with respect to its first pair of indices,

$$\mathcal{N}_{ijk} = \begin{cases} \tilde{\mathcal{A}}_{ij3} & (d=2) \\ \tilde{\mathcal{A}}_{ijk} & (d=3). \end{cases} \qquad (A.6)$$

Here $\tilde{\mathcal{A}}_{ijk}$ is the alternating symbol in dimension 3: $\tilde{\mathcal{A}}_{ijk} \in \{-1, 0, +1\}$ is 0 if at least two of the indices $(i, j, k)$ are identical, or the signature $\pm 1$ of the permutation if $(i, j, k)$ is a permutation of $(1, 2, 3)$. In dimension $d=2$, the last index $k$ of $\mathcal{N}_{ijk}$ in (A.6) is restricted to $k=3$ as in (A.4), which warrants that the contraction appearing in (A.5)$_1$, as well as the right-most indices on both sides of (A.5)$_2$, are consistent.

The equivalence in (A.5) is established in Appendix B.2.

When represented using an antisymmetric matrix $\hat{\theta}$, the infinitesimal rotation maps a vector $x \in \mathbb{R}^d$ to $x' = (I_d + \hat{\theta}) \cdot x$. When represented as a pseudo-vector $\theta$, it maps $x$ to $x' = x + \theta_3 (\mathcal{A}^T \cdot x)$ when $d=2$ (with $\mathcal{A}$ as the 2D Levi-Civita symbol, defined below Equation (7.11)), or $x' = x + \theta \times x$ when $d=3$ (with $\times$ as the 3D cross-product). Equation (A.5) warrants that the result is the same.

## A.3. Toolbox for a strain-gradient elasticity

In Appendix B, we show that the higher-order gradients $\nabla^k u$ and $\nabla^k \gamma$ of the displacement and rotation (for $k \ge 1$) can be expressed in terms of $\gamma$ and of the gradients $\nabla^j \tilde{\epsilon}$ (with $j \ge 0$) as

$$\begin{aligned} \nabla \gamma(\mathbf{X}) &= \mathcal{G}' : \nabla \tilde{\epsilon}(\mathbf{X}) \\ \nabla^2 \gamma(\mathbf{X}) &= \mathcal{G}'' : \nabla^2 \tilde{\epsilon}(\mathbf{X}) \\ \nabla u(\mathbf{X}) &= \mathcal{V} \tilde{\epsilon}' : \tilde{\epsilon}(\mathbf{X}) - \mathcal{N}' : \gamma(\mathbf{X}) \\ \nabla^2 u(\mathbf{X}) &= \mathcal{V} \tilde{\epsilon}'' : \nabla \tilde{\epsilon}(\mathbf{X}) \\ \nabla^3 u(\mathbf{X}) &= \mathcal{V} \tilde{\epsilon}''' : \nabla^2 \tilde{\epsilon}(\mathbf{X}). \end{aligned} \qquad (A.7)$$

Note that no gradient of the rotation appears in the right-hand sides above. The constant tensors $\mathcal{G}' \in T^{(n_r,d,n_r,d)}$, $\mathcal{G}'' \in T^{(n_r,d,d,n_r,d)}$, $\mathcal{V} \tilde{\epsilon}' \in T^{(d,d,n_r)}$, $\mathcal{V} \tilde{\epsilon}'' \in T^{(d,d,d,n_r,d)}$ and $\mathcal{V} \tilde{\epsilon}''' \in T^{(d,d,d,d,n_r,d)}$ are defined in Equations (B.9), (B.11), (B.13) and (B.15) in the Appendix. They are symmetric with respect to any pair of indices associated with second or third gradients, as specified in Appendix B. They depend only on the dimension $d$ of the Euclidean space.

Equation (A.7) is a reformulation of a result by [ME68] who derived expressions for $\nabla \gamma$, $\nabla^2 u$ and their successive gradients in terms of the strain $\tilde{\epsilon}$ and its gradients.

Equation (A.7) yields the Taylor expansion of the macroscopic fields $u$ and $\gamma$ about a point $\mathbf{X}_c$ as

$$\begin{aligned} \mathbf{u}(\mathbf{X}) = \mathbf{u}(\mathbf{X}_c) &+ [\mathcal{V} \tilde{\epsilon}' : \tilde{\epsilon}(\mathbf{X}_c) - \mathcal{N}' : \gamma(\mathbf{X}_c)] \cdot (\mathbf{X} - \mathbf{X}_c) + [\mathcal{V} \tilde{\epsilon}'' : \nabla \tilde{\epsilon}(\mathbf{X}_c)] \cdot \frac{(\mathbf{X} - \mathbf{X}_c)^{\otimes 2}}{2} \\ &+ [\mathcal{V} \tilde{\epsilon}''' : \nabla^2 \tilde{\epsilon}(\mathbf{X})] \cdot \frac{(\mathbf{X} - \mathbf{X}_c)^{\otimes 3}}{6} + \dots \qquad (A.8) \\ \gamma(\mathbf{X}) = \gamma(\mathbf{X}_c) &+ [\mathcal{G}' : \nabla \tilde{\epsilon}(\mathbf{X})] \cdot (\mathbf{X} - \mathbf{X}_c) + [\mathcal{G}'' : \nabla^2 \tilde{\epsilon}(\mathbf{X})] \cdot \frac{(\mathbf{X} - \mathbf{X}_c)^{\otimes 2}}{2} + \dots \end{aligned}$$
Transposing enables us to reorder the indices of a tensor in any desired order. Suppose that we wish to rewrite an expression $G_{ijkl}$ as the component $G'_{ijkl}$ of another tensor whose indices are in alphabetical order: $G'$ is a transpose of $G$, the permutation being found by noting that the levels (1, 2, 3, 4) in the original tensor $G$, 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 $G'$. This yields

$$G_{ijkl} = (G^{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 $G_{ijkl}G'_{ijkl} = G^{T_{1324l}}G'$. The usual transpose $G^T$ of a matrix $G \in T^{(n_1,n_2)}$ is a particular case of the generalized transpose, $G^T = G^{T_{21}}$.

Given a tensor $G(Y) \in T^{(n_1,n_2,\dots,n_p)}$ taking an argument $Y \in \mathbb{R}^s$, we denote as $\nabla G(Y) \in T^{(n_1,n_2,\dots,n_p,s)}$ its gradient,

$$\nabla G_{i_1 \dots i_p} = \frac{\partial G_{i_1 \dots i_p}}{\partial Y_j}.$$

By a standard convention, the index $j$ corresponding to differentiation is therefore the last one even though the symbol $\nabla$ appears first in $\nabla G$.

## A.2. Infinitesimal rotations

In terms of the dimension $d$ of the Euclidean space, we define $n_r = \frac{d(d-1)}{2}$ as the dimension of infinitesimal rotations: $n_r = 1$ in dimension $d=2$ and $n_r = 3$ in dimension $d=3$. An infinitesimal rotation is represented by a vector $\theta \in \mathbb{R}^{n_r}$, also known as a pseudo-vector. By convention, its indices are written as

$$\theta = \begin{cases} (\theta_3) & (d=2) \\ (\theta_1, \theta_2, \theta_3) & (d=3), \end{cases} \qquad (A.4)$$

i.e., the unique index $k$ of $\theta_k$ is taken to be $k=3$ and not $k=1$ when $d=2$, as $\theta$ then lives on the line perpendicular to the Euclidean plane $\mathbb{R}^2$.

The pseudo-vector $\theta$ is uniquely associated with an antisymmetric matrix $\hat{\theta} \in T^{(d,d)}$ by

$$\hat{\theta} = -\mathcal{N} \cdot \theta \Leftrightarrow \theta = -\frac{1}{2} \hat{\theta} : \mathcal{N}, \qquad (A.5)$$

where $\mathcal{N} \in T^{(d,d,n_r)}$ is a constant rank-3 tensor that is antisymmetric with respect to its first pair of indices,

$$\mathcal{N}_{ijk} = \begin{cases} \tilde{\mathcal{A}}_{ij3} & (d=2) \\ \tilde{\mathcal{A}}_{ijk} & (d=3). \end{cases} \qquad (A.6)$$

Here $\tilde{\mathcal{A}}_{ijk}$ is the alternating symbol in dimension 3: $\tilde{\mathcal{A}}_{ijk} \in \{-1, 0, +1\}$ is 0 if at least two of the indices $(i, j, k)$ are identical, or the signature $\pm 1$ of the permutation if $(i, j, k)$ is a permutation of $(1, 2, 3)$. In dimension $d=2$, the last index $k$ of $\mathcal{N}_{ijk}$ in (A.6) is restricted to $k=3$ as in (A.4), which warrants that the contraction appearing in (A.5)$_1$, as well as the right-most indices on both sides of (A.5)$_2$, are consistent.

The equivalence in (A.5) is established in Appendix B.2.

When represented using an antisymmetric matrix $\hat{\theta}$, the infinitesimal rotation maps a vector $\mathbf{x} \in \mathbb{R}^d$ to $\mathbf{x}' = (I_d + \hat{\theta}) \cdot \mathbf{x}$. When represented as a pseudo-vector $\theta$, it maps $\mathbf{x}$ to $\mathbf{x}' = \mathbf{x} + \theta_3 (\mathcal{A}^T \cdot \mathbf{x})$ when $d=2$ (with $\mathcal{A}$ as the 2D Levi-Civita symbol, defined below Equation (7.11)), or $\mathbf{x}' = \mathbf{x} + \theta \times \mathbf{x}$ when $d=3$ (with $\times$ as the 3D cross-product). Equation (A.5) warrants that the result is the same.

## A.3. Toolbox for a strain-gradient elasticity

In Appendix B, we show that the higher-order gradients $\nabla^k \mathbf{u}$ and $\nabla^k \gamma$ of the displacement and rotation (for $k \ge 1$) can be expressed in terms of $\gamma$ and of the gradients $\nabla^j \tilde{\mathbf{\epsilon}}$ (with $j \ge 0$) as

$$\begin{aligned} \nabla \gamma(\mathbf{X}) &= \mathcal{G}' : \nabla \tilde{\mathbf{\epsilon}}(\mathbf{X}) \\ \nabla^2 \gamma(\mathbf{X}) &= \mathcal{G}'' : \nabla^2 \tilde{\mathbf{\epsilon}}(\mathbf{X}) \\ \nabla \mathbf{u}(\mathbf{X}) &= \mathcal{V} \tilde{\mathbf{\epsilon}}' : \nabla \tilde{\mathbf{\epsilon}}(\mathbf{X}) - \mathcal{N} : \gamma(\mathbf{X}) \\ \nabla^2 \mathbf{u}(\mathbf{X}) &= \mathcal{V} \tilde{\mathbf{\epsilon}}'' : \nabla \tilde{\mathbf{\epsilon}}(\mathbf{X}) \\ \nabla^3 \mathbf{u}(\mathbf{X}) &= \mathcal{V} \tilde{\mathbf{\epsilon}}''' : \nabla^2 \tilde{\mathbf{\epsilon}}(\mathbf{X}). \end{aligned} \qquad (A.7)$$

Note that no gradient of the rotation appears in the right-hand sides above. The constant tensors $\mathcal{G}' \in T^{(n_r,d,n_c,d)}$, $\mathcal{G}'' \in T^{(n_r,d,d,n_c,d)}$, $\mathcal{V} \tilde{\mathbf{\epsilon}}' \in T^{(d,d,n_r)}$, $\mathcal{V} \tilde{\mathbf{\epsilon}}'' \in T^{(d,d,d,n_r,d)}$ and $\mathcal{V} \tilde{\mathbf{\epsilon}}''' \in T^{(d,d,d,d,n_r,d)}$ are defined in Equations (B.9), (B.11), (B.13) and (B.15) in the Appendix. They are symmetric with respect to any pair of indices associated with second or third gradients, as specified in Appendix B. They depend only on the dimension $d$ of the Euclidean space.

Equation (A.7) is a reformulation of a result by [ME68] who derived expressions for $\nabla \gamma$, $\nabla^2 \mathbf{u}$ and their successive gradients in terms of the strain $\tilde{\mathbf{\epsilon}}$ and its gradients.

Equation (A.7) yields the Taylor expansion of the macroscopic fields $\mathbf{u}$ and $\gamma$ about a point $\mathbf{X}_c$ as

$$\begin{aligned} \mathbf{u}(\mathbf{X}) = \mathbf{u}(\mathbf{X}_c) &+ [\mathcal{V} \tilde{\mathbf{\epsilon}}' : \nabla \tilde{\mathbf{\epsilon}}(\mathbf{X}_c) - \mathcal{N} : \gamma(\mathbf{X}_c)] \cdot (\mathbf{X} - \mathbf{X}_c) + [\mathcal{V} \tilde{\mathbf{\epsilon}}'' : \nabla \tilde{\mathbf{\epsilon}}(\mathbf{X}_c)] \cdot \frac{(\mathbf{X} - \mathbf{X}_c)^{\otimes 2}}{2} \\ &+ [\mathcal{V} \tilde{\mathbf{\epsilon}}''' : \nabla^2 \tilde{\mathbf{\epsilon}}(\mathbf{X})] \cdot \frac{(\mathbf{X} - \mathbf{X}_c)^{\otimes 3}}{6} + \dots \\ \mathbf{\gamma}(\mathbf{X}) = \mathbf{\gamma}(\mathbf{X}_c) &+ [\mathcal{G}' : \nabla \tilde{\mathbf{\epsilon}}(\mathbf{X})] \cdot (\mathbf{X} - \mathbf{X}_c) + [\mathcal{G}'' : \nabla^2 \tilde{\mathbf{\epsilon}}(\mathbf{X})] \cdot \frac{(\mathbf{X} - \mathbf{X}_c)^{\otimes 2}}{2} + \dots \end{aligned} \qquad (A.8)$$
where $x^{\otimes 2} = x \otimes x$, $x^{\otimes 3} = x \otimes x \otimes x$. In Section 3, Equation (A.8) is used to write the strain energy of a particular element, which is expressed in terms of the degrees of freedom at the adjacent nodes, into a Taylor expansion at the center of the element.

# APPENDIX B. PROOF OF THE GEOMETRIC IDENTITIES

## B.1. Mandel representation $\check{\epsilon}$ of the macroscopic strain

We introduce the following constant tensor $\mathcal{M} \in T^{(d,d,n_c)}$ depending only on the dimension $d$ of the Euclidean space,

$$ \mathcal{M} = \sum_{i=1}^{d} [\mathbf{i}_i^d \otimes \mathbf{i}_i^d] \otimes \mathbf{i}_i^{n_c} + \sum_{1 \le i < j \le d} \left[ \frac{[\mathbf{i}_i^d \otimes \mathbf{i}_j^d + \mathbf{i}_j^d \otimes \mathbf{i}_i^d]}{\sqrt{2}} \right] \otimes \mathbf{i}_{i+(j-i)}^{n_c} (d - \frac{j-i-1}{2}), \quad (B.1) $$

where $\mathbf{i}_a^b \in \mathbb{R}^b$ is the discrete Dirac vector in dimension $b$, whose components are all zero except for the $a$-th one which is equal to 1,

$$ (\mathbf{i}_a^b)_j = \begin{cases} 0 & \text{if } j \neq a, \\ 1 & \text{if } j = a, \end{cases} \quad \text{with } 1 \le j \le b. \quad (B.2) $$

The rationale behind the definition of $\mathcal{M}$ is that the quantities in square brackets form an orthonormal basis of the space of symmetric $d \times d$ matrices, and that the indices $i$ (first term) and $i + (j-i) (d - \frac{j-i-1}{2})$ (second term) provide a sequential numbering of the tensors forming this orthonormal basis, from 1 to $n_c$. The orthonormal basis of symmetric $d \times d$ matrices can be therefore be extracted as $(\mathcal{M} \cdot \mathbf{i}_k^{n_c})_{1 \le k \le n_c}$.

The following two identities are a consequence of the orthonormal character of the basis appearing in square brackets in (B.1),

$$ \mathcal{M}^{T_{231}} : \mathcal{M} = \mathbf{I}_{n_c}, \quad \mathcal{M} \cdot \mathcal{M}^{T_{231}} = ((\mathbf{I}_d \otimes \mathbf{I}_d)^{T_{1324}}) S_{12} \circ S_{34}, \quad (B.3) $$

where we recall that $\circ$ stands for the composition of symmetrization operations, see (5.1–2.3).

The Mandel representation $\check{\epsilon} \in \mathbb{R}^{n_c}$ of the strain tensor $\epsilon$ introduced in (3.7) is nothing but the decomposition of $\epsilon$ in this orthonormal basis,

$$ \check{\epsilon} = \epsilon : \mathcal{M}. \quad (B.4) $$

Conversely, the symmetric strain tensor $\epsilon$ can be reconstructed from the vector $\check{\epsilon}$ by

$$ \epsilon = \mathcal{M} \cdot \check{\epsilon}, \quad (B.5) $$

as can be shown using (B.3).

## B.2. Properties of the tensor $\mathcal{N}$

The tensor $\mathcal{N}'$ introduced in (A.6) is similar to $\mathcal{M}$ as it delivers an orthogonal (but not orthonormal) basis of antisymmetric tensors in $T^{(d,d)}$ indexed by $k$ in the form $(\mathcal{N}' \cdot \mathbf{i}_k^{n_c})$ where

* in dimension $d=2$, the only possible rotation index $k$ is limited to $k \in \{3\}$ and we set $\mathbf{i}_3^1 = (1)$,

* in dimension $d=3$, $1 \le k \le 3$ and the usual definition of $\mathbf{i}_k^{n_c}$ in (B.2) applies.

In view of this and of Equation (A.5)$_1$, $\theta$ is nothing but the components of $\hat{\theta}$ in this basis, up to a sign. The reciprocal formula (A.5)$_2$ follows from the identity

$$ \frac{1}{2} \mathcal{N}^{T_{231}} : \mathcal{N} = \mathbf{I}_{n_c}, \quad (B.6) $$

which expresses itself the orthogonal character of the underlying basis.

## B.3. Elimination of the rotation gradient

In this section, we provide a detailed justification of the identities (A.7), as well as expressions of the geometric tensors $\mathcal{G}', \mathcal{G}'', \mathcal{U}', \mathcal{U}'', \mathcal{U}'''$ appearing in the right-hand side. These tensors depend on the dimension $d$ of the Euclidean space only.

Combining the decomposition of $\nabla \mathbf{u} = \epsilon + \hat{\gamma}$ into a symmetric and antisymmetric part with (B.5) and (A.5)$_1$, we have $\nabla \mathbf{u}(\mathbf{X}) = \mathcal{M} \cdot \check{\epsilon}(\mathbf{X}) - \mathcal{N}' \cdot \gamma(\mathbf{X})$ which is nothing but (A.7)$_3$ when we identify

$$ \mathcal{U}' = \mathcal{M}. \quad (B.7) $$

Next, we repeat the classical argument yielding the rotation gradient $\nabla \gamma$ in terms of the strain gradient $\nabla \check{\epsilon}$. The rotation in (3.8)$_2$ can be written in components as $\gamma_a = -\frac{1}{2} u_{c,d} \mathcal{N}_{cda} = +\frac{1}{2} \mathcal{N}_{dca} u_{c,d}$, where a comma in subscript denotes differentiation. The gradient of rotation thus takes the form $\gamma_{a,b} = \frac{1}{2} \mathcal{N}_{dca} u_{c,db}$. Since the second gradient $u_{c,db}$ is symmetric with respect to the indices $(b, d)$, this can be rewritten as $\gamma_{a,b} = \frac{1}{2} \mathcal{N}_{dca} (u_{c,b})_d = \frac{1}{2} \mathcal{N}_{dca} (u_{c,b} + u_{b,c})_d - \frac{1}{2} \mathcal{N}_{dca} u_{b,cd}$. The last term is zero since it is a contraction over the indices $(c, d)$ of an antisymmetric tensor $\mathcal{N}_{dca}$ and a symmetric tensor $u_{b,cd}$. This yields $\gamma_{a,b} = \mathcal{N}_{dca} \epsilon_{cb,d} = \mathcal{N}_{dca} \mathcal{M}_{cbe} \check{\epsilon}_{e,d}$ by (B.5), which we can rewrite in tensor form as

$$ \nabla \gamma(\mathbf{X}) = \mathcal{G}' : \nabla \check{\epsilon}(\mathbf{X}), \quad (B.8) $$
where $x^{\otimes 2} = x \otimes x$, $x^{\otimes 3} = x \otimes x \otimes x$. In Section 3, Equation (A.8) is used to write the strain energy of a particular element, which is expressed in terms of the degrees of freedom at the adjacent nodes, into a Taylor expansion at the center of the element.

# APPENDIX B. PROOF OF THE GEOMETRIC IDENTITIES

## B.1. Mandel representation $\tilde{\epsilon}$ of the macroscopic strain

We introduce the following constant tensor $\mathcal{M} \in T^{(d,d,n_c)}$ depending only on the dimension $d$ of the Euclidean space,

$$ \mathcal{M} = \sum_{i=1}^{d} [\mathbf{i}_i^d \otimes \mathbf{i}_i^d] \otimes \mathbf{i}_i^{n_c} + \sum_{1 \le i < j \le d} \left[ \frac{[\mathbf{i}_i^d \otimes \mathbf{i}_j^d + \mathbf{i}_j^d \otimes \mathbf{i}_i^d]}{\sqrt{2}} \right] \otimes \mathbf{i}_{i+(j-i)}^{n_c} (d - \frac{j-i-1}{2}), \quad (B.1) $$

where $\mathbf{i}_a^b \in \mathbb{R}^b$ is the discrete Dirac vector in dimension $b$, whose components are all zero except for the $a$-th one which is equal to 1,

$$ (\mathbf{i}_a^b)_j = \begin{cases} 0 & \text{if } j \neq a, \\ 1 & \text{if } j = a, \end{cases} \quad \text{with } 1 \le j \le b. \quad (B.2) $$

The rationale behind the definition of $\mathcal{M}$ is that the quantities in square brackets form an orthonormal basis of the space of symmetric $d \times d$ matrices, and that the indices $i$ (first term) and $i + (j-i) (d - \frac{j-i-1}{2})$ (second term) provide a sequential numbering of the tensors forming this orthonormal basis, from 1 to $n_c$. The orthonormal basis of symmetric $d \times d$ matrices can be therefore be extracted as $(\mathcal{M} \cdot \mathbf{i}_k^{n_c})_{1 \le k \le n_c}$.

The following two identities are a consequence of the orthonormal character of the basis appearing in square brackets in (B.1),

$$ \mathcal{M}^{T_{231}} : \mathcal{M} = I_{n_c} \qquad \mathcal{M} \cdot \mathcal{M}^{T_{231}} = ((I_d \otimes I_d)^{T_{1324}}) S_{12} \circ S_{34}, \quad (B.3) $$

where we recall that $\circ$ stands for the composition of symmetrization operations, see (5.1-2.3).

The Mandel representation $\tilde{\epsilon} \in \mathbb{R}^{n_c}$ of the strain tensor $\epsilon$ introduced in (3.7) is nothing but the decomposition of $\epsilon$ in this orthonormal basis,

$$ \tilde{\epsilon} = \epsilon : \mathcal{M}. \quad (B.4) $$

Conversely, the symmetric strain tensor $\epsilon$ can be reconstructed from the vector $\tilde{\epsilon}$ by

$$ \epsilon = \mathcal{M} \cdot \tilde{\epsilon}, \quad (B.5) $$

as can be shown using (B.3).

## B.2. Properties of the tensor $\mathcal{N}$

The tensor $\mathcal{N}'$ introduced in (A.6) is similar to $\mathcal{M}$ as it delivers an orthogonal (but not orthonormal) basis of antisymmetric tensors in $T^{(d,d)}$ indexed by $k$ in the form $(\mathcal{N}' \cdot \mathbf{i}_k^{n_c})$ where

* in dimension $d=2$, the only possible rotation index $k$ is limited to $k \in \{3\}$ and we set $\mathbf{i}_3^1 = (1)$,

* in dimension $d=3$, $1 \le k \le 3$ and the usual definition of $\mathbf{i}_k^{n_c}$ in (B.2) applies.

In view of this and of Equation (A.5)$_1$, $\theta$ is nothing but the components of $\hat{\theta}$ in this basis, up to a sign. The reciprocal formula (A.5)$_2$ follows from the identity

$$ \frac{1}{2} \mathcal{N}^{T_{231}} : \mathcal{N} = I_{n_c}, \quad (B.6) $$

which expresses itself the orthogonal character of the underlying basis.

## B.3. Elimination of the rotation gradient

In this section, we provide a detailed justification of the identities (A.7), as well as expressions of the geometric tensors $\mathcal{G}', \mathcal{G}'', \mathcal{U}', \mathcal{U}'', \mathcal{U}'''$ appearing in the right-hand side. These tensors depend on the dimension $d$ of the Euclidean space only.

Combining the decomposition of $\nabla \mathbf{u} = \epsilon + \hat{\gamma}$ into a symmetric and antisymmetric part with (B.5) and (A.5)$_1$, we have $\nabla \mathbf{u}(\mathbf{X}) = \mathcal{M} \cdot \hat{\mathbf{\epsilon}}(\mathbf{X}) - \mathcal{N}' \cdot \hat{\gamma}(\mathbf{X})$ which is nothing but (A.7)$_3$ when we identify

$$ \mathcal{U}' = \mathcal{M}. \quad (B.7) $$

Next, we repeat the classical argument yielding the rotation gradient $\nabla \hat{\gamma}$ in terms of the strain gradient $\nabla \tilde{\epsilon}$. The rotation in (3.8)$_2$ can be written in components as $\gamma_a = -\frac{1}{2} u_{c,d} \mathcal{N}_{cda} = +\frac{1}{2} \mathcal{N}_{dca} u_{c,d}$, where a comma in subscript denotes differentiation. The gradient of rotation thus takes the form $\gamma_{a,b} = \frac{1}{2} \mathcal{N}_{dca} u_{c,db}$. Since the second gradient $u_{c,db}$ is symmetric with respect to the indices $(b, d)$, this can be rewritten as $\gamma_{a,b} = \frac{1}{2} \mathcal{N}_{dca} (u_{c,b})_d = \frac{1}{2} \mathcal{N}_{dca} (u_{c,b} + u_{b,c})_d - \frac{1}{2} \mathcal{N}_{dca} u_{b,cd}$. The last term is zero since it is a contraction over the indices $(c, d)$ of an antisymmetric tensor $\mathcal{N}_{dca}$ and a symmetric tensor $u_{b,cd}$. This yields $\gamma_{a,b} = \mathcal{N}_{dca} \epsilon_{cb,d} = \mathcal{N}_{dca} \mathcal{M}_{cbe} \tilde{\epsilon}_{e,d}$ by (B.5), which we can rewrite in tensor form as

$$ \nabla \hat{\gamma}(\mathbf{X}) = \mathcal{G}' : \nabla \tilde{\epsilon}(\mathbf{X}), \quad (B.8) $$
where

$$ \mathcal{G}' = (\mathcal{N}^{T_{231}} \cdot \mathcal{M})^{T_{1423}}, \qquad (B.9) $$

as announced in (A.7)$_1$.

Differentiating one more time, we get $\nabla^2\boldsymbol{\gamma}(\mathbf{X}) = \mathcal{G}' : \nabla^2\tilde{\boldsymbol{\epsilon}}(\mathbf{X}) = (\mathcal{G}' \otimes \mathbf{I}_d)^{T_{124536}} : \nabla^2\tilde{\boldsymbol{\epsilon}}(\mathbf{X})$. We enforce the natural symmetries of $\nabla^2\boldsymbol{\gamma}$ with respect to the pairs of indices (2, 3) and (5, 6) by writing

$$ \nabla^2 \boldsymbol{\gamma}(\mathbf{X}) = \mathcal{G}'' : \nabla^2 \tilde{\boldsymbol{\epsilon}}(\mathbf{X}), \qquad (B.10) $$

as announced in (A.7)$_2$, where

$$ \mathcal{G}'' = ((\mathcal{G}' \otimes \mathbf{I}_d)^{T_{124536}})^{S_{23} \circ S_{56}}, \qquad (B.11) $$

and $S_{23} \circ S_{56}$ denotes successive symmetrization with the two pairs of indices.

Differentiating (A.7)$_3$ in components and using (B.8), we obtain the second gradient of displacement as $u_{a,bc} = (\mathcal{U}'_{abe} \tilde{\boldsymbol{\epsilon}}_e - \mathcal{N}_{abd} \gamma_d)_{,c} = \mathcal{U}'_{abe} \tilde{\boldsymbol{\epsilon}}_{e,c} - \mathcal{N}_{abd} \mathcal{G}'_{dcef} \tilde{\boldsymbol{\epsilon}}_{e,f} = (\mathcal{U}'_{abe} \delta_{cf} - \mathcal{N}_{abd} \mathcal{G}'_{dcef}) \tilde{\boldsymbol{\epsilon}}_{e,f}$, which we can rewrite as

$$ \nabla^2 \mathbf{u}(\mathbf{X}) = \mathcal{U}'' : \nabla^2 \tilde{\boldsymbol{\epsilon}}(\mathbf{X}) \qquad (B.12) $$

where $\mathcal{U}''$ is given by

$$ \mathcal{U}'' = (\mathcal{U}' \otimes \mathbf{I}_d)^{T_{12435}} - \mathcal{N}' \cdot \mathcal{G}'. \qquad (B.13) $$

Since $\nabla^2\mathbf{u}$ is symmetric with respect to its last pair of indices in the left-hand side of (B.12), $\mathcal{U}''$ is automatically symmetric with respect to the pair (2, 3). Equation (B.12) proves (A.7)$_4$.

Next, we have $\nabla^3\mathbf{u} = \mathcal{U}'' : \nabla^2\tilde{\boldsymbol{\epsilon}} = (\mathcal{U}'' \otimes \mathbf{I}_d)^{T_{1235647}} : \nabla^2\tilde{\boldsymbol{\epsilon}}$. The natural symmetries by any permutation of the indices (2, 3, 4) on the one hand and (6, 7) on the other hand, are enforced by writing

$$ \nabla^3 \mathbf{u}(\mathbf{X}) = \mathcal{U}''' : \nabla^2 \tilde{\boldsymbol{\epsilon}}(\mathbf{X}) \qquad (B.14) $$

as announced in (A.7)$_5$, where

$$ \mathcal{U}''' = ((\mathcal{U}'' \otimes \mathbf{I}_d)^{T_{1235647}})^{S_{234} \circ S_{67}}, \qquad (B.15) $$

and $S_{234} \circ S_{67}$ denotes symmetrization with respect to the set of indices in subscript.

The proof of the identities (A.7) is complete.
where

$$ \mathcal{G}' = (\mathcal{N}^{T_{231}} \cdot \mathcal{M})^{T_{1423}}, \qquad (B.9) $$

as announced in (A.7)$_1$.

Differentiating one more time, we get $\nabla^2\gamma(\mathbf{X}) = \mathcal{G}' : \nabla^2\tilde{\epsilon}(\mathbf{X}) = (\mathcal{G}' \otimes \mathbf{I}_d)^{T_{124536}} : \nabla^2\tilde{\epsilon}(\mathbf{X})$. We enforce the natural symmetries of $\nabla^2\gamma$ with respect to the pairs of indices (2, 3) and (5, 6) by writing

$$ \nabla^2\gamma(\mathbf{X}) = \mathcal{G}'' : \nabla^2\tilde{\epsilon}(\mathbf{X}), \qquad (B.10) $$

as announced in (A.7)$_2$, where

$$ \mathcal{G}'' = ((\mathcal{G}' \otimes \mathbf{I}_d)^{T_{124536}})^{S_{23} \circ S_{56}}, \qquad (B.11) $$

and $S_{23} \circ S_{56}$ denotes successive symmetrization with the two pairs of indices.

Differentiating (A.7)$_3$ in components and using (B.8), we obtain the second gradient of displacement as $u_{a,bc} = (\mathcal{U}'_{abe} \tilde{\epsilon}_e - \mathcal{N}_{abd} \gamma_d)_{,c} = \mathcal{U}'_{abe} \tilde{\epsilon}_{e,c} - \mathcal{N}_{abd} \mathcal{G}'_{dcef} \tilde{\epsilon}_{e,f} = (\mathcal{U}'_{abe} \delta_{cf} - \mathcal{N}_{abd} \mathcal{G}'_{dcef}) \tilde{\epsilon}_{e,f}$, which we can rewrite as

$$ \nabla^2 \mathbf{u}(\mathbf{X}) = \mathcal{U}'' : \nabla^2 \tilde{\epsilon}(\mathbf{X}) \qquad (B.12) $$

where $\mathcal{U}''$ is given by

$$ \mathcal{U}'' = (\mathcal{U}' \otimes \mathbf{I}_d)^{T_{12435}} - \mathcal{N}' \cdot \mathcal{G}'. \qquad (B.13) $$

Since $\nabla^2\mathbf{u}$ is symmetric with respect to its last pair of indices in the left-hand side of (B.12), $\mathcal{U}''$ is automatically symmetric with respect to the pair (2, 3). Equation (B.12) proves (A.7)$_4$.

Next, we have $\nabla^3\mathbf{u} = \mathcal{U}'' : \nabla^2\tilde{\epsilon} = (\mathcal{U}'' \otimes \mathbf{I}_d)^{T_{1235647}} : \nabla^2\tilde{\epsilon}$. The natural symmetries by any permutation of the indices (2, 3, 4) on the one hand and (6, 7) on the other hand, are enforced by writing

$$ \nabla^3 \mathbf{u}(\mathbf{X}) = \mathcal{U}''' : \nabla^2 \tilde{\epsilon}(\mathbf{X}) \qquad (B.14) $$

as announced in (A.7)$_5$, where

$$ \mathcal{U}''' = ((\mathcal{U}'' \otimes \mathbf{I}_d)^{T_{1235647}})^{S_{234} \circ S_{67}}, \qquad (B.15) $$

and $S_{234} \circ S_{67}$ denotes symmetrization with respect to the set of indices in subscript.

The proof of the identities (A.7) is complete.