Study of the effect of dipole interactions on hyperthermia heating the cluster composed of superparamagnetic nanoparticles

In the present work, we investigate the effect of dipole interactions on hyperthermia heating the cluster composed of multi superparamagnetic nanoparticles via time-quantified Monte Carlo simulation. The dynamic hysteresis loop area of non-interacting particles calculated by a modified Rosensweig’s model is shown to be proportional to the field frequency. The inverse of the total number of Monte Carlo steps per field cycle is considered as a computational frequency in our modelling. By comparing the two proportionality constants gained from the simulation and from the Rosensweig’s model, respectively, the time scale of one Monte Carlo step is estimated. The shape of the cluster is characterised by treating it as an equivalent ellipsoid. When the morphology of cluster is highly anisotropic such in a chain and cylinder, dipole interactions align the moments of the particles to the morphology anisotropy axis of the cluster. The strength of such alignment depends on the magnitude of morphology anisotropy of the cluster. The alignment helps improve heating capability of the chain and cylinder clusters at the most angles between the field direction and morphology anisotropy axis. However, when the field direction is away from the axis too much, the high energy barrier will hamper the cluster to maintain the magnetization, leading to a reduced heating efficiency. Once the cluster loses its morphology anisotropy (i.e. cube), the influence of dipole interactions on hysteresis losses is reduced to the minimum; the probability to obtain an improved heating becomes very low no matter with the type of particle arrangement.


I. INTRODUCTION
Much attention has been focussed on magnetic hyperthermia due to its great potential to be a cancer treatment. [1][2][3][4][5][6] Different from other approaches developed to raise the temperature of bio-tissue, 7-10 magnetic hyperthermia heats the cancerous tissues through inductively heating pre-implanted singledomain magnetic nanoparticles. The particles convert magnet work done by the field to particles' internal energy via the mechanism called "hysteresis losses." 11 The amount of energy converted per field cycle is equal to the area of the magnetic hysteresis loop. 11,12 Consequently, the temperature of particles could be increased rapidly at a rate of 10 2 to 10 3 K/s. 13 Then, the heat is released to the surrounding tissue to treat the tumor.
Many researchers have applied large single-domain particles retaining permanent magnetic moments at room temperature for hyperthermia heating. [14][15][16][17][18][19] Triggering the hysteresis losses of large particles requires a high threshold field amplitude stronger than the coercivity field. 20,21 However, it has been suggested that the field frequency f and amplitude H 0 should be reduced to some level in order to conduct a safe hyperthermia treatment. 2,22 Therefore, the use of small single-domain particles seem to be more important. 4 When the magnetic anisotropy energy becomes comparable with the thermal energy at room temperature, the moment of individual particles will start to flip freely through Néel Relaxation in absence of an external magnetic field; thus, the permanent magnetization will disappear; and those particles are referred to as superparamagnetic nanoparticles (SMNPs). 23 However, almost null coercivity field confines the hysteresis loop area. 11 Colloidal clusters composed of multi-magnetic nanoparticles have shown great potential in magnetic hyperthermia heating. 24,25 Different from the isolated particles, the hysteresis losses of the clusters are affected by inter-particle dipole interactions. The role of dipole interactions is complex, and apparently contradictory findings have been reported. Several evidences show that the dipole interactions can improve heating performance when the shape of cluster is highly morphologically anisotropic, such as in chain or cylinder. [25][26][27] Mehdaoui et al. 25 found that dipole interactions could improve the heating performance of SMNP columns by generating an additional magnetic uniaxial anisotropy, especially when the magnetocristalline anisotropy is low and column is long. However, other results suggest that chain-like clusterings may not bring in an enhancement in heating efficiency. 17,28 Branquinho et al. 28 reported that only at high damping factor the heating efficiency of short chain clusters would increase with the chain length; but at lower damping factor formation of chains could lead to low heating efficiency, because the present magnetocrystalline anisotropy became farther away from the optimal value although dipole interactions reinforce the cluster's effective magnetic anisotropy. The role of dipole interactions is also dependent upon the relative orientation of the magnetic field with respect to the clusters. Several modelling results suggest that chain cluster will produce less heat once the field direction does not align with the chain axis. 29,26 Serantes et al. 27 have experimentally demonstrated a reduction in heating efficiency by increasing the relative orientation of field with respective to the chain formed by ferromagnetic particles. Their simulation has predicted that the loss per cycle will decrease when the shape of the cluster is changed from chain to cube. In fact, assembling particles into anisotropy-less clusters, such as in sphere, can change the heating efficiency in either positive 24 or negative way. 30 In order to understand how dipole interactions dominate the hysteresis losses of a cluster, it is necessary to develop a method to associate the influence of dipole interactions with the cluster's characteristics as well as the magnetic field orientation; and it can be applied to the cluster with any shape and particle arrangement.
Metropolis algorithm is widely used for studying time-dependent magnetization of magnetic nanoparticle system by Monte Carlo (MC) method. 16,[31][32][33][34] Unlike the kinetic Monte Carlo algorithm, 35,36 which calculates dynamic hysteresis loop within two-level approximation, the transition possibility of the Metropolis algorithm mainly depends on the energy difference between the current and attempted states. However, without a theory to associate physical time with the Monte Carlo step (MCS), a dynamic magnetic hysteresis cannot be described in a quantified manner. Nowak et al. 37 pioneered the work on relating one MCS of heat-bath MC algorithm to the time scale of Langevin dynamics in fluctuation dissipation theorem. Cheng et al. 38 used the Fokker-Planck equation to link the Metropolis MC algorithm and Langevin dynamics schemes, and gained accurate quantification at a large range of the damping factor. Melenev et al. 39 estimated the time scale of MCS for simulating dynamic magnetic hysteresis for the first time. They considered a model of non-interacting SMNPs responding to an oscillating field. The MC simulation was manipulated until the simulated hysteresis loop was in good agreement with the one produced by Brown's kinetic equation at a certain frequency. By correlating MCSs per field cycle, N MCS , with the frequency, the time scale of each MCS was obtained. One can expect an improved accuracy of the time quantification if a series of frequencies are considered.
In the present work, the effect of dipole interactions on the hysteresis losses of the cluster composed of SMNPs is studied. Three typical morphologies, namely chain, cylinder and cube, are selected to shape the cluster. Two cube clusters are fabricated to possess a simple cubic and facial centred cubic (FCC) lattice, respectively. The morphology anisotropy of the cluster is quantitatively characterised by using the reported method. 40 An approach has been proposed to quantify the directionally dependent influence of dipole interactions in terms of energy. The dynamic magnetic hysteresis of a cluster is simulated via the MC method. The hysteresis loop area of non-interacting particles calculated by Rosensweig's model 12 with an minor modification is shown to be proportional to the field frequency. In our modelling, the inverse of N MCS per cycle is considered as a computational frequency. By comparing the two proportionality constants gained from the simulation and the Rosensweig's model, the time scale of one Monte Carlo step is estimated. Our simulation results suggest that when the cluster is highly morphologically anisotropic, the dipole interactions try to align the individual magnetic moment of each particle to the axis of shape anisotropy. The strength of the alignment associates intimately with the magnitude of the morphology anisotropy. Once the cluster loses its morphology anisotropy, the influence of dipole interactions on the hysteresis loss will be reduced to the minimum and the chance of obtaining improved heating will become very low.

A. Characterization of cluster's morphology anisotropy
The morphology anisotropy of the cluster is characterised by treating it as an equivalent ellipsoid with uniform mass distribution in 3D. The ellipsoid has the same mass and principal moments of inertia as the cluster. The morphology anisotropy is defined in terms of the ratios of the length of its semi-principal axes. Each particle is divided into 280 unit volume elements. The moment of inertia tensor, T I , for a cluster composed of N discrete point masses (N = 280 × the number of particle) is described as, where (x l , y l , z l ) are the coordinates of point mass l. The original point locates at the mass centre of the cluster. Diagonalising T I works out the principal moments of inertia of the cluster, λ i (i = 1, 2, 3). The length of semi-principal axis of the equivalent ellipsoid d i is obtained by, where i, j and k equals 1 or 2 or 3, respectively, and i j k. Then, d i (i = 1, 2, 3) are assigned to A, B and C, ordered from the largest to smallest. A long d i comes with a small λ i , which suggests that the mass of the cluster is distributed close to the principal axis of rotation i. Therefore, high ratios of A to B and C indicate a heterogeneous mass distribution of the cluster with respect to different axes, indicating a high anisotropy in morphology.

B. Estimation of directionally dependent influence of dipole interactions
When a cluster is magnetically saturated, the total energy of dipole interaction where µ 0 is the vacuum permeability, M d the domain magnetization, V the magnetic volume of particle, θ i j the angle between the direction of magnetization and the line joining the centers of particle i and j, and ⃗ r i j the distance between the centers of particles. We denote it as E DISM . It is evidenced that the clusters formed by SMNPs can inherit superparamagnetism from the particles, 24,41-43 as long as T is high enough to make the moments rotate too rapidly to generate permanent magnetization. Moreover, the blocking temperature of the densest assembly of SMNPs can be much lower than room temperature. 44,45 When body temperature is considered, it can be assumed that  i < j E (i, j) D of the SMNP cluster whose magnetization is completely relaxed is independent upon its structure and shape. In this case, E DISM dominates the difference of  i < j E (i, j) D before and after the cluster is magnetized. Decreasing E DISM is expected to make magnetizing the cluster easier. In the present modelling, the 127232-4 Fu, Yan, and Roberts AIP Advances 5, 127232 (2015) direction of the cluster's magnetization is set always in parallel with the field direction. Therefore, E DISM is such a parameter which connects dipole interactions, the field direction, and the hysteresis losses.

C. Monte Carlo Simulation
All clusters are composed of mono-sized spherical magnetic nanoparticles. Each particle possesses a uniaxial magnetic anisotropy, and the orientation of the easy axis is randomly chosen in a 3D space. The magnetic properties of particle are set the same as the published data of magnetite nanoparticles: the effective magnetic anisotropy constant K eff is 9000 J/m 3 and M d 446 kA/m. 46 For the sake of simplicity, we assume that both of them are temperature-independent. The magnetic moment of each particle is defined as The energy model of the cluster system consists of three major sources, namely, anisotropy energy E A , Zeeman energy E H and dipolar interaction energy E D . 47 The uniaxial anisotropy E (i) A of each particle is given by, where ⃗ n i is the unit vector along the easy axis direction. The coupling with the applied field ⃗ H is described by, The energy of dipole interaction is given by, Therefore, the total energy of the cluster system can be expressed as, Five kinds of SMNP clusters are studied in this work. Chain cluster, cylinder cluster and cube cluster with simple cubic lattice are built with 64 particles. Cube cluster with a FCC lattice contains 63 particles. The particles of the chain cluster are lined up along the Z-axis. The cylinder cluster is created by repeating the unit cell of simple cubic lattice along the Z-axis. The anisotropy-less cluster with an imperfect lattice is fabricated by removing 16 randomly picked particles from the structure of the cube cluster with a simple cubic lattice. Usually, a SMNP is covered with a layer of organic stabiliser in practice to ensure the particle dispersity. The thickness of this stabiliser layer is set at 1 nm in the modelling, and the radius of magnetic fraction is kept at 5 nm.
Monte Carlo simulation featured with the Metropolis algorithm is carried out to reproduce the time dependent magnetization of SMNP clusters at body temperature. At the beginning of each MCS, one particle is picked randomly and the moment is directly agitated to a new direction, which is chosen inside of a spherical segment around the present direction with an aperture angleδθ. According to the reported work, 16 temperature dependence of δθ is given by δθ = (0.05k B T/2K eff V ) 0.5 in an usual reduced unit, where k B is the Boltzmann constant. This agitation is accepted with probability min [1, exp(−∆E/k B T)], where ∆E is the change of the total energy of the cluster system caused by agitation. The above procedure is repeated until all particles are agitated to complete one MCS. During the simulation, the particles only relax through Néel mechanism and the particle positions are fixed.
Before simulating cluster's hysteresis losses, 275000 MCSs are used for thermalization. After that, magnetic field is introduced and increased until reaching H 0 ; then it is decreased to −H 0 , and increased again to H 0 to finish the cycle.

III. RESULTS AND DISCUSSION
Figure 1(a) shows 3D schematics of chain cluster, cylinder cluster, and two cube clusters which have a simple cubic and FCC lattice, respectively. The Morphology anisotropy of SMNP clusters is characterised by treating it as an equivalent ellipsoid. As shown in Figure 1(b), the equivalent ellipsoids of the chain cluster and cylinder cluster exhibit similar characteristics to the morphology of clusters. Because the mass distribution of both clusters is rotationally symmetric with respect to Z-axis, the two shorter semi-principal axes B and C equal each other (Table I). The ratio A/B or A/C of chain cluster is 10 times larger than that of cylinder cluster. The longest principal axes of both ellipsoids coincide with Z-axis. In contrast, the equivalent ellipsoids of the two cube clusters are visually isotropic in geometrics (Figure 1(b)). Obviously, the finding of A/B=A/C=1 (Table I) further confirms that the mass distribution of these two cube clusters is truly isotropic, regardless of the type of lattice. Therefore, as the shape changes from a chain to a cylinder and to a cube, the morphology anisotropy decreases to 0. The longest principal axis is arbitrarily defined as the axis of morphology anisotropy (AMA). As shown in Figure 2, the orientation distribution of E DISM of chain and cylinder clusters displays typical uniaxial anisotropy. The easy axis coincides with the AMA of each cluster (red arrows). The angle-dependent E DISM obeys the correlation, E DISM = K 0 + K 1 sin 2 θ M A , where K 0 and K 1 are constants and θ M A is the angle between the unit vector of magnetization of the cluster and the AMA. The value of K 0 and K 1 are given in Table I. K 1 of the chain and cylinder cluster are positive, indicating that E DISM of both the clusters increases with θ M A . As θ M A gets close to 0 o the dipole interactions are expected to help maintain the current magnetization of the cluster. However, once θ M A approaches to 90 o , the dipole interactions will turn to compel the cluster to lose the magnetization. Moreover, K 0 of the chain cluster is smaller than that of the cylinder, but K 1 is larger. The deep energy valley brought by E DISM at 0 o makes the moments of the chain cluster more inclined to align with the AMA. In contrast, E DISM of the two cube clusters are independent of direction ( Figure 2) and fixed at 0. This is further validated by the presence of almost null K 0 and K 1 (Table I). Therefore, after the cluster loses its shape anisotropy, dipole interactions are expected to contribute little to the cluster's magnetism.
To investigate the effect of dipole interactions on hysteresis losses of the cluster, Monte Carlo approach featured with the Metropolis algorithm is carried out to generate M-H curve. At first, the simulation is validated by reproducing the DC equilibrium magnetic magnetization of non-interacting SMNPs at T ≈ 0 K, as shown in Figure 3 (left). The remanence magnetization and coercivity are 0.5 M S and 0.48 H a respectively. Both of them are in good agreement with the expected values proposed by Stoner and Wohlfarth. 16,48 Then, the temperature is increased to 310 K (body temperature), and the particle moments are allowed to interact magnetically with each other. N MCS per cycle is adjusted to be large enough to obtain a DC equilibrium magnetization. As shown in Figure 3  all of the M-H curves exhibit minute loops, demonstrating the presence of superparamagnetism, despite the cluster's shape. The DC equilibrium magnetization of non-interacting particles at 310 K is also shown in Figure 3 (right), which is calculated by the LRT. It is clear that LRT is no longer suitable for describing the magnetization of an ensemble of magnetically interacting particles. This notion is consistent with Schaller et al., 49 who used MC to simulate DC equilibrium magnetization of a sphere-like cluster at 293 K. In order to simulate AC dynamic hysteresis of SMNPs in a quantified manner, the estimation of the time scale of MCS is demanded. Rosensweig's model has been suggested to be useful in predicting heating behavior of non-interacting SMNPs. 4 In the present modelling, Rosensweig's model is modified to be suitable for describing dynamic hysteresis driven by high field intensity. The equilibrium susceptibility is set to change with field intensity, and the loop area is calculated by integrating magnetization against field. Details are given in the Supplementary Information. 52 The properties of the particle and field amplitude are set to be the same as those used in MC simulation. The value of pre-exponential factor τ 0 is chosen as 10 −9 s −1 , which corresponds to the damping factor of 0.28 (see Supplementary Information 52 for details). During the simulation, the temperature is maintained at 310 K and the value of δθ is fixed. Figure 4 shows the plot of simulated loop area in reduced units against 1/N MCS per cycle. It can be seen that the loop area is excellently proportional to 1/ N MCS per cycle. By using linear regression without the intercept term, the proportionality constant is found to be  . Dividing the second constant by the first one, we find that one MCS is comparable to 3.67 × 10 −11 s. 54800 MCSs per cycle is used for the study of the hysteresis losses of chain, cylinder and two cube clusters, which corresponds to 500 kHz. Figure 5(a) and 5(b) show M-H curves of chain and cylinder clusters obtained at different θ M A , respectively. The hysteresis loops of both clusters shrink quickly with θ M A . This is further demonstrated by continuous reduction of the loop area, denoted by A M in Figure 5(c). According to the above discussion, both E DISM of chain and cylinder clusters increase with θ M A . When the moments are required to overcome a high energy barrier to get aligned with the field, a reduced heating efficiency is understandable. The blue dashed line in Figure 5(c) presents the A M of non-interacting SMNPs as blue dash line. Before θ M A reaches 75 o , A M of chain and cylinder cluster remain higher than that of non-interacting SMNPs. If we assume that the cluster is randomly orientated with respect to the field direction in practice, there will be 80% of probability to get an enhanced heating efficiency after assembling SMNPs in the chain and cylinder cluster. In terms of the particle's properties and τ 0 chosen for the simulation, the damping factor of isolated particle is around 0.28, which is already higher than that of the particles studied in work of Branquinho et al. 28 Their work also suggests that the damping factor tend to increase with chain size. Perhaps high damping behavior might exist in our simulation and prompt dipole interactions to contribute positively to the chain cluster's hysteresis losses.
Interestingly, the chain cluster A M remains superior to the highest of the cylinder A M until θ M A reaches 60 o (Figure 5(c)), but E DISM of the latter at 0 o is already lower than the former at 30 o ( Figure 5(d)). It seems that the hysteresis losses of the chain and cylinder cluster is not completely dominated by E DISM . During the simulation, the projections of moments in [001] direction are collected simultaneously. The loop area in reduced unit, denoted by A 001 in Figure 5 Fig. S3). This indicates that the moments of both clusters have a clear tendency to get aligned with each AMA. Furthermore, the alignment strength is enlarged when the shape is changed from a cylinder to a chain. Before θ M A rises to 90 o , A 001 of the chain cluster remains higher than the other's. The growth of E DISM as θ M A makes the moments more inclined to stay at the AMA to minimize the energy of dipole interactions. Such an alignment to a certain axis should impede the demagnetization of the cluster significantly. It can be noticed that the energy valley of a chain cluster formed by E DISM is sharper than the cylinder one (see Figure 5(c)). With the aid of such strong alignment, the hysteresis losses of the chain cluster remain superior to the other at most θ M A .
However, to answer whether this alignment brought by dipole interactions is the additional uniaxial anisotropy mentioned in the work of Mehdaoui et al., 25 an investigation on the relationship between the alignment strength and θ M A will be necessary.
On the other hand, once the cluster loses its morphology anisotropy, we find that the hysteresis losses begin to react almost numbly to the change of field orientation with respect to cluster. The field direction is specified with cartesian reference frame. As shown in  Figure 6(c) compares the loop area of these two cube clusters in these four directions. The loop area of the cube cluster with a simple cubic lattice varies little with the field direction and is slightly lower than the non-interacting particles. The average loop area of the other cube cluster over four directions is exactly equal to that of the non-interacting particles. The finding of inefficient hysteresis losses caused by losing shape anisotropy is consistent with Mehdaoui et al. 25 and Saville et al. 16 E DISM of both cube clusters are directionally independent and kept at 0. The nearly null contribution of dipole interactions minimizes the difference between clusters with less morphology anisotropy and non-interacting particles.  So far, the synthesis of sphere-like SMNP clusters with a disordered structure has been widely developed. 42,50,51 Hayashi et al. 24 found that the heating efficiency was increased by 50 percent after assembling 9 nm of magnetite nanoparticles into sphere-like clusters with a disordered arrangement. However, Liu et al. 30 observed a sharp decrease in the heating ability after forming a dense spherical pack of particles by loading 6 nm of MnFe 2 O 4 into a polymer latex. It can be expected that a disordered particle arrangement might bring an impact on the cluster's hysteresis losses. Here, we intend to study how the loss of lattice perfection changes the loop area of anisotropy-less cluster. 50000 such clusters with imperfect lattice are fabricated by omitting 16 randomly picked particles from cube clusters with a simple cubic lattice. As a result, the average A/B and A/C are 1.06 and 1.14 with a relative standard deviation of 2.88% and 4.55% (Table I), respectively. It is noted the E DISM can achieve the minimum at any angle with respect to the AMA (see Supplementary Information 52 Fig. S4).
To find out the relationship between E DISM and the hysteresis losses of anisotropy-less clusters with imperfect lattice, numerical simulations are carried out with setting the field to oscillate along the Z-axis. Figure 7(a) shows the plots of the loop area against E DISM when N MCS is adjusted to 20000, 30400 and 54800, which correspond to 1300, 900 and 500 kHz, respectively. The dashed lines represent the loop area of non-interacting particles at these three N MCS . It can be seen that at 54800 MCSs per cycle the E DISM influences the loop area slightly, and the loop area remains comparable to non-interacting particles. When the N MCS per cycle is reduced to 30400, the loop area turns to decrease with E DISM . This tendency becomes to be much obvious at 20000 MCSs per cycle. In other words, only at an extremely high frequency, the dipole interactions could be expected to change the hysteresis losses, and it requires E DISM to be lower than −1.0 × 10 3 J/m 3 to produce a better heating efficiency. Figure 7(b) shows the probability distribution of E DISM when magnetization is along with [001] direction and the distribution averaged over 1000 directions uniformly distributed in 3D space. Both distributions suggest that the chance of E DISM being lower than −1.0 × 10 3 J/m 3 is quite low (< 3%). Obviously, it is difficult for a cluster with less anisotropy in shape to convert more magnetic work to heat.

IV. CONCLUSION
In this work, we investigate the effect of dipole interactions on hyperthermia heating SMNP clusters via time-quantified Monte Carlo simulation. The shape of the cluster is characterized by treating it as an equivalent ellipsoid. The morphology anisotropy is defined in terms of the ratios of the length of the ellipsoid's semi-principal axes. As cluster's shape is changed from chain to cylinder and to cube, the morphology anisotropy of the cluster decreases to 0.
According to the Rosensweig's model, the hysteresis loop area of non-interacting particles is supposed to be proportional to the field frequency. In our simulation, the inverse of total number of Monte Carlo step per field cycle is considered as a computational frequency. By comparing the two 127232-11 Fu, Yan, and Roberts AIP Advances 5, 127232 (2015) proportionality constants gained from the simulation and from the Rosensweig's model, respectively, the time scale of one Monte Carlo step is found to be 3.67 × 10 −11 s. E DISM suggests the total energy of dipole interactions after a cluster is magnetically saturated. It is gained by direct calculation of Equation (3). Decreasing E DISM is expected to make magnetizing the cluster easier. The orientation distribution of E DISM of chain and cylinder clusters displays typical uniaxial anisotropy. When the field direction is parallel to the morphology anisotropy axis of the cluster, the value of E DISM achieves the minimum; and E DISM increases with the angle between field direction and this axis. However, once the cluster loses the morphology anisotropy such as cube cluster, the orientation distribution of E DISM becomes isotropic and the value is fixed at 0, indicating that a minimized influence of dipole interactions can be expected.
The results of the present simulation suggest that the hysteresis loop area of chain and cylinder cluster decreases continuously with the increase of the angle between field direction and the morphology anisotropy axis of the cluster, due to the growth of E DISM . However, with the aid of low E DISM , the hysteresis loop area of the chain and cylinder cluster remain superior to that of noninteracting particles at the most angles. In addition, higher shape anisotropy also prompts the dipole interactions to benefit cluster's hysteresis losses via aligning the moments to the cluster's morphology anisotropy axis. Such alignment helps chain cluster to heat better than cylinder cluster, even when E DISM of the former is lower than the latter. However, once the cluster loses its shape anisotropy, it will be hard to obtain enhanced heating efficiency. The loop area of cube cluster is almost the same as that of non-interacting particles regardless of the change of the lattice type from simple cubic to FCC lattice as well as the field direction. Moreover, the loss of lattice perfection won't enlarge the probability for anisotropy-less cluster to heat better. Only at extremely high frequency the E DISM can get the chance to affect the hysteresis losses of anisotropy-less cluster with imperfect lattice.

ACKNOWLEDGMENTS
The project is supported by University of Nottingham Engineering scholarship.