GoKawiilThe optical response was then correlated with the cell potential in Fig. 2a for different regions of the particle (map shown in Fig. 2d) and as a function of cell temperature. The sharp–gradual–sharp intensity trend (with respective regimes labelled (i), (ii), (iii) in Fig. 2a top right) occurred both during lithiation and delithiation across all temperatures. The voltage hysteresis between the delithiation and lithiation cycles increases at lower temperatures, but the small step-like intensity changes appear throughout all datasets in all three regimes (i)–(iii).
Fig. 2: Optical step changes correlated with cell potential. Full size image a, Mean optical intensity curves from different particle regions compared with cell potential. Lithiation and delithiation curves are both shown, with direction denoted by black arrows, normalized to the start and end of delithiation and lithiation, respectively. Only regions B, C and D are shown for simplicity. Where appropriate, the curves have been cut before the onset of the stage 2L to 2 transition for clarity. Three regimes in the intensity trends are marked (i), (ii) and (iii). b, Differential capacity (dQ/dV) curves from the charge–discharge data at different temperatures (middle). Four processes labelled (1), (2a), (2b) and (3) can be distinguished. The onset of a peak corresponding to the dense transition (2L–2) is seen in these plots at temperatures 35 °C and below. The black arrows denote the direction of lithiation compared with delithiation. Scatter plot of all step events that occur during delithiation (top) or lithiation (bottom), detected from the mean optical intensity as a function of temperature (region B, cycles 1 and 2). Each circle represents one step event, with the x-position denoting the potential and the size denoting the intensity. c, Scatter plot of all detected step events in different regions of the particle (45 °C delithiation, cycles 1 and 2). The x-position denotes the potential, and the size denotes the intensity change. Cell potentials corresponding to processes (1), (2a), (2b) and (3) as found in b are marked with dotted lines. d, A map of particle regions, A–D, used for mean intensity calculations. Scale bar, 5 μm.
To further quantify the intensity steps, we used a step-detection algorithm to analyse the intensity traces from individual graphite regions (Supplementary Information section 3.3). This approach allowed us to build up statistics for these events as a function of temperature, cycle, voltage and region dependence. Each detected step event is plotted as a circle against the voltage at which it occurred (Fig. 2b, top, delithiation; bottom, lithiation). The size of the plotted circle encodes the event size (that is, magnitude of the step fall or rise). For simplicity, we have limited ourselves to events in region B of the particle in this plot. The distribution of the events is then correlated with a differential capacity analysis (dQ/dV) of the cell electrochemistry, which highlights underlying phase transitions or processes. Figure 2b (middle) shows the dQ/dV plots at different temperatures, with prominent peaks marked with numbers. The first peak during lithiation at 197–202 mV (1) is attributed to the 1′L → 4L transition. There are then two small peaks around 162 mV (2a) and 140 mV (2b) (more pronounced at higher temperatures, see Supplementary Fig. 11), and an additional peak (3) is visible as a shoulder of the large truncated peak (assigned to the onset of the 2L → 2 transition). During delithiation, similar trends are observed, with the smaller peaks (2a) and (2b) becoming more prominent. Based on this analysis, the two sharp drops in optical intensity shown in Fig. 2a, (i),(iii) are attributed to the differential capacity peaks (1) and (3), respectively, whereas peaks (2a,b) occur during the more gradual intensity changes (Fig. 2a, (ii)). Furthermore, the step-like events detected in Fig. 2a (top and bottom) cluster around the dQ/dV peak positions, with the events becoming more tightly clustered at higher temperature. This is closely related to the dQ/dV peak shape, which sharpens at higher temperature.
The same analysis can be extended to other regions of the particle, and a selection is shown in Fig. 2c for the 45 °C delithiation cycle. Regions A, B and C show consistent trends as discussed, but region D shows events that are more widely distributed, even at 45 °C. This was true generally, with some regions consistently showing a broad temporal spread of events at all temperatures (Supplementary Fig. 32). As shown in Extended Data Fig. 2, these regions (D and E) also experience more frequent smaller-sized events, whereas in other regions step events are more broadly distributed in size.
The step-like intensity features of the individual regions during the dilute stages are highly reminiscent of avalanches2,3: abrupt processes with a broad distribution of sizes, which occur generically in disordered systems under a slowly applied driving force. Complex experimental trends such as avalanche statistics and transition dynamics can be understood qualitatively using general theories for these phenomena. This motivated us to develop an RFIM to simulate the filling of a graphite lattice with discrete lithium occupation with attractive intralayer coupling, repulsive screened interlayer coupling36,37,38, and an on-site random potential representing idealized disorder (Fig. 3a). We considered both two-dimensional (2D) and three-dimensional (3D) models (for full description, see Supplementary Information section 2.1). Further justification and experimental comparisons can be found in Supplementary Information section 1.2. As shown in Fig. 3b, the simulated filling behaviour of this model under a slow potential ramp exhibits both avalanches (sudden jumps in fill fraction) and the sharp–gradual–sharp filling trend observed in the experiment. These arise because static disorder results in energy barriers to filling that are large enough to disrupt the dynamics of system-wide transitions between equilibrium phases at experimentally relevant timescales and temperatures. In our model, an avalanche usually corresponds to a local domain filling in a single or a few correlated layers (see Supplementary Fig. 27 and Supplementary Video 5). Our model also implies that pure staging does not exist during the reaction (Fig. 3c) because of the intrinsically nonequilibrium dynamics of avalanching systems.
Fig. 3: RFIM set-up and filling behaviour. Full size image a, Schematic of the intralayer (J) and up to third-neighbour interlayer (K1, K2 and K3) coupling terms in the 2D model, the colours indicate different distances between two Li in the c-direction. b, Mean lattice fill fraction (ϕ) compared with Monte Carlo steps (mcs) during a single ramped potential simulation (top), with disorder (\(\mathop{h}\limits^{ \sim }\) = 0.4), and temperature (\(\mathop{T}\limits^{ \sim }\) = 0.125), with total mcs = 80,000. The red circles denote time points (i)–(iii) from which the lattice pictures in c are taken. The red square inset highlights a portion of the filling curve, which demonstrates stepped (avalanche) behaviour (mcs = 45,700–47,200). c, The lattice structure shown at different time points (i)–(iii). Occupied sites are drawn in white, and unoccupied sites are coloured according to the respective c-axis separation between the occupied sites.
The model can be used to understand the avalanche size distribution and temporal spread in different experimental regions. For example, in Extended Data Fig. 3, we demonstrate that a higher disorder strength results in more frequent smaller avalanches and a broadened temporal spread, similar to those in the experimental particle regions D and E. The region-dependent trends (regions A–E in Fig. 2 or 0–9 in Fig. 4) seen in the optical measurements and the graphite disorder were explored using ex situ EBSD and atomic force microscopy (AFM) measurements (Supplementary Information section 1.1). The grain reference orientation deviation and the particle topology maps obtained from these two techniques, respectively, show similar regions or domains to those seen optically, with regions D and E showing a higher density of folds and bends and thus more static disorder. Further model–experiment comparisons can be found in Supplementary Information section 2.3.
Fig. 4: Spatio-temporal correlations during 1′L–4L transition. Full size image a, Demonstration of concurrent activity mapping. Cross-correlation calculated from every pair of pixels (5 × 5 binned) on the particle, using the time-differential signal (cycle 1 at 15 °C). The resultant cross-correlation function is summed for time lags (−5 s ≤ τ ≤ 5 s), and the resultant symmetric cross-correlation matrix is clustered into 12 blocks using spectral co-clustering (top). These 12 clusters are then mapped back into space, which produces regions on the particle that show concurrent activity (bottom). We exclude two clusters (edge pixels) in further analyses, marked with an asterisk. b, Identification of the most and least connected regions. The correlation function from the coarse-grained dataset (cycles 1 and 2 at 15–45 °C) is summed over time lags (−5 s ≤ τ ≤ 5 s). All inter-region correlations are then summed, excluding self-correlation, and mapped to the particle. c, Demonstration of inter-region, positive time-lag correlations. The correlation function from the coarse-grained dataset (cycle 1 at 45 °C) is summed over time lag (0.5 s ≤ τ ≤ 3 s) and the off-diagonal elements are interpreted as successor–predecessor relationship frequency between two regions (top). The relationship is mapped as a directed, weighted graph (bottom), in which each region (node) is connected by a weighted, directed edge defined as the difference of reciprocal values (i, j) − (j, i). The arrow points from the predecessor to the successor. Twenty edges with the highest weights are shown. d, Time-lag-dependent intraparticle correlations during delithiation and lithiation. All predecessor–successor relationships are identified from the dataset (cycles 1 and 2 at 15–45 °C) as per c. The reoccurring edges are mapped onto the particle at different lag times. The potential window corresponding to the 1′L → 4L transition was used for all analyses.
We note that this model is meant as a general qualitative framework and uses an idealized representation of quenched, uncorrelated disorder as a proxy for material properties such as turbostratic disorder39 and elastic deformation40. A discussion of model limitations and potential extensions can be found in Supplementary Information section 2.5.