Ansgar Adler^{1}, Jost M. Kollmeier^{2}, Nick Scholand^{1}, Sebastian Rosenzweig^{1}, Yong Wang^{3}, and Martin Uecker^{1,4}

^{1}Institute for Diagnostic and Interventional Radiology, (UMG) University Medical Center Göttingen, Göttingen, Germany, ^{2}Biomedical NMR, Max Planck Institute for Biophysical Chemistry, Göttingen, Germany, ^{3}Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany, ^{4}DZHK (German Centre for Cardiovascular Research), Göttingen, Germany

The lattice Boltzmann method (LBM) is a versatile numerical technique for simulating complex fluid dynamical systems and beyond. Here, we first describe an extension of the LBM to flow systems in external magnetic fields. The model is verified numerically with several benchmarks and shown to correctly predict signal changes caused by in-flow effects in a simple flow experiment.

In the conventional LBM, the distribution function $$$f(\vec{r}, t)$$$ of a fluid particle is discretized in physical space, time, and velocity space $$$f_i(\vec{r}, t) = f(\vec{r}, \vec{c}_i, t)$$$, with $$$\vec{r}$$$ the location on a spatial grid, $$$t$$$ the time step, and $$$i$$$ the index of the discretized velocity vectors $$$\vec{c}_i$$$. The statistical state of the internal spin degree of freedom of an ensemble of spin 1/2 particles, as the H$$$^1$$$ nucleus, is described by the quantum mechanical density operator $$$D$$$, which can be decomposed into Pauli matrices according to

$$ D = \frac{1}{2}\left( \sigma_0 + \vec a \cdot \vec \sigma \right) = \frac{1}{2}\left( \sigma_0 + a_x \sigma_x + a_y \sigma_y + a_z \sigma_z \right),$$

where $$$\vec \sigma$$$ is the vector of the three Pauli matrices $$$\sigma_x, \sigma_y, \sigma_z$$$, and $$$\sigma_0$$$ the identity matrix, and $$$\vec a$$$ the Bloch vector [3]. The macroscopic magnetization is related to the quantum mechanical expectation value of the spin which is given by the density matrix multiplied with the spin observable: $$$\gamma \frac{\hbar}{2} {\bf Tr} D \vec{\sigma}$$$. We propose a model of flowing particles with spin 1/2 using a matrix-valued semi-classical velocity-spin distribution $$$D(\vec r, t)$$$ and discretize it by

$$D_i(\vec r, t) := \frac{1}{2} \sum_j f_{i,j}(\vec r, t) \sigma_j, $$

where $$$j$$$ denotes the polarization direction $$${0,x,y,z}$$$. The positivity of the coefficient $$$f_i$$$ in the classical LBM corresponds to positive definiteness of the matrix $$$D_i$$$. Summing over all velocity directions recovers the density $$$\rho$$$, the bulk velocity $$$\vec u$$$, and the components of the macroscopic magnetization $$$\vec M$$$ at each lattice node according to

$$\rho = \sum_i f_{i,0}, \qquad \rho \vec u = \sum_i f_{i,0} \vec c_i, \quad \text{ and } \quad M_j = \sum_i f_{i,j} \text{ .}$$

The LBM simulates the discretized Boltzmann equation [1] in a streaming and collision step. The streaming step works on each $$$j$$$ independently and transports each component $$$f_{i,j}$$$ to the next grid point in the direction $$$\vec c_i$$$. The BGK collision step is

$$ f'_{i,j'} \leftarrow \sum_j B_{j',j} \left( f_{i,j} - \frac{f_{i,j} - f^{eq}_{i,j}}{\tau} \right), $$

with the relaxation time $$$\tau$$$, post-collision distribution $$$f'_{i,j}$$$ and equilibrium function $$$f^{eq}_{i,j}$$$. Here, we assume a common equilibrium distribution for the velocities defined by $$$\rho$$$ and $$$\vec u$$$ and independent from the magnetization as in the standard LBM. The matrix $$$B_{j,j'}$$$ implements a step of the time-discretized Bloch equations, modelling relaxation and the effect of all external fields.

Several canonical flows, e.g. the Pouiseuille flow, backward facing step flow and Karman vortex street, were considered. Numerical results with the proposed model agree well with the benchmarks. Furthermore, we showed that for zero flow the magnetization recovers the Ernst equation, i.e. steady-state magnetization [2] for a spoiled FLASH sequence. Finally, the agreement with the Bloch-Torrey equation was verified numerically.

Developing laminar pipe flow with different flow rates, controlled by the pump-level, was considered. To verify the prediction of the signal change caused by the in-flow of spins, flow through a cross-sectional slice was measured experimentally with a phase-contrast FLASH sequence. Slice thickness and flip angle were varied. We then simulated the flow and magnetization of such an experiment in a plane adapting the velocity in the center to the one determined by phase-contrast MRI. In the intersection line of the measured slice and the simulated plane, the signal intensities were compared between experiment and simulation (Fig. 1).

Fig. 2 shows the measured signal and velocity profiles for varying pump-levels in comparison to the numerical ones from simulations. Flat profiles in the center imply that the flows were developing. Considering the discrepancy in the flow profiles in the 3D experiment and the 2D simulation, analysis is restricted to a region of homogenous flow in the center.

Fig. 3 depicts the numerical and experimental mean signal levels for the region marked gray in Fig. 2 for varying flip angle and slice thickness.

[1] **Krüger T., Halim K.**, *The Lattice Boltzmann Method*, Springer-Verlag, 2017.

[2] **Liang Z. P., Lauterbur P. C.**, *Principles of Magnetic Resonance Imaging: A Signal Processing Perspective*, IEEE Press Series on Biomedical Engineering. Wiley, 1999.

[3] **Blum K.**, *Density Matrix Theory and Applications*, Plenum Press, 1981.

[4] **Puiseux T., Sewonu A., Moreno R., Mendez S., and Nicoud F.**, *Numerical simulation of 4D Flow MRI.* ISMRM Virtual Conference & Exhibition 2020, In Proc. Intl. Soc. Mag. Reson. Med. 28; 1042, 2020.

[5] **Jurczuk K., Kretowski M., Bellanger J., Eliat P., Saint-Jalmes H., Bézy-Wendling J.**,* Computational modeling of MR flow imaging by the lattice Boltzmann method and Bloch equation. *In Magnetic Resonance Imaging Vol 31; 7, 2013.

Figure 1: Experimental setup to verify the prediction of the signal change caused by the in-flow of spins. A cross-sectional slice through a small pipe placed in a channel (bottom picture) was measured using a flow sequence. The magnitude images obtained from the slice (yellow) were compared in a line (green) intersecting the simulated plane (red). This line is also indicated in the magnitude image.

Figure 2: The numerical (dashed) and experimental (solid) velocity (a) and signal (b) profiles for different pump-levels are plotted for the intersecting line inside the pipe. The region in the center is used in further comparisons (gray).

Figure 3: The mean signal values from the center region are shown in dependency of flip angle (a) and slice thickness (b) for different pump levels for numerical simulation (dashed) and experiment (solid). The error bar shows the standard deviation. The left figure is for a slice thickness of 6 mm and the right is for a flip angle of 20°.