Understanding slow compression and decompression of frictionless soft granular matter by network analysis †

We consider dense granular systems in three spatial dimensions exposed to slow compression and decompression, below, during, above and well above jamming. The evolution of granular systems under slow deformation is non-trivial and involves smooth, continuous, reversible (de)compression periods, interrupted by fast, discontinuous, irreversible transition events. These events are often, but not always, associated with rearrangements of particles and of the contact network. How many particles are involved in these transitions between two states can range from few to almost all in the system. An analysis of the force network that is built on top of the contact network is carried out using the tools of persistent homology. Results involve the observation that kinetic energy is correlated with the intensity of rearrangements, while the evolution of global mechanical measures, such as pressure, is strongly correlated with the evolution of the topological measures quantifying loops in the force network. Surprisingly, some transitions are clearly detected by persistent homology even though motion/rearrangement of particles is much weaker, i.e. , much harder to detect or, in some cases, not observed at all.


Introduction
In non-Newtonian systems, complex fluids, 1 colloidal suspensions, 2 and especially granular matter in its flowing states, 3 the system's macroscopic transport properties are not constants but depend on various state-variables such as the volume fraction and the granular temperature. 4This interdependence and the presence of energy dissipation is at the origin of many interesting phenomena like clustering, 5 shear-band formation, 6 jamming/un-jamming, 7,8 dilatancy, 9 shear-thickening 2,10,11 or shear-jamming. 8,12,13While granular gases are very well described by kinetic theory, 4 dense granular fluids and granular solids are much harder to understand, in particular since they can transit from fluid-like flowing states to static, reversibly elastic ones.How and why they achieve this is still under debate, e.g., by irreversible (plastic) deformations, [14][15][16][17][18][19] related also to creep/relaxation phenomena, 9,[20][21][22] and many other mechanisms.
In this paper, we focus on the compression (loading) and decompression (unloading) of three dimensional dense, solidlike granular packings of frictionless soft particles, which allow to study a wide range of volume fractions, resembling soft rubber, or gel particles.Goal is a better understanding of precisely how these packings respond to (de)compression.For tiny deformations (isotropic strain), the considered systems can, and mostly do, behave fully elastically and reversibly.4][25][26][27][28] In analysis of these events, in addition to classical measures (pressure, energy, coordination number), we find it insightful to consider explicitly also the force networks, describing the interparticle interactions on a scale between particles and system size.By now it is well accepted that the static and dynamic properties of these networks are closely related to mechanical properties of the whole granular sample; the examples include granular samples exposed to vibrations/tapping 29,30 or pullout of an intruder from a granular sample. 31It is therefore important to be able to describe the properties of these networks in ways that are simple enough but also precise and more detailed than the classical measures alone, and that can be applied in three spatial dimensions (3D) rather than the rich 2D network analysis of recent years, see e.g.ref. 32-38 and the references therein.This is not an easy task since the force networks are built on top of contact networks, which for a deformed system evolve dynamically together interfering with each other.
While there are various approaches that could be used to describe the considered networks, they often do not satisfy some of the above requirements.In the present work, we focus on application of persistent homology (PH), which is a technique based on computational topology allowing for a concise description of 3D weighted networks, and is therefore appropriate for description of force networks, where the role of weight is played by the magnitude of the forces between granular particles.PH allows for a significant data reduction, since it describes a complicated weighted network in terms of point clouds, known as persistent diagrams, PDs, that contain information about the connectivity of the contact/force networks, and how such connectivity depends on the force magnitude.One strength of this approach is that all force magnitudes can be considered at once.This being said, the concept of a force threshold is inherent in PD calculations and the information about the connectivity for different thresholds is readily available, see ref. 39 for in-depth discussion in the context of granular matter, and 40 for an introduction to computational homology, of which PH is an important part.
PDs have been used extensively for a number of systems involving weighted networks, ranging from social 41 to fluid dynamics 42 and material science context. 43,44While applications to granular systems are still rather limited, the approach based on PDs has been used for the systems exposed to compression, 45,46 vibrations, 30,47,48 or shear. 49In the present work, we will illustrate some of the strengths of this approach, including its ability to easily carry out analysis/computations in 3D, as well as to consider evolving networks, including computing their differences.The computations are carried out using the GUDHI library. 50efore closing this introduction, we point out that in the present work we also consider relative large volume fractions, that may lead to large particle deformations, which in turn could produce new contacts between particles.8][59] In the present paper the soft-DEM numerical approach is used to model the compaction of ''soft'' elastic grains; in this approach the change of shape of the grains is neglected.Since the focus of the present study is to understand the force networks using persistence analysis and DEM simulations, which is already a complex enough investigation, detailed analysis involving deformable particles using multi-particle contact methods will remain a subject of a future study.
The remainder of this paper is structured as follows: After introducing both DEM and PD methods in Section 2, we proceed to discuss the results.Section 3 focuses on a preliminary discussion of the observations based on classical measures (pressure, energy), involving the irreversible events occurring frequently during both compression (loading) and decompression (unloading) of the considered granular systems.In Section 4 we present the results obtained based on PDs, as well as by relating the quantities derived from PDs to the classical measures.Section 5 is devoted to summary, conclusions and outlook.

Discrete element method
The approach towards the microscopic understanding of macroscopic particulate material behaviour is the modelling of particles using the so-called discrete element method (DEM), a numerical scheme originally formulated and developed by Cundall et al. 60 DEM is a straightforward implementation to solve the equations of motion for the translational and rotational (not used here) degrees of freedom for a system of many interacting particles: where m p = (4p/3)r p r p 3 is the mass of the i-th particle with The basis of DEM consists of the force laws that relate the interaction force to the overlap of two particles.This force can be decomposed into a normal and tangential component With known normal and tangential forces acting on all particles, one can numerically integrate the equations of motion, eqn (1), and obtain the position of particles in future time-steps.Next, we will describe the simple-most force law used in this research, ignoring tangential forces.
2.1.1Normal contact law.The elementary units of granular material are mesoscopic grains that deform under contact forces, possibly induced by an externally applied stress.For realistic modelling of the forces between two particles, we relate the interaction force to their overlap d.Note that the evaluation of the inter-particle forces based on the pair overlap may not be sufficient to account for the inhomogeneous stress distribution inside the particles and possible large deformations.Within this simplification, two particles only interact if they are in contact, resulting in the normal overlap: where ñ ¼ ðx p À xq Þ=jx p À xq j is the normal unit vector pointing from particle q to particle p.The linear spring-dashpot contact force law is a further simplification in DEM: [61][62][63][64][65] with k n as the spring stiffness between two particles.In reality, particles collisions are inelastic, i.e. energy loss occurs during collisions.Here, the dissipation is related to relative velocity ) of interacting particles with the viscoelastic damping constant for normal contact viscosity g n which is a intrinsic model parameter accounting for energy dissipation.This model implies that the spherical particles do not deform during the simulation and considers binary contacts only through a single point during their collisions.2.1.2Input parameters and units.The N = 4913 particles, with particle diameters drawn from a random homogeneous size distribution with maximum to minimum width, d max p / d min in the following with a prime, e.g., r 0 p = 2000 or d 0 p = 2, are used in the simulations shown in this paper.For working with units, there are two alternatives: either one can read the numbers in chosen units ‡ or the units are chosen based on physical properties to achieve non-dimensional quantities.The latter option is adopted here, i.e., the unit of length is chosen as the mean particle diameter, x 0 u ¼ d 0 p D E ¼ 2, so that hd p i = 1 is the dimensionless diameter.The second unit is the material density, r 0 u ¼ r 0 p ¼ 2000, so that one has the dimensionless density, r = f, and thus the unit of mass, m 0 u ¼ r 0 u ðx 0 u Þ 3 , i.e., the particle mass, m p = (p/6).For the third unit one has several choices, where we adopt here the units of elastic stress, s 0 u ¼ k 0 n =d 0 p , with the linear normal contact stiffness, k 0 n = 10 5 , which yields the dimensionless stress s ¼ s 0 d 0 p =k 0 n , and results in the unit of time In the chosen units, the dimensionless linear stiffness is m 0 u ¼ 1, and the linear contact viscosity, g 0 n = 10 3 , becomes , with background viscosity, g 0 b = 10 2 , or g b = g 0 b t 0 u =m 0 u = 4 Â 10 À4 .The consequent physically relevant properties are the restitution coefficient r = exp(ÀZ n t c ) E 0.855, with damping factor Z n = g n /(2m 12 ), reduced mass, m 12 = 0.063, and contact duration, 316, all considered for a contact between the largest and the smallest particle, with the larger viscous damping time-scale, t v = 2m 12 /g n E 5, and the even larger background damping time-scale t b = 2m 12 / g b E 50.Note that this choice of units corresponds to setting t 0 u / t 0 c , so that the dimensionless time-scale collapses simulations with different stiffness, in the elastic regime. 67,68o prepare samples, first particles are randomly placed in a 3D simulation box with a periodic boundary condition at a low volume fraction, well below jamming (transition from liquid-to solid-like behavior).Then, the simulation box is compressed homogeneously from all directions, up to volume fraction f max E 0.9, and then decompressed.Simulations are carried out using an isotropic strain-rate, _ e 0 v = 1.05 Â 10 À6 , i.e., _ e v ¼ _ e 0 v t 0 u ¼ 4:2 Â 10 À7 , for both loading and un-loading.In the following sections, we will present a few typical particle simulations from loading-unloading cycles, with the goal to zoom into this slow rate data, to display the dynamics of what is going on during plastic, irreversible re-arrangement events; results from different rates were reported in ref. 69.

Persistent homology (PH)
For the present purpose, one could think of PH as a tool for describing a complicated weighted network (such as the one describing interaction forces between the particles) in the form of diagrams, so called persistence diagrams (PDs).These diagrams are obtained by filtration, i.e. thresholding the strength of the interactions between particles.For the present purpose, this strength is quantified by the normal interaction forces.Then, the simplest PD, which we refer to as b 0 PD, essentially traces how 'structures' (could be thought of as 'force chains', without attempting to define them) appear as a filtration level is decreased, or disappear as two structures merge (b 0 stands here for the 0th Betti number that essentially counts the number of components or 'chains').In three spatial dimensions (3D), there are three Betti numbers, b 0 , b 1 (counting loops), and b 2 , counting enclosed 3D structures, and therefore also three PDs.Since in the present context b 2 structures are rarely present, we will not discuss them further, and will focus on b 0 and b 1 PDs.We note, and will discuss later in this section, that PDs contain the information which is significantly more detailed than the Betti numbers, since they tell us about connectivity of the components, and not only about their numbers.
To illustrate how PDs encode the information about a force network, Fig. 1 shows a toy example of a simple network and the corresponding PDs.appears at level 4 and, by convention, never disappears (or dies), leading to the point with coordinate (4, 0) on b 0 PD.On the level 3, we note that the existing component (4) becomes larger (this is not encoded in the PD, since connectivity has not changed), and three more components appear.On level 2, two of these components connect to the existing component (4) and therefore these two components disappear, or die, leading to a double point at (3, 2).Finally, on the level 1, another component from the level 3 connects to the existing component, leading to the point (3, 1).Since at level 1 all the nodes are connected, and there is just a single component existing, b 0 PD is complete.We note that this PD contains the information about the number of components and their connectivity for all considered thresholds, which is a major difference to the Betti number, b 0 , which is threshold dependent (in this example, there is one component on level 4, four components on level 3, two components on level 2, and one component on level 1).All of this is encoded in the b 0 PD, and we leave to the reader, as simple exercise, to extract the Betti numbers listed above directly from this persistence diagram.
Turning now our attention to the loops and b 1 PD, we see in Fig. 1(a) that we need to go down to the threshold level 2 to form a loop (a loop forms when the weakest link appears; in this case we have two edges with the force 2 in the loop in the lower part of the figure).This loop has coordinates (2, 0) in b 1 PD.Then, on level 1, another loop (central) closes, leading to the (1, 0) point in b 1 PD.In the present simple toy example we do not observe disappearing or dying of loops: in a more complicated network, bigger loops 'die' when they get (completely) filled by edges of non-zero strength.In a granular system, this could correspond to the formation of crystalline zones, where eventually all neighboring particles interact by non-zero contact forces in a hexagonal network; see ref. 39 for further discussion of this point.
We conclude this brief discussion of PDs by pointing out one important feature: one can show that PDs are stable, in the sense that small changes of the underlying network lead to small changes of a PD; such a result can not be proven, e.g., for Betti numbers. 39In practical terms, and in the context of granular physics, this means that small errors in measuring forces between particles would lead only to correspondingly small errors in the measures derived from the PDs.
It should not come as a surprise that the PDs resulting from simulations of a large number of granular particles are significantly more complicated than the ones resulting from the toy example.
While there are many possible measures that could be introduced, one simple choice is based on the idea that there are two appropriate measures/numbers (i) points (generators) in a diagram, and (ii) the threshold range over which a point persists -both describing the force network.Using landscape as an analogy, the number of points in b 0 PD specifies the number of (mountain) peaks, and their lifespan (see below) describes how well developed these peaks are, or, to push the landscape analogy further, how high they are compared to the 'valleys' that surround them.The lifespan L, is introduced as the difference between birth (B) and death (D), coordinates, so L = B À D. Both measures could be combined into one by defining the total persistence, TP, as a sum of all lifespans.Then, for our toy example shown in Fig. 1, b 0 TP = 8, and b 1 TP = 3.We will be using TP later in the paper to quantify the force networks obtained from DEM simulations.
During the last decade, the software for computing PDs and additional measures has been developed by a number of research groups and made publicly available.All the results presented in the present work were computed using the GUDHI library. 50

Classical analysis of DEM results
In this section we compare the loading and unloading branches at the full range of volume fractions.Focus here is on DEM particle simulation output/information, whereas in the next section, the properties of force networks and their relation to the classical measures are discussed.
Fig. 2 shows the dimensionless pressure, p ¼ p 0 d p =k n plotted against volume fraction f.On the scale of the whole compression/decompression cycle, the differences in pressures between the branches are rather small, and the pressure curves look rather smooth.As we will, in what follows, zoom into these curves, rich structures show up to be explored deeper.The symbols in the figure show the location of some of the events which we will consider in more detail.It is worth mentioning that the (isotropic) inertial number for these simulations is of the order of 10 À5 , already for rather small p Z 10 À4 .Note that while faster compression rates lead to different results, simulations with slower deformation rates (not shown here for brevity), 69 are practically identical on this scale, confirming that the simulations shown in Fig. 2 are in the quasi-static, rateindependent regime.
Note that volume fractions in Fig. 2 go as high as 0.9.Recent studies 56 have shown that classical DEM is not accurate at such high compression, due to extreme deformations, so that particles are not spherical anymore, with additional, multiple contacts established.Different methods, e.g., DEM-FEM coupling, could be used to simulate high volume fractions, 52,70,71 however such methods are computationally costly.For simplicity, in the present work we stick to classical DEM, following the phenomenology in a qualitative way, leaving more quantitative considerations of the influence of the change of particle shape (and other aspects related to large volume fractions) to other studies 53 and future research.
Fig. 3 shows the energy ratio (kinetic to potential), E k /E p , of the system during loading and unloading.It is a known fact that the energy ratio is one way to identify jamming: it drops rapidly below unity at jamming, 7 while it sharply increases at unjamming, but at a different, larger value of f.This figure also illustrates that the process of loading and unloading is not smooth.Slow isotropic deformations results in an energy ratio base-line But the kinetic energy is also characterized by numerous events, i.e., peaks of energy ratios up to E k /E p E 10 À3 , still well below unity.This indicates that the system remains jammed even though a considerable kinetic energy shows up rapidly and decays (exponentially) fast, after each event.A closer zoom into some selected events follows below.
Next, we consider the coordination number C* which is related to the contact number per particle C = 2M/N (given the total number of contacts M) by C* = C/(1 À f r ), in which f r is the fraction of rattlers.Fig. 4 shows that the coordination number, C*, is quite similar for the loading and unloading branches, most different at the lower volume fractions, close to jamming and unjamming.On this scale, like for pressure, the tiny fluctuations, variations and differences can be seen only when zooming in to the data, as we will do in the next subsection.It must be noted that large deformations may lead to development of new contacts between particles 51,53 and that considering the influence of these newly developed contacts on the results in this paper remains to be carried out in the future study.

Loading and unloading events
To proceed further, three different windows of volume fraction for both loading and unloading branches were selected.
We proceed by comparing the loading and unloading branches at a few different selected values of f so to be able to zoom into windows of Df = AE0.0012around the events marked in Fig. 2-4 and summarized in Table 1.
Fig. 5 shows the pressures around the selected values of f at which events take place.We note that p, differs significantly between loading (L) and unloading (UL).Only the last event on the loading branch is followed by an elastic regime that is exactly, reversibly traced back during unloading from the maximum volume fraction, f max = 0.9013, confirming reversibility for small enough strain and strain rates.
All the events are characterized by (sometimes tiny, sometimes large) changes in pressure: for L we always observe pressure drops, while for UL, we also observe increases in pressure (note that the UL figures need to be viewed from right to left since f is decreasing).As we will see in what follows, other considered quantities show significant changes during the event windows considered in Fig. 5, even when the pressure changes are almost invisible.
Fig. 6 shows the energy ratio, E k /E p , for the zoom-in windows.We find that this ratio is more sensitive to a change of f  Table 1 Summary of the event windows marked in Fig. 2-4 and detailed as zooms in Fig. 5 than pressure.Note that the energy ratios during loading and unloading do not follow the same trend -even if f is the same.This is not surprising, since the sample configuration has changed.We note that each event is followed by an exponential decrease in energy ratio.The corresponding dissipation rates (data not shown) depend on the restitution coefficient as p(1 À r 2 ) and, somewhat, increase with volume fraction, due to an increased number of contacts, Fig. 7 shows the coordination number, C*, for the same data range as in Fig. 5 and 6.Similar to the energy ratio shown in Fig. 6, the coordination number plot reveals irreversible rearrangements during events, showing that the events (at least the ones considered in this figure) correspond to changes of the microstructure (coordination number) associated with plastic, irreversible deformations of the granular sample (in contrast to elastic, reversible deformations during which the coordination number varies continuously (i.e., outside of the event window, without significant jumps).We will discuss later in the paper also situations where events could occur without particle rearrangements, therefore without strong changes of the E k /E p ratio, or the coordination number.
For the events considered so far, Fig. 7 shows that an event always starts with a coordination number decrease.This is due to rearrangements of particles, which leads to an increase of the E k /E p ratio, as can be seen in Fig. 6.For an explanation why such kinetic energy (and displacement) fluctuations result in a Fig. 5 Pressure, p, plotted against volume fraction, f, from the same data as in Fig. 2, zooming into three events on each branch ( and symbols), for loading (L, red) and unloading (UL, blue).The same color/symbols are applied to Fig. 6 and 7. (a) L1, (b) L2, (c) L3, (d) UL1, (e) UL2, (f) UL3, see Table 1.
Fig. 6 Kinetic to potential energy ratio, E k /E p , plotted against volume fraction, f, from the same events as in Fig. 5.
drop of coordination number, see the discontinuous generalized distance probability density functions in ref. 72.Once an event is completed, the coordination number typically recovers close to its original, elastic trend, sometimes below (a and b), sometimes above (c-f).

Events force networks
In Fig. 8, force networks are plotted during events, for the loading branch only, with a threshold set such that the plots have similar volume fraction, even with a few holes/gaps that indicate extended islands of weak contacts only.The intensity of forces is color coded the same way in all plots, where blue, green, red range from low to high forces.
The four snapshots shown in Fig. 8 for each of the events L1 (f = 0.6735), L2 (f = 0.7427), and L3 (f = 0.899), are equally spaced within the event window, Df E 0.00023, so that biggest changes are expected between the two middle panels.However, given the rather small relative changes in pressure during events, also the force intensity change is practically invisible, and only a few modifications of edges are visible.This observation brings us to the conclusion that the visualisation of the force network is not very helpful to understand the evolution during events; further analysis of these networks based on persistent homology method will be discussed in Section 4. Before that, we discuss briefly the non-affine displacement.

Non-affine displacements during events
The main difficulty in describing the mechanics of granular materials is how to deal with disorder.The simplest approach is to ignore disorder altogether, and attempt to gain insight based on ordered, crystalline packings. 73This assumption considers a uniform strain at all scales, with the displacement field of the grains following the macroscopic deformation, socalled affine motion.However, a number of studies [74][75][76][77][78][79][80] have shown the failure of only considering affine motion in describing the mechanical behavior of granular packings.
The basic idea of elastic theories relevant to granular packings is that the macroscopic work done deforming the system is equal to the sum of the work done on the level of each graingrain contact.The latter scale has to be replaced by a suitable split to average (affine) and fluctuating (deviating, non-affine) quantities.Thus, the displacement of a particle can be written as: where x i (t) is the position of particle i at time t, the subscripts indicate the affine, non-affine reversible, and the event displacements.Note that the last two terms are different contributions to non-affine deformations, on top of the first, which accounts for the affine deformations only.The first two terms provide the smooth total displacements of particles during elastic, reversible phases, both affine plus non-affine.As we will see later, the last term accounts for the much more violent dynamics during non-reversible deformations (events); it is not smooth in time (volume fraction).
Next we focus on the non-affine displacement from before to after an event, which involves Dx i À Dx i affine .The reversible nonaffine displacement is typically much smaller than the irreversible event displacement, so that we do not attempt to separate it off (data for Dx i reversible not shown here separately).In order to visualize the evolution during the events, the non-affine displacement of all particles during the events (from f = 0.673 to 0.899) is plotted as thin lines with a stretch factor in Fig. 9.The figure shows the events from the loading branch for the system with periodic boundary conditions (grey bars).The particles with displacement above a certain threshold are plotted in color, whereas the particles with little displacement are omitted.The first and the last event are remarkable, since they show up strongest, both with a strong displacement of Fig. 7 Coordination number, C*, plotted against volume fraction, f, from the same events as in Fig. 5. many particles, covering the whole system.It was found that in the low volume fraction case, the system is not yet fully jammed and the event is the superposition of several smaller sequential events.Even for very large f, significant events can occur, even though the system is already strongly jammed/ over-compressed.Also some events show considerably less particles involved.However, in all cases the whole system is involved, indicating possible relevance of finite size effects.Fig. 10 shows an event from the unloading branch starting at f = 0.899 until f = 0.673.For unloading, the first event (left)  is very almost not showing up for the given thresholds, and also the last event (right) is rather weak.The three marked/ selected events in the middle show different characteristic patterns well above the threshold.All events involve groups of particles of order of half system size, indicating that these are almost system spanning events, with possible relation to the finite size of the system.Additionally, comparing non-affine figures for loading and unloading events, it is found that the non-affinity is more pronounced for loading than for unloading events.
4 Force network analysis based on persistent homology

Global view
The force networks, see Fig. 8 for examples, are difficult to analyze visually.As discussed in the Introduction and Section 2, persistent diagrams, PDs, provide one complementary approach to understand granular packings under slow deformations more quantitatively.Fig. 11 shows an example of b 0 and b 1 PDs for volume fraction f = 0.6739.Before proceeding with their quantification, we note a couple of issues that need to be considered when computing PDs for a complicated 3D network: (i) one needs to decide how to treat the forces across periodic boundaries: in our implementation, we treat these forces in the same way as all the other ones, resembling a system without wall-induced inhomogeneities or boundaries; (ii) one needs to decide how to treat so-called trivial loops formed by three particles in contact: in the current implementation, as a choice, these loops are ignored, what in practical terms means that larger loops can 'die' when they get completely filled by trivial loops.The interested reader is referred to ref. 39 for more details.
The PDs are essentially point clouds, and one needs to come up with appropriate measures to quantify their features that describe the corresponding contact-and force-network in granular matter.As discussed already in Section 2, a particularly simple quantification of the connectivity properties of the underlying networks is obtained when the information contained in the PDs is compressed into a simple measure, i.e., total persistence, TP, which could be considered separately for components ('chains') and loops.Fig. 12 shows TP for components (b 0 ) and for loops (b 1 ), together with the pressure (already plotted in Fig. 2) for both loading and unloading branches.TP for b 1 (TP1) closely follows the pressure, while TP for b 0 (TP0) is much more noisy.In particular, while the loading pressure is always larger than the unloading one, the same does not apply to TP0.We note that the evolution of TP0 is particularly nonuniform for large f's, suggesting significant changes in the force network.
Next, the finer features of both pressure and TP are discussed.Before that, however, we quantify the correlation between pressure and TP, by computing their correlation functions.Since the data are rather noisy, we consider the data averaged over a moving window, as described next.
The moving correlation between the two time series x(t) and y(t) is calculated by first finding the covariance between the subseries x i and y i in the moving window centered at time = t 0 , where W t 0 w is a selected window which centers at t 0 and w is the total number of data points in W t 0 w .Also, % x t 0 , % y t 0 are the means  of x i y i , respectively.The standard deviation, s t 0 w ðxÞ, of x i at time t 0 is with a corresponding expression for s t 0 w ðyÞ.Then, using the standard deviations, the appropriately normalized correlation function is r t 0 w ðx; yÞ ¼ c t 0 w ðx; yÞ s t 0 w ðxÞs t 0 w ðyÞ ; which is unity for perfect correlation, and can go down to negative unity for perfect anti-correlation.Fig. 13 shows the correlations between the pressure and total persistence TP (separately for TP0 and for TP1).We immediately observe that the correlation between pressure and TP1 is very high, in particular for large values of f.The correlation decreases for a few isolated events, which we discuss further below, but even for those events the correlation is still very strong.This result suggests that the properties of the loops in the force network are strongly correlated with the pressure. 38On the other hand, the correlation between TP0 and pressure is much more pronounced, in particular during the rapid isolated events showing as negative peaks in the correlation curve, indicating strong anti-correlation.Therefore, the topological properties of the components are much less representative of the pressure in the system.We proceed by discussing in more detail the fines structures of the evolution during selected events on the loading and unloading branches.

Analysis of the events
As already discussed in Section 3.1, and further illustrated by Fig. 13, a considerable number of events (tens or hundreds, depending on how an 'event' is defined) occur -in particular during loading.In what follows, we discuss further the evolution of the force networks for the events already considered in Fig. 5-7, with the goal of understanding such events more precisely.
We focus first on loading and consider (on a much finer scale now) pressure and total persistence, TP.Fig. 14 shows the results, where all quantities are normalized by their value at the peak of the event.Pressure is plotted as Dp = (p À p*)/p*, where p* is the reference value as marked with symbols on Fig. 2 and reported in Table 1.Total persistence, TP, is plotted as DTP = (TP À TP*)/TP*, where TP* is the total persistence of corresponding volume fraction given in Table 1.Even on this very fine scale, and during the event itself, TP1 closely follows the pressure.TP0 changes much more dramatically, showing that the topology of the force network changes during an event, but not all of these changes resemble pressure modifications.

Soft
Note that a single PD provides information the state of the system, at a given time, and therefore does not include any information about the dynamics.Another crucial feature of the PDs comes handy here: PDs live in a metric space, and therefore could be compared, see ref. 39 for a more detailed discussion of this point.This means that one can define a concept of distance between PDs, which could be thought of as a minimum (Euclidean) distance by which the points in one diagram need to be adjusted to map them exactly to each other; if the number of points in two PDs is different, the extra points are mapped orthogonally to the diagonal.
One type of distance which is often used is the degree q Wasserstein distance, W q (PD, PD 0 ), which is essentially the L q norm of the minimum distance required to map the points from PD to PD 0 (minimized over all possible mappings).The value of q puts emphasis on larger differences between the points, for q 4 1, or on smaller ones, for q o 1.We use q = 2 in the present work and refer the reader to ref. 39 for a more detailed discussion of this concept.
Fig. 15 shows the same event from Fig. 14, but from a different perspective.Here we plot the distance between the persistence diagrams, together with the relative change of the kinetic energy of the particles, scaled by the peaks, such as E Ã p are given in Table 1.Likewise other parameters, W 2 distances were scaled using of the peak value of the event such as We see that the changes in W 2 distances follow the changes in kinetic energy, although not too closely.
The fact that the W 2 distances before and after an event are similar in Fig. 15 inspires another question: is the topology of the force network before and after an event similar?To answer this question, we compute heat maps, which show the distances between all the persistence diagrams describing force networks during an event.Fig. 16 shows these heat maps (both for W 2 b 0 and W 2 b 1 distances) for event L1.This figure shows the heat maps as a (symmetric) square matrix, where the color of element (i, j) illustrates the distance between this element and the (i, i) one (here i counts rows, and j counts columns of the heat map); by definition, the elements on the diagonal are zero.We see that the distances between the elements progressively increase as we move away from the diagonal, meaning that the topological properties of the system before and after the event are significantly different.The blue square regions (top left and bottom right) show that before and after the event, the distances between the PDs are very small (compared to the changes during the event itself).The difference from before to after an event is considerable, however, as indicated by the dark red color in the off-diagonal corners.
To complete the discussion of the loading phase, Fig. 17 shows the results for events L2 and L3.We see that the general

Soft Matter
This journal is © Royal of Chemistry 2022 features of the results similar as already discussed in the context of event L1 -only these events at larger volume fractions are narrower, more localized.While TP1 closely follows the pressure, TP0 could increase or decrease during an event.Therefore, a general pattern emerges: properties of the loops that form in the force networks as characterized by TP1 (which combines the number of loops and their 'intensity' (distance from the diagonal), are closely related to the material response, as quantified by the pressure.Next, we proceed by discussing unloading.Fig. 12 already suggested that this phase is smoother than the loading one, with lesser and weaker events, at least on the coarse scale.On the fine scale, we consider the same three events as in Fig. 5-7, with the results shown in Fig. 18.The general features of the results are similar to the ones observed for loading, with TP1 closely following the pressure, and TP0 evolving in a much more dramatic fashion, either increasing or decreasing during an event.The relative changes are behaving similar to the kinetic energy ratio.Fig. 18(c) and (d) show an interesting feature at f E 0.7049, where the distances between PDs show changes in the topology of force networks, although kinetic energy is unchangedshowing that changes in the topology can occur even without any motion of the particles.Note that also the coordination number is not changing much, see Fig. 7 bottom-centre, however, with stronger fluctuations.
Further results could be used to quantify the changes occurring during events by considering in more detail the ingredients that define TP: average lifespan of generators and their numbers.However, these additional results do not show obvious common features of the events.Further research, perhaps including a more detailed analysis of the persistence diagrams may be needed to better understand the details of the force network evolution during plastic events for either loading or unloading.
To conclude, we use the fact that a significant amount of information can be condensed in a visually simple-to-digest form using heat maps, and show in Fig. 19 and 20 the heat maps for a set of ten loading and unloading events.Fig. 19 shows that the distances are becoming larger as the system progresses to larger f's during the loading, and Fig. 20 shows that the distances are becoming smaller for smaller f's while unloading.The bigger the blue and red squares, in most higher volume fraction cases, the sharper the event is localized in strain (volume change), whereas only a few events closer to jamming are not localized, diffuse, or overlapping multiple  events.features of the heat-map results do not to change significantly between loading and unloading, or between components and loops.
Finally, we comment on the aspect, which has not been discussed extensively so far: while most of the events that we analyzed involve considerable particle dynamics, i.e., a significant kinetic to potential energy ratio, there are also events which involve only rearrangements of the force networks, with much weaker particle motion (data not shown, orders of magnitude smaller energy ratios).To illustrate this feature of the results, Fig. 21 shows the distances between consecutive persistence diagrams and the energy ratio for two selected loading/unloading windows.Careful inspection shows that while most of the events do involve changes of both W 2 and /E p , there are a few events only changes in W 2 are present, that force networks can change while particles remain stationary.

Summary and conclusion
In this paper, the process of slow compression (loading) and decompression (unloading) was studied for a three dimensional granular model system of soft particles.While below jamming the system evolution is irreversible (fluid-like), one main finding is that above jamming, both loading and unloading proceed via either reversible (elastic-) or irreversible (plastic-solid-like) processes.While the former are smooth and reversible, the latter are associated with rapid events, with serious rearrangements of the particle micro-structure, their positions and contacts, with associated particle dynamics (kinetic energy).However, there are also events that are associated with the changes in force networks persistence diagrams (PDs) only, without similarly strong particles dynamics (kinetic to potential energy ratio).
The analysis of the force networks is carried out using PDs, a technique that can quantify the evolution during both reversible and irreversible processes.While PDs allow the extraction of a number of features of the underlying force networks, here we focus mostly on the most condensed measure, the total persistence, TP.It encodes information about complex weighted three dimensional (3D) networks into a single number.Even under such extreme information reduction, we find that TP provides interesting insights about the underlying force network and its evolution.In particular, we find that the changes of the mechanical pressure are strongly correlated with the TP computed for the loops that spontaneously form in evolving force networks during (de)compression.This finding suggests that the force network loops play a significant role in determining the mechanical system properties.
Further insights regarding the evolving force networks can be obtained by comparing consecutive system configurations.Such a comparison is made meaningful by computing the differences between the persistence diagrams that describe the main features of the underlying network.The difference, known as Wasserstein distance, shows that the force network experiences a significant change during the events discussed above, with strongly different force networks before and after an event.It should be emphasized that for obtaining such simple but insightful results, the data reduction provided by PH is crucial.This data reduction leads inevitably to an information loss, however we find that even under the most dramatic data compression, the emerging information still allows for reaching useful and relevant insights.
Future work should be devoted to a more quantitative analysis of strongly compressed and sheared frictional, soft granular packings -in the same spirit as present -where the   deformability is fully taken into account, and also developing new approaches to consider such deformation of real particles, far beyond the validity of the over-simplified Hertzian based models.We expect that further significant insights into the micro-mechanics of soft granular matter, using further detailed features of PDs (like generators life-spans and numbers).Particularly interesting is the influence of inter-particle friction on the force networks, which could be better understood by a simultaneous analysis of the properties of both normal and tangential force networks.
Fig. 1(a) shows a small set of particles (nodes) interacting by a force represented by the numerical values (intensity) along the edges connecting the nodes.The PDs keep track of how components appear and disappear as the force threshold is decreased from the highest value (4 in the present example) down to the lowest value (here 1).To explain how PDs are produced, let us first focus on b 0 PD.The first component (single edge) with value of 4 is the first one that

Fig. 1
Fig. 1 (a) Toy network with the strength of interaction between nodes shown by numerical values and color intensity.The corresponding persistence diagrams are (b) b 0 PD, and (c) b 1 PD.Note that in (b) the point with coordinates (3, 2) is double, since two components disappear (plotted as two slightly shifted points for visualization purposes).‡ Units could be, e.g., length, x ˆu = 0.5 Â 10 À3 m, time, t ˆu = 10 À5 s, and mass, m ˆu =

Fig. 3
Fig.3Kinetic to potential energy ratio, E k /E p , plotted against volume fraction, f, for loading (L, red) and unloading (UL, blue).The symbols indicate the location of three selected simulation event windows, for both L and UL.

Fig. 4
Fig. 4 Coordination number, C*, plotted against volume fraction, f, for loading (L, red) and unloading (UL, blue).The symbols indicate the location of three selected simulation event windows, for both L and UL.

Fig. 8
Fig. 8 Force networks during the three loading events L1, L2, and L3 (in the same scaling, where each event is represented by a row of four snapshots), with the forces/overlaps plotted only above a threshold of the overlap, d/(2R) 4 0.006, 0.04, and 0.11, respectively; all lines are color coded the same way, from small (blue), moderate (green) to large (orange) forces.Forces across the periodic boundaries (light grey boxes indicate the outside) are omitted.

Fig. 9
Fig.9Non-affine particle displacements projected to a plane, during event-windows with Df E 0.0012, plotted for a few representative events from begin (f = 0.673) to end (f = 0.899) of loading (L), i.e., volume fraction increasing from left to right (as visible from the grey bar).The thin black lines give the scaled displacement, and the circles indicate particles with displacement above a threshold of 0.1, with blue, green and red indicating moderate, large and very large displacement amplitude.

Fig. 10
Fig.10Non-affine particle displacements projected to a plane, for event-windows with Df E 0.0012, plotted for a few representative events from begin (f = 0.899) to end (f = 0.673), of unloading (UL), i.e., volume fraction decreasing from right to left (as visible from the grey bar).The thin black lines give the scaled displacement, and the circles indicate particles with displacement above a threshold of 0.1, with blue, green and red indicating moderate, large and very large displacement amplitude.

Fig. 11
Fig. 11 Examples of persistence diagrams, PDs, for (a) components (clusters), b 0 , and (b) for loops (cycles), b 1 for f = 0.6739.The axes show the force level at which a structure appears (B = birth) and disappears (D = death).The lifespan l = B À D (distance from diagonal marked as green horizontal line) illustrates the force threshold range over which a structure persists.See ESI † for animations showing the evolution from f = 0.6739 to f = 0.6750 (discussed later in Fig. 14-16).

Fig. 14
Fig. 14 Normalized values for the pressure change Dp (red solid), DTP0 (blue solid) and DTP1 (blue dashed) for event L1; all shown quantities are normalized by their value at the center of the event; for example, pressure is plotted as Dp = (p À p*)/p*, where p* is the reference value shown by J.The meaning of colors and line patterns is the same in the figures that follow.

Fig. 15
Fig. 15 DW 2 distances both for b 0 (blue solid) and b 1 (blue dashed) as well as relative change of kinetic energy DE ratio (red solid) for event L1, and the arrow is pointing to the reference value used for normalization.The meaning of colors and line patterns is the same in the figures that follow.

Fig. 16 (
Fig. 16 (a) Heat map (a symmetric matrix) of DW 2 for b 0 of the event L1, where the color of element (i, j) illustrates the normalized distance of PDs for b 0 between f = f i and f = f j , normalized by the same reference value used in Fig. 15, (b) normalized DW 2 for b 1 distance.The minimal steps between i and j are Df E 0.00014.

Fig. 21
Fig. 21 Energy ratio (red solid) and W 2 distances for b 0 (blue solid) and b 1 (blue dashed) for loading (a) and unloading (b) against volume fraction, f, for a range between 0.68 and 0.695.