Considerations on the relaxation time in shear-driven jamming

We study the jamming transition in a model of elastic particles under shear at zero temperature, with a focus on the relaxation time $\tau_1$. This relaxation time is from two-step simulations where the first step is the ordinary shearing simulation and the second step is the relaxation of the energy after stopping the shearing. $\tau_1$ is determined from the final exponential decay of the energy. Such relaxations are done with many different starting configuration generated by a long shearing simulation in which the shear varible $\gamma$ slowly increases. We study the correlations of both $\tau_1$, determined from the decay, and the pressure, $p_1$, from the starting configurations as a function of the difference in $\gamma$. We find that the correlations of $p_1$ are more long lived than the ones of $\tau_1$ and find that the reason for this is that the individual $\tau_1$ is controlled both by $p_1$ of the starting configuration and a random contribution which depends on the relaxation path length -- the average distance moved by the particles during the relaxation. We further conclude that it is $\gammatau$, determined from the correlations of $\tau_1$, which is the relevant one when the aim is to generate data that may be used for determining the critical exponent that characterizes the jamming transition.


I. INTRODUCTION
In the everyday world there are many examples of disordered systems that change from a flowing state to a rigid state due to a change of some control parameter.The examples are as different as shaving cream and grains in silos.An early suggestion was the existence of the "jamming phase diagram" [1] with the conjecture that the transition from flowing to rigid in disordered systems has the same properties whether the control parameter is density, shear stress, or temperature.It was however later shown that the shear-driven jamming transition, which takes place due to the change in density at zero temperature, is a different phenomenon than equilibrium glassy behavior [2,3].Another fundamental insight is that the shear-driven jamming transition of soft particles is perfectly sharp only in the limit of vanishing shear rate [4].
A path towards a better understanding of these systems is the study of simple models of disorderd collections of particles through computer simulations.A complication is however that the study of zero-temperature processes requires different methods than the ones used in molecular dynamics.In standard molecular dynamics the velocity of the particles automatically makes the system explore phase space and the measurement of different kinds of quantities along this trajectory gives ensemble averages of these quantities.For macroscopic particles of granular matter the thermal velocity is however negligible and ensembles of configurations need to be generated by other means.One way to achieve this is by constantly shearing the system, commonly done with Lees-Edwards boundary conditions [5], which allow for shearing of the system indefinitely with periodic boundaries in all directions.
The shear-driven jamming transition is characterized by the increase of shear viscosity with increasing density and the cleanest behavior is for a system of hard disks.For hard disks the viscosity only depends on the density, φ, and the exponent β describes the divergence of the shear viscosity at φ J , A method for simulating with hard particles has been devised and was used in Ref. [6].That method is however in practice limited to rather small systems and densities at some distance below jamming.This is so since the method relies on the diagonalization of a matrix, which has to be done each time the contact network changes.Another approach is to simulate soft particles at different shear rates γ-where the limit γ → 0 corresponds to the hard disk limit-and to try to determine the behavior in that limit by scaling analyses [4,7].Yet another method to approach the hard disk limit-a method of relevance for the present paper-is to start from configurations produced in the shearing simulations, stop the shearing and let the energy relax down towards zero energy.(Similar relaxations were first done in Ref. [8].)The final part of this relaxation turns out to be exponential to an excellent approximation and each relaxation gives a relaxation time τ 1 .We here use a notation where τ 1 , z 1 , and p 1 denote quantities from single configurations or relaxations, whereas τ , z, and p denote averages over many configurations.From the configurations we also determine z 1 which is the average number of contacts per particle in the final configuration, after the rattlers have been removed.(The rattlers are the particles with ≤ D contacts which do not contribute to the stability of the network.)It is then found that τ 1 is directly related to the distance to isostaticity, given by δz 1 ≡ z c − z 1 [9], with z c = 2d [10].(See Ref. [11] for the generalization of this expression to finite N ).The relaxation time is algebraically related to the distance to isostaticity, τ 1 ∼ δz −β/uz 1 [9].Here z c − z ∼ (φ J − φ) uz and it is commonly believed that u z = 1 [12].
We also remark that it has been claimed [13] that the determinations of τ 1 suffer from a problematic finite size dependence which in effect invalidates the method.As discussed in Appendix A other studies confirm this finite size effect but show that it is a serious problem only for certain cases and that it does not pose a problem for the simulations as in Ref. [9].
The determination of the relaxation time by starting from configurations obtained through steady shearing [9] is therefore a method that gives results relevant for the hard disk limit, that may be used for simple and clean determinations of the critical behavior.The question however remains on how to perform such analyses efficiently and the original motivation behind the present work was to find guidelines for such more efficient simulations.These investigations did however lead to a number of interesting findings both on the workings of the relaxation simulations and for the shear-driven simulations which means that these findings are the central results whereas the quest for an optimal approach in the simulations becomes secondary.
When approaching close to φ J one would better make use of big system sizes both to avoid the spurious jamming that can happen in smaller systems and to get a smaller spread in the obtained τ 1 [9].It is also desirable to perform the initial shearing simulations with small shear strain rates, since that somewhat speeds up the ensuing relaxations.The question is then how to chose the distance in terms of shear strain, γ, between successive starting configurations.The desire is to avoid getting relaxations that are strongly correlated to each other and at the same time avoid wasting simulation time on unnecessarily long shearing simulations between the successive relaxations.
Beside the practical questions there are also questions on the workings of the relaxations and one such question is the connection between the properties of the starting configuration and the relaxed configuration.Specifically we ask to what extent it is possible to predict τ 1 from the pressure p 1 of the starting configuration.The answer to this question is important both for the above stated question on the efficient use of simulations and for a better understanding of the dynamics at densities below jamming.
The organization of the manuscript is as follows: In Sec.II we describe the model and the simulations and also the scaling assumptions and the analyses employed in this work.In Sec.III we start by first examining the correlations of relaxation times τ 1 generated from configurations a distance γ apart.Because of the difficulty in getting good statistics for that quantity we examine to what extent it is possible to predict τ 1 from the pressure of the starting configuration, p 1 , and find that τ 1 is governed both by p 1 and by a random term which is related to the real space distance between the relaxed configuration and its starting configuration.We then also examine the pressure correlations and find a rich behavior, but also show that the behavior may be understood in terms of the scaling approach.In Sec.IV we discuss the implications of the findings for efficient simulations with the aim of high-precision determinations of a critical exponent.In Sec.V we give a short summary of the findings.We also include an appendix which discusses the possibility of logarithmic corrections to scaling that-if present-could make the determination of the critical divergence exceedingly difficult, and has been put forward as the reason for the discrepancy between numerically determined exponents [7,14,15] and the ones obtained from a certain theoretical approach [16,17].

A. Simulation model
For the simulations we follow O'Hern et al. [18] and use a simple model of bi-disperse frictionless disks in two dimensions with equal numbers of particles with two different radii in the ratio 1.4.We use Lees-Edwards boundary conditions [5] to introduce a time-dependent shear strain γ = t γ.We take r ij for the distance between the centers of two particles and d ij for the sum of their radii.The relative overlap then becomes δ ij = 1 − r ij /d ij .The interaction between two overlapping particles is V p (r ij ) = ǫδ 2 ij /2; we take ǫ = 1.The force on particle i from particle j is f el ij = −∇ i V p (r ij ), which means that the magnitude becomes f el ij = ǫδ ij /d ij .The simulations are performed at zero temperature.Length is measured in units of the diameter of the small particles, d s .
The total interaction force on particle i is f el i = j f el ij , where the sum extends over all particles j in contact with i.The simulations have been done with the RD 0 (reservoir dissipation) model [19] with the dissipating force f dis i = −k d v i where v i ≡ dr i /dt − y i γ x is the non-affine velocity, i.e. the velocity with respect to a uniformly shearing velocity field, y i γ x.In the overdamped limit the equation of motion is The equations of motion were integrated with the Heuns method with time step ∆t/τ 0 = 0.2.
Many quantities were measured during the shearing simulations.Two examples of relevance for the present work are pressure p, and the average magnitude of the non-affine particle velocity v = |v i | .We will below refer to earlier analyses of p in Ref. [7].The properties of v have been discussed in Ref. [20]; we here just note that different powers of v scale differently as criticality is approached, such that v 2 i and |v i | 2 diverge differently.

B. Scaling relations
Shear-driven jamming has been found to be a critical phenomenon with shear strain rate one of the relevant parameters, which means that jamming transition takes place at (φ, γ) = (φ J , 0) and that many properties are expected to scale with the distance to jamming.A general review of scaling may be found in Ref. [21].With δφ = φ − φ J the pressure is expected to scale as p(δφ, γ) = b −y/ν g(δφ b 1/ν , γb z ). ( where ν is the correlation length exponent, z is the dynamical critical exponent, and y is the scaling dimension of p.It should be noted that this expression neglects the correction to scaling term [7] which is important for precise analyses close to jamming and has gotten a new interpretation in recent works [22,23].With b = γ−1/z and the notation q = y/zν this becomes [4] To describe the deviations from the hard disk limit we note that the pressure is one of many quantities which is ∝ γ in the γ → 0 limit.For η p ≡ p/ γ Eq. ( 2) translates to and (φ J − φ)b 1/ν = 1 gives b = (φ J − φ) −ν which in turn leads to where β = zν − y.
The present work focuses on the relaxation time τ which, in the hard disk limit, behaves the same as η p [6,14], Another quantity of interest is the velocity per unit of shear strain which diverges as [20] lim with u v ≈ 1.1.For finite γ these expressions may be generalized through the same kind of relations, Note however the different nature of these quantities.The determination of the relaxation time is from two-step simulations that are done by first shearing at a certain shear strain rate and then stop the shearing and let the system relax to zero energy.The non-affine velocity per shear rate, v(φ, γ)/ γ, characterizes the flowing sheared steady state, whereas the relaxation time, τ , is from the final stage of this relaxation.

C. Analyses
To determine the relaxation time we run simulations as described above at zero temperature and fixed γ (i.e. with γ = 0) which leads to an energy decreasing towards zero; the simulations are aborted when the energy per particle is E < 10 −20 .The relaxation time is then determined from the exponential decay of the energy per particle by fitting E(t) to For each parameter set, N , φ, and γ, the starting configurations are from a long shearing simulation that gives a large number of different starting configurations with different γ.The relaxation simulations then give a set of relaxation times, τ 1 (γ).(The individual determinations are denoted by τ 1 whereas their average is denoted by τ .)On general grounds one expects two starting configurations that differ only by a small shear γ to be similar and to also lead to similar relaxation times and one of the goals of the present work is to examine how quickly this similarity decays and thus to determine the correlation shear-which is the analogous of the correlation time-for τ 1 .
The autocorrelation function of some quantity A is a measure of the correlations in the fluctuations of A, i.e. δA = A − A .The autocorrelation function is then δA(t ′ )δA(t ′ + t) , where the average is over different initial times t ′ .When comparing systems that are sheared at different shear strain rates there is a trivial shear strain rate dependence and it is therefore convenient to instead express the correlation function in terms of the change in the shearing variable, γ, We will examine this correlation function for two different quantities, A = 1/τ 1 and A = p and beside some results extracted from these correlation functions we will be able to draw some conclusions about the relaxation processes.
particles and at densities φ = 0.830 through 0.841, closely below φ J ≈ 0.8434.The rationale for determining the correlations of 1/τ 1 instead of τ 1 is that τ 1 can occasionally be very big which could mean that a few big values would dominate the correlation function.This problem is not present when instead determining the correlation of 1/τ 1 .
Fig. 1(a) shows ρ 1/τ (γ).We find that each curve to a decent approximation shows an exponential decay and that the slope of the curves increase as φ increases towards φ J .To determine the correlation shear-the analogous quantity to the correlation time-one should ideally fit the tail of ρ 1/τ (γ) to an exponential decay, but since our data are rather noisy this isn't feasible and we therefore instead determine the correlation shear γ τ from the shear strain that gives ρ 1/τ (γ τ ) = e −1 , i.e. from the crossings of the horizontal dashed line in Fig. 1(a).Fig. 1(b) shows γ τ vs φ and we note that the figure can be taken to suggest a linear behavior, γ τ ∝ φ J − φ.
The data above are determined with N = 4096 and an interesting and non-trivial question is on the finite size dependence of γ τ .One possibility would be that a larger number of particles would give more places with important changes of the contact pattern which would lead to a bigger change in τ 1 for each change in γ.Our results do however suggest that there is no real finite size dependence for our quite big systems.(For sufficiently small N one would however expect changes of all kinds of quantities.)This is illustrated by Fig. 1(c) which shows ρ 1/τ (γ) for N = 4096, 16384 particles, determined with φ = 0.838 and γ = 10 −8 .
Fig. 1(b) is for a single shear strain rate, γ = 10 −8 , and a further question is on the dependence of γ τ on the shear strain rate.To that end we have attempted the same kind of analyses for different shear strain rates but found the data to be somewhat too noisy for any safe conclusions.The basic problem is that the determinations of ρ 1/τ (γ) requires a large number relaxation times, τ 1 , which are obtained through time consuming relaxations of the system to (almost) vanishing energy.With a limited number of such τ 1 the correlation function ρ 1/τ (γ) for different φ become quite noisy which makes it difficult to achieve precise determinations of γ τ .
As an alternative route to more insight we have therefore turned to other analyses.The approach is based on the expectation that the individual τ 1 should be strongly influenced by properties-as e.g. the pressure-of the respective initial configurations, and that the correlations between such quantities should be much easier to determine since there are considerably more available data.For such an initial property we here make use of the pressure, and the approach is therefore to first turn to a comparison between relaxation time and the pressure of the corresponding initial configuration, and then, as a second step, examine these pressure correlations.As we will see these correlations turn out to be somewhat different and do not directly help answer our questions.This approach does nevertheless lead to some unexpected new insights.

B. Relaxation time and pressure
To examine the relation between τ 1 , determined from the last stage of the relaxation, and η p1 ≡ p 1 / γ, from the initial configurations (before the relaxation), Fig. 2 shows τ 1 and η p1 for γ = 10 −8 and N = 4096, at two different densities, φ = 0.836 and 0.841.Fig. 2(a) which is τ 1 vs η p1 for φ = 0.836 shows that these data are strongly correlated whereas Fig. 2(b), obtained at the higher density φ = 0.841, gives evidence for a weaker correlation.Panels shows that these quantities are strongly correlated for the lower density φ = 0.836 and more weakly correlated at φ = 0.841.The lower panels are the same data but now plotted vs γ (∼ simulation time).Panel (c) is for φ = 0.836 whereas panel (d) is for φ = 0.841.In both cases we find that τ1 closely follows ηp1-the minima and maxima of ηp1 are closely followed by similar minima and maxima in τ1-but also that the spread is bigger in panel (d) which is for the higher φ.
(c) and (d) show the same data but now plotted vs γ.For the lower φ in Fig. 2(c) it is clear that the two quantities follow each other closely but Fig. 2(d), which displays the same kind of data for the higher φ = 0.841, shows that the minima and the maxima of η p1 are reflected in τ 1 , but that there is also a big fluctuating random component to τ 1 .
To quantify this fluctuating factor we introduce which is from the ratio of the pressure of the initial configuration and the relaxation rate determine from relaxation simulations.The rational for taking the logarithm is, as shown in Fig. 3(a) and (b), that the distribution of τ 1 /η p1 is strongly skewed whereas the distribution of the logarithm of the same quantity is similar to a Gaussian distribution, as if f 1 were the sum of a number of independent random variables.If τ 1 were exactly predicted by η p1 f 1 would be a constant and to quantity the spread in f 1 we determine which is then a measure of the size of the random fluctuating factor.Fig.The reason for the strong correlation between η p1 and τ 1 at the lower densities is that the each final configuration is very close to its corresponding initial configuration.Conversely, the bigger fluctuations of f 1 at higher densities suggests the that the properties of the system often change a lot during the relaxations.As a measure of the distance between initial and final configurations in ordinary space, we introduce the relaxation path length which is the average distance moved by the particles during the relaxation, The average is here over both particles and relaxation runs performed with the same parameters, φ, γ, and N .
Figure 3(d) is a parametric plot of sdev[f 1 (φ, γ, N )] vs s(φ, γ, N ) for two system sizes, N = 4096, 16384, shear strain rate γ = 10 −8 , and densities φ = 0.834 through 0.8412.The system size dependence found in Fig. 3(c) is here taken care of through a factor of √ N and the plotted quantity is thus an algebraic dependence on s(φ, γ, N ) gives the exponent 0.49 which suggests We note that this is consistent with elementary statistics if one considers f 1 to be the average of N different terms which each is the sum of ∝ s different terms.The logarithm in the definition of f 1 lets us conclude that each small ∆s contributes a random factor to τ 1 .
To summarize this part we have found that τ 1 is directly controlled by η p1 at low densities and that the behavior at higher φ is similar but with big random fluctuations.We have also found that the size of this random contribution to τ 1 is directly related to the average distance moved by the particles during the relaxation.

C. Pressure correlations
We then turn the correlations of pressure with the aim to get an understanding of the size of the shear strain, γ p , that characterizes the decay of the pressure correlations in shear-driven simulations.
Figure 4(a) shows the correlation function ρ p (γ) for N = 4096 particles, shear strain rate γ = 10 −8 , and densities φ = 0.830 through 0.8434.We find that ρ p (γ) decays exponentially for each φ, and determine the correlation shear γ p from the condition ρ p (γ p ) = e −1 .Fig. 4(b) shows γ p vs φ for different γ.The solid dots are from the data in Figure 4(a), the other symbols are for three higher shear strain rates.
We note that the behavior of γ p (φ) as γ → 0 suggest that γ p (φ) in the hard disk limit goes approximately linearly to zero, as shown by the dashed line.This is the same kind of behavior as shown for γ τ in Fig. 1(b) and it is also consistent with the finding that a characteristic shear determined from the velocity-velocity correlation at φ J vanishes as γ → 0 [24].This is also directly related to the fact that v/ γ-the distance traveled per unit shearincreases as jamming is approached.
Another interesting observation from Fig. 4(b) is that this correlation shear, for each finite shear strain rate, γ, depends non-monotonously on φ.For each constant shear strain rate the correlation shear, γ p , first decreases towards a minimum and then increases again when φ approaches φ J from below.We will now first relate this behavior to ideas from scaling, then discuss the physical mechanisms behind this decrease and finally consider the change to an increasing trend of γ p vs φ.
From the scaling assumption in Eq. ( 3) follows that the behavior should be controlled by the combination (φ − φ J )/ γ1/zν .To test this expectation we plot γ p vs (φ − φ J )/ γ1/zν in Fig. 4(c) and we then find that the upward turns, to a decent approximation, take place at a constant (φ J − φ)/ γ1/zν .There is a deviation from that behavior for the smallest shear strain rate, γ = 10 −8 , and we attribute this to a finite size effect which could be visible when the correlation length, which increases as φ J is approached from below [25], becomes comparable to the system size.
Turning to the physical reason for the decrease of γ p with increasing φ, at low φ, we believe that this is an effect of the increasing particle velocity as φ → φ J , described by Eq. ( 7).The reasoning is that we expect two ) The correlation shear, γp, is determined as the γ for which ρp(γ) = e −1 .Panel (b) which is the correlation shear, γp vs φ for three different γ, shows a clear dependence on γ.For each γ there is a minimum at a certain φmin( γ).As γ decreases the position of this minimim moves towards φJ .Panel (c) The same data but plotted vs the scaled distance from jamming as suggested by Eq. ( 3).We note that the minima of the respective curves are on top of each other, except for the data for the lowest γ.This deviation is tentatively attributed to a finite size effect that appears as the system size becomes comparable to the correlation length.
configurations that are generated by the shearing dynamics should be substantially different-such that their respective p 1 are also different-if the particles have on average moved a certain characteristic distance, ℓ.This gives the time scale t v = ℓ/v and γ p = t v γ = ℓ γ/v, which together with Eq. ( 7) for the divergence of v/ γ becomes γ p = (ℓ/A v )(φ J − φ) uv .The dashed line in Fig. 4(b) is an approximate description of the behavior in the γ → 0 limit, assuming u v = 1.The rectilinear behavior shown there is in reasonable agreement with the behavior described by the exponent u v ≈ 1.1.
We now turn to the change in trend of γ p (φ) and we are going to argue that it goes together with a change from hard to soft particles-i.e. a change from negligible to a finite particle overlaps.When the shearing is sufficiently slow that the system has time to relax down to very small particle overlaps the system is in the hard particle limit.Since the energy relaxation is governed by the relaxation time, τ , the system should be close to the hard disk limit as long as the relaxation time is smaller than all other relevant time scales.With the velocity time scale, t v = ℓ/v, introduced above, we get the criterion τ ≪ t v , and the expectation of a change of trend when that condition is no longer fulfilled.
We now argue that the change in behavior when t v ≈ τ is consistent with the expectation from scaling discussed above, that the behavior should depend on the combination (φ − φ J )/ γ1/zν .For this comparison we first make use of the relation zν = β + y [7] to rewrite the scaling variable as From Eq. ( 7) we then find which together with Eq. ( 6) and the criterion τ ≪ t v gives which implies that the criterion for being in the hard disk limit becomes This is similar to Eq. ( 14) and also leads to the suggestion u v = y.We also note that the possibility that these two different exponents actually are the same is in agreement with the numerical values in the literature, y = 1.08±0.03[7] and u v ≈ 1.10 [20].We now also take this one step further and try to approximately express γ p in terms of τ and t v .Since τ is the time scale needed to relax energy or pressure when shearing has been stopped, it is not unreasonable to expect τ to affect the pressure relaxations also in the presence of shearing.However, since the conditions are so different it follows that a possible relation between τ and the pressure correlations could at most be an approximative one, only.
In the following, we will make use of a related but different time scale-the dissipation time, τ diss , described in Attempt to predict γp from our two time scales.The time scales are tv = ℓ/v from the average particle velocity and the dissipation time τ diss , from the initial decay of energy in a relaxation, but determined from properties measured in the shearing simulations and given by Eq. ( 18).Though the agreement with γp in Fig. 4(b) is by no means perfect, we note that this quantity has the same general behavior.
Ref. [14].This time scale is determined from the initial decay of a relaxation in contrast to τ which is determined from the final part of the relaxation.These two time scales-τ and τ diss -behave the same in the γ → 0 limit but have the opposite dependence on γ [14].The decisive advantage with τ diss before τ is however that it may be determined directly from the properties (energy and shear stress) of the shearing simulations [14], without the need for any relaxation steps.Figure 5 shows the simplest possible way to include the effect of these two different time scales, which is to assume that the effective correlation time is obtained by adding together contributions from these two time scales, such that the correlation shear is given by γ p = (t v + τ diss ) γ = (ℓ/v + τ diss ) γ.Here v is the measured average velocity, τ diss is from Eq. ( 18) and the free parameters are ℓ = 0.1, and the prefactor of τ diss which we take to be equal to unity.Though the agreement is by no means perfect we note that the curves in Figure 5 capture the general behavior of the simulation data in Figure 4(b).The conclusion is thus that the non-monotonic behavior of γ p is an effect of the increase of τ as φ approaches φ J from below.

A. Comparing the different correlation functions
In the above sections we examined correlations from two different "time" series, τ 1 (γ) and p 1 (γ) and determined the related γ τ and γ p .Since the determinations of a large numbers of τ 1 are quite computationally demanding the idea was to extract the same kind of information from p 1 of the initial configurations which would thus reduce the need for a large number of relaxation runs.
Fig. 1 suggests that γ τ from the γ dependence of τ 1 vanishes approximately linearly as φ J is approached from below whereas Fig. 4 shows that γ p has minima that depend on γ and move towards φ J as γ decreases.To explain the difference in behavior we then noted that there are two contributions to τ 1 .The first is related to the values of p 1 , and thereby η p1 , in the initial configuration and the second is due to the relaxation processes.The first contribution changes slowly whereas the second component fluctuates much more rapidly and it appears that it is this second component that is responsible for the continued decay of the correlations in Fig. 1(b) when γ p from the pressure correlations in Fig. 4(c) turn upwards.
These conclusions appear to be true when γ τ and γ p are determined as the values of γ that make the correlation functions equal to some given constant, here taken to be e −1 .The correlations that are visible in the pressure correlations ought however to be present also in ρ 1/τ (γ) and should be visible if one were able to access the tail of ρ 1/τ with sufficient precision.To see this we may consider an interval in γ where p 1 is almost a constant.In this range of γ τ 1 will vary around a certain average that depends on the average p.For a different interval in γ with another average p, τ 1 will fluctuate around another average value.From this consideration it appears that correlations that are seen in p 1 should also be present in τ 1 .Since the fluctuations in p 1 are quite small compared to the fluctuations in τ 1 it does however seem that it would be virtually impossible to verify the existence of these correlations in τ 1 .
The question is now what to conclude for efficient relaxation simulations, and what should reasonably be the distance in terms of γ between successive starting configurations.It then appears that the answer-as is often the case-depends on what the obtained data should be used for.If the goal is to get statistically independent values τ 1 for given φ and γ it is reasonable to consider the correlations of the initial configurations and it could then be reasonable to do the relaxation runs by starting at points with the separation γ ≫ γ p .If, on the other hand, the goal is to get a set of points (τ 1 , δz 1 ) for determining the exponent β/u z from τ 1 ∼ (δz 1 ) −β/uz , then a possible bias of these points towards high or small values of τ 1 and δz 1 would not be a problem, and one could then well make use of starting configurations that are quite close together in γ, but still with γ ≫ γ τ where γ τ is from the apparent correlations of 1/τ .

B. Implications for precise determinations of the critical behavior from relaxation simulations
The present results together with earlier findings [14,26] lead to some suggestions for determinations of points (τ 1 , δz 1 ) close to criticality for more precise determinations of the exponent β/u z in τ ∼ (δz) −β/uz : (i) The determinations should best be done with a big number of particles since the spread in the different quantities are ∝ 1/ √ N [14].(ii) It is preferable to do the simulations with small γ since one then has a lower energy to start with, but there is always a trade-off since at large N a small γ could give very long times for generating starting configurations with sufficient big distance in γ. [To be explicit on numbers we note that the simulation of 10 6 time units requires 2 hours when our parallel code is run on 28 cores.Considering simulations at (φ, γ) = (0.842, 10 −9 ) where we have γ τ (0.842) ≈ 0.0045 the simulation to advance γ by 2γ τ (0.842) would then require ≈ 20 hours.](iii) It would be possible to substantially speed up the relaxation simulations by using the fast minimization protocol of Ref. [27] instead of the slow simulations with steepest descent (which e.g. was used in Ref. [14]).The final part of the simulation, which is used for the actual determination of τ 1 , needs however be performed by steepest descent.(iv) It should be noted that the finite precision in the double precision variables that are typically used to store the positions may lead to artifacts for very big systems [23].This is due to two facts.First, that a larger value of a coordinate means that fewer bits are available for storing the fractional part of the position and, second, the fact that the net force on a particle, which determines the dynamics, is often considerably (i.e. a factor of τ ) smaller than the typical contact force.This problem may be taken care of with a code that stores the position in two variables, with fraction part and integer part, which ensures that large values of the position coordinates don't affect their precision.

V. SUMMARY
To summarize we have examined the correlations of both τ 1 and p 1 as a function of γ motivated both by a desire to understand the basic physical mechanisms and to answer the question on how to most efficiently determine statisticaly independent values of the relaxation time, to be used for the determination of a critical exponent.
From τ 1 (γ) and p 1 (γ) we determine the respective correlation shears, γ τ and γ p and our first conclusion is on the behavior in the hard disk limit.For the hard disk limit we conclude that our two different correlation shears both vanish essentially linearly as φ → φ J from below.We note that this is in consistent with the earlier finding that the velocity-velocity correlation at jamming vanishes as γ → 0 [24].
For γ p from the pressure correlations determined for different φ and γ we find that γ p (φ) at constant γ is a non-monotonous function which first decreases with increasing φ, reaches a minimum and then increases again as φ → φ J .We interpret this behavior as an effect of two different time scales where the first is directly related to the average non-affine particle velocity, t v = ℓ/v, and the second is the average relaxation time, τ .Close to the hard disk limit-i.e. at sufficiently low φ-the behavior is dominated by t v but at higher φ the behaviour is instead dominated by τ which diverges as φ → φ J .
Our data for γ τ -the correlation shear for τ 1 -suggests that this quantity decreases monotonously as φ → φ J and we set out to analyze this difference in behavior compared to γ p .For the hard disk limit we expect τ 1 ∝ η p1 [14], and this is also borne out by our data a low densities where the proportionality holds to a good precision for each individual relaxation.At higher densities τ 1 and η p1 do however behave very differently which is seen through τ 1 /η p1 spreading considerably around its average.To quantify this spread we determine the standard deviation in f 1 ≡ ln(τ 1 /η p1 ) and find that sdev[f 1 ] ∼ √ s, where s is the average distance moved by a particle during the relaxation.Due to the logarithm in the definition of f 1 we conclude that that relation implies that the contribution due to each small ∆s is a random factor to τ 1 .The dependence on N is in accordance with elementary statistics of N independent values, which means that this random contribution to τ 1 decreases as the system size decreases.
When it comes to efficient simulations we conclude that the distance between succcessive configurations could reasonably be taken to be γ ≈ 2γ τ (φ), with γ τ ≈ 3.25(φ J − φ), from Fig. 1(b), which means that the necessary distance in γ decreases as φ → φ J and that there is no need for any very extensive shearing simulations between successive starting configurations.such islands for bigger N , and a simple assumption of the distribution of relaxation times, follows the τ ∼ ln N dependence.The further conclusion was that τ is an ill-defined quantity because of this finite-size dependence and that it may therefore not be used to determine the critical exponent.
A later study [26] did however modify these conclusions in several ways.The suggested finite size dependence was confirmed, but it was also shown that the splitting into different islands with different relaxation times only sets in for quite big N , and is not the dominant mechanism for the finite size dependence.The dominant mechanism is instead that big random initial configurations have large density fluctuations that, to some degree, survive the relaxation process and affect the final relaxation.This effect is not present in relaxations starting from configurations obtained at steady shearing, since these starting configurations have a long pre-history of a slow shearing and therefore already have an essentially uniform density.It was also shown that it is possible to define a relaxation time which isn't plagued by the ln N -dependence, and that the problematic finite size effect is not present at the system sizes that have been used in earlier determinations of the critical behavior [9,29], and would not seem to be a problem in possible future attempts to determine the critical behavior with higher precision.

Appendix B: Logarithmic corrections to scaling
The underlying assumption in the above discussion and in the present paper is that it should be possible to determine the critical behavior through analyses of data from shear-driven simulations.An alternative view is that the data are affected by logarithmic corrections that could make such analyses very difficult or perhaps altogether unfeasible.The ground for the suggestion of logarithmic corrections to scaling is twofold: (i) The first reason is that one could expect the presence of logarithmic corrections in systems that are at the upper critical dimension of the model, and since the upper critical dimension of the jamming transition is widely believed to be d ucp = 2 [11,30], such corrections should be expected in two dimensions.(ii) The second reason to consider logarithmic corrections is as an attempt to explain the discrepancy between the values of certain critical exponents from theoretical arguments [16,17] and to the ones from simulations [7,14,15].
We do, however, not consider the evidence for the presence of logarithmic corrections to scaling to be at all compelling.(i) A shear-driven system is in many ways different from static jamming, which means that it could well be that the upper critical dimension is different from two in shear-driven systems even if it were true that it is equal to two for static jamming.(ii) A recent study [23] gives reasons to question one of the key assumptions between the theoretical approach [16,17] which is that the process that governs the divergence of the shear viscosity is "spatially extended".The new conclusion [23] is that the divergence of the viscosity is caused by the fastest particles [20] and that these particles are short range correlated, only [23].A consequence of that reasoning is that the above mentioned discrepancy between theory and simulations could be due to the failure of the recent theoretical treatment [16,17], and that the earlier approach of Ref. [6] might perhaps be correct.
We also comment on the fit of the 2D data to an expression with logarithmic correction to scaling [17] and argue that the fit is not as conclusive as it could seem.From the expectation that arguments to the logarithm function should be dimensionless, we believe that Eq. ( 12) for the relevant eigenvalue in Ref. [17] should rather be written with the additional free parameter δz 0 , λ 1 ∼ (δz) β/uz | ln(δz/δz 0 )| α (β/u z in our notation is their β), where both α and δz 0 are unknown constants.The approach in Ref. [17] amounts to tacitly assuming that δz 0 = 1, without any reason.For τ ∼ λ −1 1 the expression for τ when including corrections to scaling becomes τ ∼ (δz) −β/uz | ln(δz) − ln(δz 0 )| −α , (B1) and it is then evident that the value of δz 0 can not be absorbed in a prefactor.The fact that the value of the exponent, β/u z = 3.41, which fits the data in 3D without corrections to scaling, happens to give a good fit of the 2D data to Eq. (B1) with δz 0 = 1, then appears to be fortuitous only, since there is no reason to believe that δz 0 = 1 actually is the correct value.When taking both α and β/u z to be free parameters one expects to find that acceptable fits are obtained with a wide range of values of β/u z .This should be kept in mind when trying to assess the evidence of Ref. [17].It also seems that ordinary corrections to scaling could give an equally good fit.We finally comment on a possible and reasonable criticism of the determinations of the exponents in Ref. [7], which is that they should be regarded with quite some skepticism since the resulting exponents were obtained through a fitting of data with not too much structure to the sum of two unknown scaling functions-a main term and a correction term-which are both from polynomials with quite a number of fitting parameters.It could then seem that this complicated analysis may well give exponents that are somewhat off the correct values even without possible additional complications as the need for logarithmic corrections to scalinig.This concern is however in effect addressed in Ref. [23] since the origin of the correction term is there identified which makes the correction term known apart from a multiplicative factor.The correction term is there found to be given by const × W p (φ, γ), where W p (φ, γ) is determined from the velocity distribution.This means a great simplification of the analysis and the fact that that approach gave similar values of the critical exponents to the ones in Ref. [7], gives reason to believe that these earlier results actually were correct.

FIG. 1 .
FIG.1.Correlation function for from the inverse relaxation time.Panel (a) is the correlation function ρ 1/τ (γ) obtained with N = 4096 particles, shear strain rate γ = 10 −8 , at densities φ = 0.830 through 0.841.The correlation shear, γτ , is defined to be the value of γ for which the correlation function is ρ 1/τ = e −1 .Panel (b) which is γτ vs φ, shows that γτ decreases as φJ is approached from below in an approximately linear way, γτ ∼ φJ − φ.Panel (c) which shows ρ 1/τ (γ) for two different system sizes, N = 4096 and N = 16384, illustrates that the correlation function is independent of N .

1 FIG. 2 .
FIG.2.Plots of τ1 and ηp1 at two different densities, φ = 0.836 and 0.841.Here ηp1 are from the starting configurations whereas τ1 are from the relaxations.The question is to what extent τ1 from the relaxation may be predicted from ηp1 ≡ p1/ γ of the starting configuration.The data are for N = 4096 particles and shear rate γ = 10 −8 for the shearing simulations.Panels (a) and be which show τ1 vs ηp1 shows that these quantities are strongly correlated for the lower density φ = 0.836 and more weakly correlated at φ = 0.841.The lower panels are the same data but now plotted vs γ (∼ simulation time).Panel (c) is for φ = 0.836 whereas panel (d) is for φ = 0.841.In both cases we find that τ1 closely follows ηp1-the minima and maxima of ηp1 are closely followed by similar minima and maxima in τ1-but also that the spread is bigger in panel (d) which is for the higher φ.

1 P (f 1 )FIG. 3 .
FIG.3.Properties of f1 ≡ ln(τ1/ηp1).All data are for γ = 10 −8 .Panels (a) and (b) show that the distribution of τ1/ηp1 is strongly skewed whereas the distribution of f1 ≡ ln(τ1/ηp1) is similar to a Gaussian distribution.This is taken to suggest that f1 is a reasonable quantity for analyses.Panel (c) shows the spread in f1-denoted by sdev[f1]-vs φ for two different sizes, N = 4096 and 16384.The finite size dependence is to a good approximation given by ∼ 1/ √ N .Panel (d) shows the same data versus s, which is the average displacement during the relaxation, now compensated for the √ N dependence.A fit gives √ N sdev[f1] ∼ s 0.49 , which suggests a dependence ∼ √ s.The conclusion is that each small ∆s from the relaxation contributes a random constant to f1 and thus a random factor to τ1/ηp1.

FIG. 4 .
FIG.4.Correlation shear, γp, determined from p1(γ) measured in the shearing simulation.Panel (a) shows ρp(γ) for N = 4096 particles, shear strain rate γ = 10 −8 , and densities φ = 0.830 through 0.840.(Data for more densities are excluded in order not to clutter the figure.)The correlation shear, γp, is determined as the γ for which ρp(γ) = e −1 .Panel (b) which is the correlation shear, γp vs φ for three different γ, shows a clear dependence on γ.For each γ there is a minimum at a certain φmin( γ).As γ decreases the position of this minimim moves towards φJ .Panel (c) The same data but plotted vs the scaled distance from jamming as suggested by Eq. (3).We note that the minima of the respective curves are on top of each other, except for the data for the lowest γ.This deviation is tentatively attributed to a finite size effect that appears as the system size becomes comparable to the correlation length.
FIG.5.Attempt to predict γp from our two time scales.The time scales are tv = ℓ/v from the average particle velocity and the dissipation time τ diss , from the initial decay of energy in a relaxation, but determined from properties measured in the shearing simulations and given by Eq. (18).Though the agreement with γp in Fig.4(b) is by no means perfect, we note that this quantity has the same general behavior.