# Exploring Nanofibrous Networks with X-ray Photon Correlation Spectroscopy

Tomas Rosén<EMAIL_ADDRESS>[  HongRui He  [  Ruifu Wang  [  Korneliya
Gordeyeva  [  Ahmad Reza Motezakker  [  Andrei Fluerasu  [  L. Daniel
Söderberg  [  Benjamin S. Hsiao  [

###### Abstract

Nanofibrous networks are the foundation and natural building strategy for all
life forms on our planet. Apart from providing structural integrity to cells
and tissues, they also provide a porous scaffold allowing transport of
substances, where the resulting properties rely on the nanoscale network
structure. Recently, there has been a great deal of interest in extracting and
reassembling biobased nanofibers to create sustainable, advanced materials
with applications ranging from high-performance textiles to artificial
tissues. However, achieving structural control of the extracted nanofibers is
challenging as it is strongly dependent on the extraction methods and source
materials. Furthermore, the small nanofiber cross-sections and fast Brownian
dynamics make them notoriously difficult to characterize in dispersions. In
this work, we study the diffusive motion of spherical gold nanoparticles in
semi-dilute networks of cellulose nanofibers (CNFs) using X-ray Photon
Correlation Spectroscopy (XPCS). We find that the motion becomes increasingly
_subdiffusive_ with higher CNF concentration, where the dynamics can be
decomposed into several _superdiffusive_ relaxation modes in reciprocal space.
Using simulations of confined Brownian dynamics in combination with simulated
XPCS-experiments, we observe that the dynamic modes can be connected to pore
sizes and inter-pore transport properties in the network. The demonstrated
analytical strategy by combining experiments using tracer particles with a
_digital twin_ may be the key to understand nanoscale properties of
nanofibrous networks.

###### keywords:

Nanofibrous networks $|$ Brownian motion $|$ Anomalous diffusion $|$ Cellulose
nanofibers $|$ X-ray photon correlation spectroscopy $|$ Digital twins

FPT] Department of Fibre and Polymer Technology, Royal Institute of
Technology, 100 44 Stockholm, Sweden [WWSC] Wallenberg Wood Science Center,
Royal Institute of Technology, 100 44 Stockholm, Sweden [SBU] Department of
Chemistry, Stony Brook University, Stony Brook, NY 11794-3400, USA SBU]
Department of Chemistry, Stony Brook University, Stony Brook, NY 11794-3400,
USA SBU] Department of Chemistry, Stony Brook University, Stony Brook, NY
11794-3400, USA FPT] Department of Fibre and Polymer Technology, Royal
Institute of Technology, 100 44 Stockholm, Sweden [WWSC] Wallenberg Wood
Science Center, Royal Institute of Technology, 100 44 Stockholm, Sweden FPT]
Department of Mechanics, KTH Royal Institute of Technology, SE-100 44
Stockholm, Sweden [WWSC] Wallenberg Wood Science Center, Royal Institute of
Technology, 100 44 Stockholm, Sweden BNL] National Synchrotron Light Source
II, Brookhaven National Lab, Upton, NY 11793, USA FPT] Department of Fibre and
Polymer Technology, Royal Institute of Technology, 100 44 Stockholm, Sweden
[WWSC] Wallenberg Wood Science Center, Royal Institute of Technology, 100 44
Stockholm, Sweden SBU] Department of Chemistry, Stony Brook University, Stony
Brook, NY 11794-3400, USA

##  1 Introduction

Networks of semi-flexible nanofibers make up many key architectural elements
in nature, such as actin filaments and microtubules in the cytoskeleton of
eukaryote cells, collagen nanofibrils in connective tissues and bones, chitin
nanofibrils in exoskeletons of arthropods, fibroin nanofibrils in silk from
spiders or silkworms, and cellulose nanofibers (CNFs) in plant cell walls 1,
2, 3, 4. The primary function of the nanofibrous networks is to provide
structural integrity and resist deformation from stretching, compressing,
bending, and shearing forces. Additionally, these networks provide a natural
porous scaffold for transport of essential substances and other nanoscale
objects. Through intricate and complex hierarchical assemblies that have
evolved over millions of years, the nanoscale properties of these fibrous
network can be transferred to macroscale, thus enhancing the chances of
survival of the species they belong to 5.

Nanofibrous networks also play an important role in another modern scenario.
By extracting CNFs from the plant cell wall and reassembling them in a
controlled manner, there have been numerous examples of advanced biobased
nanomaterials made directly from biomass. These materials possess excellent
mechanical properties and functionalizability 6, 7. Some examples include
strong filaments for new textile materials 8, 9, 10, 11, membranes for water
purification 12, 13, aerogels for therapeutic delivery 14, or artificial
tissues and wound dressings in the biomedical field 15, 16. By letting the
CNFs provide a structural scaffold, other functional nanoparticles can be
introduced to make new hybrid materials that are conductive, magnetic, fire
retardant, or provide structural coloring and biosensing abilities 17, 18, 19,
20, 21, 22, 23. The hydrophilic CNFs can further be dispersed in organic
solvents, making them compatible as reinforcement agents in hydrophobic
polymer materials 24, 25.

Despite the wide range of applications for CNF materials, the nanoscale three-
dimensional structures of the CNF networks, giving rise to the material
properties, are notoriously difficult to characterize 8 and therefore to
control. The main techniques for characterizing CNFs include imaging
techniques such as atomic force microscopy (AFM) or transmission electron
microscopy (TEM), which can provide accurate descriptions of the individual
CNF morphology, but fail to capture dynamics and structures in dispersions.
Small-angle X-ray scattering (SAXS) can provide both statistical information
about cross-section distributions 26 and segmental aggregation in dispersions
27, 28. Additionally, rheological studies can address the mechanical response
from the network subject to external deformation 29. Still, there are no
characterization techniques that have been able to provide a deeper insight
into the nanoscale structures of the dynamic CNF networks, and especially the
nanoscale transport properties within the network.

By _in situ_ measurements of the diffusion of spherical gold nanoparticles
(GNPs) using X-ray photon correlation spectroscopy (XPCS), we have probed the
characteristics and nanoparticle transport properties of dilute nanofibrous
(CNF) networks in liquid. In the first step of the developed methodology, we
characterized the morphology of CNFs using AFM and SAXS to find average
nanofiber lengths, cross-sections, degree of fibrillation and segmental
aggregation. Thereafter, we characterized the diffusive motion and dimensions
of GNPs in pure liquid and at different CNF concentrations using XPCS and
SAXS. The data was analyzed to provide the mean diffusive time scales, showing
that with increasing CNF concentration and GNP size, the dynamics become: (1)
slower, (2) more non-uniform, and (3) increasingly subdiffusive.

The experimental data was parametrized as dynamic modes for the case of the
largest GNP, which were fitted to simulated results from a digital twin of the
experiments. This made it possible to extract the physical properties, such as
pore sizes and pore connectivity. Apart from providing new insights into the
CNF system used in this study, the presented methodology can be readily used
to quantify general nanoparticle transport properties of nanofibrous networks
in liquids for various applications, _e.g._ materials science or biomedicine.

Figure 1: Illustration and characterization of the samples used in this study;
two different dispersions of CNFs are used with high oxidation (HO) and low
oxidation (LO); (A) AFM images of typical CNFs (scale bar 1 $\mu$m); (B)
length distribution of CNFs obtained from AFM images; (C) Lorentz-corrected
SAXS data $Iq^{2}$ for CNFs at 0.3 wt%, where solid curves are fits of the
cross-section form factor to determine the degree of fibrillation 26 and the
dashed line shows the lower $q$-limit for this fit; the inset figure shows the
structure factor $S(q)$ of the two dispersions at low $q$; (D) three-
dimensional renderings of the two CNF networks; the higher degree of
fibrillation and less local entanglements (segmental aggregates) likely result
in smaller pores in the HO network (typical pore size illustrated with purple
semi-transparent sphere); (E) illustration of the XPCS experiment; an X-ray
beam is focused on the sample containing a mixture of CNFs and gold
nanoparticles (GNPs), where the scattering intensity $I(q,t)$ is recorded on a
detector and will be dominated by GNPs; (F) representative correlation curves
of $c(\tau)=g_{1}^{2}(\tau)$ at various $q$ for a system with 200 nm GNPs in a
HO CNF network at $0.24$ wt%;  the solid curves show the result from the
CONTIN-algorithm.

##  2 Results

###  2.1 Morphological characterization of CNF samples

Dispersions of CNFs are extracted from paper pulp typically through mechanical
high-pressure homogenization. 6. Prior to homogenization, the material is
usually oxidized (_e.g._ through carboxymethylation or TEMPO-mediated
oxidation) to introduce charged groups, thus lowering the energy requirement
for the CNF extraction due to electrostatic repulsion and allowing the CNFs to
be stably dispersed in water. The individual morphologies of CNFs are strongly
dependent on extraction methods, and specifically to the degree of oxidation
applied to the material prior to homogenization 26, 29. In this study, two
aqueous TEMPO-oxidized CNF dispersions were prepared from wood pulp with two
degrees of oxidation, which introduce different densities of charged COO-
groups on the CNF surfaces and the two dispersions are referred to as HO (high
oxidation) and LO (low oxidation), respectively. Through a recently
demonstrated solvent exchange process 25, the CNFs are stably dispersed in
propylene glycol (PG). Not only is this material interesting in terms of
making novel biocomposites, but it has previously been shown that the solvent
exchange does not affect the CNF morphology, but rather mainly causing a
reduction of Brownian motion owing to the high viscosity of PG, which will be
beneficial in order to capture the dynamics with XPCS.

Representative AFM images of the HO and LO CNFs are shown in Fig. 1A. Length
distributions are determined in Fig. 1B, where it is found as expected that
the higher oxidation treatment gives rise to shorter CNFs on average
($\bar{L}=381$ nm for HO CNFs and $\bar{L}=497$ nm for LO CNFs).
Characterization of the dispersions at 0.3 wt% using SAXS (Fig. 1C) and
fitting the Lorentz-corrected intensity curves $Iq^{2}$ with an assumption of
quantized polydispersity 26 revealed that the HO CNFs also had a higher degree
of fibrillation (71.2 % compared to 52.6 %). The inset figure in Fig. 1C shows
the structure factor $S(q)$ in the two dispersions, where a value of $S(q)=1$
would indicate no higher order structures, which is typical in dilute samples.
The exponential shape of the LO CNF indicates strong segmental aggregation,
_i.e._ locations where fibrils are locally entangled and thus in close
proximity 27, 28. The structure factor is fitted to determine the average
diameter of the segmental aggregates 28 (solid curve in inset figure), which
is found to be approximately 42 nm. The structure factor of the HO CNF instead
shows a peak at $q\approx 0.01$ Å-1 indicating that the stronger electrostatic
repulsion causes nearby CNF segments to be separated at an average spacing of
$d=2\pi/q\approx 60$ nm, similar to values found in dispersions of cellulose
nanocrystals 30.

The lower degree of fibrillation, longer fibrils and segmental aggregation
suggest that the mass distribution inside the LO dispersion is more
heterogeneous than the HO CNF, and that consequently pore sizes in the network
will be larger for a given concentration (schematically illustrated in Fig.
1D). Interestingly, despite these differences between the two dispersions, the
rheology stays practically the same, which is shown in Fig. S1 in the
supplementary information (SI). In this work, three different concentrations
of CNFs are studied, which correspond to a transition from a dilute to a semi-
dilute regime: 0.08, 0.16, and 0.24 wt%. More details of these initial
characterizations of the CNF dispersions are provided in the Materials and
Methods section.

The structure of the CNF network is probed indirectly through the motion of
GNPs that are added to the dispersion. The GNPs themselves are sterically
stabilized with a hydrophilic polymer coating to keep them stable in
electrolytes regardless of pH (more details under Materials and Methods). This
also ensures negligible electrostatic interactions with the CNFs, and that
they can be passively transported through the network.  The volume fraction of
the GNPs is $\Phi\approx 0.0001$ and no aggregation is observed through the
structure factor of the system, thus neglecting any influence of interactions
between GNPs. Three different sizes of GNPs are used here: 50, 100, and 200
nm. Using the theory by Ogston 31, 32 for pore sizes in networks of randomly
oriented rods (see Fig. S2 in SI), we find that the expected pore sizes at the
present CNF concentrations are of the same order as the sizes of the GNPs.
Thus, we expect to be able to study the effect on the GNP dynamics from an
unconfined to a confined surrounding. Three-dimensional renderings of the
approximate experimental conditions in the CNF/GNP dispersions are provided in
a supplementary movie.

###  2.2 XPCS: experiments, analysis and validation

The CNF/GNP samples are mixed and injected in quartz capillaries, which are
studied with XPCS at the CHX beamline (11-ID) at the National Synchrotron
Light Source II (NSLS-II), Brookhaven National Laboratory, USA. The principle
of XPCS is shown in Fig. 1E and described in detail in the Materials and
Methods section. In brief, a coherent X-ray beam (with wavelength $\lambda$)
is focused on the sample and the scattering intensity $I(q)$ is recorded on a
detector at various $q=(4\pi/\lambda)\sin\theta$ (with scattering angle
$2\theta$). Given the much higher electron density of gold, the scattering
intensity is completely dominated by the scattering from the GNPs. The
interference maxima on the detector, the so called _speckles_ , contain
information about spatial correlations between particles in the sample at a
length scale $2\pi/q$. By recording a sequence of speckle patterns with delay
time $\tau$, the average second-order autocorrelation of the speckle
intensities over time $g_{2}(\tau)$ at a certain $q$ can be extracted. Through
the Siegert relation, the square first-order autocorrelation is obtained by
$c(\tau)=g_{1}^{2}(\tau)=(g_{2}-g_{\infty})/\beta$, where $g_{\infty}$ is the
baseline and $\beta$ is the speckle contrast 33. As thermal fluctuations of
particle positions give rise to decaying autocorrelations of speckle
intensities, the decay of $c(\tau)$ thus captures the Brownian motion in the
sample.

The square first-order autocorrelation $c(\tau)$ is assumed to consist of
various exponential dynamic modes
$c(\tau)=\sum_{i}G^{(d)}_{i}\exp(-2\Gamma^{(d)}_{i}\tau)$, where the discrete
distribution $G^{(d)}_{i}$ of log-spaced relaxation rates $\Gamma^{(d)}_{i}$
is found through a regularized inverse Laplace transform using the CONTIN-
algorithm 34, 35 (see Fig. 1F). The mean decay rate $\bar{\Gamma}$ is found
through $\ln\bar{\Gamma}=\sum_{i}G^{(d)}_{i}\ln\Gamma^{(d)}_{i}$. Furthermore,
a polydispersity index ($PDI$) is extracted from the relaxation rate
distribution to describe the uniformity of dynamics using the method of
cumulants 36 (see details in the Materials and Methods section). It was also
noted that the data could be well fitted with a stretched exponential
function, but there are no features in the data to suggest that the decay is
truly non-exponential. More discussion and details of the various fitting
techniques is found in Figs. S3 and S4 in SI.

The potential for radiation damage was carefully assessed and in the present
systems it was found that radiation was starting to influence the dynamics at
a dose equivalent to 5 s of exposure from the non-attenuated beam (see Fig. S5
in SI for details). In all the experiments in this work, we study 2 different
attenuation levels for each sample: 0 % (full beam) and 81 % beam attenuation.
As expected, no significant attenuation-dependent trends were found in the
data since the total full-beam exposure never exceeded 1.4 s. Every
measurement is repeated twice at a different spatial location in the sample,
and we ensure that the sample is in thermal equilibrium with spatially
homogeneous dynamics (see Materials and Methods section for details).
Furthermore, any direction-dependent dynamics arising from potential gravity-
induced sedimentation were examined and determined to be negligible (see Fig
S6 in SI).

Figure 2:  Mean results from the XPCS analysis for the three sizes of GNPs;
(A) and (B) show the Brownian diffusion coefficient $D_{B}$ of GNPs in HO and
LO CNF, respectively, at various nanofiber concentrations compared to free
diffusion in the solvent (propylene glycol) $D_{0}$ assuming
$\Gamma=D_{B}q^{2}$; the solid black curve shows the zero-shear viscosity of
the CNF suspension $\eta$ compared to the viscosity $\eta_{0}$ of the solvent;
(C) and (D) shows the polydispersity index $PDI$ as function of concentration
compared to $PDI_{0}$ of free diffusion; (E) and (F) shows the dynamic scaling
exponent $n$ assuming $\Gamma=(Dq^{2})^{1/n}$, where a value of $n<1$
indicates subdiffusive dynamics. The symbol sizes correspond to the two
attenuation levels used in the study.

When GNPs of uniform sizes are moving freely in the solvent without CNFs,
$c(\tau)$ is expected to decay exponentially (low $PDI$) with
$\bar{\Gamma}=D_{0}q^{2}$ ($D_{0}$ is the Brownian diffusion constant in free
solvent), which was analyzed in initial experiments for the three GNP sizes.
The results in Figs. S7 and S8 in SI illustrate that there is some inherent
non-uniformity of dynamics in the GNP systems owing to a spread of actual
geometric sizes (found with SAXS). Also interesting to note is that the
polymer coating (with thickness $\lesssim 5$ nm) gives rise to discrepancies
between geometric and dynamic radii. This effect is naturally larger for the
smallest (50 nm) GNPs, where the mean geometric diameter is $39.9$ nm and mean
hydrodynamic diameter is $57.5$ nm. This discrepancy does not play any major
role to the study, where we investigate the _relative_ effect of the dynamics
compared to free diffusion in the solvent,  _i.e._ the geometric diameter will
stay the same but the hydrodynamic diameter will naturally change.

###  2.3 Mean dynamics of GNPs in CNF dispersions

Adding CNFs to the system and assuming a dispersion of straight rigid
nanorods, previous theoretical works 37, 38, 39 suggest that the GNP motion is
not hindered if the characteristic pore size (_mesh size_) of the network is
much larger than the GNPs, _i.e._ they will have mean diffusion constants
being the same as in free diffusion $D_{B}=D_{0}$. At higher concentrations,
when the mesh size is much smaller than the particle, the motion is determined
by the momentum transfer within the network, which on the macroscopic scale is
interpreted as an effective zero-shear viscosity $\eta$ higher than the
solvent viscosity $\eta_{0}$. Consequently, the scaling of $D_{B}/D_{0}$
should approach the scaling of $\eta_{0}/\eta$. In Fig. 2A, we find that the
expected behavior seems to be captured fairly well in the HO CNF dispersion.
Most striking is the diffusion of 50 nm GNPs with dynamics similar to free
diffusion ($D_{B}=D_{0}$) below a concentration of 0.2 wt% and then possibly
approaching the scaling of $\eta$. Similar trends are seen also for the 100
and 200 nm GNPs, where of course the threshold concentration is lower with
larger GNPs. Although the 50 nm GNPs show similar behavior in the LO CNF (see
Fig. 2B), the larger GNPs are not as affected by CNF concentration between
0.16 and 0.24 wt%.

The mean diffusion constant however does not fully reflect the dynamic
behavior in the systems. The increasing CNF concentration also results in less
uniform GNP dynamics, reflected by higher values of $PDI$ compared to the
$PDI_{0}$ in free diffusion (see Fig. 2C-D). The semi-dilute CNF network thus
seems to slow down the GNP motion unevenly, with both faster and slower moving
particles in the system.

An even stronger indication that the dynamics is more complicated is found
when studying the $q$-dependence of the relaxation rate $\bar{\Gamma}$ in Fig.
2E-F. In a system with normal Brownian motion, the mean square displacement
($MSD$) of particles is known to be linear in time, with the diffusion
constant being the proportionality constant ($MSD\varpropto D_{B}\tau$) with
unit m2/s. However, in a confined system, the dynamics can become
_subdiffusive_ with $MSD\varpropto D\tau^{n}$, where $n<1$ and the
proportionality constant gets the unit m2/sn 40. To ensure this dimension of
$D$, the relaxation rate is assumed the form $\Gamma=(Dq^{2})^{1/n}$ and it is
found in Fig. 2E-F that the GNP dynamics are indeed subdiffusive in both CNF
networks. Consequently, it can be argued that the comparison of $D_{B}$ in
Fig. 2A-B is not suitable, as these comparisons are based in an assumption of
$n=1$. Furthermore, the values of $D$ are not comparable as their unit depends
on the specific value of $n$.  For this reason, we will now investigate the
origin of this apparent subdiffusivity by studying more detailed features of
the dynamic spectra and compare with a numerical simulation of the system.

Figure 3: Illustration of the extraction of dynamic modes using a lognormal
deconvolution; (A) the resulting discrete distribution $G_{i}^{(d)}$ of
relaxation rates $\Gamma_{i}^{(d)}$ at $q=0.0025$ Å-1 as circles and the solid
black line shows the continuous distribution $G(\Gamma)$ obtained by fitting
three lognormal curves with mean values $\Gamma_{1}$, $\Gamma_{2}$ and
$\Gamma_{3}$ and magnitudes $G_{1},G_{2}$ and $G_{3}$; (B) and (C) shows the
evolution of the three dynamic modes for 200 nm GNPs with increasing $q$
(illustrated with arrows) for the HO and LO CNF networks, respectively, at
0.24 wt%.

###  2.4 Dynamic modes in the subdiffusive systems

Although various studies have been presented recently to analyze anomalous
diffusion with XPCS 41, 42, 43, a different approach is employed here, similar
to the method by Andrews et al. 35.  Since Fig. 2 displays clear trends of
slower, more polydipserse and subdiffusive dynamics with increasing CNF
concentration and GNP size, we focus our attention to the system of large (200
nm) GNPs at high CNF concentration.

Fig. 3A shows the resulting distribution $G^{(d)}_{i}$ of relaxation rates
$\Gamma^{(d)}_{i}$ from the dynamics of 200 nm GNPs in HO CNF at 0.24 wt%,
where we can distinguish three individual peaks. Through a lognormal
deconvolution, the discrete spectra is converted to a continuous function
using three modes with mean relaxation rates $\Gamma_{1-3}$ and magnitudes
$G_{1-3}$ as seen in Fig. 3B. The evolution of these dynamic modes can then be
studied individually with respect to $q$ as done in Fig. 3B-C for HO and LO
CNF at 0.24 wt%, respectively. Interestingly, the individual relaxation rates
$\Gamma_{1-3}$ do not monotonically increase with $q$. Instead, for both CNF
samples, there is a clear switching of which modes are relevant at a given
$q$-range, where the slower mode ($\Gamma_{3}$) contributes more to dynamics
at low $q$, while the faster mode ($\Gamma_{1}$) contributes more at high $q$.
The main difference between the two samples is that the contributions of the
slow and intermediate modes ($\Gamma_{3}$ and $\Gamma_{2}$) are much lower at
higher $q$.

Figure 4: Evolution of the relaxation rates of the dynamic modes $\Gamma_{1}$,
$\Gamma_{2}$ and $\Gamma_{3}$ as function of $q$ for the 200 nm GNPs; the line
width is scaled with the magnitude of the mode $G_{1},G_{2}$ or $G_{3}$, with
linear interpolation between the $q$-values (note that $G_{1}+G_{2}+G_{3}=1$
at any given $q$); the thin dashed and dotted lines indicate expected scaling
behavior for Brownian ($\Gamma\varpropto q^{2}$) and ballistic
($\Gamma\varpropto q$) dynamics, respectively; the thick black solid and
dashed lines show the free diffusion of GNPs ($\Gamma_{\text{ref}}$) and mean
relaxation rate ($\bar{\Gamma}$), respectively; (A)-(C) show the dynamic
behavior in HO CNF with increasing concentration (0.08, 0.16 and 0.24 wt%),
while (D)-(E) show the behavior in LO CNF at the same concentrations.

In Fig. 4, the switching of dynamic modes is illustrated in more detail. The
figure shows the relaxation rates as function of $q$, where the line width of
the three modes $\Gamma_{1-3}$ is scaled with their individual magnitude
$G_{1-3}$. The mean relaxation rate $\bar{\Gamma}(q)$ and the reference of
freely diffusing GNPs $\Gamma_{\text{ref}}(q)$ are plotted with thick black
dashed and solid lines, respectively. The thin dashed and dotted lines
correspond to the expected scaling for Brownian dynamics ($n=1$ and
$\Gamma\varpropto q^{2}$) and ballistic dynamics ($n=2$ and $\Gamma\varpropto
q$), respectively. There are some interesting similarities and differences to
highlight with respect to CNF concentration in the two systems. At some low
threshold concentration, there is a transition from single-mode dynamics
similar to the free diffusion ($\bar{\Gamma}\approx\Gamma_{\text{ref}}$) to
dynamics governed by two modes: (1) one faster mode $\Gamma_{1}$ that is
basically constant with $q$ but only dominating at high $q$ and (2) one slower
mode $\Gamma_{2}$ that still remains Brownian ($n=1$), but with slower rates
than the reference and dominating at lower $q$. The main difference between
the two CNFs is that the splitting of modes occurs at higher concentration in
the LO CNF within the given $q$-range. In Fig. 4D, for LO CNF at a
concentration of 0.08 wt%, the splitting of modes is clearly visible at around
$q=0.0045$ Å-1, which corresponds to length scales of around 140 nm. In Fig.
4A on the other hand, for HO CNF at the same concentration, this mode-
splitting seems to occur at higher $q$ values, but it is likely that there
exist a lower concentration where the mode-splitting is at a similar $q$-value
as for the LO CNF.

At higher concentrations, the faster constant mode $\Gamma_{1}$ remains at the
same level as for the lower concentrations, while the slower Brownian mode
$\Gamma_{2}$ decreases. However, at some concentration between 0.16 and 0.24
wt%, it stagnates and turns towards a more ballistic behavior with $n\approx
2$ ($\Gamma\varpropto q$). At this stage, a new even slower Brownian mode
appears at low $q$ ($\Gamma_{3}$). Interestingly, even though none of the
contributing modes is associated with subdiffusive dynamics ($n<1$), the mean
rate $\bar{\Gamma}$ is slightly subdiffusive, which is also the reason for the
results presented earlier in Fig. 2. _The overall subdiffusive behavior is
thus the result of superdiffusive modes contributing differently at different
$q$._

Additionally, since the trends with increasing concentration are similar for
the two CNF networks, it is clear that the LO dispersion is going through the
same stages as the HO dispersion, but generally at higher concentrations. This
is of course in line with our assumption that the mass distribution of CNFs in
the LO dispersion is more heterogenous and thus has a more open pore space.

Although the methodology provides detailed trends of the dynamic modes in
reciprocal space, it is more complicated to relate these trends to actual
physical properties of the CNF network. Instead of relying on theoretical
models based on idealized systems to describe the trends, we use a concept of
a _digital twin_ ,  _i.e._ a numerical simulation of the system, where the
same observable quantities are extracted and matched to the experiment. The
same approach was recently used to study dynamics of flowing CNFs _in situ_
using SAXS 44 and to study cross-sections of CNFs 26 with SAXS/WAXS.

Figure 5: Results from the digital twin of confined Brownian motion of 200 nm
spherical nanoparticles (NP) in propylene glycol (at room temperature); (A)
the simulation domain is divided into cells of size 50 nm, within which the
midpoint of the NP can move freely (effective pore size is thus 250 nm); the
probability for an NP to cross to another cell is determined by the
permeability $P$, where $P\rightarrow\infty$ leads to free diffusion while
$P=0$ leads to the NP being trapped inside the cell; (B) mean square
displacement ($MSD$) as function of delay time $\tau$ and permeability $P$,
where Brownian dynamics would lead to $MSD\varpropto\tau$; the shaded region
shows the approximate length scales of the XPCS-experiment, _i.e._ between
$(2\pi/q_{\text{max}})^{2}$ and $(2\pi/q_{\text{min}})^{2}$; (C) shows the
square first order autocorrelation curves $c(\tau)$ at $q=0.0045$ Å-1 at three
levels of $P$ with solid curves being the fitting result; the inset figure
shows the corresponding dynamic modes after the lognormal deconvolution;
(D)-(I) show the evolution of the relaxation rates $\Gamma_{1}(q)$ and
$\Gamma_{2}(q)$ with increasing permeability $P$ illustrated in the same way
as in Fig. 4; the black line in (D)-(F) shows the $q$-value used for the data
in (C).

###  2.5 Development of a digital twin

Subdiffusive motion of NPs in various semi-flexible networks has been
demonstrated previously through numerical simulations, indicating that it is a
natural consequence of the confined motion inside the network 40, 45, 46.
However, it is not clear from these models what such behavior would look like
in an XPCS experiment, where the dynamics are analyzed in reciprocal space
rather than real space. Therefore, in order to investigate our experimental
system, we set up numerical simulations of confined Brownian spherical
nanoparticles (NPs) of 200 nm, where their individual positions are used to
simulate speckle patterns through a Fourier transform, which in turn can be
analyzed in the same way as in the experiment. The model is practically a 2D
version of the simulation method presented by Ernst et al. 45, where NPs are
moving within a square grid providing equally sized cells, which are here
chosen to have side 50 nm (see Fig. 5A). Note that the cell size in the model
is significantly smaller than the NP size, but since it is the midpoint that
is restricted by the cell boundaries, the true pore size is actually 250 nm
(the 200 nm NP can move undisturbed $\pm 25$ nm from the middle of the cell).
When a Brownian step causes the midpoint of the NP to cross the cell boundary,
there is a probability of passing determined by a permeability $P$ (unit
s-1/2, see Materials and Methods section for details). If the NP is not
allowed to pass, it is simply bouncing off the boundary and remaining inside
the cell. Thus, if $P\rightarrow\infty$, the NP is freely diffusing without
any effect of the cell boundaries, but if $P\rightarrow 0$, the NP is
contained within the cell where it started. Fig. 5B shows the $MSD$ versus
delay time $\tau$ of this system, where $P\geq 5$ demonstrates Brownian
dynamics with $MSD\varpropto\tau$ ($n=1$). As the NP motion gets more and more
constrained by the cell boundaries around $P=1$, the system still remains
almost unaffected at small $\tau$, but is clearly subdiffusive at intermediate
$\tau$ when the particle has a ”jumping” motion from one cell to another 45,
46. When $P=0$, the value of $MSD$ is naturally limited by the size of the
cell and remains constant for all larger $\tau$.

In reciprocal space, at $q=0.0045$ Å-1, after performing the same dynamic mode
analysis, we find very similar behavior as in the experiment (see Fig. 5C).
The dynamic behavior at $P=20$ clearly is described by one single dynamic
mode, with a faster mode appearing as the permeability decreases. The entire
evolution of the splitting of modes with decreasing permeability is
illustrated in Figs. 5D-I, which looks almost identical to the trends found in
the experimental data. From being almost indistinguishable from free diffusion
at $P=20$, the dynamic behavior is described with two separate modes at lower
$P$, where the faster mode remains stationary both with $q$ and with
decreasing $P$. The slower mode quickly slows down further with lower
permeability, and with lower and lower contribution at higher $q$. Just like
in the experiments, although both modes are superdiffusive with $n>1$, the
mean relaxation rate still shows a subdiffusive behavior (thick dashed line).
Eventually, as the contribution of the slower mode diminishes, the dynamic
behavior at $P=0$, where NPs are completely contained inside the cell, is
dominated by the fast stationary mode. The two modes can thus be explained by
the two types of NP dynamics at different length scales using this simplified
model: (fast) inter-cell dynamics and (slow) macroscopic ”jumping” dynamics
between cells.

Figure 6:  Schematic illustration of the dynamics of 200 nm GNPs in the two
CNF dispersions during network formation with increasing concentration,
quantified by permeability $P$ by matching dynamic features to the digital
twin; at dilute conditions, the GNPs are moving freely within the solvent and
are unaffected by the fast Brownian CNFs; as concentration reaches a level
where the average pore sizes are smaller than the GNPs, a local cage structure
of nanofibers will restrict the GNP motion, but it has some probability of
moving due to structural changes of the dynamic network; the local structure
and dynamics would be the source of the apparent permeability $P$ and
unrestricted wiggle-room of $\pm 25$ nm; the combination of less dynamic
reconfiguration, and smaller pore sizes at higher concentrations gives rise to
a plateau of the permeability.

###  2.6 Discussion

How can the above observation be connected to the network properties of the
two CNF dispersions in the experiment and what do the simulated cell size and
permeability correspond to for the real network structure? The GNPs undergo
Brownian motion that is determined mainly by the interaction with solvent
molecules. Individual CNFs are also highly Brownian with much faster dynamics
than the GNPs, which can be interpreted on a nanometer scale as a solvent with
seemingly higher apparent viscosity, thus slowing down the motion of the GNPs.
At a higher CNF concentration, when the average pore becomes smaller than the
GNP size, a loose cage of CNFs is also formed around the individual GNPs. The
cage is naturally slightly larger than the GNP and provides some restriction
to its motion. However, the GNP can escape its local cage as the surrounding
dynamic CNF network will constantly undergo structural changes causing escape
routes for the GNP. This probability of the GNP escaping its local cage could
thus be perceived as the source of the permeability in the model. With an
increasing CNF concentration, the cage will become stronger and GNPs will
remain essentially stationary, which corresponds to the case of zero
permeability.

The main permeability-dependent feature of the dynamic modes in Fig. 5, is the
decrease of the second mode $\Gamma_{2}$ with decreasing $P$, which also seems
to correspond to the effect observed in the experimental system. By matching
the mean rate of $\Gamma_{2}$ between the experiment and the digital twin, we
can obtain quantitative values of $P$ in the experimental system as
illustrated in Fig. 6 (see details in Fig. S10 in SI). At low concentrations
the permeability is higher in the LO CNF compared to HO CNF at a given
concentration, which makes sense owing to the heterogeneous mass distribution
and therefore less direct interactions for the GNPs with the CNFs.
Interestingly, both for the HO CNF and LO CNF, we observe a slight increase of
permeability values in a certain concentration range. We hypothesize that this
is linked to the actual network formation, where the Brownian dynamics of CNFs
slow down rapidly, thus decreasing the CNF-GNP interactions despite the
smaller pore sizes. When the network is formed, the permeability $P$ again
decreases rapidly with concentration. Somewhat counter-intuitive, this plateau
of $P$ occurs at higher concentrations for LO CNF suggesting that the
heterogeneity of the dispersion causes network formation at a higher volume
fraction despite the larger CNF aspect ratios. At high concentrations, the LO
CNF will again have higher permeability than the HO CNF owing to the larger
physical pore sizes in the network.

It is natural to assume that the model cell size (or rather wiggle-room)
should also be linked to CNF concentration. However, doing the same analysis
for both larger (75 nm) and smaller (25 nm) model cell sizes reveals a
different scenario for the splitting of modes (see Fig. S9 in SI for details),
which indicates that both CNFs in our study indeed create a local surrounding
of the 200 nm GNP with an unrestricted wiggle-room of $\approx 50$ nm. We
hypothesize that this apparent wiggle-room might rather be influenced by the
flexibility of individual CNFs. The flexibility might also be influencing the
magnitude of the stationary mode ($\Gamma_{1}$). In a recent similar XPCS
study by Reiser et al. 42, the dynamics of spherical NPs in networks of
wormlike micelles also revealed a stationary mode, which was argued to be
related to the Kuhn length of the individual micelle. This could naturally not
be investigated in this work as the cell boundaries are fixed.

Although revealing much of the relevant dynamics in the experiment, there are
some features that the reduced model cannot capture. For example, the
emergence of a third dynamic mode ($\Gamma_{3}$), which suggests the real
system also includes dynamics on a higher hierarchical level. Obviously, the
”network” in the reduced model is static, while the network of CNFs is
dynamic, and it is therefore intuitive to assume that there should be an
additional dynamic mode connected to the network dynamics, although this
requires further investigation.

##  3 Conclusions

This work addressed the dynamics of sterically stabilized spherical gold
nanoparticles (GNPs) in semi-dilute dispersions of cellulose nanofibers (CNFs)
using X-ray photon correlation spectroscopy (XPCS). CNFs were prepared at two
different degrees of oxidation, which provided differences in the CNF length
distribution as well as the degree of fibrillation, which in turn had a clear
effect on the pore spaces in the network at a given concentration. From the
smallest GNPs (50 nm), it is found that their dynamics are unaffected by CNFs
up until a concentration of $\approx 0.2$ wt% in both systems, suggesting that
the undisturbed pore sizes in the semi-dilute networks were of the similar
order. However, the dynamics of the larger GNPs (200 nm) could be described by
several dynamic modes, which through a comparison with a digital twin were
found to be connected to the GNPs’ ability to escape their local surrounding
of CNFs.  This in turn could be used to quantify the effective permeability of
the GNPs in the dynamic semi-dilute CNF network and to characterize threshold
concentrations for network formation.

The fact that the relatively simple model of confined Brownian motion was able
to capture the main dynamic modes from the demonstrated experiment in this
study opens up vast possibilities for exploring nanofibrous networks by
characterizing the dynamics of tracer NPs using XPCS. The obvious next step
would be to study a wider $q$-range with the digital twin and carry out a
proper parameter study of both pore size and permeability to pinpoint how
these parameters can affect the measurable quantities. This in turn can aid in
designing new improved experiments. As a future outlook, the digital twin can
be improved using coarse-grained molecular dynamics simulations of a realistic
nanofibrous network to directly connect the dynamic modes to CNF properties
and interactions, possibly also capturing features that cannot be measured
experimentally. This could potentially also be used as training data for
machine learning, allowing physical properties to be directly obtained from
the experimental data.

Finally, we would like to emphasize that the methodology is not limited to our
experimental system but can be used more generally to other nanofibrous
networks, _e.g._ from polymeric, ceramic or carbon fiber networks.
Technically, the same methodology could also be applied to systems without a
bounding solvent, _e.g._ transport in membranes or aerogels, although it might
me more challenging/hazardous experimentally to study the unbounded NP
transport through these networks.

##  4 Materials and Methods

### Samples

The spherical gold nanoparticles (GNPs; diameters of 50 nm, 100 nm and 200 nm)
dispersed in propylene glycol were purchased from Nanopartz™ and were
functionalized with a proprietary hydrophilic polymer with a MW 1500 Da
(product name H11-PHIL-PRG). The dispersion was characterized by the vendor
using both transmission electron microscopy (TEM) and dynamic light scattering
(DLS) ensuring the quality and uniformity of the GNPs. Analysis in the present
work using SAXS further confirmed high uniformity of the GNP sizes (see Figs.
S7 in SI for details).

Never dried sulfite softwood pulp (Domsjö dissolving pulp) was provided by
RISE Research Institute of Sweden and was used for preparation of cellulose
nanofibers (CNFs). Sodium hypochlorite ($6$-$12$ wt%, Alfa Aesar), TEMPO
(2,2,6,6-tetramethyl-1-piperidinyloxy free radical, $\geq 98$ %, Alfa Aesar),
sodium hydroxide ($\geq 99.2$ %, VWR Chemicals), sodium bromide (BioUltra,
$\geq 99.5$ %, Sigma Aldrich), concentrated hydrochloric acid (VWR Chemicals,
$\geq 35$ %) were used as received. The TEMPO-mediated oxidation was performed
on never dried sulfite softwood pulp (41.6 wt%) following the protocol
described elsewhere 47, where two degrees of oxidation were achieved by adding
3 and 1.5 mmol of oxidant per gram of pulp for high oxidation (HO) and low
oxidation (LO) CNFs, respectively. To extract CNFs, treated pulp was
fibrillated by a high pressure (1600 bars) Microfluidizer (M-110EH,
Microfluidics) with a 400/200 $\mu$m (1 pass) and a 200/100 $\mu$m (4 passes)
wide chambers connected in series. The obtained 1 wt% CNF gel was further
diluted to 0.3 wt% and homogenized by Ultra-Turrax dispersing tool (IKA,
Sweden) for 10 min at 10000 rpm. Conductometric titration was performed to
estimate the total CNF charge according to protocols described elsewhere 48
and average charge densities of carboxylic (COO-) groups were estimated to be
$448\pm 16$ and $1023\pm 26$ $\mu$mol/g for the LO and HO CNF, respectively.
The solvent exchange process to propylene glycol (PG) followed the same
procedure as described elsewhere 25, and the resulting PG-CNF dispersions were
used to create dilutions of $0.1,0.2$ and 0.3 wt%. The final CNF
concentrations $0.08,0.16$ and 0.24 wt% are reached after mixing with the GNP
dispersions.

The mixing of GNPs into the CNF dispersions resulted in a final GNP
concentration of 0.2 wt%, which is comparable to the weight fraction of CNFs.
However, since the volume of a single GNP is much larger than a CNF, the
number density in the network is still quite low (see examples of 3D rendered
networks in the supplementary movie), and thus not affecting the overall
structure of the network. The samples were mixed and dispersed through vortex-
mixing (VWR Fixed -Speed Mini) and remain static for 24 hours before the
experiment. Before the actual XPCS measurements, the samples were re-dispersed
by vortex-mixing (ORIGINAL Vortex-Genie 2) and then put into an ultrasonic
bath for 5 minutes to prevent precipitation of nanoparticles. The samples were
then injected to thin-walled quartz capillary tubes for the XPCS experiments
with an approximate time of 1-2 h between injection and measurement.

### Atomic Force Microscopy (AFM)

Topography of particles was visualized with the help of an Atomic Force
Microscope (AFM, Multimode 8, Bruker, USA) in Tapping mode in air. The
original dispersion of CNFs (0.3 wt% in DI water or PG) was diluted to 0.005
wt% and homogenized for 10 min at 10000 rpm using Ultra Turrax (IKA T25D,
Sweden). Remaining aggregates were removed by centrifugation for 1 hour at
4000 rpm. Supernatant was collected and mixed on Vortex Genie 2 (Scientific
Industries Inc., USA) for 5 min right prior casting. A 20 $\mu$L drop of
(3-aminopropyl)thriethoxysilane (99%, Sigma Aldrich) was placed and kept on a
freshly cleaved mica for 30 s, whereafter it was blown away intensively by
nitrogen gas. A 20 $\mu$L drop of CNF dispersion (0.005 wt%) was placed
rapidly on functionalized mica and kept for 30 s until being blown away by
nitrogen gas in a similar manner. Obtained substrates were left to dry
overnight. Diameters and lengths of CNFs were estimated from height images
using software Nanoscope Analysis and ImageJ, acquiring above 300 and 500
measurements, respectively. With the help of software Origin 2021 statistical
analysis was performed on collected values and by fitting a Gaussian
distribution function it was estimated that average diameters (dav) for high
and low charge CNF were $2.73\pm 0.63$ and $2.56\pm 0.51$ nm, respectively.
Both CNF types demonstrated a broad distribution of lengths varying from 50 nm
to 1.6 $\mu$m. By fitting a lognormal distribution, the mean lengths $\bar{L}$
were determined.

### Small-angle X-ray scattering (SAXS)

The SAXS experiments were performed at the LIX beamline (16-ID) at the
National Synchrotron Light Source II (NSLS-II), Brookhaven National
Laboratory, USA. The CNF samples (without GNPs) are injected into liquid
sample holders with mica windows, which are scanned at five different
positions with 1 s exposure at each position, where a mean was used for
analysis. The wavelength was $\lambda=0.82$ Å and the sample-detector distance
was 3.7 m. The beam size was approximately 50 $\times$ 50 $\mu$m2. The
scattered X-ray intensity $I(q)$ was recorded on a Pilatus 1 M detector with
pixel size 172 $\times$ 172 $\mu$m2 at a range of $q=0.006-0.2$ Å-1, with
scattering vector magnitude $q=(4\pi/\lambda)\sin\theta$ (the angle between
incident and scattered light is 2$\theta$). The background scattering
intensity was subtracted using scans of DI water. The degree of fibrillation
of the CNFs is determined with the fitting method by Rosén et al. 26 applied
to the Lorentz-corrected SAXS data $Iq^{2}$ between 0.02 and 0.2 Å-1,
resulting in the form factor of the system $P(q)$. The structure factor is
found through $S(q)=I(q)/P(q)$, where the radius of segmental aggregates $R$
can be found by fitting $S(q)=A\exp[-(0.396R)^{2}q^{2}]+1$ between 0.005 and
0.05 Å-1 according to previous work28.

### Rheology

Steady shear viscosity measurements were performed on Rheometer DHR-2 (TA
Instruments, USA) using plate-on-plate geometry (diameter 25 mm) with the gap
set to 1 mm. Prior to the measurements, the CNF dispersions were degassed by
vacuum to remove air bubbles and left to rest at normal conditions for at
least overnight. The experiment consisted of two steps: sample conditioning at
25 ∘C for 5 min and flow sweep test with shear rate varying between 0.1 and
100 s-1. Data was acquired in a logarithmic manner collecting 10 points per
decade with equilibration time set to 5 s and averaging time reaching 30 s.

The approximate zero-shear viscosity $\eta$ of the CNF dispersions was
obtained by fitting the experimental curves. More details about the
rheological characterization is provided in Fig. S1 in SI.

### X-ray photon correlation spectroscopy (XPCS)

The XPCS experiments were performed at the CHX beamline (11-ID) at NSLS-II,
BNL. An X-ray beam with size 40 $\times$ 40 $\mu$m2 and wavelength
$\lambda=1.28$ Å was focused on the sample and the scattered X-rays were
collected on a detector (Dectris 3-Eiger X, pixel size 75 $\times$ 75 $\mu$m2)
at a sample-detector distance of 16 m covering $q=0.002-0.04$ Å-1. Each sample
was scanned with 1000 images at three different exposure times (1.3, 10 and
100 ms; also setting the frame rate) and four attenuation levels (0, 81, 96.4
and 99.3 %). Each scan was repeated twice, yielding a total of 24 measurements
per sample, which all were taken at different positions in the capillary. It
was however found that only the data with the lowest exposure time of 1.3 ms
(frame rate 750 Hz) could be used for analysis and the scattering intensity
using the highest attenuation levels (96.4 and 99.3 %) was too low to be used
for standard XPCS analysis. An assessment of the influence of the beam on the
sample dynamics revealed that the samples were affected with a total dose
equivalent to 5 s of exposure to the non-attenuated beam (see Fig. S5 in SI
for details), which is approximately four times higher than the maximum dose
in any experiment used for analysis. The auto-correlation $g_{2}(q,\tau)$ at a
certain delay time $\tau$ of the intensity at a certain pixel $I(q,t)$ is
determined both through the one-time correlation function (1T) obtained
through a multi-tau algorithm 49:

$g_{2}^{(1T)}(q,\tau)=\frac{\langle I(q,t)I(q,t+\tau)\rangle}{\langle
I(q,t)\rangle^{2}}$ (1)

or through the mean over all $\tau=|t_{2}-t_{1}|$ of the two-time correlation
function (2T) 50:

$g_{2}^{(2T)}(q,\tau)=\left\langle\frac{\langle
I(q,t_{1})I(q,t_{2})\rangle}{\langle I(q,t_{1})\rangle\langle
I(q,t_{2})\rangle}\right\rangle_{\tau=|t_{2}-t_{1}|}.$ (2)

The average correlations were determined in circular regions on the detector
with centers $q=0.0025,0.0035,0.0045$ and $0.0055$ Å-1, where each region had
a width of $\Delta q=0.001$ Å-1. Ideally, the two definitions above should
provide the same results, but it was found that they differ for low scattering
samples because of spurious parasitic scattering. Here, a data set was
discarded due to low scattering intensity if the average relative difference
($ARD$) between $g_{2}^{(1T)}$ and $g_{2}^{(2T)}$ was larger than 1 %.
Otherwise, a mean of the two was used as the true auto-correlation
$g_{2}(q,\tau)$. To ensure spatial homogeneity of the samples, the repeated
scans (at different positions) were used to determine the mean
$g_{2}(q,\tau)$. However, if one repeated experiment showed significantly
higher correlations than its counterpart ($ARD>0.5$ %), it was considered an
outlier and removed from the analysis. The final square first order auto-
correlation used for analysis
$c(q,\tau)=g_{1}^{2}(q,\tau)=(g_{2}(q,\tau)-g_{\infty})/\beta$ is found by
subtracting the baseline $g_{\infty}=g_{2}(q,\tau>5~{}s)$ and scaling with the
speckle contrast $\beta$ ($\beta=\beta^{*}+1-g_{\infty}$, where $\beta^{*}$ is
found through initial calibration with a test membrane).  The time-dependence
of the dynamics over the course of each experiment was also assessed with the
two-time correlation function to exclude any heterogeneous dynamics or other
time-dependent phenomena.

### XPCS analysis

Each resulting correlation curve $c(\tau)$ at a certain $q$ was assumed a form
of:

$c(\tau)=\sum_{i}G^{(d)}_{i}\exp(-2\Gamma^{(d)}_{i}\tau),$ (3)

where the discrete distribution $G^{(d)}_{i}$ is found through a regularized
inverse Laplace transform using the CONTIN-algorithm implemented in MATLAB
2020b and adapted from the code by Marino 51. Here, we assumed 50 log-spaced
modes in range $\Gamma^{(d)}_{i}\in[75,750000]$ s-1 with constraints
$G^{(d)}_{i}\geq 0\quad\forall i$, $G^{(d)}_{i}\rightarrow 0$ at the extreme
values of $\Gamma^{(d)}_{i}$ and that $\sum_{i}G^{(d)}_{i}=1$, _i.e._ that
$c(\tau\rightarrow 0)=1$. The CONTIN-algorithm requires a regularization
parameter for additional smoothing of the distribution 34, 35, which was here
chosen to 0.05. However, an assessment of the influence of the regularization
parameter (provided in Figs. S3 and S4 in SI) showed that the exact value had
little effect on the extracted dynamic modes for $RP\leq 0.1$.

The mean relaxation rate $\bar{\Gamma}$ of the decay is determined from the
distribution through:

$\ln\bar{\Gamma}=\sum_{i}G^{(d)}_{i}\ln\Gamma^{(d)}_{i}.$ (4)

The $q$-dependence of $\bar{\Gamma}$ was then used to find the Brownian
diffusion coefficient $D_{B}$ through fitting $\bar{\Gamma}(q)=D_{B}q^{2}$ and
the dynamic scaling exponent $n$ through $\bar{\Gamma}(q)=(Dq)^{1/n}$.

Through a cumulant expansion 36, the initial decay ($\tau\rightarrow 0$) of
the auto-correlation can be approximated with $\ln
c(\tau)\approx-2\bar{\Gamma}_{0}\tau+k_{2}^{2}\tau^{2}+O(\tau^{3})$, where
$\bar{\Gamma}_{0}$ is the mean relaxation rate at $\tau=0$ (note that
$\bar{\Gamma}_{0}\neq\bar{\Gamma}$). The polydispersity index ($PDI$) is
determined through $PDI=k_{2}^{2}/\bar{\Gamma}_{0}^{2}$. By doing a low $\tau$
expansion of $c(\tau)$ both from the CONTIN result and the cumulant expression
and comparing terms, it can be found that:

$\displaystyle\sum_{i}G_{i}^{(d)}$ $\displaystyle=$ $\displaystyle 1$ (5)
$\displaystyle\sum_{i}G_{i}^{(d)}\Gamma_{i}^{(d)}$ $\displaystyle=$
$\displaystyle\bar{\Gamma}_{0}$ (6)
$\displaystyle\sum_{i}G_{i}^{(d)}\left(\Gamma_{i}^{(d)}\right)^{2}$
$\displaystyle=$
$\displaystyle\left(\bar{\Gamma}_{0}^{2}+\frac{k_{2}^{2}}{2}\right)$ (7)

and the polydispersity index can thus be found from the discrete relaxation
rate distribution $G^{(d)}_{i}$:

$PDI=\frac{2\left(\sum_{i}G_{i}^{(d)}\left(\Gamma_{i}^{(d)}\right)^{2}-\bar{\Gamma}_{0}^{2}\right)}{\bar{\Gamma}_{0}^{2}}.$
(8)

The polydispersity index was seen not seen to vary significantly with $q$ and
a mean value was taken over all $q$ for comparisons in Fig. 2. The relative
quantities $D_{B}/D_{0}$ and $PDI/PDI_{0}$ in Fig. 2 were obtained by using
the values $D_{0}$ and $PDI_{0}$ from experiments of free diffusion of GNPs in
the solvent (see Fig. S8 in SI for details).

### Extraction of dynamic modes through lognormal deconvolution

The lognormal deconvolution is used to create a continuous distribution
$G^{(c)}(\Gamma)$ of the resulting discrete distribution $G^{(d)}_{i}$ and was
done by assuming the form with three dynamic modes in the continuous
distribution ($\Gamma_{1}$, $\Gamma_{2}$, and $\Gamma_{3}$):

$G^{(c)}(\Gamma)=\sum_{k=1}^{3}G^{(c)}_{k}\exp(-a_{k}[\ln(\Gamma)-\ln(\Gamma_{k})]^{2}),$
(9)

where magnitude $G^{(c)}_{k}$, relaxation rate $\Gamma_{k}$ and parameter
$a_{k}$ of each mode $k=1,2,3$ are found through fitting. The fitted magnitude
$G^{(c)}_{k}$ is then used to determine the contribution of each mode
$G^{(d)}_{k,i}$ to the original discrete distribution $G^{(d)}_{i}$, where
$G_{i}^{(d)}=\sum_{k=1}^{3}G^{(d)}_{k,i}$. The magnitudes of each mode $G_{k}$
is found by $G_{k}=\sum_{i}G^{(d)}_{k,i}$, which fulfills
$\sum_{k=1}^{3}G_{k}=1$.

### Numerical simulations of confined Brownian motion

The numerical simulation model is a 2D version of the model presented by Ernst
et al. 45 to simulate subdiffusive Brownian motion in mucus. Here a square
domain of $L=5$ $\mu$m was divided into $n=L/l$ quadratic subdomains (”cells”)
with cell size $l=50$ nm (analysis of cell sizes $l=25$ & 75 nm is provided in
Fig. S9 in SI), where $N=100$ particles are initially randomly placed in the
domain. The position of a particle, $\bm{r}(t)$ equals to integration of the
stochastic differential equation (SDE)
$\dot{\bm{r}}(t)=\sqrt{2D_{0}}\bm{\zeta}(t)$, where $\bm{\zeta}(t)$ is a
stationary Markov process satisfying $\langle\bm{\zeta}(t)\rangle=\bm{0}$. The
diffusion constant is set according to $D_{0}=k_{B}T/(6\pi\eta R)$, where
radius is $R=100$ nm, temperature is $T=293$ K, and viscosity is $0.042$ Pa s
(viscosity of PG at this temperature 52). If a particle crosses a boundary of
a cell, it will be reflected with probability $r=1-P\sqrt{\Delta\tau}$, where
$P>0$ is the membrane permeability with unit s-1/2 and $\Delta\tau=1$ ms is
the time step. The simulation is running for a total of 6000 iterations
(_i.e._ 6 s). The mean-square-displacement is calculated as the ensemble
average over all particles
$MSD(\tau)=\langle[\bm{r}(t=\tau)-\bm{r}(t=0)]^{2}\rangle$. It can be noted
that in the extreme case of $P=0$, the particle is contained within its pore
with $MSD(\tau\rightarrow\infty)=l^{2}/6$ 45.

### Simulated XPCS experiments

The simulated speckle patterns were created by a pre-defined detector with
$200\times 200$ pixels covering a range of $q_{x},q_{y}\in[-0.0065,0.0065]$
Å-1. The dynamic structure factor at a given time $S(q,t)$ was calculated
using a discrete Fourier transform of the particle positions $\bm{r}(t)$:

$S(q,t)=\left|\sum_{i=1}^{N}\exp[j\bm{q}\cdot\bm{r}_{i}(t)]\right|^{2},$ (10)

with $j=\sqrt{-1}$. The intensity $I(q,t)=P(q)S(q,t)$ is found by
multiplication of the form factor for a sphere 53:

$P(q)=\left(\frac{\sin(qR)-qR\cos(qR)}{(qR)^{3}}\right)^{2}$ (11)

using radius $R=100$ nm. The $q$-regions for evaluating intensity
autocorrelations in the simulated detector images were chosen with centers
$q=0.0025,0.003...0.006$ Å-1 and width of $\Delta q=0.005$ Å-1. The resulting
autocorrelations are then analyzed through the same steps as the experimental
data by deriving the square first order correlation
$c(\tau)=g_{1}^{2}(q,\tau)=(g_{2}(q,\tau)-g_{\infty})/\beta$, where
$g_{2}(q,\tau)=g_{2}^{(2T)}(q,\tau)$ and $\beta$ is found through
extrapolation to $\tau=0$.

The authors acknowledge the financial support from the National Science
Foundation (DMR-1808690), Wallenberg Wood Science Center (WWSC), the Alf de
Ruvo Memorial Foundation and the Hans Werthen Foundation. The authors also
acknowledge experimental assistance by C. Zhan, J. Tian, S. Chodankar and Y.
Zhang as well as helpful discussions with F. Lundell, J. Sellberg, M.
Nordenström and L. Wågberg.

The experiments were performed at the CHX beamline (11-ID) and the LIX
beamline (16-ID) of the National Synchrotron Light Source II, a U.S.
Department of Energy (DOE) Office of Science User Facility operated for the
DOE Office of Science by Brookhaven National Laboratory under Contract No. DE-
SC0012704.

The Center for BioMolecular Structure (CBMS) is primarily supported by the
National Institutes of Health, National Institute of General Medical Sciences
(NIGMS) through a Center Core P30 Grant (P30GM133893), and by the DOE Office
of Biological and Environmental Research (KP1607011).

The Supplementary Material includes: a movie showing a 3D rendering of the
experimental conditions in the experiment (available at
https://play.kth.se/media/SI_XPCS_paper_conditions/0_f8cgf96f), rheology of
the CNF dispersions, a theoretical assessment of pore sizes in random fiber
networks, assessment of various fitting methods of the XPCS data, assessment
of radiation damage, assessment of gravity-induced sedimentation, sizes and
dynamics of GNPs in pure solvent, and results from numerical simulations for
other cell sizes.

## References

  * Gardel et al. 2004 Gardel, M.; Shin, J. H.; MacKintosh, F.; Mahadevan, L.; Matsudaira, P.; Weitz, D. A. Elastic behavior of cross-linked and bundled actin networks. _Science_ 2004, _304_ , 1301–1305
  * Gitai 2005 Gitai, Z. The new bacterial cell biology: moving parts and subcellular architecture. _Cell_ 2005, _120_ , 577–586
  * Jansen et al. 2018 Jansen, K. A.; Licup, A. J.; Sharma, A.; Rens, R.; MacKintosh, F. C.; Koenderink, G. H. The role of network architecture in collagen mechanics. _Biophysical journal_ 2018, _114_ , 2665–2678
  * Ling et al. 2018 Ling, S.; Kaplan, D. L.; Buehler, M. J. Nanofibrils in nature and materials engineering. _Nature Reviews Materials_ 2018, _3_ , 1–15
  * Wegst et al. 2015 Wegst, U. G.; Bai, H.; Saiz, E.; Tomsia, A. P.; Ritchie, R. O. Bioinspired structural materials. _Nature materials_ 2015, _14_ , 23–36
  * Klemm et al. 2011 Klemm, D.; Kramer, F.; Moritz, S.; Lindström, T.; Ankerfors, M.; Gray, D.; Dorris, A. Nanocelluloses: a new family of nature-based materials. _Angewandte Chemie International Edition_ 2011, _50_ , 5438–5466
  * Li et al. 2021 Li, T.; Chen, C.; Brozena, A. H.; Zhu, J.; Xu, L.; Driemeier, C.; Dai, J.; Rojas, O. J.; Isogai, A.; Wågberg, L., et al. Developing fibrillated cellulose as a sustainable technological material. _Nature_ 2021, _590_ , 47–56
  * Rosén et al. 2020 Rosén, T.; Hsiao, B. S.; Söderberg, L. D. Elucidating the Opportunities and Challenges for Nanocellulose Spinning. _Advanced Materials_ 2020, 2001238
  * Håkansson et al. 2014 Håkansson, K. M.; Fall, A. B.; Lundell, F.; Yu, S.; Krywka, C.; Roth, S. V.; Santoro, G.; Kvick, M.; Wittberg, L. P.; Wågberg, L., et al. Hydrodynamic alignment and assembly of nanofibrils resulting in strong cellulose filaments. _Nature communications_ 2014, _5_ , 1–10
  * Mittal et al. 2018 Mittal, N.; Ansari, F.; Gowda. V, K.; Brouzet, C.; Chen, P.; Larsson, P. T.; Roth, S. V.; Lundell, F.; Wagberg, L.; Kotov, N. A., et al. Multiscale control of nanocellulose assembly: transferring remarkable nanoscale fibril mechanics to macroscale fibers. _ACS nano_ 2018, _12_ , 6378–6388
  * Walther et al. 2011 Walther, A.; Timonen, J. V.; Díez, I.; Laukkanen, A.; Ikkala, O. Multifunctional high-performance biofibers based on wet-extrusion of renewable native cellulose nanofibrils. _Advanced Materials_ 2011, _23_ , 2924–2928
  * Voisin et al. 2017 Voisin, H.; Bergström, L.; Liu, P.; Mathew, A. P. Nanocellulose-based materials for water purification. _Nanomaterials_ 2017, _7_ , 57
  * Sharma et al. 2020 Sharma, P. R.; Sharma, S. K.; Lindström, T.; Hsiao, B. S. Water Purification: Nanocellulose-Enabled Membranes for Water Purification: Perspectives (Adv. Sustainable Syst. 5/2020). _Advanced Sustainable Systems_ 2020, _4_ , 2070009
  * Rostami et al. 2021 Rostami, J.; Gordeyeva, K.; Benselfelt, T.; Lahchaichi, E.; Hall, S. A.; Riazanova, A. V.; Larsson, P. A.; Ciftci, G. C.; Wågberg, L. Hierarchical build-up of bio-based nanofibrous materials with tunable metal–organic framework biofunctionality. _Materials Today_ 2021,
  * Czaja et al. 2006 Czaja, W.; Krystynowicz, A.; Bielecki, S.; Brown Jr, R. M. Microbial cellulose—the natural power to heal wounds. _Biomaterials_ 2006, _27_ , 145–151
  * Hickey and Pelling 2019 Hickey, R. J.; Pelling, A. E. Cellulose biomaterials for tissue engineering. _Frontiers in bioengineering and biotechnology_ 2019, _7_ , 45
  * Eskilson et al. 2020 Eskilson, O.; Lindström, S. B.; Sepulveda, B.; Shahjamali, M. M.; Güell-Grau, P.; Sivlér, P.; Skog, M.; Aronsson, C.; Björk, E. M.; Nyberg, N., et al. Self-Assembly of Mechanoplasmonic Bacterial Cellulose–Metal Nanoparticle Composites. _Advanced Functional Materials_ 2020, _30_ , 2004766
  * Yao et al. 2018 Yao, J.; Ji, P.; Wang, B.; Wang, H.; Chen, S. Color-tunable luminescent macrofibers based on CdTe QDs-loaded bacterial cellulose nanofibers for pH and glucose sensing. _Sensors and Actuators B: Chemical_ 2018, _254_ , 110–119
  * Hamedi et al. 2014 Hamedi, M. M.; Hajian, A.; Fall, A. B.; Hakansson, K.; Salajkova, M.; Lundell, F.; Wagberg, L.; Berglund, L. A. Highly conducting, strong nanocomposites based on nanocellulose-assisted aqueous dispersions of single-wall carbon nanotubes. _ACS nano_ 2014, _8_ , 2467–2476
  * Amiralian et al. 2020 Amiralian, N.; Mustapic, M.; Hossain, M. S. A.; Wang, C.; Konarova, M.; Tang, J.; Na, J.; Khan, A.; Rowan, A. Magnetic nanocellulose: A potential material for removal of dye from water. _Journal of hazardous materials_ 2020, _394_ , 122571
  * Wu et al. 2012 Wu, C.-N.; Saito, T.; Fujisawa, S.; Fukuzumi, H.; Isogai, A. Ultrastrong and high gas-barrier nanocellulose/clay-layered composites. _Biomacromolecules_ 2012, _13_ , 1927–1932
  * Liu et al. 2011 Liu, A.; Walther, A.; Ikkala, O.; Belova, L.; Berglund, L. A. Clay nanopaper with tough cellulose nanofiber matrix for fire retardancy and gas barrier functions. _Biomacromolecules_ 2011, _12_ , 633–641
  * An et al. 2018 An, B. W.; Heo, S.; Ji, S.; Bien, F.; Park, J.-U. Transparent and flexible fingerprint sensor array with multiplexed detection of tactile pressure and skin temperature. _Nature communications_ 2018, _9_ , 1–10
  * Okita et al. 2011 Okita, Y.; Fujisawa, S.; Saito, T.; Isogai, A. TEMPO-oxidized cellulose nanofibrils dispersed in organic solvents. _Biomacromolecules_ 2011, _12_ , 518–522
  * Wang et al. 2019 Wang, R.; Rosen, T.; Zhan, C.; Chodankar, S.; Chen, J.; Sharma, P. R.; Sharma, S. K.; Liu, T.; Hsiao, B. S. Morphology and flow behavior of cellulose nanofibers dispersed in glycols. _Macromolecules_ 2019, _52_ , 5499–5509
  * Rosén et al. 2020 Rosén, T.; He, H.; Wang, R.; Zhan, C.; Chodankar, S.; Fall, A.; Aulin, C.; Larsson, P. T.; Lindström, T.; Hsiao, B. S. Cross-Sections of Nanocellulose from Wood Analyzed by Quantized Polydispersity of Elementary Microfibrils. _ACS nano_ 2020, _14_ , 16743–16754
  * Geng et al. 2017 Geng, L.; Peng, X.; Zhan, C.; Naderi, A.; Sharma, P. R.; Mao, Y.; Hsiao, B. S. Structure characterization of cellulose nanofiber hydrogel as functions of concentration and ionic strength. _Cellulose_ 2017, _24_ , 5417–5429
  * Rosén et al. 2021 Rosén, T.; Wang, R.; He, H.; Zhan, C.; Chodankar, S.; Hsiao, B. S. Understanding Ion-Induced Assembly of Cellulose Nanofibrillar Gels through Shear-Free Mixing and In Situ Scanning-SAXS. _Nanoscale advances_ 2021, _3_ , 4940–4951
  * Geng et al. 2018 Geng, L.; Mittal, N.; Zhan, C.; Ansari, F.; Sharma, P. R.; Peng, X.; Hsiao, B. S.; Söderberg, L. D. Understanding the mechanistic behavior of highly charged cellulose nanofibers in aqueous systems. _Macromolecules_ 2018, _51_ , 1498–1506
  * Rosén et al. 2021 Rosén, T.; Wang, R.; He, H.; Zhan, C.; Chodankar, S.; Hsiao, B. S. Shear-Free Mixing to Achieve Accurate Temporospatial Nanoscale Kinetics through Scanning-SAXS: Ion-Induced Phase Transition of Dispersed Cellulose Nanocrystals. _Lab on a Chip_ 2021, _21_ , 1084–1095
  * Ogston 1958 Ogston, A. The spaces in a uniform random suspension of fibres. _Transactions of the Faraday Society_ 1958, _54_ , 1754–1757
  * Chatterjee 2012 Chatterjee, A. P. A simple model for the pore size distribution in random fibre networks. _Journal of Physics: Condensed Matter_ 2012, _24_ , 375106
  * Madsen et al. 2010 Madsen, A.; Leheny, R. L.; Guo, H.; Sprung, M.; Czakkel, O. Beyond simple exponential correlation functions and equilibrium dynamics in x-ray photon correlation spectroscopy. _New Journal of Physics_ 2010, _12_ , 055001
  * Provencher 1982 Provencher, S. W. CONTIN: a general purpose constrained regularization program for inverting noisy linear algebraic and integral equations. _Computer Physics Communications_ 1982, _27_ , 229–242
  * Andrews et al. 2018 Andrews, R. N.; Narayanan, S.; Zhang, F.; Kuzmenko, I.; Ilavsky, J. Inverse transformation: unleashing spatially heterogeneous dynamics with an alternative approach to XPCS data analysis. _Journal of applied crystallography_ 2018, _51_ , 35–46
  * Koppel 1972 Koppel, D. E. Analysis of Macromolecular Polydispersity in Intensity Correlation Spectroscopy: The Method of Cumulants. _The Journal of Chemical Physics_ 1972, _57_ , 4814–4820
  * Tracy and Pecora 1992 Tracy, M. A.; Pecora, R. Synthesis, characterization, and dynamics of a rod/sphere composite liquid. _Macromolecules_ 1992, _25_ , 337–349
  * Kluijtmans et al. 2000 Kluijtmans, S. G.; Koenderink, G. H.; Philipse, A. P. Self-diffusion and sedimentation of tracer spheres in (semi) dilute dispersions of rigid colloidal rods. _Physical Review E_ 2000, _61_ , 626
  * Pryamitsyn and Ganesan 2008 Pryamitsyn, V.; Ganesan, V. Dynamics of probe diffusion in rod solutions. _Physical review letters_ 2008, _100_ , 128302
  * Metzler et al. 2014 Metzler, R.; Jeon, J.-H.; Cherstvy, A. G.; Barkai, E. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. _Physical Chemistry Chemical Physics_ 2014, _16_ , 24128–24164
  * Pal et al. 2021 Pal, A.; Kamal, M. A.; Zinn, T.; Dhont, J. K.; Schurtenberger, P. Anisotropic dynamics of magnetic colloidal cubes studied by x-ray photon correlation spectroscopy. _Physical Review Materials_ 2021, _5_ , 035603
  * Reiser et al. 2020 Reiser, M.; Hallmann, J.; Möller, J.; Kazarian, K.; Orsi, D.; Randolph, L.; Rahmann, H.; Westermeier, F.; Stellamanns, E.; Sprung, M., et al. Nanoscale Rigidity in Cross-Linked Micelle Networks Revealed by XPCS Nanorheology. _arXiv preprint arXiv:2010.08267_ 2020,
  * Kwaśniewski et al. 2014 Kwaśniewski, P.; Fluerasu, A.; Madsen, A. Anomalous dynamics at the hard-sphere glass transition. _Soft Matter_ 2014, _10_ , 8698–8704
  * Gowda et al. 2022 Gowda, V. K.; Rosén, T.; Roth, S. V.; Söderberg, L. D.; Lundell, F. Nanofibril Alignment during Assembly Revealed by an X-ray Scattering-Based Digital Twin. _ACS nano_ 2022, _16_ , 2120–2132
  * Ernst et al. 2017 Ernst, M.; John, T.; Guenther, M.; Wagner, C.; Schaefer, U. F.; Lehr, C.-M. A model for the transient subdiffusive behavior of particles in mucus. _Biophysical journal_ 2017, _112_ , 172–179
  * Xu et al. 2021 Xu, Z.; Dai, X.; Bu, X.; Yang, Y.; Zhang, X.; Man, X.; Zhang, X.; Doi, M.; Yan, L.-T. Enhanced Heterogeneous Diffusion of Nanoparticles in Semiflexible Networks. _ACS nano_ 2021, _15_ , 4608–4616
  * Saito et al. 2007 Saito, T.; Kimura, S.; Nishiyama, Y.; Isogai, A. Cellulose nanofibers prepared by TEMPO-mediated oxidation of native cellulose. _Biomacromolecules_ 2007, _8_ , 2485–2491
  * Katz et al. 1984 Katz, S.; Beatson, R. P., et al. The determination of strong and weak acidic groups in sulfite pulps. _Svensk papperstidning_ 1984, _87_ , 48–53
  * Lumma et al. 2000 Lumma, D.; Lurio, L.; Mochrie, S.; Sutton, M. Area Detector Based Photon Correlation in the Regime of Short Data Batches: Data Reduction for Dynamic X-ray Scattering. _Review of Scientific Instruments_ 2000, _71_ , 3274–3289
  * Fluerasu et al. 2007 Fluerasu, A.; Moussaïd, A.; Madsen, A.; Schofield, A. Slow dynamics and aging in colloidal gels studied by x-ray photon correlation spectroscopy. _Physical Review E_ 2007, _76_ , 010401
  * Marino 2007 Marino, I.-G. rilt - Regularized Inverse Laplace Transform. https://www.mathworks.com/matlabcentral/fileexchange/6523-rilt, 2007; MATLAB Central File Exchange. Retrieved June 3, 2021
  * Sagdeev et al. 2017 Sagdeev, D. I.; Fomina, M. G.; Abdulagatov, I. M. Density and viscosity of propylene glycol at high temperatures and high pressures. _Fluid Phase Equilibria_ 2017, _450_ , 99–111
  * Feigin et al. 1987 Feigin, L.; Svergun, D. I., et al. _Structure analysis by small-angle X-ray and neutron scattering_ ; Springer, 1987; Vol. 1