How the Magnetization Angle of a Linear Halbach Array Influences Particle Steering in Magnetic Drug Targeting—A Systematic Evaluation and Optimization

Abstract

The main challenge in magnetic drug targeting lies in steering the magnetic particles, especially in deeper body layers. For this purpose, linear Halbach arrays are currently in focus. However, to the best of the authors’ knowledge, the impact of the magnetization angle between two neighboring magnets in Halbach arrays has not been investigated for particle steering so far. Therefore, in this paper, a systematic numerical parameter study of varying the magnetization angle of linear Halbach arrays is conducted. This is completed by undertaking a typical magnetic drug targeting scenario, where magnetic particles have to be steered in an optimized manner. This includes the calculation of the magnetic flux density, its gradient, the total magnetic energy, and the resulting magnetic force based on a fitting function for the different Halbach constellations in the context of examining their potential for predicting the particle distribution. In general, increased magnetization angles result in an increased effective range of the magnetic force. However, as there is a trade-off between a weak force on the weak side of the array and a simple manufacturing process, a magnetization angle of ${90}^{\circ }$ is recommended. For evaluating the steering performance, a numerical or experimental evaluation of the particle distribution is mandatory.

1. Introduction

The interest in steering, separating, or trapping magnetic nanoparticles or cells, especially in biological and medical applications, has grown rapidly in recent decades [1,2,3,4]. One of these applications is so-called magnetic drug targeting (MDT) in cancer therapy [5]. This method is based on superparamagnetic iron oxide nanoparticles (SPIONs) as carriers to which the anticancer agents are attached. These SPIONs are injected into the vascular system and conveyed by the blood flow to the targeted region due to the influence of an external applied magnetic field [6]. In the case of lung cancer, the SPIONs are inhaled [7]. MDT considerably decreases the overall dose and side effects on healthy tissues since SPIONs coated with drugs are concentrated in the targeted region and maintained there due to a holding magnetic field [8]. In other words, MDT improves the efficiency of cancer therapy as it has a local cancer treatment effect [9,10]. Moreover, it reduces unwanted adverse effects, which was experimentally proven in 2013 by Tietze et al. [11]. A typical MDT scenario and its modelling are illustrated in Figure 1. An extensive review of MDT approaches, release control of the particles, and an overview of the existing clinical studies are provided by Ulbrich et al. [12]. The drug is released at the tumor region over time, or an alternating magnetic field is applied to the particles, resulting in heating, which is known as hyperthermia [13]. However, the success of such a treatment depends on the number of SPIONs reaching the tumor and being trapped there. For this purpose, most researchers think of placing permanent magnets [14,15] or electromagnets [16,17,18] outside the body. Since the usable reach of the magnetic field and especially of the magnetic force is a significant problem [19], it is also being considered to place the patient inside a magnetic resonance imaging (MRI) scanner [20], offering a strong constant field plus switchable overlay fields, or even implanting permanent magnets next to the tumor region [21,22].

No matter the specific magnetic field design, there is always a trade-off between too strong a magnetic field, which causes all SPIONs to stick to the vessel wall, or a magnetic field that is too weak to attract the SPIONs [15]. Figure 2 illustrates that, due to the background velocity flow field, the particles in the upper half of the vessel automatically take the correct outlet in a symmetrical bifurcation. Thus, the particles injected in the lower half of the vessel must be pulled into the upper part while simultaneously preventing the particles already in the upper part to become accumulated and remain stuck at the vessel wall. As the magnetic field approximately decays with the cubic distance $R^3$ [23], this worsens steering performance since the magnetic force is stronger closer to the magnet. Moreover, accumulated particles show a shielding behavior, reducing the magnetic forces [24].

To overcome this problem, many researchers use electromagnets (EMs) on two sides of a tube, which are alternately switched on and off to steer particles through a bifurcation [25,26,27] and optimize the switching algorithm by defining safety zones [28,29,30,31]. Moreover, Le et al. [32] used four EMs to generate a focused field-free point that can be shifted in place to avoid stuck particles at the vessel wall. However, EMs require an additional power supply in the kW range, with currents greater than 200 A inducing heat, therefore requiring additional cooling [33]. On the other hand, rare earth permanent magnets (PMs), like NdFeB, generate a low-cost and strong magnetic field. Thus, they are studied for attracting or steering magnetic nanoparticles too [15,34,35]. However, permanent magnets have the disadvantage that their magnetic field can only be changed by mechanical operations (e.g., rotating or changing the position). Thus, Surpi et al. [36] investigated the behavior of SPIONs under two moving PMs placed on opposite sides of a tube, and Baun et al. [37] precisely steered SPIONs in a static fluid with a cylindrical array of PMs surrounding the particle container. A more detailed review of the current state-of-the-art publications for steering magnetic nanoparticles in an MDT scenario is provided and discussion in Section 5.4.

Nevertheless, the accumulation of the SPIONs in the tumor and, thus, the efficiency of an MDT treatment depend on numerous parameters, like the flow of blood with the corresponding geometry of the vessel and tumor, the properties of the SPIONs, as well as the gradient of the magnetic field [3,38]. As we cannot adjust the blood flow and the nanoparticles have to have a diameter <200 nm to not be recognized by the body’s immune system and be able to penetrate into tissue [39], only the magnetic field can be appropriately designed. Here, a strong gradient field over a large distance is mandatory [38,40]. However, our previous study showed that, for the same magnetic effort, a longer magnetic array with a lower absolute magnetic field strength leads to the same particle attraction as a shorter array with a higher field strength [41].

In the literature, different types of arrays, mostly entirely consisting of PMs [42,43] but also of EMs [44], as well as hybrid ones [45,46], are discussed. The arrays are arranged linearly and parallel to the vessel [47] as well as coaxial structures surrounding the vessel [37]. As in medical applications placing the magnets surrounding the vessel is not possible, this paper will focus on linear arrays. In most of the literature, the magnets are arranged as a Halbach array [19], which is a promising structure since its pattern produces the strongest magnetic flux density [48]. Kang et al. [49] showed that the attraction force of a Halbach structure is stronger than when the single magnets are magnetized in alternating directions. Häfeli et al. [50] compared different array shapes and also concluded that Halbach arrays have the best particle trapping close to a surface. However, in [50], a Halbach structure with a shifted magnetization direction of 45$°$ between two neighboring magnets was investigated. In our previous publication [51], we demonstrated that the attraction force is even higher when the neighboring magnets have a shifted magnetization of 90$°$. Also, in most of the literature, 90$°$ is used for steering or attracting magnetic nanoparticles [47,49,52,53,54]. In other applications of Halbach arrays, various magnetization angles are investigated and discussed, e.g., 45$°$ by Zhang et al. [55], 30$°$ by Di Gerlando et al. [56], 18$°$ and 22.5$°$ by Shen et al. [57], and 90$°$ by Li et al. [58].

To the best of the authors’ knowledge, the optimal magnetization angle between two adjacent magnets in a Halbach array for the application of attracting magnetic nanoparticles has not been systematically analyzed so far. With this paper, we want to close this gap in the state of the art. Our key contributions are summarized:

• We systematically numerically investigate the impact of the magnetization angle in a reduced complexity 2D simulation in COMSOL Multiphysics^® 6.1 regarding the steering performance of SPIONs in a background flow through a Y-shaped bifurcation.
• We further evaluate the magnetic force for the strong and weak side of the array for different magnetization angles at its strongest position. We additionally compare the magnitude of the magnetic forces with the hydrodynamic drag force.
• Since the calculation of the magnetic gradient leads to huge errors due to the discrete mesh in COMSOL, it was determined analytically using an exponential fitting function similar to [51]. By doing so, we further analyze if the magnetic flux density can be approximated with an exponential function for other magnetization angles too.
• Since we were not able to identify standardized evaluation parameters with our comprehensive analysis of the state of the art (compare Section 5.4), the significance of various evaluation parameters, such as the magnetic field, total applied magnetic energy, or maximum gradient in the vessel, is examined regarding their prediction of particle steering.
• Based on our investigations, recommendations for the design of a linear Halbach array for steering SPIONs in MDT are derived and discussed.

Our paper is structured as follows: after the introduction in Section 1, the fundamentals and physical background of SPIONs, (linear) Halbach arrays, and their magnetic flux density are provided in Section 2. In Section 3, the simulation model is defined, and the evaluation process of this study is introduced. Afterwards, the results are presented in Section 4. Finally, in Section 5, the results are discussed regarding particle steering, including the strength and usable range of the magnetic field and force. Furthermore, recommendations for the design of a linear Halbach array are provided and the results of this paper are compared with the state of the art of particle steering in an MDT scenario. Moreover, the limitations of this study are identified. The paper is concluded in Section 6 and an outlook of future steps is provided.

2. Fundamentals and Background

2.1. Superparamagnetic Iron Oxide Nanoparticles

In MDT, superparamagnetic iron oxide nanoparticles (SPIONs) with a diameter of <250 nm are usually used as carriers for cancer drugs [39,59]. The SPIONs are coated, e.g., by lauric acid layers, to be biocompatible [60] and stable in blood [39]. Furthermore, they can be synthesized small enough not to be attacked by the body’s immune system [6,61].

Superparamagnetism exists in small ferromagnetic particles [62]. These particles consist only of single-domain magnetization, which flips randomly in time, resulting in an average magnetization to be zero [39,62]. Therefore, SPIONs have strong magnetization with no hysteretic behavior, and, consequently, the particles do not accumulate [61]. This is advantageous for medical applications since the risk of obstructing blood vessels is reduced. Lunnoo et al. [63] analyzed the size-dependent capture efficiency of SPIONs by a magnet at a desired location depending on magnetic cores made from different materials, such as Fe₃O₄, Fe₂O₃, and Fe nanoparticles with different coating layers, such as Au, PEG, and SiO₂. They focused on particle sizes ranging from 10 nm to 200 nm, which is seen as optimal for in vivo delivery [6]. Their results showed that larger particles had better capture efficiency than smaller particles. This is reasonable as the magnetic force is directly proportional to the volume of the particles [40]. The material of magnetic cores and coating layers did not significantly influence the capture efficiency.

2.2. Forces on SPIONs

Using Ampere’s law $\nabla \times \overrightarrow{H}=\overrightarrow{J}$ in the quasi-static case, and magnetostatic field dependencies $\overrightarrow{B}=\nabla \times \overrightarrow{A}$, the equation of the magnetic field can be described as [64]

$$\nabla \times (\frac{1}{{\mu }_0}\nabla \times \overrightarrow{A}-\overrightarrow{M})=0.$$

(1)

Here, ${\mu }_0=4\pi ·{10}^{-7}$ H/m corresponds to the permeability of the free space, $\overrightarrow{A}$ is the magnetic vector potential (Vs/m), and $\overrightarrow{M}$ is the volume magnetization vector (A/m). Furthermore, $\overrightarrow{J}$ is the current density (A/m²), $\overrightarrow{B}$ is the magnetic flux density (T), and $\overrightarrow{H}$ is the magnetic field strength (A/m).

The theory of magnetic flux conservation in the human body can be represented using the relation between $\overrightarrow{B}$ and $\overrightarrow{H}$:

$$\overrightarrow{B}=\left\{\begin{array}{cc}{\mu }_0{\mu }_{\mathrm{r}}\overrightarrow{H}\hfill & \mathrm{blood},\mathrm{tissue},\mathrm{and}\mathrm{air}\hfill \\ {\mu }_0(\overrightarrow{H}+\overrightarrow{M})\hfill & \mathrm{PM}\hfill \end{array}\right.$$

(2)

Here, ${\mu }_{\mathrm{r}}$ describes the relative permeability of a medium. For example, for blood and air, ${\mu }_{\mathrm{r}}$ is approximately 1 [65]. Aside from the magnetic field, the blood flow and properties of SPIONs also play an essential role in the performance of particle steering. To determine the velocity field and the pressure distribution of the background flow, the equation of mass and the Navier–Stokes equation for a Newtonian and incompressible fluid are used [63,66]

$$\begin{array}{cc}\hfill \rho \nabla ·\overrightarrow{u}& =0,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \rho (\overrightarrow{u}·\nabla \overrightarrow{u})& =-\nabla p+\eta \nabla ·(\nabla \overrightarrow{u})+\overrightarrow{F}.\hfill \end{array}$$

(4)

Thereby, $\rho$ corresponds to the mass density of the background flow (kg/m³), $\overrightarrow{u}$ is the velocity vector (m/s), p is pressure (Pa), $\eta$ is the fluid’s viscosity (Pa·s), and $\overrightarrow{F}$ is the external volume force (N/m³). For the sake of simplicity, in this paper, the medium of the background flow is assumed to be water since the flow properties of blood are much more complex. The interested reader can find detailed information regarding the mathematical modeling of magnetic nanoparticles in a blood flow for MDT in the work of Fanelli et al. [67].

To represent the particle tracing process, Newton’s second law is used to describe the balance of forces on one SPION [66]

$$\sum {\overrightarrow{F}}_{\mathrm{p}}=\frac{\mathrm{d}\left(m_{\mathrm{p}}{\overrightarrow{u}}_{\mathrm{p}}\right)}{\mathrm{d}t},$$

(5)

where $m_{\mathrm{p}}$ corresponds to the mass of the particle (kg), and ${\overrightarrow{u}}_{\mathrm{p}}$ is the velocity of the particle (m/s). ${\overrightarrow{F}}_{\mathrm{p}}$ is a force acting on the particle (N).

As the SPIONs are small in size, when one SPION is moving through the blood vessel, it mainly experiences two forces, namely the magnetophoretic force ${\overrightarrow{F}}_{\mathrm{m}}$ and the hydrodynamic drag force ${\overrightarrow{F}}_{\mathrm{d}}$ [3,68,69]. Other forces such as gravitational force ${\overrightarrow{F}}_{\mathrm{g}}$, buoyancy force ${\overrightarrow{F}}_{\mathrm{b}}$, and lift force ${\overrightarrow{F}}_{\mathrm{l}}$ can be neglected because they are significantly weaker compared to ${\overrightarrow{F}}_{\mathrm{m}}$ and ${\overrightarrow{F}}_{\mathrm{d}}$ in MDT [3,66,70]. As a result, Equation (5) can be rewritten as

$${\overrightarrow{F}}_{\mathrm{m}}+{\overrightarrow{F}}_{\mathrm{d}}=\frac{\mathrm{d}\left(m_{\mathrm{p}}{\overrightarrow{u}}_{\mathrm{p}}\right)}{\mathrm{d}t}.$$

(6)

The hydrodynamic drag force is also called the viscous drag force or only drag force [71]. This force describes the resistance when a spherical particle is moving through viscous fluid. For low-flow regimes with low Reynolds numbers Re, a laminar flow with a parabolic flow profile is obtained [3,72]. According to Stokes’ law, ${\overrightarrow{F}}_{\mathrm{d}}$ is calculated as [27,37]

$${\overrightarrow{F}}_{\mathrm{d}}=-6\pi \eta R_{\mathrm{H}}({\overrightarrow{u}}_{\mathrm{p}}-{\overrightarrow{u}}_{\mathrm{b}}),$$

(7)

where $R_{\mathrm{H}}$ is the hydrodynamic radius of the particle (m), $\eta$ corresponds to the dynamic viscosity of the fluid medium (Pa·s), and ${\overrightarrow{u}}_{\mathrm{b}}$ and ${\overrightarrow{u}}_{\mathrm{p}}$ are the velocity of the background fluid (m/s) and the particle (m/s), respectively. For sake of simplicity, in this paper, the hydrodynamic radius $R_{\mathrm{H}}$ of one SPION is set to be equal to the magnetic radius of the particle $R_{\mathrm{p}}$. In MDT, the magnitude of $|{\overrightarrow{F}}_{\mathrm{d}}|$ is in the range of ${10}^{-14}$ to ${10}^{-7}$ N and the velocity of blood $0.5$ mm/s to 40 cm/s, respectively [1].

In general, on a microscopic scale and for a constant temperature, the magnetic force ${\overrightarrow{F}}_{\mathrm{m}}$ acting on a magnetic material at a certain position $\overrightarrow{r}$ is defined by [23,37]

$${\overrightarrow{F}}_{\mathrm{m}}\left(\overrightarrow{r}\right)=\nabla \left(\overrightarrow{m}\left(\overrightarrow{r}\right)·\overrightarrow{B}\right),$$

(8)

with the magnetic flux density $\overrightarrow{B}$ acting on a material of the magnetic moment $\overrightarrow{m}$ ($\mathrm{A}·{\mathrm{m}}^2$). As the magnetic moment $\overrightarrow{m}$ inside one SPION aligns parallel to $\overrightarrow{B}$, the magnetic force can also expressed as

$${\overrightarrow{F}}_{\mathrm{m}}\left(\overrightarrow{r}\right)=\nabla \left(\overrightarrow{m}\left(\overrightarrow{r}\right)·\overrightarrow{B}\right)=\nabla \left(m·B\right).$$

(9)

Moreover, Equation (8) can be rewritten as [40]

$${\overrightarrow{F}}_{\mathrm{m}}=\nabla \left(\overrightarrow{m}·\overrightarrow{B}\right)=\left(\overrightarrow{B}·\nabla \right)\overrightarrow{m}+\left(\overrightarrow{m}·\nabla \right)\overrightarrow{B}+\overrightarrow{B}\times \left(\nabla \times \overrightarrow{m}\right)+\overrightarrow{m}\times \left(\nabla \times \overrightarrow{B}\right).$$

(10)

For SPIONs, the magnetic moment $\overrightarrow{m}$ is usually assigned to be homogeneous (no spatial dependency) within one SPION. This leads to $\left(\overrightarrow{B}·\nabla \right)\overrightarrow{m}=0$. In [40], it was stated that, for particle steering, the SPIONs can be assumed to be magnetic dipoles. This cancels $\left(\nabla \times \overrightarrow{m}\right)=0$ because $\overrightarrow{m}$ is not a function of $\overrightarrow{B}\left(\overrightarrow{r}\right)$. As in this paper, permanent magnets are used as magnetic field sources; the current density $\overrightarrow{J}$ is equal to zero. This results in $\nabla \times \overrightarrow{B}=0$. Thus, ${\overrightarrow{F}}_{\mathrm{m}}$ can be written as

$${\overrightarrow{F}}_{\mathrm{m}}=\left(\overrightarrow{m}·\nabla \right)\overrightarrow{B},$$

(11)

which is also known as Kelvin force [40]. The magnetic moment $\overrightarrow{m}$ can be deduced from the integral of the magnetization $\overrightarrow{M}$ (A/m) over the corresponding volume V (${\mathrm{m}}^3$) [23]. For SPIONs, the magnetization is usually assigned to be homogeneous within one SPION and provided by

$$M=M_{\mathrm{sat}}·\mathcal{L}\left({\displaystyle \frac{\overrightarrow{m}·\overrightarrow{B}}{kT}}\right),$$

(12)

where $M_{\mathrm{sat}}$ corresponds to the saturation magnetization (A/m), $k\approx 1.3806·{10}^{-23}$ J/K is the Boltzmann constant, and T (K) the absolute temperature [73,74]. Moreover, $\mathcal{L}$ is the Langevin function, which has a strictly monotonously increasing behavior.

If the magnetic force is acting on microscale entities, it is called magnetophoretic force [3,75]. By solving the magneto-static boundary value problem for one spherical particle, the magnetophoretic force ${\overrightarrow{F}}_{\mathrm{m}}$ can be formulated as [23,69]

$${\overrightarrow{F}}_{\mathrm{m}}=2\pi {\mu }_{\mathrm{b}}{R}_{\mathrm{p}}^3{\displaystyle \frac{{\mu }_{\mathrm{p}}-{\mu }_{\mathrm{b}}}{{\mu }_{\mathrm{p}}+2{\mu }_{\mathrm{b}}}}\nabla {\left|\overrightarrow{H}\right|}^2.$$

(13)

Thereby, ${\mu }_{\mathrm{b}}$ and ${\mu }_{\mathrm{p}}$ are the magnetic permeability of the background medium and the particle, respectively. Moreover, $R_{\mathrm{p}}$ is the radius of the particle and $\nabla \overrightarrow{H}$ is the gradient of the magnetic field excitation. Assuming a homogeneous magnetization within one SPION that is aligned parallel to the applied magnetic field, and based on Equation (11), the magnitude of the magnetophoretic force $|{\overrightarrow{F}}_{\mathrm{m}}|$ can be simplified to [51]

$$|{\overrightarrow{F}}_{\mathrm{m}}|=|\left(\overrightarrow{m}·\nabla \right)\overrightarrow{B}|=m|\nabla B|=V_{\mathrm{p}}M|\nabla B|.$$

(14)

Thus, $|{\overrightarrow{F}}_{\mathrm{m}}|$ only depends on the SPION’s volume $V_{\mathrm{p}}$, the magnetization M, and the gradient of the magnetic flux density $\nabla B$. In the application of MDT, the magnitude $|{\overrightarrow{F}}_{\mathrm{m}}|$ can obtain values from ${10}^{-25}$ to ${10}^{-11}$ N, depending especially on $\nabla B$ [1].

In MDT, ${\overrightarrow{F}}_{\mathrm{m}}$ and ${\overrightarrow{F}}_{\mathrm{d}}$ can usually be assumed to be perpendicular (compare Figure 2). Since ${\overrightarrow{F}}_{\mathrm{d}}$ is determined by the velocity flow of the blood, which is not adjustable, ${\overrightarrow{F}}_{\mathrm{m}}$ has to be maximized, which is the case for a Halbach arrangement.

2.3. Linear Halbach Array

A Halbach array is a special magnetic array in which the array’s single magnets are arranged so that the magnetic flux density is destructively and constructively superimposed on opposite sides of the array [76,77]. An illustration of a linear Halbach array arrangement is depicted in Figure 3. The magnetization pattern results in the ideal case in a one-sided magnetic field. Kang et al. [49] showed that the magnetic flux density and its gradient on the strong side of the array and, therefore, the magnetic force are stronger compared to a pattern where the magnetization direction of the single magnets point in alternating directions.

Halbach arrays were proposed by K. Halbach in 1980 for focusing charged particles in accelerator beams [77,78]. However, these arrays are found in the current literature mainly for shaping the fields of electric motors [79,80,81], but also in other various applications such as energy harvesting [82], induction-based cooking [83], refrigerators [84], or in medical applications such as magnetic resonance imaging (MRI) [85] or just as well in MDT [19].

For particle steering, most Halbach arrangements are linear and consist entirely of PMs (e.g., [43,47,49,51,52,54]). However, there are also 3D Halbach arrays consisting of EMs [86], which are combined with a core to focus SPIONs at a desired location (usually the tumor), or 3D constructs [37,87], which enclose the vessel and can precisely steer SPIONs here. However, since it is impossible to place the magnets around the blood vessel in reality, this paper focuses on linear arrays placed sideways to a vessel. The pattern of a linear Halbach array with a magnetization angle $\vartheta ={90}^{\circ }$ between two neighboring magnets is depicted in Figure 3. Hilton et al. [88] analyzed the torque ($\mathrm{N}·\mathrm{m}$) $\overrightarrow{\tau }=\overrightarrow{m}\times \overrightarrow{B}$ [23] of a Halbach array, as depicted in Figure 3, and figured out that, for a practical implementation, mechanical stabilization is crucial. Moreover, our paper also investigates other Halbach array patterns with varying magnetization angles between two neighboring magnets. An overview of different investigated array patterns is provided in Figure 4.

2.4. Magnetic Flux Density of a Linear Halbach Array

According to the literature, a Halbach array’s magnetic flux density $\overrightarrow{B}$ is computed with numerical methods or programs like COMSOL Multiphysics^® 6.1 [88,89]. Nevertheless, in our previous paper [51], we showed that the magnitude of the magnetic flux density perpendicular to a linear Halbach array with a magnetization angle of $\vartheta ={90}^{\circ }$ can be approximated with an exponential function

$$B\widehat{=}{k}_1{\mathrm{e}}^{-k_{2}y}.$$

(15)

For an array with $N=8$ magnets, this exponential approach describes B with a relative error of approximately 3.85% [51]. In Equation (15), $k_1$ and $k_2$ correspond to arbitrary fitting constants, which can be determined using a least square (LS) algorithm. In the scope of this paper, we want to investigate if the magnetic flux density of the Halbach array can also be approximated for other magnetization angles by this exponential function. Thus, the constants $k_1$ and $k_2$ are fitted based on the simulation results for both the strong and the weak side of the Halbach array for varying magnetization angles.

3. Definition of Simulation Model and Data Evaluation

3.1. Model Definition

In this article, COMSOL Multiphysics^® 6.1 is used to solve the magnetic field, fluid flow, and particle tracing equations in order to determine the steering performance in a 2D MDT scenario. To generate the magnetic flux density, a Halbach array consisting of rare-earth NdFeB PMs is utilized. The magnetic field is computed in a rectangular air-filled domain with a perfectly matched layer (PML) around it. The velocity flow consists of water and is assumed to be a Newtonian, incompressible, and single-phase laminar flow within a symmetrical Y-shaped bifurcation with an angle of ${30}^{\circ }$. One pulse of 100 SPIONs is released at the time $t=0$ s. The SPIONs are uniformly distributed across the cross-section, corresponding to the distribution in real scenarios [90]. However, 2 mm each are left out at the top and bottom of the inlet since the background velocity is equal to zero at the vessel walls, and the SPIONs, therefore, also move very slowly close to the vessel walls. By doing so, the simulation time is shortened significantly. The simulation geometry is shown in Figure 5, and the fixed simulation parameters are summarized in

Table 1. Summary of the reference simulation parameters.

Category	Symbol	Value	Unit	Label
	$B_{\mathrm{r}}$	1.2	T	Remanent flux density
	N	8	1	Number of PMs in the array
NdFeB	$R_{\mathrm{m}}$	1	cm	Radius of one single magnet
	$d_{\mathrm{m}}$	0.25	cm	Distance between single magnets
	$\vartheta$	90	$°$	Angle between magnetization direction
	$L_{\mathrm{v}}$	30	cm	Length
	$R_{\mathrm{v}}$	0.5	cm	Radius
Vessel	d	0.5	cm	Distance between magnets and vessel wall
	${\overrightarrow{u}}_{\mathrm{mean}}$	1	cm/s	Average background velocity
	I	$\left[0.2,0.8\right]$	cm	Inlet range of SPIONs
	$\rho$	2200	kg/m3	Mass density
SPIONs	$N_{\mathrm{p}}$	100	1	Number of simulated SPIONs
	$R_{\mathrm{p}}$	250	nm	Radius of one single SPION
	${\mu }_{\mathrm{r}}$	10	1	Relative permeability

. Moreover, the magnitude of the magnetic flux density $\left|B\right|$ is investigated in detail at its strongest position for the strong side of the array. This position is highlighted in Figure 5.

To implement the model, three interfaces in COMSOL are necessary:

(1)

the “Magnetic Fields, No Currents (mfnc)” interface is used for generating the magnetic flux density;

(2)

the “Laminar Flow (spf)” interface is related to the velocity flow;

(3)

the “Particle Tracing for Fluid Flow (fpt)” interface is capable of mimicking particle trajectories in the velocity flow.

In this model, particle–particle and particle–fluid interactions are neglected. It is also assumed that attracted SPIONs do not influence the magnetic field distribution.

3.2. General Fluid Mechanics Model Considerations

To estimate the behavior of a velocity flow in tubes, channels, and vessels, two characteristic numbers are often calculated, namely the Reynolds number Re and the Péclet number Pe [72,91]. The dimensionless Reynolds number $\mathrm{Re}\in {\mathbb{R}}^{+}$ can be used to predict whether the flow profile is laminar or turbulent. It is defined by

$$\mathrm{Re}={\displaystyle \frac{\mathrm{Inertial}\mathrm{force}}{\mathrm{Viscous}\mathrm{force}}}={\displaystyle \frac{u_{\max }{R}_{\mathrm{v}}}{\eta }},$$

(16)

where $u_{\max }$ corresponds to the maximum velocity in the vessel (m/s), $R_{\mathrm{v}}$ is the radius of the vessel (m), and $\eta$ is the the kinematic viscosity of the fluid’s medium (m²/s) [91,92]. For pure water, $\eta$ is equal to ${10}^{-6}$ m²/s [92]. For $\mathrm{Re}\gg 2100$, the flow is assumed to be turbulent and for $\mathrm{Re}\ll 2100$ laminar, respectively [91]. With $u_{\max }=2·u_{\mathrm{mean}}$ and the parameters in Table 1, we obtain $\mathrm{Re}=100\ll 2100$. This is in good accordance with the literature as for most blood vessels $\mathrm{Re}<500$ [91,93]. Therefore, in this paper, the velocity profile is expected to be fully laminar.

The movement of particles in a fluid is determined by both advective and diffusive flow. Which effect dominates can be estimated using the dimensionless Péclet number $\mathrm{Pe}\in {\mathbb{R}}^{+}$ [91,94] defined by

$$\mathrm{Pe}={\displaystyle \frac{\mathrm{Advective}\mathrm{transport}}{\mathrm{Diffusive}\mathrm{transport}}}={\displaystyle \frac{2·u_{\mathrm{mean}}{R}_{\mathrm{v}}}{D}}.$$

(17)

Here, $u_{\mathrm{mean}}$ is the average background velocity (m/s), $R_{\mathrm{v}}$ is again the radius of the vessel (m), and D is the diffusion coefficient (m²/s). In the case of

$$\mathrm{Pe}\gg {\displaystyle \frac{L_{\mathrm{v}}}{2R_{\mathrm{v}}}},$$

(18)

advective flow dominates [95]. According to [92], D of SPIONs is less than ${10}^{-11}$ m²/s. With Equation (17), Pe is less than ${10}^7$, which is much greater than $L_{\mathrm{v}}/\left(2R_{\mathrm{v}}\right)=30$. Therefore, advective transport is dominant for the chosen parameter setup of this paper.

3.3. Evaluation Procedure

In this study, the magnetic flux density with the resulting SPION distribution for a typical MDT scenario for different Halbach array scenarios is studied in order to evaluate the impact factors on the particle steering and find an optimized setup. By doing so, different parameters and scenarios were evaluated, which are described in detail in the following.

3.3.1. Evaluation of the Magnetic Flux Density and Its Gradient

At first, the distribution of the magnetic flux density B for the different magnetization angles $\vartheta \in \{{30}^{\circ },{45}^{\circ },{60}^{\circ },{65}^{\circ },{75}^{\circ },\dots ,{115}^{\circ },{120}^{\circ },{135}^{\circ }\}$ between the single magnets is investigated in detail. The pattern of the arrays always starts with a magnetization direction in the negative y-direction. The number of magnets was set to $N=8$, and the distance between the array and the vessel was chosen as 0.5 cm. In addition, two simulations with $\vartheta ={180}^{\circ }$ were conducted. Here, the magnetization direction was in the $\pm y$-direction and $\pm x$-direction, respectively. The magnetization pattern for these simulations is illustrated in Figure 4.

In the next step, the array’s strong and its weak side are compared by evaluating the absolute value of the magnetic flux density in the vertical direction of the array at its strongest position. For an odd number of magnets, this is always underneath the middle magnet, and, for an even number, under the middle two magnets. In this paper, the magnitude of B was evaluated for $N=8$ magnets in more detail. Thus, $\left|B\right|$ was extracted underneath the fourth magnet, as highlighted in Figure 5.

Moreover, the extracted B underneath the fourth magnet is fitted for every magnetization angle $\vartheta$ for the strong and weak side to the exponential function depicted in Equation (15). By doing so, the absolute and relative fitting errors, as well as the fitting parameters $k_1$ and $k_2$, are evaluated. From the fitting results, the gradient grad(B) is calculated for the strong and the weak side for every $\vartheta$ as the grad(B) is essential for the particle steering as the magnetic force is directly proportional to it (compare Equation (8)).

3.3.2. Evaluation of the SPION Distribution

From the perspective of magnetic drug targeting, the aim is to make the magnetophoretic force acting on particles stronger. As a result, the magnet attracts more particles to the desired region. Thus, for optimizing the Halbach array structure for a typical MDT scenario, the most important performance factor is the SPION distribution to the upper and lower branches, as well as the number of particles trapped by the magnets. In general, the stronger the force, the more particles take the upper branch or are trapped by the magnets and the less particles take the lower branch. Therefore, in this study, the SPION distribution on the upper and lower branches, as well as the ones trapped by the magnets for the geometry depicted in Figure 5, were investigated for different parameters of the linear Halbach array. The studied parameters are summarized in

Table 2. Investigated simulation parameters used in the study of the SPION distribution.

Symbol	Value	Unit	Label
N	$6\dots 12$	1	Number of PMs in the array
$\vartheta$	${30}^{\circ },{45}^{\circ },{60}^{\circ },{65}^{\circ },{75}^{\circ },\dots ,{115}^{\circ },{120}^{\circ },{135}^{\circ }$	$°$	Angle between magnetization direction
d	0.25; 0.5; 1; 1.5; 2	cm	Distance between magnets and vessel wall

.

As the simulation time of the particle study is quite long, and much computational effort is necessary, we want to check if another magnetic parameter can be conducted in order to predict the particle distribution. Thus, next to the SPION distribution, additionally, the maximum gradient in the vessel grad(B)_max and the total applied magnetic energy inside the vessel are derived for each investigated scenario. For the investigated 2D scenario in Figure 5, the total applied magnetic energy $W_{\mathrm{m}}$ (J/m) is provided by

$$W_{\mathrm{m}}=\frac{1}{2}\int \int \langle \overrightarrow{B},\overrightarrow{H}\rangle \mathrm{d}{A}_{\mathrm{vessel}}\widehat{=}\frac{1}{2}\int \int \frac{B^2}{{\mu }_0}\mathrm{d}{A}_{\mathrm{vessel}},$$

(19)

while integrating over the whole area cross-section of the vessel. Since the relative permeability of tissue and blood is approximately 1 [65], only the permeability of vacuum ${\mu }_0$ is taken into account for the medium of the background flow.

3.3.3. Evaluation of the Magnetic and Hydrodynamic Drag Force

As the movement of SPIONs through the vessel mainly depends on the forces acting on the SPIONs, the two strongest forces $F_{\mathrm{m}}$ and $F_{\mathrm{d}}$ are investigated in this paper. These forces are evaluated using the extracted magnetic flux density as depicted in Figure 5 for $N=8$ magnets and the vessel at a distance to the magnets of 0.5 cm. The magnetization angle $\vartheta$ was again varied in the range of $\vartheta \in \{{30}^{\circ },{45}^{\circ },{60}^{\circ },{65}^{\circ },{75}^{\circ },\dots ,{115}^{\circ },{120}^{\circ },{135}^{\circ }\}$. Aside from the forces, the ratio

$$F_{\mathrm{d}-\mathrm{m}}={\displaystyle \frac{|F_{\mathrm{d}}|}{|F_{\mathrm{m}}|}}$$

(20)

is analyzed inside the vessel. The evaluation of the forces $F_{\mathrm{m}}$ and $F_{\mathrm{d}}$ is described in the following:

(1)

Evaluation of the (hydrodynamic) drag force $F_{\mathrm{d}}$:

The drag force is calculated in the main branch using Equation (7). Thus, it can be assumed that the drag force only has a component in the x-direction. As the velocity $\overrightarrow{u}=u_{\mathrm{x}}\left(y\right)·{\overrightarrow{\mathrm{e}}}_{\mathrm{x}}$ is assumed to be laminar and constant in the main branch, it has a parabolic profile, which is provided by

$$u_{\mathrm{x}}\left(y\right)=u_{\max }\left[1-{\left({\displaystyle \frac{y}{R_{\mathrm{v}}}}\right)}^2\right],$$

(21)

with $u_{\max }=2u_{\mathrm{mean}}$. y corresponds to the y-coordinate, and $R_{\mathrm{v}}$ to the radius of the vessel.

(2)

Evaluation of the magnetic force $F_{\mathrm{m}}$:

The magnetic force pulls the particles towards the magnet. Thus, in this paper, $F_{\mathrm{m}}$ is assumed to have only a component in the y-direction as a consequence of the fitted exponential decay of $\left|B\right|$ in y-direction. This force is calculated according to Equation (14). Therefore, the particle’s magnetization M has to be known. It is determined by once reading out $F_{\mathrm{m},\mathrm{y}}$ and the gradient grad(B) at a fixed location in the COMSOL simulation. Then, M is derived using Equation (14). For reasons of simplification, M is assumed to be constant. In the next step, $F_{\mathrm{m}}$ is calculated using the fitted magnetic gradient for all $\vartheta$.

4. Results

4.1. Evaluation of the Magnetic Flux Density

In Figure 6, the isolines of the magnetic flux density for varying magnetization angles $\vartheta$ between the single magnets are depicted. For comparison, the color map is the same for all subfigures. It becomes immediately observable that the Halbach array independent of the applied $\vartheta$ has a strong and a weak side. In the case of Figure 6, the strong side is facing the vessel. Furthermore, the field distribution of B is symmetric for all $\vartheta$. This also holds true for an odd number of magnets. The maximum of the magnetic flux density B is always located between the first two and between the last two magnets, respectively. It can also be seen that $B_{\max }$ decreases with increasing angle $\vartheta$. Moreover, by comparing the shape of the isolines in Figure 6, it is visible that the magnetic field distribution is more compact for higher $\vartheta$. This can be seen particularly clearly in

Table 3. Summary of the maximum value of the magnetic flux density $B_{\max }$, and the distances $d_{\max }$, and $d_{\mathrm{mid}}$ between the bottom edge of the magnetic array to the isoline with $B=0.02$ T. $B_{\max }$ is always located between the first two magnets, $d_{\max }$ is the maximum distance, and $d_{\mathrm{mid}}$ was evaluated in the center of the array, so between the fourth and fifth magnets.

$\mathit{\vartheta }$	30$°$	45$°$	60$°$	65$°$	70$°$	75$°$	80$°$	85$°$	90$°$	95$°$	100$°$	105$°$	110$°$	115$°$	120$°$	135$°$
$B_{\max }$ in T	1.02	0.97	0.92	0.91	0.89	0.88	0.87	0.86	0.85	0.85	0.84	0.84	0.83	0.81	0.80	0.93
$d_{\mathrm{mid}}$ in cm	10.91	9.44	7.30	6.67	6.07	5.54	5.07	4.67	4.36	4.11	3.90	3.76	3.65	3.59	3.56	3.33
$d_{\max }$ in cm	10.91	9.44	7.30	6.67	6.18	6.04	5.97	5.81	5.64	5.43	5.14	4.87	4.64	4.54	4.47	3.87

. It shows that, for $\vartheta ={30}^{\circ }$, the distance between the array and the outer isoline $B=0.02$ T is still 10.91 cm. With increasing $\vartheta$, the maximum distance $d_{\max }$ decreases to $d_{\max }=3.87$ cm for $\vartheta ={135}^{\circ }$ and is no longer in the center of the array. Moreover, the centered distance $d_{\mathrm{mid}}$ decreases too. As a result, the field is much more homogeneous in the x-direction for higher $\vartheta$. Figure 7 shows the magnetic flux density for $\vartheta ={180}^{\circ }$ with the magnetization direction in the $\pm x$-direction and $\pm y$-direction. It reveals that the shape of the magnetic flux density for the two simulations looks quite similar. Here, the field is even more compact with $d_{\mathrm{mid}}=2.32$ cm, and the array no longer has a strong and a weak side.

The aforementioned compactness of the magnetic field distribution also influences the range of the magnetic flux density. This can be seen for chosen magnetization angles $\vartheta$ in Figure 8 on the left side. Overall, the flux density B is strongest for both the strong and the weak side of the array facing the vessel for large $\vartheta$ at the surface of the magnets but decays faster.

4.2. Evaluation of the Magnetic Gradient

In COMSOL Multiphysics^® 6.1, the magnetic flux density is calculated numerically on a discrete mesh. This can be seen clearly in Figure 8 on the left side from the hard edges in the function of the strong side. Thus, it is hardly possible to derive the gradient of the magnetic flux density directly because many irregularities and noise would be present. Therefore, in this paper, the gradient is determined by first fitting B to the exponential function in Equation (15) and analytically deriving the fitted function:

$$\nabla B\widehat{=}{\partial }_{y}B=-k_1{k}_2{\mathrm{e}}^{-k_{2}y}$$

(22)

Figure 8 on the right side depicts exemplary results for the determined gradient. To see the trends, a small (${45}^{\circ }$), middle (${90}^{\circ }$), and big (${135}^{\circ }$) magnetization angle are depicted and analyzed in more detail in the following.

Overall, the gradient in Figure 8 shows a clear difference between the weak and strong sides of the array. The gradient is stronger for the weak side of the array directly at the surface of the magnets but decays very fast with distance. For the strong side, it is vice versa. Thus, at a distance of approximately 0.5 cm, the gradient of the strong side is higher than the one of the weak side. When comparing grad(B) for the different $\vartheta$ in Figure 8, it is noticeable that, for both the strong and the weak sides, the gradient is larger for higher angles. For the strong side, the curves converge again at a distance of about 1.5 cm from the magnets’ surface, while the curves of the weak side still do not intersect at a distance of 2 cm.

4.3. Evaluation of the Fitting of the Magnetic Flux Density

As mentioned in the previous section, the gradient of the magnetic flux density was determined by fitting the magnetic flux density to the exponential approach in Equation (15). The mean and the maximum fitting error ${\epsilon }_{\mathrm{mean}}$ and ${\epsilon }_{\max }$ for all $\vartheta$ as absolute and relative errors are summarized in

Table 4. Overview of the fitting results for the different magnetization angles of the Halbach array. The absolute fitting errors $\epsilon$ are provided in mT. In brackets, the relative errors are listed in %. In the last two columns, the fitting error was only calculated inside the vessel for a distance of 0.5 cm between the magnetic array and the vessel as the magnetic force on the SPIONs is only evaluated in this range.

Magnetization Angle	Side	${\mathit{\epsilon }}_{\mathbf{mean}}$	${\mathit{\epsilon }}_{\mathbf{max}}$	${\mathit{\epsilon }}_{\mathbf{mean}\mathbf{in}\mathbf{vessel}}$	${\mathit{\epsilon }}_{\mathbf{max}\mathbf{in}\mathbf{vessel}}$
30$°$	strong	24.8 (9.4%)	120.5 (22.2%)	18.5 (7.3%)	35.7 (12.5%)
	weak	5.0 (34.5%)	17.3 (95.8%)	2.5 (23.8%)	7.6 (72.3%)
45$°$	strong	21.5 (6.6%)	118.0 (17.3%)	16.1 (5.0%)	30.6 (8.5%)
	weak	3.4 (27.5%)	20.8 (122.3%)	2.5 (26.2%)	4.0 (122.3%)
60$°$	strong	18.7 (5.4%)	80.3 (13.7%)	13.7 (3.9%)	26.5 (6.4%)
	weak	2.1 (7.6%)	7.8 (25.8%)	1.6 (5.1%)	2.9 (12.5%)
65$°$	strong	18.2 (5.3%)	80.5 (13.5%)	13.2 (3.7%)	25.9 (6.3%)
	weak	3.3 (12.8%)	15.5 (49.6%)	2.1 (6.9%)	4.3 (21.3%)
70$°$	strong	17.8 (5.2%)	86.4 (13.5%)	12.7 (3.6%)	26.0 (6.0%)
	weak	4.4 (15.2%)	28.2 (57.3%)	2.7 (8.1%)	5.5 (26.4%)
75$°$	strong	17.5 (5.2%)	88.5 (13.5%)	12.2 (3.5%)	26.1 (6.0%)
	weak	5.2 (16.2%)	35.8 (59.2%)	3.3 (8.9%)	6.5 (28.5%)
80$°$	strong	16.9 (5.2%)	91.3 (13.5%)	11.7 (3.5%)	26.0 (5.9%)
	weak	6.9 (18.3%)	25.5 (61.4%)	4.5 (11.2%)	7.7 (31.9%)
85$°$	strong	16.2 (5.2%)	86.4 (13.4%)	11.1 (3.4%)	25.0 (5.8%)
	weak	5.9 (16.1%)	31.2 (56.6%)	4.0 (9.6%)	6.8 (28.1%)
90$°$	strong	15.2 (5.1%)	81.4 (13.2%)	10.3 (3.3%)	23.7 (6.0%)
	weak	5.9 (15.1%)	25.0 (51.5%)	4.0 (9.4%)	6.5 (25.9%)
95$°$	strong	14.0 (4.9%)	81.7 (12.9%)	9.5 (3.1%)	21.6 (6.0%)
	weak	5.4 (13.3%)	24.7 (43.8%)	3.7 (8.6%)	6.5 (22.6%)
100$°$	strong	12.8 (4.7%)	81.7 (12.6%)	8.6 (2.9%)	19.4 (5.9%)
	weak	3.5 (8.9%)	10.4 (29.7%)	2.4 (5.4%)	5.1 (15.1%)
105$°$	strong	11.6 (4.5%)	75.0 (12.3%)	7.7 (2.8%)	17.2 (5.7%)
	weak	4.3 (9.1%)	19.8 (26.5%)	2.9 (6.2%)	5.9 (15.0%)
110$°$	strong	10.7 (4.3%)	65.2 (12.1%)	7.0 (2.6%)	15.3 (5.6%)
	weak	3.7 (7.3%)	16.0 (22.1%)	2.4 (4.8%)	5.1 (11.7%)
115$°$	strong	10.1 (4.2%)	52.7 (12.1%)	6.5 (2.5%)	14.2 (5.6%)
	weak	3.1 (6.0%)	10.6 (0.4%)	2.4 (3.7%)	4.2 (9.3%)
120$°$	strong	9.7 (4.3%)	40.7 (12.4%)	6.2 (2.5%)	13.7 (5.6%)
	weak	2.7 (4.8%)	8.7 (17.5%)	1.9 (2.9%)	3.4 (7.3%)
135$°$	strong	1.6 (5.1%)	53.6 (14.7%)	6.3 (2.9%)	15.1 (7.4%)
	weak	3.1 (4.4%)	34.8 (12.7%)	2.3 (3.4%)	5.0 (8.1%)
mean value	strong	15.4 (5.3%)	80.2 (13.1%)	10.7 (3.5%)	22.6 (6.6%)
	weak	4.2 (13.6%)	20.7 (45.8%)	2.8 (9.0%)	5.4 (28.6%)

for the whole evaluation domain as well as only limited to a distance inside the vessel. Furthermore, the fitting parameters $k_1$ plus $k_2$ and the relative fitting error are depicted in Figure 9.

Overall, the fitting performance for the strong side performs well with ${\epsilon }_{\mathrm{mean}}=5.3$% and ${\epsilon }_{\max }=13.1$% over the whole evaluation domain. Moreover, as the right plot in Figure 9 shows, the error is relatively constant for the strong side over the magnetization angle $\vartheta$. However, the fitting results for the weak side are, with an error of ${\epsilon }_{\mathrm{mean}}=13.6.3$% and ${\epsilon }_{\max }=45.8$%, inadequate. Especially for the small angles $\vartheta ={30}^{\circ }$ and $\vartheta ={45}^{\circ }$ as well as for angles in the range $\vartheta \in \left[{65}^{\circ },{95}^{\circ }\right]$, the relative fitting errors are insufficiently high.

Next to the relative fitting error, the left plot in Figure 9 also shows the values for fitting parameters $k_1$ and $k_2$. Overall, the amplitude parameter $k_1$ is approximately a factor of 100 higher for the strong side compared to the weak side. For the strong side, $k_1$ has its maximum at $\vartheta ={60}^{\circ }$ with $k_1=17.5$ mT and decays with increasing $\vartheta$. On the other hand, for the weak side, $k_1$ increases with increasing $\vartheta$. The exponent parameter $k_2$, in contrast, is greater for the weak side of the array. Here, $k_2$ increases with increasing $\vartheta$ for the strong side and decreases with increasing $\vartheta$ for the weak side of the array, respectively.

4.4. Evaluation of the Parameter Study on the SPION Distribution

The main evaluation parameter in particle steering is the distribution of the SPIONs to the upper and lower branches as well as the percentage of the particles trapped by the magnets. This section presents the results for the different influence parameters on the SPION distribution in detail, namely the influence of the weak vs. strong side of the array, the quantity of magnets (Figure 10), the magnetization angle $\vartheta$ between the single magnets (Figure 11), and the distance between the magnetic array and vessel (Figure 12). In addition, the results for the maximum gradient grad(B) and the total applied magnetic energy $W_{\mathrm{m}}$ inside the whole vessel are opposed.

4.4.1. Influence of Weak vs. Strong Side of the Magnet Array

A Halbach array has a strong and a weak side, resulting out of the constructive and destructive overlapping of the magnetic fields of the single magnets. Thus, fewer particles should be trapped if the weak side is applied to the vessel. This can especially be seen in Figure 11 and Figure 12, where the number of trapped particles is much higher for the strong side (blue lines). The number of SPIONs taking the lower branch is lower for the strong side.

4.4.2. Influence by Quantity of Magnets

A Halbach array consisting of more magnets means a longer array, and, therefore, the time when a magnetic force is applied to the particles while propagating underneath the array is longer. This can be seen in Figure 10a on the left side, where the SPION distribution for the strong side of the array over the number of magnets is depicted. The more magnets are used, the more SPIONs become trapped, and the less SPIONs take the lower branch.

In Figure 10b, the corresponding maximum gradient in the vessel and the magnetic energy calculated according to Equation (19) are plotted. The maximum gradient is always located below the middle magnet (odd number of magnets) or the middle two magnets (even number of magnets). The influence on the maximum gradient in the vessel is negligibly low (see axis scale in Figure 10). On the other hand, the magnetic energy in the vessel increases linearly with longer magnetic arrays.

4.4.3. Influence of Magnetization Angle between Magnets

To the best of the authors’ knowledge, in the literature, the influence of the magnetization angle $\vartheta$ has not been investigated for MDT so far. However, Figure 11a reveals that the magnetization angle $\vartheta$ between the magnets also plays an important role in influencing the particle movement. Overall, a higher magnetic force corresponds to a higher number of trapped particles and a smaller number of SPIONs taking the lower branch. In Figure 11, it can be seen that, for the strong side of the array, higher angles cause more trapped particles till a magnetization angle of ${120}^{\circ }$. The maximum number of trapped particles is approximately 1/3 of all particles (36 of 100). Interestingly, higher angles than ${120}^{\circ }$ reduce the number of trapped particles again. The minimum number of SPIONs taking the lower branch is 30% for $\vartheta \in \left[{80}^{\circ },{110}^{\circ }\right]$. For the weak side of the array, with increasing $\vartheta$, the number of trapped SPIONs increases approximately linearly, and the number of SPIONs taking the lower branch decreases approximately linearly, respectively.

Figure 11b shows the trend of the maximum gradient and the magnetic energy in the vessel. The maximum gradient starts at 32.4 T/m at ${30}^{\circ }$ and increases nearly linearly with the magnetization angle to 47.1 T/m at ${120}^{\circ }$. The magnetic energy in the vessel is low for small and high $\vartheta$. It has its maximum between $\vartheta ={60}^{\circ }$ and $\vartheta ={90}^{\circ }$.

The results for higher magnetization angles are listed in

Table 5. Comparison of the two magnetization patterns with a magnetization angle $\vartheta ={180}^{\circ }$. In the magnetization pattern, the single magnets of the array can be magnetized in the $\pm y$-direction or the $\pm x$-direction. The maximum gradient and 2D magnetic energy are evaluated inside the vessel.

Magnetization	Trapped	Upper Branch	Lower Branch	Maximum Gradient	2D Magnetic Energy
Pattern	in %	in %	in %	in T/m	in J/m
Down–up	27	34	39	45.761	20.983
Right–left	27	34	39	47.937	20.983

. It compares the particle distribution plus the maximum gradient, and magnetic energy for the two simulations with a magnetization angle $\vartheta ={180}^{\circ }$. The results reveal no influence in the SPION distribution either in the maximum gradient or magnetic energy in the vessel for arranging the magnetization direction in the $\pm x$-direction or $\pm y$-direction.

4.4.4. Influence of Distance between Magnetic Array and Vessel

Figure 12 reveals that varying the distance between the array and the vessel also influences the simulation results. Overall, a higher distance between magnets and vessel corresponds to a weaker magnetic force on the SPIONs. Thus, with increasing distance, the number of trapped SPIONs reduces for both the strong and the weak side of the array. For the strong side, the number of SPIONs taking the desired upper branch increases approximately linearly with the distance and obtains 52% at a distance of 2 cm even over the 50% distribution without an applied magnet. For the weak side of the array and at a distance of 1 cm, the SPIONs spread equally to the upper and lower branch, and no particle is trapped by the magnets anymore. Thus, at a distance of 1 cm and greater, the magnetic force does not influence the chosen parameter setup. Figure 12b shows that, with increasing distance, both the maximum gradient and the magnetic energy in the vessel decrease.

4.5. Evaluation of the Magnetic and Hydrodynamic Drag Force

As the movement of the SPIONs can be predicted when knowing the forces acting on the SPIONs, in this section, the two main forces, namely the magnetic and (hydrodynamic) drag force ($F_{\mathrm{m}}$ and $F_{\mathrm{d}}$), are investigated in more detail. Both forces are calculated for the strong and weak sides of the array underneath the fourth magnet, as illustrated in Figure 5. The resulting magnitudes of the forces are depicted in Figure 13.

4.5.1. Evaluation of $F_{\mathrm{d}}$

As aforementioned, the drag force is calculated based on Equations (7) and (21). Thus, due to the parabolic velocity profile, $F_{\mathrm{d}}$ has its maximum at the center line of the vessel with $|F_{\mathrm{d},\max }|\approx 9.4·{10}^{-11}$ N, and a minimum magnitude of 0 N at the vessel wall. This can also be seen in Figure 13. Furthermore, as the drag force is the force washing the particles within the background velocity, it only has a component in the x-direction.

4.5.2. Evaluation of $F_{\mathrm{m}}$

On the contrary, the magnetic force $F_{\mathrm{m}}$ pulls the SPIONs towards the magnet. Thus, it mainly has a component in the y-direction. Before the magnitude of the magnetic force can be calculated according to Equation (14), the magnetization M of one SPION has to be determined. For this purpose, $F_{\mathrm{m},\mathrm{y}}=8.38·{10}^{-13}$ N and grad$\left(B\right)=32$ T/m are extracted from a fixed position in the simulation. Then, $M\approx 4·{10}^5$ A/m is derived using Equation (14). This value is in good accordance with the magnetization values of SPIONs with regard to [5,96,97]. By applying $M=4·{10}^5$ A/m to Equation (14), the magnetic force reaches values from $2·{10}^{-13}$ N to $2.3·{10}^{-12}$ N inside the vessel for all magnetization angles $\vartheta$ for the strong side of the array (compare Figure 13). For the weak side of the array, the magnetic force was in the range of $4.9·{10}^{-15}$ N to $3.2·{10}^{-12}$ N for all $\vartheta$.

Figure 13 shows that the magnitude of $F_{\mathrm{m}}$ decreases with the distance to the magnets. Thus, for both sides of the array, the strongest force is at the vessel wall closer to the magnets and the weakest at the vessel wall further away from the magnets, respectively. For the strong side, $F_{\mathrm{m}}$ is much higher for greater $\vartheta$, whereas, for the weak side, the influence of $\vartheta$ is much smaller. This can also be seen in Figure 14a, where the mean and maximum values of $F_{\mathrm{m}}$ for the weak and the strong side for all $\vartheta$ are depicted. In Figure 14b, the ratio between the strong and the weak side is shown. By comparing the magnitude of $F_{\mathrm{m}}$ for the strong and the weak sides, it reveals that $F_{\mathrm{m}}$ is for all $\vartheta$ approximately double (factor of ≈1.85 for $\vartheta \in [{60}^{\circ },{100}^{\circ }]$) as high for the weak side compared to the strong side. For $\vartheta ={135}^{\circ }$, the difference is only a factor of approximately 1.5. However, as Figure 13 shows, the magnetic force decays much faster for the weak side, and at a distance of 2 cm to the magnets, $|F_{\mathrm{m}}|$ is approximately a factor of 100 stronger for the strong side. In addition, Figure 14 reveals that the maximum magnetic force inside the vessel for the weak side is for small $\vartheta$ higher than for the strong side. However, this difference becomes smaller with increasing $\vartheta$, and, for $\vartheta >{65}^{\circ }$, $F_{\mathrm{m},\max }$ is greater for the strong side again. Nevertheless, as shown on the right plot in Figure 14, the ratio between strong and weak side for the maximum value of the magnetic force is relatively small for all $\vartheta$.

4.5.3. Comparison of $F_{\mathrm{m}}$ and $F_{\mathrm{d}}$

Since the ratio of magnetic to drag force is particularly relevant for the direction of motion of the SPIONs, this ratio is illustrated in Figure 15. The ratio $F_{\mathrm{d}-\mathrm{m}}$ is equal to zero at the vessel walls as $F_{\mathrm{d}}=0$ N there. The maximum values for the strong and the weak side of the array are located at a distance between 1 cm and 1.3 cm from the magnets.

Overall, the ratio for the weak side is much higher than for the strong side as the magnetic force is weaker for the weak side. This can also be seen in Figure 14, where the mean and maximum value of $F_{\mathrm{d}-\mathrm{m}}$ for both sides of the array are depicted on the left plot, and the ratio between these two sides on the right plot. Moreover, Figure 15 reveals the influence of the magnetization angle $\vartheta$. Whereas the influence is negligible for the strong side, for the weak side and $\vartheta ={30}^{\circ }$, $F_{\mathrm{d}-\mathrm{m},\max }$ is with $\approx 1279$ very high and decreases with increasing $\vartheta$. For $\vartheta ={135}^{\circ }$, $F_{\mathrm{d}-\mathrm{m},\max }$ is only 274. The same trend can be seen in the ratio on the right plot in Figure 14 too. For $\vartheta ={30}^{\circ }$, the ratio $F_{\mathrm{d}-\mathrm{m},\max }$ is only approximately 0.19 and increases to $F_{\mathrm{d}-\mathrm{m},\max }=0.55$ for $\vartheta ={135}^{\circ }$.

5. Discussion and Limitations of This Study

This paper conducted a systematic parameter study with the focus on varying the magnetization angle for steering SPIONs in a typical MDT scenario using a linear Halbach array as a magnetic field source. In this section, the results are discussed in detail in order to build an optimized array for particle steering. In addition, the content is compared with the state-of-the-art approaches in MDT and the limitations of this work are identified.

5.1. Discussion of the Magnetic Flux Density and Its Effect on Particle Steering

In the ideal case, a Halbach array has a one-sided magnetic field caused by constructive and destructive inference on each side of the array [76,77,88]. However, as can be seen in Figure 6, this superposition effect is not ideal due to the air space between the magnets and their round shape. As shown in Section 4.1, the magnetic flux density is more compact for a higher magnetization angle $\vartheta$ (between $\vartheta ={30}^{\circ }$ and $\vartheta ={135}^{\circ }$, there is approximately a factor of 2.8 for $d_{\max }$ and 3.3 for $d_{\mathrm{mid}}$, respectively). Thus, this paper assumes that the magnetic flux density is constant in the x-direction. This assumption has already been proven and discussed in [43,51] for a Halbach array with a magnetization angle of $\vartheta ={90}^{\circ }$, where the resulting standard deviation was approximately a factor of 50 smaller than the mean magnetic flux density. Sharma et al. [98] investigated the magnitude of the magnetic force in the x- and y-directions for one single permanent magnet. They also concluded that the force in the y-direction is dominant. On the other hand, Bernad et al. [99] showed that the difference in $F_{\mathrm{m},\mathrm{x}}$ and $F_{\mathrm{m},\mathrm{y}}$ depends on the shape of the magnet. However, the SPIONs accelerate before the magnet and decelerate after the magnet due to the magnetic force. When the magnetic field source is symmetric, this effect is approximately canceled out. Nevertheless, as can be seen in Figure 13, Figure 14 and Figure 15, the magnitude of the drag force $|F_{\mathrm{d}}|$ is much higher, and, therefore, it can be assumed that the impact of the magnetic force in the x-direction can be neglected.

Figure 10 shows that the quantity of magnets in the array, and thus the length of the array, significantly influences the particle distribution. A Halbach array with a higher number of magnets means a larger area where the magnetic field can influence the particles in the vessel. Thus, the magnetic force has more time to act on the particles in the case of a more extended array. Also, the magnetic field itself becomes larger with more magnets. Consequently, more particles are attracted to the wall of the vessel where the array is placed. In our recently published letter [41], we showed that the sorting performance depends primarily on the total magnetic effort in terms of the applied magnetic energy. Here, it was shown that a more extended array with a lower magnetic field strength has the same impact on particle propagation as a shorter array with a higher field strength measured in terms of investigated magnetic energy. Figure 10b in this paper shows that the magnetic energy increases with increasing number of magnets too. On the other hand, the maximum gradient in the vessel is relatively constant because more magnets do not really influence the change in the magnetic field. Therefore, the magnetic field of the added magnets has a very low influence on the high magnetic field area between the magnets.

In [51], it was stated that, by decreasing the distance $d_{\mathrm{m}}$ between the single magnets, the magnetic flux density becomes more homogeneous, and thus, the superposition effect becomes stronger. However, because of the torque, the distance between the magnets is crucial for mechanical stability [88]. In general, mechanical stability is an essential factor that must be considered for the practical implementation of the presented steering approach. For example, the two simulations with $\vartheta ={180}^{\circ }$ show equal results regarding the isolines of the magnetic flux density (compare Figure 7), and the SPION distribution, the maximum gradient, and the energy in the vessel were equal (see Table 5). However, of course, the setup with magnetization directions down and up is mechanically stable, whereas the magnets in the array with the magnetization directions right and left repel each other.

In addition to the strength of the magnetic force, its effective range is another crucial factor. In general, there is a trade-off between the strength of the magnetic force and its effective range. This becomes particularly obvious when comparing the strong and weak sides of the array: Figure 8b shows that the gradient of B (directly proportional to the magnetic force) for the weak side of the array is stronger close to the array than for the strong side but decays much faster. A similar effect can also be seen in Figure 13. Here again, the magnetic force close to the array is stronger for the weak side than for the strong side. We can, therefore, conclude that it would be desirable to have a constant gradient over the vessel’s diameter. If the gradient closer to the vessel wall increases, this expedites the trapping of particles. Obviously, the desired range depends strongly on the particular application. Thus, Ijiri et al. [52] separated different sizes of the magnetic nanoparticles in a tube close to the surface of a 47-magnet Halbach array, which was arranged so that the weak side faces the tube. However, especially in MDT, a high range of magnetic force is essential to steer the particles in deeper body regions. Therefore, in most cases, the strong side of the Halbach array is used in MDT [46,47,49,100]. The magnetization angle $\vartheta$ between the magnets influences the shape of the magnetic field. Figure 6 shows that the field underneath the vessel becomes compressed with higher magnetization angles. Therefore, the magnetic field lines in the vessel area are closer to each other, which means that the field is stronger and can influence the particles more strongly. Thus, more particles become attracted to the upper vessel wall. This also affects the effective range of the magnetic flux density and its gradient. For larger $\vartheta$, the gradient of B (and thus the magnetic force) directly at the magnets is greater but then also decays more strongly. However, in our considered distance range, the gradient of B for large $\vartheta$ is always higher than that of smaller angles, so larger $\vartheta$ are preferable for determining a higher range. This is in good accordance with the results for the particle distribution in Figure 11 and with the mean and maximum magnitude of $|F_{\mathrm{m}}|$ in Figure 14. Here, the mean value of the magnetic force was the strongest for $\vartheta \in [{80}^{\circ },{100}^{\circ }]$.

The same trend can be seen for the weak side of the array. As suggested in our previous work, the strong and weak sides can be switched by mechanical rotation of all single magnets [43,51] or by changing the current direction of EMs [46], to wash the trapped particles off. This requires the weakest possible force (and therefore the weakest possible gradient) on the weak side. Figure 8 shows that this is the case for small magnetization angles. Thus, switching between the strong and weak side results in a trade-off between a field with a long effective range (large $\vartheta$) on the strong side and a field with a short effective range (small $\vartheta$) for the weak side. By evaluating the ratio between the strong and the weak side, as shown in Figure 14b, the best trade-off is approximately in the range of $\vartheta \in [{60}^{\circ },{100}^{\circ }]$. Since, as aforementioned, the simple production of the array is also an important design criterion, $\vartheta ={90}^{\circ }$ seems to be a good choice. Furthermore, for easy switching between the strong and the weak side, the use of EMs as in [46] is recommended.

5.2. Discussion of the Fitting Results

As the magnetic force is directly proportional to the gradient of the magnetic flux density, it is worth investigating the gradient over the vessel cross-section in detail. For determining the gradient, the magnitude of B along the evaluation line (depicted in Figure 5) was fitted to the exponential approach in Equation (15) and afterwards derived analytically as shown in Equation (22) for the strong and the weak side of the array and various magnetization angles $\vartheta$. As Figure 9 and Table 4 depict, the fitting results for the strong side were quite good (mean fitting error of ${\epsilon }_{\mathrm{mean}}=5.3$%) but insufficient for the weak side, with a mean error of ${\epsilon }_{\mathrm{mean}}=13.6$% and individual maximum errors of up to ${\epsilon }_{\max }=122.3$%. Especially for the weak side and small $\vartheta$, the fitting results were poor. The reason for this can be found in the shape of B for these cases. For example, Figure 8a shows that B of the weak side for $\vartheta ={45}^{\circ }$ at a distance of approximately 1.8 cm rises slightly again. This, of course, contradicts the assumed strictly monotonic decrease in the exponential fitting function. Since the field strength values are also very small (${\epsilon }_{\mathrm{mean}}=4.2$ mT), this results in significant relative errors. Nevertheless, as the magnitude of the magnetic field strength and the magnetic gradient are small, the magnitude of the magnetic force is small as well. Thus, the drag force dominates the magnetic force anyway. For this reason, despite the high relative errors, the calculations in this paper continue with the fitted values for the consideration of the forces. For a more accurate fitting result, as in [51], fitting functions with more degrees of freedom or analytical approaches, such as in [56,57,101], must be used.

5.2.1. Discussion of Fitting Parameter $k_1$

Fitting parameter $k_1$ corresponds to the magnetic flux density’s amplitude factor. As shown in Figure 9, $k_1$ is approximately a factor of 100 stronger for the strong side than for the weak side. By inserting the fitted magnetic flux density of Equation (22) in Equation (14), the magnetic force can be expressed as

$$F_{\mathrm{m},\mathrm{y}}=V_{\mathrm{p}}M|\nabla B|=-V_{\mathrm{p}}M{k}_1{k}_2{\mathrm{e}}^{-k_{2}y}.$$

(23)

It is obtained that the magnitude of the magnetic force is directly proportional to the factor $k_1$. Therefore, the greater $k_1$, the greater the magnitude of the magnetic force $|F_{\mathrm{m}}|$. However, as Equation (23) depicts, the exponential parameter $k_2$ also contributes to $|F_{\mathrm{m}}|$. For this reason, the results that a greater $k_1$ corresponds to a greater $|F_{\mathrm{m}}|$ cannot be seen in Figure 13, Figure 14 and Figure 15.

5.2.2. Discussion of Fitting Parameter $k_2$

As aforementioned, $k_2$ corresponds to the exponential fitting parameter. As Figure 9 shows, $k_2$ is much greater for the weak side of the array than for the strong side. By determining the gradient of B, this parameter is both in the exponent and directly proportional to the amplitude of the magnetic force $|F_{\mathrm{m}}|$ (compare Equation (23)). This means that $k_2$ influences both the range of the effective magnetic force and its magnitude. For a large magnitude of $|F_{\mathrm{m}}|$, $k_2$ should be large, whereas, for an extensive range, $k_2$ should be small. This results in a trade-off, which depends on the requested application. As a more extensive range is required for MDT, $k_2$ should tend to have lower values here. This can also be seen in the simulation results. For the strong side of the array, the maximum force $|F_{\mathrm{m}}|$ is much lower due to the small $k_2$ (Figure 9 and Figure 13), but the effective range in the y-direction is much greater in this case so that more SPIONs are deflected (Figure 11).

5.3. Discussion of the Maximum Gradient and Magnetic Energy for Predicting Particle Steering

Since the analysis of the particle distribution is very complex, this section will examine whether the evaluation of the maximum gradient in the vessel or the total applied 2D magnetic energy in the vessel can also be used to determine the particle distribution. The variation in the quantity of magnets in Figure 10 shows that the magnetic energy increases linearly with the quantity of magnets. A Halbach array with a higher number of magnets means a longer time when the magnetic field can influence the particles in the vessel. Also, the magnetic field itself becomes greater with more magnets. Consequently, more particles will be attracted to the wall of the vessel where the array is placed. This trend can be seen for the trapped SPIONs in Figure 10 and follows the same trend as the magnetic energy in the vessel. On the other hand, the maximum gradient in the vessel is relatively constant because more magnets do not influence the magnetic field change in y-direction that much. The highest gradient is between the single magnets. The magnetic field of the added magnets has a very low influence on the high magnetic field area between the magnets. The change in the total magnetic energy in the vessel is much more significant in this case because a higher number of magnets means a larger area with an applied magnetic field. Summarized, these results show that the maximum gradient is unsuitable for predicting the particle distribution. In contrast, Figure 12 shows the particle distribution over the distance between the vessel and the magnetic array. Here, it can be observed that the trend of the magnetic energy and the maximum gradient corresponds well with that of the trapped particles.

The magnetization angle $\vartheta$ between the magnets influences the shape of the magnetic field (see Figure 6). As Table 3 shows, with higher angles, the field underneath the vessel becomes compressed. Therefore, the magnetic field lines in the vessel area are closer, which means that the field is stronger and can influence the SPIONs more, resulting in more trapped particles. The corresponding SPION distribution is depicted in Figure 11. For $\vartheta \in [{80}^{\circ },{110}^{\circ }]$, the fewest particles take the lower output. This means the magnetic force is strongest here, and the effective range is highest. On the other hand, most SPIONs are trapped for $\vartheta >{95}^{\circ }$. This means that the magnetic force is strongest near the magnets. These effects can also be recognized in Figure 13. Nevertheless, the magnetic energy is highest in the range of 60–80$°$ and decreases for higher angles. The maximum gradient, on the other hand, increases approximately linearly. However, neither trend matches the SPION distribution.

Overall, the trends of the maximum gradient in the vessel and the magnetic energy in the vessel deviate from the results of the SPION distribution. Therefore, neither parameter is suitable for predicting the particle distribution. Thus, it is better to consider the total magnetic force. In this paper, the magnetic force was only analyzed at its strongest position. However, we showed in our previous publication [51] that the magnitude of the magnetic force depends strongly on the evaluated position, especially for the strong side of the array. Therefore, we recommend evaluating the particle distribution, although the simulation effort is much higher.

5.4. Comparison with the State-of-the-Art Research

In the following, the investigations of this paper are compared with the state-of-the-art research for steering magnetic nanoparticles in an MDT scenario. An overview of the corresponding studies with their field sources, evaluation parameters, and a summary of the studies’ aims is provided in

Table 6. Overview of state-of-the-art publications in the context of steering magnetic nanoparticles for MDT. The study type (simulation or measurement), study topology, and evaluation parameters regarding analyzing the steering performance are provided. Y-vessel corresponds to a Y-shaped vessel. Furthermore, the aims of the studies or the optimization parameters are summarized.

Authors	Year	Study Type		Study Topology	Magnetic Field Source	Evaluation Parameters	Aim of Study, Optimization Parameter
		Sim.	Meas.
Abolfathi et al. [102]	2020	x	x	Sym. Y-vessel	4 EMs on opposite sides	Dispersion changes in particles	Aggregating particles and reducing dispersion of particles
Cai et al. [103]	2020	x	x	Mice vessel	One EM	Particle attraction (capture rate), fluorescent imaging of mice w. particles	MDT efficiency in mice, optimizing capture rate
Hoshiaret et al. [30]	2020	x		Sym. vessel with mult. branches	3 circular arranged EMs	Particle distribution	Algorithm for steering particles through mult. bifurcations at different velocities
Kee et al. [104]	2020	x	x	Straight vessel	3D Halbach array with 9 PMs, $\vartheta ={45}^{\circ }$	$F_{\mathrm{m}}$, number of trapped particles, images of gray scale intensity	Trapping particles using 3D Halbach array
Park et al. [31]	2020	x		Sym. Y-vessel, brain vessel	6 EMs on opposite sides	B, grad(B), particle distribution, target and sticking ratio, velocity profile	Algorithm for steering particles using a haptic device
Shiriny et al. [54]	2020	x		Sym. vessel with 3 outlets	Halbach array with 3 PMs, $\vartheta ={90}^{\circ }$	H, B, and $F_{\mathrm{m}}$ along tube center line, particle distribution	Optimizing magnetophoretic separation, parameter study
Le et al. [32]	2021	x	x	Sym. Y-vessel, 3D model vessel in human brain	4 EMs on opposite sides	H, grad(H), $F_{\mathrm{m}}$, particle distribution and trajectory, images of vessel	Trapping and steering particles by shifting focal point
Nguyen et al. [17]	2021	(x) *	x	Sym. Y-vessel	9 EMs focused on region of interest	B, grad(B), particle distribution, trapping rate	Design of magn. field for particle trapping, parameter study
Sarraf et al. [105]	2021	x		Sym. Y-vessel mult. branches, unsym. vessel mult. branches	6 diff. array with each 5 PMs: linear and 3D Halbach arrays with $\vartheta \in \{0^{\circ },{90}^{\circ },{180}^{\circ }\}$	B, particle distribution	Comparing particle tracing through healthy and tumorous vessel structures
Stevens et al. [47]	2021	x	x	Straight tube	Lin. Halbach arrays (diff.  number of magnets), $\vartheta ={90}^{\circ }$	B, grad(B), magn. recovery of particles, number of trapped particles via	Immunomagnetic enrichment process, magn. recovery of cells
Bernad et al. [99]	2022	(x) *	x	Unsym. vessel	Different single PMs (varying size, material)	$B,H$ at fixed pos. (all components), $F_{\mathrm{m}}$, camera images of tube cross-section, particles deposition and target efficiency, particle deposition length and thickness	Investigation of diff. PMs and PM positions for particle steering in MDT
Chakrabarty et al. [25]	2022	x		Sym. Y-vessel	2 EMs on opposite side of vessel	H, grad(H), trajectory of one particle based on analytical calc. of $F_{\mathrm{m}}$ and $F_{\mathrm{d}}$	Steering particles, avoiding trapping
Chakrabarty et al. [16]	2022	x	x	Sym. Y-vessel	2 double coils on opposite sides of vessel	H (diff. cross-sections), $F_{\mathrm{m}}$, images of vessel, particle trajectory	Steering particles and avoiding stiction using time dept. fields
Hussain et al. [106]	2022		x	Straight 8 parallel lanes	2 EMs on opposite sides (Helmholtz coils)	B at fixed pos., velocity and translation velocity of particles, images of particles and fluorescence images along tubes	Controlling nanoparticle flow for simultaneous multichannel testing
Liu et al. [45]	2022	(x) *		Hybrid array Straight tube	6 EMs and 2 PMs	Magn. energy and $F_{\mathrm{m}}$ in horiz. and vert. dir. at fixed pos. and distance; no particles evaluated	Generating $F_{\mathrm{m}}$ in propagation dir. of particles for washing trapped particles out
Camargo et al. [107]	2023	x		Topology of blood vessel	One rectangular PM	Trajectory of particles, particle distribution	Parameter study, varying distance magnet to vessel, varying nanoparticle size
Durme et al. [28]	2023	x		Sym. Y-vessel mult. branches	4 EMs on opposite sides	Ferrofluid concentration over cross-section, particle trajectory, particle distribution	Optimizing particle steering, parameter study
Patel et al. [42]	2023	(x) *	x	In real mice	Structure with 4 PMs	B, $F_{\mathrm{m}}$, iron concentration at end of tube, fluorescent images of mice w. particles	MDT efficiency in tubes and mice, optimizing capture rate
Surpi et al. [36]	2023	(x) *	x	Straight tube	2 PMs on opposite side of vessel	H, images of tube, velocity of particles, magn. energy, particle energy, magnetic moment of particles	Controlling particle velocity in a straight tube
Thalmayer et al. [46]	2023	x		Sym. Y-vessel	Lin. hybrid Halbach array with 3 EMs and 4 PMs, $\vartheta ={90}^{\circ }$	B, $B_{\max }$ and grad${\left(B\right)}_{\max }$ in vessel, particle distribution	Investigating hybrid Halbach array and strength of its current for particle steering
Zhou et al. [108]	2023	x		Artificial vascular vessels, sym. Y-vessel mult. branches	Time-varying artificial magnetic field	Particle trajectory, particle distribution, particle speed	Steering particles using stochastic algorithms
This work	2024	x		Sym. Y-vessel	Halbach arrays consisting of PMs, 6 to 12 magnets, $\vartheta \in [{30}^{\circ },{135}^{\circ }]$	B, grad($B)$ at strongest pos., max. grad(B), magn. energy in vessel, $F_{\mathrm{m}}$, $F_{\mathrm{d}-\mathrm{m}}$, SPION distribution	Steering SPIONs, parameter study, investigate potential of diff. parameters for particle steering incl. strength and range of $F_{\mathrm{m}}$

* no simulations for particle steering.

. The steering results are not included in Table 6 as the performance depends significantly on the analyzed geometry, the overall setup, and parameters such as the investigated particle size. As these parameters differ considerably in the analyzed publications, comparing the steering performance is not meaningful.

As the simulation time of the particle study is quite long, and much computational effort is necessary, many researchers evaluated other magnetic parameters in order to predict the particle distribution. In Table 6, an overview of the evaluated parameters is provided. Here, e. g., the magnitude of the magnetic field $|\overrightarrow{H}|$ and magnetic flux density $|\overrightarrow{B}|$ are evaluated [16,25,31,36,42,54,99,105,106]. However, as Equation (8) depicts, the magnetic force is directly proportional to the gradient of B. Therefore, in other papers, the gradient of H or B is also analyzed [17,25,31,32,47]. Other researchers evaluated the magnetic force $F_{\mathrm{m}}$ [32,42,45,54,99,104]. Nevertheless, the gradients and $F_{\mathrm{m}}$ are usually analyzed at predefined evaluation positions, but, as shown in our previous paper [51], the magnitude of $F_{\mathrm{m}}$ strongly depends on the evaluation position. Another important point is that the actual attraction force also depends on the magnetic force distribution across the vessel cross-section. As shown in Figure 13, for example, the magnitude of $|F_{\mathrm{m}}|$ on the weak side is initially enormous but then decays quickly, which in turn means that the particles in the lower branch, which are to be pulled upwards, no longer experience any significant attraction force. Moreover, since the movement of the particles depends on the equilibrium of forces (Equation (5)), the magnitude of the drag force $|F_{\mathrm{d}}|$ must also be taken into account. However, this is not the case in many studies [32,42,45,54,99,104]. For this reason, the reliability of these findings should also be treated with care. In other studies, the magnetic energy is considered [36,45]. However, the results of this paper have shown that a statement about the particle distribution via magnetic energy should also be interpreted with caution (compare Figure 11).

In most state-of-the-art publications [17,25,28,30,31,32,46,47,102,103,105,106,108], the particle distribution is only evaluated for a fixed distance between the magnetic field source and the vessel. However, in this paper, we have shown that the distance between the magnetic field source and vessels has a significant impact on particle steering (compare Figure 12). Furthermore, in most publications (including our paper), the distance between the magnetic field source and the vessel is within only a few centimeters (or less), quite short. Surpi et al. [36] varied the distance up to 0.5 cm, Kee et al. [104] and Shiriny et al. [54] up to 1 cm, Patel et al. [42] up to 1.5 cm, Sarraf et al. [105] up to 2.5 cm, Bernad et al. [99] up to 3 cm, and Chakrabarty et al. [16] plus Liu et al. [45] up to 4 cm. Only Camargo et al. [107] investigated greater distances up to 7.5 cm. In MDT, this short distance between magnetic field source and vessel enables particle steering only very close to the body surface. Based on this, it can be concluded that steering magnetic nanoparticles in deeper body regions is still an open research question, which we addressed by highlighting the trade-off between the strength and usable range of the magnetic force.

Overall, in this paper, we conclude that, for the evaluation of particle steering, the particle distribution or the distribution based on ferrofluid concentrations must be evaluated. When measurements are conducted, camera images are often taken and used to extract the particle distribution, for example by using gray-scale images of the tube [16,32,36,99,104,106] or fluorescence [42,54,103,106]. In other publications, the particles are collected at the end of the tubes, and the particle concentration is determined from these samples [17,42,99]. This shows that the precise detection of particle concentration, especially during particle propagation, is also an open research topic. However, for steering magnetic nanoparticles in MDT, it is crucial to know where and when they are located. This information can then be fed back to the magnetic field via a control loop system [109,110].

Another clearly visible point in Table 6 is that the used magnetic source varies too. While simpler magnets are used, especially in publications in the medical context [42,103,107], many researchers aim to avoid particle trapping by placing two or more EMs on opposite sides of the vessel [16,25,28,31,32,102,106]. However, this may work under laboratory conditions, but, for an actual MDT treatment in the human body, it is not possible to arrange the magnets in this way. For this reason, many other authors [47,54,104,105], including those of this work, prefer using (linear) Halbach arrays from one side for particle steering.

5.5. Limitations of This Study

It is worth mentioning that the proposed study has some limitations: In order to reduce the simulation time and simplify the evaluation, a 2D geometry was investigated. However, reality is represented by a much more complex 3D geometry of real vessels. As the distance between the magnets and the upper boundary of the vessel is not constant for a cylindrical vessel in a 3D geometry, deviations in the particle distribution when using a 3D model may occur. Furthermore, in the simulations, the SPIONs were considered separately as individual particles. This means that the interaction of the SPIONs with each other was not considered. However, it has been shown in measurements that, when a magnetic field is applied, the particles form chain-like structures that are orientated in the parallel direction to the magnetic field lines [111]. This changes the effective permeability of the particle solution, resulting in an increased magnitude of the magnetic force [18,111]. Moreover, the interaction between the local particle concentration and the magnetic field was not considered. However, in [24], we have shown that the presence of magnetic particles changes the distribution of the magnetic field, resulting in a change in the magnitude of the magnetic gradient up to a factor of 2. In addition, the magnetization within the particles was assumed to be constant. Nevertheless, as, e.g., shown in [39,61], the actual structure of SPIONs is much more complex. Also, more effects like temperature changes or Brownian and Néel relaxation are not included in the simulation model. The impact of these effects is discussed in detail by Kolhatkar et al. [112]. Furthermore, in this study, as particle–particle interaction was not considered, only a limited number of 100 particles was taken into account. However, in actual MDT therapy, approximately ${10}^{10}$ particles per mL [9] with injection volumes of approximately 5 mL to 10 mL [12,113] are injected into the patient. Furthermore, the particles are injected slowly over several minutes in the artery via a catheter [114], whereas, in this study, all the particles were released simultaneously within one pulse. As demonstrated in [15,72], the parabolic velocity profile results in velocity dispersion of the particle pulse, which is assumed to be less evident when the particles are injected continuously.

Therefore, for simulating magnetic nanoparticles in a magnetic drug targeting scenario, models like ferrofluid or inhomogeneous media are also discussed [18,28,35,115,116]. Furthermore, there are also already approaches for integrating particle–particle interaction into simulation models based on measurement results [18] or particle red blood cells (RBCs) interaction [117]. To address these points, it is planned to expand and improve our simulation model in future studies.

Another limitation of this paper is that, for the magnetic force ${\overrightarrow{F}}_{\mathrm{m}}$, only its component in the y-direction was considered. However, of course, the magnetic force also has a component in the x-direction, which accelerates and decelerates the particles before and behind the magnet, respectively [98,99]. Nevertheless, as this component is much smaller, it can be assumed that the impact is negligibly small. Moreover, in this study, the background velocity was constant and set to water. However, the rheological behavior of blood is different, and the blood flow is pulsed. The interested reader can find a review on how to model the pressure and flow properties in vessels in the work of Zhou et al. [118]. Modeling of magnetic nanoparticles in blood flow was investigated numerically by Fanelli et al. [67], and the magnetic accumulation of SPIONs under arterial flow conditions in blood was measured by Hennig et al. [59].

6. Conclusions and Outlook

This paper conducted a systematic parameter study with a focus on varying the magnetization angles of linear Halbach arrays for steering SPIONs in a typical MDT scenario. The study was conducted numerically using COMSOL Multiphysics^® 6.1. In order to optimize the steering of SPIONs for future MDT therapies, the short usable range of the magnetic force is a limiting factor in addition to the precise adjustability of the magnetic field. For the strong side of the array, the results showed that large magnetization angles lead to greater magnetic forces. For the weak side, on the other hand, a weak magnetic force is desired for washing the SPIONs out of the vessel. This is the case for small magnetization angles. The reason for that can be found in the approximately exponential manner of the magnetic field and the resulting derivative. Either the gradient (and thus the magnetic force) decays quickly, or it remains at a high level for a greater distance. This results in a trade-off between a field with a long effective range (large angle) on the strong side and a field with a short effective range (small angle) on the weak side. Since simple manufacturing must be ensured for bringing this to practical implementation, we recommend a magnetization angle of $\vartheta ={90}^{\circ }$. In order to be able to precisely adjust the magnetic field for targeted therapy, we recommend to use hybrid arrays with electromagnets in addition to permanent magnets, as already presented in [45,46]. In future applications, it is planned to integrate a long array into a catheter and then insert the catheter near the tumor for precise steering of drug-loaded particles. Another challenge is to measure the position of the particles simultaneously to the steering process in order to adjust the magnetic field with the help of control loops [109,110]. For the authors’, ultrasonic (US) methods seem particularly promising here [119,120]. Furthermore, we concluded in this paper that particle steering cannot be predicted by simple quantities such as mean or maximum gradient or magnetic energy, which is common in many state-of-the-art publications, as we have shown in an extensive state-of-the-art analysis in Table 6. Thus, a numerical or experimental evaluation of the particle distribution is mandatory.

To optimize the steering or capture efficiency in magnetic drug targeting, parameters such as magnetic field intensity, particle size, and the magnetic permeability of the particles should be taken into consideration in future work. Also, the size and shape of the individual magnets can be investigated for the magnetic array. For further improving the simulation model, the magnetic nanoparticles should be modeled as an inhomogeneous medium including particle–particle interaction rather than individually.