# Supplementary Information:
Non-affinity of liquid networks and bicontinuous mesophases

Michael S. Dimitriyev<EMAIL_ADDRESS>Department of Polymer Science and
Engineering, University of Massachusetts, Amherst, MA 01003, USA  Department
of Materials Science and Engineering, Texas A&M University, College Station,
TX 77843  Xueyan Feng  Department of Macromolecular Science, State Key
Laboratory of Molecular Engineering of Polymers, Fudan University, Shanghai,
China 200437  Edwin L. Thomas  Department of Materials Science and
Engineering, Texas A&M University, College Station, TX 77843  Gregory M.
Grason<EMAIL_ADDRESS>Department of Polymer Science and Engineering,
University of Massachusetts, Amherst, MA 01003, USA

##  S1 Non-affinity of pre-stressed, equilibrium networks

We start with nodes in a reference configuration specified by positions
$\\{\mathbf{r}^{(m)}\\}$, assumed to be in an equilibrium configuration so
that $\delta H[\\{\mathbf{r}^{(m)}\\}]=0$. Moreover, we assume that every
strut in this configuration has length
$\ell_{0}=|\mathbf{r}^{(n)}-\mathbf{r}^{(m)}|$ for all node neighbors $m$ and
$n$. After an affine deformation, specified by deformation matrix $\Lambda$,
the node positions are given by $\\{\Lambda\mathbf{r}^{(m)}\\}$. The energy of
this configuration is given by
$H\left[\left\\{\Lambda\mathbf{r}^{(m)}\right\\}\right]$. Next, the network
will relax to a new equilibrium configuration. Without some constraints on
node positions, the network would simply relax to its non-deformed
configuration. These constraints come in two forms: (i) specification of
“boundary” nodes that do not relax [1]; (ii) specification of fixed lattice
vectors that constrain relate the positions of nodes in all copies of the unit
cells [2, 3]. In the following we will only consider nodes that relax (i.e.
non-boundary nodes) with the idea that fixed nodes impose boundary conditions
to the resulting equations. To obtain the equilibrium configuration, consider
the change in energy $\Delta H$ between an equilibrium configuration and the
affine-deformed configuration, given by

$\begin{split}\Delta H&\equiv
H\left[\left\\{\Lambda\mathbf{r}^{(m)}+\mathbf{u}^{(m)}\right\\}\right]-H\left[\left\\{\Lambda\mathbf{r}^{(m)}\right\\}\right]\\\
&=\sum_{mn}A_{mn}\bigg{\\{}h\left(\left|\Lambda\mathbf{r}^{(n)}+\mathbf{u}^{(n)}-\Lambda\mathbf{r}^{(m)}-\mathbf{u}^{(m)}\right|\right)-h\left(\left|\Lambda\mathbf{r}^{(n)}-\Lambda\mathbf{r}^{(m)}\right|\right)\bigg{\\}}\,,\end{split}$
(S1)

where relaxed node positions are displaced by vectors $\\{\mathbf{u}^{(m)}\\}$
from affine-deformed node positions. Next, we expand $\Delta H$ to quadratic
order in displacement:

$\begin{split}\Delta
H=\sum_{mn}A_{mn}\bigg{\\{}&\frac{1}{2}(\mathbf{u}^{(m)}-\mathbf{u}^{(n)})^{T}\mathbf{D}^{(mn)}\,(\mathbf{u}^{(m)}-\mathbf{u}^{(n)})\\\
&-\mathbf{f}^{(mn)}(\Lambda)\cdot(\mathbf{u}^{(m)}-\mathbf{u}^{(n)})+\mathcal{O}(u^{3})\bigg{\\}}\,,\end{split}$
(S2)

where $\mathbf{D}^{(mn)}$ is a “nodal” dynamical matrix that weighs the cost
of extending the strut joining node $m$ to $n$ in various directions and
$\mathbf{f}^{(mn)}(\Lambda)$ is a force acting along that strut, which depends
on the applied affine deformation $\Lambda$ of the unit cell. Since each node
of DP has pairs of anti-parallel struts, the net force on a given node is 0,
so DP has a completely affine response.

In the linear-response limit, we expand $\Delta H$ to leading order in the
symmetric strain

$\bm{\varepsilon}\equiv\frac{1}{2}\left(\bm{\Lambda}^{T}\bm{\Lambda}-\mathbbm{1}\right)$
(S3)

and find that the dynamical matrix is given by

$\mathbf{D}^{(mn)}\approx\frac{\tau(\ell_{0})}{\ell_{0}}\left(\mathbbm{1}-(1-\eta)\hat{\mathbf{R}}^{(mn)}\otimes\hat{\mathbf{R}}^{(mn)}\right)$
(S4)

and the strut force is given by

$\mathbf{f}^{(mn)}\approx\tau(\ell_{0})\,(1-\eta)\,\varepsilon^{(mn)}\,\hat{\mathbf{R}}^{(mn)}$
(S5)

where $\tau(\ell_{0})=h^{\prime}(\ell_{0})$ is the residual tension of the
struts in the undeformed configuration,
$\varepsilon^{(mn)}\equiv\hat{R}^{(mn)}_{i}\,\varepsilon_{ij}\,\hat{R}^{(mn)}_{j}$
is a projection of the strain tensor onto the strut joining nodes $m$ and $n$,
and $\eta\equiv\frac{\ell_{0}\,\tau^{\prime}(\ell_{0})}{\tau(\ell_{0})}$ is a
parameter that characterizes the tension-length constitutive relation.
Finally, the linear response $u\sim\varepsilon$ can be determined by
minimizing $\Delta H$ with respect to the non-affine displacements of each
node within the unit cell of a given network. This yields, for each node $m$,
an equilibrium condition,

$\sum_{n\in\mathcal{N}_{m}}\left[\mathbf{D}^{(mn)}(\mathbf{u}^{(m)}-\mathbf{u}^{(n)})-\mathbf{f}^{(mn)}\right]=0\,,$
(S6)

where $\mathcal{N}_{m}$ is the set of neighbors of node $m$. Given that the
nodal dynamical matrix $\mathbf{D}$ and the affine force $\mathbf{f}$ have the
residual stress $\tau(\ell_{0})$ as an overall scale factor, we find that (i)
the following results only hold for residually stressed networks (for which
$\tau(\ell_{0})\neq 0$) and (ii) the non-affine response is completely
characterized by the parameter $\eta$.

###  S1.1 Trihedral and tetrahedral graphs

First we solve Eq. S6 for the case of the trihedral and tetrahedral graphs
shown in Fig. 2 of the main text. In these cases, there is only one central
node that relaxes to a new equilibrium position and the remaining 3
(trihedral) and 4 (tetrahedral) nodes are “boundary nodes.” It is useful to
introduce coordinates to the nodes. The central node $\mathbf{r}^{(0)}$ will
be located at the origin for both graphs. For the trihedral graph, the
boundary nodes will be located at

$\begin{split}\mathbf{r}^{(1)}&=\frac{\ell_{0}}{\sqrt{2}}\left(1,\,-1,\,0\right)\\\
\mathbf{r}^{(2)}&=\frac{\ell_{0}}{\sqrt{2}}\left(0,\,1,\,-1\right)\\\
\mathbf{r}^{(3)}&=\frac{\ell_{0}}{\sqrt{2}}\left(-1,\,0,\,1\right)\end{split}$
(S7)

For the tetrahedral graph, the boundary nodes will be located at

$\begin{split}\mathbf{r}^{(1)}&=\frac{\ell_{0}}{\sqrt{3}}\left(1,\,1,\,1\right)\\\
\mathbf{r}^{(2)}&=\frac{\ell_{0}}{\sqrt{3}}\left(1,\,-1,\,-1\right)\\\
\mathbf{r}^{(3)}&=\frac{\ell_{0}}{\sqrt{3}}\left(-1,\,1,\,-1\right)\\\
\mathbf{r}^{(4)}&=\frac{\ell_{0}}{\sqrt{3}}\left(-1,\,-1,\,1\right)\end{split}$
(S8)

For both graphs, we can rewrite Eq. S6 for the central node as

$\mathbf{D}\,\mathbf{u}=\mathbf{f}\,,$ (S9)

where we suppress the superscript on the displacement vector
$\mathbf{u}\equiv\mathbf{u}^{(0)}$. Then the scale of the non-affine
displacement is given by

$\frac{|\mathbf{u}|^{2}}{\ell_{0}^{2}}=\frac{1}{\ell_{0}^{2}}\mathbf{f}^{T}(\mathbf{D}^{-1})^{T}\mathbf{D}^{-1}\mathbf{f}\,.$
(S10)

For the trihedral network, the dynamical matrix $\mathbf{D}$ is given by

$\begin{split}\mathbf{D}&=\frac{3\tau(\ell_{0})(1+\eta)}{2\ell_{0}}\left\\{\mathbbm{1}+\frac{1-\eta}{1+\eta}\hat{\mathbf{N}}\otimes\hat{\mathbf{N}}\right\\}\end{split}$
(S11)

where $\hat{\mathbf{N}}$ is the direction normal to the plane of the trihedral
graph (here, $\hat{\mathbf{N}}=(1,1,1)/\sqrt{3}$). For the tetrahedral
network, the dynamical matrix is given by

$\mathbf{D}=\frac{8\tau(\ell_{0})(1+\frac{\eta}{2})}{3\ell_{0}}\mathbbm{1}\,.$
(S12)

We will focus on the case of uniaxial deformation, where

$\bm{\varepsilon}=\varepsilon\,\hat{\mathbf{d}}\otimes\hat{\mathbf{d}}\,.$
(S13)

Consequently, for the trihedral graph,

$\frac{|\mathbf{u}|^{2}}{\ell_{0}^{2}}=\frac{1}{4}\left(\frac{1-\eta}{1+\eta}\right)^{2}\,|\hat{\mathbf{d}}\times\hat{\mathbf{N}}|^{4}\,\varepsilon^{2}=\Gamma_{\rm
tri.}(\hat{\mathbf{d}})\,\varepsilon^{2}$ (S14)

and for the tetrahedral graph,

$\frac{|\mathbf{u}|^{2}}{\ell_{0}^{2}}=\frac{1}{6}\left(\frac{1-\eta}{1+\frac{\eta}{2}}\right)^{2}\,\left[1-\hat{d}_{x}^{4}-\hat{d}_{y}^{4}-\hat{d}_{z}^{4}\right]\,\varepsilon^{2}=\Gamma_{\rm
tet.}(\hat{\mathbf{d}})\,\varepsilon^{2}\,.$ (S15)

Writing $\Gamma(\hat{\mathbf{d}})=\Gamma^{\rm max}\psi(\hat{\mathbf{d}})$, we
have

$\begin{split}\Gamma^{\rm tri.}_{\rm
max}&=\frac{1}{4}\left(\frac{1-\eta}{1+\eta}\right)^{2}\\\ \psi^{\rm
tri.}(\hat{\mathbf{d}})&=|\hat{\mathbf{d}}\times\hat{\mathbf{N}}|^{4}\end{split}$
(S16)

and

$\begin{split}\Gamma^{\rm tet.}_{\rm
max}&=\frac{1}{9}\left(\frac{1-\eta}{1+\frac{\eta}{2}}\right)^{2}\\\ \psi^{\rm
tet.}(\hat{\mathbf{d}})&=\frac{3}{2}\left[1-\hat{d}_{x}^{4}-\hat{d}_{y}^{4}-\hat{d}_{z}^{4}\right]\end{split}$
(S17)

Finally, we determine the average $\overline{\Gamma}$ by integrating
$\Gamma(\hat{\mathbf{d}})$ over the 2-sphere,

$\overline{\Gamma}=\Gamma_{\rm max}\oint_{\mathbbm{S}^{2}}\frac{{\rm
d}\Omega}{4\pi}\,\psi(\hat{\mathbf{d}})\,,$ (S18)

yielding

$\begin{split}\overline{\Gamma}^{\rm tri.}&=\frac{8}{15}\Gamma^{\rm tri.}_{\rm
max}\\\ \overline{\Gamma}^{\rm tet.}&=\frac{2}{5}\Gamma^{\rm tet.}_{\rm
max}\end{split}$ (S19)

###  S1.2 DG and DD networks

Triply-periodic networks do not possess boundary nodes; instead, network
equilibrium is determined by the position of nodes within a given unit cell.
In this case, Eq. S6 must be solved with periodic boundary conditions. A
single DG network has initial node positions given by

$\begin{split}\mathbf{r}^{(1)}&=\frac{1}{8}\left(1,1,1\right)\\\
\mathbf{r}^{(2)}&=\frac{1}{8}\left(3,7,5\right)\\\
\mathbf{r}^{(3)}&=\frac{1}{8}\left(7,5,3\right)\\\
\mathbf{r}^{(4)}&=\frac{1}{8}\left(5,3,7\right)\\\
\mathbf{r}^{(5)}&=\frac{1}{8}\left(5,5,5\right)\\\
\mathbf{r}^{(6)}&=\frac{1}{8}\left(1,7,3\right)\\\
\mathbf{r}^{(7)}&=\frac{1}{8}\left(7,3,1\right)\\\
\mathbf{r}^{(8)}&=\frac{1}{8}\left(3,1,7\right)\end{split}$ (S20)

and edge relations given by

$\begin{split}\langle
1,6\rangle&:\mathbf{r}^{(1)}\to\mathbf{r}^{(6)}+(0,-1,0)\\\ \langle
1,7\rangle&:\mathbf{r}^{(1)}\to\mathbf{r}^{(7)}+(-1,0,0)\\\ \langle
1,8\rangle&:\mathbf{r}^{(1)}\to\mathbf{r}^{(6)}+(0,0,-1)\\\ \langle
2,5\rangle&:\mathbf{r}^{(2)}\to\mathbf{r}^{(5)}+(0,0,0)\\\ \langle
2,6\rangle&:\mathbf{r}^{(2)}\to\mathbf{r}^{(6)}+(0,0,0)\\\ \langle
2,8\rangle&:\mathbf{r}^{(2)}\to\mathbf{r}^{(8)}+(0,1,0)\\\ \langle
3,5\rangle&:\mathbf{r}^{(3)}\to\mathbf{r}^{(5)}+(0,0,0)\\\ \langle
3,6\rangle&:\mathbf{r}^{(3)}\to\mathbf{r}^{(6)}+(1,0,0)\\\ \langle
3,7\rangle&:\mathbf{r}^{(3)}\to\mathbf{r}^{(7)}+(0,0,0)\\\ \langle
4,5\rangle&:\mathbf{r}^{(4)}\to\mathbf{r}^{(5)}+(0,0,0)\\\ \langle
4,7\rangle&:\mathbf{r}^{(4)}\to\mathbf{r}^{(7)}+(0,0,1)\\\ \langle
4,8\rangle&:\mathbf{r}^{(4)}\to\mathbf{r}^{(8)}+(0,0,0)\end{split}$ (S21)

where $\langle mn\rangle$ denotes a pair of joined nodes, with strut joining
node $\mathbf{r}^{(m)}$ to $\mathbf{r}^{(n)}+\Delta\mathbf{e}^{(mn)}$, where
$\Delta\mathbf{e}^{(mn)}$ is an additional translation vector that accounts
for struts that pass through a unit cell boundary. Similarly, a single DD
network has initial node positions given by

$\begin{split}\mathbf{r}^{(1)}&=\frac{1}{8}\left(1,1,1\right)\\\
\mathbf{r}^{(2)}&=\frac{1}{8}\left(3,3,3\right)\\\
\mathbf{r}^{(3)}&=\frac{1}{8}\left(1,5,5\right)\\\
\mathbf{r}^{(4)}&=\frac{1}{8}\left(5,1,5\right)\\\
\mathbf{r}^{(5)}&=\frac{1}{8}\left(5,5,1\right)\\\
\mathbf{r}^{(6)}&=\frac{1}{8}\left(3,7,7\right)\\\
\mathbf{r}^{(7)}&=\frac{1}{8}\left(7,3,7\right)\\\
\mathbf{r}^{(8)}&=\frac{1}{8}\left(7,7,3\right)\end{split}$ (S22)

and edge relations given by

$\begin{split}\langle
1,2\rangle&:\mathbf{r}^{(1)}\to\mathbf{r}^{(2)}+(0,0,0)\\\ \langle
1,6\rangle&:\mathbf{r}^{(1)}\to\mathbf{r}^{(6)}+(0,-1,-1)\\\ \langle
1,7\rangle&:\mathbf{r}^{(1)}\to\mathbf{r}^{(7)}+(-1,0,-1)\\\ \langle
1,8\rangle&:\mathbf{r}^{(1)}\to\mathbf{r}^{(8)}+(-1,-1,0)\\\ \langle
2,3\rangle&:\mathbf{r}^{(2)}\to\mathbf{r}^{(3)}+(0,0,0)\\\ \langle
2,4\rangle&:\mathbf{r}^{(2)}\to\mathbf{r}^{(4)}+(0,0,0)\\\ \langle
2,5\rangle&:\mathbf{r}^{(2)}\to\mathbf{r}^{(5)}+(0,0,0)\\\ \langle
3,6\rangle&:\mathbf{r}^{(3)}\to\mathbf{r}^{(6)}+(0,0,0)\\\ \langle
3,7\rangle&:\mathbf{r}^{(3)}\to\mathbf{r}^{(7)}+(-1,0,0)\\\ \langle
3,8\rangle&:\mathbf{r}^{(3)}\to\mathbf{r}^{(8)}+(-1,0,0)\\\ \langle
4,6\rangle&:\mathbf{r}^{(4)}\to\mathbf{r}^{(6)}+(0,-1,0)\\\ \langle
4,7\rangle&:\mathbf{r}^{(4)}\to\mathbf{r}^{(7)}+(0,0,0)\\\ \langle
4,8\rangle&:\mathbf{r}^{(4)}\to\mathbf{r}^{(8)}+(0,-1,0)\\\ \langle
5,6\rangle&:\mathbf{r}^{(5)}\to\mathbf{r}^{(6)}+(0,0,-1)\\\ \langle
5,7\rangle&:\mathbf{r}^{(5)}\to\mathbf{r}^{(7)}+(0,0,-1)\\\ \langle
5,8\rangle&:\mathbf{r}^{(5)}\to\mathbf{r}^{(8)}+(0,0,0)\\\ \end{split}$ (S23)

Note that the DD skeletal graph detailed above is generated from the symmetry
group $Fd\overline{3}m$ and contains 8 copies of the standard DD unit cell
generated from the $Pn\overline{3}m$ space group. The above form allows for a
straightforward representation of single networks, similar to the
representation used for a the DG skeletal graph.

Next, it is convenient to represent $\\{\mathbf{u}^{(m)}\\}$ as a $3N_{\rm
nodes}$-dimensional vector

$\mathbf{U}=\left(\begin{array}[]{c}u^{(1)}_{x}\\\ u^{(1)}_{y}\\\
u^{(1)}_{z}\\\ u^{(2)}_{x}\\\ \vdots\\\ u^{(N_{\rm
nodes})}_{z}\end{array}\right),$ (S24)

where $N_{\rm nodes}$ is the number of nodes per unit cell (8 for each of the
single networks above). Each $3\times 3$ nodal dynamical matrix
$\mathbf{D}^{(mn)}$ can then be folded into a $3N_{\rm nodes}\times 3N_{\rm
nodes}$ dynamical matrix $\mathbbm{D}$ and the sum over forces at each node,
$\sum_{n\in\langle mn\rangle}\mathbf{f}^{(mn)}$, can be written as a $3N_{\rm
nodes}$-dimensional vector $\mathbf{F}$. As a result, the equilibrium node
displacement vector $\mathbf{U}$ satisfies a modified form of Eq. S6, namely

$\mathbbm{D}\mathbf{U}=\mathbf{F}\,.$ (S25)

Since the energy of the network is invariant under continuous, rigid
translations (i.e., $\mathbf{u}^{(m)}={\rm const.}$ for all $m$), the
dynamical matrix $\mathbbm{D}$ possesses three zero-eigenvalues and is
therefore not invertible. To solve Eq. S25, we restrict the possible solution
set of $\mathbf{U}$ to exclude rigid translations. We do this via the _Moore-
Penrose pseudoinverse_ $\mathbbm{D}^{+}$, which acts as the inverse of
$\mathbbm{D}$ in the complement to its nullspace. Thus, solutions are given by

$\mathbf{U}=\mathbbm{D}^{+}\mathbf{F}$ (S26)

and the non-affine response is given by

$\frac{\langle|\mathbf{u}|^{2}\rangle}{\ell_{0}^{2}}=\frac{\left|\mathbf{U}\right|^{2}}{N_{\rm
nodes}\ell_{0}^{2}}=\frac{1}{N_{\rm
nodes}{\ell_{0}^{2}}}\mathbf{U}^{T}\left(\mathbbm{D}^{+}\right)^{T}\mathbbm{D}^{+}\mathbf{U}\,.$
(S27)

We use Mathematica to compute the pseudoinverse $\mathbbm{D}^{+}$ and the non-
affine response to strain components $\varepsilon_{ij}$ via Eq. S27, yielding

$\begin{split}\frac{\langle|\mathbf{u}|^{2}\rangle_{\rm
DG}}{\ell_{0}^{2}}&=\frac{3}{12}\left(\frac{1-\eta}{1+3\eta}\right)^{2}\left[\varepsilon_{xx}^{2}+\varepsilon_{yy}^{2}+\varepsilon_{zz}^{2}-\varepsilon_{xx}\varepsilon_{yy}-\varepsilon_{yy}\varepsilon_{zz}-\varepsilon_{zz}\varepsilon_{xx}\right]+\frac{1}{4}\left(\frac{1-\eta}{1+\eta}\right)^{2}\left[\varepsilon_{xy}^{2}+\varepsilon_{yz}^{2}+\varepsilon_{zx}^{2}\right]\\\
\frac{\langle|\mathbf{u}|^{2}\rangle_{\rm
DD}}{\ell_{0}^{2}}&=\frac{1}{3}\left(\frac{1-\eta}{2+\eta}\right)^{2}\left[\varepsilon_{xy}^{2}+\varepsilon_{yz}^{2}+\varepsilon_{zx}^{2}\right]\end{split}$
(S28)

from which we determine the non-affinity components.

##  S2 Comparisons with strain-stiffening polymers

Figure S1:  (a) Effective force law parameter $\eta$ for the Freely-Jointed
Chain (FJC) model (solid) and the Worm-Like Chain (WLC) model (dashed) as a
function of chain stretch relative to chain length $\lambda$. (b) Non-affinity
estimates from linear response theory for the trihedral graph (purple) and the
tetrahedral graph (orange) for both chain models.

Polymer networks strain-stiffen due to the finite extensibility of polymer
chains when strongly stretched. Two models that describe polymers of finite
extent are the Freely-Jointed Chain (FJC) model of flexible polymers and the
Worm-Like Chain (WLC) model [4, 5]. The force $f$ needed to stretch a chain to
end-separation $r$ is well-approximated within both models as

$\frac{fb}{k_{\rm
B}T}=\frac{2r}{Nb}+\frac{1}{p}\left(1-\frac{r}{Nb}\right)^{-p}-\frac{1}{p}\,,$
(S29)

where $b$ is the statistical segment length, $N$ is the polymerization index,
$k_{\rm B}T$ is the thermal energy scale, and $p=1$ (FJM) or 2 (WLC). Using
these force laws as constitutive models for strut tension, we can estimate the
non-affine response of finitely-extensible polymer networks from the
characteristic parameter $\eta(r)=rf^{\prime}/f$. Defining $\lambda\equiv
r/(Nb)$ as the stretch fraction relative to total chain length, the
characteristic parameter is given by

$\eta(\lambda)=\frac{1+\frac{1}{2}\left(1-\lambda\right)^{-p-1}}{1+\frac{1}{2p\lambda}\left[\left(1-\lambda\right)^{-p}-1\right]}\,,$
(S30)

which only depends on $\lambda$ and $p$. As shown in Fig. S1(a), the parameter
$\eta\geq 1$, correctly categorizing these models as strain-stiffening; as
$\lambda\to 1$, $\eta\to\infty$. Consequently, as $\lambda\to 1$, the expected
non-affine response approaches values similar to those expected by LNT, as
shown in Fig. S1(b).

##  S3 Symmetries of linear response

The non-affinity tensor $\bm{\Gamma}$ provides the quadratic-order strain
coupling in the linear-response regime of deformation. Much like the elastic
modulus tensor of linear elasticity, $\bm{\Gamma}$ possesses a variety of
symmetries, including

  * •

Minor symmetries: $\Gamma^{ijkl}=\Gamma^{jikl}=\Gamma^{ijlk}$

  * •

Major symmetry: $\Gamma^{ijkl}=\Gamma^{klij}$

  * •

Symmetries related to the crystal system

Due to the symmetry of the strain tensor $\varepsilon_{ij}=\varepsilon_{ji}$
and the minor symmetries of $\bm{\Gamma}$, it is useful to represent both
using Voigt notation $\varepsilon_{ij}\mapsto\varepsilon_{I}$ and
$\Gamma^{ijkl}\mapsto\Gamma^{IJ}$, where $I,J,\dots\in\\{1,2,\dots,6\\}$
represents, in order, the set of indices $\\{xx,yy,zz,yz,xz,xy\\}$. Setting
$\\{\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}\\}=\\{\varepsilon_{xx},\varepsilon_{yy},\varepsilon_{zz}\\}$
and
$\\{\varepsilon_{4},\varepsilon_{5},\varepsilon_{6}\\}=\\{2\varepsilon_{yz},2\varepsilon_{xz},2\varepsilon_{xy}\\}$
allows for the identity
$\langle|\mathbf{u}|^{2}\rangle/\ell_{0}^{2}=\Gamma^{ijkl}\varepsilon_{ij}\varepsilon_{kl}=\Gamma^{IJ}\varepsilon_{I}\varepsilon_{J}$.
The cubic crystal system is invariant under $\pi/2$ rotations of each axis;
much like the elastic modulus tensor, $\Gamma$ thus consists of three non-
zero, independent constants (reduced from 21 for a full triclinic system),

$\Gamma^{IJ}=\left\\{\begin{array}[]{ccc|ccc}A&B&B&0&0&0\\\ B&A&B&0&0&0\\\
B&B&A&0&0&0\\\ \hline\cr 0&0&0&C&0&0\\\ 0&0&0&0&C&0\\\
0&0&0&0&0&C\end{array}\right\\}$ (S31)

where, in order for $\langle|\mathbf{u}|^{2}\rangle$ to be positive
semidefinite ($\langle|\mathbf{u}|^{2}\rangle\geq 0$), the constants $A$, $B$,
and $C$ must be such that ${\rm det}\,\bm{\Gamma}=(A-B)^{2}(A+2B)C^{3}\geq 0$.
Thus, $\langle|\mathbf{u}|^{2}\rangle$ is given by

$\frac{\langle|\mathbf{u}|^{2}\rangle}{\ell_{0}^{2}}=A\left(\varepsilon^{2}_{xx}+\varepsilon^{2}_{yy}+\varepsilon^{2}_{zz}\right)+2B\left(\varepsilon_{xx}\varepsilon_{yy}+\varepsilon_{yy}\varepsilon_{zz}+\varepsilon_{zz}\varepsilon_{xx}\right)+4C\left(\varepsilon_{xy}^{2}+\varepsilon_{yz}^{2}+\varepsilon_{zx}^{2}\right)\,.$
(S32)

Defining the deviatoric strain tensor (the traceless part of the strain
tensor) as

$\hat{\varepsilon}_{ij}\equiv\varepsilon_{ij}-\frac{1}{3}\delta_{ij}\varepsilon_{kk}\,,$
(S33)

we have

$\frac{\langle|\mathbf{u}|^{2}\rangle}{\ell_{0}^{2}}=A\left(\hat{\varepsilon}^{2}_{xx}+\hat{\varepsilon}^{2}_{yy}+\hat{\varepsilon}^{2}_{zz}\right)+2B\left(\hat{\varepsilon}_{xx}\hat{\varepsilon}_{yy}+\hat{\varepsilon}_{yy}\hat{\varepsilon}_{zz}+\hat{\varepsilon}_{zz}\hat{\varepsilon}_{xx}\right)+\frac{A+2B}{3}\left({\rm
tr}\,\bm{\varepsilon}\right)^{2}+4C\left(\varepsilon_{xy}^{2}+\varepsilon_{yz}^{2}+\varepsilon_{zx}^{2}\right)\,.$
(S34)

For the class of mechanical networks considered here, a force-balance is
unchanged under a uniform re-scaling of space; consequently, $A+2B=0$. The
fact that this relationship implies that ${\rm det}\,\bm{\Gamma}=0$ merely
reflects the dilational symmetry of such mechanical networks. Thus, we can
write

$\begin{split}\frac{\langle|\mathbf{u}|^{2}\rangle}{\ell_{0}^{2}}&=\frac{A}{2}\left((\hat{\varepsilon}_{xx}-\hat{\varepsilon}_{yy})^{2}+(\hat{\varepsilon}_{yy}-\hat{\varepsilon}_{zz})^{2}+(\hat{\varepsilon}_{zz}-\hat{\varepsilon}_{xx})^{2}\right)+4C\left(\varepsilon_{xy}^{2}+\varepsilon_{yz}^{2}+\varepsilon_{zx}^{2}\right)\\\
&=3A\sum_{i}\frac{\hat{\varepsilon}_{ii}^{2}}{2}+2C\sum_{i\neq
j}\varepsilon_{ij}^{2}\end{split}$ (S35)

and identify the first term as relating to volume-preserving extensional
strain $\varepsilon_{\rm ext}$, where

$\varepsilon_{\rm
ext}\equiv\sqrt{\sum_{i}\frac{\hat{\varepsilon}_{ii}^{2}}{2}}=\sqrt{\frac{\hat{\varepsilon}_{xx}^{2}+\hat{\varepsilon}_{yy}^{2}+\hat{\varepsilon}_{zz}^{2}}{2}}\Bigg{(}=\sqrt{\frac{(\hat{\varepsilon}_{xx}-\hat{\varepsilon}_{yy})^{2}+(\hat{\varepsilon}_{yy}-\hat{\varepsilon}_{zz})^{2}+(\hat{\varepsilon}_{zz}-\hat{\varepsilon}_{xx})^{2}}{6}}\Bigg{)}$
(S36)

(here, equality with the term in the parenthesis is ensured since ${\rm
tr}\,\hat{\bm{\varepsilon}}=0$) and the second term as relating to the shear
strain $\varepsilon_{\rm shear}$, where

$\varepsilon_{\rm shear}\equiv\sqrt{\sum_{i\neq
j}\varepsilon_{ij}^{2}}=\sqrt{2\left(\varepsilon_{xy}^{2}+\varepsilon_{yz}^{2}+\varepsilon_{zx}^{2}\right)}\,.$
(S37)

Thus, we have $3A=\Gamma_{\rm ext}$ and $2C=\Gamma_{\rm shear}$.

Finally, to plot multiple strain paths on the same axis, we use the above
result that to linear order, the response can be represented as a function of
$\varepsilon_{\rm ext}$ and $\varepsilon_{\rm shear}$, to project the full
deformation path on a linearized subspace spanned by $\\{\varepsilon_{\rm
ext},\varepsilon_{\rm shear}\\}$. We can then introduce polar angle $\alpha$
via

$\tan\alpha\equiv\frac{\varepsilon_{\rm shear}}{\varepsilon_{\rm ext}}$ (S38)

and radial measure $\varepsilon_{r}$ via

$\varepsilon_{r}^{2}\equiv\varepsilon_{\rm ext}^{2}+\varepsilon_{\rm
shear}^{2}=\sum_{i}\frac{\hat{\varepsilon}_{ii}^{2}}{2}+\sum_{i\neq
j}\varepsilon_{ij}^{2}=\frac{1}{2}{\rm
tr}\left(\hat{\bm{\varepsilon}}^{2}\right)\,,$ (S39)

where the final term in the above equality is used to plot non-affinity,
length variation, and volume variation data on a common axis. Note that the
above arguments showing that the linear non-affine response
$\langle|\mathbf{u}|^{2}\rangle$ depends only on volume-preserving extensional
and shear strain measurements holds for the linear response of strut length
variations $\langle\Delta\ell^{2}\rangle$ and nodal volume variations
$\langle\Delta V^{2}\rangle$ as both measures are similarly expected to be
invariant under uniform dilations.

##  S4 SCF procedure

Self-consistent field (SCF) theory calculations were performed using the open-
source PSCF software [6]. Using parameters relevant to the experimental data,
we set A-block fractions as $f_{\rm A}=0.4$ for DG and $f_{\rm B}=0.42$ for DD
and assumed equal statistical segment lengths for both blocks. Finally, we
chose $\chi N=30$ as representative of the intermediate-segregation regime.

In order to calculate equilibrium density fields along simulated deformation
paths, we first calculate the equilibrium fields assuming a cubic unit cell
and constraining Fourier modes to obey the relevant space group; unit cell
volume was allowed to fluctuate. Using the equilibrium unit cell dimensions
and fields obtained from this method as input, we then relax the space group
constraint, impose a fixed volume (determined from the result of the first
calculation), and calculate the resulting equilibrium structure. We then
iteratively adjust the unit cell dimension (at fixed volume), and use the
equilibrium fields produced by the prior calculation as a starting point for
each new calculation. In this way, we are able to sweep out a deformation path
at fixed volume.

In general, the deformed, triclinic unit cell can be parameterized by 3
lattice constants $a$, $b$, and $c$, as well as 3 angles $\alpha$, $\beta$,
and $\gamma$. Adopting the parameterization used in PSCF, the each triclinic
unit cell can be transformed from a cubic unit cell by the deformation matrix

$\bm{\Lambda}=\frac{1}{a_{0}}\left(\begin{array}[]{ccc}a&0&0\\\
b\cos\gamma&b\sin\gamma&0\\\
c\cos\alpha\sin\beta&c\sin\alpha\sin\beta&c\cos\beta\end{array}\right)\,,$
(S40)

where $a_{0}$ is the lattice constant of the (undeformed) cubic unit cell. We
parameterize extensional deformations are by the 1D family of deformation
matrices $\bm{\Lambda}_{\rm ext}(s)$, where

$\bm{\Lambda}_{\rm ext}(s)=\left(\begin{array}[]{ccc}1+s&0&0\\\
0&\frac{1}{\sqrt{1+s}}&0\\\ 0&0&\frac{1}{\sqrt{1+s}}\end{array}\right)$ (S41)

and shear deformations by the 1D family $\bm{\Lambda}_{\rm shear}(s)$, where

$\bm{\Lambda}_{\rm shear}(s)=\left(\begin{array}[]{ccc}1&0&0\\\ s&1&0\\\
0&0&1\end{array}\right)\,.$ (S42)

To linear order, $\varepsilon_{\rm ext}=\frac{\sqrt{3}}{2}s$ and
$\varepsilon_{\rm shear}=0$ for the extensional deformation and
$\varepsilon_{\rm shear}=\frac{s}{\sqrt{2}}$ and $\varepsilon_{\rm
ext}=\mathcal{O}(s^{2})$ for the shear deformation.

Figure S2:  Deformation paths #1, #2, and #3 for DG (solid curves), projected
onto the space spanned by the linear combination of extensional strain
$\varepsilon_{\rm ext}$ and shear strain $\varepsilon_{\rm shear}$. These
curves pass through points representing the three triclinic DG grains from
experimental data. Dashed lines represent the linearization of these paths in
the small-strain regime. Also shown is the polar representation of this space,
where $\alpha$ is the polar angle and $\sqrt{({\rm
tr}\,\hat{\bm{\varepsilon}}^{2})/2}=\sqrt{\varepsilon_{\rm
ext}^{2}+\varepsilon_{\rm shear}^{2}}$ is the radial coordinate. The three
linearized paths have polar angles $\alpha_{1}\approx 66^{\circ}$,
$\alpha_{2}\approx 49^{\circ}$, and $\alpha_{3}\approx 58^{\circ}$. Contours
represent level sets of constant non-affine response
$\langle|\mathbf{u}|^{2}\rangle/\ell_{0}^{2}$, in the linear response
approximation, as obtained from SCF measurements.

In order to simulate a more general, volume-preserving, deformation path, we
take advantage of the exponential map

$\bm{\Lambda}(t)=e^{t\,\bm{\lambda}}$ (S43)

for some generating matrix $\lambda$. Due to the identity

${\rm det}\,\bm{\Lambda}(t)={\rm det}\,e^{t\,\bm{\lambda}}=e^{t\,{\rm
tr}\,\bm{\lambda}}\,,$ (S44)

it is guaranteed that if the generating matrix is traceless (i.e. ${\rm
tr}\,\bm{\lambda}=0$) then volume preserved at each point along the
deformation (i.e. ${\rm det}\,\bm{\Lambda}(t)=1$). In order to determine an
appropriate generating matrix $\bm{\lambda}$, we select a “target” deformation
matrix $\bm{\Lambda}_{1}=\bm{\Lambda}(t=1)$. Then the generating matrix is
found by the logarithm map, $\bm{\lambda}={\rm ln}\,\bm{\Lambda}_{1}$. We
apply this method to simulate deformation “trajectories” that originate at the
cubic configuration (where $\bm{\Lambda}(0)=\mathbbm{1}$) and pass through the
experimental point at $t=1$ (where $\bm{\Lambda}(1)=\bm{\Lambda}_{1}$). As
shown in Fig. S2, the resulting deformation paths are slightly curved, when
projected onto the space spanned by the linear combination of extensional
strain $\varepsilon_{\rm ext}$ and shear strain $\varepsilon_{\rm shear}$.
This curvature arises from non-linear strain coupling that occurs for large
deformations.

##  S5 Measuring non-affine response from SCF calculations

SCF calculations output density fields for the A- and B-block monomers. Using
these density fields and the skeletonization process outlined in ref. [7], we
extract skeletal graphs for each relaxed network structure. We then compare
these skeletal graphs to affine deformations of the skeletal graphs of the
cubic structures. Each node $m$ in the skeletal relaxed graph can be obtained
from a corresponding node in affinely-deformed skeletal graph by a
displacement vector $\mathbf{u}^{(m)}=\mathbf{r}^{(m)}-\mathbf{r}_{\rm
aff}^{(m)}$. Given that the SCF calculations of the deformed networks involve
the full set of Fourier modes within a single triclinic unit cell, periodic
boundary conditions allow for free translation of the density field within the
unit cell; we find that indeed there is a slight translation of the density
field with respect to the usual choice of origin. We therefore add a
translation vector $\mathbf{t}$ so that the non-affine displacement vectors
are given by

$\mathbf{u}^{(m)}=\mathbf{r}^{(m)}-(\mathbf{r}_{\rm aff}^{(m)}+\mathbf{t})\,.$
(S45)

Then we determine a measure of the total displacement by the average displaced
distance,

$\langle|\mathbf{u}|^{2}\rangle=\sum_{m=1}^{N_{\rm
nodes}}\frac{\left|\mathbf{u}^{(m)}\right|^{2}}{N_{\rm nodes}}\,,$ (S46)

where $N_{\rm nodes}$ is the number of nodes within the unit cell,
corresponding to the multiplicity of the Wyckoff position: $N_{\rm nodes}=16$
for DG, where nodes sit at the 16b Wyckoff positions; $N_{\rm nodes}=2$ for
DD, where nodes sit at the 2a Wyckoff positions. Finally, we minimize
$\langle|\mathbf{u}|^{2}\rangle$ with respect to the translation vector
$\mathbf{t}$ to “subtract out” the residual affine deformation that comes from
translation, leaving a measure of non-affine response.

###  S5.1 Determination of non-affinity components

In order to measure $\Gamma_{\rm ext}$ and $\Gamma_{\rm shear}$ from this
process, we perform a pair of sweeps along deformation paths
$\bm{\Lambda}_{\rm ext}(s)$ and $\bm{\Lambda}_{\rm shear}(s)$ involving (to
linear order) only extensional strain and shear strain components,
respectively. Using the method outlined above, we extract a pair of non-affine
responses, $\langle|\mathbf{u}|^{2}\rangle_{\rm ext}/\ell_{0}^{2}$ and
$\langle|\mathbf{u}|^{2}\rangle_{\rm shear}/\ell_{0}^{2}$ as functions of
individual strain components $\varepsilon_{\rm ext}$ and $\varepsilon_{\rm
shear}$, respectively. Next, we fit these data to functions of the form
$\langle|\mathbf{u}|^{2}\rangle/\ell_{0}^{2}=\Gamma\,\varepsilon^{2}+C$, where
$C$ is an offset that represents the additional non-affine response
originating from numerical error at $\varepsilon=0$. However, given that we
only seek to quantify the _linear response_ (i.e. the lowest-order, quadratic
contribution to $\langle|\mathbf{u}|^{2}\rangle$), the fit must be taken over
a suitably chosen region of data. This is achieved by choice of a cut-off
strain $\varepsilon_{\rm cut}$; each non-affinity component $\Gamma$
determined through the quadratic fit is thus a function of $\varepsilon_{\rm
cut}$.

Given that the non-affinity measurement depends on the choice of
$\varepsilon_{\rm cut}$, our accepted value involves a degree of error from
this choice. To estimate this error, we consider the variation of
$\Gamma(\varepsilon_{\rm cut})$ over a range of $\varepsilon_{\rm cut}$
values. This range is given by the interval $\left[\varepsilon^{<}_{\rm
cut},\varepsilon^{>}_{\rm cut}\right]$; the lower limit, $\varepsilon^{<}_{\rm
cut}$, is chosen to include enough data to provide a reasonable parabolic fit,
given the noise present in SCF calculations; we choose $\varepsilon^{<}_{\rm
cut}=0.02$. We calculate $\Gamma(\varepsilon_{\rm cut})$ over this range
($\varepsilon_{\rm cut}\in\left[\varepsilon^{<}_{\rm cut},\varepsilon^{>}_{\rm
cut}\right]$); this function is made continuous using linear interpolation.
Next, we consider the statistics of $\Gamma(\varepsilon_{\rm cut})$ via the
mean and variance,

$\begin{split}\langle\Gamma\rangle&\equiv\frac{1}{\varepsilon^{>}_{\rm
cut}-\varepsilon^{<}_{\rm cut}}\int_{\varepsilon^{<}_{\rm
cut}}^{\varepsilon^{>}_{\rm cut}}{\rm d}\varepsilon_{\rm
cut}\,\Gamma(\varepsilon_{\rm cut})\\\
\langle\delta\Gamma^{2}\rangle&\equiv\frac{1}{\varepsilon^{>}_{\rm
cut}-\varepsilon^{<}_{\rm cut}}\int_{\varepsilon^{<}_{\rm
cut}}^{\varepsilon^{>}_{\rm cut}}{\rm d}\varepsilon_{\rm
cut}\,\Gamma^{2}(\varepsilon_{\rm cut})-\langle\Gamma\rangle^{2}\end{split}$
(S47)

Using the variance in the fitting parameter as a measure of sensitivity to
choice of $\varepsilon_{\rm cut}$, we consider the choice of $\varepsilon_{\rm
cut}^{>}$ that minimizes the error $\sqrt{|\langle\delta\Gamma^{2}\rangle|}$
(at fixed $\varepsilon_{\rm cut}^{<}=0.02$); the mean fit
$\langle\Gamma\rangle$ and error are shown as a function of $\varepsilon_{\rm
cut}^{>}$ in Fig. S3(a). For DG, we find that $\varepsilon_{\rm
cut}^{>}\approx 0.04$ seems to be a reasonable upper limit on the fitting
range, with values of the quadratic fitting parameter $\Gamma^{\rm DG}_{\rm
shear}$ saturating beyond this point. We adopt the mean values of the fitting
parameter over this interval, $\langle\Gamma\rangle$, for both data sets as
the accepted results, leading to estimates of the non-affinity of $\Gamma^{\rm
DG}_{\rm ext}\approx 0.19\pm 0.01$ and $\Gamma^{\rm DG}_{\rm shear}\approx
0.15\pm 0.01$. As shown in Fig. S3(b), these accepted values seem to fit the
data quite well over a broad range of $\varepsilon\lesssim 0.1$.

Figure S3:  (a) Non-affinity measurements, as characterized by the mean
$\langle\Gamma\rangle$ and standard deviation of the fitting parameter
$\Gamma(\varepsilon_{\rm cut})$, where the cut-off strain for the fit,
$\varepsilon_{\rm cut}$ is sampled over the interval $[0.02,\varepsilon_{\rm
cut}^{>}]$, as a function of the upper limit of the cut-off strain,
$\varepsilon_{\rm cut}^{>}$. Non-affinity parameter fit statistics are shown
in orange for extensional strain and blue for shear strain. (b) Non-affine
response $\langle|\mathbf{u}|^{2}\rangle/\ell_{0}^{2}$ measured along the
simulated extensional (orange) and shear deformation paths (blue) are shown
from SCF calculations (points) and from the mean quadratic fit (solid curves).
(c) The non-affine response of the three SCF-simulated deformation paths
passing through the three triclinic crystal grains from experimental data are
shown as a function of total strain $\sqrt{({\rm
tr}\,\hat{\bm{\varepsilon}}^{2})/2}$ (solid curves), focusing on the low-
strain regime. By comparison, linear response predictions obtained from
weighted combinations of response to extensional and shear deformation paths
are shown as dashed curves. Deviation of SCF results from the linear response
theory suggest that the linear regime is valid for strains up to $\lesssim
5\%$.

###  S5.2 Estimating departure from linear response

Using symmetry arguments, we showed that the linear response can be
represented in terms of a linear combination of extensional $\varepsilon_{\rm
ext}$ and shear $\varepsilon_{\rm shear}$ strains. However, departures from
this approximation are expected (and seen) to occur for larger deformations.
We can attribute departures from the linear theory presented here to several
factors:

  * •

Higher-order corrections to the fixed-volume constraint ${\rm
tr}\,\bm{\varepsilon}=0+\mathcal{O}({\rm tr}\,\bm{\varepsilon}^{2})$ (see e.g.
[8]).

  * •

Increasing importance of higher-order strain couplings (e.g. cubic terms of
the form
$\Gamma_{(3)}^{ijklmn}\varepsilon_{ij}\varepsilon_{kl}\varepsilon_{mn}$)
resulting in “mixing” between strain terms that prevents factoring into powers
of $\varepsilon_{\rm ext}$ and $\varepsilon_{\rm shear}$.

  * •

Dilational modes (i.e. ${\rm tr}\,\bm{\varepsilon}$) in SCF calculations and
experiments that, upon relaxation to thermodynamic equilibrium, spontaneously
break network symmetry (i.e. those symmetries ensured by space group of the
original network), leading to non-affine deformations.

The first two items listed above require a broader non-linear response
extension to the work presented here. The third item would represent a
breakdown of a key assumption, namely the stability of network phases (with
corresponding space group symmetries intact) under uniform dilations of space.
Non-affine response under dilations can be potentially incorporated into a
mechanical network model by including forces that scale differently under re-
scaling transformations of space; these may, for example, lead to “buckling”
deformations of the network (sometimes seen in stiff filament networks [9]).
In our example system of block copolymer melts, the assumed incompressibility
of the melt requires that dilations be accompanied by a flux of molecules
into/out of the network, breaking the assumed fixed mass constraint.

To probe the extent of the linear response regime, we examine the deviation of
SCF calculations from the expectations of our linear response theory.
Equilibrium calculations of the extensional and shear strain deformations of
DG predict that the non-affine response to extensional strain is larger than
the response to shear strain, $\Gamma_{\rm ext}^{\rm DG}>\Gamma_{\rm
shear}^{\rm DG}$. For all other deformation paths, we find it useful to turn
to the “polar projection” of the deformation paths, characterized by their
radius $\sqrt{({\rm tr}\,\hat{\bm{\varepsilon}}^{2})/2}$ and angle $\alpha$,
given by $\tan\alpha=\varepsilon_{\rm shear}/\varepsilon_{\rm ext}$. For this
analysis, we compute $\alpha$ from the first non-zero strain points along each
deformation path. For the deformation paths corresponding to the three DG
deformation paths (labeled $\\#1$, $\\#2$, and $\\#3$), we find
$\alpha_{1}\approx 66^{\circ}$, $\alpha_{2}\approx 49^{\circ}$, and
$\alpha_{3}\approx 58^{\circ}$. In these polar coordinates, our linear
response theory predicts that the non-affine response of each of these paths
can be computed from the extremal paths via

$\frac{\langle|\mathbf{u}|^{2}\rangle}{\ell_{0}^{2}}=\left(\Gamma^{\rm
DG}_{\rm ext}\cos^{2}\alpha+\Gamma^{\rm DG}_{\rm
shear}\sin^{2}\alpha\right)\frac{{\rm
tr}\,\hat{\bm{\varepsilon}}^{2}}{2}\equiv\Gamma^{\rm DG}_{\alpha}\,\frac{{\rm
tr}\,\hat{\bm{\varepsilon}}^{2}}{2}\,,$ (S48)

where the polar angle $\alpha$ weights the two non-affinity parameters in the
calculation of the non-affinity of for each path. Since
$0^{\circ}<\alpha_{2}<\alpha_{3}<\alpha_{1}<90^{\circ}$, as shown in Fig. S2,
we expect that $\Gamma^{\rm DG}_{2}>\Gamma^{\rm DG}_{3}>\Gamma^{\rm DG}_{1}$.
As shown in Fig. S3(c), this is indeed the case for small values of
$\sqrt{({\rm tr}\,\hat{\bm{\varepsilon}}^{2})/2}$. Around $\sqrt{({\rm
tr}\,\hat{\bm{\varepsilon}}^{2})/2}\approx 0.05$, path $\\#3$ starts
overtaking path $\\#2$ in scale of non-affine response, fully overtaking path
$\\#2$ around $\sqrt{({\rm tr}\,\hat{\bm{\varepsilon}}^{2})/2}\lesssim 0.1$.
Moreover, the predicted linear, non-affine response shown by the dashed curves
in Fig. S3(c), seems to follow the SCF data well for $\sqrt{({\rm
tr}\,\hat{\bm{\varepsilon}}^{2})/2}\approx 0.05$. Thus, we estimate that the
linear response regime is well-characterized for strain values
$\varepsilon\lesssim 5\%$, with increasing non-linear response for higher
strain values; as shown in Fig. 3(c) of the main text, non-linear effects
rapidly become important for $\varepsilon\gtrsim 10\%$.

##  S6 Tomographic imaging of PS-PDMS

The polystyrene-b-polydimethylsiloxane (PS-PDMS) diblock copolymer was
synthesized by sequential anionic polymerization of styrene and
hexamethylcyclotrisiloxane. For the DG samples, the average molecular weight
is 43.5 $\text{kg mol}^{-1}$ for PS and 29.0 $\text{kg mol}^{-1}$ for PDMS,
with an overall composition of 40% PDMS (by volume) and a molecular weight
polydispersity index of 1.04. The sample studied was cast slowly (over the
course of one week) from a 10 wt% solution (2 ml) in toluene. After drying,
the sample was heated to $60^{\circ}$C for 3 days in a vacuum in order to
remove any residual solvent. For the DD samples, the average molecular weight
is 51 $\text{kg mol}^{-1}$ for PS and 35 $\text{kg mol}^{-1}$ for PDMS, with
an overall composition of 42% PDMS (by volume) and a molecular weight
polydispersity index of 1.05. The sample studied was slowly cast over 2 weeks
from a 5 wt % solution (2 mL) in chloroform without further thermal annealing.

Given that the samples were heated to a temperature less than the glass
temperature for polystyrene, the experiments reflect a structure that is not
completely relaxed. This may account for some of the discrepancy between the
experimental data and SCF results and calls for a more in-depth study of
structural relaxation under different processing conditions, including
annealing protocols.

Data acquisition was performed using “slice-and-view” scanning electron
microscopy (SVSEM), as detailed in ref. [10]. A gallium-ion (Ga+) beam was
used to mill each sample surface. Imaging of the free surface was performed
using a 1-keV electron beam, enabling 2D density maps of PDMS relative to PS
with a planar resolution of $\sim 3{\rm nm}/{\rm pixel}$. After imaging the
free surface, a new sample surface was milled, resulting in a z-stack of 2D
SEM images, where each z-slice was $\sim 3{\rm nm}$ thick. Details of the
location of the reconstructions are provided upon reasonable request. 3D
reconstruction of the density fields and skeletal graphs were then performed,
following procedures in refs. [7, 10]. In brief, the 3D density fields
obtained from SVSEM were binarized based on block volume fraction of the
sample and smoothed using a Gaussian filter. Initial skeletonization was
performed using a watershed algorithm; extraction of the skeletal network
involved iterative merging of voxels in the skeletal set, until graphs with
straight edges and 3- and 4-coordinated nodes were obtained. Finally, a graph
optimization was performed, wherein node positions were varied until a cost
functional, given by the integral of the PDMS volume fraction along every
strut in the graph was maximized. This helps ensure that the skeletal graph is
a good representation of the “center” of the PDMS domains.

We additionally employed the method outlined in ref. [10] to extract Bragg
peaks from the 3D Fourier transform of the density field data. Since deformed
crystals transform the positions of the expected Bragg peaks of the
corresponding cubic crystal, we used these data to extract affine deformation
matrices $\bm{\Lambda}$ for each of the samples. As $\bm{\Lambda}$ extracted
in this way incorporates isometries of Euclidean space, which are irrelevant
for describing deformations of the unit cell, we use a QR-decomposition to
extract the triangular matrix part of the affine deformation, allowing for
representation of $\bm{\Lambda}$ in the form given by Eq. S40. The extracted
deformation matrices for each of the samples is listed below:

$\displaystyle\bm{\Lambda}^{\rm DG}_{1}$
$\displaystyle=\left(\begin{array}[]{ccc}1.033&0&0\\\ 0.035&1.065&0\\\
-0.079&0.249&0.909\end{array}\right)$ (S52) $\displaystyle\bm{\Lambda}^{\rm
DG}_{2}$ $\displaystyle=\left(\begin{array}[]{ccc}1.071&0&0\\\
0.201&0.852&0\\\ 0.018&0.091&1.097\end{array}\right)$ (S56)
$\displaystyle\bm{\Lambda}^{\rm DG}_{3}$
$\displaystyle=\left(\begin{array}[]{ccc}1.015&0&0\\\ -0.119&1.137&0\\\
0.266&-0.056&0.866\end{array}\right)$ (S60) $\displaystyle\bm{\Lambda}^{\rm
DD}_{1}$ $\displaystyle=\left(\begin{array}[]{ccc}1.041&0&0\\\
-0.078&0.929&0\\\ -0.023&-0.003&1.034\end{array}\right)$ (S64)
$\displaystyle\bm{\Lambda}^{\rm DD}_{2}$
$\displaystyle=\left(\begin{array}[]{ccc}0.994&0&0\\\ -0.006&1.165&0\\\
-0.108&0.113&0.863\end{array}\right)$ (S68)

In order to extract the non-affine displacement, we used the same fitting
method employed in SCF calculations. Since the 3D-reconstructions of each
grain contains several unit cells, each unit cell was considered an
independent measurement of non-affine displacement, allowing for statistical
estimates of error on each non-affinity measurement. We measured the non-
affine displacements of each node and averaged over identical nodes in
adjacent unit cells to obtain a single-node average non-affine displacement
measure $\overline{\left|\mathbf{u}^{(m)}\right|^{2}}$ for node $m$, as well
as a corresponding standard deviation for the single node, $\sigma^{(m)}$. The
non-affine displacement averaged over a unit cell,
$\langle|\mathbf{u}|^{2}\rangle$, was then found by averaging over
$\overline{\left|\mathbf{u}^{(m)}\right|^{2}}$, weighted by each node’s
individual standard deviation, i.e.

$\begin{split}\langle|\mathbf{u}|^{2}\rangle&=\frac{1}{\sigma_{\rm
tot}^{2}}\sum_{m}\frac{\overline{\left|\mathbf{u}^{(m)}\right|^{2}}}{(\sigma^{(m)})^{2}}\\\
\sigma^{-2}_{\rm tot}&=\sum_{m}(\sigma^{(m)})^{-2}\end{split}$ (S69)

where $\sigma_{\rm tot}$ is the overall standard deviation of the averaged
measurements. Note that this same method of extracting data from multiple unit
cells within each sample was used to calculate error on strut length variance
and node volume variance.

###  S6.1 Comparison between DD experiments, SCF computations, and LNT

Here, we report comparisons between measurements of DD grains, simulated
deformation paths in SCF, and predictions from LNT. The results are summarized
in Fig. S4. As shown in Fig. S4(a), while the two grains have significantly
different deformations, their respective simulated deformation paths have very
similar trajectories in the linear regime, making it difficult to
differentiate their linear response behavior. Indeed, both the SCF and LNT
predictions show nearly-overlapping results for the non-affine displacement
$\langle|\mathbf{u}|^{2}\rangle/\ell_{0}^{2}$ (shown in Fig. S4(b)) and the
r.m.s. length variations
$\sqrt{\langle\Delta\ell^{2}\rangle}/\langle\ell\rangle$ (shown in Fig.
S4(c)). Note that for sample #1, we observed a slight gradient to the lattice
spacing, so to minimize the effect of this gradient on the calculation of the
non-affine displacement, we excluded points a distance greater than
$2\ell_{0}$ from the center of the 3D reconstruction; inclusion of further
points increased the measured value of $\langle|\mathbf{u}|^{2}\rangle$ by an
order of magnitude. Moreover, due to relatively large number of 3-coordinated
nodes observed in sample #2, we set the cutoff on minimum graph edge length
during the reconstruction to be large enough to yield only 4-coordinated
nodes. As discussed in Sec. S8 below, deformations of DD, particularly those
that lead to large extensional strain, result in splitting of 4-coordinated
nodes into pairs of 3-coordinated nodes. Since sample #2 has a considerably
larger value of $\varepsilon_{\rm ext}$ than sample #1, the increase in node
splitting is consistent with this observation.

Figure S4:  (a) Deformation paths #1 and #2 for DD (solid curves), projected
onto the space spanned by the linear combination of extensional strain
$\varepsilon_{\rm ext}$ and shear strain $\varepsilon_{\rm shear}$. These
curves pass through points representing the two triclinic DD grains from
experimental data. (b) Non-affine displacement response
$\langle|\mathbf{u}|^{2}\rangle/\ell_{0}^{2}$ of experimental data (points),
SCF calculations (curves), and LNT predictions (dashed curves). (c) r.m.s.
strut length variations for the same data sets.

##  S7 Nodal volume measurements in DG

Calculations of nodal volume for 3D reconstructions of experiments and SCF
calculations were performed by associating each voxel in the 3D A-block
density field to the nearest node, according to Euclidean distance. Given the
volume $\delta V$ of each voxel and local A-block volume fraction $\phi_{v}$
at voxel $v$, the total A-block volume $V^{(m)}$ associated to node $m$ is
given by

$V^{(m)}=\delta V\,\sum_{v\in\mathcal{N}_{m}}\phi_{v}\,,$ (S70)

where $\mathcal{N}_{m}$ is the set of voxels associated to node $m$. For the
experimental measurements, we only calculated the volumes of nodes that are a
suitable distance away from the sample boundary, to reduce error. This
amounted to nodes that lie close enough to the center of the data set; we
found that only calculating nodal volume for nodes within a distance of
$4\ell_{0}$ of the center was sufficient for both reducing boundary error and
obtaining a suitable number of corresponding nodes to obtain error estimates
via the method outlined in Sec. S6.

For LNT comparisons, we modeled the volume of each node as the volume of the
Voronoi cell generated by the node. In the cubic configuration, all nodes are
equivalent, so each Voronoi cell is equivalent. In each deformed
configuration, node displacements lead to variations in Voronoi cell volumes.
We find that the r.m.s. volume variance computed in this way is consistently
less than the volume variance computed from SCF and experiment. This, in part,
is due to variations in network domain thickness that are not captured in the
Voronoi model. Moreover, the polygonal faces of Voronoi cells only roughly
approximate the domain boundaries of network phases, which are better
approximated as minimal surfaces (see e.g. [11, 12]).

Figure S5:  (a) Single DG skeletal graph showing the 8 nodal regions and pairs
of equivalent nodes represented by the same color. (b)-(d) Volume variation
with respect to the mean $\Delta V/\langle V\rangle$ for DG grains #1 - #3,
where red and blue data points represent contributions from the two separate
networks. Subsets of identical nodes are labeled by (i)-(iv).

###  S7.1 Nodal polydispersity in DG

While all nodal regions within an undeformed, cubic DG, DD, or DP crystal are
equivalent and can be identified as “mesoatomic” building blocks, symmetry-
breaking deformations can result in different deformations amongst the various
nodes. Our analysis of the trihedral graphs, as reported in Eq. S16, shows
that the response of a nodal region depends on the orientation of the node. In
the DG skeletal graphs, each graph has 8 nodes, and of those nodes, pairs
share the same orientation, as shown in Fig. S5(a). These orientation-sharing
nodes additionally have the same neighbors: they are identical under a half-
translation of the unit cell (i.e., by $\left(1,1,1\right)/2$). Since the
second skeletal graph is obtained by a translation and inversion of the first,
it shares the same node orientations; thus, out of 16 nodes, there are 4
equivalent subgroups. Moreover, since the triclinic deformations considered
here preserve inversion symmetry, these 4 nodes subsets should be preserved in
triclinic DG networks. One metric of nodal geometry is provided by the nodal
volume; variance in nodal volume indicates some degree of deformation (note
that nodes may deform without changing volume, however). The variation in
nodal volume amongst the 8 different nodes from each skeletal graph (denoted
by red and blue points) obtained from the 3 experimental data sets is shown in
Fig. S5(b)-(d). We see that patterns in the range of volume variations for
each of the 8 indexed nodes consistently fall into 4 equivalent nodes between
the two networks, labeled (i)-(iv). We find that the two networks in each data
set are nearly identical according to volume variance, with slight
differences, notably for sample #2. While further studies are needed to better
quantify nodal shape variations, we can nonetheless identify persistent sets
of identical nodes within DG networks. By contrast, DD all nodes are
equivalent under translation and inversion and do not have this same mode of
nodal polydispersity. Instead, we find that DD nodes can undergo “splitting”
into pairs of 3-coordinated nodes, as described in Sec. S8.

##  S8 Node splitting in double diamond (DD)

Figure S6:  (a) DD cubic unit cell from SCF calculations, along with skeletal
graphs. (b) Relaxed DD structure after imposed tetragonal deformation
($\varepsilon_{\rm ext}\approx 16\%$), showing that the 4-coordinated split
into pairs of 3-coordinated nodes. (c) View of 8 unit cells of the tetragonal-
DD network, showing the full structure of the double-network.

In order to calculate the non-affinity of the deformed DD networks, depicted
in Fig. S6(a), we imposed a constraint on the edge length cutoff for the
skeletonization algorithm that ensured the skeleton graphs consist of
4-coordinated nodes. However, if we relax the edge cutoff to standard values
used for constructing DG networks, we find that above a modest value of
extensional strain, $\varepsilon_{\rm ext}\gtrsim 5\%$, the 4-coordinated
nodes split into pairs of 3-coordinated nodes. As the strain is increased, the
length of the strut joining the new nodes increases, eventually reaching a
length that is about $20\%$ of the initial strut length of the DD near strain
values of $\varepsilon_{\rm ext}\sim 16\%$, depiced in Fig. S6(b),(c). Since
the new nodes are close to planar and have strut angles close to
$120^{\circ}$, we expect that such networks are well-described by LNT. As the
equilibrium networks of LNT minimize their total length at constant volume,
this suggests that such “topological transitions” of the graph structure can
be rationalized by the mathematics of length-minimizing, “Steiner” networks.
Recent work [3] has shown that the skeletal network corresponding to the DG,
the srs network (named for the strontium silicide crystal structure), is the
globally length-minimizing, triply-periodic Steiner network. We recognize that
the skeletal networks of the tetragonally-deformed DD closely resemble the ths
network (named for the thorium silicide crystal structure), which include
Steiner networks with lengths slightly larger than the srs. The ths networks,
also called (10,3)-b [13], consist of planar “zig-zag” strut patterns that
resemble the skeletal structure of the O70 (orthorhombic $Fddd$ phase) [14,
15]. Interestingly, the SCF calculations showed a continuous, monotonic free-
energy _decrease_ from the cubic state, suggesting that under conditions of
fixed volume (i.e. fixed total mass), the cubic DD is weakly unstable to
tetragonal deformations; at $\varepsilon\approx 16\%$, we found a slight free-
energy decrease of $\sim-0.1\%$ relative to the cubic-DD. SCF calculations
starting for the tetragonal-DD and deforming back to the cubic-DD showed the
same free-energy change, indicating that this transition is indeed continuous.
We leave the exploration of the node-splitting mechanics and thermodynamics
for future study.

## References

  * Ivanov and Tuzhilin [1994] A. Ivanov and A. A. Tuzhilin, Minimal Networks: The Steiner Problem and Its Generalizations (CRC Press, 1994).
  * Alex and Grosse-Brauckmann [2019] J. Alex and K. Grosse-Brauckmann, Periodic networks of fixed degree minimizing length (2019), arXiv:1911.01792 [math.CO] .
  * Alex and Grosse-Brauckmann [2023] J. Alex and K. Grosse-Brauckmann, Periodic steiner networks minimizing length, Discrete & Computational Geometry 10.1007/s00454-023-00576-z (2023).
  * Marko and Siggia [1995] J. F. Marko and E. D. Siggia, Stretching DNA, Macromolecules 28, 8759 (1995).
  * Rubinstein and Colby [2003] M. Rubinstein and R. Colby,  _Polymer Physics_ (Oxford University Press, 2003).
  * Arora  _et al._ [2016] A. Arora, J. Qin, D. C. Morse, K. T. Delaney, G. H. Fredrickson, F. S. Bates, and K. D. Dorfman, Broadly Accessible Self-Consistent Field Theory for Block Polymer Materials Discovery, Macromolecules 49, 4675 (2016).
  * Prasad  _et al._ [2018] I. Prasad, H. Jinnai, R.-M. Ho, E. L. Thomas, and G. M. Grason, Anatomy of triply-periodic network assemblies: characterizing skeletal and inter-domain surface geometry of block copolymer gyroids, Soft Matter 14, 3612 (2018).
  * Warner and Terentjev [2007] M. Warner and E. Terentjev,  _Liquid Crystal Elastomers_ , International Series of Monographs on Physics (OUP Oxford, 2007).
  * Conti and MacKintosh [2009] E. Conti and F. C. MacKintosh, Cross-linked networks of stiff filaments exhibit negative normal stress, Phys. Rev. Lett. 102, 088102 (2009).
  * Feng  _et al._ [2019] X. Feng, C. J. Burke, M. Zhou, H. Guo, K. Yang, A. Reddy, I. Prasad, R.-M. Ho, A. Avgeropoulos, G. M. Grason, and E. L. Thomas, Seeing mesoatomic distortions in soft-matter crystals of a double-gyroid block copolymer, Nature 575, 175–179 (2019).
  * Reddy  _et al._ [2022] A. Reddy, M. Dimitriyev, and G. Grason, Medial packing and elastic asymmetry stabilize the double-gyroid in block copolymers, Nature Communications 13, https://doi.org/10.1038/s41467-022-30343-2 (2022).
  * Dimitriyev  _et al._ [2023] M. S. Dimitriyev, A. Reddy, and G. M. Grason, Medial packing, frustration, and competing network phases in strongly segregated block copolymers, Macromolecules 56, 7184 (2023).
  * Wells [1977] A. F. Wells,  _Three Dimensional Nets and Polyhedra_ (Wiley, New York, 1977).
  * Bailey  _et al._ [2002] T. S. Bailey, C. M. Hardy, T. H. Epps, and F. S. Bates, A noncubic triply periodic network morphology in poly(isoprene-b-styrene-b-ethylene oxide) triblock copolymers, Macromolecules 35, 7007 (2002).
  * Tyler and Morse [2005] C. A. Tyler and D. C. Morse, Orthorhombic $fddd$ network in triblock and diblock copolymer melts, Phys. Rev. Lett. 94, 208302 (2005).