The following description of the optical characteristics is similar to ALICE 2 facility used on MIR experiments (see (Lecoutre et al., 2014; Garrabos et al., 2015; Lecoutre et al., 2015; Mota et al., 2015; Nikolayev et al., 2015; Garrabos et al., 2016; Durieux et al., 2017; Garrabos et al., 2018) and references therein). ALI insert has a modular optical design with a “source optical box” containing the laser, different filters, and photodiodes, the “thermostat box” that includes the sample cell unit, and a “collecting optical box” that contains the CCD and additional photodiodes (Lecoutre et al., 2009). The complete optical scheme is detailed in (Marcout et al., 1994), and the optical performances are precisely analyzed for other experiments, such as turbidity measurements (Lecoutre et al., 2009).

The fluid sample cell is visualized through light transmission normal to the windows using LED illumination with a spectrum centered around 660 nm. ALI instrumentation works with a wide field of view of 10 × 10 mm² object image at 30 µm resolution. It also records small field of view images of approximately 1 × 1 mm² object image from the center of the DOC with a microscope objective within Δx = 3.5 µm resolution. All images analyzed in this study were of the small field-of-view type obtained with the optical microscope at the focal plane in the middle of fluid layer, centered on the optical axis of the fluid sample. We analyzed N _img = 1,922 recorded with a sampling rate of f _s = 11.5 frames per second (the time interval between two successive images is t _s = 1/f _s ≈ 87 ms). All images I ( r ⃗ , t ) are digitized with a spatial vector given by r ⃗ = Δ x ( n x , n y ) , where n _x and n _y are integer numbers comprised between 1 and the image size M = 1,024 pixels.

Figure 2A1 shows an original image and its corresponding power spectrum in the Fourier space (q _x , q _y ) = q _min(m _x , m _y ), with m _x , m _y integers comprised between − M _w /2 and M _w /2, where M _w is the Fourier window size, and q _min = 2π/(MΔx) = 2π/W in which W ≈ 1 mm is the side of the small field-of-view image (see Figure 2A2). In general, the Fourier window size M _w is a power of two. In this study, the Fourier window size M _w was the same as the image size M = 1,024. The power spectrum has azimuthal symmetry, as illustrated by the “ring” surrounding the DC component at the center of Figure 2A2. Although all spectra have symmetry, the azimuthal symmetry is not perfect for all spectra. For some IMFs (see Figures 2C2–D2), the vertical and horizontal anisotropy of the objects present in these IMFs leads to a rounded square shape of the power spectra for some wavenumbers. Due to the azimuthal symmetry of the power spectra (see Figures 2A2–E2), only the azimuthal averages of the power spectra versus the magnitude of the wavenumber q = q x 2 + q y 2 = q min m x 2 + m y 2 were considered for the DDM analysis (see the red continuous line insets in Figures 2A2–E2). The azimuthal average is calculated over thin rings with nearly the same wave vector modulus q = q x 2 + q y 2 . Each ring has a width of 1 unit in the discrete Fourier space. Since the DC component of the power spectrum near q = 0 is orders of magnitude larger than the rest of the power spectrum, and because it does not contribute any valuable information to the subsequent analysis, we removed it from the insets. The DC component leakage into adjacent wavenumbers and significantly affects the residual image’s power spectrum shown in Figure 2E2. The wavenumber resolution in Fourier space is q _min = 2π/W ≈ 62.8 cm⁻¹, where W ≈ 1 mm is the side of the small field-of-view image. At small wave vectors, the azimuthal averaging statistic is relatively poor, but for large wavenumbers, the number of independent samples in a single wavevector ring q = q x 2 + q y 2 = q min m x 2 + m y 2 increases proportionally to π × m, in which m is the number of the channel which varies from 1 to M w / 2 , where M _w is the window size of Fourier transform. Since our images are square, m = m _x = m _y . At the maximum wave vector m ≈ M w / 2 ≈ 724, which results in over 2,000 statistically independent samples for the respective value of the wavevector q. The wavevector can be expressed either in Fourier space integers 1 ≤ m ≤ M _w for an M _w -point Fourier transform or in corresponding cm units by using the Fourier resolution formula q = m × q _min.

Without repeating the details of the almost two decades old Differential Dynamic Microscopy (DDM) method (see Supplementary Section S5.1 for a brief description and references (Cerbino and Trappe, 2008; Giavazzi et al., 2009)), we only mention novel results based on this revolutionary approach to the dynamic light scattering technique. The DDM is a spectral, Fourier-based image analysis method. The two-dimensional intensity images I ( r ⃗ , t ) are first normalized by the average image intensity to account for light source intensity fluctuations during the measurement, i.e, i ( r ⃗ , t ) = I ( r ⃗ , t ) / I ( r ⃗ , t ) ̄ r ⃗ . The fundamental quantity in DDM is the normalized image difference δ i ( r ⃗ , Δ t , t ) = i ( r ⃗ , t + Δ t ) − i ( r ⃗ , t ) , where Δt is a delay time between frames (Cerbino and Trappe, 2008). The advantage of considering the image difference is a significant reduction in static scattering coming from dust, scratches on the sample cell, or other static imperfections of the setup (Kessler et al., 2020).

The DDM also uses a moving average of the power spectrum of the fluctuating image δ i ( q ⃗ , Δ t , t ) = | F F T [ δ i ( r ⃗ , Δ t , t ) ] | 2 over different reference times t and the same delay Δt to reduce optical noise. Such a moving average is called the image structure function, i.e., C m ( q ⃗ , Δ t ) = δ i ( q ⃗ , Δ t , t ) ̄ t . The reason for the temporal average is to increase “the statistical accuracy of the data” (Cerbino et al., 2017). The number of reference times over which the average is computed N _avg is generally smaller than the correlation time of fluctuations such that fast relaxation processes are not averaged out. In practice, one uses a recursive selection process in which we start with a high enough number of averaged images to reduce the noise and obtain good statistics, e.g., N _avg = 32. Once the correlation time of fluctuations is determined, one can change N _avg and repeat the estimation of the correlation time. By varying the delay time Δt during image fluctuation calculations, one can extract the power spectrum both over the wavenumbers q and elapsed time Δt between frames. The dependence of the power spectrum on the elapsed time Δt between frames contains information on the correlation time of fluctuations. The two-dimensional image structure function probes the sample dynamics in different directions in the q ⃗ plane. Whenever the image structure function bears a circular symmetry, the azimuthal averaging, i.e., C m q , Δ t = C m q ⃗ , Δ t ̄ q ⃗ is used to obtain the one-dimensional image structure function. The one-dimensional image structure function C _m (q, Δt) “rises until saturating when the images are totally decorrelated” as seen in Figure 3A (Lu et al., 2012).

Examples of azimuthal averages of power spectra, such as those shown with the continuous red line in Figures 2A2–E2 insets, versus the delay time Δt for fixed wavenumbers q are shown in Figure 3A. The (one-dimensional) time-dependent structure functions C _m (q, Δt) describes how the spectral power changes with the delay time Δt between images for a fixed wavenumber q, and it is given by (Cerbino and Trappe, 2008; Giavazzi et al., 2009; Lu et al., 2012): C m q , Δ t = 2 A q 1 − G q , Δ t + B q , (2) where A(q) is an amplitude term that contains information about the static scattering from the sample and the optical system’s transfer function, B(q) is the background contribution to the time-dependent structure functions, and G(q, Δt) is “equivalent” (Lu et al., 2012) of the normalized ISF from the traditional dynamic light scattering experiments (see Figure 3B). For particular optical setups, such as Schlieren and Shadowgraph, the analytical form of the optical transfer function T(q) is known (see (Trainoff and Cannell, 2002; Croccolo et al., 2006; Cerbino and Trappe, 2008; Giavazzi et al., 2009; Oprisan and Leilani Payne, 2013; Croccolo et al., 2016a)).

The effect of DDM moving average on the structure factor and fluctuation correlation time. As mentioned above briefly, DDM relies on the temporal moving average of power spectra to reduce as much as possible the ubiquitous image noise (see (Vailati and Giglio, 1998; Croccolo, 2006; Croccolo et al., 2006; Croccolo et al., 2007; Cerbino and Vailati, 2009; Giavazzi et al., 2009; Vailati et al., 2011; Ortiz de Zárate et al., 2014; Cerbino et al., 2015; Croccolo et al., 2016a; Giavazzi et al., 2016a; Bataller et al., 2016; Croccolo et al., 2016b; Giavazzi et al., 2016b; Cerbino et al., 2017) and references therein) and to compute the structure factor and correlation times over averaged fluctuation images. The number of images over which denoising average is performed depends on the particular conditions of each experiment. As mentioned above, the general recommendations are (see (Vailati and Giglio, 1998; Croccolo, 2006; Croccolo et al., 2006; Croccolo et al., 2007; Cerbino and Vailati, 2009; Giavazzi et al., 2009; Vailati et al., 2011; Ortiz de Zárate et al., 2014; Cerbino et al., 2015; Croccolo et al., 2016a; Giavazzi et al., 2016a; Bataller et al., 2016; Croccolo et al., 2016b; Giavazzi et al., 2016b; Cerbino et al., 2017) and references therein) that the number of images averaged needs to be small enough to avoid averaging out fast-changing processes or completely de-correlating the images in a batch, and, at the same time, the number of images in an averaging batch must be large enough to provide good noise filtering. We based our selection of the number of images in an averaging batch on the delay time between frames Δt it takes for the time-dependent structure functions C _m (q, Δt) to begin saturating. As shown in Figure 3A, the time-dependent structure functions C _m (q, Δt) starts saturating around 1–3 s, corresponding to 20–70 consecutive images. In other words, the correlation between fluctuations becomes negligible after 20–70 successive images. In this study, we performed the N _avg = 32-image temporal moving average of power spectra since it was on the lower end of the above saturation range of time-dependent structure functions C _m (q, Δt). Such a value for the number of images in an averaging batch still captures short correlation times down to 1×10^–2 ms while producing decent image denoising (see Figure 3D). We found that a lower size of the moving average did not change the amplitude factor A(q) fitting nor the correlation time estimate (not shown). On the upper limit side, a N _avg = 64 images moving average batch could also work, but that would be towards the upper limit of the saturation time for the time-dependent structure functions C _m (q, Δt) (see Figure 3A) and has the potential of filtering out shorter saturation times. Additionally, such a large averaging batch size leads to fewer data points at low correlation times in Figure 3D and, therefore, a less precise estimation of the effective diffusion coefficient D. Since time-dependent structure functions C _m (q, Δt) were averaged over N _avg = 32 images, the image index and the corresponding recording time reported on all figures referring to C _m (q, Δt) or any quantity derived from it, such as the amplitude A(q), the ISF G(q, Δt), the correlation time of fluctuations τ(q), and the diffusion coefficient D, are for the first image in the batch.

DDM data correction to account for the optical transfer function T(q). The amplitude A(q) = S(q) × T(q) from Eq. 2 contains information on both the structure factor S(q) and the transfer function T(q) of the optical setup. As shown by Cerbino and Vailati (Cerbino and Vailati, 2009), “in the limit of coherently illuminated” optical field and “weakly scattering objects,” the transfer function is T ( q ) = sin 2 ( q 2 ⁡ z λ 4 π ) , which has increasingly fast oscillations as the wavevector q increases. As noted by Cerbino and Vailati (Cerbino and Vailati, 2009), “in practice, for sufficiently large wave vectors q, the function T(q) may lose its oscillatory character and in some cases there is no need for data correction.” According to Cerbino and Vailati (Cerbino and Vailati, 2009), for wave vectors larger than the “crossover wave vector” q cross = min π 2 e cyl λ , π L sensor z λ , “there is no need for data correction” due to the transfer function T(q). In our experiment, the sample thickness is e _cyl = 4.115 mm, the sensor size L _sensor is given by sensor’s resolution Δx = 3.5 µm and the number M = 1,024 of pixels, i.e., L _sensor = Δx × M = 3.584 mm. The wavelength of the LED used was λ = 660 nm, and the estimated distance z for the microscopic view images is the order of 1 mm. The visualization distance z is measured between the sample and the plane imaged onto the sensor (Croccolo, 2006). With the above data, the crossover wavenumber due to sample thickness is q _cross ≈ 850 cm⁻¹. Note that this estimation of the crossover wave number q _cross due to the sample’s thickness still holds for all values of z up to 20 cm, which is well within the range of our microscope setup. For wavenumbers q larger than the crossover value q _cross ≈ 850 cm⁻¹, we used the amplitude A(q) and the structure factor S(q) interchangeably. To be on the conservative side, we only estimated the slope of the amplitude A(q) for wavenumbers larger than q _cross to ensure that there is no need for T(q) correction (see Figure 3). Based on the above estimate of the crossover wavenumber, the amplitude A(q) was also fitted with Eq. 3 in the large wavenumber range where there is no need for data correction (Cerbino and Vailati, 2009).

Another potential issue when estimating the slope and fitting the amplitude A(q) with Eq. 3 is the temporal moving average performed by DDM. In general, image averaging is equivalent to a lowpass filtering that attenuates high (spatial) frequencies and could change the slope of the structure factor or slightly shifts it along the wavenumber axis. We checked that the slope of the amplitude A(q) shown in Figure 3C with a denoising DDM average of N _avg = 32 did not change compared to the slope of A(q) for individual images. The variance for the slopes obtained without averaging was larger than when averaging over batches of N _avg = 32 images.

The experimental structure factor S(q) obtained by fitting the one-dimensional image structure function with Eq. 2 can be approximated with S q = S 0 1 + q ξ n = S 0 1 + q / q corr n , (3) where S ₀ is the structure factor value at very low (q → 0) wavenumbers, ξ is the correlation distance of fluctuations, q _corr = 1/ξ is the correlation wavenumber given in Table 1, and n is a power-law exponent usually close to 2 for critical fluctuations (see also Figure 3A in Ref. (Oprisan et al., 2012; Oprisan et al., 2021a) for structure factors of equilibrium fluctuations near the critical point in other systems). The curve described in Eq. 3 is shown with a green dashed line in Figure 3C on top of the amplitude A(q) experimental data. Finally, we notice that the slope of the amplitude A(q) shown with a continuous black line changes with the wavenumber range considered. For example, near the peak of A(q), the slope, i.e., the power-law exponent n from Eq. 3, approaches zero, whereas, at large wavenumbers, it approaches −4. The existence of multiple power-law exponents for the structure factor was attributed to the fractal nature of the observed dynamics using scattered light (Sorensen, 2001). One should not give too much weight to data points at large wavenumbers as they represent data at the lower end of CCD camera resolution, which is around q _max = 2π/3.5 µm ≈ 18 ,000 cm⁻¹ in our experiment. Other factors that limit the wavenumber range are poor signal-to-noise ratio, inadequate sampling frequency, and recording duration. For a more detailed analysis of different limiting factors, see (Giavazzi et al., 2017).

The saturation of time-dependent structure functions C _m (q, Δt) shown in Figure 3A and the dependence of the amplitude A(q) on the wavenumber q shown in Figure 3C provide insight into the characteristic relaxation time of the density fluctuations (Oprisan et al., 2012). The corresponding ISF shown in Figure 3B exhibits an exponential relaxation for delay times shorter than 3 s. The slope of the ISF shown in Figure 3B gives the inverse of the correlation time of fluctuations shown in Figure 3D. At large wavenumbers, the log-log plot of correlation time is always linear with a slope of −2 (see the continuous black line in Figure 3D).

Theoretically, the relaxation time τ of critical fluctuations near critical temperature T _c should obey Kawasaki (Kawasaki, 1970) formula: τ − 1 ≈ D q 2 , (4) when D is the thermal diffusivity coefficient. Such results are exemplified in Figure 4B, where the dashed green curve corresponds to a smooth spline interpolation of the correlation time of fluctuations in the original image (see the solid black square in Figure 4B) to highlight the inflection point at critical wavenumber q _c . In this case, the straight line with a slope of −2 (see the continuous white line in Figure 3) gives the effective diffusion coefficient D.

Density fluctuations near critical points have a wide range of sizes limited only by the boundaries of the enclosing container (Beysens, 1986; Beysens et al., 1987; Beysens and Garrabos, 2000; Onuki, 2002; Barmatz et al., 2007; Midya and Das, 2017). Critical density fluctuations near critical points show fractal patterns (Guenoun et al., 1989; Schaefer et al., 1989; Antoniou et al., 1998; Antoniou et al., 2000; Oprisan et al., 2021a; Oprisan et al., 2021b). Nonequilibrium fluctuations can also lead to fractal patterns (Vailati et al., 2011). The fractal structures are usually described by a spectrum of fractal dimensions, which are power-law exponents associated with different spatial or temporal scales present in the data. How would a fluctuating image near the critical point look if we could break it into disjoint spatial scales, like decomposing white light into narrow-band, monochromatic waves? We noticed in the previous Section 3.1 that the amplitude A(q) factor in Eq. 2, which coincides with the structure factor S(q) above a crossover wavenumber that eliminates the need for optical transfer function correction of the data, scales as A(q) ∝ q ⁻² (see Figure 3C). We also noticed that the amplitude A(q) slope changes for different wavenumber ranges. Could different power-law exponents required over different wavenumber ranges reflect the multifractal nature of critical fluctuations? What would such a multiscale analysis over disjoint wavenumber ranges reveal about the behavior of the fluctuations’ correlation time? The DDM-based results obtained in Section 3.1 from the original images indicate that the correlation time decreases with the wavenumber (see Figure 3D), i.e., the correlation time increases with the characteristic size of fluctuations. What happens to the slope of the structure factor over different spatial scales? What happens to the critical wavenumber q _c that marks the transition to a different fluctuation regime? Where is the critical wavenumber q _c located on images that only capture fluctuation over a limited wavenumber range? Do the scaling laws governing each spatial scale change when decomposing fluctuations into narrow wavenumber ranges? How are the relaxation times of fluctuations at each spatial scale related to the dynamics of fluctuations in the original image? We decomposed the original images over disjoint wavenumber ranges using the BEMD method to answer such questions. We then determined the correlation time of fluctuations and critical wavenumbers q _c using DDM as described in Section 3.1.

We used BEMD to break the original image (Figure 2A1) into three IMFs (Figures 2B1–D1) and a residual background (Figure 2E1). We performed the data analysis described in the previous subsection on all IMFs. To compare side-by-side the structure factors of the original image and the three IMFs shown in Figures 2A1–D1, we normalized each of them by their respective maximum values and plotted all four curves in Figure 4A. We notice that the amplitude factor A(q) of the original image (solid black squares in Figure 4A), which is proportional to the structure factor S(q) for wavenumbers larger than q _cross = 850 cm⁻¹, shows a power-law behavior with an exponent n close to −2 (see the yellow dashed line in Figure 4A). The theory of critical fluctuations predicts a power-law exponent equal to −2 (Beysens, 1986; Beysens and Garrabos, 2000; Domb et al., 2001).

As the BEMD decomposes the original image into orthogonal IMFs by separating the spatial scales of fluctuations, it only retains spatial structures larger than a characteristic size starting from the smallest (IMF1), then intermediate (IMF2) and large (IMF3). This means that each IMF has both a lower and upper boundary on the wavenumbers that it can capture. While there is no set rule for defining the exact boundaries of a given IMF order, we used the peak q _peak of the structure factor criterion, i.e., the lower wavenumber for a given IMF is the peak of the structure factor. The upper limit of the wavenumbers for a given IMF order is the lower boundary of the precedent IMF order. For example, q _maxIMF2 = q _minIMF1 . Specifically, IMF1 only contains small-size objects from the original images, which means that its corresponding wavenumbers will start around the peak of the IMF1’s structure factor, which is q _minIMF1 ≈ 8,500 cm⁻¹ and expand up to the larges wavenumbers captured by the original images. IMF2 will only capture intermediate size objects, which means the wavenumbers over which its results are accurate cover q _minIMF2 ≈ 5,000 cm⁻¹ to q _maxIMF2 ≈ q _minIMF1 . Ideally, if the IMFs are truly orthogonal decompositions of the original image, there should be no leakage between them and, therefore, the wavenumber ranges for each IMF should be disjoint. In such an ideal case, we would expect that the smallest wavenumber in IMF1 q _minIMF1 will be close to the largest wavenumber in IMF2 q _maxIMF2 and so on. Since we have a bit of leakage between IMFs, the wavenumber ranges may overlap, most likely due to image noise. The IMF3 will only capture large-size fluctuations, which means the wavenumbers over which its results are accurate cover q _minIMF3 ≈ 3,000 cm⁻¹ to q _maxIMF3 ≈ q _minIMF2 . We note that the range of wavenumbers for IMF3 is well above the crossover wavenumber q _cross = 850 cm⁻¹ determined with the criterion from (Cerbino and Vailati, 2009).

Suppose we retain only the wavenumbers larger than the peak location of IMF1 in the structure factor of the original image (see point A1 in Figure 4A). In that case, the structure factor of the original image is determined mainly by the structure factor of IMF1 and IMF2. This is because, as seen in Figure 4A, for wavenumbers larger than those indicated by the peak A1 of IMF1, the structure factor of IMF2 is one order of magnitude smaller than IMF1, and IMF3 is at least three orders of magnitude smaller than IMF1. At the same time, the slope of the original image’s structure factor is intermediate between the slope of IMF1, which is on average −2.1 ± 0.6, and the slope of IMF2, which is on average −4.3 ± 0.3. The solid red rectangle marks the wavenumber range for IMF1 on the wavenumber axis in Figure 4A. The solid blue rectangle marks the corresponding wavenumber range for IMF2 on the wavenumber axis in Figure 4A.

Finally, at small wavenumbers, between the peak of the IMF2 (blue vertical dashed line in Figure 4A) and IMF3 (green dashed vertical line in Figure 4A), the structure factor of the original image lies between points A2 and A3 in Figure 4A. The structure factor of the original image is determined mainly by IMF2 and IMF3 as the structure factor for IMF1 is more than one order of magnitude smaller than the other two. Since the slope of the original image’s structure factor swings from positive (at low wavenumbers) to negative (at higher wavenumbers), it is a mix of IMF2, which has a positive slope, and IMF3, which has a very abrupt negative slope of −5.6 ± 1.1. The fact that the original image’s negative slope of the structure factor is not as steep as for IMF3 shows that the positive slope of IMF2 contributed significantly to the original image’s structure factor. The peak wavenumber decreases for successive IMFs, as shown in Figure 4A.

A similar picture emerges when analyzing the correlation time of fluctuation in the original image and the corresponding IMFs (see Figure 4B). We did not normalize the correlation times but retained the actual values in seconds as obtained from the intermediate scattering function (see Figure 3B) with Supplementary Appendix SEq. 7. We repeated the correlation time fitting for all IMFs following the procedure detailed in the previous section and determined the corresponding D and q _c . As we notice from Figure 4B, the peak location of the correlation time (q _c ) for IMF1 (see point B1 in Figure 4B) determines the range of wavenumbers over which it overlaps with the correlation time of the original image (see the red rectangle marked IMF1 on the wavenumber axis in Figure 4B). Similarly, the IMF2 correlation time starts overlapping the correlation time of the original image at the onset of intermediate-range wavenumber (point B2 in Figure 4B). Finally, the IMF3 wavenumber range extends from low wavenumbers marked by point B3 in Figure 4B, and the overlap does not extend farther than B1. At large wavenumbers, the correlation time for IMF3 flattens out (solid green inverted triangles in Figure 4B), suggesting that we reached the limit of maximum acceptable wavenumber for the respective IMF. Similar behavior can be observed for IMF2 (blue solid upright triangles in Figure 4B). However, as seen in Figure 4B, the flattening out of the correlation time of IMF3 (solid inverted triangles in Figure 4B) is determined by the fact that the results past q _maxIMF3 are not accurate since there is no information in IMF3 regarding objects at such wavenumbers.

In addition, the correlation time of fluctuations (see Figure 4B) shows that the maximum correlation time τ _max increases with IMF’s order. In other words, large wavenumber fluctuations, i.e., small physical size fluctuations, have short correlation times. For example, the smallest spatial scale of IMF1 has a maximum correlation time of 0.2 s, the IMF2 reaches 2 s, and the largest spatial scale of IMF3 reaches a maximum correlation time of 10 s (see Figure 4B). This relationship between the maximum correlation time and the physical size of the fluctuations confirms the fractal nature of fluctuations (see also (Oprisan et al., 2021a)). Most of the focus of scaling laws near criticality was on spatial criticality, i.e., deriving the universal laws for correlation length of fluctuations (Moldover et al., 1979; Beysens, 1986; Beysens et al., 1990; Wilkinson et al., 1998; Barmatz et al., 2007). This study shows that near the critical point, the correlation time of fluctuations also separates in multiple domains, similar to the spatial separation of fluctuations.

Once we fit the correlation time of fluctuations (Figure 4B) with Eq. 4, one gets the effective diffusion coefficient D shown in Figure 5A. The critical wavenumber q _c is the inflection point of the correlation time (see the points B1, B2, and B3 in Figure 4B), and they are shown in Figure 5B. Based on the examples shown in Figure 4B, we expect that the correlation time slopes for the original images and their IMFs are almost identical within the experimental variability. Since the slopes of the correlation times at large wavenumber are virtually identical, we anticipate from Figure 4B that the effective diffusion coefficient is the same across the original images and all IMFs. Indeed, as we see from Figure 5A, the effective diffusion coefficients computer for each batch of N _avg = 32 images are very similar within experimental errors. We also notice from Figure 5A example and Table 1 that the effective diffusion coefficient variance increases with the IMF’s order.

The continuous black line that fits the effective diffusion coefficients for the original image in Figure 5A is given by equation D = (3.66 ± 0.06) 10^–8 + (−5.9 ± 0.6) 10^–11 × t (in cm s⁻¹), where t is the temporal variable shown in Figure 5. The above linear fit’s statistics show a χ ² = 5.47 × 10^–18 and an adjusted coefficient of determination R ² = 0.588. Given the significant data variance of the average diffusion coefficient D (see Table 1), even for the original image, the slight negative slope of −1.18×10^–10 cm² s⁻² might not be statistically significant. To determine that the effective diffusion coefficient decreases over time, as the fitting line for the original images in Figure 5A suggests, future studies need to include more experimental data.

We also extracted the value of the critical wavenumber q _c , i.e., the wavenumber for which the correlation time reached the peak in Figure 4B and plotted it in Figure 5B. As expected, the critical wavenumber decreases with the IMF’s order. The horizontal lines in Figure 5B show the average for the corresponding IMFs.

To conclude, using the BEMD decomposition over a narrow wavenumber range, we obtained structure factors with a single power-law exponent, e.g., IMF1 exponent is −2.1 ± 0.6, IMF2 exponent is −4.3 ± 0.3, and IMF3 exponent is −5.6 ± 1.1. We previously reported that depending on the wavenumber range, one gets different power-law exponents for the structure factor in the original image. Using BEMD we can trace those exponents to different IMFs.

To find the values of the correlation length wavenumber q _corr = 1/ξ, we fitted each structure factor with Eq. 3. Although the correlation length wavenumber q _corr seems to change over time slightly (data not shown), we estimated its average value and standard deviation in Table 1. Using q _corr = (4,088 ± 653) cm⁻¹ from the original image, we found that the corresponding correlation length ξ = 1/q _corr ≈ (2.4 ± 0.2) µm. This means q _c × ξ ≤ 1, which shows that we are in the hydrodynamic regime (Onuki, 2002), which justifies the use of the Kawasaki (Kawasaki, 1970) formula shown in Eq. 4 to determine the effective diffusion coefficient.

The experimentally estimated correlation length of ξ = 1/q _corr ≈ (2.4 ± 0.2) µm allowed us to determine the corresponding temperature distance from T _c by using the scaling law for correlation length (Domb et al., 2001): ξ = ξ + T − T c T c − ν , (5) where ξ ₊ = 1.892×10^–10 m, T _c = 318.733 K and the universal exponent is ν = 0.6304 (Moldover et al., 1979; Lecoutre et al., 2009). Our computations show that this 0.2 mK temperature quench stepped through T _c , which seems to be right in the middle of the quench with the first image in the 1922 series 73–121 µK above T _c .

If we assume that the same scaling law given by Eq. 5 applies to the IMFs, we can estimate their “temperatures.” The range of temperatures is as follows: IMF1 823–943 µK above T _c , IMF2 186–217 µK above T _c , and for IMF3 64–75 µK above T _c (see Table 1. Because the BEMD breaks the original image into different wavenumber ranges, the “temperature” of the first IMF will always be larger than the second, and so on.