To include quadrupole field in the model of quantitative susceptibility mapping (QSM).

Magnetic field, $$$H$$$, can be written in terms of scale potential as $$$H=-\nabla\Phi$$$, based on governing magnetostatic equations, $$$ \nabla\times H=0 $$$ and $$$ \nabla\cdot H=-\nabla\cdot M $$$. $$$H$$$ field is expanded up to 2nd order as follows.$$ \Phi(x)=-\frac{1}{4\pi}\int\frac{\nabla'\cdot M(x')}{|x-x'|}d^{3}x'\approx\frac{1}{4\pi}\left(\frac{m\cdot x}{r^{3}}+\sum_{i=1}^{3}\sum_{j=1}^{3}\frac{1}{2}Q_{ij}\frac{x_{i}x_{j}}{r^{5}}\right) $$where, $$$ Q_{ij}=\int d^{3}x'\left(3M{}_{i}x_{j}'+3M_{j}x_{i}'-2\delta_{ij}M\cdot x'\right) $$$.

The first and second term corresponds to dipole and quadrupole interaction, respectively. Only the z component is of interest, $$$ H_{z}=-\partial_{z}\Phi $$$.$$ H_{z}(r)=d*\chi+\sum_{i=1}^{3}g_{i}*Q_{i} $$where $$$ d=\frac{1}{4\pi}\frac{3(\hat{m}\cdot x)^{2}-r^{2}}{r^{5}} $$$, $$$ g_{1}=-\frac{3}{4\pi}\frac{x\left(x^{2}+y^{2}-4z^{2}\right)}{r^{7}} $$$, $$$ g_{2}=-\frac{3}{4\pi}\frac{y\left(x^{2}+y^{2}-4z^{2}\right)}{r^{7}} $$$, $$$ g_{3}=-\frac{3}{4\pi}\frac{z\left(3x^{2}+3y^{2}-2z^{2}\right)}{r^{7}} $$$, $$$ \chi\equiv\int d^{3}x'M_{z} $$$, $$$ Q_{1}\equiv\int d^{3}x'M_{z}x' $$$, $$$Q_{2}\equiv\int d^{3}x'M_{z}y' $$$, and $$$ Q_{3}\equiv\int d^{3}x'M_{z}z' $$$. The ratio of quadrupole moment to susceptibility corresponds to the center of susceptibility within a voxel: $$$ CS_{1}=Q_{1}/\chi $$$, $$$ CS_{2}=Q_{2}/\chi $$$, and $$$ CS_{3}=Q_{3}/\chi $$$.

The corresponding Fourier Transform of $$$H_{z}(r)$$$ is $$ H_{z}(k)=D\cdot\chi+\sum_{i=1}^{3}G_{i}\cdot Q_{i} $$.

Where $$$ D=\frac{1}{3}-\frac{k_{z}^{2}}{k^{2}} $$$, $$$ G_{1}=-\frac{i\pi k_{x}}{k^{2}}\left(k^{2}-5k_{z}^{2}\right) $$$, $$$G_{2}=-\frac{i\pi k_{y}}{k^{2}}\left(k^{2}-5k_{z}^{2}\right) $$$, and $$$ G_{3}=-\frac{i\pi k_{z}}{k^{2}}\left(3k^{2}-5k_{z}^{2}\right) $$$.

The formulation above is based on the assumption that the main field is along $$$z$$$. If the main field is along $$$\hat{n}=[n_{1},\,n_{2},\,n_{3}]^{T} $$$, dipole ($$$D$$$), quadrupole function ($$$G_{i}$$$), and quadrupole moment ($$$Q_{i}$$$) are changed correspondingly as follows. $$ D(k)\rightarrow D(Rk),G(k)\rightarrow G(Rk),Q\rightarrow RQ $$

where $$$R$$$ is rotation matrix, $$$ \hat{z}=R\hat{n} $$$.

$$ R=\left[\begin{array}{ccc}n_{3}+\frac{n_{2}^{2}}{1+n_{3}} & -\frac{n_{1}n_{2}}{1+n_{3}} & -n_{1}\\-\frac{n_{1}n_{2}}{1+n_{3}} & n_{3}+\frac{n_{1}^{2}}{1+n_{3}} & -n_{2}\\n_{1} & n_{2} & n_{3}\end{array}\right] $$

Multipole orientation measurements are needed because of four unknown parameters.

$$ f_{j}(k)=D(R^{j}k)\cdot\chi(k)+\sum_{i=1}^{3}G_{i}(R^{j}k)\cdot R^{j}Q_{i}(k) $$

where $$$f_{j}(k)$$$ is the local field in k-space, and $$$R^{j}$$$ denotes the $$$j^{'}$$$th rotation.

A 1% agarose phantom was prepared which contains four balloons with different concentrations of gadolinium solution. For data acquisition, all experiments were conducted on a 3.0T clinical Siemens scanner with multi-echo gradient echo sequence. The resolution was set to 1 × 1 × 1 mm, TR to 27 ms, bandwidth to 260 Hz/pixel, flip angle to 15º. Four TEs were used: 6.2, 11.5, 16.9, and 22.2 ms. The phantom were rotated with respect to the magnetic field. Eight different rotations were acquired at the COSMOS optimal angles (0, 60, 120 and 270º) in the xy and yz planes¹: $$$ \hat{n} = $$$ [0 -1 0], [-0.804 -0.574 -0.028], [-0.912 0.373 -0.004], [0.986 0.156 0.001], [0.008 -0.004 -1.007], [0.830 -0.040 -0.474], [0.812 -0.048, 0.515], and [-0.603 0.038 0.807]. For the post-processing, phase unwrapping and background field removal² were used to obtain local field ($$$f_{i}$$$). COSMOS result was calculated with the multiple orientation data. The susceptibility ($$$\chi$$$) and quadrupole moment ($$$Q_{i}$$$) were calculated solving the last equation of $$$f_{j}(k)$$$ with Conjugate Gradient method (residue 0.05). Fig. 1 shows the susceptibilities of four balloons based on COSMOS (Fig. 1a) and based on our quadrupole model (Fig. 1b). The COSMOS susceptibility result were 2.7732 $$$\pm$$$ 0.1195 ppm, 1.5241 $$$\pm$$$ 0.0368 ppm, 0.7794 $$$\pm$$$ 0.0186 ppm, and 0.4375 $$$\pm$$$ 0.0182 ppm, respectively. Our quadrupole model susceptibility result were 2.7504 $$$\pm$$$ 0.2272 ppm, 1.3689 $$$\pm$$$ 0.0973 ppm, 0.6472 $$$\pm$$$ 0.0581 ppm, and 0.3387 $$$\pm$$$ 0.0663 ppm, respectively, which agreed with the dipole model within ~0.15ppm. Fig. 2 shows quadrupole moment in the unit of ppm×mm (Fig. 2a-c) and the ratio of quadrupole moment to dipole moment in the unit of mm, which corresponds to the center of susceptibility (Fig. 2d-f). The quadrupole effect was noticeable at the boundaries where Gd only fills part of the voxels. The color in Fig. 2d-f indicates the degree of how far away the center of susceptibility locates from the voxel center, e.g. white indicates the susceptibility center locates further than the voxel center along each direction.

In this study, we established the theoretical frame work to calculate quadrupole moment, and verified it with phantom data. With this method, we quantitatively measure both susceptibility and non-uniform distribution of magnetization inside voxel, e.g. center of susceptibility distribution in a voxel.

In the inclusion of quadrupole moments, the number of variables substantially worsens the condition of the inverse problem. Multiple orientations may be sampled to estimate parameters in the quadruple model, but misregistation among orientations are inevitable and cause artifacts and errors in the reconstructed parameter maps. Alternative approach using regularization to address the ill posedness may alleviate misregistration issue and allow translation into in vivo applications.

In conclusion, the magnetic source to field forward problem may include quadruple in additional to dipole, and quantitative susceptibility mapping of dipole and quadruple moments is feasible.