Modeling the Magnetization Dynamics for Large Ensembles of Immobilized Magnetic Nanoparticles in Multi-dimensional Magnetic Particle Imaging

inst1]Hannes Albers
inst2,inst3]Tobias Knopp
inst2,inst3]Martin Möddel
inst2,inst3]Marija Boberg
inst1,inst4]Tobias Kluth

[inst1]organization=Center for Industrial Mathematics, University of Bremen,
addressline=Bibliothekstrasse 5,

[inst2]organization=Section for Biomedical Imaging, University Medical Center Hamburg-Eppendorf,
addressline=Lottestrasse 55,

[inst3]organization=Institute for Biomedical Imaging, Hamburg University of Technology,
addressline=Am Schwarzenberg-Campus 1,

[inst4]organization=Center for Optimization and Approximation, University of Hamburg,
addressline=Bundesstrasse 55,

Magnetic nanoparticles (MNPs) play an important role in biomedical applications including imaging modalities such as magnetic resonance imaging (MRI) and magnetic particle imaging (MPI). The latter one exploits the non-linear magnetization response of a large ensemble of magnetic nanoparticles to magnetic fields which allows determining the spatial distribution of the MNP concentration from measured voltage signals. The image-to-voltage mapping is linear and described by a system matrix. Currently, modeling the voltage signals of large ensembles of MNPs in an MPI environment is not yet accurately possible, especially for liquid tracers in multi-dimensional magnetic excitation fields. As an immediate consequence, the system matrix is still obtained in a time consuming calibration procedure. While the ferrofluidic case can be seen as the typical setting, more recently immobilized and potentially oriented MNPs have received considerable attention. By aligning the particles magnetic easy axis during immobilization one can encode the angle of the particle's magnetic easy axis into the magnetization response providing a sophisticated benchmark system for model-based approaches. In this work we address the modeling problem for immobilized and oriented MNPs in the context of MPI. We investigate a model-based approach where the magnetization response is simulated by a Néel rotation model for the particle's magnetic moments and the ensemble magnetization is obtained by solving a Fokker-Planck equation approach. Since the parameters of the model are a-priori unknown, we investigate different methods for performing a parameter identification and discuss two different models: One where a single function vector is used from the space spanned by the model parameters and another where a superposition of function vectors is considered. We show that our model can much more accurately reproduce the orientation dependent signal response when compared to the equilibrium model, which marks the current state-of-the-art for model-based system matrix simulations in MPI.

Néel rotation immobilized magnetic nanoparticles particle magnetization modeling magnetic particle imaging

0000 1111

0000 1111

§ INTRODUCTION

Magnetic nanoparticles (MNP) play an important role in various biomedical applications. They can be used for diagnosis in imaging applications and therapy <cit.>. Coupled with a drug they can be used for targeted drug delivery <cit.>, e.g. to treat a vessel occlusion. In hyperthermia applications the particles are heated locally by applying an oscillating magnetic field, which can be used in tumor therapy. The treatment using MNP usually requires a good localization of the particles, which requires tomographic imaging techniques. One indirect way is to use magnetic resonance imaging (MRI) where the MNP shorten the $T_2^*$ relaxation time and in turn can be identified by a negative contrast mechanism. However, in practice it can be challenging to differentiate tissue and particle signal using MRI. In 2005 an alternative method called magnetic particle imaging (MPI) was proposed that exploits the non-linear magnetization response of magnetic nanoparticles to external magnetic fields <cit.>. Since these particles saturate even for low field strength in the region of only few mT$\mu_0^{-1}$, MPI does provide positive contrast without any background signal from tissue.

MPI has shown to be a promising tool for a variety of medical application including vascular <cit.> and oncological applications <cit.>. In particular MPI shows great potential in the area of intervention imaging since one can use one tracer to visualize the vessel tree using a blood pool tracer <cit.> and at the same time track medical devices such as guide wires, catheters <cit.>, and stents <cit.>. In order to do so the particles are usually immobilized such that they can be permanently mounted on the device, e.g., by using a lacquer or by directly integrating the MNP into the polymer of a catheter. Quite recently, it has been shown that MPI can not only be used for tracking the position of MNP but also to determine the orientation <cit.>. This is done by immobilizing the particles within a strong static magnetic field, which due to the MNP's magnetic anisotropy results in a particle ensemble with parallel aligned magnetic easy axis, which influences the magnetization response during excitation. Using multi-contrast image reconstruction techniques <cit.> it is then possible to simultaneously determine the position and orientation of a sample with immobilized easy axis aligned particles <cit.>.

Here, the influence of the spatial orientation (parallel to the yellow line) of the easy axis of immobilized magnetic nanoparticles with parallel aligned magnetic easy axes on a 2D Lissajous type MPI system matrix is illustrated. Three selected frequency components with mixing factors $(m_x, m_y) \in \{ (1,1), (3,3), (5,5)\}$ for the $x$-channel are shown. On the left, system matrices measured for $\vartheta = \SI{0}{\degree},\SI{30}{\degree}, \dots, \SI{150}{\degree}$ are shown. On the right a fit of the equilibrium model (see M0 in Section <ref>) is shown for particle diameter $D=\SI{22}{nm}$. The later is a good approximation for the system matrix components with $0^\circ$ and $90^\circ$ easy axis orientation only. In all figures the complex values of the system matrix components are colored using a 2D colormap shown on the right, where saturation encodes the amplitude (absolute value) and color encodes the phase of a complex number. Amplitudes are normalized for each frequency component as indicated by the a.u. to account for their wide range. This normalization remains the same along the different orientations and the equilibrium model.

The image reconstruction problem in MPI consists of solving an inverse problem that is linear in the particle concentration <cit.>. The linear system obtained after discretization of the continuous imaging equation is ill-conditioned and requires proper regularization for a stable solution. Both the single-contrast and the multi-contrast imaging scenario can be formulated as a linear system and solved the same way using direct or iterative solvers. A prerequisite for image reconstruction is that the MPI system matrix $S\in \mathbb{C}^{K\times N}$ is known. The state-of-the art method for determining the system matrix $S$ is a fully data-driven calibration approach, where $S$ is measured column by column by measuring a small delta sample at each image voxel position. This procedure yields a physically accurate but noisy $S$ but it requires calibration times up to several days <cit.> and it is also very memory demanding, particularly when considering 3D multi-patch field of views.

For this reason there is a series of research works aiming to replace this measurement-driven approach by a model-driven approach. In <cit.> and <cit.> it has been shown that the simple equilibrium particle magnetization model (see, e.g., <cit.> for a definition) can be used to approximate the behavior of fluid MNP. The fully data-driven calibration method however, still outperformed these approaches for multi-dimensional MPI imaging protocols in terms of reconstruction image quality. In part this is due to the fact that the simplified equilibrium model relies on the assumption that the magnetization instantaneously follows the external magnetic field, i.e. no relaxation is present. In general, these assumptions are not fulfilled, especially when sufficiently fast changing fields are considered. In <cit.> it has been shown that a more sophisticated magnetization model yields much more accurate system matrices for fluid MNP that allows to describe relaxation effects that can be observed in real experiments.

In case of fast changing magnetic fields the dynamic behavior of the MNPs' magnetic moments is affected by Brownian (magnetic moment rotation due to rotation of the whole particle) and Néel (internal rotation of magnetic moment) mechanisms (see also <cit.>). For MPI, these new modeling approaches were initially suggested in <cit.> in terms of stochastic ordinary differential equations and the first solutions using the Fokker-Planck equation for Brownian rotation in a one-dimensional setup were presented by Yoshida and Enpuku <cit.> followed by several works in this direction, e.g., see <cit.> for further reading and <cit.> for a stronger focus on MPI excitation patterns in the context of Brownian and Néel rotations. The Néel rotation was exploited to formulate the already mentioned sophisticated magnetization models for ferrofluids in <cit.> where the authors also suggested an oriented easy axis model case which turned out to be not suitable for the fluid case. In particular, the Néel rotation mechanism is further influenced by the orientation of the particle's easy axis, which significantly affects the magnetization behavior, e.g., in the immobilized case for one-dimensional excitation, see <cit.>. Even if the immobilized case is less complex when compared to the ferrofluid case where a coupling between Brownian and Néel rotation needs to be taken into account, the modeling problem of the magnetization of large ensembles of MNPs has not been addressed for multi-dimensional excitation patterns so far and is the major focus of the present work.

§.§ Problem Statement

In this work we consider immobilized magnetic nanoparticles with parallel aligned magnetic easy axes and aim at modeling MPI system matrices for different spatial orientations of said axis using the Néel rotation model B2 derived in <cit.>. For orientation imaging <cit.>, one needs to determine several system matrices $S^{\vartheta_l}$, $l=1,\dots, L$, where each matrix has been acquired under a different orientation angle $\vartheta_l$. Within this work we focus on easy axis alignments which lie within a plane. While for single-contrast imaging the data-driven calibration approach for determining the system matrix can still be handled, in the case of multi-contrast imaging it becomes too tedious since the calibration time is increased by a factor of $L$. For example, if the measurement of a single system matrix requires one day, as reported in <cit.>, the multi-contrast calibration would require $L$ days in which the MPI system is occupied. A problem that is simply not an issue if a model-based approach is used to obtain the system matrices.

Fig. <ref> shows a subset of measured MPI system matrices measured with an immobilized and easy-axis aligned sample at different orientations. We note that the image resolution is chosen small ($11 \times 11$) such that the calibration time could be kept short. The measured matrices serve as the ground truth in the present work and the details of the experimental parameters are outlined in Section <ref>. The image shows three selected matrix rows that have been reshaped to represent the signal that a certain spatial position contributed to the voltage signal in Fourier space within a tomographic imaging experiment. This representation is very well known in MPI since the matrix rows for the equilibrium model can be represented using tensor products of Chebyshev polynomials <cit.> for Lissajous-type excitations. What is much more interesting now and has not been studied is the influence of parallel aligned magnetic easy axes on the structure of the MPI system matrices. Independent of the orientation of this axis one can see that a wave-hill structure is preserved in this case. However, what is also visible is that the wave hills are merging along certain directions which are global for each easy axis orientation. When changing the orientation, the merging angle changes accordingly, which provides the contrast mechanism for orientation imaging. Another observation is an increased amplitude when the orientation is closer to $45^\circ$ or $135^\circ$. Here, it is conspicuous that these angles are the orientations of the drive field vector in which a fast change of the length of the applied magnetic field vector is performed in the used cosine Lissajous-type applied magnetic field.

What one can also see is that the equilibrium model is not expressive enough to describe the observed effect at all and thus, this model is not suitable for modeling the system matrices in this particular multi-contrast experiment. The purpose of the present work is to use the model developed in <cit.> and calculate model-based system matrices that can describe the easy-axis dependent effects on the MPI system matrix shown in Figure <ref>. This is not only an important prerequisite for the development of more efficient calibration procedures but also for future investigations of different applied fields and their optimization with respect to multi-contrast imaging for magnetic easy axis orientation.

§ MATERIALS AND METHODS

Motivated by a general MPI imaging experiment where the particles are located in the MPI scanner and the change of the particle magnetization is measured using induction coils, we outline the general setting in the following. Let $\Omega \subset \R^3$ be the field-of-view (below we further restrict $\Omega$ to be the $xy$-plane). A voltage induced in a receive coil with sensitivity profile $p: \R^3 \rightarrow \R^3$ in m$^{-1}$ is given by

\begin{align}
\tilde{v} (t) &=- \mu_0 \int_\Omega  c(r) p(r)^\trans \frac{\partial}{\partial t}\bar{m}(r,t) \ \d r   \label{eq:complete-problem}
\end{align}

in V where $c: \Omega \rightarrow \R^+_0$ in  is the concentration of the magnetic nanoparticles and $\bar{m}: \R^3 \times [0,T] \rightarrow \R^3$ in $10^{-3}$ is the molar mean magnetic moment.
The mean magnetic moment $\bar{m}(r,t)$ depends on the applied magnetic field $H: \R^3 \times [0,T] \rightarrow \R^3$ in , which is usually a $T$-periodic function represented on the time interval $[0,T]$. In MPI a static selection field $H_\mathrm{SF}: \R^3 \rightarrow \R^3$ is combined with a dynamic drive field $H_\mathrm{DF}: \R^3 \times [0,T]\rightarrow \R^3$, i.e., $H(r,t)=H_\mathrm{SF}(r) + H_\mathrm{DF}(r,t)$. Note that the notation is motivated by the MPI point of view. Effectively, the selection field evaluated at discrete points encodes a variety of offset fields taken into account. The induced signal is analogously filtered to remove the main contribution of the direct feed-through. In good approximation it is represented by a convolution with a $T$-periodic filter kernel $a:\R \rightarrow \R$ yielding the signal $v=a\ast\tilde{v}$. Therefore, $v$ is $T$-periodic as well and can be expanded into a Fourier series with coefficients

\begin{align}
\hat{v}_k &= - \hat{a}_k \frac{\mu_0}{T} \int_\Omega  c(r) p(r)^\trans \int_{0}^{T}  \left( \frac{\partial}{\partial t}\bar{m}(r,t) \right){\textrm{e}}^{-2\pi i tk/T} \ \d t \ \d r \label{eq:MPI_model_fourier}
\end{align}

for $k\in\N_0$ and the conjugate complex $\hat{v}_{-k}=\overline{\hat{v}_k}$.
This formulation is commonly used, since the signal at the excitation frequencies is blocked using an analog band-stop filter prior to the signal digitization. With $s_k(r) := - \hat{a}_k\frac{\mu_0}{T} \int_{0}^{T} p(r)^\trans \frac{\partial}{\partial t}\bar{m}(r,t) {\textrm{e}}^{-2\pi i tk/T} \ \d t$
we can bring this into the standard notation

\begin{align}
\hat{v}_k &=\int_\Omega  c(r) s_k(r) \ \d r,  \label{eq:system_function_model}
\end{align}

where $s_k: \Omega \subset \R^3 \rightarrow \C$ is the system function, e.g., as previously discussed in Figure <ref>.
When measuring the system function at discrete positions $r_\nu \in \Omega, \nu=1,\dots,N$, and with a certain sampling rate in time, one obtains the MPI system matrix $S = \left( s_k(r_\nu) \right)_{k=0,\dots,K; \nu=1,\dots,N}$ storing the Fourier transformed measurements in the respective columns. In this work we  consider both modeled system matrices as well as measured system matrices that are obtained by a hybrid approach <cit.> exploiting the focus field in the MPI scanner.
In case of multiple receive channels, the corresponding system matrices can be stacked to arrive at a joint linear system of equation; see for example <cit.> for a formal definition.

§.§ Particle magnetization model

One important aspect in the model equation (<ref>) is the mean magnetic moment $\bar{m}$ of the ensemble of nanoparticles.
Over-simplified models like the equilibrium model (see, e.g., <cit.>) do not take into account the particles' anisotropy.
As we have the immobilized and oriented particle configuration in mind, we need a more sophisticated model and  exploit a Néel rotation model taking into account an uniaxial anisotropy.
For this we follow the Fokker-Planck equation approach for the Landau-Lifshitz-Gilbert equation as already described in <cit.>.
We determine the mean magnetic moment via the probability density function $f: \Omega \times \mathbb{S}^2 \times [0,T] \rightarrow \R_0^+$ which is the solution to the corresponding Fokker-Planck equation where $\mathbb{S}^2$ is the surface of the sphere in $\R^3$.
The mean is then given by

\begin{equation}\label{eq:mean_magmom_FP}
\bar{m}(r,t) = m_0 \int_{\mathbb{S}^2} m f(r,m,t) \ \d m
\end{equation}

where $f$ is the solution to the following specific case of a convection-diffusion equation on the sphere

\begin{equation}
\label{eq:FP-general}
\frac{\partial}{\partial t} f = \mathrm{div}_{\mathbb{S}^2}\left(\frac{1}{2\tau} \nabla_{\mathbb{S}^2} f \right) - \mathrm{div}_{\mathbb{S}^2}\left(b f\right)
\end{equation}

where $\tau >0$ is the relaxation time constant and the (velocity) field $b:\mathbb{S}^2 \times \R^3 \times \mathbb{S}^2 \rightarrow \R^3$ given by

\begin{align}
&b(m,H,n) = p_1 H \times m + p_2 (m\times H) \times m \notag \\ & \quad+ p_3 (n\cdot m) n \times m + p_4 (n\cdot m) (m\times n) \times m \label{eq:convection}
\end{align}

where $p_i\geq 0$, $i=1,\hdots,4$, are physical constants and $n\in \mathbb{S}^2$ is the easy axis of the particles.
Differentiation in terms of gradient $\nabla_{\mathbb{S}^2} $ and divergence $ \mathrm{div}_{\mathbb{S}^2}$ is considered with respect to the surface $\mathbb{S}^2$.
The Euclidean scalar product of $\R^3$ is denoted by $n\cdot m$.

A more detailed derivation of the Fokker-Planck equation from the Langevin equation can, for example, be found in <cit.>.
A Néel rotation including anisotropy is then given by $p_1=\tilde{\gamma}\mu_0$, $p_2=\tilde{\gamma}\alpha \mu_0$, $p_3=2\tilde{\gamma}\frac{K^\text{anis}}{M_\mathrm{S}}$, and $p_4=\alpha p_3$ with $\tau=\frac{V_\mathrm{C} M_\mathrm{S}}{2 k_\mathrm{B} T_\mathrm{B} \tilde{\gamma} \alpha}$ and $\tilde{\gamma}=\frac{\gamma}{1+ \alpha^2}$. Here, $V_\mathrm{C}$ is the core volume of the nanoparticles depending on the core diameter $D$, $M_\mathrm{S}$ denotes the saturation magnetization, $\alpha$ is the damping parameter and $k_\mathrm{B}$ and $T_\mathrm{B}$ denote the Boltzman constant and the temperature, respectively. The uniaxial anisotropy constant is denoted by $K^\text{anis}$.

We note that the parabolic partial differential equation in (<ref>) has no dependence on derivatives with respect to the spatial variable $r$.
It can thus be considered as parametric with respect to $r$, respectively the constant offset fields encoded in the selection field.

The equation is solved numerically by using a finite volume method, which is conservative by design in the sense that the integral over the solution does not change over time. This makes sense, since the solution of the Fokker-Planck equation is a probability density. The sphere is segmented into spherical triangles, and considering the probability flow over the edges leads to a system of ordinary differential equations. This system is subsequently solved with a multistep integrator. Here, we exploited the computational toolbox provided in <cit.>, where more details on the implementation can be found.

§.§ Data-driven physically motivated models

Given an appropriate model for the mean magnetic moment $\bar{m}$ it is still unlikely that single scalar parameters such as particle diameter and anisotropy constant, which would represent a perfect monodisperse tracer, will be sufficient for the proper prediction of the measured signal in the MPI setting.
Additionally, in an MPI measurement we do not have direct access to the mean magnetic moment such that we need to take into account uncertainties in the involved quantities describing the system function $s_k$ in (<ref>).
We thus formulate the following physical model-based but still data-driven approaches for the system functions.

The general setting is based on the following assumptions:

(i) The analog filter $\hat{a}$ is not given with necessary accuracy.

(ii) The sensitivity profile $p$ of the receive coil is homogeneous and it is assumed to be given.

(iii) The monodisperse Néel rotation model $\bar{m}_{(D,K^\text{anis},\vartheta)}$ including uniaxial anisotropy constant $K^\text{anis}$, particle diameter $D$, and easy axis $n(\vartheta)=(\cos (\vartheta),\sin(\vartheta),0)^\trans\in S^2$ according to (<ref>), (<ref>), (<ref>), i.e., $\bar{m}_{(D,K^\text{anis},\vartheta)}=\bar{m}$ with $\bar{m}$ from (<ref>), is exploited in the data-based approaches in the following.

We particularly consider easy axes in the $xy$-plane for a given angle $\vartheta$ with respect to the $x$-axis.

With these assumptions, predefined discrete parameter sets $\mathcal{P}_\text{diam},\mathcal{P}_\text{anis}, \mathcal{P}_\text{angle}$, and a frequency index set $I_\text{freq}\subset\{0,1,\hdots,K\}$  we consider the following data-driven models describing the system functions in (<ref>)

M0 Even if the equilibrium model <cit.> is not able to express the dependence on the easy axis orientation $\vartheta$, we include it as a reference to the state of the art in MPI and also for the sake of completeness: For individual $D\in \mathcal{P}_\text{diam}$ and $\vartheta\in \mathcal{P}_\text{angle}$ identify $\hat{a}=(\hat{a}_k)_{k\in I_\text{freq}}$ in

\begin{equation}
s^{(\vartheta)}_k(r)= \hat{a}_k \hat{\phi}_{D;k}(r), \ k\in I_\text{freq}
\end{equation}

for given  $\hat{\phi}_{D;k}(r)=\frac{1}{m_0 T}\int_{0}^{T} p^\trans \frac{\partial}{\partial t}\bar{m}^\text{equ}_{D}(r,t) {\textrm{e}}^{-2\pi i tk/T} \ \d t$ from measurements $S^{\mathrm{calib}}_{\vartheta}$ (measured system matrix for one predefined easy axis angle $\vartheta$; see Section <ref>).
$\bar{m}^\text{equ}_{D}$ is explicitly given by the equilibrium model , i.e,

\begin{equation}
\bar{m}^\text{equ}_{D}(r,t) = m_0 \mathcal{L}_\beta(| H(r,t)|) \frac{H(r,t)}{| H(r,t)|}
\end{equation}

where $\mathcal{L}_\beta: \R \rightarrow \R$ is given in terms of the Langevin function by the following:

\begin{equation}
\label{eq:Langevin}
\mathcal{L}_\beta(z) =  \left( \coth( \beta z) - \frac{1}{ \beta z} \right)
\end{equation}

for $m_0,\beta >0$. Both $m_0=V_\mathrm{C} M_\mathrm{S}$ and $\beta=\frac{\mu_0 V_\mathrm{C} M_\mathrm{S}}{k_\mathrm{B} T_\mathrm{B}}$ with $V_\text{C} = \frac{1}{6}\pi D^3$ depend on the particle core diameter $D$ in a cubic fashion.

M1 For individual $D\in \mathcal{P}_\text{diam}$, $K^\text{anis}\in \mathcal{P}_\text{anis}$, and $\vartheta \in \mathcal{P}_\text{angle}$ identify $\hat{a}=(\hat{a}_k)_{k\in I_\text{freq}}$ in

\begin{equation}
s^{(\vartheta)}_k(r)= \hat{a}_k \hat{\phi}_{(D,K^\text{anis},\vartheta);k}(r), \ k\in I_\text{freq}
\end{equation}

for given

\begin{equation}
\hat{\phi}_{(D,K^\text{anis},\vartheta);k}(r)=\frac{1}{m_0 T}\int_{0}^{T} p^\trans \frac{\partial}{\partial t}\bar{m}_{(D,K^\text{anis},\vartheta)}(r,t) {\textrm{e}}^{-2\pi i tk/T} \ \d t
\end{equation}

from measurements $S^{\mathrm{calib}}_{\vartheta}$ (measured system matrix for one predefined easy axis angle $\vartheta$; see Section <ref>).

M2 In the absence of an accurate estimate of the transfer function and aiming for a polydisperse model, as a preparatory step we generate an improved dictionary of functions such that we follow the following heuristic two step procedure:

* For any $D\in \mathcal{P}_\text{diam}$, $K^\text{anis}\in \mathcal{P}_\text{anis}$, and $\vartheta \in \mathcal{P}_\text{angle}$ obtain an $\hat{a}_{(D,K^\text{anis},\vartheta)}$ according to the procedure outlined for M1 which gives us the preprocessed auxiliary functions

\begin{equation}\label{eq:psi_def}
\psi_{(D,K^\text{anis},\vartheta);k}(r)= \hat{a}_{(D,K^\text{anis},\vartheta);k} \hat{\phi}_{(D,K^\text{anis},\vartheta);k}(r), \ k\in I_\text{freq}.
\end{equation}

* We then obtain a weighted linear combination with respect to the physical parameters to take the polydisperse characteristic of the tracer into account.

Here, we aim for finding weights $w_{D,K^\text{anis}}\geq 0$ such that

\begin{equation}
s^{(\vartheta)}_k(r)=\negthickspace\!\sum_{D\in\mathcal{P}_\text{diam}} \sum_{K^\text{anis}\in\mathcal{P}_\text{anis}}  w_{D,K^\text{anis}} \psi_{(D,K^\text{anis},\vartheta );k}(r), \ k\in I_\text{freq}
\end{equation}

approximates the measurements $S^{\mathrm{calib}}_{\vartheta}$ for any $\vartheta \in \mathcal{P}_\text{angle}$.
In contrast to M1 we obtain global weights not depending on $\vartheta$ but only depending on $D$ and $K^\text{anis}$.

Remark: Note that the fitted $\hat{a}$/$\hat{a}_{(D,K^\text{anis},\vartheta)}$ are parallel to the transfer function in (<ref>) but its absolute also includes some scalar factors (such as $\mu_0$, $m_0$, etc.).

Technical details regarding the parameter identification

We exploit the discrete nature in space given by the discrete points $\{r_\nu\}_{\nu=1,\hdots,N} \subset \Omega$ of the problem setup and introduce the following variables and transformations:

* $S^{\mathrm{model}} = (s_k^{(\vartheta)}(r_\nu))_{k\in I_\text{freq},\nu=1,\hdots,N, \vartheta \in \mathcal{P}_\text{angle}} \in \C^{|I_\text{freq}| \times N \times |\mathcal{P}_\text{angle}|}$ being a tensor including system matrices for the respective model for different orientations $\vartheta$. With $S^\mathrm{model}_\vartheta$/$S^\mathrm{calib}_\vartheta \in \C^{|I_\text{freq}| \times N}$ we denote the respective system matrix for easy axis orientation $\vartheta$.

* $\Phi_{(D,K^\text{anis},\vartheta)}=(\hat{\phi}_{(D,K^\text{anis},\vartheta);k}(r_\nu))_{k\in I_\text{freq},\nu=1,\hdots,N} \in \C^{|I_\text{freq}| \times N}$ being an element of the dictionary of system matrices according to the description of M1. $\Phi_D$ analogously for M0.

* $\Psi_{(D,K^\text{anis},\vartheta)}=({\psi}_{(D,K^\text{anis},\vartheta);k}(r_\nu))_{k\in I_\text{freq},\nu=1,\hdots,N} \in \C^{|I_\text{freq}| \times N}$ being an element of the dictionary of system matrices according to (<ref>)

* $T^\mathrm{mat2vec}: \C^{K\times N} \rightarrow \C^{NK}$ concatenates the columns of the input matrix in one vector.

* $T^\mathrm{diag}: \C^K \rightarrow \C^{K\times K}$ generates a diagonal matrix with the input vector as the diagonal.

Model fitting in any of the previously outlined model approaches is performed by doing a (possibly weighted) least squares fit to the data.

M0/M1 For any $\vartheta \in \mathcal{P}_\text{angle}$ and

\begin{equation*}
q:=\left\{ \begin{array}{ll} D \in \mathcal{P}_\text{diam} & \text{for M0} \\ (D,K^\text{anis},\vartheta) \in \mathcal{P}_\text{diam} \times \mathcal{P}_\text{anis} \times \{\vartheta\} & \text{for M1}  \end{array} \right.
\end{equation*}

we obtain

\begin{align}
\hat{a}_q &= \mathrm{arg\ min}_{\hat{a}\in \C^{|I_\text{freq}|}} \| T^\mathrm{diag}(\hat{a}) \Phi_q - S^{\mathrm{calib}}_{\vartheta} \|_\mathrm{F}^2 \notag \\ &\left(= \mathrm{arg\ min}_{\hat{a}\in \C^{|I_\text{freq}|}} \| T^\mathrm{mat2vec}(T^\mathrm{diag}(\hat{a}) \Phi_q - S^{\mathrm{calib}}_{\vartheta}) \|_2^2\right)
\end{align}

explicitly by (see, e.g.,<cit.>)

\begin{equation}
\hat{a}_{q;k} = \frac{S^{\mathrm{calib}}_{\vartheta;k} \overline{\Phi_{q;k}}^T}{\| \Phi_{q;k} \|_2^2}
\end{equation}

for any $k \in \{1, \hdots,|I_\text{freq}|\}$ where $S^{\mathrm{calib}}_{\vartheta;k}$ and $\Phi_{q;k}$ denote the respective $k$-th row.

The first step is given in terms of  $\Psi_q=T^\mathrm{diag}(\hat{a}_q) \Phi_q,\ q \in  \mathcal{P}_\text{diam} \times \mathcal{P}_\text{anis} \times \mathcal{P}_\text{angle}$ where $\hat{a}_q$ are obtained as outlined in M1.
In the second step, we utilize for any $\vartheta \in \mathcal{P}_\text{angle}$ and $q^{(\vartheta)}_j:=(D_j,K_{\mathrm{anis};j},\vartheta) \in \mathcal{P}_\text{diam} \times \mathcal{P}_\text{anis} \times \{ \vartheta \}$, $j=1,\hdots,J:=|\mathcal{P}_\text{diam}|\cdot |\mathcal{P}_\text{anis}|$, (such that $\{q_j^{(\vartheta)}\}_{j=1,\hdots,J}=\mathcal{P}_\text{diam} \times \mathcal{P}_\text{anis} \times \{\vartheta\}$), we obtain $w\in \R_+^J$ by solving

\begin{equation}\label{eq:LSfitM2}
\min_{w\in \R^J_+} \sum_{\vartheta \in \mathcal{P}_\text{angle}} \| M_{\vartheta}w - T^\mathrm{mat2vec}(\Sigma_{\vartheta} S^{\mathrm{calib}}_{\vartheta}) \|_2^2
\end{equation}

with dictionary matrices
$$M_{\vartheta}:=\left[T^\mathrm{mat2vec}(\Sigma_{\vartheta}\Psi_{q_1^{(\vartheta)}}), \hdots, T^\mathrm{mat2vec}(\Sigma_{\vartheta} \Psi_{q_J^{(\vartheta)}}) \right]$$

for $\vartheta \in \mathcal{P}_\text{angle}$ which includes the simulations at different positions in each column. Each column corresponds to one physical parameter tuple $q_j^{(\vartheta)}$. $\Sigma_{\vartheta} \in \R^{|I_\text{freq}| \times |I_\text{freq}|}$ is a diagonal weighting matrix with respect to the individual frequency components allowing a weighted least squares fit.

Throughout this work, the weighting was chosen such that for each frequency, the measurement is normalized with respect to the infinity-norm over the offset fields:

$$\Sigma_{\vartheta} = T^\mathrm{diag}\left(\left( \| S^{\mathrm{calib}}_{\vartheta;k}\|_\infty^{-1}\right)_{k\in I_\text{freq}}\right).$$

For numerical treatment, the least-squares functional is expanded, taking into account that the data is complex while the weights are real. Thus, we obtain a so-called quadratic program with linear constraints of the form $\min_w \frac{1}{2}w^T A w + f^Tw$ for a positive semidefinite real matrix $A$ and a real vector $f$, with the constraints $w_i \geq 0$ for $i=1,\dots,J$. We solve this problem by using a standard interior-point algorithm, as described in <cit.>, among many other works.

However, since a rather high resolution of parameters was considered for the data fitting in model M2, many matrix columns of the resulting matrices $M_{\vartheta}$ differ only slightly. This leads to an ill-conditioned problem and as a result, the weight vector $w$ is not interpretable, since it might simply reflect noise in the measurements. This makes it necessary to employ regularization techniques in order to obtain weights that better reflect the real-world parameter distribution. We chose to employ the truncated singular value decomposition (TSVD) or reduced-rank regularization approach which is motivated by the equivalent formulation of (<ref>) as

\begin{equation}\label{eq:LSfitM2_alternative}
\min_{w\in \R^J_+}  \| Mw - \underset{=:y}{\underbrace{(T^\mathrm{mat2vec}(\Sigma_{\vartheta} S^{\mathrm{calib}}_{\vartheta}))_{\vartheta\in \mathcal{P}_\text{angle}}}} \|_2^2,
\end{equation}

where the matrix $M\in \C^{|I_\text{freq}|\cdot |\mathcal{P}_\text{angle}| \times |\mathcal{P}_\text{diam}| \cdot |\mathcal{P}_\text{anis}|}$ consists of the stacked $M_{\vartheta}, \vartheta \in \mathcal{P}_\text{angle}$. For this we compute the SVD, i.e., $M = U \Sigma V^\ast$ for unitary matrices $U, V$ and a real diagonal matrix $\Sigma$. Choosing a low-rank approximation by the SVD is particularly motivated by the following relationship

\begin{align}
\| M w -y \|_2 &= \| U^\ast M w - U^\ast y  \|_2 \notag \\
%    & \leq \|(U^\ast M - U_r^\ast M) w \|_2 + \| U_r^\ast M w - U^\ast y  \|_2 \notag \\
&\approx \| U_s^\ast M w - U_s^\ast y \|_2, \label{eq:approx_low_rank}
\end{align}

where $U_s := (U_{ij})_{i=1,\hdots,|\mathcal{P}_\text{angle}|\cdot|I_\text{freq}|; j=1,\hdots,s}$ includes the first $s$ columns of $U$, corresponding to the $s$ largest singular values, i.e., $s$ being the regularization parameter.

We thus obtain the regularized minimization problem

\begin{align*}
\min_{w \in \mathbb{R}^J_+} \| U_s^\ast M w - U_s^\ast y\|_2^2.
\end{align*}

Note that the global error measure (see Section <ref>) $\varepsilon_\mathrm{global}(S^\mathrm{model},S^\mathrm{calib}) \sim \| M w -y \|_2$. Due to the fast decay of singular values, using a low-rank approximation thus has only little influence on the global error measure but stabilizes the parameter reconstruction.

§.§ Distance Measure

We exploit the distance measure as introduced in <cit.>.
For quantifying the accuracy of a model
we express the frequency index $k$ in terms of the mixing orders $m_x$, $m_y$ <cit.> fulfilling

\begin{align} \label{Eq:mix}
k(m_x,m_y) &= m_x N^{\text{dens}} + m_y (N^{\text{dens}}+1)
\end{align}

Here, $N^{\text{dens}}$  is the parameter that controls the density of the Lissajous sampling pattern and for each $k$ the factors $m_x$ and $m_y$ with the smallest absolute value fulfilling (<ref>) are considered.
Now with $S^{\text{calib}}_{\vartheta;k(m_x,m_y)}$ being the measured and $S^{\text{model}}_{\vartheta;k(m_x,m_y)}$ being the modeled $k(m_x,m_y)$-th system matrix row (also named frequency component) for a given angle $\vartheta$ we consider the normalized root mean square error (NRMSE)

\begin{align} \label{Eq:error1}
\varepsilon(S^{\text{calib}}_{\vartheta;k}, S^{\text{model}}_{\vartheta;k}) &= \frac{\frac{1}{\sqrt{N}}\Vert  S^{\text{calib}}_{\vartheta;k}- S^{\text{model}}_{\vartheta;k} \Vert_2}{\Vert  S^{\text{calib}}_{\vartheta;k} \Vert_\infty},
\end{align}

where we decided to use the infinity norm for normalization.
In addition to this per-row metric we consider global measures over $k \in I_\text{freq} =\{k(m_x,m_y)\,|\, m_x=0,\dots,M_x$, $m_y=0,\dots,M_y \}\setminus \{ (0,0),(0,1),(1,0)\}$ where $M_x,M_y\in\N$ are upper bounds for the mixing order and $\vartheta \in \mathcal{P}_\text{angle}$ with $\mathcal{P}_\text{angle}$ being a discrete set of orientations.
In particular we consider the

* per frequency NRMSE measure

\begin{align}
&\varepsilon^\text{freq}_{k}(S^{\text{calib}}, S^{\text{model}}) \nonumber \\ &= \frac{1}{\sqrt{|\mathcal{P}_\text{angle}|}} \left( \sum_{\vartheta \in \mathcal{P}_\text{angle}} \varepsilon(S^{\text{calib}}_{\vartheta;k}, S^{\text{model}}_{\vartheta;k})^2 \right)^{\frac12}
\end{align}

* and the global NRMSE measure

\begin{align}
&\varepsilon^\text{global}(S^{\text{calib}}, S^{\text{model}}) \notag \\ &= \frac{1}{\sqrt{|I_\text{freq}|}} \left( \sum_{k\in I_\text{freq}} \varepsilon^\text{freq}_{k}(S^{\text{calib}}, S^{\text{model}})^2 \right)^{\frac12}  \label{eq:globalMeasure}
\end{align}

We note that in both measures the quadratic sums have been additionally performed along the receive channel dimension, which we omitted in this work to simplify the notation.

§.§ Experimental Setup

Experimental data is acquired using a pre-clinical MPI system (Bruker, Ettlingen, Germany) that is equipped with three send and three receive channels. The scanner is operated with a 2D excitation in the $xy$-plane (horizontal) with 12.59 excitation field strength in $x$-direction and  13.55 in $y$-direction. These two values have been retrospectively determined using a calibration coil and thus do not exactly match the ones used at the scanner console. The drive-field frequencies are $f_x=2.5/102$ MHz and $f_y=2.5/96$ MHz leading to a repetition time of  652.8 and a density of the Lissajous trajectory of $N^{\text{dens}} = 16$.

For signal reception we use a 3D pick-up coil while taking only the two channels in $x$- and $y$-direction into account. The analog signal is band-pass filtered and digitized at a sampling rate of 2.5.

The system matrices are sampled on a grid of size $11 \times 11$ at $N=121$ equidistant measurement points with different offset values ranging from
-14.35 to 14.35 in $x$-direction and -14.71 to 14.71 in $y$-direction. Both intervals have been calibrated using a  3-channel gaussmeter with a three-axis Hall probe (model 460, Lake Shore, Westerville, USA). Somewhat uncommonly we report the sampled grid using offset values since we measure the matrix with a static sample at the center and move the sampling trajectory using a 3D offset field that can be adjusted using the scanner. This method is called hybrid system matrix <cit.> and it is advantageous for the present study as spatial field imperfections can be ignored since the sample is always measured at the same position. Since the offset field in a real MPI imaging experiment would be translated to space, we consider the $x$- and $y$-offset field dimensions to be spatial dimensions with a virtual and perfect selection field applied. The overscan <cit.> of the calibration data is 13.9 % in the $x$-direction and 8.6 % in the $y$-direction. All measurements have been done with an active gradient of -1.0;-1.0;2.0 in $x$-, $y$-, and $z$-direction, which, however, is only enabled to suppress background signals that are much stronger with disabled selection field. This should not be confused with the virtually applied selection field gradient, which can have any value.
The calibration measurements are block averaged with a factor of 5000 amounting to a measurement time of about 3.3 per position.

For system matrix calibrations a single cylindrical delta sample of immobilized magnetic nanoparticles with magnetic easy axis in parallel to a specified direction is designed. As shown in Fig. <ref>, the delta sample is assembled from two parts, which are 3D printed using a Form 3 and clear resin (Formlabs Inc., Somerville, USA). The top part (carrier) has a 3 diameter and 3 height cavity to hold magnetic nanoparticles. The bottom part (connector) is interchangeable and can be assembled with the carrier. The assembled delta sample can be mounted onto the robot arm of our MPI system such that it has a defined orientation inside the $xy$-plane of the MPI system, which depends on the connector chosen. Initial preparation of the delta sample is done by assembling the carrier with a 0-connector. The carrier is then filled with a mixture of 20 fluidal perimag (Micromod Partikeltechnologie GmbH, Rostock, Germany) with an iron concentration of 89Fe and about 21 sodium alginate powder. At this point the nanoparticles can still change their spatial orientation by rotation. The sample is then immediately put in-between two cylindrical neodymium magnets with a 0.4 axial field, which was measured using the same gaussmeter as used for the MPI scanner offset calibration. The magnetic easy axes of the nanoparticles align in parallel to the external field <cit.> and slowly lose their mobility due to drying of the mixture, which fixes the orientation of their easy axes permanently. Placed into our MPI system the delta sample with the 0-connector has a $\vartheta = \SI{0}{\degree}$ angle between the parallel aligned magnetic easy axis and the $y$-axis of our MPI system.

In this preparation process the mixing of perimag and sodium alginate powder is most error prone, which is why only a single carrier is prepared and different orientations inside our MPI system are realized by disassembling carrier and connector and reassembling the former with a different suitable connector. For our measurements a series of connectors was printed, allowing us to measure 12 system matrices at angles

\begin{equation}
\mathcal{P}_\text{angle} := \left\lbrace \SI{0}{\degree}, \SI{15}{\degree},\dots,\SI{165}{\degree}\right\rbrace
%\left\lbrace l \cdot 15^\circ \;|\; l=0,\dots, 11 \right\rbrace
\end{equation}

resulting in the measured system matrix tensor $S^\mathrm{calib} \in \C^{K\times N \times |\mathcal{P}_\text{angle}|}$.

A schematic overview of delta sample (top) and measurement (bottom) preparation is shown. The delta sample is assembled from two parts (top left). Its cylindrical depression is filled with a brownish alginate-tracer mixture (top center). Immediately, the sample is placed in a magnetic field (red sided permanent magnets) to align the magnetic nanoparticles magnetic easy axis in parallel and immobilize particles by drying (top right). The sample has multiple connectors, which can be selected (bottom left). Each connector can be assembled with the top part of the sample (bottom center) and attached to the robot arm of our MPI system (bottom right), providing a connector specific easy axis orientation inside the MPI system.

§.§.§ Simulation Setting

For the parameter optimization we selected for each model parameter a finite set of values and performed simulations with them. For the particle diameter $D$ we chose

\begin{equation}
\mathcal{P}_\text{diam} = \left\lbrace \SI{14}{nm}, \SI{16}{nm}, \dots, \SI{28}{nm} \right\rbrace,
\end{equation}

i.e., we chose an increment of $\SI{2}{nm}$. For the particle anisotropy constant $K^\text{anis}$ we chose

\begin{equation}
\mathcal{P}_\text{anis} = \left\lbrace \SI{100}{Jm^{-3}}, \SI{200}{Jm^{-3}}, \dots, \SI{9000}{Jm^{-3}} \right\rbrace,
\end{equation}

i.e., an increment of $\SI{100}{Jm^{-3}}$. Finally for the angles $\mathcal{P}_\text{angle}$ we took the same angles as used for the measurements.

§ RESULTS

§.§ Parameter Optimization

Parameter optimization for the model M1. The bottom shows the global NRMSE $\varepsilon^\text{global}$ combined for $x$- and $y$-channel as a function of the particle diameter $D$ and the particle anisotropy $K^\text{anis}$. The minimum of the discrete error map is reached at ($D=\SI{22}{nm}$, $K^\text{anis}=\SI{1100}{Jm^{-3}}$). For selected parameters indicated by white circles, the top part of the figure shows one row ($m_x=3$, $m_y=3$) of the modeled system matrices for the $x$-channel and the angle $\vartheta = 30^\circ$ in comparison with the measured system matrix shown in the last column. The yellow arrow indicates the easy axis for the considered angle. To highlight the influence of a change of the particle anisotropy, the central part of the figure shows a profile through the merging wave hills of the system matrix pattern. The profile is indicated as a white arrow in the considered system matrix pattern shown in the top part of the figure.

The discrete parameter optimization for the equilibrium model M0 resulted in the diameter $D=\SI{22}{nm}$ as the optimal value minimizing (<ref>). For the more interesting case of M1, the optimization found the pair $(D,K^\text{anis})=(22\,\text{nm}, 1100\, \text{Jm}^{-3} )$ to be the global minimizer of (<ref>).

More details on the influence of the particle diameter and the anisotropy constant on the resulting system matrix are illustrated in Fig. <ref> for a certain frequency component ($m_x=3$, $m_y=3$). In the top part of the figure the simulated system matrices and on the lower part the values of the error metric (<ref>) are shown. First of all one can see that the error functional in the bottom part is a rather smooth curve without oscillations. Here, one can also see that the dimensions $D$ and $K^\text{anis}$ are not uncorrelated, which means that a larger particle diameter requires a smaller anisotropy constant to result in similar error values. Independent of these qualitative findings there is a minimum of the global error at $0.2049$. Selected system matrices shown above reveal the influence of the parameters $D$ and $K^\text{anis}$. Decreasing $D$ for fixed $K^\text{anis}$ leads to an increase of a spatial blurring with respect to the offset field dimension, whereas this blurring has an isotropic nature.
As the diameter influences the diffusion coefficient $\frac{1}{2\tau}$ in (<ref>) only, when decreasing the diameter (resulting in larger $\frac{1}{2\tau}$), a larger offset field is necessary to obtain a similar effect. The blurring thus also includes a zooming effect to some degree. But due the delicate interplay of drive field, anisotropy, and offset field in the field vector $b$ in (<ref>), zooming is not the only effect which contribute to the blurring such that we stick to the term isotropic blurring as it is caused by the isotropic diffusion coefficient $\frac{1}{2\tau}$ .

In contrast,  increasing $K^\text{anis}$ for fixed $D$ leads to an anisotropic blurring, which is oriented perpendicularly to the easy axis indicated as a yellow line.
As $\tau$ is independent of $K^\text{anis}$, this is mainly caused by the structure of $b$ (see (<ref>)) and it becomes apparent for large $K^\text{anis}$.
In the case that the offset field is parallel to the easy axis, the easy axis component and the offset field component of $b$ align such that we can observe the pattern along the easy axis in the field of view. The pattern in this direction is caused by the change of the offset field strength as it contributes effectively to the dominating anisotropy effect.
If the offset field is perpendicular to the easy axis, an effective contribution by the offset field to the behavior of $b$ is small such that a larger offset field strength would be required to change the magnetization behavior significantly. This anisotropic blurring effect can be seen in the top part of Fig. <ref>, particularly in the right column.
When comparing the simulated data with the measurement shown on the upper right, one can see that $(22\,\text{nm}, 1100 \,\text{Jm}^{-3} )$ yields the most similar looking pattern, which, as we discussed before, is also the value chosen based on the error metric.

Considering M2, the TSVD-regularized (see Section <ref>) optimization was performed for different regularization parameters $s$. Here, the regularization parameter takes the form of the number of singular values that were taken into account, i.e., the rank of the low-rank approximation of $M$. Lower numbers correspond to more aggressive rank reduction, and thus to stronger regularization. The values of the NRMSE are displayed in Table <ref>. We observe that the rank can be reduced quite aggressively before the results start to deteriorate. Only if less than 3% of the largest singular values are considered, the error starts growing significantly. For that reason, we chose this value as the best balance between accuracy and interpretability. In other words, the 720 parameter combinations considered can be reduced to 21 effective structures by removing redundant information, resulting in a model more robust to measurement inaccuracies and noise.

The obtained weight vector $w$ for M2 is illustrated in Fig. <ref>.
Here we can see  similarities to the optimal parameter of model M1 as this $(D,K^\text{anis})$ tuple got a large weight value. But we can also observe smaller contributions by surrounding particle diameters and different anisotropy values. For instance, the two smallest simulated anisotropy constants $K^\text{anis} = \text{\SIlist{100;200}{Jm^{-3}}}$ are assigned nonzero weights for $D = \SI{24}{nm}$. The weights for different anisotropy constants may be interpreted in part as a hint towards the shape distribution of the particles, since shape anisotropy plays a large role in the overall anisotropic behavior <cit.>. In this context, small anisotropies hint towards more spherical particles while larger anisotropies may stand for ellipsoid shapes. There are also small contributions by the largest simulated diameter for various anisotropy constants. This may reflect larger particles or even clusters. In general, the role of clusters and chains in this consideration is unclear. We expect that the polydisperse nature as well as agglomeration effects lead to the observed structures. Consequently, the weights should serve as a hint at best for the true distribution of nanoparticle properties.

% of SVs (Matrix rank $s$)      NRMSE $\varepsilon^\text{global}$
100 (720)     0.1994

5 (36)     0.1995

4  (28)     0.1995

3 (21)     0.1996

2  (14)     0.2005

1  (7)     0.2034

0.5 (3)     0.2178

Global NRMSE $\varepsilon^\text{global}$ for the truncated SVD approach for different numbers of singular values considered and the corresponding matrix ranks.

Calculated weights $w_{D,K^\text{anis}}$ for model M2 determined by using the dictionary summarized in (<ref>).

§.§ Qualitative Comparison

Comparison of the modeled MPI system matrices for models M0 to M2 with the measured system matrix for the $x$-channel. The left part of the figure shows a comparison for a fixed frequency component ($m_x=3$, $m_y=3$) and various angles in $30^\circ$ steps. The right part of the figure shows the data for a fixed angle of $60^\circ$ and several frequency components at mixing orders $(m_x, m_y) \in \{(1,1),(1,3),(3,3),(5,3),(5,5)\}$. The orange dashed box indicates the data that is present in both parts of the plot. In addition to the matrix patterns itself, the complex-valued difference $D_{k,\vartheta} = S^\text{model}_{k,\vartheta} - S^\text{calib}_{k,\vartheta}$ is shown. All data uses the complex colormap shown in Fig. <ref> whereas the same windowing is used for each frequency index $k$.

We proceed with a qualitative comparison of the optimized models, M0 to M2. Fig. <ref> shows parts of the modeled and measured system matrices. Since the matrices are too large such that one could show them at once, we show selected parts sampled along different axes. On the left part of the figure the frequency component for $m_x=3$, $m_y=3$ is again selected but this time every second measured angle is shown. Here one can see that model M0 is not capable of modeling the effect of the anisotropic and angle dependent blurring and just provides a mean model with no merge of wave hills and no intensity variations for the different angles. In contrast, M1 much better resembles the measurement for the chosen frequency component and one can clearly see both features: An angle dependent anisotropic blurring and an angle dependent change in intensity. This can also be seen in the different frequency components for a fixed angle, which are shown on the right hand side of the figure. To highlight the differences between the measurement and the model, Fig. <ref> also shows difference maps that highlight the angle dependent deviation for M0 and the much lower overall deviation for the models M1 and M2. When switching to the model M2, differences to model M1 are less pronounced. For $(m_x=1,m_y=3)$ one can see that the first pixel row has a sightly higher value for M2, which better resembles the intensity for the measured frequency component.

§.§ Quantitative Comparison

Differences between the measured and modeled system matrices for models M0–M2. On the top the error $\varepsilon^\text{freq}_{k}$ is illustrated for mixing orders $m_x = 0,\dots, 9$ and  $m_y = 0,\dots, 9$, while the individual errors along the angle and channel dimensions have been combined. The values for $(m_x,m_y) \in \{ (0,0), (0,1),(1,0) \}$ have been omitted since they do not contain any particle signal in the measurements (transfer function close to zero). The lower plot shows a statistical analysis of the error $\varepsilon$ as a box plot. Here, the statistical analysis has been performed along the frequency and channel dimensions, while for each angle an individual analysis has been performed.

After the qualitative comparison of the modeled MPI system matrices with the measured one we next carry out a quantitative comparison. Here, we consider the distance measures introduced in Section <ref>. We first focus on the frequency dependent measure $\varepsilon_k^\text{freq}$ where the data for all angles and receive channels has been combined into a single number per frequency index $k$, respectively mixing order pair $(m_x,m_y)$. The resulting errors are shown in the upper part of Fig. <ref>. The first finding is that the models M1 and M2 result in considerably lower differences to the measurements than the equilibrium model, which confirms the results of the qualitative comparison. When comparing the errors for M1 and M2, one can hardly see a difference, which underlines the observation that we made in the qualitative comparison. To highlight that there are smaller differences we plot the difference $\left\vert\varepsilon^\text{freq,M1}_k- \varepsilon^\text{freq,M2}_k\right\vert$ in the upper right part of the figure. One can see that model M2 shows a slight improvement that is two orders of magnitude lower than the actual error. The improvements is highest when one of both mixing orders is 1.

Next we investigate the deviation in dependence of the considered easy axis angle. Here, we perform a statistical analysis on the differences $\varepsilon$ and set up a box-whisker plot for each angle and each model in the lower part of Fig. <ref>. First of all, one can see again that models M1 and M2 outperform the equilibrium model M0 and show much lower error. Interestingly, one can also see an angular dependency of the median error, which varies between $0.141$ and $0.255$ for model M1 and between $0.234$ and $0.305$ for model M0. Finally, the global NRMSE for the three considered models is summarized in Table <ref>, which wraps up our quantitative comparison.

M0     M1     M2
NRMSE $\varepsilon^\text{global}$      0.2976     0.2049     0.1996

Global NRMSE $\varepsilon^\text{global}$ between the modeled system matrices for M0 to M2 and the measured system matrices.

§ DISCUSSION AND OUTLOOK

Simulating the particle magnetization in magnetic particle imaging remains a very challenging task, which is due to the complex interplay of the particles on a micromagnetic level and a large possible parameter space involving size and anisotropy which is not exactly known for a set of more than billions of particles that is typically located in an image voxel and contributes a signal during a typical MPI experiment. One thus needs to consider ensembles which then only can describe the physical phenomenon in a mean sense.

In this work we progressed in the development of tools for describing MPI experiments based on physical models. Our approach uses in its core a magnetization model exploiting the Fokker-Plank equation, making the assumption that particle-particle interaction as well as agglomeration effects can be neglected. The physical model is then enriched by some data-driven components to bring the model in alignment with experimental data. First of all we still use a fitted transfer function, which simplifies model-based approaches tremendously since it allows to account for amplitude and phase differences that can be caused by an absent / inaccurate transfer function of the receive chain or by imperfectly estimated particle and scanner parameters. Secondly, we use a data-driven approach for obtaining the particle parameters (diameter and anisotropy constant). Here we used two different approaches. One simply pics the parameters that best describe the experiment while the other allows for linear combination of simulations obtained for different parameter sets. This dictionary approach is very flexible but it also makes physically non-plausible solutions possible if no physical priors are integrated, which we did not in this study. Our results indicate that the dictionary approach only slightly improves the result over the single parameter set simulation in our setting.
A natural extension of M2 is a joint reconstruction of one transfer function for each receive channel only and the weight vector simultaneously. On the one hand this models the physical relationship more properly, but on the other hand its solution requires usage of more sophisticated nonlinear programming tools. This direction of research remains immediate future work.

In our experiment we used an immobilized particle setup with oriented easy axis, which marks both an interesting and challenging setup for model studies since it allows to change the orientation as a parameter and to compare if the model can describe the changes in the measured signal. As expected, the equilibrium model can not describe changes in easy axis orientation, whereas both our models M1 and M2 were successful in describing the observed phenomenons which manifest themselves as an anisotropic blurring of the system matrix patterns in perpendicular direction to the magnetic easy axis.
Of course, this might also be an interplay between easy axis orientation and the Lissajous pattern.
Future works should include the study of different excitation patterns and might also include excitation optimization for multi-contrast reconstruction with respect to easy axis orientations.

We note that there was some angular dependency in the model accuracy, which cannot be explained by the model itself. One reason might be that the positioning of the sample or the stability of the MPI scanner did slightly change leading to a higher deviation for, in particular, the $150^\circ$ case.

Our focus in this work was not on the deployment of a simulation-based model, which would require additional work. First of all we used a fully measured system matrix for parameter estimation (transfer function, particle diameter and particle anisotropy), which would not be available in practice since the modeled system matrix aims to replace the measurement. One route that has already been exploited for the equilibrium model in <cit.> is to rely on a low-resolution based system matrix, which has the advantage that its calibration scans can be done in few minutes and one still has reference data to compare the model with. This can be combined with compressed sensing based system matrix recovery <cit.> or simultaneous high resolution system matrix and concentration reconstruction <cit.>. In particular, motivated by the structure of the reconstructed weight vector in this work, a sparsity-promoting approach with respect to the dictionary obtained by the physically-driven model we developed can be highly beneficial to reduce the amount of necessary calibration data. We also did not deploy the modeled system matrices for image reconstruction yet since this is its own topic and needs more sophisticated phantoms in case of immobilized samples with aligned easy axis. Subject of a future study will be the development of an immobilized phantom with differently oriented samples in the FOV such that the evaluation can be done based on the multi-contrast reconstruction outlined in <cit.>.

We note that the considered model-system with axis-aligned easy axis is not just of theoretical interest but has clear application in medical imaging scenarios in particular in the field of interventional imaging where the instruments will be coated with immobilized particles that can be additionally orientated.
This would enable the determination of position and orientation of a medical instrument such as a catheter, guide wire, or stent. As discussed in <cit.>, a large number of calibration measurements is necessary for image reconstruction in this scenario, which will require model-based approaches for system matrix calibration since the calibration time on a fine spatial grid would be prohibitive in case of 3D excitation, which would be employed in a clinical environment. Our work contributes to this goal and provides a platform for future research where the goal will be to remove the remaining minor but still observable model-measurement mismatch that our model still has.

§ ACKNOWLEDGEMENTS

H. Albers and T. Kluth acknowledge funding by the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft) - project 426078691.


# Resolving a discrepancy between experimental and theoretical lifetimes in
atomic negative ions

Tomas Brage<EMAIL_ADDRESS>Division of Mathematical Physics,
Department of Physics, Lund University, Sweden  Jon Grumer
<EMAIL_ADDRESS>Division of Mathematical Physics, Department of
Physics, Lund University, Sweden

###### Abstract

Recently the lifetime of the excited ${}^{2}P_{1/2}$-state of S- was measured
to be $503\pm 54$ s (Bäckström et al. Phys. Rev. Lett. 114, 143003 (2015)).
The earlier theoretical lifetime of $436$ s was clearly outside the
experimental error bars. To investigate this discrepancy we have performed
systematic and large-scale multiconfiguration Dirac-Hartree-Fock calculations
for this system. After including a careful treatment of correlation and
relativistic effects, we predict a well-converged value of $492$ s for this
lifetime, with an uncertainty considerably less than 1%, thereby removing the
apparent conflict between theory and experiment. We also show that this result
corresponds to the non-relativistic limit in the $LS$-approximation for the M1
transition within this ${}^{2}P$ term. We also demonstrate the usefulness of
the latter approach for ${}^{2}P$ transitions in O-, Se- and Te-, as well as
for analogous M1 transitions within ${}^{2}D$ terms in Ni- and Pt- ions.

###### pacs:

31.15.A-, 31.15.-p, 32.70.Cs

##  I Introduction

Negative atomic ions remain an interesting and challenging problem to both
experimental and theoretical studies Bäckström  _et al._ (2015); Andersen
(2004); Pegg (2004); Andersen  _et al._ (1999). In the theoretical models we
are faced with strong correlation contributions, since these system are far
from dominated by a central, nuclear potential, but are rather defined by
electron-electron interactions Fischer (1989). On the experimental side, the
spectroscopic identification of these systems is extremely challenging since
they only have a few bound states - often only the ground state. In cases
where more than one state is bound, we can observe the bound-bound spectrum of
these ions. However, we are then restricted to single- and few-line
spectroscopy which, except in very few cases Bäckström  _et al._ (2015);
Walter  _et al._ (2014), only involve states of the same parity. To support
the experimental determination, accurate theoretical predictions of transition
properties are therefore needed.

Halogen-like negative ions have recently been investigated Bäckström  _et al._
(2015); Blondel  _et al._ (2006); Andersson  _et al._ (2006). These have
ground states of the form $np^{5}$ ${}^{2}P_{3/2}$, and the $np^{5}$
${}^{2}P_{1/2}$ as the only excited states. This opens up a bound-bound decay
channel ${}^{2}P_{1/2}$ $\rightarrow$ ${}^{2}P_{3/2}$, which is dominated by a
forbidden, magnetic dipole (M1) transition in all the cases investigated here.

A lingering problem Blondel  _et al._ (2006); Andersson  _et al._ (2006) in
these studies seemed to be a disagreement between experimental Bäckström  _et
al._ (2015) and theoretical lifetimes from multiconfiguration Hartree-Fock
calculations Fischer  _et al._ (2007), which casts some doubt on the
identification of this transition in S- (where $n=3$). In this work we
therefore present a systematic and elaborate theoretical investigation of the
${}^{2}P$ transition in S- based on the relativistic multiconfiguration Dirac-
Hartree-Fock (MCDHF) method, and compare with other available experimental and
theoretical results. In addition we also present results for this transition
in some other negative ions (O-, Se- and Te-) as well as for the analogous M1
transition in ${}^{2}D$ terms of Ni- and Pt- ions.

##  II Theoretical Method

In the MCDHF method the atomic eigenstate is represented by an atomic state
function (ASF),

$\Psi(\Gamma;J)=\sum_{i=1}^{N}c_{i}\Phi(\gamma_{i}J)\;,$ (1)

where $\Phi(\gamma_{i}J)$ are configuration state functions (CSFs),
constructed as angular-momentum-coupled, antisymmetrized products of one-
electron Dirac-orbitals. The sum over CSFs in (1) is determined by
substitutions of one-electron orbitals from one or more reference states, in
this case $1s^{2}\ldots\,3p^{5}$ ${}^{2}P_{1/2,3/2}$, to an active set of
correlation orbitals. This is increased systematically while monitoring the
convergence in order to build effectively complete wavefunctions and obtain
accurate results for the physical quantities of interest - which in this case
are fine-structure separations and transition rates.

The construction of the correlation models also relies on the definition of
two parts of the electron cloud of the ions - the core and the valence
electrons. In the present approach the core is defined as all subshells with
$n<3$, while the valence subshells are defined as those with $n=3$. Depending
on the physical quantity under investigation, the most important correlation
is often between electrons in the valence subshells, while core-valence and
core-core correlation generally is of less importance in atoms and positive
ions (note that core-core effects play an important role when calculating e.g.
isotope shifts Filippin  _et al._ (2016)). In spite of the fact that
correlation, including core-valence and also core-core, is more important in
negative ions, the partition of the electron subshells into two parts is still
a useful concept. The correlation models employed in the present work are
defined in the following section.

Details on the MCDHF method and its non-relativistic counterpart can be found
in the recently published review Fischer  _et al._ (2016) and book Froese
Fischer  _et al._ (1997) by Fisher et al. as well as in the comprehensive book
by Grant Grant (2007). In the present work we use the MCDHF method as
implemented in the most recent version of the Grasp2k code by Jönsson et al.
Fischer  _et al._ (2007).

##  III Results and Discussions

Table 1:  Energies ($E$) and rates ($A$) for the ${}^{2}P_{1/2}\rightarrow^{2}P_{3/2}$, transition in S-, together with lifetime ($\tau$) of the ${}^{2}P_{1/2}$ for different correlation models with systematically enlarged active sets of orbitals. The $\delta E$ values show the relative difference between our ab initio energies and the experimental value $E_{exp}$ reported by Blondel et al. Blondel  _et al._ (2006). This value for the energy and values for the lifetime from other calculations and experiments are given in the bottom of the table. The various model acronyms are defined in the text. The _adjusted_ values for $A$ and $\tau$ are rescaled to the experimental energy of Blondel et al. | ab initio |  | adjusted
---|---|---|---
| $E$ | $\delta E$ | $A$ |  | $A$ | $\tau$
Model | (meV) | (%) | (ms-1) |  | (ms-1) | (s)
DF | 63.3024 | 5.59 | 2.3930 |  | 2.0327 | 491.97
DFBQ | 60.6233 | 1.12 | 2.1019 |  | 2.0327 | 491.96
VV3SD3 | 56.9078 | -5.08 | 1.7386 |  | 2.0327 | 491.96
VV3SD4 | 59.0948 | -1.43 | 1.9469 |  | 2.0327 | 491.96
VV3SD5 | 59.3941 | -0.95 | 1.9756 |  | 2.0327 | 491.95
VV3SD6 | 59.5337 | -0.70 | 1.9905 |  | 2.0327 | 491.97
VV3SD7 | 59.5715 | -0.63 | 1.9943 |  | 2.0327 | 491.96
VV3SD8 | 59.5994 | -0.59 | 1.9971 |  | 2.0327 | 491.96
CV2SD4 | 61.8510 | 3.17 | 2.2322 |  | 2.0327 | 491.96
CV2SD5 | 60.6721 | 1.20 | 2.1070 |  | 2.0327 | 491.96
CV2SD6 | 60.1256 | 0.29 | 2.0505 |  | 2.0327 | 491.96
CV2SD7 | 60.4645 | 0.86 | 2.0854 |  | 2.0327 | 491.96
VV3SDT3 | 56.6467 | -5.81 | 1.6985 |  | 2.0327 | 491.95
VV3SDT4 | 58.6979 | -2.09 | 1.9079 |  | 2.0327 | 491.96
VV3SDT5 | 58.9217 | -1.72 | 1.9298 |  | 2.0327 | 491.96
$E_{exp}$ |  $59.9507\pm 0.0004$ Blondel  _et al._ (2006)
$\tau_{exp}$ |  |  $503\pm 54$ Bäckström  _et al._ (2015)
$\tau_{analytical}$ |  | $491.93^{a}$
$\tau_{others}$ |  |  $437$ Andersson  _et al._ (2006)
aEq. (6) with $\sigma_{exp}$ from Blondel  _et al._ (2006).

The result from various correlation models are presented in Tab. 1 and
compared with other experimental and theoretical results. The first column
defines the theoretical method and correlation model where the DF acronym
refers to Dirac-Fock (i.e. when a single CSF defines the ASF in Eq. (1)) while
DFBQ to a model which includes also Breit and leading radiative QED
corrections (vacuum polarization and self-energy) Jönsson  _et al._ (2013).
The following models represent systematically enlarged correlation models. VV3
denotes models with valence-valence correlation within $n=3$ subshells while
CV2 adds additional core-valence correlation contributions with $n=2$
subshells. SD implies single and double substitutions from the reference set
(the $np^{5}\,{}^{2}P_{J}$ CSFs where $n=3$ for S-) and SDT refers to single,
double and triple substitutions. The subscript number in the end of the model
acronym designation denotes the maximum principle quantum number, $n_{max}$,
of the orbitals in the active set. The models represent valence-valence and
core-valence correlation, and includes single, double and triple substitutions
to an active set of orbitals.

We present two results for the rates - first the ab initio, $A_{ab-initio}$,
and second an adjusted value, $A_{adjusted}$, where we rescale the rate with
experimental energies according to

$A_{adjusted}=\left(\frac{E_{exp}}{E_{ab-initio}}\right)^{3}A_{ab-initio}$ (2)

In this equation $E_{exp}$ and $E_{ab-initio}$ are experimental and
theoretical ab initio transition energies, respectively. This approach is
based on the non-relativistic limit, where the line strength, $\mathcal{S}$,
is energy independent and the rate is simply $A\propto
E^{3}\times\mathcal{S}$.

It is clear that the adjusted rate, $A_{adjusted}$, and corresponding lifetime
is converging fast and does not change as a function of neither the type of
correlation model nor its size, $n_{max}$. Even if it might be argued that the
transition energy is not fully converged, it is clear that there is little
chance that any additional contributions to the transition rate could affect
the final value of $A_{adjusted}$.

From a straightforward analytical and non-relativistic model for the M1
transition in the $LS$-approximation, it becomes clear why the convergence of
our results is so fast. The M1 is a forbidden transition, but expected Brage
_et al._ (2009); Grumer  _et al._ (2014a), since it is not induced by small
corrections to the wavefunction. More importantly, the rate of the M1
transition in this model, does not depend on the radial part of the
wavefunctions, but is rather given by Cowan (1981)

$\displaystyle A_{analytical}$ $\displaystyle=$
$\displaystyle\frac{64\pi^{2}e^{2}a_{0}^{2}(\alpha/2)^{2}\sigma^{3}}{3h(2J^{\prime}+1)}\left|\langle
LSJ||{\bf J}^{(1)}+\left(g_{s}-1\right){\bf
S}^{(1)}||L^{\prime}S^{\prime}J^{\prime}\rangle\right|^{2}$ (3)
$\displaystyle\approx$ $\displaystyle\frac{2.6973\times
10^{-11}\sigma^{3}}{2J^{\prime}+1}\left|\langle LSJ||{\bf
J}^{(1)}+\left(g_{s}-1\right){\bf
S}^{(1)}||L^{\prime}S^{\prime}J^{\prime}\rangle\right|^{2}$

in a general case of fine-structure, M1 transition, where the primed
quantities represent the initial, upper state. In our specific case with a
${}^{2}P_{1/2}\rightarrow^{2}P_{3/2}$ transition, we get

$\displaystyle A_{analytical}$ $\displaystyle\approx$
$\displaystyle\frac{2.6973\times
10^{-11}\sigma^{3}}{2}\left|\langle^{2}P_{3/2}||{\bf
J}^{(1)}+\left(g_{s}-1\right){\bf S}^{(1)}||^{2}P_{1/2}\rangle\right|^{2}$ (4)

where $\sigma=E/hc$ is the wavenumber for the transition (in cm-1). It can be
shown that

$\left|\langle^{2}P_{3/2}||{\bf J}^{(1)}+\left(g_{s}-1\right){\bf
S}^{(1)}||^{2}P_{1/2}\rangle\right|^{2}\approx\frac{4}{3}$ (5)

in pure $LS$-coupling. Finally, we find that

$A_{analytical}\approx 1.7982\times 10^{-11}\sigma^{3}.$ (6)

By using the experimental transition energy $\sigma=59.9507$ meV$/hc$ =
483.535 cm-1 Blondel  _et al._ (2006) we get $A_{analytical}=2.0328\times
10^{-3}$ s-1, in perfect agreement with our adjusted result and experiment.

Table 2:  Experimental wavenumbers for the ${}^{2}P_{1/2}\rightarrow^{2}P_{3/2}$ transitions ($\sigma_{exp}$), together with lifetimes for the upper ${}^{2}P_{1/2}$ level in four different negative ions. Lifetimes are reported from experiment ($\tau_{exp}$), Eq. (6) using $\sigma_{exp}$ values ($\tau_{analytical}$), and from calculations ($\tau_{calc}$). In addition, we also present values for ${}^{2}D_{3/2}\rightarrow\,^{2}D_{5/2}$ M1 transitions in two selected negative ions. Ion |  $\sigma_{exp}$ (cm${}^{-1})$ |  $\tau_{analytical}$ (s) |  $\tau_{exp}$ (s) |  $\tau_{calc}$ (s)
---|---|---|---|---
${}^{2}P$-terms
O- | 177.084 Andersen  _et al._ (1999) | $1.00\times 10^{4}$ |  |  $7.24\times 10^{3}$ Andersson  _et al._ (2006)
S- | 483.535 Blondel  _et al._ (2006) | 492 |  $503\pm 54$ Bäckström  _et al._ (2015) | 437 Andersson  _et al._ (2006), 492a
Se- | 2278.15 Andersson  _et al._ (2006) | 4.70 |  $4.78\pm 0.18$ Bäckström  _et al._ (2015) | 4.92 Andersson  _et al._ (2006)
Te- | 5005.36 Scheer  _et al._ (1998) | 0.443 |  $0.463\pm 0.008$ Bäckström  _et al._ (2015) | 0.454 Andersson  _et al._ (2006)
${}^{2}D$-terms
Ni- | 1484.8 Kamińska  _et al._ (2016) | 18.88 |  $15.1\pm 0.4$ Kamińska  _et al._ (2016) |
Pt- | 9740.9 Thøgersen  _et al._ (1996) | 0.0669 |  |  $0.071$ Thøgersen  _et al._ (1996)
aour MCDHF results presented in Tab. 1

This simple model can of course also be used for ${}^{2}P$ transitions in
other negative ions. Andersson et al. Andersson  _et al._ (2006) reported the
lifetimes of the single bound excited states in Se- ($4p^{5}\;P_{1/2}$) and
Te- ($5p^{5}\;P_{1/2}$) using a room temperature magnetic ion storage ring.
Recently, a breakthrough in ion-beam storage technology at cryogenic
temperatures was reported Schmidt  _et al._ (2013); Thomas  _et al._ (2011).
This new technology allows the storage of beams of negative ions for much
longer times than what has been possible before, which brings measurements of
atomic lifetimes of the order of tens of minutes within reach. This enabled
Bäckström et al. to measure the lifetime of the $3p^{5}\;{}^{2}P_{1/2}$ level
to be $503\pm 54$ s and to measure the lifetimes of the corresponding states
in Se- and Te- with considerably higher accuracy Bäckström  _et al._ (2015)
than in earlier studies.

One could also apply this analytical model to systems with different terms.
The $3d^{9}\,4s^{2}\,{}^{2}D_{3/2}\rightarrow 3d^{9}\,4s^{2}\,{}^{2}D_{5/2}$
transition in e.g. Ni- was recently observed by Kamińska et al. Kamińska  _et
al._ (2016). An anologous ${}^{2}D$ transition for $5d^{9}\,6s^{2}$
configurations was observed in Pt- by Thøgersen et al. Thøgersen  _et al._
(1996) (note that no experimental value for the lifetime was reported in this
work). In the ${}^{2}D$ case the rate is given by

$\displaystyle A_{analytical}$ $\displaystyle\approx$
$\displaystyle\frac{2.6973\cdot
10^{-11}\sigma^{3}}{4}\left|\langle^{2}D_{5/2}||{\bf
J}^{(1)}+\left(g_{s}-1\right){\bf S}^{(1)}||^{2}D_{3/2}\rangle\right|^{2}$ (7)

where the reduced matrix element is

$\left|\langle^{2}D_{5/2}||{\bf J}^{(1)}+\left(g_{s}-1\right){\bf
S}^{(1)}||^{2}D_{3/2}\rangle\right|^{2}\approx\frac{3}{5}$ (8)

for pure $LS$-coupled states. From Tab. 2 it is clear that the fact that the
rates are well-represented by the analytical expression of Eq. (6), is not
unique for S-, but a general case for at least the halogen-like transitions.

There are two possible sources of the uncertainty in the predicted lifetimes -
the experimental uncertainty in the measured transition energy and the
deviation from the $LS$-approximation for the negative ions studied here. The
former is for S- estimated to be less than 0.01 ‰ while the latter, according
to the stability of the adjusted value in Tab. 1, has an uncertainty of less
than 0.1 ‰. The uncertainty of the predicted lifetime is therefore well within
1 % for S-, but increasing with nuclear charge for the other systems in Tab.
2, due to an increasing deviation from the $LS$-approximation for heavier
ions.

##  IV Conclusions

We can conclude that the rate of the fine structure M1 transitions are
straightforward to compute for the negative ions considered here, while an
accurate determination of the wavelengths is dependent on a careful treatment
of correlation. We therefore recommend the use of the measured lifetime of
these upper levels as a support for the spectroscopic identification of the
corresponding fine-structure transitions. A similar method was proposed and
used for highly ionized species in what was labelled ”single-line
spectroscopy” for silver-like ions, e.g. Tungsten Grumer  _et al._ (2014b).

Measurements of the lifetimes of the upper states involved in these
transitions could be used to confirm their spectroscopic identification, This
could be an important tool since atomic anions often have just a few bound
states, and therefore only a few photon signals from bound-bound transitions,
thereby making traditional spectroscopic analyses impossible. In addition to
this, in cases where the lifetime, but not the transition energy, is known,
the latter can be determined from the former by using the simple energy
dependence of the lifetime.

##  V Acknowledgments

We gratefully acknowledge Dag Hanstorp, Henrik Cederquist and Henning T.
Schmidt for helpful comments and fruitful discussions. This work was supported
by the Swedish Research Council (VR) under Contract No. 2015-04842.

## References

  * Bäckström  _et al._ (2015) E. Bäckström, D. Hanstorp, O. M. Hole, M. Kaminska, R. F. Nascimento, M. Blom, M. Björkhage, A. Källberg, P. Löfgren, P. Reinhed, S. Rosén, A. Simonsson, R. D. Thomas, S. Mannervik, H. T. Schmidt, and H. Cederquist, Phys. Rev. Lett. 114, 143003 (2015).
  * Andersen (2004) T. Andersen, Phys. Rep. 394, 157 (2004).
  * Pegg (2004) D. J. Pegg, Rep. Prog. Phys. 67, 857 (2004).
  * Andersen  _et al._ (1999) T. Andersen, H. K. Haugen, and H. Hotop, J. Phys. Chem. Ref. Data 28, 1511 (1999).
  * Fischer (1989) C. F. Fischer, Phys. Scr. 40, 25 (1989).
  * Walter  _et al._ (2014) C. W. Walter, N. D. Gibson, D. J. Matyas, C. Crocker, K. A. Dungan, B. R. Matola, and J. Rohlén, Phys. Rev. Lett. 113, 063001 (2014).
  * Blondel  _et al._ (2006) C. Blondel, W. Chaibi, C. Delsart, and C. Drag, J. Phys. B: At. Mol. Opt. Phys. 39, 1409 (2006).
  * Andersson  _et al._ (2006) P. Andersson, K. Fritioff, J. Sandström, G. Collins, D. Hanstorp, A. Ellmann, P. Schef, P. Lundin, S. Mannervik, P. Royen, K. C. Froese Fischer, F. Österdahl, D. Rostohar, D. J. Pegg, N. D. Gibson, H. Danared, and A. Källberg, Phys. Rev. A 73, 032705 (2006).
  * Fischer  _et al._ (2007) C. F. Fischer, G. Tachiev, G. Gaigalas, and M. R. Godefroid, Comp. Phys. Comm. 176, 559 (2007).
  * Filippin  _et al._ (2016) L. Filippin, M. Godefroid, J. Ekman, and P. Jönsson, Phys. Rev. A 93, 062512 (2016).
  * Fischer  _et al._ (2016) C. F. Fischer, M. Godefroid, T. Brage, P. Jönsson, and G. Gaigalas, J. Phys. B: At. Mol. Opt. Phys. 49, 182004 (2016).
  * Froese Fischer  _et al._ (1997) C. Froese Fischer, T. Brage, and P. Jönsson,  _Computational atomic structure: an MCHF approach_ (Inst. of Physics Publishing, Bristol, 1997).
  * Grant (2007) I. P. Grant,  _Relativistic Quantum Theory of Atoms and Molecules: Theory and Computation (Springer Series on Atomic, Optical, and Plasma Physics)_ (Springer-Verlag New York, Inc., 2007).
  * Jönsson  _et al._ (2013) P. Jönsson, G. Gaigalas, J. Bieroń, C. F. Fischer, and I. Grant, Comp. Phys. Comm. 184, 2197 (2013).
  * Brage  _et al._ (2009) T. Brage, M. Andersson, and R. Hutton, AIP Conf. Proc. 1125, 18 (2009).
  * Grumer  _et al._ (2014a) J. Grumer, T. Brage, M. Andersson, J. Li, P. Jönsson, W. Li, Y. Yang, R. Hutton, and Y. Zou, Phys. Scr. 89, 114002 (2014a).
  * Cowan (1981) R. D. Cowan,  _The theory of atomic structure and spectra_ , Vol. 3 (University of California Press, 1981).
  * Scheer  _et al._ (1998) M. Scheer, R. C. Bilodeau, and H. K. Haugen, J. Phys. B: At. Mol. Opt. Phys. 31, L11 (1998).
  * Kamińska  _et al._ (2016) M. Kamińska, V. T. Davis, O. M. Hole, R. F. Nascimento, K. C. Chartkunchand, M. Blom, M. Björkhage, A. Källberg, P. Löfgren, P. Reinhed, S. Rosén, A. Simonsson, R. D. Thomas, S. Mannervik, P. A. Neill, J. S. Thompson, H. T. Schmidt, H. Cederquist, and D. Hanstorp, Phys. Rev. A 93, 012512 (2016).
  * Thøgersen  _et al._ (1996) J. Thøgersen, M. Scheer, L. D. Steele, H. K. Haugen, and W. P. Wijesundera, Phys. Rev. Lett. 76, 2870 (1996).
  * Schmidt  _et al._ (2013) H. T. Schmidt, R. D. Thomas, M. Gatchell, S. Rosén, P. Reinhed, P. Löfgren, L. Brännholm, M. Blom, M. Björkhage, E. Bäckström, J. D. Alexander, S. Leontein, D. Hanstorp, H. Zettergren, L. Liljeby, A. Källberg, A. Simonsson, F. Hellberg, S. Mannervik, M. Larsson, W. D. Geppert, K. G. Rensfelt, H. Danared, A. Paál, M. Masuda, P. Halldén, G. Andler, M. H. Stockett, T. Chen, G. Källersjö, J. Weimer, K. Hansen, H. Hartman, and H. Cederquist, Rev. Sci. Instrum. 84, 055115 (2013).
  * Thomas  _et al._ (2011) R. D. Thomas, H. T. Schmidt, G. Andler, M. Björkhage, M. Blom, L. Brännholm, E. Bäckström, H. Danared, S. Das, N. Haag, P. Halldén, F. Hellberg, A. I. S. Holm, H. A. B. Johansson, A. Källberg, G. Källersjö, M. Larsson, S. Leontein, L. Liljeby, P. Löfgren, B. Malm, S. Mannervik, M. Masuda, D. Misra, A. Orbn, A. Paál, P. Reinhed, K.-G. Rensfelt, S. Rosén, K. Schmidt, F. Seitz, A. Simonsson, J. Weimer, H. Zettergren, and H. Cederquist, Rev. Sci. Instrum. 82, 065112 (2011).
  * Grumer  _et al._ (2014b) J. Grumer, R. Zhao, T. Brage, W. Li, S. Huldt, R. Hutton, and Y. Zou, Phys. Rev. A 89, 062511 (2014b).