# Elastic Response of Wire Frame Glasses. I. Two Dimensional Model

David A. King<EMAIL_ADDRESS>Cavendish Laboratory, University of Cambridge,
J. J. Thomson Ave., Cambridge CB3 0HE, UK  Masao Doi  Centre of Soft Matter
and its Applications, Beihang University, Beijing 100191, China  Erika Eiser
Cavendish Laboratory, University of Cambridge, J. J. Thomson Ave., Cambridge
CB3 0HE, UK

###### Abstract

We study the elastic response of concentrated suspensions of rigid wire frame
particles to a step strain. These particles are constructed from infinitely
thin, rigid rods of length $L$. We specifically compare straight rod-like
particles to bent and branched wire frames. In dense suspensions the wire
frames are frozen in a disordered state by the topological entanglements
between their arms. We present a simple, geometric method to find the scaling
of the elastic stress with concentration in these glassy systems. We apply
this method to a simple 2D model system where a test particle is placed on a
plane and constrained by a random distribution of points with number density
$\nu$. Two striking differences between wire frame and rod suspensions are
found: 1) The linear elasticity per particle for wire frames is very large,
scaling like $\nu^{2}L^{4}$, whereas for rods it much smaller and independent
of concentration. 2) Rods always shear thin but wire frames shear harden for
densities less than $\sim\sqrt{K/k_{B}TL^{4}}$, where $K$ is the bending
modulus of the particles. The deformation of wire frames is found to be
important even for small strains, with the proportion of deformed particles at
a particular strain, $\gamma$, being given by $(\nu L^{2})^{2}\gamma^{2}$. Our
results agree well with a simple simulation of the 2D system.

††preprint: AIP/123-QED

##  I Introduction

Concentrated suspensions of rigid macro-molecules are a prototypical soft
matter system, constituting a rich field of study for both experimentalists
and theoreticians. The great variety of non-Newtonian flow behaviours
displayed by these systems has been a particular interest. Being able to
predict and engineer these behaviours based on the shape of the suspended
particles is a problem of broad scientific an technological importance.

The interactions between the particles are generally what lead to the
pronounced non-Newtonian behaviour. The strongest effects are observed for
concentrated suspensions, where the number density of particles, $\rho$, is
large such that $\rho L^{3}\gg 1$, with $L$ being a typical length-scale of
the particles. Hydrodynamic and excluded volume interactions will always be
present and have an effect on the rheology, but for certain particle shapes
and concentrations a different kind of interaction dominates. These are
‘kinetic’ or ‘topological’ constraints, which originate from the fundamental
property that no two particles can cross. This is the case when the particles
are very long compared to their width and the particles are rarely in close
contact, such that $\rho V_{exc}\ll 1$, where $V_{exc}$ is the excluded volume
of the particles. As the excluded volume is irrelevant, all equilibrium
properties are known; the system is equivalent to an ideal gas. On the other
hand, the non-equilibrium properties, such as the rheology, are severely
affected by the non-crossing condition. This is often referred to as the
‘entanglement effect’.

We will focus on how kinetic constraints affect the rheology of the suspension
depending on the shape of suspended particles. The concentration range
considered is, $1/L^{3}\ll\rho\ll 1/V_{exc}$. This ensures that the kinetic
constraints are the dominant interaction and that there is no long range
orientational order in the system. Each of the particle shapes we will
consider is constructed from infinitely thin, rigid rods. We consider the
joints between the rods to be effectively rigid and unaffected by thermal
fluctuations. Such particle shapes can be divided into three different
classes; straight, bent and branched. Straight particles are simply rods, see
Figs.(1a) & (2c). Bent particles are two rods joined at some angle, see
Figs.(1b) & (2a). Branched particles are those where three or more component
rods meet at a point. These can be planar stars as in Fig.(2a) or fully three
dimensional shapes, for instance the shapes in Fig.(2d) with an extra rod
piercing the plane. We refer to the bent or branched particles as ‘wire
frames’.

Wire frame particles can be realised practically using DNA origami
techniquesSeeman (1982); Rothemund (2006). Arms made of double stranded DNA
sequences can be engineered so that they rigidly link at precise angles to
form ‘nanostars’ Bomboi  _et al._ (2019); Biffi  _et al._ (2013); Xing  _et
al._ (2018). Modelling the double stranded DNA arms as infinitely thin,
perfectly rigid rods is reasonable, due to their large aspect ratio ($\sim
20$) and stiffness (persistence length $\sim 390$Å)Gross  _et al._ (2011).
Being able to predict the different rheological responses depending on the
shape of the nano-stars is very important for the design of these systems as
functional nano-materials.

The kinetic constraints in these systems can also be responsible for a glass
transition. Suspensions of 3D crosses, particles consisting of three mutually
perpendicular, infinitely thin rods joined at their midpoints, have been shown
to have glassy dynamical behaviour in the range of densities we are
consideringVan Ketel, Das, and Frenkel (2005). In this system the
translational diffusion coefficient goes to zero with an exponential
dependence on the density. The majority of particles are rendered effectively
immobile by the kinetic constraints imposed by their surroundings and the
system is frozen in an isotropic, disordered state.

The geometry of the suspended particles must influence this behaviour. In
Fig.(1) we sketch two examples of the systems we are considering. A dense
suspension of rods is shown in Fig.(1a) and a dense glassy state of L-shaped
wire frames is shown in Fig.(1b). In each panel a particular test particle is
highlighted in blue. The rod may still diffuse in this concentration regime by
the reptation mechanism Doi (1975); Doi and Edwards (1986). However, the wire
frames cannot. This is because the reptation process relies on the rods’
ability to diffuse freely along their length even though their transverse
motion is severely hindered. If a wire frame diffuses along the length of one
of its constituent rods, then any one of the others will quickly become
entangled, blocking its motion.

Figure 1:  Sketches of the systems considered. (a) A dense suspension of rod
like particles. A test rod is shown in cyan. The motion of the test rod is
restricted by the surrounding red rods, but it can still diffuse along its
length, so the suspension is still in a fluid state. (b) L-shaped particles in
dense suspension. The cyan test L-shape is trapped by the surrounding red
L-shapes. If it moves along the length of one of its legs, the other becomes
entangled. The system is frozen in a glassy state.

It is interesting to compare the rheology of glassy states of wire frames,
Fig.(1b), to dense fluid states of rods, Fig.(1a). For flexible polymers, it
is well known that branched star polymers display significantly different flow
behaviour from linear chains, with extremely slow stress relaxation Doi and
Edwards (1986); Gennes (1975); Pearson and Helfand (1984); Doi and Kuzuu
(1980); Milner and McLeish (1997). This is as a result of the branching making
reptation impossible. It is reasonable then to expect an equally pronounced
difference in behaviour between rigid branched particles and rods. In fact, it
has been observed in simulationsHeine, Petersen, and Grest (2010) that the
viscosity for dense suspensions of branched 3D crosses is significantly larger
than for rods and that the viscosity has a much stronger dependence on
concentration Petersen, Lane, and Grest (2010).

Evidently, we cannot use the standard approach based on the tube model Doi and
Edwards (1978, 1986), where the reptation mechanism is used in a hydrodynamic
formalism. Instead we must use a new approach which determines the rheological
properties of these glassy states directly from the constraints on each
particle. In this two part series of papers, we restrict our attention to a
simple question which gives useful insight into the problem. We do not address
in detail the time (or frequency) dependence of the rheological response here,
but rather focus on its instantaneous magnitude. Our hope is that the method
presented here will form a useful basis for more detailed studies in the
future.

We begin with some simple definitions before stating the problem explicitly.
The rheological properties of a suspension are expressed by the constitutive
equation relating the stress tensor to the applied rate strain tensor, $E(t)$.
The stress will consist of two parts; the elastic stress, which we call
$\sigma$, and the viscous stress, $\sigma_{V}$. The viscous stress is related
to the energy dissipation in the system and can generally be written as a
function of time, $\sigma_{V}(t)=\mathcal{V}:E(t)$, where the 4th rank tensor
$\mathcal{V}$ can be a function of the deformation history of the suspension.
The elastic stress is related to the change in the free energy per unit
volume, $F$, by the virtual work principle;

$\delta F=\sigma:E\delta t,$ (1)

where $\delta F\equiv F(E)-F(0)$ is the change in the free energy, calculated
as the difference between the free energy in a system deformed by $E$ and that
in the undeformed state. If the applied strain is small enough, the elastic
stress will depend on it linearly,

$\sigma(t)=\int_{-\infty}^{t}dt^{\prime}\
\mathcal{E}(t-t^{\prime}):E(t^{\prime}),$ (2)

where the time dependent, fourth rank tensor $\mathcal{E}$ is the elastic
response function.

To probe the rheology of these systems, we consider the stress response to a
step strain. In this case, the rate of strain tensor is given by,

$E(t)=\kappa\delta(t),$ (3)

where we have defined the strain tensor, $\kappa$. We always take the strain
to be simple $x-y$ shear for which $\kappa$ has only one non vanishing
component, $\kappa_{xy}\equiv\gamma$. The extension to other strains is
straightforward. For step strain, the viscous stress must also be proportional
to a delta function, and so can be taken to vanish for all practical
timescales.

The virtual work principle for the elastic stress can be re-written,

$\sigma=\frac{\partial F}{\partial\kappa},$ (4)

and for small strains will take the form,

$\sigma(t)=\mathcal{E}(t):\kappa=\gamma\mathcal{E}_{xy}(t),$ (5)

hence the choice of step strain allows us to directly probe the elastic
response function.

Typically the stress will start at a particular value, $S_{0}$, and then
decay. The initial timescale for this decay, $\tau_{0}$, is expected to be the
same for all concentrations. This is because immediately after the
deformation, the particles are free to move without hindrance from the
surroundings. After a time, $\tau_{\text{int}}$, the particles have moved
sufficiently to begin interacting with their neighbours and the decay
timescale will increase to $\tau(\rho)$, due to the entanglement effect. This
new decay rate is a function of the particle density, since for denser systems
the entanglement effect is stronger. Therefore a reasonable approximation for
the time dependent elastic stress is,

$\sigma(t)=\begin{cases}S_{0}\ e^{-t/\tau_{0}}\ ,&\text{for}\
t<\tau_{\text{int}}\\\ S_{0}\
e^{\tau_{\text{int}}(1/\tau(\rho)-1/\tau_{0})}e^{-t/\tau(\rho)}\ ,&\text{for}\
t>\tau_{\text{int}}\end{cases}.$ (6)

In concentrated suspensions, it is known that $\tau_{0}\ll\tau$ and hence we
can ignore the initial rapid decay of the stress and take $\sigma(t)$ to be
approximated by,

$\sigma(t)\approx S_{0}\
e^{-\tau_{\text{int}}/\tau_{0}}e^{-t/\tau(\rho)}\equiv\sigma_{0}(\rho)e^{-t/\tau(\rho)},$
(7)

where we have defined $\sigma_{0}(\rho)$ as the initial value of the stress
measured in response to a step strain. Throughout these papers, we refer to
$\sigma_{0}(\rho)$ as simply ‘the elastic stress’.

From this discussion we see that there are two main parts to the stress
response. Its initial magnitude, $\sigma_{0}(\rho)$, and $\tau(\rho)$, the
dominant timescale of the subsequent decay. This timescale is generally very
long, so the elastic stress persists for a long period after the initial
strain. There are many interesting questions relating to this timescale in the
case when reptation is not possible. In particular, is there a critical
concentration at which this timescale diverges? Or is the relaxation a single
or stretched exponential, as could be expected?Edwards and Vilgis (1986).
While we do not address these problems here, we hope that the methods we
introduce may be of use in their resolution.

In this paper (paper I), we present a simple geometric method for determining
the scaling of the elastic stress, $\sigma_{0}(\rho)$. This method is
applicable to a wide range of wire-frame shapes and is outlined in the next
section. In section III we introduce a two dimensional model system to which
we apply this method. We find a striking difference between straight, rod-like
particles and bent or branched particles. In section IV, we extend our
treatment to the non-linear elasticity and show that there is a critical
density, related to the bending modulus of the individual particles, where the
behaviour changes from shear hardening to shear softening. This is in contrast
to the result for rods, where the suspension always shear thins. Our
theoretical results are found to be consistent with those of a simulation of
the 2D system. The purpose of this paper is to introduce the method in the 2D
model, where the calculation can be performed simply and exactly. A full
treatment of the problem in three dimensions will be given in paper II.

##  II General Method: Linear Response

To determine the elastic stress from the virtual work principle (4), we need
the free energy as a function of strain. The free energy is, of course,
$F=U-TS$, where $U$ is the internal energy, $S$ is the entropy, and $T$ is the
temperature. We consider the system in the absence of an external potential,
therefore the strain only changes the internal energy if the particles
themselves are deformed. This will be shown to contribute only to the non-
linear elastic response, which is dealt with in section IV. In this section,
we outline a simple, geometric method for calculating $\sigma_{0}$ from the
entropy change per particle caused by the strain in the linear regime.

In the glassy states, the motion of a test particle is impaired by the other
particles in the system, with its centre of mass confined to a ‘cage’ and only
certain orientations accessible at any given moment. Generally, it is only the
Brownian rotations of the particle which contribute to the stress, not the
motion of its centre of mass. This allows us to treat the centre of mass as
fixed and focus only on the constraints placed on the particle’s orientation.

Let us define the orientation of the particle as $\Gamma$. In two dimensions,
$\Gamma$ is represented by the angle, $\phi$, that the particle makes to a
reference axis and in three dimensions, three Euler angles specify $\Gamma$.
The surrounding particles constrain the test particle to a region in the
configuration space. The size of this region depends on $\Gamma$ and the
configuration of the surrounding particles $C$, and is written
$\Omega(\Gamma,C)$. The test particle will explore all of $\Omega$ over a
timescale of roughly $\tau_{int}$, which is generally very short. Therefore
the particle essentially instantly equilibrates inside $\Omega$. The entropy
associated with this particle then follows from the Boltzmann definition,

$S(\Gamma,C)=k_{B}\log\Omega(\Gamma,C).$ (8)

When the step strain is applied, the surrounding particles will move, and so
the accessible volume of configuration space will instantaneously change from
$\Omega$ to $\tilde{\Omega}(\kappa)$. The change in entropy caused by this
transformation alters the free energy and hence yields an elastic stress in
the system. This method is similar to that introduced by Edwards to study
flexible polymers with topological constraints and crosslinksEdwards (1967,
1969); Edwards and Freed (1969). As the surrounding particles rearrange due to
their Brownian motion, the test particle can explore more of the configuration
space, so $\Omega$ begins to increase and the stress relaxes. We are
interested in the initial value of the elastic stress, $\sigma_{0}$, so we
only need to consider the instantaneous change after the step strain.

If there are $\rho$ particles per unit volume, the free energy per unit volume
immediately after the step strain is,

$F(\kappa)=-k_{B}T\rho\int
d\Gamma\psi(\Gamma)\Big{\langle}\log\tilde{\Omega}(\Gamma,C\
;\kappa)\Big{\rangle},$ (9)

where $\psi(\Gamma)$ is the orientational distribution function, which may be
taken to be uniform since the system in a disordered state, and
$\langle\cdots\rangle$ denotes averaging over all configurations, $C$, of the
constraints. Differentiating this with respect to $\kappa$ will give the
stress immediately after the step strain, $\sigma_{0}$.

The advantage of equation (9) for the free energy is that $\tilde{\Omega}$ can
be calculated from purely geometric considerations, without detailed knowledge
of the dynamics of the particles. In three dimensions the geometry required
will become very complicated, especially when defining the constraints placed
on the test particle in order to calculate $\tilde{\Omega}$. In two
dimensions, however, the geometry is simple enough that the problem may be
solved exactly.

##  III 2D Model and Simulation Details

The calculation of $F(\kappa)$, becomes quite complex in three dimensions.
Therefore, here we first consider a simpler, two dimensional model.

In this model we only consider planar wire frames, as shown in Fig.(2) and
each particle is fixed parallel to the $x$-$y$ plane. The wire frame particles
are infinitely thin so no two lie in the same plane and as such, each particle
forms its own independent, 2D system. The rotation of each particle in the
plane is constrained by many point-like obstacles which it cannot cross, as
shown in Fig.(3a). The number of obstacles per unit area is, $\nu$, and in the
concentrated limit, $\nu L^{2}\gg 1$. These points play the role of the other
wire frame particles in three dimensions, preventing a test wire frame from
freely rotating from one orientation to another. We ignore the translational
part of the entropy and only consider the rotation of the test particle. The
obstacles are also assumed to be fixed and are displaced affinely when the
deformation is applied. This model system can be thought of as a cross section
of a three dimensional system taken in the plane of the test particle. The
obstacle points are then the surrounding particles which intersect this plane.

The particles we consider are all formed of rigid rods of the same length,
$L$, referred to as ‘legs’, which will be indexed by $l$. Each leg has the
same length, $L$, and the legs all meet at the same point which is fixed at
the origin. The direction of the $l$th leg is given by the unit vector,
$\textbf{e}_{l}$, running parallel to it and pointing away from the origin.
The orientation of the particle is given by the angle $\phi$ of the $l=0$ leg
to the $x$-axis, so that its unit vector is,
$\textbf{e}_{0}\equiv(\cos\phi,\sin\phi)$.

Figure 2:  Sketches of the particle shapes considered. (a & b) An example of a
bent particle, an L-shape with arms of length $L$. (c) A rod of length $2L$.
(d) An example of a branched particle, a ‘5-Star’. In each panel the unit
direction vectors of the legs, $\textbf{e}_{l}$, are shown and the orientation
of the $l=0$ leg is set at an angle $\phi$ from the $x$-axis, shown as a
dashed line. In panels (a) and (b), the L-shape is constrained to lie between
the two dotted lines. In (a) the angles, $\eta_{l}$ and $\theta_{l}$, from
each leg to the nearest anti- and clockwise constraint respectively, are
shown. In (b) the angles which the constraints make to the $x$-axis are shown;
$\alpha_{l}$ for the closest anti-clockwise constraint to leg $l$ and
$\beta_{l}$ for the closest clockwise constraint. The points representing the
constraints (see Fig.(3)) are omitted from this sketch for clarity. The angle
from the $l=0$ to $l$th leg is $\Delta_{l}$. Shown in (a), this is $\pi/2$ for
the L-shape. For the 5-star in (d) this is $2l\pi/5$.

Imagine rotating the particle clockwise through a full turn. As the particle
is rotated one particle leg will eventually collide with one of the
constraining points. Let us define the vector $\textbf{a}_{l}$ as the position
vector of the constraint which the $l$th leg first crosses during this
clockwise rotation. In a similar way we define, $\textbf{b}_{l}$ as the
position vector of the corresponding point for the anti-clockwise rotation.
The important feature of these vectors is their angle to the $x$-axis defined
in terms of the components of the vectors as,

$\begin{split}\tan\alpha_{l}=\frac{a^{l}_{y}}{a^{l}_{x}}\ \ \ \text{and}\ \ \
\tan\beta_{l}=\frac{b^{l}_{y}}{b^{l}_{x}}.\end{split}$ (10)

Figure 3:  A sketch of the 2D model system for rods. The rod is shown as a
solid line at an angle $\phi$ from the dashed $x$-axis. Panel (a) shows the
constraining points surrounding the rod with the two points closest to the rod
highlighted in red. These two points constrain the orientation to the range of
angles, $\Omega$. In (b) the step strain has been applied, indicated by the
blue arrows, and the constraining points closest to the rod have moved, with
their new positions shown in the darker, solid colour and their previous
positions in a lighter colour. This changes the range of accessible angles to
$\tilde{\Omega}$.

Further, we define the angles indicated in Fig.(2a),

$\eta_{l}=\phi+\Delta_{l}-\alpha_{l}\ \ \ \text{and}\ \ \
\theta_{l}=\beta_{l}-\phi-\Delta_{l},$ (11)

where $\Delta_{l}$ is the angle between the $l=0$ leg and the $l$th leg. The
particle is constrained by the points which are closest to it, so the range of
accessible angles is,

$\Omega=\min_{l,m}\ \eta_{l}+\theta_{m}.$ (12)

To obtain the free energy as given in equation (9), an average needs to be
taken over the distribution of constraints, which means averaging over all the
angles, $\eta_{l}$ and $\theta_{m}$. The probability distribution,
$P(\eta_{l})$ and $P(\theta_{m})$, can be determined from a simple argument.
For any given area of size $a$, the average number of points inside is the
same, $\nu a$, and independent of the number found in a different area.
Therefore, the points are distributed according to a Poisson distribution.
When the particle rotates through an angle $\eta_{l}$, say, the $l$th leg
sweeps out an area of size $L^{2}\eta_{l}/2$. The probability that there are
no constraining points in this area is, $P_{0}(\eta_{l})=\exp(-\nu
L^{2}\eta_{l}/2)$. By definition the $l$th leg firsts collide with a
constraint after rotating clockwise through an angle $\eta_{l}$, so it follows
that $P(\eta_{l})=-P_{0}^{\prime}(\eta_{l})$. So we find,

$P(\eta_{l})=\frac{\nu L^{2}}{2}e^{-\nu L^{2}\eta_{l}/2},$ (13)

and exactly the same form for $P(\theta_{m})$.

When the shear strain is applied, the constraining points all move and the
angles $\alpha_{l}$, $\beta_{l}$, $\eta_{l}$ and $\theta_{l}$ all change. We
assume the the constraining points are displaced affinely when the shear is
applied, therefore the position vectors, $\textbf{a}_{l}$ and
$\textbf{b}_{l}$, transform according to the rule,

$\textbf{a}_{l}\to\tilde{\textbf{a}_{l}}=(\mathbb{I}+\kappa)\cdot\textbf{a}_{l},$
(14)

with the same expression for $\tilde{\textbf{b}}_{l}$. For the simple $x$-$y$
shear we consider, the new angles can be determined from,

$\tan\tilde{\alpha}_{l}=\frac{a^{l}_{y}}{a^{l}_{x}+\gamma
a^{l}_{y}}=\frac{\tan\alpha_{l}}{1+\gamma\tan\alpha_{l}},$ (15)

and an analogous expression for $\tilde{\beta}_{l}$. We focus on small strains
with $\gamma\ll 1$. Expanding (15) in this limit, we find,

$\tilde{\alpha}_{l}=\alpha_{l}-\gamma\sin^{2}\alpha_{l}+\gamma^{2}\sin^{3}\alpha_{l}\cos\alpha_{l}+\mathcal{O}(\gamma^{3}),$
(16)

and similarly for $\beta_{l}$. The terms to second order in $\gamma$ will
contribute to the linear elasticity. For the non-linear response, more terms
must be retained. The angles, $\tilde{\eta}_{l}$ and $\tilde{\theta}_{l}$ are
then found from,

$\tilde{\eta}_{l}=\phi+\Delta_{l}-\tilde{\alpha}_{l}\ \ \ \text{and}\ \ \
\tilde{\theta}_{l}=\tilde{\beta}_{l}-\phi-\Delta_{l},$ (17)

which leads to the transformed range of accessible angles,

$\tilde{\Omega}=\min_{l,m}\ \tilde{\eta}_{l}+\tilde{\theta}_{m}.$ (18)

This model forms the basis of a simple numerical simulation we use for
comparison to our analytic results. This computes the accessible angle and the
free energy numerically by the following simple procedure. First, a particle
is placed in a random orientation in the plane. Then the obstacles are placed
at random on a circle surrounding the particle; the number of obstacles is
proportional to the density, $\pi\nu L^{2}$. The two closest obstacles to the
particle are then found and their positions determine the accessible angle.
The positions of the obstacles on the circle are then all moved according to
(15). The logarithm of the ratio of the accessible angle before and after the
transformation gives the free energy for that configuration. These steps are
repeated for many different configurations of the particle and obstacles and
the average free energy over these realisations is calculated. This
essentially calculates (9) numerically using a Monte-Carlo integration for the
averages.

###  III.1 Rods

The first shape we consider is a rod of length $2L$, which we can think of as
two legs of length $L$ which are anti-parallel. The unit vectors giving the
orientations of these legs are therefore,
$\textbf{e}_{0}=(\cos\phi,\sin\phi)=-\textbf{e}_{1}$. Due to the symmetry of
the rod and the shear flow considered, it does not matter which leg of the rod
is constrained. This means we only need to consider the angles $\eta_{0}$ and
$\theta_{0}$, from which we will drop the subscripts in this section. Using
equations (16) and (17), we can determine the transformed angles
$\tilde{\eta}$ and $\tilde{\theta}$,

$\tilde{\eta}=\eta+\gamma\sin^{2}(\phi-\eta)-\gamma^{2}\sin^{3}(\phi-\eta)\cos(\phi-\eta),$
(19)
$\tilde{\theta}=\theta-\gamma\sin^{2}(\phi+\theta)+\gamma^{2}\sin^{3}(\phi+\theta)\cos(\phi+\theta).$
(20)

The transformed range of accessible angles is then calculated from
$\tilde{\Omega}=\tilde{\eta}+\tilde{\theta}$,

$\begin{split}&\tilde{\Omega}=\Omega+\gamma\sin(\eta+\theta)\sin(\eta-\theta-2\phi)-\frac{\gamma^{2}}{8}\bigg{[}\sin\big{(}4(\eta-\phi)\big{)}\\\
&+\sin\big{(}4(\theta+\phi)\big{)}-4\cos(\eta-\theta-2\phi)\sin(\eta+\theta)\bigg{]}.\end{split}$
(21)

For convenience let us define the functions, $f(\eta,\theta;\phi)$ and
$g(\eta,\theta;\phi)$ such that,

$\tilde{\Omega}=\Omega(1+\gamma f+\gamma^{2}g/2).$ (22)

To second order in $\gamma$, the entropy of the test particle for this
realisation is,

$S=S_{0}+k_{B}\gamma f+\frac{k_{B}\gamma^{2}}{2}(g-f^{2}),$ (23)

where $S_{0}=k_{B}\log\Omega$ is the entropy of the undeformed state. Now the
free energy to the same order is according to (9),

$\begin{split}F=F_{0}-\frac{k_{B}T}{2\pi}\rho\int_{0}^{2\pi}d\phi\
\bigg{\langle}\gamma
f+\frac{\gamma^{2}}{2}(g-f^{2})\bigg{\rangle},\end{split}$ (24)

where $F_{0}$ is the free energy of the un-deformed state which is an
irrelevant constant we may drop. The distribution function for $\phi$ is taken
to be uniform and the average of the constraints is defined as,

$\big{\langle}\cdots\big{\rangle}=\int_{0}^{\pi}d\theta
P(\theta)\int_{0}^{\pi}d\eta P(\eta)\big{(}\cdots\big{)}.$ (25)

It is straightforward to show that the integrals over $\phi$ of $f$ and $g$
vanish, leaving,

$\begin{split}F&=\frac{k_{B}T}{8}\rho\gamma^{2}(\nu
L^{2})^{2}\int_{0}^{\pi}d\theta d\eta\ e^{-\nu
L^{2}(\theta+\eta)/2}\frac{\sin^{2}(\eta+\theta)}{(\eta+\theta)^{2}}\\\
&\times\frac{1}{2\pi}\int_{0}^{2\pi}d\phi\
\sin^{2}(\eta-\theta-2\phi).\end{split}$ (26)

The $\phi$ average is simply $1/2$. To compute the integrals over $\eta$ and
$\theta$, we make the substitutions $x=\nu L^{2}\eta$ and $y=\nu L^{2}\theta$.
We can then use the fact that $\nu L^{2}\gg 1$ to replace the upper
integration limits by $\infty$.

$F=\frac{k_{B}T}{16}\rho\gamma^{2}\int_{0}^{\infty}dxdy\
e^{-(x+y)/2}\frac{\sin^{2}\big{(}\frac{x}{\nu L^{2}}+\frac{y}{\nu
L^{2}}\big{)}}{\big{(}\frac{x}{\nu L^{2}}+\frac{y}{\nu L^{2}}\big{)}^{2}}.$
(27)

This can be calculated for $\nu L^{2}\gg 1$, see appendix A. The final answer
is

$F\approx\frac{k_{B}T\gamma^{2}}{4}\rho\bigg{(}1-\frac{32}{3(\nu
L^{3})^{2}}\bigg{)}.$ (28)

This directly leads to the linear elastic stress using (4),

$\sigma_{0}(\nu)\sim\frac{k_{B}T}{2}\rho\gamma\bigg{(}1-\frac{32}{3(\nu
L^{2})^{2}}\bigg{)},$ (29)

The important features of this expression are that it is approximately
independent of density and it agrees with a calculation based on a reptation
theory.

In the case of rod shaped particles, it is possible to compute the elastic
stress by an alternative argument. The main details of this textbook
derivation are left to appendix B; here we summarise the main results. The
relevant component of the elastic stress is found to be,

$\sigma_{xy}=\frac{k_{B}T}{2}\rho\gamma\begin{cases}e^{-4D_{r}^{0}t}\
,&\text{for}\ t<\tau_{\text{int}}\\\
e^{4\tau_{\text{int}}(D_{r}(\nu)-D_{r}^{0})}e^{-4D_{r}(\nu)t}\ ,&\text{for}\
t>\tau_{\text{int}}\end{cases}.$ (30)

where $D_{r}^{0}$ and $D_{r}(\nu)$ are the rod’s rotational diffusion
coefficients in the absence and presence of constraints respectively. The
latter is estimated using the reptation argument such that, $D_{r}(\nu)\sim
D_{r}^{0}(\nu L^{2})^{-2}$. As the average angle between constraints is
$\sim(\nu L^{2})^{-1}$, the time taken for the rod to interact with the
constraints is approximated by, $D_{r}^{0}\tau_{\text{int}}\sim(\nu
L^{2})^{-2}$. This leads to the initial elastic stress,

$\sigma_{0}(\nu)\sim\frac{1}{2}k_{B}T\rho\gamma e^{-b/(\nu
L^{2})^{2}}\sim\frac{1}{2}k_{B}T\rho\gamma\bigg{(}1-\frac{b}{(\nu
L^{2})^{2}}\bigg{)},$ (31)

where $b$ is a positive constant. This is precisely the same form as in (29),
found by our geometric argument. The reptation argument does not predict the
value of $b$ exactly and, due to the approximations used, we do not expect the
value of $32/3$ in (29) to be precise either. However, the prefactor of $1/2$
and the scaling with density found by our new geometric method agree precisely
with those found using the well known techniques.

###  III.2 Wire Frame Particles

We now consider the simplest variant to a rod, a bent particle with two legs
joined at an angle $\chi$. In the framework outlined previously we have
$\Delta_{l}=l\chi$ in (11) and $l=0,1$. The free energy is then calculated by
essentially the same procedure as for rods, the only difference being that we
need to deal with different possible combinations of constraints. This is done
by weighting the different contributions by their appropriate probabilities,
such that the free energy is,

$\begin{split}F=&-\frac{k_{B}T}{2\pi}\rho\int_{0}^{2\pi}d\phi\sum_{l,m=0,1}\int_{0}^{\pi}d\eta_{l}P(\eta_{l})\int_{0}^{\pi}d\theta_{m}P(\theta_{m})\\\
&\int_{\eta_{l}}^{\pi}d\eta
P(\eta)\int_{\theta_{m}}^{\pi}P(\theta)\log(\tilde{\eta}_{l}+\tilde{\theta}_{m}).\end{split}$
(32)

Making use of the same substitutions as before, $x_{j}=\nu L^{2}\eta_{j}$ and
$y_{j}=\nu L^{2}\theta_{j}$, the integrals on $\eta$ and $\theta$ can be
evaluated simply, with each giving a factor of $e^{-\eta_{l}/2}$ and
$e^{-\theta_{m}/2}$ respectively. The expansion in (16), along with the
definitions (17), are then used to expand the logarithm to second order in
$\gamma$, from which it is straightforward to show that the only contribution
which will not vanish when integrated over $\phi$ is,

$\begin{split}&\frac{k_{B}T\gamma^{2}}{16\pi}\rho\int_{0}^{2\pi}d\phi\sum_{l,m=0,1}\int_{0}^{\infty}dx_{l}\int_{0}^{\infty}dy_{m}\frac{(\nu
L^{2})^{2}e^{-(x_{l}+y_{m})}}{(x_{l}+y_{m})^{2}}\\\
&\times\bigg{[}\sin^{2}\Big{(}\phi+l\chi-\frac{x_{l}}{\nu
L^{2}}\Big{)}-\sin^{2}\Big{(}\phi+m\chi+\frac{y_{m}}{\nu
L^{2}}\Big{)}\bigg{]}^{2}.\end{split}$ (33)

At this stage we can drop the subscripts on $x$ and $y$, and take the sums and
$\phi$ integral to give,

$\begin{split}F=\frac{k_{B}T\gamma^{2}}{16}\rho\int_{0}^{\infty}&dxdy\frac{(\nu
L^{2})^{2}e^{-(x+y)}}{(x+y)^{2}}\\\
&\times\bigg{[}1-\cos^{2}\chi\cos\bigg{(}\frac{2(x+y)}{\nu
L^{2}}\bigg{)}\bigg{]}.\end{split}$ (34)

At first, this integral appears to diverge because of the inclusion of the
point $x=y=0$. However, this point can be omitted for physical reasons. In any
real system there will be a finite excluded volume. This leads to the integral
being cut off at small values, making it finite. The qualitative behaviour of
the integral can be understood by ‘pre-averaging’ the constraints so that the
angles, $\eta$ and $\theta$, are replaced by their average values,
$\langle\eta\rangle=\langle\theta\rangle=2(\nu L^{2})^{-1}$. This sets $x\sim
y\sim 1$ in (34). Expanding the resulting expression in powers of $\nu L^{2}$
gives the free energy as,

$\begin{split}&F=\frac{k_{B}T\gamma^{2}}{2}\rho\Big{[}a\ (\nu
L^{2})^{2}(1-\cos^{2}\chi)+b\ \cos^{2}\chi\Big{]},\end{split}\ $ (35)

where $a$ and $b$ are positive constants. In appendix C, we treat the excluded
volume cut off more rigerously, and find that the scaling found in (35) has
only logarithmic corrections which depend on the exact cutoff. Therefore,
without losing qualitative accuracy, we take the coefficients $a$ and $b$ to
be unknown fitting parameters.

The stress follows easily as,

$\sigma_{0}(\nu)=\rho\Big{[}a\ (\nu L^{2})^{2}(1-\cos^{2}\chi)+b\
\cos^{2}\chi\Big{]}k_{B}T\gamma.$ (36)

This can show a concentration dependence very different from that of rods,
depending on the size of $\chi$. For this different behaviour to be seen we
require,

$(\nu L^{2})^{2}(1-\cos^{2}\chi)\gtrsim\cos^{2}\chi.$ (37)

After some simple re-arrangement this is,

$\tan\chi\gtrsim\frac{1}{\nu L^{2}},$ (38)

which, since $\nu L^{2}\gg 1$, implies that if the particle is bent through an
angle much larger than $\chi_{c}\sim(\nu L^{2})^{-1}$, the dominant scaling of
the stress will be $\propto(\nu L^{2})^{2}$. This means that even modestly
bent particles ($\chi\sim\mathcal{O}(1)$) will have a significantly different
elastic response than rods.

In particular if we consider ‘L-shaped’ particles, where $\chi=\pi/2$ the
elastic stress is simply,

$\sigma_{0}(\nu)\propto k_{B}T\rho(\nu L^{2})^{2}\gamma.$ (39)

This is also observed for other symmetric, N-legged shapes with the legs all
joined at the same angle, $2\pi/N$. Examples of these particles are Y-shapes
and crosses have been fabricated in practice from double stranded DNA Xing
_et al._ (2018); Biffi  _et al._ (2013). It is interesting to note that these
special, symmetric shapes have no elasticity in dilute suspension for
fundamental symmetry reasonsKing, Doi, and Eiser (2020). Suspensions of these
particle shapes therefore will be most sensitive to concentration changes.
This has important implications for designing these DNA systems as functional
nano-materials, where it is desirable to have highly tuneable, concentration
or connectivity dependent properties.

Equations (39) and (29) show that wire frames behave in a strikingly different
way compared to rods. As $\nu L^{2}\gg 1$, we predict that the elastic stress
for wire frames is significantly larger than that of rods, and very sensitive
to the concentration of constraining points. These results are qualitatively
consistent with the simulations, see Fig(4). Here the simulation data for the
stress per particle is plotted as a function of density for rods, Y-shapes
(N=3 stars) and X-shapes (N=4 stars) and compared to the theoretical
predictions (note that the data for rods has been multiplied by 20 for ease of
visualisation). The stress is calculated in the simulations by fitting the
free energy to a quadratic for strains between $\gamma=0.02$ and $0.1$.

There is, however, a quantitative difference between the predicted scaling and
that of the data. The stress for branched particles is predicted to scale as
$\nu^{2}$, but the fitted exponent for the data is noticeably smaller, $\sim
1.4$. In the simulations particles which bend are removed from the system, i.e
the bending modulus $K=0$. As will be shown later (section IV), this leads to
strong shear thinning in these systems even for relatively small strains,
leading to the apparently smaller stress than predicted by the linear theory.

Figure 4:  A log-log plot of the elastic stress per particle $\sigma_{0}$ as a
function of density, $\nu$, for rods, Y and X shaped particles. The points
represent the simulation data; circles for rods, triangles for Y-shapes and
squares for X-shapes. Note that the data for rods has been multiplied by 20.
The dashed line is the theoretical result for rods (multiplied by twenty) and
the solid line the scaling theory for branched or bent particles, with the
intercept fitted. There is clear qualitative agreement between the theory and
the data. For rods the stress is small and approximately independent of
density, whereas for branched particles the stress is much larger and
increases rapidly with the density. On the other hand there is a quantitative
discrepancy. The predicted scaling is $\propto\nu^{2}$, but the fitted
exponent for the data is smaller $\sim 1.4$. This can be explained by
deviations from the linear elasticity at very small strains in these systems.

##  IV General Method: Non-Linear Response

We now extend our approach to study the non-linear elasticity of these
systems. To calculate this we seek an expression for the free energy accurate
to order $\gamma^{4}$. In addition to expanding the transformed entropy,
$\log\tilde{\Omega}$, to this order, we need to introduce another physical
process; the deformation of the particles themselves. For example, as sketched
in Fig.(5), an L shaped particle is forced to bend when the surrounding
constraining points are displaced by the shear deformation. This situation
cannot occur for straight rods, but can for general wire frame particles.

Figure 5:  A sketch of the bending mechanism for an L-Shaped wire frame. (a)
An undeformed configuration. The clockwise rotation of the L-shaped particle
is constrained by the point, $P_{0}$, which hinders the leg, $OP_{0}$. The
anti-clockwise rotation is constrained by the point $P_{1}$ hindering leg
$OL_{1}$. The L-shape can only rotate between the orange line, $OP_{0}$, and
the green line, $OP_{1}$. The accessible configuration space is the angle
between these two lines, $\Omega$. (b) When the step strain is applied, the
points $P_{0}$ and $P_{1}$ are displaced to $\tilde{P}_{0}$ and
$\tilde{P}_{1}$ respectively. The L-shape must now lie between the lines
$O\tilde{P}_{0}$ and $O\tilde{P}_{1}$. The line $O\tilde{P}_{1}$ is now to the
right of $O\tilde{P}_{0}$, therefore the L-shape cannot keep its original
shape and $\tilde{\Omega}<0$. (c) To satisfy the constraints, the particle
must deform from its original state. Assuming its two legs stay straight, the
angle between them becomes $\pi/2-|\tilde{\Omega}|$.

To demonstrate this, consider two lines at angles $\alpha$ and $\beta$ to the
$x$-axis with $\alpha<\beta$. These represent the constraints placed on a
particular particle; the particle is trapped between these two lines. The
particle will need to bend if, after the deformation, these lines cross. For
this to happen we must be able to find a pair of angles, $\alpha\neq\beta$,
for which $\tilde{\alpha}=\tilde{\beta}$. Initially, take the particle to be a
rod. This means that both lines will transform in the same way, according to
(15), which can be re-arranged to,

$\cot\tilde{\alpha}=\cot\alpha+\gamma.$ (40)

Subtracting this from the equivalent expression for $\tilde{\beta}$, and
insisting that the transformed angles need to be the same we find that
$\cot\alpha=\cot\beta$, which implies, $\alpha=\beta+n\pi$. Due to the
symmetry of the particle shape, rotation of one of the constraining angles by
an integer multiple of $\pi$ results in the same constraint. Therefore, for
rods, $\tilde{\alpha}=\tilde{\beta}$ if and only if $\alpha=\beta$ and hence,
bending is impossible.

For L-shaped particles, if the constraints are placed on different legs, the
transformation rule for one of the angles is different. Taking $\beta$ to be a
constraint on one leg and $\alpha$ to be on the other, then $\beta$ will
transform as if it were rotated by $\pi/2$ with respect to $\alpha$. Hence
$\tilde{\beta}$ is given by,

$\cot\tilde{\beta}=\gamma-\tan\beta.$ (41)

Again, we subtract this from (40) and insist $\tilde{\alpha}=\tilde{\beta}$
which yields, $\cot\alpha=\tan\beta$. This has solutions,
$\beta=\alpha-\pi/2+n\pi$. Clearly it is possible to find valid solutions with
$\alpha\neq\beta$ and so for these particles, bending is possible for certain
configurations.

To deal with the possibility of bending we introduce the function,
$P(\gamma;\Gamma)$. This is the proportion of particles with orientation
$\Gamma$ not bent at strain $\gamma$. The general form of this function can be
determined from simple considerations. When the system is not strained, none
of the particles are bent, hence $P(0;\Gamma)=1$. The behaviour of the system
is symmetric under a change in sign of $\gamma$, which implies,
$P(\gamma;\Gamma)=P(-\gamma;\Gamma)$. The size of strain is taken to be small
so that this function is expanded to give,

$P(\gamma;\Gamma)=1-\gamma^{2}/\gamma_{c}^{2}(\Gamma),$ (42)

where we have defined the critical strain $\gamma_{c}(\Gamma)$, whose scaling
with density will be determined in the following section. This critical strain
should be interpreted as the strain above which all particles in orientation
$\Gamma$ are bent. This strain also indicates roughly the limit of accuracy of
the treatment presented here, as above this strain additional processes not
taken into account will need to be considered, such as the non linear
elasticity or potential breakage of the particles themselves.

We may also consider the orientational average of (42),

$P(\gamma)=1-\gamma^{2}/\gamma_{c}^{2}.$ (43)

This should be interpreted as the total fraction of particles which have not
bent at strain $\gamma$. We will use this for comparison to simulations
through the scaling of $\gamma_{c}$ with the density.

When a particle bends, it stores some elastic energy. This must contribute to
the free energy of the system and $P(\gamma;\Gamma)$ can be used to determine
this bending contribution. The probability that a particle first bends at a
strain $\gamma$ is $-\partial P/\partial\gamma$. If the applied strain is
$\gamma$, but a particle in orientation $\Gamma$ first bent at
$\gamma^{\prime}<\gamma$, then the particle must bend by an angle
$\vartheta=(\gamma-\gamma^{\prime})\Delta(\Gamma)$, where $\Delta$ depends on
the particle shape. If the bending modulus of the particle is $K$, the energy
associated with this bending is, $K\vartheta^{2}/2$. The bending contribution
to the free energy from one particle in orientation $\Gamma$ is given by;

$F_{B}(\gamma;\Gamma)=\frac{K}{2}\rho\Delta^{2}(\Gamma)\int d\gamma^{\prime}\
(\gamma-\gamma^{\prime})^{2}\bigg{(}-\frac{\partial
P(\gamma;\Gamma)}{\partial\gamma}\bigg{\lvert}_{\gamma^{\prime}}\bigg{)},$
(44)

where all the possible strains at which the particle could first bend have
been summed over, weighted by the appropriate probability. We can now write an
expression for the free energy per unit volume which includes both the
entropic contribution and the bending contribution.

$\begin{split}&F(\gamma)=-k_{B}T\rho\int d\Gamma\
\psi(\Gamma)P(\gamma;\Gamma)\big{\langle}\log\tilde{\Omega}(\Gamma;\gamma)\big{\rangle}\\\
&+\frac{K}{2}\rho\int d\Gamma\
\psi(\Gamma)\Delta^{2}(\Gamma)\int_{0}^{\gamma}d\gamma^{\prime}\
(\gamma-\gamma^{\prime})^{2}\bigg{(}-\frac{\partial
P}{\partial\gamma}\bigg{\lvert}_{\gamma^{\prime}}\bigg{)}.\end{split}$ (45)

The orientation needs to be averaged over the distribution function
$\psi(\Gamma)$, which for our case can be taken to be uniform. The first term
in the above equation represents the entropic contribution, $F_{S}(\gamma)$.
The factor of $P$ needs to be included since a particle which bends is forced
to be in one configuration by the constraints, and as such does not contribute
entropically.

The form of the free energy (45) is the basis of our treatment of the non-
linear elasticity in these systems. We will use this to determine the free
energy to order $\gamma^{4}$ first for rods, then for L-shaped particles. The
latter shape is taken for ease, but we would find the same features for any
bent or branched particle.

###  IV.1 Rods

The extension to non-linear elasticity for rods is straightforward. As
discussed previously, there is no possibility for bending, so
$P(\gamma;\Gamma)=1$, and all that needs to be done is to expand
$\tilde{\Omega}$ and the entropy to order $\gamma^{4}$. After lengthy but
straightforward algebra, the first and third order contributions to the free
energy can be shown to vanish after integrating over $\phi$. The second order
contribution is that found previously in (28). The fourth order contribution
is, after averaging on $\phi$,

$\bigg{\langle}\frac{\gamma^{4}\sin^{3}(\eta+\theta)}{32(\eta+\theta)^{4}}\Big{[}4(\eta+\theta)\cos(\eta+\theta)+\big{(}2(\eta+\theta)^{2}-3\big{)}\Big{]}\bigg{\rangle}.$
(46)

The averages over $\eta$ and $\theta$ are then taken in the same way as
before, and the free energy to 4th order in $\gamma$ is obtained,

$\frac{F(\gamma)}{\rho k_{B}T}=\frac{\gamma^{2}}{4}\bigg{(}1-\frac{c}{(\nu
L^{2})^{2}}\bigg{)}-\frac{\gamma^{4}}{32}\bigg{(}1-\frac{d}{(\nu
L^{2})^{4}}\bigg{)},$ (47)

where $c$ and $d$ are positive constants. The elastic stress is, dropping the
subdominant terms for $\nu L^{2}\gg 1$,

$\sigma_{0}=k_{B}T\rho\bigg{(}\frac{1}{2}\gamma-\frac{1}{8}\gamma^{3}\bigg{)}.$
(48)

This shows that the suspension of rods shear thins, as the co-efficient of the
cubic term is negative. This is consistent with the results of a more detailed
theory for this system based on the reptation model Doi and Edwards (1986),
and serves as a useful reference for comparison to the behaviour of the
L-shapes.

###  IV.2 L-Shapes

For L-shapes, bending is possible since configurations can be found where
$\tilde{\Omega}\leq 0$. To compute $P(\gamma;\phi)$ for these shapes, we first
find the probability that a particle with orientation $\phi$ has bent,

$P_{B}(\gamma;\phi)=1-P(\gamma;\phi).$ (49)

This is the probability that $\tilde{\Omega}\leq 0$ which, from the previous
discussion, happens only when the constraints are placed on different legs.
Therefore there are two possibilities for bending to occur,

$\begin{split}\text{(I)}\ \ &\tilde{\eta}_{1}<\tilde{\eta}_{0}\ \ ,\ \
\tilde{\theta}_{0}<\tilde{\theta}_{1}\ \ \text{and}\ \
\tilde{\eta}_{1}+\tilde{\theta}_{0}\leq 0,\\\ \text{(II)}\ \
&\tilde{\eta}_{0}<\tilde{\eta}_{1}\ \ ,\ \
\tilde{\theta}_{1}<\tilde{\theta}_{0}\ \ \text{and}\ \
\tilde{\eta}_{0}+\tilde{\theta}_{1}\leq 0.\end{split}$ (50)

The probability of bending is then the sum of the integrals of the constraint
probabilities over each of these regions. We give the calculation for region
(I) in detail, the steps are the same for (II) so we will simply state the
result. The integral over region (I) is formally,

$\begin{split}P_{B}(\gamma;\phi)=&\int_{\tilde{\eta}_{1}+\tilde{\theta}_{0}\leq
0}d\eta_{1}d\theta_{0}\ P(\eta_{1})P(\theta_{0})\\\
&\times\int_{\tilde{\eta}_{0}>\tilde{\eta}_{1}}d\eta_{0}\
P(\eta_{0})\int_{\tilde{\theta}_{1}>\tilde{\theta}_{0}}d\theta_{1}\
P(\theta_{1}).\end{split}$ (51)

To make progress we require expressions for the transformed constraining
angles, which are deduced straightforwardly from previous definitions,

$\begin{split}&\tilde{\eta}_{0}=\eta_{0}+\gamma\sin^{2}(\phi-\eta_{0})-\gamma^{2}\sin^{3}(\phi-\eta_{0})\cos(\phi-\eta_{0}),\\\
&\tilde{\eta}_{1}=\eta_{1}+\gamma\cos^{2}(\phi-\eta_{1})+\gamma^{2}\cos^{3}(\phi-\eta_{1})\sin(\phi-\eta_{1}),\\\
&\tilde{\theta}_{0}=\theta_{0}-\gamma\sin^{2}(\phi+\theta_{0})+\gamma^{2}\sin^{3}(\phi+\theta_{0})\cos(\phi+\theta_{0}),\\\
&\tilde{\theta}_{1}=\theta_{1}-\gamma\cos^{2}(\phi+\theta_{1})-\gamma^{2}\cos^{3}(\phi+\theta_{1})\sin(\phi+\theta_{1}).\end{split}$
(52)

Let us define the function $A(\theta_{0})$ as the value of $\theta_{1}$ where
$\tilde{\theta}_{1}=\tilde{\theta}_{0}$. The function $\tilde{\theta}_{1}$ in
(52) is monotonic in $\theta_{1}$ for small values of $\gamma$, so it follows
that if $\theta_{1}>A(\theta_{0})$, then
$\tilde{\theta}_{1}>\tilde{\theta}_{0}$. We also define $B(\eta_{1})$ as the
value of $\eta_{0}$ where $\tilde{\eta}_{1}=\tilde{\eta}_{0}$. These functions
are found to be,

$\begin{split}A(\theta_{0})&=\theta_{0}+\gamma\cos(2\theta_{0}+2\phi)\\\
&+\frac{1}{2}\gamma^{2}\big{(}\sin(2\theta_{0}+2\phi)-\sin(4\theta_{0}+4\phi)\big{)},\end{split}$
(53a) $\begin{split}B(\eta_{1})&=\eta_{1}+\gamma\cos(2\eta_{1}-2\phi)\\\
&-\frac{1}{2}\gamma^{2}\big{(}\sin(2\eta_{1}-2\phi)+\sin(4\eta_{1}-4\phi)\big{)}.\end{split}$
(53b)

The integrals over $\theta_{1}$ and $\eta_{0}$ in (51) are now,

$\int_{B(\eta_{1})}^{\pi}d\eta_{0}\
P(\eta_{0})\int_{A(\theta_{0})}^{\pi}d\theta_{1}\ P(\theta_{1}).$ (54)

However it can be shown that $A(0)<0$ and $B(0)<0$ when,

$\frac{\pi+\gamma}{4}<\phi<\frac{3\pi+\gamma}{4}\ \ \text{and}\ \
\frac{5\pi+\gamma}{4}<\phi<\frac{7\pi+\gamma}{4}.$ (55)

In which case, the lower limits on the integrals in (54) are both replaced by
zero.

In a similar fashion, we define the function $\eta_{1}(\theta_{0})$ where
$\tilde{\eta}_{1}+\tilde{\theta}_{0}=0$,

$\begin{split}\eta_{1}(\theta_{0})&=-\theta_{0}-\gamma\cos(2\theta_{0}+2\phi)\\\
&-\frac{1}{2}\gamma^{2}\big{(}\sin(2\theta_{0}+2\phi)-\sin(4\theta_{0}+4\phi)\big{)}.\end{split}$
(56)

This defines a line in the $\theta_{0}$, $\eta_{1}$, plane, which is
approximately a straight line between the points $(0,\eta_{1}(0))$ and
$(\theta_{0}(\phi),0)$, with gradient negative one, as shown in Fig.(6). The
shaded region in Fig.(6) enclosed by the axes and this line is the region
where $\tilde{\Omega}\leq 0$. For certain orientations this region vanishes,
and bending is not possible. For the region to exist, it is necessary for
$\eta_{1}(0)>0$, which from (56) occurs for orientations satisfying (55).
Therefore, the contribution to $P_{B}$ from region (I) is determined from the
integral,

$\int_{0}^{\theta_{0}(\phi)}d\theta_{0}\int_{0}^{\eta_{1}(\theta_{0})}d\eta_{1}\int_{0}^{\pi}d\eta_{0}d\theta_{1}\
P(\eta_{0})P(\eta_{1})P(\theta_{0})P(\theta_{1}),$ (57)

with $\phi$ satisfying (55) and,

$\theta_{0}(\phi)=-\gamma\cos 2\phi-\frac{\gamma^{2}}{2}\big{(}\sin 2\phi+\sin
4\phi\big{)}.$ (58)

Figure 6:  A sketch of the $\theta_{0},\eta_{1}$ plane. The line
$\eta_{1}(\theta_{0})$ with $\gamma=0.1$ and $\phi=\pi/2$, is shown. From (56)
this must be approximately a line with gradient $-1$, which intersects the two
axes at $(0,\theta_{0}(\phi))$ and $(\eta_{1}(0),0)$. The small variations
about this line are too small to be noticeable on this figure. The shaded
region is where $\tilde{\Omega}<0$.

Using the substitutions $x_{l}=\nu L^{2}\eta_{l}$ and $y_{l}=\nu
L^{2}\theta_{l}$ as before, the $\eta_{0}$ and $\theta_{1}$ integrals will
each evaluate to one and we are left with,

$\frac{1}{4}\int_{0}^{\nu L^{2}\theta_{0}(\phi)}dy_{0}\int_{0}^{\nu
L^{2}\eta_{1}(y_{0}/\nu L^{2})}dx_{1}e^{-(x_{1}+y_{0})/2}.$ (59)

This integral is evaluated in appendix D. The result to lowest order in
$\gamma$, is,

$1-e^{-\nu L^{2}\theta_{0}(\phi)/2}-\frac{1}{2}\nu
L^{2}\theta_{0}(\phi)e^{-\nu L^{2}\eta_{0}(0)/2}.$ (60)

Using the definitions of $\theta_{0}(\phi)$ and $\eta_{0}(0)$, this is
expanded to second order in $\gamma$ to give the contribution to $P_{B}$ from
region (I),

$\text{(I)}=\ \ \begin{cases}0,&\text{if}\
\phi\notin\text{(\ref{allowedphi})}\\\ \Big{(}\frac{\nu
L^{2}}{2}\Big{)}^{2}\gamma^{2}\cos^{2}2\phi\ ,&\text{if}\
\phi\in\text{(\ref{allowedphi})}.\end{cases}$ (61)

The contribution from region (II) is calculated by following the same
procedure, and is found to be,

$\text{(II)}=\ \ \begin{cases}0,&\text{if}\
\phi\in\text{(\ref{allowedphi})}\\\ \Big{(}\frac{\nu
L^{2}}{2}\Big{)}^{2}\gamma^{2}\cos^{2}2\phi\ ,&\text{if}\
\phi\notin\text{(\ref{allowedphi})}.\end{cases}$ (62)

The bending probability is simply the sum of these to contributions,

$P_{B}(\gamma;\phi)=\Big{(}\frac{\nu
L^{2}}{2}\Big{)}^{2}\gamma^{2}\cos^{2}2\phi$ (63)

Taking the orientational average of this, we obtain the proportion of un-bent
particles,

$P(\gamma)=1-(\nu L^{2})^{2}\gamma^{2}/8.$ (64)

From which the scaling of $\gamma_{c}$ can be extracted,

$\gamma_{c}\sim\frac{1}{\nu L^{2}}.$ (65)

Even though this was derived explicitly for L-shapes, this scaling is expected
for any bent or branched particles. In Fig.(7) this scaling is compared to
simulation data for Y and X-shaped particles and we see excellent agreement.

Figure 7:  A log-log plot of the critical strain, $\gamma_{c}$, against
dimensionless density, $\pi\nu L^{2}$. The points are data taken from the 2D
simulations, triangles for Y-shaped particles and squares for X-shapes. The
solid line is the scaling predicted theoretically (with the unknown intercept
fitted to the Y-shape data). We see excellent agreement between the predicted
scaling and the data.

The entropic term will need to be expanded to fourth order as it was for rods
but with $P(\gamma,\phi)$ included before the orientational average is taken
so that the integrals are only taken over the region with $\tilde{\Omega}>0$.
The exact calculation is possible albeit long winded. Fortunately the scaling
can be found using the same pre-averaging method as in section III.2 (up to
logarithmic corrections, see appendix C). We replace the angles $\eta$ and
$\theta$ by their average values, and pre-average $\log\tilde{\Omega}$ and
$P(\gamma,\phi)$ over $\phi$ in (45). After several lines of algebra, we find,

$\begin{split}&\langle\log\tilde{\Omega}\rangle=-\frac{(\nu
L^{2})^{2}}{(16)^{2}}\gamma^{2}-\frac{(\nu
L^{2})^{4}}{(16)^{4}}\gamma^{4}.\end{split}$ (66)

The pre-averaged entropic free energy is then found by multiplying this by
(64),

$\begin{split}&\frac{\bar{F}_{S}(\gamma)}{\rho k_{B}T}=\frac{(\nu
L^{2})^{2}}{(16)^{2}}\gamma^{2}+\frac{(\nu
L^{2})^{4}}{(16)^{4}}\gamma^{4}-\frac{(\nu
L^{2})^{4}}{8(16)^{2}}\gamma^{4}.\end{split}$ (67)

Hence, the general form of the free energy to fourth order for bent and
branched particles is,

$\begin{split}&\frac{F_{S}(\gamma)}{\rho k_{B}T}=\frac{a}{2}(\nu
L^{2})^{2}\gamma^{2}-\frac{b}{4}(\nu L^{2})^{4}\gamma^{4},\end{split}$ (68)

where $a$ and $b$ are positive constants.

We now compute the bending contribution to the free energy, which requires us
to determine the angle through which the particle has bent at a particular
strain. At the strain when the particle starts to deform, $\gamma^{\prime}$,
the constraining points are positioned on different legs at $\eta=\theta=0$.
Therefore from (11), $\alpha=\phi$ and $\beta=\phi+\pi/2$, so that the angle
between the two legs is $\chi=\beta-\alpha=\pi/2$. Straining the system
further to $\gamma>\gamma^{\prime}$ changes this angle to
$\tilde{\chi}=\tilde{\beta}-\tilde{\alpha}$. This process bends the particle
through an angle, $\vartheta=\tilde{\chi}-\pi/2$. Using previous results, it
is straightforward to show that,

$\vartheta=-(\gamma-\gamma^{\prime})\cos
2\phi-\frac{1}{2}(\gamma-\gamma^{\prime})^{2}\sin 2\phi.$ (69)

Only the first term is needed for the bending free energy to fourth order,
which is given by

$\begin{split}F_{B}(\gamma)&=\frac{1}{4}K(\nu
L^{2})^{2}\rho\int_{0}^{2\pi}\frac{d\phi}{2\pi}\cos^{4}2\phi\int_{0}^{\gamma}d\gamma^{\prime}(\gamma-\gamma^{\prime})^{2}\gamma^{\prime}\\\
&=\frac{K(\nu L^{2})^{2}}{32}\rho\gamma^{4}.\end{split}$ (70)

Finally, we obtain an expression for the free energy for wire frames accurate
to order $\gamma^{4}$, which includes both the entropic and bending
contributions,

$\frac{F(\gamma)}{\rho k_{B}T}=\frac{a}{2}(\nu
L^{2})^{2}\gamma^{2}+\frac{b}{4}(\nu L^{2})^{2}\bigg{(}\frac{K}{k_{B}T}-c(\nu
L^{2})^{2}\bigg{)}\gamma^{4}.$ (71)

The constants $a,b$ and $c$ are all positive and depend on the particle
geometry chosen. The stress immediately follows,

$\frac{\sigma_{0}}{\rho k_{B}T}=a(\nu L^{2})^{2}\gamma+b(\nu
L^{2})^{2}\bigg{(}\frac{K}{k_{B}T}-c(\nu L^{2})^{2}\bigg{)}\gamma^{3}.$ (72)

This should be compared to the stress for rods (48). Not only is the stress
for the bent particles significantly larger and more sensitive to
concentration than for rods, the sign of the cubic term is not necessarily
negative. This means that while the rod system will always be shear thinning,
the bent particle system can shear harden. If the elastic modulus of the
particles themselves is large enough,

$K\gtrsim(\nu L^{2})^{2}k_{B}T,$ (73)

the cubic term is positive and the response is shear hardening. Alternatively
this condition can be seen as a critical density above which the behaviour
transitions from shear hardening to shear softening,

$\nu_{c}\sim\sqrt{\frac{K}{k_{B}TL^{4}}}.$ (74)

These conditions are non-trivial, because, while we have assumed $K\gg
k_{B}T$, so the particles are very rigid, we also take $\nu L^{2}\gg 1$. If
$K/k_{B}T\sim(\nu L^{2})^{p}$ for any $p\geq 1$, our results are valid but
only for $p>2$ is shear hardening predicted.

This behaviour has an explanation at the level of the model presented here.
When a particle starts to deform, its orientation is completely determined by
the surroundings, and as such cannot contribute to the entropic free energy.
This effect is captured by the $P(\gamma)$ factor in the first term of
equation (45). As the applied strain is increased, more and more particles
begin to bend, so fewer and fewer contribute entropically. This deficit leads
to the shear thinning behaviour of (68). If, at a given strain, the bending
contribution is not sufficient to make up this deficit, the total stress will
be shear thinning. Therefore, there is some critical value of the bending
modulus which must be exceeded to see a shear hardening response.

Another difference between the behaviour of the wire frames and rods is the
relative size of the non-linear term in the stress. For rods both the co-
efficients of $\gamma$ and $\gamma^{3}$ in (48) are of order one, implying
non-linear effects are only important for larger strains. On the other hand,
the coefficient of $\gamma^{3}$ in (72) is $\mathcal{O}((\nu L^{2})^{4})$ but
the coefficient of $\gamma$ is only $\mathcal{O}((\nu L^{2})^{2})$. This shows
that the non-linear effects can begin to be important at very small strains.
This can explain the difference between the simulation data and the linear
theory in Fig.(4). In the simulations $K=0$, so shear softening is expected,
reducing the measured stress compared to what is predicted from the linear
theory.

##  V Discussion

We have discussed the elastic stress response to step strain in a dense
suspension of rod-like, bent and branched particles using a simple geometrical
method. The principle of this method is to determine the volume of
configuration space accessible to a particular particle, given the constraints
placed on it by its surroundings, and how this volume changes when the
constraints are moved by the applied strain. The accessible volume of
configuration space is related to the entropy of the particle, and so when it
changes this leads to a change in the free energy of the system. The stress
calculated from this free energy should be interpreted as that measured
immediately after the step strain is applied.

In this paper we illustrate this method on a simple 2D model system. Here a
test particle is surrounded by point-like constraints which it cannot cross.
These points play the role of the other particles in three dimensions by
providing kinetic constraints. The test particle is free to rotate about its
fixed centre until it collides with the constraints. The constraints are taken
to transform affinely with the applied strain. We focus specifically on the
difference between straight rods and bent or branched wire frames, using
L-shaped particles as an example.

The crucial difference between rods and the other particles is that the
constraints on rotation can be placed on different legs. The consequence of
this is that the range of angles accessible to an L-shape after the strain is
applied can be zero or even negative. This means the particle must have
deformed to satisfy the constraints. This must contribute to the free energy.
This is taken into account by introducing the function, $P(\gamma)$,
interpreted as the proportion of particles which have not bent at a strain
$\gamma$. We find that, to lowest order in $\gamma$, this is given by
$P=1-(\gamma/\gamma_{c})^{2}$, where $\gamma_{c}$ is the critical strain above
which most particles have bent. We determined that $\gamma_{c}\sim(\nu
L^{2})^{-1}$, which agrees well with the scaling found from simple simulations
of the 2D model. With the inclusion of the bending mechanism, the free energy
was found to $\mathcal{O}(\gamma^{4})$ and used to discuss the magnitude of
the linear elastic stress as well as the nature of the non-linear response for
different particle shapes. The linear elastic stress for rods calculated by
this method agrees precisely with that found from a calculation based on a
reptation model.

This simple model shows two interesting differences between rods and wire
frame particles:

1) The linear elastic stress per particle for the wire frames scales
proportional to $(\nu L^{2})^{2}$ whereas for rods it is roughly independent
of density. This is a significant difference. As we are taking $\nu L^{2}\gg
1$, the stress for wire frames is much larger than for rods and much more
sensitive to concentration. This stronger scaling is found for any particle
bent through an angle $\sim\mathcal{O}(1)$ and so the rheology of the system
is highly sensitive to the particle shape. This is in good agreement with
simulation data.

2) There is a critical density beyond which the elastic stress for branched
particles changes from shear hardening to shear softening. The value of this
critical density depends on the elastic modulus of the particle,
$\nu_{c}L^{2}\sim(K/k_{B}T)^{1/2}$. This is also very different from the
behaviour for rods, where the system is always shear thinning.

We hope the simple model presented here may be a useful starting point for
investigating other interesting features of wire frame systems, such the decay
timescale for the stress in the system of branched particles. For such
particles, the reptation process which allows rods to diffuse in a densely
constrained system is no longer possible. Hence we may expect to find very
different time dependence of the stress and the 2D model provides a framework
for investigating this where exact calculations are possible.

###### Acknowledgements.

We gratefully acknowledge Prof. Daan Frenkel for many important and insightful
discussions. D.A.K. acknowledges financial support from the UK Engineering and
Physical Sciences Research Council Ph.D. Studentship award No. 1948692.

##  Appendix A Evaluation of rod free energy integral

The free energy for rods to second order in $\gamma$ is given by the integral,

$F=\frac{k_{B}T}{16}\rho\gamma^{2}\int_{0}^{\infty}dxdy\
e^{-(x+y)/2}\frac{\sin^{2}\big{(}\frac{x}{\nu L^{2}}+\frac{y}{\nu
L^{2}}\big{)}}{\big{(}\frac{x}{\nu L^{2}}+\frac{y}{\nu L^{2}}\big{)}^{2}}.$
(75)

To compute this in the limit $\nu L^{2}\to\infty$, we expand the integrand in
the variable, $y/(\nu L^{2})$, which yields

$\frac{k_{B}T}{16}\rho\gamma^{2}\int_{0}^{\infty}dxdy\
e^{-(x+y)/2}\bigg{[}\frac{\sin^{2}(\frac{x}{\nu L^{2}})}{(\frac{x}{\nu
L^{2}})^{2}}+C\Big{(}\frac{x}{\nu L^{2}}\Big{)}\frac{y}{\nu L^{2}}\bigg{]},$
(76)

where $C(a)\equiv 2\sin(a)(a\cos(a)-\sin(a))/a^{3}$. Integrating over $y$
gives,

$\frac{k_{B}T}{8}\rho\gamma^{2}\int_{0}^{\infty}dx\
e^{-x/2}\bigg{[}\frac{\sin^{2}(\frac{x}{\nu L^{3}})}{(\frac{x}{\nu
L^{3}})^{2}}+\frac{2}{\nu L^{2}}C\Big{(}\frac{x}{\nu L^{3}}\Big{)}\bigg{]}.$
(77)

Expanding the integrand in powers of $x/(\nu L^{3})$ we finally obtain the
free energy as given in equation (28) of the main text,

$\begin{split}F&\approx\frac{k_{B}T}{8}\rho\gamma^{2}\int_{0}^{\infty}dx\
e^{-x/2}\bigg{(}1-\frac{x(4+x)}{3(\nu L^{2})^{2}}\bigg{)}\\\
&\approx\frac{k_{B}T\gamma^{2}}{4}\rho\bigg{(}1-\frac{32}{3(\nu
L^{3})^{2}}\bigg{)}.\end{split}$ (78)

##  Appendix B Calculation for rods based on reptation

Here we follow a textbook procedure to find the elastic stress for rods. Under
$x$-$y$ shear flow, a rod in two dimensions will rotate with angular velocity
$\dot{\phi}=-\dot{\gamma}\sin^{2}\phi$, about its centre. If the rotational
diffusion co-efficient is $\mathcal{D}_{r}$, then in the presence of this flow
the orientational distribution function, $\psi$ satisfies the Smoluchowski
equation,

$\frac{\partial\psi}{\partial
t}=\mathcal{D}_{r}\frac{\partial^{2}\psi}{\partial\phi^{2}}+\dot{\gamma}\frac{\partial}{\partial\phi}\bigg{(}\sin^{2}\phi\
\psi\bigg{)}.$ (79)

The elastic stress is found from the change in the free energy according to
the virtual work principle (1). The free energy is written,

$F=k_{B}T\rho\int d\phi\ \psi\log\psi.$ (80)

When the shear flow is applied, the distribution function changes and so does
the free energy. The change in the distribution function is determined from
the Smoluchowski equation (79). For a step shear of magnitude $\gamma$ applied
over a very short time $\delta t$, the shear rate is
$\dot{\gamma}=\gamma/\delta t$ and the change in the distribution function is,

$\delta\psi=\gamma\frac{\partial}{\partial\phi}\bigg{(}\sin^{2}\phi\
\psi\bigg{)},$ (81)

while the change in the free energy is,

$\delta F=k_{B}T\rho\int d\phi\ \delta\psi\ \big{(}\log\psi+1\big{)}.$ (82)

Substituting (81) into (82) and integrating by parts we find,

$\delta F=-k_{B}T\rho\gamma\int d\phi\ \sin^{2}\phi\
\frac{\partial\psi}{\partial\phi},$ (83)

one final integration by parts yields the elastic stress,

$\sigma_{xy}=k_{B}T\rho\gamma\big{\langle}\sin 2\phi\big{\rangle}.$ (84)

We now obtain an equation for the average $\langle\sin 2\phi\rangle$. By
definition, the time derivative of this is,

$\frac{\partial}{\partial t}\big{\langle}\sin 2\phi\big{\rangle}=\int d\phi\
\sin 2\phi\ \frac{\partial\psi}{\partial t}.$ (85)

Using the Smoluchowski equation and integrating by parts it is easy to show
that,

$\frac{\partial}{\partial t}\big{\langle}\sin
2\phi\big{\rangle}=-4\mathcal{D}_{r}\big{\langle}\sin
2\phi\big{\rangle}-2\gamma\big{\langle}\sin^{2}\phi\ \cos 2\phi\big{\rangle}.$
(86)

To calculate the linear elastic stress, the average in the final term can be
replaced by its equilibrium value, $\big{\langle}\sin^{2}\phi\ \cos
2\phi\big{\rangle}=-1/4$. This allows (86) to be solved easily and using (84)
we find,

$\sigma_{xy}=\frac{1}{2}k_{B}T\rho\gamma e^{-4\mathcal{D}_{r}t}.$ (87)

So far, we have not mentioned the constraints placed on the rod which will
clearly have an effect on the dynamics, altering the diffusion constant in the
Smoluchowski equation.

The constraints effectively force the rod to move along its length in a
channel of width $a\sim 1/(\nu L)$ and length approximately $L$. For a short
time, the rod is unaffected by these constraints, since it will not have
diffused far enough to be hindered by them. This is the case while it has not
rotated through and angle of about $a/L\sim 1/(\nu L^{2})$. If the diffusion
constant without constraints is $D_{r}^{0}$ then this is true for times up to
of order $\tau_{\text{int}}\sim 1/(D_{r}^{0}\nu^{2}L^{4})$, as quoted in the
main text. During this time period then, the stress decays according to (87),
but with $\mathcal{D}_{r}=D_{r}^{0}$.

For times longer than $\tau_{\text{int}}$, the rods orientation can only
change via the reptation process. When the rod has diffused roughly $L/2$
along its length, it can disengage from its original channel, and rotate by an
angle of order $a/L$. The time scale for this process is approximately
$\tau\sim L^{2}/D_{t}^{0}$, where $D_{t}^{0}\sim D_{r}^{0}L^{2}$ is the
translational diffusion constant without constraints. Therefore the rotational
diffusion coefficient estimated for the reptation process is,

$D_{r}(\nu)\sim\frac{a^{2}}{L^{2}\tau}\sim\frac{a^{2}D_{t}^{0}}{L^{4}}\sim\frac{D_{r}^{0}}{(\nu
L^{2})^{2}},$ (88)

as given in the main text. Hence for times longer than $\tau_{\text{int}}$,
the stress decays as in (87) but with $\mathcal{D}_{r}=D_{r}(\nu)$.

Putting all of these pieces together and ensuring that the stress is
continuous, we recover equation (30) of the main text.

##  Appendix C A more rigorous treatment of integrals

Here we treat the divergent integral (34) more rigorously with a small lower
cutoff. This cutoff will ultimately be set by the excluded volume of the legs
of the wire frame shapes. If the width of the legs is $b$, then the
probability of finding a constraining point with angles less than or
approximately $b/L$ must be zero due to the excluded volume effect. We define
the small quantity, $\varepsilon=b/L$, and the reduced density $p=\nu L^{2}$
for convenience. Thus the dominant term in (34) is,

$p^{2}\int_{\varepsilon}^{\pi}d\theta\int_{\varepsilon}^{\pi}d\eta\frac{e^{-p(\theta+\eta)}}{(\theta+\eta)^{2}}.$
(89)

The integral is taken over a square with side length $\pi$ with a smaller
square with side length $\varepsilon$ taken out around the origin. To make
progress without changing the scaling we approximate this region as a quarter
circle of radius $\pi$ with a small quarter circle radius $\varepsilon$
removed. The integral is then written in polar co-ordinates, by defining the
vectors $\textbf{r}=(\theta,\eta)$ and $\textbf{d}=(1,1)$ and the angle
between them as $\varphi$,

$p^{2}\int_{\varepsilon}^{\pi}\frac{dr}{r}\int_{0}^{\pi/4}d\varphi\frac{e^{-\sqrt{2}pr\cos\varphi}}{\cos^{2}\varphi}.$
(90)

Transforming the radial integral to the variable, $R=pr$, we have, for large
$p$,

$p^{2}\int_{p\varepsilon}^{\infty}\frac{dR}{R}\int_{0}^{\pi/4}d\varphi\frac{e^{-\sqrt{2}R\cos\varphi}}{\cos^{2}\varphi}.$
(91)

The $R$ integral can then be evaluated as an incomplete gamma
functionAbramowitz and Stegun (1964),

$p^{2}\int_{0}^{\pi/4}\frac{d\varphi}{\cos^{2}\varphi}\
\Gamma(0,\sqrt{2}p\varepsilon\cos\varphi).$ (92)

Since we are taking the density to be in the range, $1/L^{2}\ll\nu\ll
1/V_{exc}$, we may take the limit $\nu\varepsilon\to 0$, giving

$-p^{2}\int_{0}^{\pi/4}\frac{d\varphi}{\cos^{2}\varphi}\
\Big{[}\gamma_{em}+\log(\sqrt{2}p\varepsilon\cos\varphi)\Big{]},$ (93)

where $\gamma_{em}$ is the Euler-Mascheroni constant. The integral can now be
taken exactly using integration by parts for the second term,

$\begin{split}&-p^{2}\Big{[}\tan\varphi\bigg{(}\gamma_{em}+1+\log(\sqrt{2}p\varepsilon\cos\varphi)\bigg{)}-\varphi\Big{]}\Big{\lvert}_{0}^{\pi/4}\\\
&=p^{2}\big{(}|\log p\varepsilon|-C\big{)},\end{split}$ (94)

where we have used $p\varepsilon\ll 1$, to make the sign more obvious and
defined the positive, order unity constant, $C$.

Hence we find only logarithmic corrections to the $\nu^{2}$ scaling discussed
in the main text. This does not change the qualitative features of our
results.

##  Appendix D Evaluation of Bending Probability Integrals

To get the contribution to the bending probability from region (I), see (50),
we need to evaluate the integral,

$\frac{1}{4}\int_{0}^{\nu L^{2}\theta_{0}(\phi)}dy_{0}\int_{0}^{\nu
L^{2}\eta_{1}(y_{0}/\nu L^{2})}dx_{1}e^{-(x_{1}+y_{0})/2}.$ (95)

Where the functions, $\theta_{0}(\phi)$ and $\eta_{1}(y/\nu L^{2})$ are
defined in equations (58) and (56) of the main text respectively. Taking the
$x_{1}$ integral gives,

$\frac{1}{2}\int_{0}^{\nu
L^{2}\theta_{0}(\phi)}dy_{0}\bigg{(}e^{-x_{1}/2}-e^{x_{1}/2-\nu
L^{2}\eta_{1}(y_{0}/\nu L^{2})/2}\bigg{)}.$ (96)

Integrating the first term in the integrand is straightforward whereas the
second integral is schematically,

$\int_{0}^{c\gamma+d\gamma^{2}}dx\ e^{\nu L^{2}\gamma a(x/\nu L^{2})+\nu
L^{2}\gamma^{2}b(x/\nu L^{2})}.$ (97)

The integrand can be expanded in powers of $x/\nu L^{2}$, because over the
whole integration range $x/\nu L^{2}\sim\mathcal{O}(\gamma)$,

$\int_{0}^{c\gamma+d\gamma^{2}}dx\ e^{\nu L^{2}\gamma a(0)+\nu
L^{2}\gamma^{2}b(0)}\Big{[}1+x\big{(}\gamma
a^{\prime}(0)+\gamma^{2}b^{\prime}(0)\big{)}\Big{]},$ (98)

clearly the integral of the second term in the square brackets will be
$\mathcal{O}(\gamma^{3})$, and so only the first term contributes to second
order. Therefore and the result of the integral (96) is,

$1-e^{-\nu L^{2}\theta_{0}(\phi)/2}-\frac{1}{2}\nu
L^{2}\theta_{0}(\phi)e^{-\nu L^{2}\eta_{0}(0)/2},$ (99)

as given in equation (60) of the main text.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Seeman (1982) N. C. Seeman, “Nucleic acid junctions and lattices,” Journal of Theoretical Biology (1982), 10.1016/0022-5193(82)90002-9.
  * Rothemund (2006) P. W. Rothemund, “Folding DNA to create nanoscale shapes and patterns,” Nature 440, 297–302 (2006).
  * Bomboi  _et al._ (2019) F. Bomboi, D. Caprara, J. Fernandez-Castanon, and F. Sciortino, “Cold-swappable DNA gels,” Nanoscale 11, 9691–9697 (2019).
  * Biffi  _et al._ (2013) S. Biffi, R. Cerbino, F. Bomboi, E. M. Paraboschi, R. Asselta, F. Sciortino, and T. Bellini, “Phase behavior and critical activated dynamics of limited-valence DNA nanostars,” Proceedings of the National Academy of Sciences 110, 15633–15637 (2013).
  * Xing  _et al._ (2018) Z. Xing, A. Caciagli, T. Cao, I. Stoev, M. Zupkauskas, T. O’Neill, T. Wenzel, R. Lamboll, D. Liu, and E. Eiser, “Microrheology of DNA hydrogels,” Proceedings of the National Academy of Sciences 115, 8137–8142 (2018).
  * Gross  _et al._ (2011) P. Gross, N. Laurens, L. B. Oddershede, U. Bockelmann, E. J. Peterman, and G. J. Wuite, “Quantifying how DNA stretches, melts and changes twist under tension,” Nature Physics 7, 731–736 (2011).
  * Van Ketel, Das, and Frenkel (2005) W. Van Ketel, C. Das, and D. Frenkel, “Structural arrest in an ideal gas,” Physical Review Letters 94, 8–11 (2005).
  * Doi (1975) M. Doi, “Rotational relaxation time of rigid rod-like macromolecule in concentrated solution,” Journal de Physique 36, 607–611 (1975).
  * Doi and Edwards (1986) M. Doi and S. F. Edwards,  _The theory of polymer dynamics_ (Oxford University Press, 1986).
  * Gennes (1975) P. G. D. Gennes, “Reptation of stars,” Journal de Physique 36, 1199–1203 (1975).
  * Pearson and Helfand (1984) D. S. Pearson and E. Helfand, “Viscoelastic Properties of Star-Shaped Polymers,” Macromolecules 17, 888–895 (1984).
  * Doi and Kuzuu (1980) M. Doi and N. Y. Kuzuu, “Rheology of star polymers in concentrated solutions and melts,” Journal of Polymer Science: Polymer Letters Edition 18, 775–780 (1980).
  * Milner and McLeish (1997) S. T. Milner and T. C. McLeish, “Parameter-free theory for stress relaxation in star polymer melts,” Macromolecules 30, 2159–2166 (1997).
  * Heine, Petersen, and Grest (2010) D. R. Heine, M. K. Petersen, and G. S. Grest, “Effect of particle shape and charge on bulk rheology of nanoparticle suspensions,” Journal of Chemical Physics 132 (2010), 10.1063/1.3419071.
  * Petersen, Lane, and Grest (2010) M. K. Petersen, J. M. D. Lane, and G. S. Grest, “Shear rheology of extended nanoparticles,” Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 82, 1–4 (2010).
  * Doi and Edwards (1978) M. Doi and S. F. Edwards, “Dynamics of rod-like macromolecules in concentrated solution. Part 1,” Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics 74, 560–570 (1978).
  * Edwards and Vilgis (1986) S. F. Edwards and T. Vilgis, “The dynamics of the glass transition,” Physica Scripta 1986, 7–16 (1986).
  * Edwards (1967) S. F. Edwards, “Statistical mechanics with topological constraints: I,” Proceedings of the Physical Society 91, 513–519 (1967).
  * Edwards (1969) S. F. Edwards, “Theory of cross-linked polymerized material,” Journal of Physics C: Solid State Physics 2, 1–13 (1969).
  * Edwards and Freed (1969) S. F. Edwards and K. F. Freed, “The entropy of a confined polymer. I,” Journal of Physics A: General Physics 2, 145–150 (1969).
  * King, Doi, and Eiser (2020) D. A. King, M. Doi, and E. Eiser, “Particle shapes leading to Newtonian dilute suspensions,” Phys. Rev. E 102, 32615 (2020).
  * Abramowitz and Stegun (1964) M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ (Dover Publications, New York City, 1964).