The samples were prepared as in Ref. [37]. Shortly, sterically stabilized polymethylmethacrylate (PMMA) spheres of diameters σ _l = 3.10 μm (polydispersity 0.07) and σ _s = 0.56 μm (polydispersity 0.13) are dispersed in a refractive index and density-matching mixture of cis-decalin/cycloheptylbromide to form binary mixtures with different mixing ratios and total volume fractions. The use of a density-matching solvent avoids sedimentation effects on the time- and length-scales probed in our experiments [37]. The smaller particles are fluorescently labeled with nitrobenzoxadiazole. Upon salt (tetrabutylammoniumchloride) addition, this system presents hard-sphere–like interactions [51, 52].

Confocal microscopy experiments were performed by using a Nikon A1R-MP confocal scanning unit mounted on a Nikon Ti-U inverted microscope with a 60× Nikon Plan Apo oil immersion objective (NA = 1.40). We acquired 10⁴ images with 512 × 512 pixels (corresponding to 107.5 × 107.5 μm²) at 30 fps by focusing on a plane at a depth of approximately 30 μm from the coverslip. The confocal images were acquired with the maximum pinhole diameter allowed by the microscope, 255 μm. The experiments were performed at a temperature of T = 293 ± 2 K.

A DDM analysis of the confocal images was performed as described in Refs. [20, 36]. In brief, we calculated the image structure function D ( q , Δ t ) = ⟨ | I ̂ ( q , t + Δ t ) − I ̂ ( q , t ) | 2 ⟩ , where I ̂ ( q , t ) is the spatial Fourier transform of the image intensity distribution I(x,Δt) at time t and the symbol ⟨⋅⟩ indicates a time average over the initial time t. The image structure function is connected to the collective ISF f(q,Δt) through the relation: D q , Δ t = A q 1 − R f q , Δ t + B q , (1) where R [ ⋅ ] is the real part, B(q) is an additive term accounting for the noise in the detection chain, and A(q) accounts for the static scattering amplitude [20, 22]. In most cases, including this study, ISF is a real quantity, that is, R [ f ( q , Δ t ) ] = f ( q , Δ t ) . In addition, if the dynamics is isotropic, the image structure function is expected to exhibit circular symmetry, and one can consider its azimuthal average D(q,Δt), where q = |q|. Reliably extracting f(q,Δt) from D(q,Δt) requires properly estimating A(q) and B(q). In principle, one could use the fact that f(q,Δt → 0) → 1 and obtain the camera noise contribution B(q) as the limit for Δt → 0 of D(q,Δt). In practice, this procedure is prevented by the fact that Δt is bounded from below by the inverse of the frame rate. We, thus, fit a polynomial function of the second degree to D(q,Δt) over a narrow time interval including the smallest accessible lag-times (typically the first five) and estimate B(q) as the intercept of the best fitting curve [53]. Once B(q) is estimated, we could, in principle, obtain A(q) as the long-time limit of D(q,Δt) − B(q) since f(q,Δt → ∞) → 0. Unfortunately, the complete decorrelation of f(q,Δt) cannot be observed for all q -values within the acquisition time interval. To overcome this limitation, we used the procedure introduced in Ref. [54]. We first compute a background image I ₀(x) as the minimum projection of each pixel of the entire stack in time I ₀(x) = min{I(x,t)} _t . Considering the lowest intensity value for each pixel during the acquisition time allows to efficiently remove the contribution of the rapidly moving small particles and hence a fairly accurate reconstruction of the background intensity, which in the case of glassy samples is mostly given by the auto-fluorescence of the large spheres (see Figure 5). We then estimate A(q) as the time-averaged power-spectrum of the individual images after the subtraction of the background image A ( q ) = ⟨ | I ̂ ( q , t ) − I 0 ̂ ( q ) | 2 ⟩ − B ( q ) . Once A(q) and B(q) are computed, we obtain the collective ISF f(q,Δt) by inverting Eq. 1.

As shown explicitly in Section 3.4, the collective ISF f( q ,Δt) can be written as the sum of its self-part f _self(q,Δt) and its distinct part f _dist(q,Δt). In particular, f _self(q,Δt) coincides with the Fourier-transform of the self-van Hove function, that is, the PDF P(Δx,Δt) of particle displacements [55] f self q , Δ t = ∫ d 2 Δ x P Δ x , Δ t e − j q ⋅ Δ x . (2) As long as P(Δx,Δt) is isotropic, f _self(q,Δt) depends only on q = |q| and can be expanded in cumulants for q → 0 [28, 56]. − log f self q , Δ t = 〈 | Δ x Δ t | 2 〉 4 q 2 + 1 2 α 2 Δ t 〈 | Δ x Δ t | 2 〉 4 q 2 2 + ⋯ , (3) where MSD(Δt) = ⟨|Δx(Δt)|²⟩ is the 2D mean square displacement and α ₂ (Δt) is the 2D non-Gaussian parameter α 2 Δ t = 1 2 〈 | Δ x Δ t | 4 〉 〈 | Δ x Δ t | 2 〉 2 − 1 . (4) Inspection of Eq. 3 suggests a strategy to simultaneously access MSD and the non-Gaussian parameter by focusing on the low-q behavior of the self ISF for a fixed time delay Δt, which we implement in our ConDDM analysis. We consider only the q-values and time delays Δt, for which f _self(q,Δt) has relaxed by less than 25% (see Figure 1). For each time delay Δt, this condition identifies as an upper bound q _u for q. The lower bound q _c ≃ 0.9 μm⁻¹ is a cutoff wave vector which is introduced to avoid the artifacts related to the confocal optical sectioning, as described in detail in Section 3.1. If the interval [q _c,q _u] contains at least fifteen data points, we fit a function of the form a + cq ² + dq ⁴ to log(f _self(q,Δt)), obtaining the parameters a, c, and d of the best fitting curve. Finally, we estimate the 2D MSD and non-Gaussian parameter for the considered time delay as MSD = −4c and α ₂ = 2d/c ², respectively. This procedure can be considered a generalization of the one proposed in Ref. [5].

We first consider a one-component sample, where only the small particles (of diameter σ _s = 0.56 μm) are present, suspended at a very low volume fraction (ϕ _s < 0.001) in a density-matching solvent. In this condition, the particles essentially do not interact, and each of them performs an independent Brownian motion. Importantly, in the absence of interactions, one has f(q,Δt) = f _self (q,Δt). Supplementary Movie SM1 shows a 10-s-long portion of a longer image acquisition: the fluorescent signal of a single particle is too low compared to the background noise to enable reliable tracking and subsequent reconstruction of the particle trajectories in real space; in contrast, the excellent stability of the fluorescence signal over time and the absence of appreciable bleaching allow for a quantitative analysis with ConDDM [20, 22, 36].

Representative normalized ISFs obtained from the ConDDM analysis are shown in Figure 2A. Rescaling the horizontal axis with q ² (Figure 2B) shows that a simple exponential model f(q,Δt) = e ^−Γ(q)Δt provides an excellent fit to the data over a wide range of wave vectors. In particular, the obtained relaxation rate Γ(q) displays a quadratic dependence on q for q ≳ 0.9 μm⁻¹, whereas for lower values of q, it deviates significantly from this behavior and tends to saturate to a constant value Γ_c ≃ 0.1 s⁻¹. As discussed in Ref. [36], this is due to the fact that the thickness l _c of the optical section of the confocal microscope sets a characteristic wave vector q _c = 1/l _c below which, similar to fluorescence correlation spectroscopy, the relaxation of the ISF is dominated by the fluctuation in the number of particles within the optical section; such relaxation occurs with a q-independent characteristic rate Γ c ≃ D 0 / l c 2 , corresponding to the inverse of the diffusion time across the optical section. Since, at least as far as the particle motion is isotropic, the value of q _c (≃ 1 μm⁻¹ in our experiments) is determined only by the properties of the optical system, we will restrict our attention to a q-range that is bounded from below by q _c.

In this q-range, the experimentally determined relaxation rate Γ(q) can be fitted by a homogenous quadratic function D _0,a q ², which provides the estimate D _0,a = 0.28 ± 0.03 μm²/s for the diffusion coefficient of the particles (Figure 2B, inset). This value can be compared with the result obtained from a linear fit 4D _0,bΔt of the MSD (see Methods), which provides the estimate D _0,b = 0.29 ± 0.02 μm²/s (Figure 2C). The non-Gaussian parameter α ₂(Δt), shown in Figure 2D, points to essentially Gaussian dynamics with deviations occurring only for very small time delays. Such deviations are not to be interpreted as non-Gaussian dynamics as they probably originate from the large uncertainty associated with the fact that α ₂(Δt) is obtained as the ratio of two quantities that become indistinguishably small for Δt → 0 (see Section 2.4).

Once the dynamics of a dilute suspension of small particles diffusing freely is known, we proceed to study these small particles in a dense matrix of larger particles (of diameter σ _l ≃ 5σ _s and volume fraction ϕ = 0.60, see Supplementary Movie SM2). The volume fraction of the small particles is ϕ _s = 0.006. The ISFs obtained with ConDDM at different q-values are shown in Figure 3A. Some curves display a nonzero plateau for t → ∞, which suggests non-ergodicity arising from the inability of the system to fully explore the configuration space during the observation time. As proposed in Ref. [37], this can be easily understood as a consequence of the fact that the small particles can freely move only in the voids left by the large particles, which are almost completely frozen in their positions. Here, we substantiate this picture with a simple model (see Section 3.4) according to which the expression f q , Δ t = 1 − α q f self q , Δ t + α q (5) holds. α(q) is a time-independent non-ergodicity parameter, which reflects the matrix structure (see Eq. 21 for an approximate analytical expression of α(q)). Eq. 5 shows that 1) the collective ISF f (q,Δt) is fully determined by its self-part f _self (q,Δt) and 2) the presence of a kinetically arrested matrix introduces a q-dependent plateau α(q).

The situation in which small particles can diffuse in the void space of a slowly evolving glassy matrix is compatible with the diffusing diffusivity scenario, in which particles can change their diffusivity by exploring different local environment and whose analytical predictions [16, 18, 19] can be cast in reciprocal space for a diffusion coefficient D that fluctuates with a characteristic correlation time τ ₀. By introducing Γ(q) = Dq ² and β = 1 + 2 Γ ( q ) τ 0 , the self-ISF can be expressed as f self q , Δ t = 4 β e − β − 1 Δ t / τ 0 1 + β 2 − 1 − β 2 e − 2 β Δ t / τ 0 n / 2 (6) where n is a free parameter defining the dimensionality of the stochastic process of D. The mean square displacement MSD can be then evaluated as follows M S D Δ t = − 2 ∂ 2 f self q , Δ t ∂ q 2 q = 0 = 2 n D Δ t , (7) showing a linear dependence with time Δt. Inserting Eq. 6 with n = 2 in Eq. 5 provides a model that can be fit to our experimental collective ISF. The fit is performed at fixed q, the fitting parameters being Γ, τ ₀, and α. The results show that the ISFs are well-described by the model (Supplementary Figure S2A) but point out that our experiments were too short to capture the characteristic time τ ₀ of the fluctuating diffusive coefficient (Supplementary Figure S2B), which is expected to be a q-independent quantity.

For our particular experiments, it is, thus, useful to focus on the two limiting regimes (for n = 2):

• Γ(q)τ ₀ ≫ 1, in which the density fluctuations relax faster than the characteristic fluctuation time of the diffusion coefficient, and one has

We note that the second regime corresponds to ordinary diffusion, where the PDF of 2D particle displacements corresponding to the exponentially relaxing ISF in Eq. 9 is a Gaussian function P Δ x , Δ t = 1 4 π D Δ t exp − | Δ x | 2 4 D Δ t . (10) On the other hand, the PDF corresponding to the first regime, which can be obtained as the 2D Fourier transform of the self-ISF in Eq. 8, is P Δ x , Δ t = 1 2 π D Δ t K 0 | Δ x | D Δ t , (11) where K ₀ is the modified Bessel function of the second kind of order 0 [57]. As it can be appreciated from Supplementary Figure S1 in ESI, where representative curves are shown, in this case, the PDF displays strongly non-Gaussian, exponential-like tails.

The results shown in ESI Supplementary Figure S2, suggest adopting Eq. 8 as a model for the self-ISF. Such a model is in excellent agreement with our experiments (Figure 3). In particular, the collective ISFs (Figure 3A) collapse well on a single master curve when converted to the self-ISF and plotted as a function of the rescaled time q ²Δt (Figure 3B). The inset of Figure 3B shows the extracted q-dependent relaxation rate, which can be fitted by Γ(q) = D _60,a q ² to provide a diffusion coefficient D _60,a = (11 ± 1) × 10^–2 μm²/s. This diffusive scaling along with the Lorentzian-like fitting model f _self(q,Δt) = [1 + Γ(q)Δt]⁻¹ indicates that although the small particles exhibit Brownian diffusion, the probability distribution of the displacements is no longer Gaussian.

To further support our findings, we show in Figures 3C,D the particle MSD and the non-Gaussian parameter extracted from the self-ISFs, respectively. The MSD displays a linear scaling with a diffusion coefficient D _60,b = (10.5 ± 0.5) × 10^–2 μm²/s, which is in agreement with D _60,a. Moreover, the non-Gaussian parameter α ₂ now is significantly different from zero and hence points to non-Gaussian dynamics. α ₂ was obtained, according to Eq. 3, by considering the first two terms of the cumulant expansion of the ISF. In principle, the expansion could be extended to the calculation of cumulants (or moments) of arbitrary order, which can then be used to retrieve the full PDF of particle displacements. In practice, the noise on the data rapidly makes this procedure numerically extremely unstable, rendering the determination of higher-order cumulants very challenging. An alternative strategy to retrieving direct-space information consists in adopting a suitable analytical model for the ISFs (like the ones in Eqs 6–8), determining the free parameters via a fitting procedure, and estimating the PDF as the (analytically or numerically evaluated) Fourier transform of this model function. According to the aforementioned discussion, the application of this procedure to our data results in a prediction for the PDF of particle displacements in the form of Eq. 11, which is characterized by exponential-like tails (see Supplementary Figure S1). This result is in slight disagreement with the findings of Ref. [10], where the motion of diluted small particles diffusing in suspensions of larger spheres is investigated. In that case, the PFDs of particle displacements were found to show consistent deviations from Gaussianity and robust diffusive scaling but without a clear indication of an exponential tail. On the one hand, this discrepancy can originate from the fact that in Ref. [10] the large particles were always in a fluid-like state, and thus a clear-cut separation of the time scales associated with tracer diffusion and structural rearrangement of the matrix was not present. Moreover, compared to Ref. [10], in our experiments, the unbalance between the size of the small and large particles is less pronounced. As a consequence, in the present case, the small particles move within the relatively narrower “voids” formed by the large spheres, and this is expected to enhance hydrodynamic coupling and sensitivity to the local environment.

When the volume fraction ϕ of the large particles is increased, we expect the dynamics of the small particles to change. This expectation is only partially confirmed by a ConDDM analysis of the movies acquired for ϕ = 0.60, ϕ = 0.61, and ϕ = 0.62 (Figure 4). The relative volume fraction of the small particles is x _s = ϕ _s/ϕ = 0.01 for all cases, meaning that ϕ _s is slightly different for the three samples. Most of the features outlined in the previous section for ϕ = 0.60 are also observed for the higher volume fractions. In particular, the self-ISFs for ϕ = 0.61 and 0.62 show a scaling similar to the one of ϕ = 0.60 when plotted as a function of the rescaled time q ²Δt (Figure 4A). The relaxation rates obtained from the best fitting curves using the model in Eqs 5, 8 for the three different concentrations are shown in Figure 4B, in which we observe a modest change of the diffusion coefficient for this narrow range of volume fractions. This result is confirmed by MSDs (Figure 4D) and non-Gaussian parameters (Figure 4E), which also show minor changes as a function of the volume fraction.

The most evident difference between the experiments performed at different ϕ is found in the magnitude of α(q) which for ϕ = 0.61 and ϕ = 0.62 is significantly larger than for ϕ = 0.60 (Figure 4C). The fact that α(q) is significantly different from zero even in the high-q limit implies that the small particles are not free to diffuse, even not over small distances. As discussed in more detail in Section 3.4, this suggests the presence of an immobile fraction (of about 6% and 10% for ϕ = 0.61 and ϕ = 0.62, respectively) of small particles stuck in the voids of the matrix. The appearance of an immobile fraction of particles upon increasing the volume fraction can be qualitatively appreciated in Figure 5, where we show the minimum projection of the image sequence for ϕ = [0.60, 0.61, 0.62] and in which the stuck small particles appear as bright spots. This effect was not observed in previously published molecular dynamics simulations [37], which indicated that small-particle localization is not to be expected for size ratio δ < δ _c ≃ 0.35, for all the investigated volume fractions. Therefore, the observed increase in the number of stuck particles may be attributed to system aging or particle aggregation, which would both prevent particles from freely exploring the matrix.

As demonstrated in Refs. [37, 39], the binary mixtures considered in Section 3.2 are in a “single-glass” state, that is, the dynamics of the large particles is kinetically arrested, whereas the small particles are not arrested and remain free to diffuse over arbitrarily large distances. Therefore, the self-dynamics of the small particles is expected to be ergodic. In this condition, the glassy nature of the matrix formed by the large particles induces a decoupling between the self and the collective ISF of the small particles: while f _self(q,Δt) can completely relax to zero for Δt → ∞; the f(q,Δt) displays a nonzero plateau, usually referred to as the non-ergodicity factor α(q), which reflects the spatial modulations in the collective dynamics imposed by the partially frozen structure of the large spheres [39].

While a detailed theoretical discussion of these aspects is beyond the scope of this study (see, e.g., Refs. [39, 49] for a thorough treatment), we present here a simplified calculation aimed at 1) providing a simple physical picture for the emergence of a finite non-ergodicity factor, 2) providing a prediction for the q-dependent non-ergodicity factor to be compared, at least qualitatively, to the experimental results, and 3) showing that, at least for highly diluted tracers, the self-part of the ISF can be reliably extracted from the collective ISF.

We start with the definition of the (unnormalized) collective ISF F(q,Δt) of the small particles F q , Δ t = 1 N s ⟨ ρ s ̂ q , t + Δ t ρ s ̂ q , t * ⟩ , (12) where ρ s ( x , t ) = ∑ n δ ( x − x n ( s ) ) is the number density of small particles and N _s = ∫d x ρ _s(x,t) is the total number of particles within the considered volume. In general, the ISF can be decomposed into its self and its distinct part [55] F q , Δ t = f self q , Δ t + f dist q , Δ t , (13) where f self q , Δ t = 1 N s ⟨ ∑ n e − j q ⋅ x n s t + Δ t − x n s t ⟩ , (14) and f dist q , Δ t = 1 N s ⟨ ∑ n ≠ m e − j q ⋅ x n ( s ) t + Δ t − x m s t ⟩ . (15) Under the hypotheses that 1) the motion of the small particles takes place in the presence of a completely frozen matrix of larger particles and 2) both the volume fraction ϕ _s of the small particles and the ratio δ = σ _s/σ _l between the diameters of the small and large particles are very small, the space accessible to the small particles essentially coincides with the portion left free by the large ones; since the effective volume fraction ϕ _s/(1 − ϕ _l) of the small particles within the voids is very small, positional correlations between distinct small particles are negligible. The time-averaged density distribution ρ ̄ s ( x ) = ⟨ ρ s ( x , t ) ⟩ of the small particles can, thus, be assumed to be uniform over the whole accessible region ρ ̄ s x = N s V 0 1 − ϕ l 1 − p l x * ∑ n δ x − x n l . (16) In this expression, V ₀ indicates the total volume, p _l(x) is the characteristic function of a sphere of diameter σ _l, the symbol * indicates the convolution operation, and the prefactor N _s/V ₀ (1 − ϕ _l) is determined by imposing the normalization condition ∫ d x ρ ̄ s ( x ) = N s . Since the motion of two small particles can be assumed to be completely uncorrelated, the distinct part of the ISF (Eq. 15) is expected to be time-independent and, up to an additive term of order 1/N _s , it is given by f dist q ∼ 1 N s 〈 ∑ n e − j q ⋅ x n s t + Δ t 〉 〈 ∑ m e + j q ⋅ x m s t 〉 = 1 N s ρ ̄ ̂ s q ρ ̄ ̂ s q * = 1 N s | ρ ̄ ̂ s q | 2 . (17) Using Eq. 16, the power spectrum | ρ ̄ ̂ s ( q ) | 2 of the equilibrium density distribution of the small particles can be expressed in terms of the form factor P l ( q ) = ( 1 / V l 2 ) | p ̂ l ( q ) | 2 and of the static structure factor S l ( q ) = ( 1 / N l ) | ρ ̂ l ( q ) | 2 of the large particles as follows | ρ ̄ ̂ s q | 2 = N s V 0 1 − ϕ l 2 V l 2 P l q N l S l q . (18) Plugging this expression in Eq. 17 we finally get f dist q = x s δ 3 ϕ l 1 − ϕ l 2 P l q S l q , (19) where δ is the previously introduced ratio between the radius of the small and large particles (δ = 0.18 in the present case) and x _s = ϕ _s /ϕ is the previously introduced relative fraction of the small particles (x _s = 0.01 in our case). Combining Eq. 19 with Eq. 13, we obtain the following expression for the normalized ISF f(q,Δt) = F(q,Δt)/F(q, 0) f q , Δ t = 1 − α q f self q , Δ t + α q , (20) where we have introduced the “non-ergodicity factor” α q = 1 + δ 3 x s 1 − ϕ l ϕ l 2 1 P l q S l q − 1 . (21) Eqs 20, 21 show that even in the absence of explicit interactions between the moving particles, the presence of fixed obstacles introduces a q-dependent, non-ergodicity factor in the normalized ISF.

In order to compare the theoretical expression in Eq. 21 with the experimentally determined α(q), we need an estimate for the static structure factor S _l(q) of the large particles. To this end, we have used the code described in Ref. [58] and freely available at https://github.com/VasiliBaranov/packing-generation to generate different configurations of randomly closed packed hard-spheres from which we evaluated the three-dimensional (3D) static structure factor. As a consistency check, we have also calculated from the same data a two-dimensional (2D) structure factor, obtained by considering the centers of the particles included in slices of height Δz = 1 μm, similar to the thickness l _c of the optical section of our confocal microscope. This function can be directly compared to the experimentally determined 2D structure-function of the large spheres, which was obtained by performing the minimum projection of the image sequence (see Figure 6A) and by locating the positions of the centers of the large particles using the particle-tracking software freely available at https://github.com/dsseara/microrheology [59]. As it can be appreciated from Figure 6B, the experimental and the simulated 2D static structure functions are in excellent agreement, confirming the reliability of our simulation. In Figure 6C, the theoretical non-ergodicity factor calculated inserting in Eq. 21 the simulated 3D S _l(q) and the theoretical form factor P _l(q) for a sphere of radius σ _l is compared with the one obtained by fitting Eq. 5 to the collective ISF. A different estimate for α(q), obtained by evaluating the ISFs for the large delay time Δt = 10⁴ s (as carried out in Ref. [39]), is also reported. We observe that despite the rather crude approximations leading to Eq. 21, the overall agreement between the theoretical prediction and experimental data is qualitatively good, and only for qσ _l ≲ 4 (corresponding to q ≲ 1.3 μm⁻¹), a more systematic deviation is observed. This could be due to the fact that while the simulated α(q) is computed considering a fully three-dimensional matrix, for low q, the experimentally estimated α(q) could be affected by the finite thickness of the confocal optical section.

Based on Eq. 21, we expect that α(q → ∞) → 0 because P _l(q → ∞) → 0. In Section 3.3, however, we have seen that, at least for the samples with ϕ > 0.60, α(q) does not decay to zero. This was attributed to the presence of a small fraction ϕ _s,i of immobile small particles. The simplest way to account for these particles is to introduce the change α(q) → α(q) + (1 − α(q))ϕ _s,i. This change captures the observed behavior of α(q) for large q (see Figure 4C) but not for wave vectors smaller than the structural peak of the large particles qσ _l < 6, suggesting that further contributions should be incorporated in Eq. 16. This case illustrates the need for extensions of Eq. 5, which is strictly valid only in the idealized case of infinitely diluted and infinitely small particles diffusing in a perfectly frozen matrix of immobile large particles. For example, if the finite size of the particles becomes relevant, the expression given in Eq. 16 must be modified in order to account for the additional excluded volume. Similar changes need to be introduced whenever the positional fluctuations of the large particles or the interactions between the small particles need to be taken into account.