Hybrid Molecular and Spin Dynamics Simulations for Ensembles of Magnetic Nanoparticles for Magnetoresistive Systems

The development of magnetoresistive sensors based on magnetic nanoparticles which are immersed in conductive gel matrices requires detailed information about the corresponding magnetoresistive properties in order to obtain optimal sensor sensitivities. Here, crucial parameters are the particle concentration, the viscosity of the gel matrix and the particle structure. Experimentally, it is not possible to obtain detailed information about the magnetic microstructure, i.e., orientations of the magnetic moments of the particles that define the magnetoresistive properties, however, by using numerical simulations one can study the magnetic microstructure theoretically, although this requires performing classical spin dynamics and molecular dynamics simulations simultaneously. Here, we present such an approach which allows us to calculate the orientation and the trajectory of every single magnetic nanoparticle. This enables us to study not only the static magnetic microstructure, but also the dynamics of the structuring process in the gel matrix itself. With our hybrid approach, arbitrary sensor configurations can be investigated and their magnetoresistive properties can be optimized.


Introduction
The well-known giant magnetoresistance (GMR) effect is commonly used to design extremely sensitive sensors that respond to external magnetic fields with changes in their electrical resistance. The GMR effect was originally discovered in magnetic multilayer systems [1,2]. Later it was also found for systems that contain magnetic grains in metallic matrices [3,4]. Recently, it could be shown that systems made of magnetic particles that are immersed in conductive gel matrices show GMR behavior as well. Besides achieving very high GMR effect amplitudes, the use of gel matrices opens up the perspective of printable, low-cost, magnetoresistive sensor devices such as sensors for the detection of biomolecules [5][6][7]. The sensor characteristics directly depend on the microscopic behavior of the magnetic particles in the gel matrix, i.e., on the dynamics of the magnetic moments of the particles in combination with their local positions and motion in the gel. Since such microscopic information cannot be directly extracted from experiments, numerical simulations have to be taken into account. In order to simulate the magnetodynamics of an ensemble of magnetic nanoparticles that are embedded in a conductive gel matrix, two different simulation methods have to be combined in order to simultaneously consider both the magnetic and translational degrees of freedom of the magnetic particles. In order to do so, one has to consider the time scales of the two classes of degrees of freedom as well as a method to control the temperature. Here, we present a hybrid approach that couples molecular dynamics and classical spin dynamics methods so as to simulate the structuring process of spherical magnetic nanoparticles in a viscous surrounding medium. In contrast to other hybrid molecular and spin dynamics algorithms presented in the literature [8][9][10][11][12][13][14], a purely classical approach is followed here. Moreover, the known methods follow completely different approaches and, therefore, instead of a coupling of the spin and mechanical degrees of freedom, couplings of spin and other degrees of freedom are realized. For example, in [8,9], a hybrid method is presented that covers the spin and electronic degrees of freedom whereas in [13] a coupling between spin and lattice dynamics is introduced. In addition to that, we follow a coarse grain approach, whereas the other methods pursue an atomistic model. Our hybrid molecular and spin dynamics approach has been introduced in [6] as one part of a simulation tool chain. Therein, a basic description is given without computational details. Here, we provide a complete and comprehensive derivation of our approach. Starting from the general approach, we briefly describe both numerical methods separately. Subsequently, we present the theoretical foundations and introduce our hybrid algorithm. Finally, we show one exemplary simulation run for didactic purposes.

General Approach
In this paper we assume magnetically interacting nanoparticles that are free to move in a viscous medium. In order to simulate their behavior one has to solve the following sets of coupled equations of motion: Here, the first and second set of equations is used to compute the velocity and position of each particle in the ensemble. The force describes the viscous drag of the surrounding fluid. The term ℋ is the magnetic potential due to the dipole-dipole interaction among all particles which depends on the positions of the particles as well as on the orientations of their magnetic moments . The last term of the first equation describes the repulsive force due to the Weeks-Chandler-Andersen potential which is necessary in order to simulate hard sphere dynamics. The third set of equations calculates the orientations of the magnetic moments . Here, denotes an effective magnetic field which is generated by the magnetic interactions between the particles. The last term is the so-called Landau-Lifshitz damping which is non-conservative and used to find the appropriate magnetic low energy configurations for a given particle arrangement (see below). Here, describes the friction and can be adjusted in order to reach a sufficient numerical stability. In summary, there are three important potential or force contributions. These contributions are described in detail in Section 2.1.2.
Although the equations of motion (1-3) are fully coupled by the dipole-dipole potential ℋ one can show that under certain conditions (see below) it is not necessary to solve these equations simultaneously but to separate their calculation leading to a unidirectional consecutive approach. In this approach, the magnetic configuration, i.e., the orientation of the magnetic moments, for a given initial nanoparticle arrangement is calculated first. This magnetic configuration creates magnetic forces that act on the particles inducing their subsequent motion which in turn changes the orientation of the magnetic moments and so on, resulting in a hybrid approach. The advantage of our hybrid approach is that the calculation of the magnetic low energy configuration and the calculation of the resulting translational motion can be executed separately by two different simulation methods which are explained in the following.

Molecular Dynamics Equations of Motion
Classical molecular dynamics (MD) provides a method to simulate the motion of atoms, molecules or particles due to interactions between the particles or interactions with external fields by solving the following set of equations of motion: In contrast to ab initio methods which are based on quantum mechanical principles, classical MD is based on classical potential formulations from which the according forces on the particles are derived. MD can be applied to a broad range of fields, e.g., theoretical physics, physical chemistry, biophysics and materials science. In this work, we have extended the open-source code HOOMD-blue [15][16][17] with the necessary customizations as given below.

Integration of the Equations of Motion
Particles in classical MD systems obey the laws of Newtonian mechanics. There exists a multitude of standard algorithms to integrate the Newtonian equations of motion (4,5), i.e., the Euler algorithm, different types of Runge-Kutta schemes and Verlet-type methods [18]. Here, we use one of the most important algorithms for the integration of the equations of motion, namely the Velocity Verlet algorithm which is a variant of the well-known Leapfrog algorithm. In contrast to other integration algorithms, positions and velocities are calculated for the same steps in time which is a great advantage of the Velocity Verlet algorithm. Moreover, this algorithm exhibits good energy conservation and is time-reversible. Equations (6) and (7) show how the positions and velocities are calculated for each particle [18,19]:

Force Calculation
The forces that act on the particles due to particle-particle or particle-field interactions have to be calculated in every time step of an MD simulation. The calculation of particle-particle interactions is extremely time-consuming because it is of order . In some cases the computational effort can be reduced by means of neighbor list or cell computation approaches [18]. In our case we assume that the magnetic nanoparticles solely interact by magnetic dipole-dipole interaction which exhibits particular long-ranged characteristics. In order to preserve this long-range behavior, methods to reduce the computational effort cannot be applied because all of these approaches utilize a potential cut-off. Therefore, the full calculation of all pairwise contributions is used here which provides correct forces for short and long distances between the particles [18,19]. In order to model magnetic nanoparticles that move due to their magnetic interactions while being immersed in a viscous matrix, three contributions to the net force must be considered: Magnetic Dipole-Dipole Interaction Every magnetic particle in the system under consideration creates a magnetic dipole field which exerts a torque on the dipole field of every other particle. The pairwise potential contributions can be calculated according to the following Equation (8): Therein, and represent the magnetic moments of two particles i and j while describes the unit vector that connects the centers of the particles. The distance between the centers of the particles is given by . The magnetic dipole-dipole interaction causes a force that acts along . Depending on the relative orientations of the magnetic moments of the particles, it can be attractive or repulsive. An example of the distance dependence of the magnetic dipole-dipole interaction energy of two interacting cobalt particles with diameters of 10 nm is depicted in Figure 1. It is evident that a well-defined cut-off radius does not exist due to the long-range character of the interaction.

Figure 1.
Magnetic dipole-dipole energy curves for the case of two interacting cobalt particles with diameters of 10 nm. The magnetic dipole-dipole energy is plotted vs. the center-to-center distance of the particles. For the case of particles with magnetic moments aligned in parallel, negative energy values are obtained (black curve), resulting in particle attraction. In contrast to this, for an antiparallel alignment of the magnetic moments, positive energy values are obtained (red curve), resulting in particle repulsion. It becomes apparent that the magnetic dipole-dipole potential contribution cannot be cut off due to its long-range character over a distance of several particle diameters even for small particles.

Hard Particle Approach
In general, point masses are used in MD simulations. However, this means for the particle positions that they can become identical or that distances can become infinitesimally small, i.e., → 0. In such cases the dipole-dipole potential Equation (8) diverges which would lead to numerical instabilities and unphysical behavior. In order to avoid this, potential functions for molecular dynamics simulations of hard particles have to be used which are based on a combination of attractive and repulsive contributions. The so-called Weeks-Chandler-Andersen (WCA) potential provides a force-shifted formulation of the well-known Lennard-Jones potential which provides the required short-range repulsion [20]. Above a distance of two particle diameters ≥ 2 , the WCA potential is cut off. Below this distance, a Lennard-Jones type repulsion is modeled: Here, the parameter σ is given in units of distance and ε is given in units of energy. Figure 2 illustrates the WCA potential for the case σ = ε = 1. (red) and −4 × (black). Here, the respective potential ℋ which is given in units of energy is plotted as a function of the distance between two fictitious particles given in units of distance.

Stokes Drag
To take into account the surrounding medium that affects the motion of the particles due to its viscosity, a third force contribution must be included. Stokes' law describes the force that acts on a small spherical particle with a radius and with a given velocity that is embedded in a viscous medium with kinematic viscosity . The resulting drag force acts in the direction opposite to the particle's motion and it is proportional to its velocity and the viscosity of the medium. = 6 (10)

Classical Spin Dynamics Simulations
Classical spin dynamics is a proven method for the numerical analysis of the static and dynamic thermodynamic properties of microscopic and mesoscopic ensembles of classical magnetic moments. Here, we use a spin dynamics algorithm that involves solving the Landau-Lifshitz equation. A Langevin approach is utilized to model the contact to an external heat bath. Moreover, the particles are assumed to be homogeneously magnetized and therefore modeled by a single magnetic moment. A detailed description of the algorithm is presented in [21].
Spin Equations of Motion: A magnetic moment that is subject to an external magnetic field precesses around the external field direction. In a classical picture this precession is given by: where ℋ is the Hamiltonian that contains the magnetic interactions, i.e., exchange, Zeeman, anisotropy and dipole-dipole interactions. In this paper we only consider the dipole-dipole interaction. For systems of magnetic nanoparticles that are exposed to a homogeneous external magnetic field, the Zeeman contribution has to be included as well [21]. It is convenient to use the partial derivative of the Hamiltonian H, i.e., ℋ = which is the effective magnetic field caused by all interactions that act on the spin under consideration. Equation (11) gives rise to a conservative dynamics. In order to be able to calculate low energy configurations, Equation (11) needs to be extended by a damping term which allows the system to relax towards its ground state. This is commonly done by adding a damping term as proposed by Landau and Lifshitz: In doing so, the spins relax until they are oriented parallel to the effective field direction. The damping force which is proportional to the positive damping constant , acts in the direction of the effective field [21]. By solving Equation (12) it is not guaranteed that the system's ground state is found. In fact, for more complex systems the algorithm "gets stuck" in the local energy minimum that is first reached in phase space and there is no way to get out of this. In order to overcome this problem one can extend Equation (12) by a heat bath coupling. Here, we use a Langevin approach in which the heat bath acts on the system by means of stochastic forces. Hence, the classical spins fluctuate due to forces that are driven by the temperature of the heat bath and "kick" the system out of local minima. The fluctuations are introduced into the equations of motion by the additional term: with ( ) representing the undirected fluctuations that have a white noise characteristic. As a result one obtains the following stochastic differential equation of motion for the spin degrees of freedom: In order to integrate Equation (14) we have used a fourth-order Runge-Kutta scheme as proposed by Milstein and Tretyakov [21,22].

Hybrid Molecular and Spin Dynamics Simulations
As pointed out in Section 2, the simulation of interacting magnetic nanoparticles in a viscous medium generally requires solving the system of coupled differential Equations (1)-(3) for the translational and magnetic degrees of freedom. However, this is only necessary if the time scales of the dynamics for the translational degrees of freedom and magnetic degrees of freedom are of the same order. If this is not the case the magnetic and translational degrees of freedom can be treated separately, leading to a simple coupling scheme reminiscent of the method proposed by Dünweg and Ladd for the numerical simulation of colloidal dispersions [23]. In order to validate this assumption for our problem, the relaxation times of the two types of degrees of freedom are determined in the following.

Translational Relaxation Time
In order to determine the time scale of the translational relaxation of a particle of spherical shape, the configurational relaxation time [23] must be evaluated. The configurational relaxation time, or Brownian relaxation time, is the time that the particle needs to diffuse across its own radius and can be evaluated according to: Here, represents the kinematic viscosity of the matrix material, is the particle diameter, is the Boltzmann constant and represents the temperature [24].

Magnetic Relaxation Time
The magnetization of a colloidal magnetic particle relaxes due to two different mechanisms. First, the relaxation can occur via rotation of the particle in the liquid. The characteristic time associated with this mechanism is called Brownian rotational diffusion time and can be calculated according to Equation (16) with representing the particle volume: The second mechanism is the so-called Néel rotation where the magnetization of the particle relaxes due to the coherent rotation of atomic magnetic moments within the particle. For a uniaxial, single-domain, ferromagnetic particle, the associated time can be calculated by: where, represents the magnetic anisotropy constant of the particle material and is a characteristic frequency [25]. As shown in Figure 3, for small particles the Néel relaxation dominates whereas for larger particles, the Brown relaxation dominates. In addition to that, a transition regime can be identified around particle diameters of about 20 nm. Because both mechanisms contribute to the relaxation behavior of ensembles of magnetic particles, an effective relaxation time can be defined according to the following Equation [26]:

Comparison of Translational and Magnetic Relaxation Times
In order to check whether the translational and the magnetic degrees of freedom can be separated, the time scales have been evaluated for cobalt particles with diameters up to 100 nm that are immersed in agarose gel matrices with a common concentration. In our current algorithm we do not consider a difference between Brownian and Néel rotation but will address this in further studies. Thus, the relevant time scales are set by the configurational relaxation time and both the Néel rotation diffusion time and the Brownian relaxation time that are combined to the effective relaxation time .
For we have chosen a characteristic frequency of = 10 Hz and an anisotropy constant of = 10 J • m [6,27]. For we have used = 0.017 Pa • s for the kinematic viscosity of the gel matrix which corresponds to a 2% agarose gel matrix as it is used in [5,7]. The results are shown in Figure 3. Here, it becomes apparent that the Néel relaxation time is very small for small particle diameters and rapidly increases with increasing particle diameter. In contrast to that, the Brown relaxation time slowly increases with increasing particle diameter. By combining both relaxation times according to Equation (18), the relevant effective magnetic relaxation time is obtained. A comparison of this effective relaxation time and the mechanical, configurational relaxation time shows that a separation of the time scales exists for all particle diameters except for very small diameters up to around 3 nm. The smallest difference between the time scales for larger particles can be identified around the transition regime for particle diameters of about 20 nm. Nevertheless, even at this point, the time scales are clearly separated by at least one order of magnitude. Within our algorithm, a difference between Brown and Néel relaxation cannot be addressed and hence, both mechanisms result in a change of the magnetic moment orientation without changing the orientation of the particle itself. Hence, the configurational relaxation time must simultaneously be compared to the magnetic relaxation times that are combined to the effective magnetic relaxation time (orange) according to Equation (18). There is a significant difference between the magnetic and mechanical relaxation over all particle diameters except for very small particle and one small window around the particle diameter of 20 nm.

Hybrid Molecular and Spin Dynamics Coupling Scheme
In our hybrid coupling scheme which is depicted in Figure 4, we start from an initial random magnetic moment configuration . In principle, the magnetic moments are three-dimensional vectors. However, due to their initial coplanar arrangement the dipole-dipole interaction leads to a coplanar orientation of the magnetic moments as well. The initial particle positions are taken from experimental data, e.g., microscopic images. Using and the magnetic moment orientations are relaxed towards a low energy state by means of spin dynamics (SD) simulations. Hence, the magnetic moment orientations are updated while the particle positions are fixed. This procedure leads to a new set of data, i.e., and which is used as input data for the subsequent MD simulation. As presented in Section 2.1, the MD algorithm first calculates the net force that acts on every particle according to Equation (1). In a second step, the equations of motion are integrated over a time step according to Equations (6) and (7) leading to the new positions while the magnetic moment orientations are kept unchanged. The new set of data and is used as input data for the subsequent SD simulation. As shown in Figure 4, this process is repeated until a predefined number of total time steps is reached.

The Role of Temperature
In general, both methods MD and SD work in the canonical ensemble, i.e., the temperature is kept constant during the simulation. In order to avoid conflicts between the two heat baths, the temperature of the SD part of the simulation is set to 0 K whereas the MD simulations are performed at room temperature. As explained in [28], it is common practice for the simulation of macromolecules in solution to couple separate thermostats to the subsystems. This practice is based on ab initio methods such as the Car-Parrinello approach [29] where the nuclear and electronic subsystems are separated dynamically. Whereas the slow nuclei are connected to a "physical" temperature, the fast electronic degrees of freedom are assumed to be linked to a "fictitious" temperature. According to [30], the "cold" electronic subsystem is assumed to be close to its instantaneous minimum energy. In analogy to that, here, we assume that the "cold" magnetic degrees of freedom are close to their minimum energy configuration. Hence, it is justified to set the temperature of the SD simulation to 0 . However, as pointed out in Section 2.2. we use a finite spin bath temperature during the SD simulation step in order to find configurations of lowest energy. . Schematic representation of the hybrid algorithm that couples spin dynamics (SD) and molecular dynamics (MD) simulations in order to calculate trajectories of magnetic nanoparticles that interact via magnetic dipole-dipole interaction. Starting from a random initial configuration , the first magnetic low energy state is calculated by means of SD. Therein, the particle positions are retained whereas the orientations of the magnetic moments are changed to . The new orientations together with the fixed particle positions are then passed over to the MD algorithm. Within MD, the magnetic moment orientations are maintained while the particles are moved forward for a small step in time due to the forces that act on the particles. The resulting particle positions are then passed over to SD again in order to calculate the magnetic low energy state of the new configuration. This procedure is repeated until a predefined number of time steps is reached.
Here, slowly cooling down from a high initial temperature to 0 K during the simulation leads in most cases to configurations of lowest energy. If this method is not applicable we use more sophisticated strategies that are based on experimental demagnetization protocols which involve rotating and damped external magnetic fields for finding lowest energy configurations in spin ice systems [27].

Calculation of Qualitative GMR Curves
We use our hybrid SD-MD algorithm in order to predict the particle arrangements in the liquid state of the gel which corresponds to the experimental preparation stage. After the particle-gel mixture is put on a substrate, the structuring process with or without an external magnetic field takes place. Afterwards, the gel matrix is dried out and the particle structure is preserved. For the application in a sensor device, the actual measurement takes places in the solid state of the gel. Therefore, the magnetoresistive properties in the solid state of the gel are of great importance. In order to determine the magnetoresistive properties in a qualitatively way we use the SD algorithm alone which corresponds to the solid state of the gel with mechanically frozen particles. By means of further spin dynamics simulations, magnetization curves can be calculated. Finally, qualitative GMR curves can be calculated from the magnetization curves by: Therein, is the GMR effect given in percent, is the GMR effect amplitude that has to be extracted from experimental data or further numerical investigation [6,31], is the magnetization and represents the material's saturation magnetization. The resulting GMR curves lead to a qualitative estimate of the magnetoresistive properties of a particular particle-gel combination which provides useful information for the actual sensor development.

Application: Magneto-Dynamics of Interacting Magnetic Nanoparticles in Gel Matrices for the Efficient Design of Magnetoresistive Systems
The hybrid SD-MD algorithm can be used for the investigation of the magneto-dynamics of magnetic nanoparticles that are used in combination with conductive gel matrices to build nanoparticular magnetoresistive sensor devices. As presented in [6], this type of sensor can be used for the detection of biomolecules. Therefore, the GMR sensor is covered with antibodies that are chosen for a specific analysis purpose, i.e., the biomolecules under consideration in a sample solution bind to the sensor surface. In order to utilize the GMR effect for the measurement of the biomolecule concentration, magnetic marker particles have to be used. The binding of the marker particles to the biomolecules is realized by using antibodies as well. As a result, the stray field of the magnetic marker particles changes the initial configuration of the magnetic particles and therefore the magnetoresistive properties. Thus, a change in the electrical resistance of the structure can be measured. The sensitivity of such a sensor is related to its magnetoresistive properties which are defined by the microstructure, i.e., the arrangement of the magnetic nanoparticles. Here, we show simulations of the structuring process of an ensemble of magnetic nanoparticles in a viscous matrix. In addition to that, we show how the magnetoresistive properties can be obtained from simulated magnetization curves.

Model System
In general, the hybrid SD-MD algorithm presented in Section 3 can be applied to arbitrarily large systems and is only limited by computational resources. In order to demonstrate its applicability we have modeled a section of a real system made of cobalt nanoparticles [5,7] with a broad size distribution as shown in Figure 5. It can be considered as a small region of a magnetoresistive sensor structure. Here, 50 cobalt particles with diameters in the range between 10 nm and 80 nm are distributed over an area of 350 nm × 350 nm. The particles are immersed in a 4% agarose gel matrix, i.e., we assume a viscous surrounding medium with a kinematic viscosity of = 0.11 Pa•s [7]. The size distribution (b) is taken from experimental studies on particle-based GMR sensors [5,7].

Structuring Process
The structuring process of the cobalt nanoparticle ensemble has been simulated by our hybrid SD-MD method. In order to reach a stable end configuration, a MD run with a total number of 3 × 10 7 time steps and time step length = 1 × 10 −14 (in reduced Lennard-Jones units) has been carried out with one SD run every 10 4 time steps. Every SD run took 1 × 10 4 time steps with a step length of = 1 × 10 −15 s. Due to significant changes of the magnetic moment orientations in the beginning of the simulation, the structural changes are much larger in the beginning and become slower in the end of the simulation. The simulation results are shown in Figure 6. In agreement with previous theoretical and experimental studies [5][6][7], the self-assembly process of the cobalt particles leads to a configuration that consists of parallel aligned areas, antiparallel chains and particle islands that contain magnetic moment vortices.

Estimation of GMR Properties
As presented above, the hybrid SD-MD algorithm considers the particle-gel systems in the liquid states of the gel. By switching off the MD algorithm, the magnetization curve of the system in the solid state of the gel can be calculated by means of spin dynamics simulations alone. The resulting magnetization curves and the GMR calculated using Equation (19) are shown in Figure 7. Furthermore, we show in Figure 7, the corresponding magnetization and GMR curve of the initial particle configuration in order to emphasize the difference of the magnetic properties of the initial configuration and the final configuration that results from the hybrid molecular and spin dynamics trajectory calculation. The GMR effect amplitude is assumed to be 20% in agreement with previous experimental investigations on similar structures [5][6][7].   Figure 6c (black line) and, for comparison, the magnetization curve of the initial configuration shown in Figure 6a (black dots). A rough estimate of the GMR properties of the system has been calculated by means of equation 12 (initial configuration with blue dots and final configuration with blue line) for an effect amplitude of 20% according to previous experimental investigations [5][6][7].
The GMR curves can only be considered a first rough estimate for the magnetoresistive properties of the complete structure since we have only calculated a small section of it. This section contains a few rather large cobalt particles with large magnetic moments that dominate the total magnetization calculation. As a result of this the curves are rather bumpy which would disappear if one simulates larger sections. Nevertheless, even by simulating small sections of a sensor structure different particle materials, size distributions and gel matrices can be investigated in order to optimize parameters for future GMR sensor systems.

Conclusions/Outlook
GMR sensor devices based on magnetic nanoparticles immersed in conductive gel matrices show promising features for the development of disposable low-cost devices for the detection of biomolecules.
In order to optimize the magnetoresistive properties of such sensors it is necessary to know detailed information about the magnetic microstructure. Since this information is not accessible by experimental techniques we provide a novel simulational approach based on hybrid spin dynamics and molecular dynamics methods that allows one to simulate the structuring process of the particles in the liquid state of the gel. Moreover, we have shown that our approach can be used to characterize the GMR properties of a sensor structure. Future enhancements of the hybrid algorithm will include inhomogeneous magnetic fields for the structuring process as well as non-spherical particles such as nanorods and nanocubes.