Elastic Response of Wire Frame Glasses. I. Two Dimensional Model
EElastic Response of Wire Frame Glasses. I. 2D Model
Elastic Response of Wire Frame Glasses. I. Two Dimensional Model
David A. King, a) Masao Doi, and Erika Eiser Cavendish Laboratory, University of Cambridge, J. J. Thomson Ave., Cambridge CB3 0HE,UK Centre of Soft Matter and its Applications, Beihang University, Beijing 100191, China (Dated: 5 February 2021)
We study the elastic response of concentrated suspensions of rigid wire frame particles to a step strain. These particlesare constructed from infinitely thin, rigid rods of length L . We specifically compare straight rod-like particles to bentand branched wire frames. In dense suspensions the wire frames are frozen in a disordered state by the topologicalentanglements between their arms. We present a simple, geometric method to find the scaling of the elastic stress withconcentration in these glassy systems. We apply this method to a simple 2D model system where a test particle isplaced on a plane and constrained by a random distribution of points with number density ν . Two striking differencesbetween wire frame and rod suspensions are found: 1) The linear elasticity per particle for wire frames is very large,scaling like ν L , whereas for rods it much smaller and independent of concentration. 2) Rods always shear thin butwire frames shear harden for densities less than ∼ (cid:112) K / k B T L , where K is the bending modulus of the particles. Thedeformation of wire frames is found to be important even for small strains, with the proportion of deformed particles ata particular strain, γ , being given by ( ν L ) γ . Our results agree well with a simple simulation of the 2D system. I. INTRODUCTION
Concentrated suspensions of rigid macro-molecules are aprototypical soft matter system, constituting a rich field ofstudy for both experimentalists and theoreticians. The greatvariety of non-Newtonian flow behaviours displayed by thesesystems has been a particular interest. Being able to predictand engineer these behaviours based on the shape of the sus-pended particles is a problem of broad scientific an technolog-ical importance.The interactions between the particles are generally whatlead to the pronounced non-Newtonian behaviour. Thestrongest effects are observed for concentrated suspensions,where the number density of particles, ρ , is large such that ρ L (cid:29)
1, with L being a typical length-scale of the particles.Hydrodynamic and excluded volume interactions will alwaysbe present and have an effect on the rheology, but for certainparticle shapes and concentrations a different kind of inter-action dominates. These are ‘kinetic’ or ‘topological’ con-straints, which originate from the fundamental property thatno two particles can cross. This is the case when the parti-cles are very long compared to their width and the particlesare rarely in close contact, such that ρ V exc (cid:28)
1, where V exc is the excluded volume of the particles. As the excluded vol-ume is irrelevant, all equilibrium properties are known; thesystem is equivalent to an ideal gas. On the other hand, thenon-equilibrium properties, such as the rheology, are severelyaffected by the non-crossing condition. This is often referredto as the ‘entanglement effect’.We will focus on how kinetic constraints affect the rheol-ogy of the suspension depending on the shape of suspendedparticles. The concentration range considered is, 1 / L (cid:28) ρ (cid:28) / V exc . This ensures that the kinetic constraints are thedominant interaction and that there is no long range orienta-tional order in the system. Each of the particle shapes we a) Electronic mail: [email protected] will consider is constructed from infinitely thin, rigid rods.We consider the joints between the rods to be effectively rigidand unaffected by thermal fluctuations. Such particle shapescan be divided into three different classes; straight, bent andbranched. Straight particles are simply rods, see Figs.(1a) &(2c). Bent particles are two rods joined at some angle, seeFigs.(1b) & (2a). Branched particles are those where three ormore component rods meet at a point. These can be planarstars as in Fig.(2a) or fully three dimensional shapes, for in-stance the shapes in Fig.(2d) with an extra rod piercing theplane. We refer to the bent or branched particles as ‘wireframes’.Wire frame particles can be realised practically using DNAorigami techniques . Arms made of double stranded DNAsequences can be engineered so that they rigidly link at preciseangles to form ‘nanostars’ . Modelling the double strandedDNA arms as infinitely thin, perfectly rigid rods is reasonable,due to their large aspect ratio ( ∼
20) and stiffness (persistencelength ∼ . Being able to predict the different rheologi-cal responses depending on the shape of the nano-stars is veryimportant for the design of these systems as functional nano-materials.The kinetic constraints in these systems can also be re-sponsible for a glass transition. Suspensions of 3D crosses,particles consisting of three mutually perpendicular, infinitelythin rods joined at their midpoints, have been shown to haveglassy dynamical behaviour in the range of densities we areconsidering . In this system the translational diffusion coeffi-cient goes to zero with an exponential dependence on the den-sity. The majority of particles are rendered effectively immo-bile by the kinetic constraints imposed by their surroundingsand the system is frozen in an isotropic, disordered state.The geometry of the suspended particles must influence thisbehaviour. In Fig.(1) we sketch two examples of the systemswe are considering. A dense suspension of rods is shown inFig.(1a) and a dense glassy state of L-shaped wire frames isshown in Fig.(1b). In each panel a particular test particle ishighlighted in blue. The rod may still diffuse in this concen- a r X i v : . [ c ond - m a t . s o f t ] F e b lastic Response of Wire Frame Glasses. I. 2D Model 2 FIG. 1. Sketches of the systems considered. (a) A dense suspensionof rod like particles. A test rod is shown in cyan. The motion ofthe test rod is restricted by the surrounding red rods, but it can stilldiffuse along its length, so the suspension is still in a fluid state. (b)L-shaped particles in dense suspension. The cyan test L-shape istrapped by the surrounding red L-shapes. If it moves along the lengthof one of its legs, the other becomes entangled. The system is frozenin a glassy state. tration regime by the reptation mechanism . However, thewire frames cannot. This is because the reptation process re-lies on the rods’ ability to diffuse freely along their lengtheven though their transverse motion is severely hindered. If awire frame diffuses along the length of one of its constituentrods, then any one of the others will quickly become entan-gled, blocking its motion.It is interesting to compare the rheology of glassy states ofwire frames, Fig.(1b), to dense fluid states of rods, Fig.(1a).For flexible polymers, it is well known that branched starpolymers display significantly different flow behaviour fromlinear chains, with extremely slow stress relaxation . Thisis as a result of the branching making reptation impossible. Itis reasonable then to expect an equally pronounced differencein behaviour between rigid branched particles and rods. Infact, it has been observed in simulations that the viscosityfor dense suspensions of branched 3D crosses is significantlylarger than for rods and that the viscosity has a much strongerdependence on concentration .Evidently, we cannot use the standard approach based onthe tube model , where the reptation mechanism is used ina hydrodynamic formalism. Instead we must use a new ap-proach which determines the rheological properties of theseglassy states directly from the constraints on each particle. Inthis two part series of papers, we restrict our attention to a sim-ple question which gives useful insight into the problem. Wedo not address in detail the time (or frequency) dependence ofthe rheological response here, but rather focus on its instanta-neous magnitude. Our hope is that the method presented herewill form a useful basis for more detailed studies in the future.We begin with some simple definitions before stating theproblem explicitly. The rheological properties of a suspensionare expressed by the constitutive equation relating the stresstensor to the applied rate strain tensor, E ( t ) . The stress willconsist of two parts; the elastic stress, which we call σ , andthe viscous stress, σ V . The viscous stress is related to theenergy dissipation in the system and can generally be writtenas a function of time, σ V ( t ) = V : E ( t ) , where the 4th ranktensor V can be a function of the deformation history of the suspension. The elastic stress is related to the change in thefree energy per unit volume, F , by the virtual work principle; δ F = σ : E δ t , (1)where δ F ≡ F ( E ) − F ( ) is the change in the free energy, cal-culated as the difference between the free energy in a systemdeformed by E and that in the undeformed state. If the ap-plied strain is small enough, the elastic stress will depend onit linearly, σ ( t ) = (cid:90) t − ∞ dt (cid:48) E ( t − t (cid:48) ) : E ( t (cid:48) ) , (2)where the time dependent, fourth rank tensor E is the elasticresponse function.To probe the rheology of these systems, we consider thestress response to a step strain. In this case, the rate of straintensor is given by, E ( t ) = κδ ( t ) , (3)where we have defined the strain tensor, κ . We always takethe strain to be simple x − y shear for which κ has only onenon vanishing component, κ xy ≡ γ . The extension to otherstrains is straightforward. For step strain, the viscous stressmust also be proportional to a delta function, and so can betaken to vanish for all practical timescales.The virtual work principle for the elastic stress can be re-written, σ = ∂ F ∂ κ , (4)and for small strains will take the form, σ ( t ) = E ( t ) : κ = γ E xy ( t ) , (5)hence the choice of step strain allows us to directly probe theelastic response function.Typically the stress will start at a particular value, S , andthen decay. The initial timescale for this decay, τ , is expectedto be the same for all concentrations. This is because imme-diately after the deformation, the particles are free to movewithout hindrance from the surroundings. After a time, τ int ,the particles have moved sufficiently to begin interacting withtheir neighbours and the decay timescale will increase to τ ( ρ ) ,due to the entanglement effect. This new decay rate is a func-tion of the particle density, since for denser systems the entan-glement effect is stronger. Therefore a reasonable approxima-tion for the time dependent elastic stress is, σ ( t ) = (cid:40) S e − t / τ , for t < τ int S e τ int ( / τ ( ρ ) − / τ ) e − t / τ ( ρ ) , for t > τ int . (6)In concentrated suspensions, it is known that τ (cid:28) τ andhence we can ignore the initial rapid decay of the stress andtake σ ( t ) to be approximated by, σ ( t ) ≈ S e − τ int / τ e − t / τ ( ρ ) ≡ σ ( ρ ) e − t / τ ( ρ ) , (7)lastic Response of Wire Frame Glasses. I. 2D Model 3where we have defined σ ( ρ ) as the initial value of the stressmeasured in response to a step strain. Throughout these pa-pers, we refer to σ ( ρ ) as simply ‘the elastic stress’.From this discussion we see that there are two main parts tothe stress response. Its initial magnitude, σ ( ρ ) , and τ ( ρ ) , thedominant timescale of the subsequent decay. This timescaleis generally very long, so the elastic stress persists for a longperiod after the initial strain. There are many interesting ques-tions relating to this timescale in the case when reptation isnot possible. In particular, is there a critical concentration atwhich this timescale diverges? Or is the relaxation a single orstretched exponential, as could be expected? . While we donot address these problems here, we hope that the methods weintroduce may be of use in their resolution.In this paper (paper I), we present a simple geometricmethod for determining the scaling of the elastic stress, σ ( ρ ) .This method is applicable to a wide range of wire-frameshapes and is outlined in the next section. In section III weintroduce a two dimensional model system to which we applythis method. We find a striking difference between straight,rod-like particles and bent or branched particles. In sec-tion IV, we extend our treatment to the non-linear elasticityand show that there is a critical density, related to the bend-ing modulus of the individual particles, where the behaviourchanges from shear hardening to shear softening. This is incontrast to the result for rods, where the suspension alwaysshear thins. Our theoretical results are found to be consistentwith those of a simulation of the 2D system. The purpose ofthis paper is to introduce the method in the 2D model, wherethe calculation can be performed simply and exactly. A fulltreatment of the problem in three dimensions will be given inpaper II. II. GENERAL METHOD: LINEAR RESPONSE
To determine the elastic stress from the virtual work prin-ciple (4), we need the free energy as a function of strain. Thefree energy is, of course, F = U − T S , where U is the internalenergy, S is the entropy, and T is the temperature. We con-sider the system in the absence of an external potential, there-fore the strain only changes the internal energy if the particlesthemselves are deformed. This will be shown to contributeonly to the non-linear elastic response, which is dealt within section IV. In this section, we outline a simple, geometricmethod for calculating σ from the entropy change per parti-cle caused by the strain in the linear regime.In the glassy states, the motion of a test particle is impairedby the other particles in the system, with its centre of massconfined to a ‘cage’ and only certain orientations accessibleat any given moment. Generally, it is only the Brownian ro-tations of the particle which contribute to the stress, not themotion of its centre of mass. This allows us to treat the centreof mass as fixed and focus only on the constraints placed onthe particle’s orientation.Let us define the orientation of the particle as Γ . In twodimensions, Γ is represented by the angle, φ , that the particlemakes to a reference axis and in three dimensions, three Eu- ler angles specify Γ . The surrounding particles constrain thetest particle to a region in the configuration space. The sizeof this region depends on Γ and the configuration of the sur-rounding particles C , and is written Ω ( Γ , C ) . The test particlewill explore all of Ω over a timescale of roughly τ int , whichis generally very short. Therefore the particle essentially in-stantly equilibrates inside Ω . The entropy associated with thisparticle then follows from the Boltzmann definition, S ( Γ , C ) = k B log Ω ( Γ , C ) . (8)When the step strain is applied, the surrounding particles willmove, and so the accessible volume of configuration spacewill instantaneously change from Ω to ˜ Ω ( κ ) . The change inentropy caused by this transformation alters the free energyand hence yields an elastic stress in the system. This methodis similar to that introduced by Edwards to study flexible poly-mers with topological constraints and crosslinks . As thesurrounding particles rearrange due to their Brownian motion,the test particle can explore more of the configuration space,so Ω begins to increase and the stress relaxes. We are inter-ested in the initial value of the elastic stress, σ , so we onlyneed to consider the instantaneous change after the step strain.If there are ρ particles per unit volume, the free energy perunit volume immediately after the step strain is, F ( κ ) = − k B T ρ (cid:90) d Γ ψ ( Γ ) (cid:68) log ˜ Ω ( Γ , C ; κ ) (cid:69) , (9)where ψ ( Γ ) is the orientational distribution function, whichmay be taken to be uniform since the system in a disorderedstate, and (cid:104)· · · (cid:105) denotes averaging over all configurations, C ,of the constraints. Differentiating this with respect to κ willgive the stress immediately after the step strain, σ .The advantage of equation (9) for the free energy is that˜ Ω can be calculated from purely geometric considerations,without detailed knowledge of the dynamics of the particles.In three dimensions the geometry required will become verycomplicated, especially when defining the constraints placedon the test particle in order to calculate ˜ Ω . In two dimensions,however, the geometry is simple enough that the problem maybe solved exactly. III. 2D MODEL AND SIMULATION DETAILS
The calculation of F ( κ ) , becomes quite complex in threedimensions. Therefore, here we first consider a simpler, twodimensional model.In this model we only consider planar wire frames, asshown in Fig.(2) and each particle is fixed parallel to the x - y plane. The wire frame particles are infinitely thin so notwo lie in the same plane and as such, each particle formsits own independent, 2D system. The rotation of each particlein the plane is constrained by many point-like obstacles whichit cannot cross, as shown in Fig.(3a). The number of obstaclesper unit area is, ν , and in the concentrated limit, ν L (cid:29) L , referred to as ‘legs’, which will be indexed by l . Each leg has the same length, L , and the legs all meet at thesame point which is fixed at the origin. The direction of the l th leg is given by the unit vector, e l , running parallel to it andpointing away from the origin. The orientation of the particleis given by the angle φ of the l = x -axis, so that itsunit vector is, e ≡ ( cos φ , sin φ ) .Imagine rotating the particle clockwise through a full turn.As the particle is rotated one particle leg will eventually col-lide with one of the constraining points. Let us define thevector a l as the position vector of the constraint which the l thleg first crosses during this clockwise rotation. In a similarway we define, b l as the position vector of the correspondingpoint for the anti-clockwise rotation. The important feature ofthese vectors is their angle to the x -axis defined in terms of thecomponents of the vectors as,tan α l = a ly a lx and tan β l = b ly b lx . (10)Further, we define the angles indicated in Fig.(2a), η l = φ + ∆ l − α l and θ l = β l − φ − ∆ l , (11)where ∆ l is the angle between the l = l th leg.The particle is constrained by the points which are closest toit, so the range of accessible angles is, Ω = min l , m η l + θ m . (12)To obtain the free energy as given in equation (9), an averageneeds to be taken over the distribution of constraints, whichmeans averaging over all the angles, η l and θ m . The proba-bility distribution, P ( η l ) and P ( θ m ) , can be determined froma simple argument. For any given area of size a , the averagenumber of points inside is the same, ν a , and independent ofthe number found in a different area. Therefore, the pointsare distributed according to a Poisson distribution. When theparticle rotates through an angle η l , say, the l th leg sweepsout an area of size L η l /
2. The probability that there are noconstraining points in this area is, P ( η l ) = exp ( − ν L η l / ) .By definition the l th leg firsts collide with a constraint af-ter rotating clockwise through an angle η l , so it follows that P ( η l ) = − P (cid:48) ( η l ) . So we find, P ( η l ) = ν L e − ν L η l / , (13)and exactly the same form for P ( θ m ) . FIG. 2. Sketches of the particle shapes considered. (a & b) Anexample of a bent particle, an L-shape with arms of length L . (c) Arod of length 2 L . (d) An example of a branched particle, a ‘5-Star’. Ineach panel the unit direction vectors of the legs, e l , are shown and theorientation of the l = φ from the x -axis, shownas a dashed line. In panels (a) and (b), the L-shape is constrained tolie between the two dotted lines. In (a) the angles, η l and θ l , fromeach leg to the nearest anti- and clockwise constraint respectively,are shown. In (b) the angles which the constraints make to the x -axisare shown; α l for the closest anti-clockwise constraint to leg l and β l for the closest clockwise constraint. The points representing theconstraints (see Fig.(3)) are omitted from this sketch for clarity. Theangle from the l = l th leg is ∆ l . Shown in (a), this is π / l π / When the shear strain is applied, the constraining pointsall move and the angles α l , β l , η l and θ l all change. We as-sume the the constraining points are displaced affinely whenthe shear is applied, therefore the position vectors, a l and b l ,transform according to the rule, a l → ˜ a l = ( I + κ ) · a l , (14)with the same expression for ˜ b l . For the simple x - y shear weconsider, the new angles can be determined from,tan ˜ α l = a ly a lx + γ a ly = tan α l + γ tan α l , (15)and an analogous expression for ˜ β l . We focus on small strainslastic Response of Wire Frame Glasses. I. 2D Model 5 FIG. 3. A sketch of the 2D model system for rods. The rod isshown as a solid line at an angle φ from the dashed x -axis. Panel (a)shows the constraining points surrounding the rod with the two pointsclosest to the rod highlighted in red. These two points constrain theorientation to the range of angles, Ω . In (b) the step strain has beenapplied, indicated by the blue arrows, and the constraining pointsclosest to the rod have moved, with their new positions shown in thedarker, solid colour and their previous positions in a lighter colour.This changes the range of accessible angles to ˜ Ω . with γ (cid:28)
1. Expanding (15) in this limit, we find,˜ α l = α l − γ sin α l + γ sin α l cos α l + O ( γ ) , (16)and similarly for β l . The terms to second order in γ will con-tribute to the linear elasticity. For the non-linear response,more terms must be retained. The angles, ˜ η l and ˜ θ l are thenfound from,˜ η l = φ + ∆ l − ˜ α l and ˜ θ l = ˜ β l − φ − ∆ l , (17)which leads to the transformed range of accessible angles,˜ Ω = min l , m ˜ η l + ˜ θ m . (18)This model forms the basis of a simple numerical simulationwe use for comparison to our analytic results. This computesthe accessible angle and the free energy numerically by thefollowing simple procedure. First, a particle is placed in arandom orientation in the plane. Then the obstacles are placedat random on a circle surrounding the particle; the number ofobstacles is proportional to the density, πν L . The two clos-est obstacles to the particle are then found and their positionsdetermine the accessible angle. The positions of the obsta-cles on the circle are then all moved according to (15). Thelogarithm of the ratio of the accessible angle before and afterthe transformation gives the free energy for that configuration.These steps are repeated for many different configurations ofthe particle and obstacles and the average free energy overthese realisations is calculated. This essentially calculates (9)numerically using a Monte-Carlo integration for the averages. A. Rods
The first shape we consider is a rod of length 2 L , which wecan think of as two legs of length L which are anti-parallel.The unit vectors giving the orientations of these legs are there-fore, e = ( cos φ , sin φ ) = − e . Due to the symmetry of the rod and the shear flow considered, it does not matter whichleg of the rod is constrained. This means we only need toconsider the angles η and θ , from which we will drop thesubscripts in this section. Using equations (16) and (17), wecan determine the transformed angles ˜ η and ˜ θ ,˜ η = η + γ sin ( φ − η ) − γ sin ( φ − η ) cos ( φ − η ) , (19)˜ θ = θ − γ sin ( φ + θ ) + γ sin ( φ + θ ) cos ( φ + θ ) . (20)The transformed range of accessible angles is then calculatedfrom ˜ Ω = ˜ η + ˜ θ ,˜ Ω = Ω + γ sin ( η + θ ) sin ( η − θ − φ ) − γ (cid:20) sin (cid:0) ( η − φ ) (cid:1) + sin (cid:0) ( θ + φ ) (cid:1) − ( η − θ − φ ) sin ( η + θ ) (cid:21) . (21)For convenience let us define the functions, f ( η , θ ; φ ) and g ( η , θ ; φ ) such that,˜ Ω = Ω ( + γ f + γ g / ) . (22)To second order in γ , the entropy of the test particle for thisrealisation is, S = S + k B γ f + k B γ ( g − f ) , (23)where S = k B log Ω is the entropy of the undeformed state.Now the free energy to the same order is according to (9), F = F − k B T π ρ (cid:90) π d φ (cid:28) γ f + γ ( g − f ) (cid:29) , (24)where F is the free energy of the un-deformed state which isan irrelevant constant we may drop. The distribution functionfor φ is taken to be uniform and the average of the constraintsis defined as, (cid:10) · · · (cid:11) = (cid:90) π d θ P ( θ ) (cid:90) π d η P ( η ) (cid:0) · · · (cid:1) . (25)It is straightforward to show that the integrals over φ of f and g vanish, leaving, F = k B T ργ ( ν L ) (cid:90) π d θ d η e − ν L ( θ + η ) / sin ( η + θ )( η + θ ) × π (cid:90) π d φ sin ( η − θ − φ ) . (26)The φ average is simply 1 /
2. To compute the integrals over η and θ , we make the substitutions x = ν L η and y = ν L θ .We can then use the fact that ν L (cid:29) ∞ . F = k B T ργ (cid:90) ∞ dxdy e − ( x + y ) / sin (cid:0) x ν L + y ν L (cid:1)(cid:0) x ν L + y ν L (cid:1) . (27)lastic Response of Wire Frame Glasses. I. 2D Model 6This can be calculated for ν L (cid:29)
1, see appendix A. The finalanswer is F ≈ k B T γ ρ (cid:18) − ( ν L ) (cid:19) . (28)This directly leads to the linear elastic stress using (4), σ ( ν ) ∼ k B T ργ (cid:18) − ( ν L ) (cid:19) , (29)The important features of this expression are that it is approxi-mately independent of density and it agrees with a calculationbased on a reptation theory.In the case of rod shaped particles, it is possible to computethe elastic stress by an alternative argument. The main detailsof this textbook derivation are left to appendix B; here wesummarise the main results. The relevant component of theelastic stress is found to be, σ xy = k B T ργ (cid:40) e − D r t , for t < τ int e τ int ( D r ( ν ) − D r ) e − D r ( ν ) t , for t > τ int . (30)where D r and D r ( ν ) are the rod’s rotational diffusion coef-ficients in the absence and presence of constraints respec-tively. The latter is estimated using the reptation argumentsuch that, D r ( ν ) ∼ D r ( ν L ) − . As the average angle betweenconstraints is ∼ ( ν L ) − , the time taken for the rod to inter-act with the constraints is approximated by, D r τ int ∼ ( ν L ) − .This leads to the initial elastic stress, σ ( ν ) ∼ k B T ργ e − b / ( ν L ) ∼ k B T ργ (cid:18) − b ( ν L ) (cid:19) , (31)where b is a positive constant. This is precisely the same formas in (29), found by our geometric argument. The reptationargument does not predict the value of b exactly and, due tothe approximations used, we do not expect the value of 32 / / B. Wire Frame Particles
We now consider the simplest variant to a rod, a bent parti-cle with two legs joined at an angle χ . In the framework out-lined previously we have ∆ l = l χ in (11) and l = ,
1. The freeenergy is then calculated by essentially the same procedure asfor rods, the only difference being that we need to deal withdifferent possible combinations of constraints. This is doneby weighting the different contributions by their appropriateprobabilities, such that the free energy is, F = − k B T π ρ (cid:90) π d φ ∑ l , m = , (cid:90) π d η l P ( η l ) (cid:90) π d θ m P ( θ m ) (cid:90) πη l d η P ( η ) (cid:90) πθ m P ( θ ) log ( ˜ η l + ˜ θ m ) . (32) Making use of the same substitutions as before, x j = ν L η j and y j = ν L θ j , the integrals on η and θ can be evaluatedsimply, with each giving a factor of e − η l / and e − θ m / respec-tively. The expansion in (16), along with the definitions (17),are then used to expand the logarithm to second order in γ ,from which it is straightforward to show that the only contri-bution which will not vanish when integrated over φ is, k B T γ π ρ (cid:90) π d φ ∑ l , m = , (cid:90) ∞ dx l (cid:90) ∞ dy m ( ν L ) e − ( x l + y m ) ( x l + y m ) × (cid:20) sin (cid:16) φ + l χ − x l ν L (cid:17) − sin (cid:16) φ + m χ + y m ν L (cid:17)(cid:21) . (33)At this stage we can drop the subscripts on x and y , and takethe sums and φ integral to give, F = k B T γ ρ (cid:90) ∞ dxdy ( ν L ) e − ( x + y ) ( x + y ) × (cid:20) − cos χ cos (cid:18) ( x + y ) ν L (cid:19)(cid:21) . (34)At first, this integral appears to diverge because of the inclu-sion of the point x = y =
0. However, this point can be omittedfor physical reasons. In any real system there will be a finiteexcluded volume. This leads to the integral being cut off atsmall values, making it finite. The qualitative behaviour of theintegral can be understood by ‘pre-averaging’ the constraintsso that the angles, η and θ , are replaced by their average val-ues, (cid:104) η (cid:105) = (cid:104) θ (cid:105) = ( ν L ) − . This sets x ∼ y ∼ ν L gives thefree energy as, F = k B T γ ρ (cid:104) a ( ν L ) ( − cos χ ) + b cos χ (cid:105) , (35)where a and b are positive constants. In appendix C, we treatthe excluded volume cut off more rigerously, and find that thescaling found in (35) has only logarithmic corrections whichdepend on the exact cutoff. Therefore, without losing qualita-tive accuracy, we take the coefficients a and b to be unknownfitting parameters.The stress follows easily as, σ ( ν ) = ρ (cid:104) a ( ν L ) ( − cos χ ) + b cos χ (cid:105) k B T γ . (36)This can show a concentration dependence very different fromthat of rods, depending on the size of χ . For this differentbehaviour to be seen we require, ( ν L ) ( − cos χ ) (cid:38) cos χ . (37)After some simple re-arrangement this is,tan χ (cid:38) ν L , (38)which, since ν L (cid:29)
1, implies that if the particle is bentthrough an angle much larger than χ c ∼ ( ν L ) − , the domi-nant scaling of the stress will be ∝ ( ν L ) . This means thatlastic Response of Wire Frame Glasses. I. 2D Model 7even modestly bent particles ( χ ∼ O ( ) ) will have a signifi-cantly different elastic response than rods.In particular if we consider ‘L-shaped’ particles, where χ = π / σ ( ν ) ∝ k B T ρ ( ν L ) γ . (39)This is also observed for other symmetric, N-legged shapeswith the legs all joined at the same angle, 2 π / N . Examples ofthese particles are Y-shapes and crosses have been fabricatedin practice from double stranded DNA . It is interesting tonote that these special, symmetric shapes have no elasticity indilute suspension for fundamental symmetry reasons . Sus-pensions of these particle shapes therefore will be most sen-sitive to concentration changes. This has important implica-tions for designing these DNA systems as functional nano-materials, where it is desirable to have highly tuneable, con-centration or connectivity dependent properties.Equations (39) and (29) show that wire frames behave in astrikingly different way compared to rods. As ν L (cid:29)
1, wepredict that the elastic stress for wire frames is significantlylarger than that of rods, and very sensitive to the concentrationof constraining points. These results are qualitatively consis-tent with the simulations, see Fig(4). Here the simulation datafor the stress per particle is plotted as a function of densityfor rods, Y-shapes (N=3 stars) and X-shapes (N=4 stars) andcompared to the theoretical predictions (note that the data forrods has been multiplied by 20 for ease of visualisation). Thestress is calculated in the simulations by fitting the free energyto a quadratic for strains between γ = .
02 and 0 . ν , but the fitted exponent forthe data is noticeably smaller, ∼ .
4. In the simulations par-ticles which bend are removed from the system, i.e the bend-ing modulus K =
0. As will be shown later (section IV), thisleads to strong shear thinning in these systems even for rel-atively small strains, leading to the apparently smaller stressthan predicted by the linear theory.
IV. GENERAL METHOD: NON-LINEAR RESPONSE
We now extend our approach to study the non-linear elastic-ity of these systems. To calculate this we seek an expressionfor the free energy accurate to order γ . In addition to ex-panding the transformed entropy, log ˜ Ω , to this order, we needto introduce another physical process; the deformation of theparticles themselves. For example, as sketched in Fig.(5), anL shaped particle is forced to bend when the surrounding con-straining points are displaced by the shear deformation. Thissituation cannot occur for straight rods, but can for generalwire frame particles.To demonstrate this, consider two lines at angles α and β tothe x -axis with α < β . These represent the constraints placedon a particular particle; the particle is trapped between thesetwo lines. The particle will need to bend if, after the defor-mation, these lines cross. For this to happen we must be ableto find a pair of angles, α (cid:54) = β , for which ˜ α = ˜ β . Initially, FIG. 4. A log-log plot of the elastic stress per particle σ as a func-tion of density, ν , for rods, Y and X shaped particles. The points rep-resent the simulation data; circles for rods, triangles for Y-shapes andsquares for X-shapes. Note that the data for rods has been multipliedby 20. The dashed line is the theoretical result for rods (multiplied bytwenty) and the solid line the scaling theory for branched or bent par-ticles, with the intercept fitted. There is clear qualitative agreementbetween the theory and the data. For rods the stress is small and ap-proximately independent of density, whereas for branched particlesthe stress is much larger and increases rapidly with the density. Onthe other hand there is a quantitative discrepancy. The predicted scal-ing is ∝ ν , but the fitted exponent for the data is smaller ∼ .
4. Thiscan be explained by deviations from the linear elasticity at very smallstrains in these systems. take the particle to be a rod. This means that both lines willtransform in the same way, according to (15), which can bere-arranged to, cot ˜ α = cot α + γ . (40)Subtracting this from the equivalent expression for ˜ β , and in-sisting that the transformed angles need to be the same wefind that cot α = cot β , which implies, α = β + n π . Due tothe symmetry of the particle shape, rotation of one of the con-straining angles by an integer multiple of π results in the sameconstraint. Therefore, for rods, ˜ α = ˜ β if and only if α = β andhence, bending is impossible.For L-shaped particles, if the constraints are placed on dif-ferent legs, the transformation rule for one of the angles isdifferent. Taking β to be a constraint on one leg and α to beon the other, then β will transform as if it were rotated by π / α . Hence ˜ β is given by,cot ˜ β = γ − tan β . (41)Again, we subtract this from (40) and insist ˜ α = ˜ β whichyields, cot α = tan β . This has solutions, β = α − π / + n π .Clearly it is possible to find valid solutions with α (cid:54) = β and sofor these particles, bending is possible for certain configura-tions.To deal with the possibility of bending we introduce thefunction, P ( γ ; Γ ) . This is the proportion of particles with ori-entation Γ not bent at strain γ . The general form of this func-tion can be determined from simple considerations. When thesystem is not strained, none of the particles are bent, hence P ( Γ ) =
1. The behaviour of the system is symmetric un-der a change in sign of γ , which implies, P ( γ ; Γ ) = P ( − γ ; Γ ) .lastic Response of Wire Frame Glasses. I. 2D Model 8 FIG. 5. A sketch of the bending mechanism for an L-Shaped wire frame. (a) An undeformed configuration. The clockwise rotation ofthe L-shaped particle is constrained by the point, P , which hinders the leg, OP . The anti-clockwise rotation is constrained by the point P hindering leg OL . The L-shape can only rotate between the orange line, OP , and the green line, OP . The accessible configuration space isthe angle between these two lines, Ω . (b) When the step strain is applied, the points P and P are displaced to ˜ P and ˜ P respectively. TheL-shape must now lie between the lines O ˜ P and O ˜ P . The line O ˜ P is now to the right of O ˜ P , therefore the L-shape cannot keep its originalshape and ˜ Ω <
0. (c) To satisfy the constraints, the particle must deform from its original state. Assuming its two legs stay straight, the anglebetween them becomes π / − | ˜ Ω | . The size of strain is taken to be small so that this function isexpanded to give, P ( γ ; Γ ) = − γ / γ c ( Γ ) , (42)where we have defined the critical strain γ c ( Γ ) , whose scalingwith density will be determined in the following section. Thiscritical strain should be interpreted as the strain above whichall particles in orientation Γ are bent. This strain also indicatesroughly the limit of accuracy of the treatment presented here,as above this strain additional processes not taken into accountwill need to be considered, such as the non linear elasticity orpotential breakage of the particles themselves.We may also consider the orientational average of (42), P ( γ ) = − γ / γ c . (43)This should be interpreted as the total fraction of particleswhich have not bent at strain γ . We will use this for compari-son to simulations through the scaling of γ c with the density.When a particle bends, it stores some elastic energy. Thismust contribute to the free energy of the system and P ( γ ; Γ ) can be used to determine this bending contribution. The prob-ability that a particle first bends at a strain γ is − ∂ P / ∂ γ . Ifthe applied strain is γ , but a particle in orientation Γ firstbent at γ (cid:48) < γ , then the particle must bend by an angle ϑ =( γ − γ (cid:48) ) ∆ ( Γ ) , where ∆ depends on the particle shape. If thebending modulus of the particle is K , the energy associatedwith this bending is, K ϑ /
2. The bending contribution to thefree energy from one particle in orientation Γ is given by; F B ( γ ; Γ ) = K ρ ∆ ( Γ ) (cid:90) d γ (cid:48) ( γ − γ (cid:48) ) (cid:18) − ∂ P ( γ ; Γ ) ∂ γ (cid:12)(cid:12)(cid:12)(cid:12) γ (cid:48) (cid:19) , (44)where all the possible strains at which the particle could firstbend have been summed over, weighted by the appropriateprobability. We can now write an expression for the free en-ergy per unit volume which includes both the entropic contri- bution and the bending contribution. F ( γ ) = − k B T ρ (cid:90) d Γ ψ ( Γ ) P ( γ ; Γ ) (cid:10) log ˜ Ω ( Γ ; γ ) (cid:11) + K ρ (cid:90) d Γ ψ ( Γ ) ∆ ( Γ ) (cid:90) γ d γ (cid:48) ( γ − γ (cid:48) ) (cid:18) − ∂ P ∂ γ (cid:12)(cid:12)(cid:12)(cid:12) γ (cid:48) (cid:19) . (45)The orientation needs to be averaged over the distributionfunction ψ ( Γ ) , which for our case can be taken to be uniform.The first term in the above equation represents the entropiccontribution, F S ( γ ) . The factor of P needs to be included sincea particle which bends is forced to be in one configuration bythe constraints, and as such does not contribute entropically.The form of the free energy (45) is the basis of our treatmentof the non-linear elasticity in these systems. We will use thisto determine the free energy to order γ first for rods, then forL-shaped particles. The latter shape is taken for ease, but wewould find the same features for any bent or branched particle. A. Rods
The extension to non-linear elasticity for rods is straight-forward. As discussed previously, there is no possibility forbending, so P ( γ ; Γ ) =
1, and all that needs to be done is to ex-pand ˜ Ω and the entropy to order γ . After lengthy but straight-forward algebra, the first and third order contributions to thefree energy can be shown to vanish after integrating over φ .The second order contribution is that found previously in (28).The fourth order contribution is, after averaging on φ , (cid:28) γ sin ( η + θ ) ( η + θ ) (cid:104) ( η + θ ) cos ( η + θ ) + (cid:0) ( η + θ ) − (cid:1)(cid:105)(cid:29) . (46)The averages over η and θ are then taken in the same way asbefore, and the free energy to 4th order in γ is obtained, F ( γ ) ρ k B T = γ (cid:18) − c ( ν L ) (cid:19) − γ (cid:18) − d ( ν L ) (cid:19) , (47)lastic Response of Wire Frame Glasses. I. 2D Model 9where c and d are positive constants. The elastic stress is,dropping the subdominant terms for ν L (cid:29) σ = k B T ρ (cid:18) γ − γ (cid:19) . (48)This shows that the suspension of rods shear thins, as the co-efficient of the cubic term is negative. This is consistent withthe results of a more detailed theory for this system based onthe reptation model , and serves as a useful reference for com-parison to the behaviour of the L-shapes. B. L-Shapes
For L-shapes, bending is possible since configurations canbe found where ˜ Ω ≤
0. To compute P ( γ ; φ ) for these shapes,we first find the probability that a particle with orientation φ has bent, P B ( γ ; φ ) = − P ( γ ; φ ) . (49)This is the probability that ˜ Ω ≤ η < ˜ η , ˜ θ < ˜ θ and ˜ η + ˜ θ ≤ , (II) ˜ η < ˜ η , ˜ θ < ˜ θ and ˜ η + ˜ θ ≤ . (50)The probability of bending is then the sum of the integralsof the constraint probabilities over each of these regions. Wegive the calculation for region (I) in detail, the steps are thesame for (II) so we will simply state the result. The integralover region (I) is formally, P B ( γ ; φ ) = (cid:90) ˜ η + ˜ θ ≤ d η d θ P ( η ) P ( θ ) × (cid:90) ˜ η > ˜ η d η P ( η ) (cid:90) ˜ θ > ˜ θ d θ P ( θ ) . (51)To make progress we require expressions for the transformedconstraining angles, which are deduced straightforwardlyfrom previous definitions,˜ η = η + γ sin ( φ − η ) − γ sin ( φ − η ) cos ( φ − η ) , ˜ η = η + γ cos ( φ − η ) + γ cos ( φ − η ) sin ( φ − η ) , ˜ θ = θ − γ sin ( φ + θ ) + γ sin ( φ + θ ) cos ( φ + θ ) , ˜ θ = θ − γ cos ( φ + θ ) − γ cos ( φ + θ ) sin ( φ + θ ) . (52)Let us define the function A ( θ ) as the value of θ where˜ θ = ˜ θ . The function ˜ θ in (52) is monotonic in θ for smallvalues of γ , so it follows that if θ > A ( θ ) , then ˜ θ > ˜ θ . Wealso define B ( η ) as the value of η where ˜ η = ˜ η . Thesefunctions are found to be, A ( θ ) = θ + γ cos ( θ + φ )+ γ (cid:0) sin ( θ + φ ) − sin ( θ + φ ) (cid:1) , (53a) B ( η ) = η + γ cos ( η − φ ) − γ (cid:0) sin ( η − φ ) + sin ( η − φ ) (cid:1) . (53b)The integrals over θ and η in (51) are now, (cid:90) π B ( η ) d η P ( η ) (cid:90) π A ( θ ) d θ P ( θ ) . (54)However it can be shown that A ( ) < B ( ) < π + γ < φ < π + γ π + γ < φ < π + γ . (55)In which case, the lower limits on the integrals in (54) are bothreplaced by zero.In a similar fashion, we define the function η ( θ ) where˜ η + ˜ θ = η ( θ ) = − θ − γ cos ( θ + φ ) − γ (cid:0) sin ( θ + φ ) − sin ( θ + φ ) (cid:1) . (56)This defines a line in the θ , η , plane, which is approximatelya straight line between the points ( , η ( )) and ( θ ( φ ) , ) ,with gradient negative one, as shown in Fig.(6). The shadedregion in Fig.(6) enclosed by the axes and this line is the re-gion where ˜ Ω ≤
0. For certain orientations this region van-ishes, and bending is not possible. For the region to exist, itis necessary for η ( ) >
0, which from (56) occurs for orien-tations satisfying (55). Therefore, the contribution to P B fromregion (I) is determined from the integral, (cid:90) θ ( φ ) d θ (cid:90) η ( θ ) d η (cid:90) π d η d θ P ( η ) P ( η ) P ( θ ) P ( θ ) , (57)with φ satisfying (55) and, θ ( φ ) = − γ cos 2 φ − γ (cid:0) sin 2 φ + sin 4 φ (cid:1) . (58)Using the substitutions x l = ν L η l and y l = ν L θ l as before,the η and θ integrals will each evaluate to one and we areleft with,14 (cid:90) ν L θ ( φ ) dy (cid:90) ν L η ( y / ν L ) dx e − ( x + y ) / . (59)This integral is evaluated in appendix D. The result to lowestorder in γ , is,1 − e − ν L θ ( φ ) / − ν L θ ( φ ) e − ν L η ( ) / . (60)Using the definitions of θ ( φ ) and η ( ) , this is expanded tosecond order in γ to give the contribution to P B from region(I), (I) = , if φ / ∈ (55) (cid:16) ν L (cid:17) γ cos φ , if φ ∈ (55) . (61)lastic Response of Wire Frame Glasses. I. 2D Model 10 FIG. 6. A sketch of the θ , η plane. The line η ( θ ) with γ = . φ = π /
2, is shown. From (56) this must be approximately a linewith gradient −
1, which intersects the two axes at ( , θ ( φ )) and ( η ( ) , ) . The small variations about this line are too small to benoticeable on this figure. The shaded region is where ˜ Ω < The contribution from region (II) is calculated by followingthe same procedure, and is found to be,(II) = , if φ ∈ (55) (cid:16) ν L (cid:17) γ cos φ , if φ / ∈ (55) . (62)The bending probability is simply the sum of these to contri-butions, P B ( γ ; φ ) = (cid:16) ν L (cid:17) γ cos φ (63)Taking the orientational average of this, we obtain the propor-tion of un-bent particles, P ( γ ) = − ( ν L ) γ / . (64)From which the scaling of γ c can be extracted, γ c ∼ ν L . (65)Even though this was derived explicitly for L-shapes, thisscaling is expected for any bent or branched particles. InFig.(7) this scaling is compared to simulation data for Y andX-shaped particles and we see excellent agreement.The entropic term will need to be expanded to fourth orderas it was for rods but with P ( γ , φ ) included before the orienta-tional average is taken so that the integrals are only taken overthe region with ˜ Ω >
0. The exact calculation is possible albeitlong winded. Fortunately the scaling can be found using thesame pre-averaging method as in section III B (up to logarith-mic corrections, see appendix C). We replace the angles η and θ by their average values, and pre-average log ˜ Ω and P ( γ , φ ) over φ in (45). After several lines of algebra, we find, (cid:104) log ˜ Ω (cid:105) = − ( ν L ) ( ) γ − ( ν L ) ( ) γ . (66) FIG. 7. A log-log plot of the critical strain, γ c , against dimensionlessdensity, πν L . The points are data taken from the 2D simulations,triangles for Y-shaped particles and squares for X-shapes. The solidline is the scaling predicted theoretically (with the unknown interceptfitted to the Y-shape data). We see excellent agreement between thepredicted scaling and the data. The pre-averaged entropic free energy is then found by multi-plying this by (64),¯ F S ( γ ) ρ k B T = ( ν L ) ( ) γ + ( ν L ) ( ) γ − ( ν L ) ( ) γ . (67)Hence, the general form of the free energy to fourth order forbent and branched particles is, F S ( γ ) ρ k B T = a ( ν L ) γ − b ( ν L ) γ , (68)where a and b are positive constants.We now compute the bending contribution to the free en-ergy, which requires us to determine the angle through whichthe particle has bent at a particular strain. At the strain whenthe particle starts to deform, γ (cid:48) , the constraining points arepositioned on different legs at η = θ =
0. Therefore from(11), α = φ and β = φ + π /
2, so that the angle between thetwo legs is χ = β − α = π /
2. Straining the system further to γ > γ (cid:48) changes this angle to ˜ χ = ˜ β − ˜ α . This process bendsthe particle through an angle, ϑ = ˜ χ − π /
2. Using previousresults, it is straightforward to show that, ϑ = − ( γ − γ (cid:48) ) cos 2 φ − ( γ − γ (cid:48) ) sin 2 φ . (69)Only the first term is needed for the bending free energy tofourth order, which is given by F B ( γ ) = K ( ν L ) ρ (cid:90) π d φ π cos φ (cid:90) γ d γ (cid:48) ( γ − γ (cid:48) ) γ (cid:48) = K ( ν L ) ργ . (70)Finally, we obtain an expression for the free energy for wireframes accurate to order γ , which includes both the entropiclastic Response of Wire Frame Glasses. I. 2D Model 11and bending contributions, F ( γ ) ρ k B T = a ( ν L ) γ + b ( ν L ) (cid:18) Kk B T − c ( ν L ) (cid:19) γ . (71)The constants a , b and c are all positive and depend on theparticle geometry chosen. The stress immediately follows, σ ρ k B T = a ( ν L ) γ + b ( ν L ) (cid:18) Kk B T − c ( ν L ) (cid:19) γ . (72)This should be compared to the stress for rods (48). Not onlyis the stress for the bent particles significantly larger and moresensitive to concentration than for rods, the sign of the cubicterm is not necessarily negative. This means that while therod system will always be shear thinning, the bent particlesystem can shear harden. If the elastic modulus of the particlesthemselves is large enough, K (cid:38) ( ν L ) k B T , (73)the cubic term is positive and the response is shear hardening.Alternatively this condition can be seen as a critical densityabove which the behaviour transitions from shear hardeningto shear softening, ν c ∼ (cid:114) Kk B T L . (74)These conditions are non-trivial, because, while we have as-sumed K (cid:29) k B T , so the particles are very rigid, we also take ν L (cid:29)
1. If K / k B T ∼ ( ν L ) p for any p ≥
1, our results arevalid but only for p > P ( γ ) factor in the first term of equation(45). As the applied strain is increased, more and more parti-cles 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 tomake up this deficit, the total stress will be shear thinning.Therefore, there is some critical value of the bending moduluswhich must be exceeded to see a shear hardening response.Another difference between the behaviour of the wireframes and rods is the relative size of the non-linear term in thestress. For rods both the co-efficients of γ and γ in (48) are oforder one, implying non-linear effects are only important forlarger strains. On the other hand, the coefficient of γ in (72)is O (( ν L ) ) but the coefficient of γ is only O (( ν L ) ) . Thisshows that the non-linear effects can begin to be important atvery small strains. This can explain the difference betweenthe simulation data and the linear theory in Fig.(4). In thesimulations K =
0, so shear softening is expected, reducingthe measured stress compared to what is predicted from thelinear theory.
V. DISCUSSION
We have discussed the elastic stress response to step strainin a dense suspension of rod-like, bent and branched parti-cles using a simple geometrical method. The principle of thismethod is to determine the volume of configuration space ac-cessible to a particular particle, given the constraints placedon it by its surroundings, and how this volume changes whenthe constraints are moved by the applied strain. The accessi-ble volume of configuration space is related to the entropy ofthe particle, and so when it changes this leads to a change inthe free energy of the system. The stress calculated from thisfree energy should be interpreted as that measured immedi-ately after the step strain is applied.In this paper we illustrate this method on a simple 2D modelsystem. Here a test particle is surrounded by point-like con-straints which it cannot cross. These points play the role ofthe other particles in three dimensions by providing kineticconstraints. The test particle is free to rotate about its fixedcentre until it collides with the constraints. The constraintsare taken to transform affinely with the applied strain. We fo-cus specifically on the difference between straight rods andbent or branched wire frames, using L-shaped particles as anexample.The crucial difference between rods and the other particlesis that the constraints on rotation can be placed on differentlegs. The consequence of this is that the range of angles ac-cessible to an L-shape after the strain is applied can be zero oreven negative. This means the particle must have deformed tosatisfy the constraints. This must contribute to the free energy.This is taken into account by introducing the function, P ( γ ) ,interpreted as the proportion of particles which have not bentat a strain γ . We find that, to lowest order in γ , this is givenby P = − ( γ / γ c ) , where γ c is the critical strain above whichmost particles have bent. We determined that γ c ∼ ( ν L ) − ,which agrees well with the scaling found from simple simu-lations of the 2D model. With the inclusion of the bendingmechanism, the free energy was found to O ( γ ) and used todiscuss the magnitude of the linear elastic stress as well as thenature of the non-linear response for different particle shapes.The linear elastic stress for rods calculated by this methodagrees precisely with that found from a calculation based on areptation model.This simple model shows two interesting differences be-tween rods and wire frame particles:1) The linear elastic stress per particle for the wire framesscales proportional to ( ν L ) whereas for rods it is roughlyindependent of density. This is a significant difference. Aswe are taking ν L (cid:29)
1, the stress for wire frames is muchlarger than for rods and much more sensitive to concentration.This stronger scaling is found for any particle bent throughan angle ∼ O ( ) and so the rheology of the system is highlysensitive to the particle shape. This is in good agreement withsimulation data.2) There is a critical density beyond which the elastic stressfor branched particles changes from shear hardening to shearsoftening. The value of this critical density depends on theelastic modulus of the particle, ν c L ∼ ( K / k B T ) / . This islastic Response of Wire Frame Glasses. I. 2D Model 12also very different from the behaviour for rods, where the sys-tem is always shear thinning.We hope the simple model presented here may be a use-ful starting point for investigating other interesting features ofwire frame systems, such the decay timescale for the stress inthe system of branched particles. For such particles, the rep-tation process which allows rods to diffuse in a densely con-strained system is no longer possible. Hence we may expectto find very different time dependence of the stress and the2D model provides a framework for investigating this whereexact calculations are possible. ACKNOWLEDGMENTS
We gratefully acknowledge Prof. Daan Frenkel for manyimportant and insightful discussions. D.A.K. acknowledgesfinancial support from the UK Engineering and PhysicalSciences 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 γ is given by theintegral, F = k B T ργ (cid:90) ∞ dxdy e − ( x + y ) / sin (cid:0) x ν L + y ν L (cid:1)(cid:0) x ν L + y ν L (cid:1) . (A1)To compute this in the limit ν L → ∞ , we expand the inte-grand in the variable, y / ( ν L ) , which yields k B T ργ (cid:90) ∞ dxdy e − ( x + y ) / (cid:20) sin ( x ν L )( x ν L ) + C (cid:16) x ν L (cid:17) y ν L (cid:21) , (A2)where C ( a ) ≡ ( a )( a cos ( a ) − sin ( a )) / a . Integrating over y gives, k B T ργ (cid:90) ∞ dx e − x / (cid:20) sin ( x ν L )( x ν L ) + ν L C (cid:16) x ν L (cid:17)(cid:21) . (A3)Expanding the integrand in powers of x / ( ν L ) we finally ob-tain the free energy as given in equation (28) of the main text, F ≈ k B T ργ (cid:90) ∞ dx e − x / (cid:18) − x ( + x ) ( ν L ) (cid:19) ≈ k B T γ ρ (cid:18) − ( ν L ) (cid:19) . (A4) Appendix B: Calculation for rods based on reptation
Here we follow a textbook procedure to find the elasticstress for rods. Under x - y shear flow, a rod in two dimen-sions will rotate with angular velocity ˙ φ = − ˙ γ sin φ , about itscentre. If the rotational diffusion co-efficient is D r , then in the presence of this flow the orientational distribution function, ψ satisfies the Smoluchowski equation, ∂ ψ∂ t = D r ∂ ψ∂ φ + ˙ γ ∂∂ φ (cid:18) sin φ ψ (cid:19) . (B1)The elastic stress is found from the change in the free energyaccording to the virtual work principle (1). The free energy iswritten, F = k B T ρ (cid:90) d φ ψ log ψ . (B2)When the shear flow is applied, the distribution functionchanges and so does the free energy. The change in the distri-bution function is determined from the Smoluchowski equa-tion (B1). For a step shear of magnitude γ applied over a veryshort time δ t , the shear rate is ˙ γ = γ / δ t and the change in thedistribution function is, δ ψ = γ ∂∂ φ (cid:18) sin φ ψ (cid:19) , (B3)while the change in the free energy is, δ F = k B T ρ (cid:90) d φ δ ψ (cid:0) log ψ + (cid:1) . (B4)Substituting (B3) into (B4) and integrating by parts we find, δ F = − k B T ργ (cid:90) d φ sin φ ∂ ψ∂ φ , (B5)one final integration by parts yields the elastic stress, σ xy = k B T ργ (cid:10) sin 2 φ (cid:11) . (B6)We now obtain an equation for the average (cid:104) sin 2 φ (cid:105) . By defi-nition, the time derivative of this is, ∂∂ t (cid:10) sin 2 φ (cid:11) = (cid:90) d φ sin 2 φ ∂ ψ∂ t . (B7)Using the Smoluchowski equation and integrating by parts itis easy to show that, ∂∂ t (cid:10) sin 2 φ (cid:11) = − D r (cid:10) sin 2 φ (cid:11) − γ (cid:10) sin φ cos 2 φ (cid:11) . (B8)To calculate the linear elastic stress, the average in thefinal term can be replaced by its equilibrium value, (cid:10) sin φ cos 2 φ (cid:11) = − /
4. This allows (B8) to be solved easilyand using (B6) we find, σ xy = k B T ργ e − D r t . (B9)So far, we have not mentioned the constraints placed on therod which will clearly have an effect on the dynamics, alteringthe diffusion constant in the Smoluchowski equation.The constraints effectively force the rod to move along itslength in a channel of width a ∼ / ( ν L ) and length approx-imately L . For a short time, the rod is unaffected by theselastic Response of Wire Frame Glasses. I. 2D Model 13constraints, since it will not have diffused far enough to behindered by them. This is the case while it has not rotatedthrough and angle of about a / L ∼ / ( ν L ) . If the diffusionconstant without constraints is D r then this is true for timesup to of order τ int ∼ / ( D r ν L ) , as quoted in the main text.During this time period then, the stress decays according to(B9), but with D r = D r .For times longer than τ int , the rods orientation can onlychange via the reptation process. When the rod has diffusedroughly L / a / L . The timescale for this process is approximately τ ∼ L / D t , where D t ∼ D r L is the translational diffusion constant without con-straints. Therefore the rotational diffusion coefficient esti-mated for the reptation process is, D r ( ν ) ∼ a L τ ∼ a D t L ∼ D r ( ν L ) , (B10)as given in the main text. Hence for times longer than τ int , thestress decays as in (B9) but with D r = D r ( ν ) .Putting all of these pieces together and ensuring that thestress 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 rigorouslywith a small lower cutoff. This cutoff will ultimately be setby the excluded volume of the legs of the wire frame shapes.If the width of the legs is b , then the probability of finding aconstraining point with angles less than or approximately b / L must be zero due to the excluded volume effect. We define thesmall quantity, ε = b / L , and the reduced density p = ν L forconvenience. Thus the dominant term in (34) is, p (cid:90) πε d θ (cid:90) πε d η e − p ( θ + η ) ( θ + η ) . (C1)The integral is taken over a square with side length π with asmaller square with side length ε taken out around the origin.To make progress without changing the scaling we approxi-mate this region as a quarter circle of radius π with a smallquarter circle radius ε removed. The integral is then writtenin polar co-ordinates, by defining the vectors r = ( θ , η ) and d = ( , ) and the angle between them as ϕ , p (cid:90) πε drr (cid:90) π / d ϕ e −√ pr cos ϕ cos ϕ . (C2)Transforming the radial integral to the variable, R = pr , wehave, for large p , p (cid:90) ∞ p ε dRR (cid:90) π / d ϕ e −√ R cos ϕ cos ϕ . (C3)The R integral can then be evaluated as an incomplete gammafunction , p (cid:90) π / d ϕ cos ϕ Γ ( , √ p ε cos ϕ ) . (C4) Since we are taking the density to be in the range, 1 / L (cid:28) ν (cid:28) / V exc , we may take the limit νε →
0, giving − p (cid:90) π / d ϕ cos ϕ (cid:104) γ em + log ( √ p ε cos ϕ ) (cid:105) , (C5)where γ em is the Euler-Mascheroni constant. The integral cannow be taken exactly using integration by parts for the secondterm, − p (cid:104) tan ϕ (cid:18) γ em + + log ( √ p ε cos ϕ ) (cid:19) − ϕ (cid:105)(cid:12)(cid:12)(cid:12) π / = p (cid:0) | log p ε | − C (cid:1) , (C6)where we have used p ε (cid:28)
1, to make the sign more obviousand defined the positive, order unity constant, C .Hence we find only logarithmic corrections to the ν scal-ing discussed in the main text. This does not change the qual-itative features of our results. Appendix D: Evaluation of Bending Probability Integrals
To get the contribution to the bending probability from re-gion (I), see (50), we need to evaluate the integral,14 (cid:90) ν L θ ( φ ) dy (cid:90) ν L η ( y / ν L ) dx e − ( x + y ) / . (D1)Where the functions, θ ( φ ) and η ( y / ν L ) are defined inequations (58) and (56) of the main text respectively. Takingthe x integral gives,12 (cid:90) ν L θ ( φ ) dy (cid:18) e − x / − e x / − ν L η ( y / ν L ) / (cid:19) . (D2)Integrating the first term in the integrand is straightforwardwhereas the second integral is schematically, (cid:90) c γ + d γ dx e ν L γ a ( x / ν L )+ ν L γ b ( x / ν L ) . (D3)The integrand can be expanded in powers of x / ν L , becauseover the whole integration range x / ν L ∼ O ( γ ) , (cid:90) c γ + d γ dx e ν L γ a ( )+ ν L γ b ( ) (cid:104) + x (cid:0) γ a (cid:48) ( ) + γ b (cid:48) ( ) (cid:1)(cid:105) , (D4)clearly the integral of the second term in the square bracketswill be O ( γ ) , and so only the first term contributes to secondorder. Therefore and the result of the integral (D2) is,1 − e − ν L θ ( φ ) / − ν L θ ( φ ) e − ν L η ( ) / , (D5)as given in equation (60) of the main text.lastic Response of Wire Frame Glasses. I. 2D Model 14 DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request. N. C. Seeman, “Nucleic acid junctions and lattices,” Journal of TheoreticalBiology (1982), 10.1016/0022-5193(82)90002-9. P. W. Rothemund, “Folding DNA to create nanoscale shapes and patterns,”Nature , 297–302 (2006). F. Bomboi, D. Caprara, J. Fernandez-Castanon, and F. Sciortino, “Cold-swappable DNA gels,” Nanoscale , 9691–9697 (2019). 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 , 15633–15637 (2013). Z. Xing, A. Caciagli, T. Cao, I. Stoev, M. Zupkauskas, T. O’Neill, T. Wen-zel, R. Lamboll, D. Liu, and E. Eiser, “Microrheology of DNA hydrogels,”Proceedings of the National Academy of Sciences , 8137–8142 (2018). P. Gross, N. Laurens, L. B. Oddershede, U. Bockelmann, E. J. Peterman,and G. J. Wuite, “Quantifying how DNA stretches, melts and changes twistunder tension,” Nature Physics , 731–736 (2011). W. Van Ketel, C. Das, and D. Frenkel, “Structural arrest in an ideal gas,”Physical Review Letters , 8–11 (2005). M. Doi, “Rotational relaxation time of rigid rod-like macromolecule in con-centrated solution,” Journal de Physique , 607–611 (1975). M. Doi and S. F. Edwards,
The theory of polymer dynamics (Oxford Uni-versity Press, 1986). P. G. D. Gennes, “Reptation of stars,” Journal de Physique , 1199–1203(1975). D. S. Pearson and E. Helfand, “Viscoelastic Properties of Star-Shaped Poly-mers,” Macromolecules , 888–895 (1984). M. Doi and N. Y. Kuzuu, “Rheology of star polymers in concentrated solu-tions and melts,” Journal of Polymer Science: Polymer Letters Edition ,775–780 (1980). S. T. Milner and T. C. McLeish, “Parameter-free theory for stress relaxationin star polymer melts,” Macromolecules , 2159–2166 (1997). D. R. Heine, M. K. Petersen, and G. S. Grest, “Effect of particle shape andcharge on bulk rheology of nanoparticle suspensions,” Journal of ChemicalPhysics (2010), 10.1063/1.3419071. M. K. Petersen, J. M. D. Lane, and G. S. Grest, “Shear rheology of ex-tended nanoparticles,” Physical Review E - Statistical, Nonlinear, and SoftMatter Physics , 1–4 (2010). M. Doi and S. F. Edwards, “Dynamics of rod-like macromolecules in con-centrated solution. Part 1,” Journal of the Chemical Society, Faraday Trans-actions 2: Molecular and Chemical Physics , 560–570 (1978). S. F. Edwards and T. Vilgis, “The dynamics of the glass transition,” PhysicaScripta , 7–16 (1986). S. F. Edwards, “Statistical mechanics with topological constraints: I,” Pro-ceedings of the Physical Society , 513–519 (1967). S. F. Edwards, “Theory of cross-linked polymerized material,” Journal ofPhysics C: Solid State Physics , 1–13 (1969). S. F. Edwards and K. F. Freed, “The entropy of a confined polymer. I,”Journal of Physics A: General Physics , 145–150 (1969). D. A. King, M. Doi, and E. Eiser, “Particle shapes leading to Newtoniandilute suspensions,” Phys. Rev. E , 32615 (2020). M. Abramowitz and I. A. Stegun,