Quantifying reversibility in a phase-separating lattice gas: an analogy with self-assembly. James Grant and Robert L. Jack ∗ Department of Physics, University of Bath, Bath, BA2 7AY We present dynamic measurements of a lattice gas during phase separation, which we use as an analogyforself-assemblyofequilibriumorderedstructures. Weusetwoapproachestoquantifythe degree of ‘reversibility’ of this process: firstly, we count events in which bonds are made and bro- ken;secondly,weusecorrelation-responsemeasurementsandfluctuation-dissipationratiostoprobe reversibility during different time intervals. We show how correlation and response functions can be related directly to microscopic (ir)reversibility and we discuss time-dependence and observable- dependence of these measurements, including the role of fast and slow degrees of freedom during 2 assembly. 1 0 PACSnumbers: 2 n I. INTRODUCTION are typically reversible on short time scales, even though a they change macroscopically on long timescales. J Recentwork[19,20,22,23]hassoughttoquantifythis 3 Self-assembly [1–3] is the spontaneous formation of qualitative idea by making dynamic measurements of re- 1 complex ordered equilibrium structures from simpler versibility. We will compare two approaches in a simple component particles. The range of possible structures ] lattice gas model where particles assemble into clusters t includes novel crystals [4–7], viral capsids [8, 9], and col- f and the effectiveness of self-assembly is identified by the o loidal molecules [10–12]. Self-assembly is viewed as a presence (or absence) of under-coordinated particles in s promising alternative technology to current fabrication the clusters. Firstly, we count individual bonding and . at techniques, offering a bottom-up approach to design and unbonding events and compare their relative frequencies manufacture [13]. For example, experimental work has m by defining dynamical flux and traffic observables [23]. usedDNAtogenerateabackboneforpotentialnanofab- Then we use measurements of responses to an applied - rication [14, 15]; potential applications of artificial viral d fieldtoprobe‘reversibility’,usingfluctuation-dissipation n capsidsincludeinertvaccinesanddrugdelivery[16]; and theory [24–26] combined with measurements of out-of- o the potential range of self-assembled structures available equilibrium correlation functions. We relate measure- c through control of particle shape and interactions have ments of fluctuation dissipation ratios (FDRs) [26] to [ been discussed extensively [2, 9, 10]. Self-assembly has the reversibility of self-assembly, extending the analysis 2 even been considered for the potential regeneration of of [19, 22]. In particular, we derive a formula for the v human organs and tissues [17]. However a key challenge response that elucidates its relation to microscopic re- 8 for design and control of self-assembly processes is that versibility and to measurements of flux and traffic [23]. 6 even in systems where the equilibrium states are known, We discuss how these results affect the idea [19, 22] that 0 the conditions which lead to effective assembly remain 6 measurements of reversibility on short time scales might poorlyunderstoodandthereremainsnogeneraltheoret- . be used to predict long-time behaviour. 0 ical approach to support experimental advances. 1 1 A major step towards such a theoretical approach has 1 been the recognition of the importance of reversibility II. THE MODEL : in the assembly process [9, 18–21]. As a system evolves v fromaninitiallydisorderedstateitmakesbondsbetween Xi the constituent particles. If structures form which are A. Definition r not typical of equilibrium configurations then they must a anneal before the system arrives at equilibrium. When We use an Ising lattice gas as a simple model system bonds between particles are too strong, thermal fluctu- for self-assembly. While this model appears much sim- ations are insufficient to break incorrect bonds before pler than more detailed models of crystallisation [4, 7] additional particles aggregate, and the system becomes or viral capsid self-assembly [9, 20], previous work has trapped in long-lived kinetically frustrated states. To shown that it mirrors many of the physical phenomena avoid this problem and achieve effective self-assembly, thatoccurinmodelself-assemblingsystems,particularly bonds must be both made and broken as the system kinetic trapping [21, 23, 27]. evolves in time. In this sense, self-assembling systems At low temperatures, the equilibrium state of the lat- ticegasconsistsoflargedenseclustersofparticles,butef- ficientassemblyoftheseclustersrequiresreversiblebond- inginordertoavoidkinetictrapping(particularlyaggre- ∗Electronicaddress: [email protected] gation into ramified fractal structures [28]). 2 ulations begin from a random arrangement of particles, a) b) 0.75 1 eqbm and we simulate dynamics at fixed bond strength (cid:15) /T. b FLUID binodal The dynamics are based on the ‘cleaving’ algorithm [30] 0.75 0.5 t = 106 which makes cluster moves in accordance with detailed T/!b VAPOUR + LIQUID n4(t) 0.5 balance and ensures physically realistic diffusion. We propose clusters of particles to be moved by picking a 0.25 0.25 t = 104 seed particle at random: each neighbour of that particle isaddedtotheclusterwithprobability1 exp( λ(cid:15) /T) b 0 0 where λ is a parameter that determines t−he rela−tive like- 0 0.25 0.5 0.75 1 0 0.2 0.4 0.6 0.8 " T/!b lihood of cluster rearrangement and cluster motion. We takeλ=0.9(seebelow). Theclustertobemovedisbuilt FIG. 1: a) Exact phase diagram for the lattice gas [29]. At up by recursively adding neighbours of those particles in high temperature the system equilibrates in a single fluid the cluster: each possible bond is tested once. phase while below the binodal it separates into high and low In addition, a maximum cluster size n is chosen max densityphases. b)Plotofyieldn4(t)againstreducedtemper- for each proposed move, and the move may be accepted atureT/(cid:15) ,fort=104 and106 MCS.Forthisrangeoftimes, b only if the size of the cluster to be moved is less than theyieldismaximisedatT/(cid:15) ≈0.35. Theequilibriumyield b n . We choose n to be a real number, greater than is shown (labelled ‘eqbm’): as t → ∞ the yield approaches max max unity, distributed as P(n > n) = 1/n2. The result thisresult. ThebinodalatT/(cid:15) =0.547isshownasavertical max b line. is that clusters of n particles have a diffusion constant D 1/n, consistent with Brownian dynamics: see [30] ∝ for full details. Themodelconsistsofalattice, eachsiteofwhichmay Having generated a cluster, one of the four lattice di- containatmostoneparticle,andtheenergyofthesystem rections is chosen at random, and we attempt to move is the cluster a single lattice spacing in that direction, re- (cid:15) jecting any moves which cause particles to overlap. Fi- E = b np (1) nally, the proposed cluster move is accepted with prob- − 2 p ability P , which depends on its energy change ∆E. (cid:88) a For ∆E = 0, we take a Metropolis formula P = wherethesummationisoverallparticles, (cid:15)b istheinter- min 1,exp(cid:54)[(λ 1)∆E/T]) while for ∆E = 0 weatake action strength, and np the number of occupied nearest P ={α with α−=0.9 a con}stant (see below). (Note that neighbours of particle p. We consider N = 1638 parti- a the definition of the energy in this model means that cles on a square lattice of dimension d = 2 and linear ∆E/(cid:15) is an integer, which ensures that P is monotonic size L = 128. The relevant system parameters are the b a in ∆E for all values of (cid:15) that we consider: cases where density ρ=N/Ld 0.1 and the bond strength (cid:15) /T (or b equivalently, the re≈duced temperature T/(cid:15) ). Thbe phase ∆E/(cid:15)b isnon-integerarediscussedinSec.IIIC.)Simula- b tion time, t, is measured in Monte Carlo sweeps (MCS), diagram is shown in Fig. 1, indicating the single fluid with 1MCS corresponding to N attempted MC moves; phase and the region of liquid-gas coexistence which oc- the time associated with a single attempted MC move is curs when sinh4((cid:15) /2T) > (1 (2ρ 1)8) 1[29]. At the b − − − therefore δt = 1/N. We note that the constants α and density ρ=0.1 at which we work, the binodal is located λ within this algorithm may be chosen freely between 0 at T/(cid:15) =0.547, and below this temperature the system b and 1, while still preserving detailed balance. Our MC separates into high and low density phases at equilib- algorithm mimics diffusive cluster motion most closely rium (the critical point for the model is at ρ = 0.5, c when α and λ are close to unity. However, when evalu- T /(cid:15) = 0.567). For the analogy with self-assembly that c b ating response functions (see below), the MC transition we are considering, we are interested in behaviour at ratesshouldbedifferentiablefunctionsofanyperturbing fairly low densities, as in [9, 10, 20–22]. All of the data fields: this may not be the case if either α and λ is equal that we show is taken at ρ = 0.1 but the same qual- to unity. We therefore take α=0.9 and λ=0.9. itative behaviour is found for lower densities too. (At higher densities, we find that clusters of particles start to percolate through the system. We do not consider this regime since it would corresponds to gelation in the B. Effectiveness of self-assembly self-assembling systems, and this is not relevant for the systems we have in mind.) In studying the lattice gas we have in mind a sys- We use a Monte Carlo (MC) scheme to simulate the temself-assembling(orphaseseparating)intoanordered diffusive motion of particles as they assemble. Our cen- phase such as a crystal, coexisting with a dilute fluid. tralassumptionsare(i)thatclustersofnparticlesdiffuse In the ordered phase, particles have a specific local co- withratesproportionalto1/n,and(ii)thatbond-making ordination environment. In the lattice gas model, we is diffusion-limited (that is, bond-making rates depend take this local environment to be that a particle has all weakly on the bond strength (cid:15) /T while bond-breaking neighbouring sites occupied. We define the ‘assembly b rates have an Arrhenius dependence on (cid:15) /T.) Our sim- yield’,n (t),astheproportionofparticleswiththemax- b 4 3 clusters. We include in Fig. 1 the equilibrium behaviour Temperature for systems of this size, which we obtain by running dy- T/!b=0.15 T/!b=0.35 T/!b=0.5 namicalsimulationsstartingfromafullyphase-separated stateinwhichasinglelargeclustercontainsallparticles. As t , the yield must approach its equilibrium re- t = 104 → ∞ sult: the laws of thermodynamics state that if we wait long enough then the highest quality assembly (or the Time lowest energy final state) will be at the lowest tempera- ture. t = 106 III. MEASURING REVERSIBILITY A. Flux-Traffic Ratio FIG.2: Snapshotsofconfigurationsatthreerepresentative We now turn to measurements of reversibility and temperatures. At (cid:15) /T = 0.15 (low temperature), long-lived b fractal clusters survive to form a stable gel-like structure, their relation to kinetic trapping effects and the non- which we identify as kinetically frustrated. At (cid:15) /T = 0.5 monotonic yield shown in Fig. 1. We follow [23] in con- b (a higher temperature) and for long times, the system fully sidering the net rate of energy changing events at a mi- phaseseparatesintoalargeclustersurroundedbyadilutegas. croscopiclevel,theflux,inproportiontothetotalrateof However,thelargeclusterisnotclose-packedandalargepro- energychangingevents, thetraffic(seealso[31–33]). We portion of particles remain in small clusters: this is therefore countanevent(or‘kink’)wheneveraparticlechangesits poor assembly. At (cid:15)b/T = 0.35 (intermediate temperature) number of neighbours: the number of times that parti- the majority of the particles are in large clusters having few cle i increases its number of bonds between times t and defects, which we identify as good assembly. t+∆t is K+(t,∆t), with K (t,∆t) the number of times p p− the particle decreases its number of bonds. We then av- erage and normalise by the time interval imal number of 4 bonds to neighbouring particles. In 1 Fig. 1b we plot the average yield at two fixed times. We k (t,∆t)= K (t,∆t) . (2) observe a maximum at a reduced temperature that we ± ∆t(cid:104) p± (cid:105) denote by T∗/(cid:15)b: the value of T∗ depends only weakly Hereandthroughout,averages ... runoverthestochas- on time t. The system shows three qualitative kinds of tic dynamics of the system an(cid:104)d o(cid:105)ver a distribution of behaviour. To illustrate this, we show snapshots from initial configurations where particle positions are chosen our simulations in Fig. 2, at three representative bond independently at random. In the limit ∆t 0, then strengths. For strong bonds, there is a decrease in yield the k (t,∆t) converge to the rates for bon→d-breaking ± as the system becomes trapped in kinetically frustrated (+) and bond-making ( ) events, which we denote by states. For weak bonds, the system is already close to k (t). The flux, f(t) =−k+(t) k (t) is the net rate of ± − equilibrium at t = 106 but the final state has a rela- bond-makingandthetrafficis−thetotalrateofallevents tivelylowvalueofn4(t). Inthelanguageofself-assembly, τ(t)=k+(t)+k−(t). (In contrast to [23], we focus here we interpret this as a poor quality product. Based on on rates for making and breaking bonds and not on the Figs. 1 and 2, we loosely identify the kinetically trapped total numbers of bonds made or broken.) regime as T/(cid:15) < 0.2 and the regime of weak bond- We also define a flux-traffic ratio b ing and poor assembly as T/(cid:15) 0.5. However, we b ≥ f(t) emphasise that these regimes are separated by smooth Q(t)= (3) crossoversandnotsharptransitions. Inthefollowing,we τ(t) use T/(cid:15)b = 0.15,0.35,0.5 as representative state points which provides a dimensionless measure of the instanta- for kinetic trapping, good assembly and poor assembly, neous reversibility of a system, consistent with the qual- respectively. itativedescriptionofreversibilityduetoWhitesides[18]. Thus, despite the simplicity of the lattice gas, it cap- The inverse of the flux-traffic ratio 1/Q(t) is the number tures the kinetic trapping effects and the non-monotonic of energy changing events per net bonding event. In a yield observed in more physically realistic model sys- system at equilibrium f(t) = 0 and so Q(t) = 0 also. tems [9, 10, 19–23]. We emphasise that the kinetically For a system quenched to T =0, our MC dynamics does trapped states we find are closely related to diffusion- notpermitbondstobebrokensoparticlesneverincrease limited aggregates [28], while the ‘assembled’ states are their energy: hence k (t) = 0 and Q(t) = 1. In systems − compact clusters. The changes in cluster morphology on that have been quenched, we expect to see a value be- varyingbondstrengtharediscussedfurtherin[27]butin tween these two limits, 0 < Q(t) < 1: the smaller the thisarticleweconcentrateonthedynamicalreversibility flux-traffic ratio the ‘more reversible’ the system, while of assembly [18–20, 22] and not on the structure of the large flux-traffic ratios permit more rapid assembly. 3 Temperature Inthelimit∆t→0,thenthek±(t,∆t)convergetothe ratesforbond-breaking(+)andbond-making( )events, T=0.15 T=0.35 T=0.5 whichwedenotebyk±(t). Theflux,f(t)=k+(−t)−k−(t) isthenetrateofbond-makingandthetrafficisthetotal t = 104 rate of all events τ(t) = k+(t)+k−(t). (In contrast to [23], we focus here on rates for making and breaking bonds and not on the total numbers of bonds made or Time broken.) Wealsodefineaflux-trafficratio t = 106 f(t) Q(t)= (3) τ(t) whichprovidesadimensionlessmeasureoftheinstanta- neousreversibilityofasystem,consistentwiththequal- FIG.2: RLJ:putthetimesandtemperatureetconthefig- itativedescriptionofreversibilityduetoWhitesides[18]. ure. Then we can shorten the caption Snapshots of config- urations at t = 104(top) and 106(bottom) at temperatures Theinverseoftheflux-trafficratio1/Q(t)isthenumber increasingfromlefttoright,T =0.15(kineticallytrapped), of energy changing events per net bonding event. In a T =0.35(goodassembly)andT =0.5(poorassembly). At system at equilibrium f(t) = 0 and so Q(t) = 0 also. low temperature long-lived leggy clusters survive to form a ForasystemquenchedtoT =0,ourMCdynamicsdoes singlegel-likestructureatlatertimes,identifiedaskinetically notpermitbondstobebrokensoparticlesneverincrease frustrated. Athightemperaturethesystemfullyphasesepa- theirenergy: hencek−(t)=0andQ(t)=1. Insystems ratesbythelatertimebutthelargeclusterhasmanydefects that have been quenched, we expect to see a value be- andalargeproportionofparticlesremaininsmallclusters, tween these two limits 0 < f(t) < 1: the smaller the poor assembly. At intermediate temperatures the majority flux-traffic ratio the ‘more reversible’ the system while oftheparticlesareinlargeclustershavingfewdefects,good largeflux-trafficratiospermitmorerapidassembly. assembly. In Fig. 3 we plot the flux f(t), traffic τ(t), and their ratioQ(t)atarangeofbondstrengths,2< /T <10. b obtainbyrunningdynamicalsimulationsstartingfroma In all cases, the flux decreases towards zero as the the 4 fully phase-separated state. As t , the yield must system evolves towards equilibrium. On this logarith- approach its equilibrium result: th→e ∞laws of thermody- mic scale, the temperature-dependence of the flux ap- namicsstatethatifwewaitlongenoughthenthehighest pearsquiteweak,althoughthedifferencebetweendiffer- a) b) c) qatuatlhiteyloawsseesmtbtelym(poerrtahtuerleo.wIenstFeinge.rg2ywfiensahloswtatsen)awpsilhlobtes eRnLtJb:ocnodnssitdreernigntshe1st0mw0iatyhbtiemuepdteopean n4doerndteirnotefgmraatgenditfluudxe?. 100 100 oafssceomnbfilgyu,rkaitnieotniscitnrdapicpaitnivge,aonfdtheefferecgtiivmeesseolff-alosswe-mqubalyl.ity Galwivaeynstlheaststthhean‘in1tth0e-e1grmataexdimflaulxn’uFm( tb2)e=rof0tpdotssfi(btle)b[2o3n]diss 10-1 10-1 (4inthiscase),it10i-s2clearthatf(t)t−1 atlongtimes. T/" T/" RLJ:showthisasdottedline? 0 10-2 b 10-2 b III. MEASURINGREVERSIBILITY Incontrafs(tt)tot1h0e-3fluTx,/t"hbetraffics1how1s0a3tla1r0g6eva!r(iat-) 00.4.55 Q(t) 0.01.51 tionwithbondstrength.0.W5 eakerbondsaremoreeasily 10-3 0.4 10-3 0.2 A. Flux-TrafficRatio broken and resul1t0i-n4 mo00r..e43traffic. Combining flux and 0.03.53 00.2.53 traffic, the ratio1Q0-(5t) is0.2less than 0.1 throughout the 10-4 0.02.52 10-4 00.3.54 good assembly regime, a0n.1d decreases approximately as 0.15 0.45 δtRfoLrJ:1aMmCusminogv∆e.tfWoreannoiwndteufirnnitetopemrieoadsuorfetmimenetsanodf t−1 showthis! wh1i0le-6goodassemblyistakingplace. The 10-5 0.1 10-5 0.5 regimeofkinetictrap10p0ing1i0s1cha1r0a2cte1r0is3ed1b0y4lar1g0e5rval- 100 101 102 103 104 105 100 101 102 103 104 105 reversibilityandtheirrelationtokinetictrappingeffects uesofQ(t)andweakertimedependtence. Similarresults t t and the non-monotonic yield shown in Fig. 1. We fol- werefoundinRef.[23],whereasimilarratiodenotedby low [23] in considering the net rate of energy changing M˜ wasusedtocomparetotalnumbersofbond-breaking events at a microscopic level, the flux, in proportion to thetotalrateofenergychangingevents,thetraffic. We aFndIGb.on3d:-mPalkoitnsgoefvetnhtse. flCuxomfp(atr)edantodMt˜r,affithce τra(tti)o, and their ratio Q(t). The flux varies little across the temperature range countanevent(or‘kink’)wheneveraparticlechangesits Qtc(aotn)ndsdinedopeterneoddnswiotsnhlihyleisotntohryteh:ettrhsatiffisatwceioldlfebpteheenudsseysfsutmlemiunchamtamtkiimonrgee strongly on the temperature. In the kinetic trapping regime T/(cid:15)b (cid:46) 0.2, number of neighbours: the number of times that parti- the ratio Q remains relatively large throughout the trajectories, while the good assembly and high-temperature regimes are contact with other measures of reversibility to be dis- cleiincreasesitsnumberofbondsbetweentimestand associated with reversible evolution and small values of Q. t+∆tisKp+(t,∆t),withKp−(t,∆t)thenumberoftimes cduusrsinedgbopeltoiwm.alTahsesermatbioly,1/inQd(itc)atisinign(tahseirnan[2g3e])10th−at10t0h0e theparticledecreasesitsnumberofbonds. Wethenav- system makes many steps forwards and backwards be- erageandnormalisebythetimeinterval foreIintaFchigie.ve3sawseingplleosttetphoefflneutxprfo(gtre)s,sttroawffiarcdsτt(hte), and their controlled, to avoid kinetic trapping effects. While flux 1 arssaetmiobleQd(stt)ataet. aInrtahnisgesenosfe,bgoonoddsatsrseemnbgltyhrse,q2uir<es(cid:15) /T < 10. and traffic observables are not readily measured except k±(t,∆t)= ∆tKp±(t,∆t). (2) siIgnniaficllanctasreevse,rtsihbeilifltyu,xasdaercgureedasaebsovteo.wardszeroastbhesystem in simple simulation models, we now show how correla- evolves towards equilibrium. On this logarithmic scale, tion and response functions can be used to reveal simi- the temperature-dependence of the flux appears quite lar information (see also [19, 22]). These functions can weak, although the difference between different bond be measured without the requirement to identifyspecific strengths may be up to an order of magnitude. Given bondingandunbondingevents;insomecasescorrelations that the ‘integrated flux’ (t) = tdtf(t) [23](Fig.3a andresponsescanalsobecalculatedexperimentally[34– inset) is always less than tFhe maxim0al (cid:48)num(cid:48)ber of possi- 36]. ble bonds (4 in this case), it is cle(cid:82)ar that f(t) (cid:46) t−1 at At equilibrium, fluctuation dissipation theorems long times. (FDTs) [37] allow the response of a system to an exter- In contrast to the flux, the traffic shows a large varia- nal perturbation to be calculated from correlation func- tion with bond strength. Weaker bonds are more easily tionsmeasuredintheabsenceoftheperturbation. Away broken and result in more traffic. Combining flux and from equilibrium, such comparisons may be used for a traffic, the ratio Q(t) is less than 0.1 throughout the classification of aging phenomena [26] and may also be good assembly regime, and decreases approximately as usefulforestimatingthedegreeof(ir)reversibilityinself- t 1 while good assembly is taking place. The regime of assembling systems [19, 22]. − kinetic trapping is characterised by larger values of Q(t) We first consider the general case for measurement of and weaker time dependence. Similar results were found aresponsefunctioninasystemwithMCdynamics. Fig. inRef.[23],whereasimilarratiodenotedbyM˜ wasused 4 illustrates our notation and the procedure used. The to compare total numbers of bond-breaking and bond- system is initialised at t = 0 and evolves up to time making events. Compared to M˜, the ratio Q(t) depends w, when a perturbing field of strength h is switched on. only on the state of the system at time t and not on its A single MC move is attempted with the field in place, history: this will be useful in making contact with other after which the field is switched off. The system evolves measuresofreversibilitytobediscussedbelow. Theratio withunperturbeddynamicsuptoafinaltime tatwhich 1/Q(t) is in the range 10 1000 during optimal assem- an observable A is measured. In general, the response − bly, indicating (as in [23]) that the system makes many function dependson thetwo times t and w and givesthe steps forwards and backwards before it achieves a single changeintheaveragevalueofA,inresponsetothesmall step of net progress towards the assembled state. In this field h. As indicated in Fig. 4, we use Greek letters to sense, good assembly requires significant reversibility, as represent configurations of the system. (Recall that δt is argued above. the time associated with a single MC move.) The MC scheme may be specified through the proba- bilities P0(ν µ) for transitions from configuration µ B. ‘Flux relation’ between correlation and response to ν in a sin←gle attempted move. However, it is con- functions venient to work with transition rates W0(ν µ) = ← (δt) 1P0(ν µ). We also define ρ (µ) as the proba- − w ← Indesigningandcontrollingself-assemblyprocesses, it bility that the system is in the specific configuration µ would be useful if reversibility could be measured and at time w, and the propagator G0(ν µ) as the prob- t ← 5 conjugate observable to the field h. It is convenient to h(t) define a connected correlation function w w+!t t C(t,w)= δA(t)δB(w) (7) (cid:104) (cid:105) R(t,w) (here and throughout we use the notation δO =O O −(cid:104) (cid:105) µ " # t for all observables O). We also define ∂ S(t,w) C(t,w). (8) ≡ ∂w [The MC dynamics evolves in discrete time steps: the interpretation of the time-derivative is discussed in Ap- pendix A.] At equilibrium, one has the FDT Req(t,w)=Seq(t,w). (9) FIG.4: Procedureformeasuringtheimpulseresponse. The system is initialised in a random configuration. After a wait- Out of equilibrium, we define ing time, w, a perturbation is applied for a single MC move duringwhichthesystemattemptsamovefromconfiguration ∆(t,w) S(t,w) R(t,w) (10) µ to ν. The field is then switched off and the system allowed ≡ − to evolve until time t, at which point its configuration is de- as adeviation betweentheresponse ofthe actualsystem noted by γ. The typical response of the system is indicated, andtheresponseofanequilibratedsystemwiththesame togetherwithsnapshotsofshowinghowparticlesmightmove correlation functions. throughthesystem. (Theclusterthatmovesinthestepwhile Various out-of-equilibrium expressions for response the perturbation is applied is highlighted in red.) functions may be derived: see for example [33, 38, 39] wheredifferentrepresentationsaredefinedandcompared with each other. In Appendix A1, we give a proof of (9) ability that a system initially in state µ will evolve to and we derive the general relation state ν over a specified time period t. The superscripts 0 indicate that no perturbation is being applied to the system: we use a superscript h if the field is applied. [It ∆(t,w)= [A(γ)−(cid:104)A(t)(cid:105)]Gu(γ ←ν)F(ν,µ)Jw0(ν,µ) follows that ρ (µ) = G0(µ κ)ρ (κ) where ρ (κ) γνµ is the probabilwity of initiκal cwondit←ion κ0in the dynam0ical (cid:88) (11) simulations.] (cid:80) where Withthesedefinitions,theaveragevalueoftheobserv- ∂ able A at time t is F(ν,µ)=T lnWh(µ ν) (12) ∂h ← (cid:104)A(t)(cid:105)= A(γ)G0u(γ ←ν)Ph(ν ←µ)ρw(µ), (4) and γµν (cid:88) J0(µ,ν)=W0(ν µ)ρ (µ) W0(µ ν)ρ (ν) (13) where A(γ) is the value of A in configuration γ and we w ← w − ← w introduce u t w δt for compactness of notation. istheprobabilitycurrentbetweenconfigurationsµandν ≡ − − (Since we use Greek letters to represent configurations, attimew. Wewillseebelowthattheprobabilitycurrent there should be no confusion between A(γ) and A(t) , J0 plays a central role in quantifying irreversibility and (cid:104) (cid:105) the former being a property of configuration γ and the deviations from equilibrium. latter a time-dependent average.) The definition of the The formula (11) is general for discrete time Markov (impulse) response is processes obeying detailed balance (the generalisation to Markov jump processes in continuous time is straight- T ∂ A(t) forward). We emphasise that (11) is a representation of R(t,w) (cid:104) (cid:105), (5) ≡ δt ∂h (6), valid whenever detailed balance holds. As such, it is mathematically equivalent to several other representa- where the derivative is evaluated at h=0. Hence, tions that have been derived elsewhere: for example, the ‘asymmetry’functiondiscussedbyLippielloetal.[39,40] ∂ R(t,w)=T A(γ)G0(γ ν) Wh(ν µ)ρ (µ). playsarolesimilarto∆(t,w)inthecaseofLangevinpro- u ← ∂h ← w γνµ cesses. The purpose of (11) is to make contact between (cid:88) (6) reversibility and deviations from FDT in self-assembling We assume that the system obeys detailed balance systems, as we now discuss. Further details of the rela- with respect to an energy function Eh =E0 hB where tion between (11) and previous analyses [33, 38–42] are − E0 is the energy of the unperturbed system and B is the discussed in Sec. IIIG below. 6 a) T/! =0.15 b) T/! =0.35 c) T/! =0.5 b b b 1 1 1 S˜∆t(t,w ) S˜∆t(t,w ) S˜∆t(t,w ) 0.75 R˜∆t(t,w ) 0.75 R˜∆t(t,w ) 0.75 R˜∆t(t,w ) ∆˜∆t(t,w ) ∆˜∆t(t,w ) ∆˜∆t(t,w ) 0.5 0.5 0.5 0.25 0.25 0.25 0 0 0 1000 1250 1500 1000 1250 1500 1000 1250 1500 t t t FIG. 5: Data for S˜∆t(t,w), R˜∆t(t,w) and ∆˜∆t(t,w) at three representative temperatures. The data are obtained at fixed w=103 MCS,whilethetimetisvaried. ThedeviationbetweencorrelationS(t,w)andresponseR(t,w)issmallinbothgood assembly and poor assembly regimes, but significant in the kinetically frustrated regime. We emphasise that the correlation and response functions are associated with perturbations to particle bond strengths: that is, A=B =np in the definitions of S and R. C. Energy correlation and response functions in Weuseastraightforwardgeneralisationofthe‘no-field’ the lattice gas method [41, 42] to allow efficient measurement of the re- sponse: see Appendix A2 for details. To attain good statistics, we consider responses in which the perturbing We now turn to response functions in the lattice gas. fieldhactsnotjustforoneMCmovebutforatimeinter- We consider how a single particle in the lattice gas re- val∆t=10MCsweeps. Sinceweareworkingatleading sponds if the strength of its bonds are increased. To this order in h, the response to such a perturbation is simply end, we write the energy in the presence of perturbing fieldsh asEh = [1(cid:15) h ]n wherethesumrunsover p p 2 b− p p (∆t/δt) 1 all particles, as in (1). We measure the response of n δt − (cid:80) (cid:104) p(cid:105) R∆t(t,w) R(t,w+jδt). (14) to the field h , by taking the observables A and B of the ≡ ∆t p j=0 previoussectionbothequaltothenumberofbondsn for (cid:88) p a specific particle p. Thus R(t,w) = T ∂(cid:104)np(cid:105). Given the The relevant correlation function in this case is δt ∂hp C(t,w)= δn (t)δn (w) where δn (t)=n (t) n (t) energyfunctionEh,thereisstillconsiderablefreedomto (cid:104) p p (cid:105) p p −(cid:104) p (cid:105) asusual. Itisconvenienttodefinenormalisedcorrelation choose the MC rates Wh while preserving detailed bal- and response functions ance. In [33, 39, 40], it was mathematically convenient to 1 take Wh(ν µ) = W0(ν µ)eh[B(µ) B(ν)]/2. Here S˜∆t(t,w)= [C(t,w+∆t) C(t,w)] (15) − (w) − ← ← ∆t our dynamics we are motivated by the central assump- N 1 tions (i) and (ii) of Sec. IIA, which indicate that rates R˜ (t,w)= R (t,w) (16) ∆t ∆t (w) forbond-makingandclusterdiffusionshoulddependonly N∆t weakly on perturbations to the particle bond strengths. with (w) = C(w +∆t,w +∆t) C(w +∆t,w) so We therefore include all h-dependence in the probabil- N∆t − that S˜ (w+∆t,w)=1. We also define ity Ph of accepting an MC move, as follows. If the ∆t a change in the unperturbed energy for an MC move is ∆˜ (t,w)=S˜ (t,w) R˜ (t,w). (17) ∆E0 and the change in the perturbation is ∆V = ∆t ∆t − ∆t h ∆n then for ∆E0 = 0 we take take Ph = − p p p (cid:54) a In all cases, the subscript ∆t indicates that these func- min(1,e(λ 1)∆E0/T ∆V/T) while for ∆E0 = 0 we take tions generalise the R(t,w), S(t,w) and ∆(t,w) of the (cid:80) − − Ph = 2α/(1 + e∆V/T). It is easily verified that this previous section, while the tilde indicates normalisation. a choice is compatible with detailed balance and reduces For small ∆t, the dependence of these functions on ∆t is to the unperturbed probabilities P if ∆V =0. Further, weak, but numerical calculation of the response is easier a coupling the field h only to the acceptance probability for larger ∆t. p ensures that the perturbation affects the rates for bond- InFig.5weshownumericalresultsforcorrelationand breaking and bond-making but does not affect the rates response functions. At T/(cid:15) = 0.5, clusters of parti- b for diffusion of whole clusters. cles are growing in the system, which is far from global 7 equilibrium. However, as found in [19, 22] for several ∆˜ (t+∆t,t) a) other self-assembling systems, the deviations from FDT ∆t 1 are small because the time-evolution is nearly reversible T/! b (recall Fig. 3). On the other hand, at T/(cid:15) = 0.15, 0.15 b 0.2 the particles are aggregating in disordered clusters and 0.75 0.25 0.3 the time evolution is far from reversible with unbonding 0.35 events being rare. 0.4 0.5 0.45 The potential utility of this result was discussed 0.5 in [19, 22]: it means that straightforward measurements on short time scales can be used to predict long-time 0.25 assembly yield, by exploiting links between correlation- response measurements and reversibility. In what fol- 0 lows, we explore in more detail how the deviation func- 0 500 1000 1500 t tion ∆ couples to microscopic irreversibility during as- sembly. ∆˜ (t+∆t,t) b) ∆t 1 D. Interpretation of ∆(w+δt,w) as a flux 0.75 T/! b In this section, we concentrate on quantities that can 0.5 be written in the form 0.5 0.45 0.4 0.35 Z(t)= z(ν,µ)J0(ν,µ) (18) 0.3 t 0.25 0.25 µν 0.2 (cid:88) 0.15 where J0 is the probability current, defined in (13). 0 w In contrast to straightforward one-time observables like 0 0.25 0.5 0.75 1 A(t) = A(µ)ρ (µ), we will show that observables Q(t) (cid:104) (cid:105) µ t like Z(t) are currents, time-derivatives and fluxes: they (cid:80) measure deviations from instantaneously reversible be- haviour, measured at time t. FIG. 6: a) The deviation from FDT measured immediately To this end, we generalise the analysis of “kinks” in afterperturbation,∆˜∆t(t+∆t,t),∆t=10. Attemperatures Sec. IIIA by defining Kνµ(t,∆t) as the average num- above T∗, ∆˜∆t(t+∆t,t) decreases quickly with t, while at ber of MC transitions(cid:104)from state(cid:105)µ to state ν between low temperatures it remains significant throughout the early times t and t+∆t. The associated kink rate is kνµ(t)= times. Intheconstructionof∆(t,w)wetaketheobservables (δt) 1 Kνµ(t,δt) . It follows from the definition of the to be A = B = np as discussed in the main text. (b) Para- ks t−ha(cid:104)t (cid:105) metricplotof∆˜∆t(t+∆t,t)againsttheflux-trafficratioQ(t), for0<t<1500. Therelationshipbetweenthesetwoquanti- kνµ(t) kµν(t)=J0(ν,µ) (19) tiesdependsveryweaklyontemperature,indicatingtheclose − t relation between them. Thus, the probability current J0 gives the difference t between the probabilities of forward and reverse MC transitions between states µ and ν, which clarifies its structure of (18) together with this choice of z ensures connection to irreversibility. At equilibrium, one has that f(t) acquires negative contributions from MC tran- J0(ν,µ)=0forallstatesµandν: thus,J0 measuresde- sitionswhereparticlepexperiencesadecreaseinitsnum- t t viations from equilibrium. However, measuring (or even ber of bonds, as required.) Similarly the particle current representing) J0 is a near-impossible task in a system in exclusion processes is obtained by taking z(ν,µ) = 1 t where the number of possible states µ is exponentially for transitions µ ν where a particle hops to the right, → large in the system size. and zero for all other pairs. Clearly, these observables Instead, one makes a specific choice of the z(ν,µ), in monitor the breaking of time-reversal symmetry, as is which case (18) shows that Z(t) is a projection of the the case for all observables Z(t). current J0 onto the observables of the system. For ex- To relate deviations from FDT to irreversible events, t ample, if z(ν,µ) = A(ν) then Z(t) = ∂ A(t) is a time we first consider the case t = w +δt, so that the per- ∂t(cid:104) (cid:105) derivative that clearly vanishes at equilibrium (or in any turbation h is applied for just one MC move and then steady state). Other choices for z(ν,µ) become relevant the response is measured immediately. In this case (11) inout-of-equilibriumsettings. Forexample,thefluxf(t) reduces to inSec.IIIAisobtainedbysettingz(ν,µ)=1ifthetran- sitionµ ν involvesanincreaseinthenumberofbonds ∆(w+δt,w)= z (ν,µ)J0(ν,µ) (20) → ∆ w of particle p, with z(ν,µ) = 0 for all other µ,ν. (The ν,µ (cid:88) 8 as in (18), with χ˜(t,w) a) z (ν,µ)=[A(ν) A(w+δt) ]F(ν,µ). (21) 1 ∆ −(cid:104) (cid:105) FDT T/!b 0.5 The function F(ν,µ) was defined in (12) and measures 0.8 0.45 0.4 the effect of the perturbation h on the transition rate 0.35 from ν to µ. The factor A(ν) A(w +δt) indicates 0.6 0.3 −(cid:104) (cid:105) 0.25 whether configuration ν has a high or a low value of the 0.2 observable A, compared to the average of A at the mea- 0.4 0.15 0.1 surement time w+δt. For the case considered in Sec. IIIC and Fig. 5, where 0.2 the perturbation is coupled to the energy of particle p, one has A(ν)=n (ν) and 0 p 0 0.2 0.4 0.6 0.8 1 F(ν,µ)=[n (ν) n (µ)]Θ(E0(µ) E0(ν)) (22) C˜(t,w) p p − − whereΘ(x)isthestepfunction,definedsuchthatΘ(0)= FIG. 7: Integrated correlation-response plot showing the 21. The step function appears because if a proposed MC responseχ˜(t,w)andcorrelationC˜(t,w). Thecorrelationand move results in a decrease in the bare energy, the ac- response functions are associated with observables A = B = ceptance probability does not depend on the perturbing np, as in previous figures. The data is shown for a fixed fieldhp. Clearly,F(ν,µ)isfiniteonlyifparticlepchanges t = 104 MCS and and 0 ≤ w < t. The × symbol indicates its energy between states µ and ν: the same is true for thepointswherew=103: impulseresponsesassociatedwith z (µ,ν). this time are shown in Fig. 5. ∆ Thus, (20) indicates that the deviation ∆(t,w) mea- sured at t = w + δt reflects the imbalance between since the system forgets its memory of the configuration rates of bond-making and bond-breaking for particle p. ν and“regressesbacktothemean”. (Thatis,thepropen- In this respect it is similar to the flux f(w): however, sities for different fixed initial conditions all converge to the weights given to different bond-making and bond- the same average value at long times.) Out of equilib- breakingprocessesdifferbetween∆(w+δt,w)andf(w), rium, this may not be the case: if a configuration ν has due to the different forms of z(ν,µ) in the two cases. a more ordered structure than another configuration ν Thus, while these quantities reflect similar physics, the (cid:48) thentheseconfigurationsmayhavedifferentpropensities details of their behaviour is different. To see the rela- for order even at much later times. Loosely, a system in tionship between the immediate response ∆(w,w +δt) ν has a ‘head start’ along the route to assembly, and the and the flux f(t) from Sec. IIIA, we present Fig. 6 system does not forget the memory of this head start as which shows that the two quantities have similar time assembly proceeds. andtemperature-dependenceandthereforerevealsimilar Fig.5showsthatforthecorrelationandresponsefunc- information about the (ir)reversibility of the dynamical tions considered here, the deviation ∆(t,w) decays only bonding and unbonding processes. very slowly with time t. The t-dependence of this func- tion comes through a weighted sum of propensities: it E. Time intervals for reversibility is therefore clear that these propensities do not regress quickly to the mean. In contrast, one may write the cor- relation as We now consider the deviation ∆(t,w) for t>w+δt. From(11),weseethat∆(t,w)canbewrittenintheform S(t,w)= [ A(ν) A(t) ][n (ν) n (µ)] (18) if we take Pu −(cid:104) (cid:105) p − p νµ (cid:88) z(ν,µ)=[PuA(ν)−(cid:104)A(t)(cid:105)]F(ν,µ). (23) ×W0(ν ←µ)ρw(µ) (25) where whose t-dependence also comes from the propensity. A(ν)= A(γ)G0(γ ν) (24) Fig. 5 shows that this function does decay quite quickly Pu u ← witht,althoughitdoesnotreachzeroonthetimescales γ (cid:88) considered here. is the propensity [43] of observable A for configuration ν Our conclusion is the following. For configurations ν after a time u. That is, A(ν) is the value of A that is that are typical at time w, the propensity A(ν) has a Pu Pu obtained by an average over the dynamics of the system, ‘fast’ contribution that decays quickly with time u, as for a fixed initial condition ν and a fixed time u. (Recall well as a ‘slow’ contribution that depends weakly on u u=t w δt.) and reflects the ‘head start’ of ν along the route to the − − Comparing(23)with(21),thedifferenceisthereplace- assembled product. The fast contribution reflects rapid ment of A(ν) by the propensity. At equilibrium, one ex- bond-making and bond-breaking events that do not lead pects A(ν) A(t) todecaytowardszeroasuincreases, directly to assembly while the ‘slow’ contribution is a Pu −(cid:104) (cid:105) 9 non-equilibrium effect that measures assembly progress χ˜(t,w) and also dominates the deviation ∆(t,w) that we have 1 T/! defined here. By contrast S(t,w) picks up contributions FDT b 0.5 from both slow and fast contributions. The utility of 0.8 0.45 0.4 the FDR is that for suitable w, it couples to the irre- 0.35 versiblebond-makingbehaviourthatismostrelevantfor 0.6 0.3 0.25 self-assembly. 0.2 To illustrate this balance of fast and slow degrees of 0.4 0.15 0.1 freedom, and also to make contact with previous stud- ies [19, 22, 26] of FDRs, we define an integrated (and 0.2 normalised) response 0 1 t 0 0.2 0.4 0.6 0.8 1 χ˜(t,w)= dw(cid:48)R(t,w(cid:48)) (26) C˜(t,w) C(t,t) (cid:90)w and we also normalise the correlation as FIG.8: Integratedcorrelation-responseplotasitedependent perturbation. Thatis,wetaketheobservablesA=B =ρi as C(t,w) describedSec.IIIF,andincontrasttopreviousfigureswhere C˜(t,w)= C(t,t) (27) A = B = np. The data is shown for fixed t = 104 with ×symbolsindicatingw=103,asinFig.7. Therelaxationof It is conventional to display these correlations func- theobservableρi isdiffusiveanddependsweaklyonthebond tions in a parametric form [26]. Specifically, making a strength, so the dependence on T/(cid:15)b is much weaker than in parametricplotofχ˜(t,w)againstC˜(t,w),forfixedtand Fig. 7 and the deviations from FDT behaviour are smaller. varying w, the gradient is X(t,w). It has been empha- − sised several times [44–46] that plotting data for fixed w the fluctuations of fast (reversible) degrees of freedom and varying t is not in general equivalent to this proce- regressbacktothemeanwhiletheslow(irreversible)ones dureandmaygivemisleadingresults,especiallyifC(t,t) remain significant. For this reason, the FDR X(t,w) is has significant time-dependence, as it does in these sys- a more sensitive measure of irreversibility than the flux- tems. Plotting data at fixed t as in Fig. 7, it is ap- traffic ratio discussed in Sec. IIIA. parent that the high temperature systems are the most reversible, with X(t,w) 1. To understand the con- ≈ nection with Fig. 5, note that the time interval t w increases from right to left in Fig. 7: when t w is la−rge F. Responses to different perturbations − then the deviation ∆(t,w) has a larger fractional contri- bution to S(t,w), so that the curves are steepest (most Weremarkthattheobservablesusedinmeasuringcor- reversible) when t w is small and least steep (least re- relation and response functions may strongly affect the − versible) when t w is large. We find that parametric results. Mostrelevantistheextenttowhichtheresponse − plots depend weakly on the fixed values of t used but function couples to the slow, irreversible degrees of free- we do not show this data, for brevity: see for example dom and the fast, reversible ones. To illustrate this, we Ref. [22] for an analysis of t-dependence in a system of compare the results of the previous section with a differ- crystallising particles. ent response function. We take the observables A and B Physically, we attribute the w-dependence of X(t,w) of Sec. IIIB both equal to the occupancy ρ = 0,1 of a i to the fact that fast degrees of freedom (like dimer for- specific site i on the lattice. The effect of the perturba- mation and breakage in the vapour phase) quickly relax tionontheMCdynamicsisthesameasthatdescribedin to a quasiequilibrium state [27] where probability cur- Sec.IIIC,exceptthattheperturbedenergyEh =E0+V rents J0 associated with this motion are small. It is with V = h ρ , where the sum runs over sites of the i i i only much slower degrees of freedom (like the gradual lattice. (cid:80) growth/assembly of large clusters) for which probability Unlike the bond numbers n which change only when p currents remain significant on long time scales, and lead bonds are made and broken, the site occupancies ρ also i to long-time contributions to ∆(t,w). change as clusters diffuse through the system. When Compared to the simple flux-traffic ratio Q(t) consid- bonds are strong, the site correlation functions decay eredinSec.IIIA,thenormaliseddeviation1 X(t,w)= faster than the bond correlation functions considered − ∆(t,w)/S(t,w) has the same effect of comparing a gen- above. Therefore, since the site response is coupling to eralisedflux∆(t,w)withanormalisationS(t,w)thatre- fastdegreesoffreedom,weexpecttosee‘morereversible’ flectsthetrafficinthesystem. However, theeffectofthe behaviour: this is borne out by Fig. 8 where the FDRs time difference t w is that while ∆(t,w) decays slowly aremuchclosertotheequilibriumvalue(unity)thanthe − with time t then S(t,w) decays quite quickly. Thus, as FDRs shown in Fig. 7. This is consistent with recent re- t increases, the deviation 1 X(t,w) becomes increas- sults of Russo and Sciortino [46] who observed a similar − ingly sensitive to the deviations from reversibility, since effect, and with earlier studies of observable-dependence 10 ofFDRs[44,45,47](seehowever[48]inwhichsomefast where (w) measures the h-dependence of the amount T degrees of freedom do not relax to quasiequilibrium). of dynamical activity (traffic) between times w and w+ Similar correlation and response measurements in an δt. The second term is therefore associated with traffic Ising lattice gas were also considered by Krząkała [49], while the first term is a correlation between A(t) and but the dynamical scheme used in that work means that the ‘excess entropy production’ at time w: it is therefore clusters of particles do not diffuse: in that case the related to a flux. However, neither of the correlation site degrees of freedom ρ are not able to reach quasi- functions in (28) may be written in the form of (18) so i equilibrium and one does observe significant deviations they do not vanish at equilibrium. The key point is that from FDT for this observable. inequilibrated(reversible)systems,thetermsin(28)are We end this section by briefly considering these corre- both non-zero and equal to each other. So while (28) lation and response functions for large t. In this limit, does involve a separation of terms symmetric and anti- standard arguments [26] indicate that both correlation symmetric under time reversal, the resulting correlation and response functions have contributions from well- functions are not fluxes in the form of (18) and do not separated fast and slow sectors, with the contributions provide the same direct measurement of the deviation fromthefastsectorobeyingFDT.Incoarseningsystems from reversibility given by ∆(t,w) in (11). like this one then there is no response in the slow sec- In another relation to irreversible behaviour, the anal- tor. Thus, for large enough times, the parametric plot ysis of Refs. [39, 40] identifies an ‘asymmetry’ measure- has a segment with X = 1 for large C and small χ; ment which in our notation is for smaller C then X = 0 so the parametric plot has a 1 (t,w)=R(t,w) [S(t,w)+S(w,t)] (29) plateau where the value of χ(t,w) depends only on equi- A − 2 libriumpropertiesofthesystem. However,weemphasise Itfollowsfromthefluctuationdissipationtheorythatthis that the approach to this large-t limit is very slow and quantityvanishesinasystemwithtime-reversalsymme- applies only after all clusters in the system are annealed try. Comparison with (10) also illustrates the similar- into compact clusters. This limit is therefore irrelevant ity between (t,w) and the deviation ∆(t,w). However, for analysing the kinetic trapping phenomena on which A (t,w)differsfrom∆(t,w)sincethedeviation∆(t,w)is we focus in this article. A aprojectionofJ0 andhencevanishesifthesystemisbe- w having reversibly at time w, regardless of what happens at later times. On the other hand, (t,w) vanishes only G. Relation to previous analyses of if all trajectories between times wAand t are equiproba- non-equilibrium response ble with their time-reversed counterparts. That is, time- reversalsymmetrymustholdthroughouttheintervalbe- We have emphasised throughout this article that our tween w and t, and not just at time w. main purpose is to use correlation and response mea- surementstomeasurereversibilityinassemblingsystems. Many other studies have considered such measurements IV. OUTLOOK inavarietyofothercontexts–herewemakeconnections between our methodology and some other results from We have analysed the extent to which ideas of re- the literature. versible bond-making are applicable to a phase separat- One recent area of interest has been the use of ‘no- ing lattice gas, which we used as a simplified model for field’ measurements of response functions. Here, the aim a self-assembling system. The data of Figs. 1, 3 and 7 is to develop formulae for the response that can be eval- further support the idea that this system is relevant for uated in Monte Carlo simulations without the introduc- studying self-assembly, since similar results have been tion of any direct perturbation [39–42, 45]. In fact, we presentedformoredetailedmodelsofself-assemblingsys- use such a straightforward generalisation of the method tems [9, 19, 20, 22, 23]. of [41] to measure responses in this paper, as discussed The flux/traffic measurements of Sec. IIIA provide an in Appendix A2. intuitive picture of reversibility, while the picture based Our analysis in Sec. IIIB is concerned not so much on the correlation-response formalism is more subtle. with measurement of the response, but with the phys- However,thecentralresult(11)demonstratesanexplicit ical interpretation of deviations from FDT. Our result link between measurements of response functions and ir- is therefore in the same spirit as [19, 24–26, 33, 38, 40]. reversibility of bonding. Also, Figs. 5 and 7 do indicate In particular, Baiesi et al [33] relate the FDT to ‘flux’ that information about both short-time reversible and and ‘traffic’ observables that separate reversible and ir- long-time irreversible behaviour can be obtained by con- reversible behaviour, but they use a different approach sidering the behaviour of correlation and response func- to ours. Full details are given in Appendix. A3 but in tions. our notation, their central result is It has been argued that correlation-response measure- mentsonshorttimescalesmightbeusedtopredictlong- 1 ∂ term assembly yield, and even to control assembly pro- R(t,w)= A(t)B(w) + A(t) (w) (28) 2 ∂w(cid:104) (cid:105) (cid:104) T (cid:105) cesses. The work presented here places this objective (cid:20) (cid:21)