Dynamic magnetic resonance scattering

Dynamic light scattering is a popular technique to determine the size distribution of small particles in the sub micrometer region. It operates in reciprocal space, by analyzing the signal fluctuations with the photon auto correlation function. Equally, pulsed field gradient magnetic resonance is a technique generating data in the reciprocal space of the density distribution of an object. Here we show the feasibility of employing a magnetic resonance imaging system as a dynamic scattering device similar to dynamic light scattering appliances. By acquiring a time series of single data points from reciprocal space, analogue to dynamic light scattering, we demonstrate the examination of motion patterns of microscopic particles. This method allows the examination of particle dynamics significantly below the spatial resolution of magnetic resonance imaging. It is not limited by relaxation times and covers a wide field of applications for particle or cell motion in opaque media. Dynamic light scattering is a method used to examine the dynamics and the size distribution of submicrometer particles and operates in reciprocal space. The authors apply the fundamentals of this technique to magnetic resonance imaging in order to examine particle motion below the spatial resolution in short measurement time.

T he principles of magnetic resonance (MR) build the foundation of a variety of imaging applications 1 . Generally, images are not acquired directly in the spatial domain but in the reciprocal space, the k-space. To reconstruct an image, a sufficient number of data points in k-space has to be sampled. Thus imaging particle motion in real time (e.g., particles at a size of several micrometers dispersed in a fluid) is hardly feasible due to the need of high spatial and temporal resolution. In this paper we show that it is possible to quantify statistical parameters of particle motion if the temporal evolution of only a single k-space point is studied. This is made possible by transferring the principles known from dynamic light scattering (DLS) to MR.
In the following, a brief introduction is given to the very basics of motion sensitive magnetic resonance measurements. For a more in-depth treatment of the topic the reader is referred to the common textbooks 2,3 . During a standard MR experiment the magnetization vector of the sample, initially colinear with the main magnetic field B 0 = (0,0,B z ), is tipped into the transverse plane by applying a resonant radio-frequency pulse. The resulting transverse magnetization precesses around B 0 with the typical Larmor frequency ω = γB 0 , where γ is the gyromagnetic ratio. Omitting relaxation effects, the measurement signal can be written as the sum of all local precessing magnetization vectors: where I(r) is proportional to the local spin density and represents the MR image. Using constant magnetic field gradients G = ∇B z , the Larmor frequency ω can be varied linearly for each spatial dimension. Eq. (1) can then be written as: where the vector k is defined as: Hence, the MR image can be obtained from the Fourier transform of the measurement signal, which is acquired in kspace 4 . To perform the Fourier transformation, the measurement signal must be sampled for a sufficient range in k-space by varying either the gradient amplitude, or the gradient lobe duration of a gradient pulse. In this way the k-space can be filled line by line as shown in the simple example of a 2D-gradient echo sequence in Fig. 1. The signal amplitude is controlled by transversal relaxation processes, such as T 2 and T Ã 2 relaxation, and the speed of rebuilding longitudinal magnetization, i.e. the T 1 relaxation. Minimum acquisition time in k-space is thus limited by constraints from relaxation times and the need for sufficient signal-to-noise ratio. Imaging small particle motion in real time (e.g., particles at a size of several micrometers dispersed in a fluid) can thus be very challenging, if not impossible. In contrast to optical methods, snapshot-like fast acquisition methods are not possible, which for example would allow for single particle tracking in real time.
However, MR acquisition in k-space allows the manipulation and quantification of the signal phase which opens up a range of possibilities to quantify motion. A very common and fundamental technique is the so-called Pulsed Gradient Spin Echo (PGSE) technique as shown in Fig. 2a 2,5 . The PGSE technique uses a Hahn-Spin-Echo extended with two short symmetric gradient pulses, one before and one after the 180°radio-frequency-pulse (which effectively inverts the polarity of the first gradient pulse) 6 . Since the motion encoding gradients are symmetric they do not influence the spatial encoding; i.e. spatial and motion encoding can be combined, thus building the fundamental technique of MR diffusion imaging. In the PGSE signal, diffusing particles lead to a signal attenuation. In the narrow gradient pulse approximation the amplitude of the center of the Spin-Echo depends on the diffusion coefficient D, the gradient pulse spacing Δ and the q-value and can be described as: The q-value specifies the motion encoding capability of the gradient pulses and is defined as: where δ is the gradient pulse duration and G its amplitude. Eq. (5) corresponds to the definition of the k-value from Eq. (3) 7 . However, the definition of a q-value is only meaningful in the narrow gradient pulse approximation and in its application with symmetric motion-encoding gradients. In that case, the q-value on its own spans a reciprocal space to the space of displacements R = r(t + Δ) − r(t), which can be seen when using the propagator formalism to describe the center echo amplitude: The averaged propagator P is the Fourier transform of the MR signal when acquired in q-space just as the MR image is the Fig. 1 Gradient-Echo sequence: Pulsed field gradients generate k-space encoding according to Eq. (3). For simplicity slice-selection is omitted. A magnetic resonance (MR)-image can be reconstructed by applying a 2D-Fourier transform to the measurement data, which is acquired in k-space Fourier transform of the MR signal when acquired in k-space 2,8 .
Sampling the MR signal in both spaces, i.e., using k-space and qspace encoding, respectively, allows to quantify diffusion and flow in a pixel-wise manner. However, this leads to a long measurement time since both k-space and q-space have to be sampled sufficiently.
Recent publications were able to show, that by breaking the symmetry between the two diffusion-encoding gradients the diffusion experiment can be shifted from being a scattering experiment to being an imaging experiment, enabling to define the shape of the diffusion boundaries which are in general not accessible by a pure scattering experiment 9,10 . Classical resolution limits in MR imaging, as defined by a limited magnitude of the maximum k-space vector, can thus be circumvented.
In this paper, we show that ensemble motion can also be quantified entirely without using motion encoding gradients (as done with q-space encoding), by using Dynamic MR Scattering (DMRS). The phase-sensitive acquisition of a time series of Structure functions for a constant particle drift: Stucture functions for brownian motion: k-space encoded signals already allows extraction of the parameters of the ensemble motion. k-space encoding signifies a harmonic spatial modulation. When particles move, i.e. when the spin-density distribution changes, the phase of the k-space signal also changes in time. The higher the k-value, the more sensitive the signal phase is to translational motion. For a certain k-space position the fluctuating signal can then be evaluated by calculating the autocorrelation function. This approach generalizes the idea of interpreting dynamic data acquisition in k-space as a scattering experiment, without the necessity of image reconstruction in analogy to concepts used in dynamic light scattering (DLS) and differential dynamic microscopy [11][12][13] .

Results
Theoretical description of dynamic magnetic resonance scattering. To illustrate the principle, we assume a sample made of a number of identical point-like particles dispersed in a viscous fluid. It is irrelevant whether the source for the signal is found in the particles or in the fluid, provided that the particles generate sufficient signal contrast to the fluid background to create variations in the spatial signal distribution. In the following discussion, for convenience, the particles are assumed to generate the signal instead of the surrounding fluid (which generally would be the case for MR). The spatial distribution of the signal density I(r,t) for N particles can then be written as the convolution 12 : δ½r À r j ðtÞ : ð7Þ I 0 represents the local spatial distribution of the particle signal analogously to the scattering potential of a single particle as known from optics. The fluctuating k-space signal for the moving particles, which is generated by applying magnetic field gradients, can be derived from (7) by applying the spatial Fourier transform: where S 0 (k) is the Fourier transform of I 0 (r). The autocorrelation function of the k-space signal corresponds to the field correlation function in dynamic light scattering: Γðk; tÞ ¼ hSðk; t′ÞS Ã ðk; t′ þ tÞi t′ For particles moving independently from each other the temporal averaging as denoted by the angular brackets in (9) vanishes for correlations between displacements of different particles: Assuming N identical particles, the autocorrelation function of the k-space signal becomes: The normalized version of the field correlation function for uncorrelated and identical particles can then be written as: In DLS experiments, g(k,t) is accessible through the Siegert relation 12 . With the phase-sensitive MR measurements, the field correlation function can be directly calculated from the time course of a single k-space point. When assuming identical particles, Eq. (12) directly corresponds to the signal attenuation in a PGSE experiment given in Eq. (6). That is to say, instead of dephasing the measured signal by applying diffusion-encoding gradients and subsequently exploiting its effect on the signal attenuation, one can simply repeat single-point k-space sampling and evaluate the temporal correlation in k-space.
The temporal distance between the diffusion gradients Δ in Eq. (6) corresponds to t in the field correlation function (Eq. (12)). Albeit in PGSE experiments, Δ must not significantly exceed T 2 (or T 1 in the case of stimulated echo experiments), t in Eq. (12) can virtually be extended to infinity with DMRS experiments, since there is no theoretical upper limit for the repetition time (T R ) in a MR experiment. This enables to also quantify very slow ensemble motion.
From dynamic light scattering, the course of the field correlation function Eq. (12) is well known for Brownian motion as an exponential decay 12 : with the diffusion coefficient D and for a constant particle drift as an oscillating function: with the drift velocity v. In cases where most of the signal energy is contained in the static background signal, it is more favorable to examine functions describing the statistics of the signal differences in order to cancel out the background signal 14 . One of these functions is the averaged mean square of the signal differences in k-space also known as the structure function: The time course of the structure function is closely related to the field correlation function. The normalized version of the structure function can be written as: where <½f denotes the real part of f. The structure function allows to extract the sample dynamics from the fluctuating kspace signal, even in cases of strong background signal. In case of a constant translational motion F(k, t) will be an oscillating function as g(k, t). Determining the oscillation angular frequency ω = v⋅k gives access to the velocity v. In case of Brownian motion, the exponential decay from Eq. (13) is also transfered to Eq. (16).
Fitting an exponential decay 1 − e −t/τ to the structure function with the relaxation time τ = 1/Dk 2 allows to determine the diffusion constant.
Simulations. To demonstrate the feasibility of the proposed technique, MR measurements are simulated using a Bloch simulation combined with a simulation for samples producing particle drift and Brownian motion. First, the signal for a sample based on glass spheres dispersed in water (ρ = 2500 kg m −3 , Ø = 100 μm) is simulated. According to Stokes Law the constant gravitational drift velocity of the particles is set to v = 8.2 mm s −1 . The MR sequence producing k-space data comprises an excitation pulse followed by gradient pulse during which different data points are sampled in the reciprocal space as shown in Fig. 2b. The MR sequence is then repeated 100 times with a repetition time T R = 10 ms such that the measurement of one experiment covers a time window of 1 s. Although sampling the time-course of just a single k-space data point would be sufficient, 128 sampling points are acquired (each one represents a different k-space position) with a sampling frequency of 100 kHz. The structure function for each single k-space value contains the full statistical information about the particle motion as shown in Fig. 2c, d. Figure 3a shows the time course of the normalized structure function for three different k-space values. The oscillating progression of the normalized structure function clearly reflects the translational movement of the particles at a constant velocity along the spatial encoding gradient. Identifying the peaks in the discrete spectra as shown in Fig. 3b, allows to determine the drift velocity. The velocities determined from the peak show a good agreement with the predefined drift velocity of 8.2 mm s −1 of the particles. Due to the limited spectral resolution, determining the peak position is less accurate than directly fitting a sinusoid. Figure 3c shows the results for the decorrelation time τ = 1/(v⋅k) plotted as a function of k, when applying sinusoidal fits to the different structure functions. The drift velocity is found to be v = 8.18 mm s −1 by fitting the decorrelation time as a function of k.
In a second step we simulate Brownian motion of 1 μm particles dispersed in water (viscosity η = 1.0 ⋅ 10 −3 N s m −2 ; temperature T = 300 K) using a 3D-random walk. According to the Stokes-Einstein-relation the diffusion coefficient is set to D = 4.39 ⋅ 10 −13 mm 2 s −1 . Figure 3d shows the structure function plotted for three different k-values and fitted to the exponential decay as given in Eqs. (13) and (16). The results of the simulation show a good agreement with the given diffusion coefficient (D 1 = MR measurements. In order to demonstrate the feasibility of DMRS, we examine sedimenting glass microspheres (ρ = 2500 kg m −3 , Ø = 100 μm) on a Bruker 17.6 T scanner equipped with a 1 T m −1 gradient insert (see Fig. 4a). The temporal evolution of the k-space-data is sampled using a 1D fast gradient echo sequence 15 . The direction of the frequencyencoding (k x -encoding) is chosen to be aligned with the drift direction. Figure 4b shows a photograph of the sedimenting glass spheres. For reference, a sequence of MR images is acquired, additionally. Due to the limited spatial resolution and primarily due to the limited temporal resolution the reconstructed MR images reveal no details about the particle size and distribution, as shown in Fig. 4c.
In contrast, with the DMRS-data, i.e., with the time series of the k-space data, the imprint of the ensemble motion becomes visible, which can be seen in the corresponding structure functions and their spectra as shown in Fig. 4d, e. Figure 4f shows the evaluation of the decorrelation time τ = 1/(v⋅k) for different k-values. Fitting this plot results in a drift velocity of v = (12 ± 1) mm s −1 . By analyzing high-resolution video frames of the sedimenting glass spheres, the mean drift-velocity is found to be v = (11.4 ± 0.9) mm s −1 , which is in good agreement with the MR experiments.
In a second experiment, the velocity distribution of rising airbubbles in a water tube is examined (see Fig. 4g). DMRS-k-space data are acquired again using a 1D fast gradient-echo sequence (i.e. without phase-encoding) on a Magnetom Skyra-3T (Siemens, see Fig. 4h). With video analysis and single air-bubble tracking, we could identify a range of velocities between 150 and 250 mm s −1 . In Fig. 4i, j, the corresponding structure functions and its spectra are plotted for three different k x -values. For each k x -value, the velocity at the peak-value of the spectrum is given. Determining the decorrelation time τ = 1/(v⋅k) for different k-values (see Fig. 4k) allows to specify a mean value for the most prominent drift velocity of the bubbles (v = (189 ± 7) mm s −1 ). The width of the spectral distribution shows a good agreement with the range of velocities that were found by video analysis. Both examples show the capability of the presented method to examine particle dynamics. The modulation of the field correlation function reflects the dynamics of the particles and is governed by the magnitude of the corresponding k-vector.

Discussion
Contrary to the PGSE technique, a contrast between the particles and the background is needed. This means that the field of applications can be found primarily in analyzing the dynamics of heterogenous systems, such as colloidal suspensions or emulsions, by directly examining the dynamics of particles or droplets. In respect thereof DMRS can be considered as a complement to the existing MR methods probing translational dynamics, which are widely based on the Stejskal-Tanner technique.
To improve the contrast between continuous phase and dispersed phase, it is also conceivable to use contrast agents such as ultrasmall superparamagnetic iron oxid (USPIOS), thus even particles at cell size (i.e., several micrometers) can provide sufficient contrast to the background. It is important to mention that it is not necessary to spatially resolve the particle distribution, as is done for instance with molecular imaging. Furthermore, the particle dynamics must provide sufficient signal fluctuation in order to be captured with the autocorrelation function. For Brownian motion, for instance small particles (<1 μm), theoretically can be examined easier with less gradient performance. The slower the particle dynamics, the higher the gradient strength necessary to capture the ensemble motion. However, the required gradient strength can be reduced by extending T R , i.e., the temporal spacing between the k-space sampling points. In case of a strong background signal or a dense particle distribution, it is also important that the receiver chain is able to resolve the k-space signal fluctuation. In that case, it might be useful to only sample the relevant higher k-space sampling points in order to suppress the high amplitude signal from the k-space center.
An interesting field of application might be the characterization of emulsion systems in natural as well as synthetic products. In this regard, PGSE-NMR is widely used to determine the dropletsize distribution (DSD) [16][17][18] . The examinations are mainly based on the quantification of the self-diffusion of the dispersed liquid within the droplet wall. However, the droplets themselves also undergo Brownian motion, leading to additional spin dephasing during the encoding gradients, which might hamper the accuracy of the DSD examination. DMRS is not sensitive to self-diffusion but primarily relies on a spatial and temporal signal contrast between the continuous and the dispersed phase, albeit the gradient performance is not needed to spatially resolve particles or droplets. Thus, mathematical models evaluating the DSD can be significantly simplified since self-diffusional processes are not accounted for by the measurement process.
Additionally, the ability to study the DSD with low-cost/lowfield MR bench-top spectrometers has led to the widespread use of PGSE-NMR in suitable emulsions and colloidal suspensions when optical methods cannot be used, as in opaque media. However, bench-top MR systems have to rely on the possibility to separate signals from the continuous and the dispersed phase by relaxation or diffusion-weighting to isolate each phase for self-diffusion measurement. In our approach, this clear separation is needed only to the extent that sufficient contrast to noise is provided to map the particle dynamics, making our presented application particularly interesting for low-cost bench-top systems.
When using PGSE-NMR, the range of droplets able to be sized is limited by the duration of the temporal separation Δ (see Fig. 2a) of the encoding gradient 18 . PGSE-NMR allows quantification of droplet sizes between 0.2 and 100 μm, which limits the range of applications for industrial bench-top MR systems 17 . To circumvent this limitation sophisticated MR techniques have been proposed, such as the use of dipolar demagnetization fields 19 . With DMRS, the inter-pulse separation time Δ from PGSE-NMR is translated to the repetition time T R , which has no upper limit dictated by the relaxation times T 2 (in case of a Spin-Echo-measurements) or T 1 (in case of the use of stimulated echoes). Thus, DMRS also allows to monitor the dynamics at a low diffusion constant D, since there is no theoretical upper limit for the repetition time T R , as long as the measurement environment can be kept stable for the whole measurement duration. Summarized, the use of DMRS to monitor chemical and physical processes in multi-phasic liquids could help to characterize and control corresponding industrial processes, e.g. hydrocracking or emulsion ripening/flocculation/coagulation/aggregation/etc.
A further field of possible applications is the investigation of granular flow and sedimentation processes. These processes have already been investigated using MR techniques. For example Seymour used motion sensitizing gradients to study stochastic and deterministic properties of the fluid motion in a granular flow 20 . In contrast to this method, the presented method, provides means to study the motion of non MR-sensitive solid particles itself, which are carried along by the fluid and show different dynamic behavior. Although this information is also encoded in the flow propagator it seems more difficult to extract this information from the data. In comparison, the proposed method can be specifically tailored to yield this information. Thus, our proposed method could provide a simple way to complete the dynamic picture of granular flow processes. For slower processes like sedimentation direct imaging has also been proposed 21 . However, this requires a sufficiently slow sedimentation process and hardware capable of spatially resolving the solid components as well as their motion. The presented approach eases these conditions. As imaging is not necessary, even processes on the time scale of a few T R or below can be resolved. Furthermore, the ensemble dynamics of particles smaller than the spatial resolution of the acquired k-space line can be investigated. This may reduce hardware requirements and thus make these experiments more feasible. These examples are not limited to the application of solid non-MR-sensitive objects in a fluid, but can be used to study these properties in any multiphasic system where the phases can be distinguished by MR methods, thus generating the necessary contrast between the phases to detect the temporal changes. This can be chemical shift, e.g. oil droplets, or different nuclei, e.g. perfluor carbon emulsions.
Another area that may profit from the presented method is the characterization of multi-phase flow in trickle or fixed bed reactors which is important for reactor efficiency [22][23][24][25][26] . For example, changing from bubbly or trickle flow to pulsed flow could result in an observable change/loss of the spatio-temporal signal-correlation and thus lead to a change in the apparent relaxation of the correlation function. Again, this does not require an air-liquid interface if both phases can be separated by MR and the reactor bed can be distinguished from all other phases. If multiple k-space points are acquired for instance, a combination with the Bayesian analysis by Holland et al. might provide a better characterization of the transporting liquid or the transport pathways in the reactor bed 27 . Also, a combination with T Ã 2 weighted imaging as described by Arbabi et al. is possible if data is sampled at different echo times 28 . This may also provide a possible way of creating different contrasts for different phases.

Methods
Derivation of the field correlation function. The spatio-temporal dynamics of a homogenous system of particles can be captured with the Van Hove correlation function 29 : δðR À r n ð0Þ þ r m ðtÞÞ The function splits into two parts. A self correlation part G S (R, t) giving the probability that a particle moved a distance R in time t: δðR À r n ð0Þ þ r n ðtÞÞ and a distinct part giving the probability to find at time t a particle j at a distance R from a place where there was a different particle i at time t = 0: In case of statistically independent particles G D vanishes and G(R, t) = G S (R, t). Further for identical particles G can be written as: Calculating the spatial Fourier transform of G S (R, t): Gðk; tÞ ¼ Z d 3 Re ikÁR hδðR À r n ð0Þ þ r n ðtÞÞi ¼ e ikÁðr n ðtÞÀr n ð0ÞÞ shows that the Fourier transform of the space time correlation function equals the first order field correlation function from Eq. (12): In case of Brownian motion superimposed with a constant translational motion such as sedimentation for example, the particle flux J can be described by Fick's first law of diffusion with an additional term describing the constant translational motion at the velocity v: Further, if the number of particles is constant in time, the conversation of particles can be expressed by the continuity equation Substituting the flux J from Eq. (23) yields The probability density function G S (R, t) obeys Fick's law in the same manner as the particle concentration c so that it is reasonable to assume that G S (R, t) is also the solution to Eq. (25) 11 : Calculating the spatial Fourier transform and using Eq. (22) yields: ∂ ∂t gðk; tÞ À ½ik Á vgðk; tÞ ¼ Àk 2 Dgðk; tÞ : Assuming g(k, 0) = 1 gives the solution for the field correlation function g(k, t): gðk; tÞ ¼ e ikÁvt e Àk 2 Dt : ð28Þ The first part represents the oscillating influence of the constant translational motion on the correlation function. The second part causes an exponential decay due to the uncorrelated Brownian motion.
MR simulations. In general, the spin density distribution is not found to be the exact Fourier transform of the k-space data in MR. This only holds true for a completely static spin density distribution. To demonstrate, that the principles of DLS can still be applied in MR, we have simulated moving particles in a spatial encoding gradient field, as it can be found in typical MR imaging sequences. The phase evolution of the particles was calculated as follows: with k Ã i ¼ γ R G i tdt. k * incorporates the effect of first order phase accumulation due to constant velocities during spatial encoding.
For the simulations, we assumed point-like, non-interacting particles. The simulation calculates the individual phase accumulation of the transverse magnetization for each particle according to Eq. (29). The k-space signal was generated by adding the signal contribution of each particle. We simulated a fully spoiled 1D gradient-echo-sequence with gradients applied only in readout direction.
Simulation of particle drift. A total of 500 particles (Ø = 1 μm) randomly positioned in a sample volume of the size of 3.0 × 3.0 × 2.0 cm 3 filled with water (viscosity η = 1.0 ⋅ 10 −3 N s m −2 ) were examined under the constant drift velocity of v ¼ d 2 gðρ glas À ρ H 2 O Þ=ð18ηÞ ¼ 8:2 mm s À1 with g = 9.81 m s −2 , ρ H 2 0 ¼ 1000 kg m À3 and ρ glas = 2500 kg m −3 . The direction of the resultant gravitational drift was set to coincide with the gradient direction. The readout-gradient strength was set to G x = 0.08 T m −1 . The magnetization was assumed to always start with the same amplitude at the sequence entrance and hence T 1 relaxation effects were omitted. Data acquisition was performed as shown in Fig. 2b with a sampling frequency of ν sampling = 100 kHz and 128 sampling points during each readoutgradient. Each sampling point represents a different k and k * -space value respectively. The sequence was repeated N t = 100 times at a repetition time of T R = 10 ms such that the measurement of one experiment covers a time window of 1 s. The whole experiment was repeated N av = 100 times and the k-space signal differences where averaged according to Eq. (15).
Simulation of Brownian motion. The sample size was chosen to be (3.0 × 3.0 × 2.0) μm 3 with 5000 dissolved particles in water. The particle diameter was set to 1 μm. The Brownian motion was simulated with a 3D-random walk with a time resolution of τ = 10 −5 s and a 1D stepsize of hxi ¼ ffiffiffiffiffiffiffiffiffiffiffi 2DΔτ p . According to the Stokes-Einstein-relation D was calculated to be D = 4.39 ⋅ 10 −13 m 2 s −1 . To adapt the MR sequence to the slow diffusion process, the following changes were made to the sequence: G x = 0.23 T m −1 , repetition time T R = 1 s, sampling frequency ν sampling = 10 kHz, number of sequence repetitions N t = 100, number of averages: N av = 250. MR measurement at 3 T. The experiments to study the velocity distribution of rising air-bubbles in water were conducted with a Magnetom Skyra-3T (Siemens) with a maximum gradient field strength of 0.045 T m -1 . The temporal evolution of the k-space-data was acquired using a FLASH-sequence without k y -encoding, i.e., without phase encoding. Under each readout-gradient 80 k x -sampling-points were acquired (actual gradient strength: 0.73 mT m −1 ; gradient duration: 1.3 ms). The whole experiment was averaged N av = 383 times. Further parameters were: flip-angle α = 15°, slice-thickness 20 mm; field-of-view (x): 0.5 m; resolution (x): 6.25 mm.
Simulations of particles undergoing Brownian motion and translation simultaneously. Additionally to the results presented in the main manuscript, particles undergoing simultaneously Brownian motion and translational motion were simulated. Particles with a diameter of 1 μm (D = 4.39 ⋅ 10 −13 m 2 s −1 ) and a constant drift velocity of 1 μm s −1 were examined. The same pulse sequence was applied as for the determination of the pure Brownian motion. The results for the structure functions and its spectra can be seen in Fig. 5a, b for three selected k-values. The structure function is now a damped harmonic oscillation. The oscillating part represents the translational motion, whereas the exponential decay is caused by the loss of spatial correlation due to the Brownian motion.
Fðv; D; k; tÞ ¼ 1 À Ae ÀDk 2 t cos v 0 Ákt ð Þ; ð30Þ Figure 5b, c shows the results when fitting a damped cosine wave as given in Eq. (30) to the k-dependant structure functions. Both, the diffusion coefficient and the drift velocity, can be extracted correctly from the structure functions. For high k-values the fast loss of correlation of the k-space signal in time blurs the characteristics of the structure functions and leads to increased imprecision of the fitted parameters as can be seen Fig. 5c, d. Suitable k-values to probe the drift can be found in the range of k low < k < k high with k low ≈ 2π/(T A v) and k high ≈ f s π/v, where T A is the temporal sampling window for the field correlation function, f s = 1/T R is the sampling frequency and v the drift velocity. The lower boundary is the spectral resolution due to the finite T A while the upper boundary is the Nyquist sampling theorem. For Brownian motion k should be on the order of magnitude of a few k low % ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1=DT A p , with the diffusion coefficient D. This ensures sufficient dynamics in the structure function for parameter fitting. For high k-values the temporal sampling bandwidth has to be elevated accordingly such that the dynamic range of the exponential decay is sampled sufficiently, otherwise the exponential fit increasingly delivers biased diffusion coefficients as can be seen in Fig. 5b for high k-values. For a given sampling rate, one finds k high % ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1=DT R p . In Fig. 5e, f, the correlation times τ are plotted as functions of k, for the particle drift (τ = 1/(v⋅k) see Fig. 5e) and the Brownian motion (τ = 1/Dk 2 see Fig. 5f), respectively. Unstable fits (i.e., when the relative error for the fit parameter exceeded 20 %) where excluded from the plot.
As shown in Fig. 5a, the structure functions change significantly with increasing k. For low k-values, nearly a pure oscillation is observed. With large k the attenuating part dominates. It can be expected that for high k-values, only an exponential curve is observable. This indicates that, in limiting cases, the structure function is either dominated by the diffusion process or by the translation. For simplicity in the following v⋅k = vk is assumed. Rescaling Eq. (30) by k = k T κ and τ = k T vt shows Unstable fit results i.e. those exceeding a relative error for the fit parameter of 20 % where excluded from the plot with k T = v/D as the k-value where the relaxation rate of the damping equals the oscillation frequency. If κ 2 τ ) κ τ ) κ ) 1 (large k-values), damping dominates the time evolution of F(κ, τ) and F κ)1 ðκ; τÞ % 1 À Ae Àκ 2 τ ð32Þ holds, which is the structure function for Brownian motion (Eqs (13) and (16) in the manuscript). The opposite case of weak damping is found if κ 2 τ ( κ τ ) κ ( 1. In this case for short times τ ( 1=κ the damping can be neglected and one finds Fðκ; τÞ % 1 À Acos κτ ð Þ; ð33Þ which corresponds to Eqs (14) and (16) of the manuscript for a constant drift. For τ≳1=κ observable damping of the oscillation would be present. For unscaled times this translates to T A ( D= κ 2 v 2 ð Þ¼1= k 2 D ð Þ. The experimental application of these limits shall be investigated exemplarily for the 1 μm particles from the simulation above. In this case one finds k T ≈ 2 × 10 6 m −1 . Starting with the fast damping limit, k-values significantly higher than k T have to be realized. With k ¼ γ R Gdt and assuming a pulse duration of 1 ms a gradient strength of G ) 7:5T m À1 would be necessary to generate k ) k T . Even if extending the gradient-pulse duration, such high k-values can only be realized with highly specialized MR hardware, e.g., dedicated diffusion probes. However, the weak-damping limit with k ( k T in the scenario of the drifting 1 μm particles can easily be met with current MR imaging systems as demonstrated in the experiment.

Data availability
The data that supports the findings of this study are available upon request from the corresponding author.