{'fraction_non_alphanumeric': 0.033515715091886056, 'fraction_numerical': 0.011924332016586109, 'mean_word_length': 4.44383598236944, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Equilibrium and dynamic pleating of a crystalline bonded network

We describe a phase transition that gives rise to structurally non-trivial states in a two-dimensional ordered network of particles connected by harmonic bonds. Monte Carlo simulations reveal that the network supports, apart from the homogeneous phase, a number of heterogeneous"pleated"phases, which can be stabilised by an external field. This field is conjugate to a global collective variable quantifying"non-affineness", i.e.~the deviation of local particle displacements from local affine deformation. In the pleated phase, stress is localised in ordered rows of pleats and eliminated from the rest of the lattice. The {\em kinetics} of the phase transition is unobservably slow in molecular dynamics simulation near coexistence, due to very large free energy barriers. When the external field is increased further to lower these barriers, the network exhibits rich dynamic behaviour: it transforms into a {\em metastable} phase with the stress now localised in a {\em disordered} arrangement of pleats. The pattern of pleats shows ageing dynamics and slow relaxation to equilibrium. Our predictions may be checked by experiments on tethered colloidal solids in dynamic laser traps.


I. INTRODUCTION
Fabricating complex shapes by folding or pleating a two-dimensional elastic manifold has recently emerged as a viable technological paradigm, applicable over a large range of length scales, from microns to nanometers [1][2][3][4]. A number of these innovative ideas are equally applicable for atomic crystals or for larger assemblies involving functionalized colloidal particles joined together using polymer tethers [5] or micron-sized lipid droplets [6,7]. To make such attempts feasible, one needs efficient ways to control local structural properties, preferably in a reversible way. Thus, microscopic understanding of the underlying thermodynamics and kinetics of these or similar local shape changes would be valuable.
In this paper, we study in detail equilibrium and dynamic aspects of a transition in a tethered network of colloidal particles -from a homogeneous to a heterogeneous phase containing an ordered arrangement of pleats. Within these pleats, the network penetrates itself, with parallel rows of vertices folding back to completely or nearly overlapping adjacent rows. These complex structures arise here spontaneously as a result of an underlying, equilibrium "first-order" phase transformation [8].
We show that pleating excitations are induced and controlled by a novel external field conjugate to a collective variable [9] defined as follows. In earlier work, it was demonstrated that local displacements of particles in a crystal away from their ideal positions may be decomposed into affine and non-affine components [10][11][12]. External stress couples to the affine part of the displacements. In analogy, one imagines an external field, which is conjugate to a global collective coordinate X, mea-suring non-affineness, to be defined explicitly later. This "non-affine" field, h X , is realisable experimentally for colloidal solids using dynamic laser traps [11,13]. For small h X > 0, non-affine displacements, and consequently X, are enhanced in a controlled manner computable within a linear response framework. In a non-bonded conventional crystal, this leads to the creation of defects [11]. The present paper is devoted to an analysis of the consequences of large positive h X in a connected network much beyond the linear response regime.
We consider a periodic lattice of point vertices in two dimensions (2d) that are connected to their nearest neighbours by harmonic springs. Crystalline 2d networks have been studied extensively in the past [14][15][16], partly because of their biological significance (e.g. as a simple model for the spectrin network in red blood corpuscles) [17,18]. Note that such a network has nontrivial properties, such as the presence of pleats, only when the bond length is random or exceeds the nearest neighbour distance set by the density [16] by a critical amount. In this case the network suffers an instability and catastrophic collapse. This is in stark contrast to the equilibrium first order transition [8] in the presence of h X (and zero external stress) to an ordered pleated state that we describe, we believe for the first time, in the present work. For this transition, we obtain the relative free energies of the pleated phases as well as individual free energy barriers in quantitative detail -a first but essential step towards efficient control over their formation. We use two distinct particle-based simulation techniques for this study.
Monte Carlo (MC) simulations [19] in combination with sequential umbrella sampling (SUS) [20] are used to compute the free energy land-arXiv:1612.00574v1 [cond-mat.soft] 2 Dec 2016 scape as a function of X at various values of h X . We show that SUS-MC simulations are able to detect metastable phases which are inaccessible by conventional MC. Free energies of interfacial structures between pleated and unpleated regions of the lattice can also be investigated and energy barriers for the formation of the product phase determined. This advanced sampling method (SUS-MC) allows us to not only identify this as a first order phase transition but also to characterise the properties of the coexisting states separated by interfaces. This latter point is remarkable because we obtain explicitly an interface between an inhomogeneous state (i.e. the pleated state) and the normal crystal. Locating this interface is nontrivial and requires the computation of local stresses.
We use molecular dynamics (MD) simulations [21,22] at constant particle number, total area and temperature (NAT) to reveal the kinetic aspects of this transition. In the MD simulation, the transformation shows features that are different from the thermodynamic first-order transition, as seen in SUS-MC. The transformation occurs at a larger value of h X and leads to a metastable phase. We show, in this work, the strong interplay between kinetics and thermodynamic phase behaviour. The equilibrium pleated states are heterogeneous phases which are quite difficult to realise via a kinetic pathway. Conversely, the nature of the metastable states obtained in the kinetic transition can only be understood if one is able to identify the underlying phase transition through SUS-MC.
We have carried out the simulations for both the pure, non-self avoiding network as well as a model where the vertices of the network are decorated with finite-sized colloidal particles. On a qualitative level, the findings for both models are similar.
The rest of the paper is organised as follows. Section II introduces the local and global non-affine parameters, χ(R) and X, as well as the field h X conjugate to the latter. Moreover, the details of the model solids and simulation methodology are presented here. Then, in Sec. III, a comprehensive exposition of our results for the network solid is given. We report results for analytic calculations for small h X , ground states and finite temperature phases and the dynamical transition in this model. We end the paper with a summary and conclusions as well as an outlook for future work.

II. MODELS, FORMALISM AND SIMULATION DETAILS
In this section, we commence our discussion by first introducing the non-affine field h X and the model Hamiltonians followed by a description of the simulation methodologies used. A. The model Hamiltonian with the non-affine field Consider a reference configuration of N particles where the particle with index i (i = 1, ..., N ) is located at position R i . A displacement of particle i from its position on the reference lattice is given by u i = r i − R i , with r i the instantaneous position of the particle. Now within a neighborhood Ω around particle i, we define relative atomic displacements ∆ j = u j − u i with particle j = i ∈ Ω. The "best fit" [9] local affine deformation D is the one that minimizes j [∆ j − D(R j − R i )] 2 with the nonaffinity parameter χ(R i ) > 0 being the minimum value of this quantity. The minimisation procedure amounts to projecting [10] ∆ i onto a non-affine subspace defined by a projection operator P such that χ(R i ) = ∆ T P∆ where ∆ is the column vector constructed out of the ∆ i . In the projector P = I−R(R T R) −1 R T , the N d×d 2 elements of R are given by R jα,γγ = δ αγ R jγ (here, the central particle i is taken to be at the origin). This projection formalism is perfectly general and can be carried through for any lattice in any dimension.
Once the dynamical matrix D(q) = v −1 BZ u q u −q is obtained (v BZ is the volume of the Brillouin zone), we can calculate the ensemble average of the non-affine parameter χ using a coarse graining procedure outlined in [10,11]. We include this in brief here for completeness.
We define the coarse grained correlations, C iα,jβ = ∆ iα ∆ jβ where as before the Roman indices denote particles and the Greek ones denote coordinates. The particles i, j both belong inside the neighborhood Ω of particle 0. Substituting the definition of the displacement differences ∆ iα = u α i − u α 0 we obtain [10], The ensemble average of χ is then given by Tr PCP, where P is the projection operator. The non-trivial eigenvalues σ µ of PCP and their corresponding eigenvectors are the non-affine modes of the lattice. For example, in the triangular lattice there are 8 such modes when Ω corresponds to the nearest neighbour shell. The orthogonal subspace, i.e. the affine displacements, are spanned by the non-trivial eigenvectors of (1 − P)C(1 − P) and correspond to the usual volumetric, uniaxial and shear strains together with local rotations. The non-affine field does not affect the statistics of the affine part of the displacements to linear order in h X . Space-time correlation functions of both affine and non-affine variables can also be obtained using a similar procedure [10,11]. In order to selectively excite lattice distortions that enhance non-affine displacements, we introduce an extended microscopic Hamiltonian [11] involving the thermodynamic conjugate variables h X and X, with X = N −1 N i χ(R i ). In analogy to conjugate variables like stress-strain or pressure-volume, we add the product of h X and X to the Hamiltonian: (1) Here, H 0 represents the Hamiltonian of any standard solid.
The second term in Eq. (1), −N h X X, involves appropriate Cartesian components P ij of the projection operator P that couple to the relevant displacement differences. Note that the P ij are constant parameters that only depend on the position vectors of the reference configuration, R i . The size of the coarse-graining volume Ω surrounding particle i is set by the range of the interaction. In the rest of this paper, we take Ω as the nearest neighbour shell. A positive value of the non-affine field h X enhances non-affine distortions of the lattice [10,11], namely, fluctuations in particle positions projected onto a subspace spanned by those eigen-distortions of Ω that cannot be represented as combination of affine deformations.
We describe below our network model. The reference lattice structure {R i } i=1,...,N is an ideal triangular lattice. The Hamiltonian of this model is that of a standard network of point non self-avoiding vertices connected by harmonic bonds [14], with p i the momentum, m the mass, r i the instantaneous position, and R i the reference position of vertex i as before. The length scale is set by the lattice parameter l, the energy scale by Kl 2 , and the time scale by m/K. We use those as our units in the following, effectively setting l = m = K = 1. A dimensionless inverse temperature is given by β = Kl 2 /k B T , with k B the Boltzmann constant.
For some of the calculations reported here, we attach finite sized repulsive particles with every vertex. The Hamiltonian is therefore augmented to The interaction potential for a pair of particles, separated by a distance r, is (3) for r r c = 2 1 6 r 0 and v WCA = 0 for r > r c . We use φ = 1 and r 0 = 0.6l, respectively.

B. Simulation Details
We perform molecular dynamics simulations as well as Monte Carlo in combination with sequential umbrella sampling of the regions of configuration space that are otherwise inaccessible using simple MC or MD simulation techniques.
Molecular Dynamics. MD simulations in the canonical ensemble, i.e. at constant number N of vertices, area A and temperature T , were done using a leapfrog algorithm, coupling the system to a Brown and Clarke thermostat [21,22]. The size of the systems ranges from N = 100 to N = 40000 vertices. Typically, unless otherwise stated, we used an MD time step of δt = 0.002 and inverse temperature β = 200. These parameters are the same even when repulsive, WCA particles are attached to the vertices. Typically, the solid was held for 5 × 10 5 MD steps followed by the collection of data for further 5 × 10 5 MD steps at an interval of 1000 steps. Sequential Umbrella Sampling. Standard Metropolis Monte Carlo [19,22] is inefficient in sampling systems with free energy barriers separating different regions in configuration space. SUS is an advanced sampling technique [19,20,22] that ensures good sampling of the entire range of pertinent states. Our implementation of SUS-MC in the N AT ensemble involves dividing the range of the relevant order parameter, X, into small windows to be sampled successively starting at X = 0. Histograms denoted by H(n) keep track of how often each value of X within the nth window is realized, with H(n) L and H(n) R representing the left and right boundaries of the nth histogram. Now, a predetermined number of MC moves are attempted per window. MC moves resulting in X within the chosen window are accepted or rejected using the conventional Metropolis criterion and the relevant histograms are modified accordingly. Any moves leading to values of X outside the chosen window are rejected with appropriate modification of the histograms at the boundaries to ensure detailed balance [20]. Finally, the un-normalized relative probability distribution of X can be computed using, The SUS-MC runs, for the network model, were done for systems with N = 900 at β = 200 and density ρ = 1.1547 (which corresponds to our choice l = 1 for the lattice parameter). For systems with N = 900, we considered the range between X = 0.0 and X = 1.0 and divided this into 500 sampling windows with 8×10 7 MC moves attempted in each window. In each MC move, maximal particle displacements of 0.2 l along the x and y directions are allowed. Apart from simulations with N = 900 vertices, systems with N = 100 and N = 400 were studied using SUS-MC. As the variance of the parameter X is proportional to inverse of the system size [11], the range of X for sampling was chosen accordingly, viz. 0.0 ≤ X ≤ 3.0 for as a function of hX , corresponding to the two degenerate non-affine modes at the lowest energy (filled squares). The solid curve is a linear fit through the last three data points. The lattice distortions corresponding to these modes are illustrated in the inserted sketches. Note that in each case nearest neighbours a, b (red lines) move apart while next nearest neighbours, c, d (blue lines) come closer (see text). N = 100, 0.0 ≤ X ≤ 1.5 for N = 400, and 0.0 ≤ X ≤ 1.0 for N = 900. Also for N = 100 and N = 400, the range of X was divided into 500 sampling windows with 8×10 7 MC trial moves in each window.

A. Crystal properties in the presence of small hX
We first study the properties of our model solid in the low temperature limit, where a harmonic approximation becomes exact even in the presence of the non-affine field h X . Indeed, the complete low-T statistical mechanics of the system can be obtained as long as the periodic crystalline phase is stable. Therefore we begin by first presenting analytic results for finite values of h X [10][11][12] in the ideal crystal. These calculations provide an estimate of the critical value, h 0 X , at which the crystal becomes unstable under the application of the non-affine field.
In order to understand how crystal properties are affected by the non-affine field, h X , we calculate the eigenvalues of the dynamical matrix [23,24] and the corresponding eigenvectors. The dynamical matrix is obtained from the Fourier transform of the Hessian D(R, R ) = ∂ 2 H/∂u(R)∂u(R ). The full expression for D in the presence of a non-affine field h X can be worked out in a straightforward manner and one obtains, Note that the leading order term in A X (q) is of order O(q 4 ), so that h X does not contribute to the speed of sound or to elastic constants.
Increasing h X enhances χ locally, which is associated with a softening of certain phonon modes. At small q, the transverse phonon modes are softened showing that the solid becomes nearly unstable to large wavelength shear modes. We find that most of this contribution comes from the two softest non-affine modes, i.e. the eigenvectors of PCP with the smallest eigenvalues σ µ (µ = 1, 2). These two eigenvalues are identical and the corresponding modes are illustrated in the sketches of Fig. 1: distortions with respect to the reference configuration (open circles) are generated where nearest neighbours a, b move away from each other while next-nearest neighbours c, d come closer tending to nucleate a dislocation dipole [11]. Note that any linear combinations of these modes are also degenerate and so the direction of the modulation wavenumber, with magnitude χ −1/2 , varies in space pointing randomly in all directions consistent with crystal symmetry.
In Fig. 1, we plot the reciprocal of the largest nonaffine eigenvalue (σ µ with µ = 1, 2) as a function of h X . As h X increases, this eigenvalue vanishes as h X − h 0 X pointing to an underlying saddle-node bifurcation point beyond which a crystalline solid cannot exist [10]. Similar behavior is shown by all the eigenvalues (not shown) corresponding to the non-affine modes [11]. As we can estimate from the plot, the eigenvalues σ µ with µ = 1, 2 diverge around h 0 X = 0.072 -the limit of stability of the crystal in a non-affine field.

B. Pleated configurations at T = 0
Having considered the ideal crystalline ground state of the non-self-avoiding triangular network, we now ask whether additional low-energy configurations exist. Indeed, we show that configurations containing one or more pleats, as illustrated in Fig. 2, are also possible. Within a pleat, two rows of vertices overlap producing a band of twice the local stiffness. Note that the pleated configuration remains two-dimensional. While similar pleated structures have been reported [14][15][16] for such networks under compression or with disorder, to the best of our knowledge, the existence of these states for regular crystalline networks has never been commented upon before. As explained in Fig. 2a., in a pleated state no bond is stretched or compressed. A pleated row of vertices does not destroy local crystalline order and can be distinguished only by a high value of local χ, pointing out that a finite non-affine displacement is necessary to produce a pleat. The displacement u however becomes a multi-valued function of the coordinates at the location of the pleats with each of the two values corresponding to the two distinct "leaves" of the pleat. Since the lattice is stiffer at the pleats these regions are also regions of enriched stress.
A positive h X encourages the formation of pleats. In a solid of size L, internal strains ε int ∼ 1/L, however, need to be introduced to fit the the pleated lattice back into the box, making configurations with a large number of pleats energetically unfavourable. Consider, for the moment, only pleats of full rows, all in the same direction, then having p pleats requires a strain of p/L elsewhere, hence elastic energy ∼ L 2 (p/L) 2 ∼ p 2 . However, we gain an energy ∼ h X pL from the field term as χ is increased to O(1) in an area of order pL. Equating the two predicts p ∼ h X L or p/L ∼ h X , i.e. a finite fraction of pleats that increases linearly with h X . At h X = 0 one has competition between energy ∼ p 2 and entropy ∼ p ln L and thus expects p ∼ ln L, hence p/L → 0 as L → ∞.
In Fig. 2b, we illustrate this with two pleats, one horizontal and the other tilted at an angle of 60 • . We create these configurations by shifting rows of vertices either downwards or to the left by amounts λ 1 and λ 2 respectively. Here λ 1 = λ 2 = 2 lattice spacings. The local χ is a quadratic function of both λ 1 and λ 2 . Note that internal strains and the term proportional to h X in H determine the relative stability of the pleated configurations. The network responds to h X by either increasing the density of the pleats or by introducing side branches at 60 • (or equivalently 120 • ) to the horizontal. Large h X favours a large number of pleats and in that regime several configurations may have the same average non-affineness X and, at the same time, be degenerate in energy. In the next section we show that these pleated configurations of the crystalline network survive at non-zero temperatures and lead to interesting phase behaviour.
C. The phase transition at T > 0 In this section we study in detail the phase transition from an un-pleated solid to one with pleats at finite tem- a. Schematic illustration of a pleat in the triangular lattice as an origami fold. A "valley" fold at I − II and a "mountain" fold at III − IV bring the vertices at I, II and 1, 2 as well as III, IV and 3, 4 together to obtain the final two-dimensional configuration shown on the right. Note that the pleated configuration has two layers less and 7 overlapped bonds together with 5 overlapped vertices. Such an operation can be continued along a crystalline row of vertices creating a complete pleat. Note further that none of the bonds in the pleat are either stretched or compressed. b. An originally 18 × 32 lattice with one horizontal and one 60 • pleat of amplitudes λ1 = λ2 = 2 lattice spacings. The vertices are shown as filled circles coloured according to the value of local χ. The pleats do not destroy local crystalline order but have large local χ values (color bar).
peratures. This transition can be located by SUS-MC which gives at a given value of h X the probability P (X) to find the system in a state with a certain value of X. The logarithm of this probability is directly related to the corresponding free energy, F (X) = −k B T ln P (X) + C (with C a constant). Figure 3a displays − ln P (X) for the system with N = 30 × 30 lattice sites at different values of h X . All the distributions exhibit a first minimum around X = 0.05 consistent with our results for the ideal network in the harmonic approximation [11]. A saddle point also appears at X saddle ≈ 0.1 which is almost independent of h X . For X > X saddle , the function − ln P (X) first decreases until a second minimum is reached. At higher values of X more minima and shoulders can be discerned in the distributions; at high values of h X these lie below the first two minima and thus correspond to the stable states in this regime.
In Fig. 3b, P (X) at h X = 0.025 for various system sizes N is plotted. In a conventional first-order phase transition [8], P (X) sharpens with increasing system size such that −N −1 ln(P (X, N )) approaches its thermodynamic limit. In our system, while P (X) does become narrower with increasing N , detailed features of P (X) depend nontrivially on the system size. Indeed for larger system sizes, additional minima appear that correspond to domain configurations not possible for smaller sizes. One of the great advantages of the SUS-MC method is that configurations that contribute to P (X) in each range are directly available.
Snapshots for different values of X at h X = 0.030 (as marked in Fig. 3a by black dots) are displayed in Fig. 3c. Here each vertex is represented as a filled circle and coloured according to the local value of log 10 χ. While the configuration corresponding to the first minimum in − ln(P (X)) is a homogeneous crystal, the one corresponding to the second minimum is an inhomogeneous phase where a band of vertices with a high positive χ value percolates through the crystal. The latter band does not simply consist of two lattice rows, as the snapshot may suggest. While the two rows fit perfectly into the hexagonal structure, each of the two rows consist of a pair of overlapping rows, i.e. a pleat. At higher values of X, in addition to the horizontal pleat, side branches form at an angle of 60 • with a geometry similar to that used in the T = 0 calculation. On further increasing X, the side branches cross the main band such that they also percolate through the system. Is there a finite temperature phase transition from a homogeneous crystalline phase to an inhomogeneous phase with a pleat of "non-affine" vertices? At such a transition, the probability distribution P (X) would exhibit two peaks with the area under both peaks being equal [19]. As one can infer from Fig. 3a, this happens for a value of h X between 0.025 and 0.030. To obtain an estimate of the coexistence value for h X , one can use histogram reweighting and deduce from a reference dis- a. − ln P (X) is plotted as a function of X at hX = 0.027, where the crystal and the inhomogeneous phase with a non-affine stripe of particles are at coexistence. b. Configuration snapshots with log 10 χ colormap corresponding to the states indicated in a. While the snapshots C1, C2, C3 represent instantaneous configurations, the snapshots M1 and M2 were obtained from a time average over 10 4 subsequent configurations.
(4) Eq. (4) may now be used to determine the distribution P (X, h coex X ) for which the area under the peaks corresponding to the coexisting phases is equal. We accomplish this using an iterative procedure where successive refinements of h coex X are obtained from SUS-MC simulations at a previously estimated values of h coex X . From this procedure we find h coex X = 0.027. We check that the final P (X, h coex X ) shows two peaks with equal area under them as required for co-existence (see Fig. 4a). States in the two-phase region (e.g. those points marked as C1, C2, and C3 in the figure) are also visible. Here, the two phases are at coexistence and separated from each other by an interfacial region. As the snapshots for the states C1, C2, and C3 show the pleat does not percolate through the system in the two-phase region and remains as a droplet terminated at two opposite vertices by the presence of the homogeneous crystal phase. Local χ is a convenient collective variable useful for characterising pleats. However, the thermodynamic variable which ought to show the interface between the homogeneous crystal and the inhomogeneous pleated phase is the space dependent, complex amplitude of the appropriate vertex density modulation. This is computationally difficult to obtain. Fortunately, as we show below, the local (uniaxial) stress distribution works just as well.
In Fig. 5 we demonstrate that the two phases at coexistence (Fig. 4) differ from each other with respect to the distribution of local stresses [22] at coexistence, P (σ), where σ = σ xx − σ yy , the deviatoric or uniaxial stress component. The distribution for the inhomogeneous phase (M2) has its maximum at the negative value σ = −4 × 10 −5 and is asymmetric, with a pronounced excess contribution at positive values of σ. The latter asymmetry reflects the localization of stresses in the pleat with a high value of log 10 χ (see snapshot corresponding to M2 in Fig. 4). In contrast, the distribution corresponding to M1 is symmetric and its maximum is at the positive value σ = 2 × 10 −5 . The corresponding snapshot indicates a homogeneous distribution of stresses. Thus, the two phases at M1 and M2 differ with respect to the average stress value and the localization of stress. In the M2 phase, the average stress value in the pleated region where the stress is localized is similar to that of the homogeneous M1 phase. All this is reflected in the snapshots corresponding to the two-phase region (C1, C2, and C3). Now the interfaces between the coexisting phases are clearly visible (dashed red lines in the plot of configuration C1, C2 and C3). While in the M1 phase stress is distributed throughout the crystal, in the M2 phase, it is eliminated from the un-pleated part and concentrated mostly at the pleat. There are two interfaces due to periodic boundary conditions and the amount of the two phases at a given state is controlled by the lever rule. Consequently, the free energy decreases linearly from state C1 to C3 (see Fig. 4), although the total area of the interfaces is constant for these three states. In the periodically repeated system, the M2 phase becomes a vertical stripe of undistorted network punctuated by a parallel array of horizontal pleats at regular intervals. Thermal undulations along the interfaces are also visible. In summary, Fig. 5 demonstrates that the order parameter of the transition from the homogeneous crystal to one with stress localisation within pleats can be connected to distribution of local stresses. We end this section with a comment on P (σ) from our simulations. Fig.5 suggests that the σ → −σ symmetry is broken at ε = 0. To obtain a complete description one would need to sample all degenerate, globally rotated copies of the crystal and also evaluate the full stress tensor σ ij . This is however not necessary for our purpose here.

D. Dynamical transition and plastic deformation
In the last section we studied the properties of the pleated configurations showing that they form the stable equilibrium phase beyond a first order transition from an un-pleated to a pleated phase, which occurs in our network solid at h X = h coex X ≈ 0.027 when β = 200. How can such configurations form dynamically? We turn now to study this kinetic transition.
We have mentioned before that pleated configurations imply a multi-valued displacement field. Obviously, such a configuration cannot be represented as a linear combination of hydrodynamic phonon fluctuations of the net- work. Pleated configurations therefore need to form by the nucleation and growth of non-hydrodynamic, localized droplets. Incomplete pleated regions surrounded by a strained network have been observed and described in detail in the last subsection (see for example the C1 configuration in Fig. 4). For values of h X at which pleated states become globally stable, such droplet configurations, which lie in the saddle region in-between minima corresponding to un-pleated and pleated states cost extremely high free energy. Such high barriers prevent the equilibrium transition from occurring in MD simulations.
A dynamical transition at T > 0 to a pleated configuration is possible only when h X becomes sufficiently large so that the lattice is close to being (but not quite!) locally unstable and the free energy barrier is substantially reduced. In Fig. 6b, we have plotted h X against X obtained from SUS-MC. The X values were obtained by a histogram reweighting method. Together with these results, we have also plotted results from MD simulations of the same network where X now represents an average over the MD simulation time. The MD and the SUS-MC results both show a jump in X at the pleating transition. However, the transition in MD occurs at a much larger value of h X showing that for a large range of h X , the un-pleated state remains metastable. Note also that the value of X in the un-pleated state just before the dynamical transition is roughly equal to the value at the saddle point as shown in Fig. 3. A Lindemann like criterion [8,25] viz. X = X L = X saddle just below the transition in the crystal phase is thus operative at this kinetic transition. Finite size effects in the MD simulations roughly follow those in the MC consistent with the shift of the position of the saddle point to smaller X values (see Fig. 3b) with increasing N . Configurations obtained just after the transition, at h X = 0.06, are plotted in Fig. 6b as both local χ and σ maps. While pleated regions of higher local stress similar to the SUS-MC results are also seen here, the arrangement of the pleats is disordered. Close examination of the configurations also suggests that some of the pleated regions are amorphous. This may be understood as follows. As soon as the thermal energy required to cross the free energy barrier is available, the solid begins to form local pleats. Since many equivalent pleated states are equally stable at these high values of h X deep within the equilibrium phase boundary, the solid locally chooses between the several degenerate pleated states and relaxes, typically, to the nearest metastable free energy minimum. Further relaxation to the true equilibrium ground state, however, now needs large scale rearrangements of the network. As a consequence, the pleated solid shows ageing dynamics.
To show this, we compute the overlap function Q(t) = N −1 N i w(|r i (t) − r i (t w )|) for all vertices with local χ ≥ χ cut = 1. The weight function w(x) is zero or 1 depending on whether x > a or x ≤ a, with a being some predetermined length, smaller than the lattice spacing l; we choose a = 0.1l. For a fixed value of h X = 0.06, starting from an initial crystalline structure, the system is allowed to relax for a "waiting" time t w , before Q(t) is computed. To obtain good statistics, Q(t) is averaged over many independent runs.
In Fig. 7a, we plot Q(t) for five values of t w spanning three orders of magnitude. In each case Q(t) shows an initial rapid decrease to a plateau value, then slow relaxation in the plateau followed by an eventual escape away from the plateau at large times. Q(t) depends on both t and the waiting time t w implying that the system shows ageing behaviour. For a process where the system relaxes quickly to a steady state configuration, Q(t) at long times always has a non-zero limiting value. To illustrate this, we plot in Fig. 7a Q(t) for all the vertices in the crystalline network at a lower h X = 0.02. The long time behaviour of Q(t) in this case is very different, saturating to a constant value. On the other hand, for a system showing complex dynamics requiring long range (and time consuming) particle rearrangements Q(t) relaxes to a plateau first but then eventually decays to zero, on a timescale that grows with t w . To distinguish between these behaviours, we plot Q(t) against the scaled time t * = (t − t w )/t w in Fig. 7b. The data collected over all three decades of t w collapse on a single curve which decays to zero at large values of t * showing that the phase transition kinetics bears strong resemblance to relaxation in a complex landscape.
All the results described so far correspond to the pleating transition in a non-self avoiding, 2d triangular crystalline network. If one thinks of the vertices as colloidal particles, then one might construct such a system experimentally and apply the non-affine field h X using dynamic laser traps in the fashion described in detail in Ref. [11] Of course, colloidal particles are self-avoiding and will have an excluded volume. This should alter the properties of the pleats. Does it also suppress the pleating transition completely? We now show that while the detailed configuration of the pleats are affected because particles cannot overlap, this does not change any of the equilibrium or dynamic results substantially. While a complete overlap is impossible, particles occupy positions determined by a compromise between the bonding and the non-bonding, hard-core repulsion.
We summarise these results in Fig.8 a.-c. The vertices of the network now have, in addition to the harmonic bonded interaction, a non-bonding interaction, which we modelled using the purely repulsive WCA potential in (2).
In Fig. 8a, we plot the probability distribution P (X) for a single value of h X = 0.054, close to but lower than the value at the equilibrium transition. The nature of this curve is very similar to its counterpart for the non-self avoiding lattice (see Fig. 3). However, the pleated states are now somewhat destabilised with respect to the homogeneous lattice. This is only to be expected because the hard core repulsion that prevents particle overlaps now makes pleated states with no bond stretching impossible. This is clear from plots of configurations in Fig. 8b where particles are seen to come close to one another without overlaps. In general, then, the position of the equilibrium transition is shifted to higher values of h X . On the other hand, the position of the saddle point X saddle is virtually unchanged. As a consequence, the location of the dynamical transition does not shift too much so that the equilibrium and dynamical transition points are now closer to each other.
The qualitative nature of the dynamics is similar and the dynamical transition for the same system size occurs again for slightly larger values of h X ≈ 0.06. In Fig. 8 we plot the overlap function Q(t) for this system at a value of h X after the transition. Similar ageing dynamics as seen in the network with point vertices is obtained.

IV. DISCUSSION AND CONCLUSIONS
In this paper we study a crystalline network consisting of a lattice of particles (with and without hard core repulsion) permanently connected to their nearest neighbours by harmonic springs. We show that this undergoes a phase transition from a homogeneous to a pleated phase provided that non-affine displacements are artificially enhanced using an external field. It is interesting to note that the pleating transition in the network, where a homogeneous phase with a uniform stress distribution gives rise to an inhomogeneous phase consisting of an ordered arrangement of pleats where stress is concentrated, has an analogy in the physics of Type II superconductors [26]. Here stress plays the rôle of the magnetic field and the elimination of stress from the un-pleated parts of the network is a manifestation of a "stress Meissner effect". The possibility of such an effect had been described in the past for crystals which have irregular modulations [27]. In these crystals, the order parameter is modulated with a space dependent amplitude and phase producing a disordered structure. The non-affine modes discussed in Section IIIA (see Fig. 1) do produce a similar modulation. However, remarkably, the product state is ordered, contrary to the predictions of Ref. [27]. Strong correlation effects between the order parameter modulations [11] may be responsible for this departure. Nevertheless, this ordered pleated state should be viewed as a classical version of the Abrikosov vortex lattice. The arrangement of pleats now performs the same function as the vortices.
There are several ways in which the calculations described here may be extended. Firstly, we believe that external stress would have significant effects on the transition observed ere. Specifically, compressive stress, should decrease the value of h coex X , perhaps even to h coex X = 0, where the equilibrium transition occurs before the network becomes locally unstable [16]. A similar decrease of h coex X is possible for uniaxial or shear stresses. Such stresses should also introduce anisotropy making it possible to design specific pleating morphologies. Pleating of networks may also have some implications for plastic deformation in these systems. Preliminary investigations by us do point to such a possibility and these results will be published elsewhere.
We have confined ourselves, in this paper, to 2d models where experiments can be performed to check all our predictions with available technology [11]. As detailed in Ref. [11], a feedback-loop may be set up where local particle configurations in a colloidal crystal may be used to compute non-affine forces which are then administered using laser traps positioned on-the-fly [13]. We expect such experiments to yield metastable structures. The enumeration of equilibrium pleated configurations together with their relative free energies and individual barriers and transition states as obtained using techniques elaborated in this paper would, we believe, be useful to analyse the results of these future experiments.
In principle all our calculations can also be extended in a straightforward fashion to higher dimensions as pointed out in [10], though experiments on colloids then become more difficult. Irrespective of this, the analogs of the pleated phase in higher dimensions should certainly be interesting.