Frequency-dependent magnetic susceptibility of magnetic nanoparticles in a polymer solution: a simulation study

Abstract

Magnetic composite materials i.e. elastomers, polymer gels, or polymer solutions with embedded magnetic nanoparticles are useful for many technical and bio-medical applications. However, the microscopic details of the coupling mechanisms between the magnetic properties of the particles and the mechanical properties of the (visco)elastic polymer matrix remain unresolved. Here we study the response of a single-domain spherical magnetic nanoparticle that is suspended in a polymer solution to alternating magnetic fields. As interactions we consider only excluded volume interactions with the polymers and hydrodynamic interactions mediated through the solvent. The AC susceptibility spectra are calculated using a linear response Green–Kubo approach, and the influences of changing polymer concentration and polymer length are investigated. Our data is compared to recent measurements of the AC susceptibility for a typical magnetic composite system [Roeben et al., Colloid Polym. Sci., 2014, 2013–2023], and demonstrates the importance of hydrodynamic coupling in such systems.

---

1 Introduction

Magnetic nanoparticles (MNPs) have been connected to many technical applications and are also promising candidates for novel bio-medical treatments. They can be manufactured by splitting larger magnetic objects into smaller particles (e.g. by grinding) or synthesized using chemical processes, allowing for particles of well-defined shape, size, and material properties.^1–9 As the human body mass is transparent to static magnetic fields, MNPs that are injected into a body can be moved by application of inhomogeneous external fields. This allows an externally controlled accumulation of MNPs in specific places, e.g. in tumors. Loading the particles with therapeutics before injection, they can be used as a targeted drug delivery system – an important step toward locally constrained chemotherapy.^10–14 Another promising approach to cancer treatment based on MNPs is the local heating of cancer cells, so-called ‘hyperthermia’.^15–20 Typical bio-medical applications of MNPs will be situated in a polymer-rich, i.e., complex cellular environment. The use of MNPs in a technical context has mainly two aspects. On one hand, they can be embedded in a polymeric gel or suspension to create a soft composite material. The coupling between the magnetic properties of the particles and (visco)elastic properties of the polymer matrix allows for the dynamic control of some aspects of the composite. This includes, e.g., changing their shape, motion, or elastic properties using external fields, which makes them interesting candidates for a variety of applications.^21–28 On the other hand, MNPs can be viewed as probes that are used to asses the local behaviour of a polymeric environment. A well established experimental technique for that is the so-called AC susceptometry, where a (sinusoidal modulated) external magnetic field is applied to a sample. Frequency-dependent susceptibility spectra are obtained by simultaneous measurement of the probes' magnetization response. Besides the insight gained into the magnetic behaviour of the composite as a whole, there exist a number of theoretical models that connect the susceptibility to frequency-dependent mechanical properties. Using, e.g., extended Debye models,^29–31 the Gemant–DiMarzio–Bishop model,^32–37 or the theoretical model derived by Raikher et al.,³⁸ one may infer the frequency-dependent local viscosity and even elastic moduli on the particle scale, a method dubbed magnetic nanorheology.

A good understanding of the coupling mechanisms is key, when it comes to tailoring the properties of MNPs to a specific bio-medical and technical use case. Unfortunately, often the details of the coupling between the particles and the polymeric environment are not fully known and moreover difficult to (directly) determine experimentally. Here, computer simulations are a powerful tool, allowing the study of well-defined model systems. Aspects of the computational models can be easily modified, where it might be difficult to do so experimentally. Computational models exist on different scales of modelling. While continuum descriptions are well-suited to model MNP–polymer composites on a macro-level,^39–43 they typically lack details on the length scale of single polymers. This, however, will be necessary to study how individual changes to the polymer matrix affect its coupling to the embedded MNPs. In this work, we therefore present a particle-based molecular dynamics (MD) simulation model for MNPs immersed in a polymer solution. As explicit resolution of the polymers comes at a high computational cost – in particular at high concentrations and large sample sizes – we focus on the efficient modelling of such systems. The model, which will be presented in Section 2, uses coarse-grained polymers, that are coupled to a computationally efficient lattice-Boltzmann hydrodynamics solver; MNPs are modelled as so-called ‘raspberry particles’,^44,45 in order to couple translations as well as rotations hydrodynamically to the environment. We use our model to simulate a typical nanorheological system. As a representative example, we use the system of Roeben et al.⁴ throughout this publication. The parameters required for simulation using our model are specified in Section 3. To compare to the experiment, we need to obtain AC susceptibility spectra from the simulations, which can be done in two ways. While the susceptibility can be directly measured for specific frequencies (Section 4.2), our main method is based on Green–Kubo linear response theory (Section 4.3). It allows us to obtain a complete susceptibility spectrum from a single (however more expensive) simulation. With this approach, we demonstrate that our model reproduces the experimental trends. Using the Green–Kubo approach, we run simulations to illustrate how changes to the polymer matrix affect the resulting susceptibility. Higher values for both, the polymer volume fraction and the polymers' chain length, are found to independently cause a shift of the susceptibility spectra towards lower frequencies (Section 5), which is in accordance with the experimental observations of Roeben et al.⁴ We end with a discussion of our model and the results obtained our the simulations of the experimental system.

2 Computational model

In this section we present an efficient MD model for simulations of typical nanorheological systems, i.e. (magnetically blocked) MNPs immersed in polymer suspension. Such systems consist of three main components that have to be included in the model, namely the polymers, the fluid, and the magnetic particles.

In MD simulations, parts of the system are represented by interacting MD particles, which are propagated through space and time by time-discrete integration of Newton's equation of motion. For the polymers, the so-called bead-spring formalism is used. Every polymer chain is explicitly resolved using MD beads, which are inter-connected by an attractive potential – here, we use a simple harmonic interaction of the form

where r is the distance between two consecutive MD beads within a polymer chain, r₀ is the potential's equilibrium distance and the spring constant k determines the stiffness of the bond. The MD beads themselves interact via a purely repulsive Weeks–Chandler–Andersen (WCA) potential⁴⁶based on the well-known 6-12 Lennard-Jones (LJ) pair potential,with ε the depth of the LJ potential well and σ the distance at which φ_LJ = 0. The repulsive interaction effectively models excluded volume. Note that typically every MD bead within the polymer will represent more than one ‘real-world’ monomer unit. This coarse-graining approach allows for significantly larger system sizes or simulation times as it reduces the number of interacting particles.

To include the hydrodynamic coupling between MNPs and close-by polymers in our model, hydrodynamics are solved using the lattice-Boltzmann (LB) algorithm.^47,48 This very efficient grid-based hydrodynamics solver can be straight-forwardly parallelized, has good scaling behaviour, and there also exist fast implementations using graphics cards. Their computational performance make LB algorithms well-suited for the mesoscale. On the small length-scales of our system, thermal fluctuations are an important part of modelling. The LB fluid can be thermalized by adding stochastic fluctuations to the stress tensor, which guarantees local mass and momentum conservation.^49–51 To couple MD particles to this fluid, a dissipative friction force is used, which also has to be thermalized. Following Ahlrichs and Dünweg,⁵² the coupling force is given by

Its strength is determined by the friction parameter γ and the difference of fluid velocity u⃑_fluid and velocity of the particle u⃑_part. Fluctuations are incorporated via the stochastic force of zero mean andwhere δ is the Kronecker delta.

In typical experimental setups, the MNPs are at least one order of magnitude larger than the polymers' persistence length λ_p (e.g., in the experimental system which serves as a basis for the parameter choices in Section 3,⁴r_h/λ_p ≈ 18.9, with hydrodynamic radius r_h of the magnetic particle). To properly capture the hydrodynamic behaviour of the MNPs in the model, it does not suffice to just use a single coupling point per particle. With such an approach every MNP would couple to only a single lattice site per time step, with no rotational coupling. In addition, the hydrodynamic impact of the MNP's shape would be neglected. Thus, to obtain the physically expected hydrodynamic behaviour, the so-called ‘raspberry model’ is used for the MNPs, which means homogeneously filling their shape with fluid coupling points that are rigidly bound to a central MD bead.⁴⁴ For preparation of the raspberry we follow the procedure described by Fischer et al.,⁴⁵ which we will only briefly sketch here. For in-depth information please consult the referenced publication. Creating the raspberry particle happens in a two-step process: resolving the surface and then homogeneously filling its volume, which is necessary to obtain consistent translational and rotational diffusivities D_t and D_r.^45,53 The preparation process uses a MD approach. To resolve the surface of the raspberry, (mobile) MD beads are put into an empty simulation box. A harmonic potential (see eqn (5)) that is shifted by the desired raspberry radius r₀ = r_rasp, is used to force the beads onto a spherical shell. To homogenize the surface density, a purely repulsive WCA potential (see eqn (2)) is added between the MD beads. By slowly increasing the harmonic spring constant k to an extremely high target value, MD beads end up forming a homogeneous surface shell. After correcting small deviations of the beads' desired distances from the raspberry centre (|r⃑| = r_rasp), they are rigidly connected to the centre position and the harmonic bonds removed. Now the raspberry can be filled by adding MD beads inside the hollow shell, which interact again via a repulsive WCA potential. To prevent numerical instabilities, the resulting forces are initially capped. The positions of the MD beads are allowed to evolve using MD (with a Langevin thermostat) while slowly raising the force cap. This procedure typically leads to a homogeneous distribution of MD beads, which is tested by tracking the shell-wise particle density, as well as the geometric centre of all raspberry bead positions (it should coincide with the raspberry centre). Once these conditions are fulfilled, the MD beads are rigidly bound to the raspberry centre.⁴⁵ All beads of the resulting raspberry-like shaped rigid body now act as fluid coupling points. They couple to the LB fluid via the same friction force as the other MD beads (see eqn (4)). The equation of motion, however, is only integrated for the body's central bead, which carries the mass M_rasp and moment of inertia I = 2M_raspr_rasp²/5 of the particle, as well as its magnetic moment m. Fig. 1 shows a cut through the raspberry particle used for this publication. Imposing a meaningful potential between the central bead and polymer MD particles, one can model respective interactions – as we will see in Section 3, we choose a purely repulsive WCA potential, here, but more involved interaction potentials would also be possible. Note that for a spherical particle these interactions can only lead to translations. Coupling of the rotational behaviour happens only via hydrodynamics which are thus essential to our model.

Fig. 1 Cut through the spherical raspberry particle used for this publication. It consists of homogeneously distributed fluid coupling points (blue) forming a rigid body. Newton's equation of motion is integrated only for its centre (black), which carries the particle's mass, as well as its magnetic moment (indicated by the yellow arrow).

3 Model parameters

To establish a close relationship between experiment and modelling, we set up simulations based on the experimental system used by Roeben et al.⁴ and Hess et al.⁵⁴ for nanorheological measurements. It consists of spherical, magnetically blocked single-domain cobalt ferrite (CoFe₂O₄) nanoparticles immersed in an aqueous polyethylene glycol (PEG) solution. The CoFe₂O₄ particles act as probes for investigating the nanorheological properties of the polymer solution. Their magnetic moment is small, and they are used at very high dilutions (3.6 × 10⁻²% volume fraction⁴). Therefore, dipolar interactions between the MNPs do not play a major role (for an estimate of their relative strength see the calculation of λ in Section 5). Consequently, we neglect dipolar interactions in our model. This allows us to simulate only a portion of the system containing a single magnetic particle. Periodic boundary conditions are used to eliminate boundary effects, but with a very large simulation box to minimize hydrodynamic interactions between the MNP and its periodic images.

All simulations for this publication were performed using the MD software package ESPResSo.⁵⁵ This software uses the concept of simulation units, that allows the user to freely choose the mass, length, and energy scales of the system, which in combination determine the simulated time scale. For reference, the values of important quantities describing the system are summarized in Table 1 using both, SI and simulation units, respectively.

Table 1 Overview of the parameters we use for the simulation unit system, for some characteristic system quantities

Quantity	Value (SI units)	Value (simulation units)
Particle radius rrasp	7.2 nm	4 [x]
Thermal energy	k B·300 K	1 [E]
Solvent density ρw	1 × 10^3 kg m^−3	1 [m]/[x]^3
Kinematic viscosity ν	8.9 × 10^−8 m^2 s^−1	1.86 [x]^2/[t]
Mass unit	5.8 × 10^−24 kg	1 [m]
Time unit	67.4 ps	1 [t]

---

The length scale in our model is prescribed by the hydrodynamic radius of the raspberry particle r_rasp, which corresponds to the experimental value of r_h = 7.2 nm.⁴ As simulation domain, we use a cubical box with side length L_box = 10r_rasp. The energy scale is given by the thermal energy of the system. As the experiments are performed at T_exp = 25 °C, we use T = 300 K in our model, prescribing the energy scale 1 [E] = k_B·300 K.

For the fluid, we make the following choices: the size of the LB grid cells is set to a_grid = 1 [x] = r_rasp/4 yielding [x] = 1.8 nm; LB fluid density is set to 1 [m]/[x]³ corresponding to water density ρ_w = 1 × 10³ kg m⁻³, thereby also defining the mass scale [m] = 5.8 × 10⁻²⁴ kg. For the magnetic particle, we use the bulk density of CoFe₂O₄, ρ_CoFe2O4 = 5.3 × 10³ kg m⁻³.⁵⁶ We make a simplifying assumption concerning the fluid viscosity by choosing its value as one tenth the experimental kinematic viscosity of water, ν_exp = 0.1ν_w = 8.9 × 10⁻⁸ m² s⁻¹,⁵⁷ in simulation units ν = 1.86 [x]²/[t]. This effectively decreases the decay time of the magnetic moment of the MNP through rotational motion, which reduces simulation time by roughly one order of magnitude. Per every MD time step Δt = 0.01 [t], the LB fluid field is updated. The simulated time scale is fully determined by the prescribed length, energy, and mass scales yielding [t] = 67.4 ps.

In the experiment, PEG is the polymer species being used. Experimentally, the mean monomer length in PEG is 2.78 Å,⁵⁸ the chain's persistence length in aqueous solution λ_p ≈ 3.8 Å.^59,60 To accurately capture the experimental size ratio between magnetic particles and persistence length of PEG would be computationally exceedingly demanding, limiting investigations to small polymer concentrations or simulation times. Our aim is, however, to end up with a computationally model that can handle polymer volume fractions up to as high as ∼25%. In addition, simulations of such systems have to capture both, the relaxation time of the polymer suspension and the relatively long self-diffusion time of the magnetic particle, making coarse-graining of the polymers indispensable. While this may change the strength of some effects in the system, general trends are expected to remain unaffected – for results and discussion, please see Section 4. As water is a good solvent for PEG, excluded volume has to be considered for the polymers, which is accounted for using a purely repulsive WCA interaction between monomer beads (see eqn (2)). The diameter of the monomer beads is given by the choice of the WCA parameter, here σ_mono = 1 [x]; we set ε_mono = 1 [E]. The equilibrium length of the harmonic potential which connects monomer beads within each polymer chain (eqn (5)) is set to r₀ = σ_mono, the spring constant k = 200 [E]/[x]². Finally, to prevent polymers from penetrating the ‘solid’ MNP, an additional WCA potential is required between the monomer beads of the polymers and the raspberry particle. In principle, this interaction could be set between monomer beads and the coupling points of the raspberry. Since, in general, the coupling points do not smoothly resolve the raspberry surface, any remaining surface roughness introduces additional rotational coupling. As detailed in Section 2, using the raspberry model the main concern is to obtain consistent translational and rotational diffusive behaviour for the MNP, which can be challenging in itself and does in no way guarantee a smooth resolution of the MNP's surface. To avoid artefacts in the rotational coupling, we therefore prescribe the repulsive WCA potential between the monomer beads and the centre or the raspberry particle. The parameters for this mixed interaction are ε_mix = 1 [E] and σ_mix = (r_rasp + σ_mono/2). The raspberry used for this publication consists of 1400 coupling points (≈5.25 per LB unit cell), the fluid-coupling friction parameter is γ = 15 [m]/[t]. A visualization of the raspberry is given in Fig. 1, a representative snapshot of the readily set-up system is shown in Fig. 2.

Fig. 2 Visualization of the simulation model described in Section 2 using the parameters of Section 3. Raspberry coupling points are colored blue, (coarse-grained) polymers green. Here, each polymer chain consists of 15 beads. For visibility reasons, the polymer volume fraction in the shown setup is only ϕ = 0.01. Typical volume fractions used for our simulations are significantly higher, up to ϕ = 0.25.

4 Obtaining AC susceptibility spectra from simulations

The AC susceptibility is the differential response of a system's magnetization M(t) to an oscillating external magnetic field H(t), i.e. χ = dM(t)/dH(t). The result is a complex value χ = χ′ − iχ′′.

4.1 For reference: Debye model

Before introducing the methods for measuring AC susceptibility spectra in the simulation, we briefly present the result of the theoretical Debye model.⁶¹ It was first derived by Debye (1929), originally to describe the relaxation of the rotational polarization of molecules, but generally captures the rotational relaxation process of any dipole. The model makes a few assumptions. Namely, it describes the magnetization of non-interacting ideal dipoles that couple to a fluid with frequency-independent dynamic viscosity η = νρ via rotational Stokes friction ζ_r = 8πηr³. The magnetic susceptibility analytically evaluates towith χ₀ the zero-frequency susceptibility and τ the characteristic decay time of the magnetic moment.

For a raspberry particle in a LB fluid, as used in our model, these assumptions are fulfilled. We can therefore use the Debye model as a reference curve to benchmark our measurement techniques against. In our model we consider magnetically blocked particles only and thus do not allow internal (Néel-)relaxation of the magnetic moment. Instead, decorrelation of the orientation of the dipole moment only happens via (rotational) Brownian diffusion. As a result we have

with Boltzmann's constant k_B, temperature T, and the hydrodynamic volume V_h = 4πr_h³/3.⁶¹

4.2 Direct measurement

The straightforward approach of measuring the AC susceptibility of the simulation is by faithfully reproducing the experimental measurement technique i.e. measuring per frequency f: for any chosen angular frequency ω = 2πf, a time-dependent alternating magnetic field

H⃑(t) = H⃑₀ sin(ωt), (8)

with amplitude H⃑₀ is explicitly applied to the simulation box, causing a torque on the magnetic particle. The magnetization response of the system is measured via the component of the particle's magnetic moment that points into field direction.

For weak coupling to the external field, the magnetization M⃑(t) will also have a sinusoidal shape

M⃑(t) = M⃑_max sin(ωt − ϑ), (9)

but with a frequency-dependent phase shift ϑ. However, this signal will be superimposed by thermal noise. To bring out the signal, averaging over several periods is necessary, as will be discussed shortly. The phase shift is then used to obtain the complex magnetic susceptibility χ = M(t)/H(t) = χ′ − iχ′′, where the real part χ′, that relates to reversible magnetization, is in phase with the external field. The imaginary component χ′′ is related to losses and can be computed as the out-of-phase part of the magnetization response. We getA visualization of the analysis is shown in Fig. 3. Exemplary results of such direct measurements for a MNP in a Newtonian fluid are shown in Fig. 4. For reference, the corresponding Debye model result is also shown. The positions of imaginary susceptibility peaks match quite well. The same is true for the real part of the susceptibility in that region. There are, however, some deviations from the curve predicted by Debye. This can be attributed to the non-linear nature of the Langevin magnetization law.^62,63 Typically, the so-called Langevin parameter is used to characterize the magnetic interaction strength. The Langevin parametercompares the magnetic interaction energy between the particle's dipole m and the field to the thermal energy. Here, μ₀ is the magnetic constant. In order to minimize non-linear effects in the susceptibility spectra, one would like to choose α ≪ 1. However, choosing a low α simultaneously increases the relative importance of thermal fluctuations in the system – the signal-to-noise ratio drops significantly, which means much longer simulation times are necessary for comparable statistics. Minimizing non-linear effects when using direct measurements thus comes at a high computational cost. Another obvious shortcoming of this technique is that each simulation yields the complex susceptibility only for the chosen frequency ω. To sample complete susceptibility spectra with high resolution is therefore hardly feasible. However, direct measurements may come in handy when observing interesting behaviour at specific frequencies.

Fig. 3 Direct measurement of the susceptibility: an external magnetic AC field H (frequency f = 1 291 549 Hz) is applied to the system and its magnetization response M(t) is measured. Due to Brownian diffusion this signal will be noisy and has thus to be averaged over several periods, 64 here. The mean magnetization (orange) is then fitted (green, for functional form see eqn (9)) to obtain the phase shift ϑ.

Fig. 4 AC susceptibility spectra measured from the simulation of a MNP in a Newtonian fluid without any polymers. Blue data points are results of direct measurements (Langevin parameter α = 2). For reference, the theoretical Debye model is shown (green). The spectrum obtained using the Green–Kubo approach (red) matches the theoretical model quite accurately.

4.3 Green–Kubo approach

In this section, we discuss an approach that allows one to obtain the full AC susceptibility spectrum from a single simulation, assuming that the applied external magnetic field is small. In general, the dynamic response of a magnetic suspension to an applied magnetic field is non-linear. However, for small fields, α ≪ 1, linear response can be assumed. In the framework of linear response theory, the so-called Green–Kubo relations are a general type of fluctuation–dissipation theorems that connect physical quantities such as polarization or complex conductivity to time correlation functions of associated dynamic variables.⁶⁴ Also for the magnetic susceptibility one such Green–Kubo relation can be derived. It readswhere V is the system's volume. The integral part of this equation is the Fourier–Laplace transform of the time. Applying this relation, a full AC susceptibility spectrum can be directly obtained from one single steady-state simulation (without an explicit external magnetic field).

As this approach solely relies on the auto-correlation of the magnetic moment as an input, one has to make sure that this quantity is well sampled even for large lag times. As a consequence, simulation runs have to be significantly longer than with the direct method. At this point, the simplifications in the used parameter set (decreased viscosity and water density for the magnetic particle) significantly decrease computation time by speeding up the decorrelation of the magnetic moment. Still, in practice, dealing with the integral part can be challenging, as thermal noise in the long-time tail of the auto-correlation function of the magnetic moment is amplified by the frequency. For infinite runs noise in the long-time tail would cancel, leaving a pure decay signal. A common approach to mitigate the impact of tail noise is thus to fit an analytic expression for the tail – a typical choice would be a simple exponential decay,

with 〈M⃑(0)M⃑(0)〉_fit the zero-frequency limit of the magnetization correlation and characteristic decay time τ_fit. Note that depending on the system, more involved functional descriptions may be required to properly capture the decay behaviour, such as the expression for an isolated particle in a viscoelastic medium derived by Ilg and Evangelopoulos.³⁴ Deviations from a pure exponential decay behaviour are also expected for interacting and polydisperse ferrofluid systems.⁶⁶ For our setup, the decay behaviour of the long-time tail is actually well captured using the simple exponential decay. For the well-sampled low-frequency region the numerical data is kept, whereas the fitted function is used for the long-time tail. Between these two regions a smooth transition is made over a certain frequency range. The merged data is given by

〈M⃑(t)M⃑(0)〉 = f_t(t)〈M⃑(t)M⃑(0)〉_num + [1 − f_t(t)]f_fit(t), (15)

where f_t(t) is a linear transition function over the interval [τ_l,τ_u],For a visualization of this procedure see Fig. 5. Here, we decided to choose the same time range for fitting and the linear transition between numerical data and analytical tail. As lower boundary (τ_l), we use the viscous time of our system τ_r² = ρ_wr_rasp²/η = 0.58 ns (the time it takes fluid momentum to diffuse by one colloidal radius).^45,65 As upper limit (τ_u), we choose the Brownian relaxation time of the Newtonian system τ_B = 100.8 ns, see eqn (7).

Fig. 5 Representative example for the time auto-correlation function of the magnetization obtained from our simulation model (blue solid line). Here, a MNP in a Newtonian fluid was used. To get rid of the Brownian noise in the long-time tail, an analytic exponential function is fitted to the data, the fit region is the indicated grey area. The same frequency range is used to linearly transition from the numerical correlation data to the analytical tail fit (see eqn (15) and (16)). The resulting correlation function is shown by the orange dashed line.

The integral in eqn (13) can then be efficiently computed via the Fourier–Laplace-transform of resulting auto-correlation function of the magnetization. The result of such a Green–Kubo susceptibility measurement on a magnetic particle in a Newtonian fluid is shown in Fig. 4. It perfectly matches the expected theoretical Debye relaxation reference curve.

5 Results and discussion

Using our modelling approach (Section 2) with the parameters of Section 3 – based on the experimental system of Roeben et al.⁴ – we set up simulations to study the signatures that changes to the polymeric surroundings generate in the measured susceptibility spectra. Here, we will focus on the effects of (1) varying the polymer volume fraction and (2) the polymer chain length.

To obtain full AC susceptibility spectra, we use the Green–Kubo approach described in Section 4.3, which by default operates in the linear response limit. The linear response assumption also holds in the experiment as we can see from the estimated experimental Langevin parameter. With field amplitude H_0,exp = 0.4 kA m⁻¹,⁴ temperature T_exp = 25 °C, and magnetic moment of the particle m = 3.6 × 10⁻¹⁹ A m²,⁵⁴ we obtain

which confirms that the experiments do indeed take place in the linear response limit. For magnetically blocked particles the magnetic moment scales with the particles' core volume, m ∼ r_rasp³ (assuming constant thickness of surface coating if present). As the used value was actually measured for a CoFe₂O₄ particle with hydrodynamic radius r_h = 10.7 nm (instead of our r_h = 7.2 nm),⁵⁴ we even overestimate the Langevin parameter α_exp, here.

With these parameters, we also estimate the relative strength of inter-particle dipolar interactions using the dipolar interaction parameter

which compares the per-particle dipolar energy for two touching particles in head-to-tail configuration to the thermal energy. We obtain λ_exp ≈ 0.25, a value that indicates rather small dipolar interactions (for particles to show notable structure, one would expect a λ larger than ∼2). Particle–particle interactions become important for sufficiently high densities of MNPs and strength of dipolar interactions.^66–69 Here, however, the low dipolar interaction parameter and volume fraction (<1%) imply that particle–particle interactions do not significantly influence the susceptibility. This justifies our approach of studying a single magnetic particle in its surrounding.

We start with simulations of varying polymer volume fraction. The polymer volume fraction refers only to the polymer solution, i.e., it ignores the volume of the magnetic particles. In general, the polymer volume fraction is thus given by

with number of polymers N_poly, monomer units per polymer N, and the radius per monomer unit σ (to calculate its excluded volume). The box volume is given by V_box, N_rasp is the number of raspberries with radius r_rasp inside the box. For the parameters used in our simulations, see Section 3. To study the effect of increasing polymer concentrations, we choose a constant chain length of N = 15 and systematically vary the volume fraction ϕ ∈ {0.05, 0.15, 0.25}. The magnetic moment of the particle is recorded every 0.1 time unit (every 6.7 ps). The total number of samples per simulations run is N_samples = 2 × 10⁷. The resulting susceptibility spectra are shown in Fig. 6. With increasing polymer fraction, the susceptibility curves experience a shift towards lower frequencies. As the Debye shape of the curve is preserved for all cases, these shifts signify an increase of the effective viscosity around the MNP. Even at the highest concentration, no entanglement of polymers is expected – as will be discussed shortly, both polymer length and volume fraction would have to be considerably higher.

Fig. 6 The AC susceptibility spectra obtained from simulations with polymers of chain length N = 15 at different volume fractions. Increasing the polymer volume fraction leads to a shift towards lower frequencies. Solid lines indicate the real, dashed lines the imaginary part of the susceptibility.

Let us now turn to the simulations at a constant volume fraction ϕ = 0.25 where the polymer chain length is systematically varied N ∈ {1, 5, 15, 25}. Results are shown in Fig. 7. Again, a shift of the measured spectra towards lower frequencies is observed while maintaining the Debye shape, where larger shifts occur with increasing chain length.

Fig. 7 AC susceptibility spectra for magnetic particles in solution with volume fraction of polymers ϕ = 0.25. The chain length of the polymers N is varied. The longer the chains are, the more the spectra get shifted towards lower frequencies. Solid lines indicate the real, dashed lines the imaginary part of the susceptibility. The inset shows the corresponding fitted relaxation times τ_fit.

It can be seen from the figure, that already the addition of unconnected monomers causes a significant shift in frequency. In our computational model, the monomer beads are coupled to the fluid via the thermalized friction force of eqn (4), thus slowing down fluid dynamics. This registers as an increased effective fluid viscosity, causing a higher relaxation time for the systems' magnetization. An overview of the relaxation times obtained via fitting the long-time tail (Section 4.3) can be found in Table 2. By increasing both, the length of polymers and their concentration, this purely hydrodynamic effect is intensified through the increased coupling to the polymer matrix. The observed behaviour matches the experimental observations of Roeben et al.⁴ In the experiment, PEG solutions with different chain lengths and polymer concentrations were studied. As in our simulation results, for low polymer fractions (depending on the chain length) the experimentally measured susceptibility spectra are simple Debye-type curves with a single relaxation time, which directly depends on the effective viscosity of the polymer suspension. With increasing PEG concentration and chain length, however, the experimentally measured spectra broaden and some measurements also suggest the existence of a second relaxation process.⁴ The emerging deviation from Newtonian behaviour is most noticeable above the overlap concentration of the respective solution. This indicates entanglement of polymer chains as a main cause. The broadening of Debye spectra for spherical MNPs in increasingly elastic environments was also observed by other experimental groups.⁵

Table 2 Overview of the zero-frequency magnetization correlations and decay times fitted using the procedure described in Section 4.3, for simulations at different polymer volume fraction ϕ and number of monomers per chain N

ϕ	N	〈M⃑(0)M⃑(0)〉fit	τ fit/τB
0.00	0	1.00	1.010
0.05	15	0.99	1.089
0.15	15	0.99	1.185
0.25	1	0.99	1.216
0.25	5	0.99	1.240
0.25	15	0.98	1.333
0.25	25	0.98	1.373

---

To get a feeling for the polymer lengths required to observe entanglement in the simulation model, the work of Kremer and Grest⁷⁰ on polymer melts may serve as an indication. Using a FENE potential (instead of harmonic) to connect the MD beads within polymers and a slightly different parameter set, they found an entanglement length of N_e = 35 at a volume fraction of ϕ = 0.45.⁷⁰ In our simulations we are well below these values, with maximum N = 25 and ϕ = 0.25. Despite the deviations in modelling, we thus expect that polymer chains would have to be significantly longer and used at higher volume fractions before seeing effects of entanglement.

As already described, our model makes several simplifying assumptions such as the use of coarse-grained polymers and using one tenth the viscosity of water for our fluid. Our results are therefore not expected to quantitatively reproduce the experimental findings. To get an idea of the relative strength of the effects in the experiment and simulation, let us make a rough comparison, neglecting polymer entanglement effects. In the simulations we found that the use of polymers of length N = 25 at a polymer volume fraction of ϕ = 0.25 leads to an increase of the relaxation time by ∼37%. Experimentally, the use of 25 vol% PEG increased the dynamic viscosity and thus the relaxation time (see eqn (7), we are in the Debye-like Newtonian limit) by roughly a factor of 2 for triethylene glycol and about a factor of 3 for PEG1000.⁴ The effect observed in the simulations is thus roughly one order of magnitude smaller than in the experiment, possibly due to the small length of our polymers (compared to PEG1000), and the neglect of any specific interactions of the polymers with the MNPs.

As a key result, our simulations show the importance of hydrodynamics for magnetic composite systems. Since all pair interactions between the MNP and surrounding monomer beads just affect the particle's translational behaviour, hydrodynamic interactions are the only possible source of rotational coupling. Our simulations demonstrate that this is already sufficient to qualitatively reproduce the experimentally observed shifts of the AC susceptibility spectra towards lower frequencies for both, higher polymer volume fractions and longer chains, respectively.

6 Conclusions

In this work, we presented a simulation model for magnetic nanoparticles in polymer suspensions and showed how frequency-dependent AC susceptibility spectra can be obtained from the simulations. Our modelling approach is based on molecular dynamics, coupled to an efficient lattice-Boltzmann hydrodynamics solver. While still costly, we managed to bring the computational complexity down to a feasible level. To do so, we made some simplifying assumptions, such as a certain level of coarse-graining regarding the polymers. Using our model, we ran simulations with parameters based on the experimental system of Roeben et al.,⁴ to study the signatures that changes to the polymeric surroundings produce in the resulting susceptibility spectra. For this, we simulated a single particle in its polymeric environment (using periodic boundary conditions). This is justified because of the low concentration of particles as well as small dipolar interaction strength in the experimental system. We presented two approaches to obtain susceptibility spectra from our simulations. Using direct measurement one may obtain the susceptibility for a specific frequency, which, however, is quite costly when it comes to sampling entire spectra over a large range of frequencies. We therefore introduced a second approach based on Green–Kubo linear response theory. To benchmark this method with our model, we showed that the simulation of a particle in a Newtonian fluid reproduces the curve expected from Debye relaxation theory. We used the Green–Kubo approach throughout our simulation study of a spherical nanoparticle in a polymer suspension. Matching the experimental observations,⁴ we found that the AC susceptibility spectra are shifted towards lower frequencies,⁴ when either increasing the polymer volume fraction or the polymer chain length for a fixed polymer concentration. In our model, hydrodynamic interactions provide the sole source of direct coupling between the polymeric environment and the rotational behaviour of the particle, i.e. the Brownian relaxation of the systems' magnetic moment. Finding that this is already sufficient to qualitatively reproduce the experimental trends highlights the essential role of hydrodynamics for the behaviour of spherical magnetic nanoparticles in complex environments.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through the priority program SPP 1681 (Grant HO 1108/23-3).

Notes and references

1. S. Laurent, D. Forge, M. Port, A. Roch, C. Robic, L. Vander Elst and R. N. Muller, Chem. Rev., 2008, 108, 2064–2110.
2. J. A. López Pérez, M. A. López Quintela, J. Mira, J. Rivas and S. W. Charles, J. Phys. Chem. B, 1997, 101, 8045–8047.
3. M. Faraji, Y. Yamini and M. Rezaee, J. Iran. Chem. Soc., 2010, 7, 1–37.
4. E. Roeben, L. Roeder, S. Teusch, M. Effertz, U. K. Deiters and A. M. Schmidt, Colloid Polym. Sci., 2014, 292, 2013–2023.
5. H. Remmer, E. Roeben, A. M. Schmidt, M. Schilling and F. Ludwig, J. Magn. Magn. Mater., 2017, 427, 331–335.
6. S. Sun, H. Zeng, D. B. Robinson, S. Raoux, P. M. Rice, S. X. Wang and G. Li, J. Am. Chem. Soc., 2004, 126, 273–279.
7. R. Sen, P. Jain, R. Patidar, S. Srivastava, R. Rana and N. Gupta, Mater. Today: Proc., 2015, 2, 3750–3757.
8. A. Günther, P. Bender, A. Tschöpe and R. Birringer, J. Phys.: Condens. Matter, 2011, 23, 325103.
9. C. Märkert, B. Fischer and J. Wagner, J. Appl. Crystallogr., 2011, 44, 441–447.
10. R. Tietze, J. Zaloga, H. Unterweger, S. Lyer, R. P. Friedrich, C. Janko, M. Pöttler, S. Dürr and C. Alexiou, Biochem. Biophys. Res. Commun., 2015, 468, 463–470.
11. R. Tietze, S. Lyer, S. Dürr, T. Struffert, T. Engelhorn, M. Schwarz, E. Eckert, T. Göen, S. Vasylyev and W. Peukert, et al. , Nanomedicine, 2013, 9, 961–971.
12. S. Goodwin, C. Peterson, C. Hoh and C. Bittner, J. Magn. Magn. Mater., 1999, 194, 132–139.
13. C. Alexiou, R. Jurgons, R. Schmid, A. Hilpert, C. Bergemann, F. Parak and H. Iro, J. Magn. Magn. Mater., 2005, 293, 389–393.
14. C. Alexiou, R. Jurgons, C. Seliger, O. Brunke, H. Iro and S. Odenbach, Anticancer Res., 2007, 27, 2019–2022.
15. E. A. Périgo, G. Hemery, O. Sandre, D. Ortega, E. Garaio, F. Plazaola and F. J. Teran, Appl. Phys. Rev., 2015, 2, 041302.
16. U. M. Engelmann, A. A. Roeth, D. Eberbeck, E. M. Buhl, U. P. Neumann, T. Schmitz-Rode and I. Slabu, Sci. Rep., 2018, 8, 13210.
17. R. Hergt, S. Dutz, R. Müller and M. Zeisberger, J. Phys.: Condens. Matter, 2006, 18, S2919.
18. A. Aqil, S. Vasseur, E. Duguet, C. Passirani, J. P. Benoît, R. Jérôme and C. Jérôme, J. Mater. Chem., 2008, 18, 3352–3360.
19. L. Lao and R. Ramanujan, J. Mater. Sci.: Mater. Med., 2004, 15, 1061–1064.
20. H. S. Huang and J. F. Hainfeld, Int. J. Nanomed., 2013, 8, 2521.
21. C. Gollwitzer, A. Turanov, M. Krekhova, G. Lattermann, I. Rehberg and R. Richter, J. Chem. Phys., 2008, 128, 164709.
22. D. Szabó, G. Szeghy and M. Zrínyi, Macromolecules, 1998, 31, 6541–6548.
23. R. Ramanujan and L. Lao, Smart Mater. Struct., 2006, 15, 952.
24. S. Odenbach, Arch. Appl. Mech., 2016, 86, 1–11.
25. T. Mitsumata and S. Ohori, Polym. Chem., 2011, 2, 1063–1067.
26. S. Monz, A. Tschöpe and R. Birringer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 021404.
27. Z. Varga, G. Filipcsei and M. Zrínyi, Polymer, 2006, 47, 227–233.
28. T. I. Volkova, V. Böhm, T. Kaufhold, J. Popp, F. Becker, D. Y. Borin, G. V. Stepanov and K. Zimmermann, J. Magn. Magn. Mater., 2017, 431, 262–265.
29. F. Ludwig, A. Guillaume, M. Schilling, N. Frickel and A. Schmidt, J. Appl. Phys., 2010, 108, 033918.
30. V. L. Calero-DdelC, D. I. Santiago-Quiñonez and C. Rinaldi, Soft Matter, 2011, 7, 4497–4503.
31. M. Feyen, E. Heim, F. Ludwig and A. M. Schmidt, Chem. Mater., 2008, 20, 2942–2948.
32. E. A. DiMarzio and M. Bishop, J. Chem. Phys., 1974, 60, 3802–3811.
33. A. Gemant, Trans. Faraday Soc., 1935, 31, 1582–1590.
34. P. Ilg and A. E. Evangelopoulos, Phys. Rev. E, 2018, 97, 032610.
35. S. Havriliak Jr. and S. J. Havriliak, J. Polym. Sci., Part B: Polym. Phys., 1995, 33, 2245–2252.
36. R. Diaz-Calleja, A. Garcia-Bernabe, M. Sanchis and L. Del Castillo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 051505.
37. K. Niss, B. Jakobsen and N. B. Olsen, J. Chem. Phys., 2005, 123, 234510.
38. Y. L. Raikher, V. Rusakov, W. Coffey and Y. P. Kalmykov, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 031402.
39. O. V. Stolbov, Y. L. Raikher and M. Balasoiu, Soft Matter, 2011, 7, 8484–8487.
40. K. A. Kalina, P. Metsch and M. Kästner, Int. J. Solids Struct., 2016, 102, 286–296.
41. P. Metsch, K. A. Kalina, C. Spieler and M. Kästner, Comput. Mater. Sci., 2016, 124, 364–374.
42. A. M. Menzel, Soft Matter, 2017, 13, 3373–3384.
43. M. Puljiz and A. M. Menzel, Phys. Rev. E, 2019, 99, 053002.
44. V. Lobaskin and B. Dünweg, New J. Phys., 2004, 6, 54.
45. L. P. Fischer, T. Peter, C. Holm and J. de Graaf, J. Chem. Phys., 2015, 143, 084107.
46. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237.
47. G. R. McNamara and G. Zanetti, Phys. Rev. Lett., 1988, 61, 2332–2335.
48. T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer, Cham, 2017.
49. A. J. C. Ladd, J. Fluid Mech., 1994, 271, 285–309.
50. B. Dünweg and A. J. Ladd, Advanced Computer Simulation Approaches for Soft Matter Sciences III, Springer, 2009, pp. 89–166.
51. L. D. Landau and E. M. Lifshitz, Course of theoretical physics. vol. 6: Fluid mechanics, Pergamon Press, London, 1959.
52. P. Ahlrichs and B. Dünweg, J. Chem. Phys., 1999, 111, 8225–8239.
53. S. T. Ollila, C. J. Smith, T. Ala-Nissila and C. Denniston, Multiscale Model. Simul., 2013, 11, 213–243.
54. M. Hess, E. Roeben, P. Rochels, M. Zylla, S. Webers, H. Wende and A. M. Schmidt, Phys. Chem. Chem. Phys., 2019, 21, 26525–26539.
55. F. Weik, R. Weeber, K. Szuttor, K. Breitsprecher, J. de Graaf, M. Kuron, J. Landsgesell, H. Menke, D. Sean and C. Holm, Eur. Phys. J.-Spec. Top., 2019, 227, 1789–1816.
56. A. Shinde, Int. J. Innov. Technol. Exp. Eng., 2013, 3, 64–67.
57. J. Kestin, M. Sokolov and W. A. Wakeham, J. Phys. Chem. Ref. Data, 1978, 7, 941–948.
58. F. Oesterhelt, M. Rief and H. Gaub, New J. Phys., 1999, 1, 6.1–6.11.
59. J. E. Mark and P. J. Flory, J. Am. Chem. Soc., 1965, 87, 1415–1423.
60. F. Kienberger, V. P. Pastushenko, G. Kada, H. J. Gruber, C. Riener, H. Schindler and P. Hinterdorfer, Single Mol., 2000, 1, 123–128.
61. P. J. W. Debye, Polar molecules, Chemical Catalog Company, Incorporated, 1929.
62. S. Odenbach, Colloids Surf., A, 2003, 217, 171–178.
63. R. Weeber, M. Hermes, A. M. Schmidt and C. Holm, J. Phys.: Condens. Matter, 2018, 30, 063002.
64. R. Kubo, J. Phys. Soc. Jpn., 1957, 12, 570–586.
65. R. Zwanzig and M. Bixon, J. Fluid Mech., 1975, 69, 21–25.
66. A. O. Ivanov and P. J. Camp, Phys. Rev. E, 2018, 98, 050602.
67. J. O. Sindt, P. J. Camp, S. S. Kantorovich, E. A. Elfimova and A. O. Ivanov, Phys. Rev. E, 2016, 93, 063117.
68. A. O. Ivanov and O. B. Kuznetsova, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 041405.
69. V. S. Mendelev and A. O. Ivanov, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 051502.
70. K. Kremer and G. S. Grest, J. Chem. Phys., 1990, 92, 5057–5086.

---