Nonequilibrium Dynamics of a Magnetic Nanocapsule in a Nematic Liquid Crystal

Colloidal particles in nematic liquid crystals show a beautiful variety of complex phenomena with promising applications. Their dynamical behaviour is determined by topology and interactions with the liquid crystal and external fields. Here, a nematic magnetic nanocapsule reoriented periodically by time-varying magnetic fields is studied using numerical simulations. The approach combines Molecular Dynamics to resolve solute–solvent interactions and Nematic Multiparticle Collision Dynamics to incorporate nematohydrodynamic fields and fluctuations. A Saturn ring defect resulting from homeotropic anchoring conditions surrounds the capsule and rotates together with it. Magnetically induced rotations of the capsule can produce transformations of this topological defect, which changes from a disclination curve to a defect structure extending over the surface of the capsule. Transformations occur for large magnetic fields. At moderate fields, elastic torques prevent changes of the topological defect by tilting the capsule out from the rotation plane of the magnetic field.


Introduction
Nematic liquid crystals are anisotropic fluids usually formed of rodlike molecules that exhibit no positional order but self-organised orientational order [1]. The order parameter S and the so-called director fieldn describe the degree of order and the average molecular orientation of the nematic state, respectively. When doped with colloidal particles (CPs), nematic liquid crystals become enormously interesting due to fundamental and technological reasons [2][3][4]. CPs induce strong deformations onn that depend on the colloidal shape and size and on the anchoring conditions that colloidal surfaces impose on liquid crystal molecules [5]. The interplay between long-range orientational order and superficial alignment inevitably creates topological defects around the CP-regions where molecular ordering is frustrated and S is strongly depressed. The typical radius of these defects is r c 10 nm [6]. However, director deformations extend to cover micrometric scales creating long-range interactions that are exclusive for CPs dispersed in nematic solvents [7]. This promotes self-assembly of CPs and facilitates the fabrication of composite structures with novel designed optical and mechanical properties [8].
To a lesser degree, the effects of driving nematic colloids away from equilibrium have also been investigated. Under specific conditions regarding size, confinement, and anchoring strength, microscopic spheres with homeotropic anchoring, for whichn is perpendicular to the colloidal surface, are surrounded by ringlike topological defects known as Saturnrings [5,21,22]. Flow past these CPs distorts and moves defect rings along the relative flow [23][24][25][26][27][28][29]. In the limit of vanishing inertial contributions, rings around nanospheres move upstream and collapse into a point defect of the hedgehog type [30]. Effects of flow could be used as means for particle selection [24] or self-assembly control [30].
The motion of CPs in nematic liquid crystals can be controlled by micropost arrays [31,32], patterned substrates [33], or undulating walls [34,35]. The latter produce gradients of the director field that attract CPs to particular locations of minimum elastic energy, thus giving rise to a lock-and-key mechanism. As spherical and ellipsoidal CPs move towards the wavy walls, their accompanying Saturn-rings are displaced contrary to the motion of the CPs and suffer a transition to a dipole defect induced by the elastic energy landscape [34]. Morphological changes of topological defects have been also reported for spherocylindrical CPs in a nematic host phase under shear flow [36]. In this case, the relative rotation betweenn and the CP induces hydrodynamic torques on the former, which periodically broaden and extend the accompanying ring defect over the cylindrical portion of the CP. The high elastic energy associated with this transformation is paid by the externally imposed shear flow.
Optical tweezers, provide important additional degrees of freedom to control the motion of microparticles in liquid crystals [37,38], whereas magnetic fields can be used for colloid-molecule coupling and rotatory manipulation [39,40]. A delicately adjusted relation between fluid-imposed, magnetic, and gravitational forces has been used for levitating magnetic nanowires in nematic liquid crystals [41]. Magneto-optic manipulations permit generation of nontrivial topological defects around individual spherical CPs and small colloidal chains [42]. The effect of a rotatory magnetic field on a peanut-shaped CP with a magnetic moment perpendicular to its long symmetry axis has been investigated [16]. The CP rotates as the magnetic field does, while the director field at far distances is unperturbed since the diamagnetic anisotropy of the liquid crystal is very small. The combined effect of magnetic and elastic torques on these CPs is used to estimate their magnetic dipole moment and to control their self-assembly along small chains [16]. Understanding the dynamics of anisotropic nematic colloids under magnetic fields is an important scientific challenge with additional applications in the use of nanoparticles for tuning the Freedericksz transition threshold and orientation relaxation times in liquid crystal cells [43].
Here, magnetically induced rotations of an anisotropic CP in a nematic liquid crystal are studied. The CP is considered to impose homeotropic anchoring conditions. This problem is analysed from a numerical perspective based on a hybrid algorithm combining Molecular Dynamics (MD) and Nematic Multiparticle Collision Dynamics (N-MPCD) [29,36]. This allows us to simulate the nematohydrodynamic behaviour of the solvent and the dynamics of the topological defects surrounding the rotating CP, which are found to be a disclination for bent rings. In addition, the numerical approach allows us to vary relevant parameters, including the frequency and strength of the magnetic field and the order parameter of the solvent, over ranges that permit identification of different rotatory regimes departed from equilibrium. The analysis is restricted to time-varying magnetic fields B(t) that rotate at small and moderate frequencies on a plane containing the director field at far distances from the CPn 0 . At small frequencies, the CP and its defect ring simultaneously follow the rotation of B. At higher frequencies, rotation of the CP is retarded with respect to B by a measurable time lapse. This behaviour is modelled in terms of a stochastic linear equation of the Langevin-type, from which relaxation of the CP towards the instantaneous equilibrium state dictated by B is analysed and a good agreement with numerical results is found. As the capsule rotates at small magnetic fields, it is also inclined out from then 0 − B plane when its long symmetry axis approachesn 0 to avoid large increments of the elastic energy. For large fields, the companion defect changes its form from an elongated Saturn ring to an extended structure that surrounds most of the surface of the capsule. This transformation occurs as an attempt of the solvent to release the elastic energy accumulated, as capsule and defect are driven through nonequilibrium configurations. Defect transformations are present for small rotational Reynolds numbers Re (r) ∼ 10 −2 and small Ericksen numbers Er ∼ 10 −2 , suggesting that they could be experimentally reproduced.

Methods
A magnetic nanocapsule with spherocylindrical shape is placed into a nematic cell of thickness L 3 and homeotropic anchoring conditions, as depicted in Figure 1, where L and σ indicate the length of the cylindrical portion and the radius of the capsule, respectively. Two auxiliary frames are used to describe the system. The first one is a laboratory frame spanned by unit vectors {ê 1 ,ê 2 ,ê 3 }. The second one is attached to the nanocapsule and is spanned by vectors {ê 1 ,ê 2 ,ê 3 }, whereê 3 coincides with the long symmetry axis. Hereafter, the orientation of the capsule will be specified by the angle θ(t) = cos −1 (ê 3 (t) ·ê 3 ). The magnetic moment of the CP is µ = µê 3 (t), and the external magnetic field changes periodically according to B(t) = B[sin(ω 0 t)ê 1 + cos(ω 0 t)ê 3 ]. For illustrative purposes only, small elongated particles are used in Figure 1 to represent the director field within the cell. The capsule promotes homeotropic anchoring on its surface, whereas anchoring on the confining plates imposes the conditionn 0 =ê 3 . It is assumed that the diamagnetic anisotropy of the solvent is very small. Thus, B(t) does not have any effect on the orientation of the director. To study the response of the CP to B(t), numerical simulations are implemented. The employed algorithm has a twofold scope. It integrates Newton's laws of motion of the solvent-solute system over a small time-step and incorporates slower hydrodynamic modes of the solvent by mimicking collective interactions between its microscopic degrees of freedom. This approach has been successfully used to simulate topological defects around spherical and spherocylindrical nematic colloids in equilibrium and nonequilibrium conditions imposed by flow [29,36]. It also reproduces equilibrium states of CPs under static magnetic fields [44]. For self-containment, and considering that the method is rather novel, complete details of the algorithm are given below.

Molecular Dynamics
On the perspective of the mesoscopic method to be used at the large time-scale, the solvent is simulated by N point particles carrying a unit orientation vector. Positions, velocities, and orientations of solvent particles are taken as continuous functions of t and denoted as r i , v i , andû i , respectively, for i = 1, 2, . . . , N. These particles interact with the capsule through a repulsive potential [45], where is the interaction strength and i is the minimum distance from r i to the line segment defining the longitudinal symmetry axis of the cylindrical portion of the capsule. Forces on solvent and solute particles are given by where R is the centre-of-mass position of the CP. The equations of motion of solvent particles are resolved over a small time-step δt according to the MD velocity-Verlet algor ithm [46] where m is the mass of the solvent particles. Similar equations are used to integrate R and the centre of mass velocity of the capsule, whereas its rotational dynamics is integrated in terms of the transformation matrix between the laboratory and moving coordinate systems A and the angular velocity and angular momentum of the CP, Ω and l, respectively. The former satisfiesê α = A ·ê α , and is updated according to the second order expansion It has been shown in reference [47] that time derivatives of A and Ω are related This permits use of the Euler equationṡ to calculate all terms on the right-hand side of Equation (4). In Equations (6)-(8), J 1 , J 2 , and J 3 are the principal moments of inertia of the capsule (J 1 = J 2 ); and κ = A · k. Notice that to sustain the evolution of A, a rule for updating Ω is also required. This is given indirectly by writing Ω in terms of the angular momentum in the laboratory system, Ω = J · A −1 · l, and using a velocity-Verlet rule to update l, namely,

Nematic Multiparticle Collision Dynamics
The mesoscopic behaviour of the solvent will be simulated through a coarse-grained approach that permits reproducing the density, orientation, and flow fields around the CP. This approach is inspired in Multiparticle Collision Dynamics (MPCD), a very useful technique for the simulation of diverse systems in soft condensed matter, including colloids and polymer solutions in isotropic liquids [48,49]. MPCD has been recently generalised to nematic liquid crystals [50][51][52]. The approach used in this paper corresponds to the specific method proposed by Shendruk and Yeomans in reference [51], in which collective collisions between solvent particles take place in localised space regions at regular time intervals of size ∆ t δ t. These collisions are conducted by stochastic collision operators that modify orientations and velocities of the solvent particles as described below.
To apply such operators, the simulation box is subdivided into adjacent cubic regions with volume a 3 , hereafter referred to as collision cells. The number of solvent particles within each collision cell N c (t) is different and changes with time. Particles in the same cell are used to approximate mesoscopic fields such as velocity, order parameter, and director, v c , S c , and n c , respectively. The former is defined as while S c and n c are estimated as the largest eigenvalue and the corresponding eigenvector of the order parameter tensor [45] where I is the identity matrix. Some other auxiliary fields that are also calculated are the centre of mass position the inertia tensor and the velocity and orientation gradients ∇v c and ∇n c , respectively, which are approximated from finite differences of fields v c and n c over the grid of collision cells. N-MPCD considers a mean-field interaction between solvent particles located in the same collision cell. Specifically, the Maier-Saupe potential [53] quantifies the tendency of orientationsû i to align with the director produced by all particles located in the same celln c . Parameter U in Equation (14) tunes the strength of the meanfield potential and is referred to as the nematicity of the solvent. Shendruck and Yeomans have shown that N-MPCD fluids exhibit nematic behaviour for U 5k B T, with the Boltzmann constant k B and temperature T, and that they exhibit an isotropic phase otherwise. The stochastic collision operator for orientations takes new vectorsû i from the canonical distribution In turn, the collision operator for velocities samples new v i vectors from the Maxwell velocity distribution using the Andersen thermostat [48] v where ξ i is the randomly sampled velocity ξ c = ∑ N c j=1 ξ j and ∆ L c is the angular momentum generated by the velocity change. The third and fourth terms on the right-hand side of Equation (15) guarantee conservation of linear and angular momentum, respectively.
An important feature of nematic liquid crystals is that director and flow fields are strongly coupled. Nonuniform flow induces director's reorientation, which might cause hydrodynamic motion known as backflow [1]. N-MPCD cares about this issue through the implementation of rules that recreate the two-way coupling. To simulate reorientation by flow, solvent particles are considered to behave as slender rods with tumbling parameter λ. Orientation changes over the time interval ∆ t are given by Jeffery's equation [54]: where C c and D c are the vorticity and strain rate matrices, i.e., the antisymmetric and symmetric parts of ∇v c , respectively. Moreover, χ HI is a control parameter restricted to the interval [0, 1] that tunes the relaxation time for solvent particles relative to ∆ t. The dynamic response of the director field in the presence of flow is dictated by λ. For λ 1, nematic particles perform a continuous rotation known as tumbling. In contrast, for λ 1, a static alignment under flow is achieved [51].
Backflow is taken into account by translating the angular momentum generated by reorientation of solvent particles within each collision cell ∆ L c ori into linear momentum. With this aim, the reorientation process is considered to be overdamped; thus, ∆ L c ori can be written in the form where γ R is a heuristic viscous rotation coefficient andû j (t + ∆ t) andû j (t) correspond to orientations after and before the N-MPCD collision event, respectively. ∆ L c ori is transferred as part of ∆ L c in Equation (15). It is worth noticing that other versions of N-MPCD resolve backflow differently. In particular, they use a method based on theoretical formulations of nematohydrodynamics [50,52].
To quantify the elastic energy stored in the nematic medium, the Frank free energy density is considered under the one-constant approximation, which has been proven to be valid for the present implementation of N-MPCD [36,51,55]. The total elastic energy of the nematic solvent is where K is the elastic constant and ∑ c indicates summation over all collision cells. By using Equation (18), it is assumed that defect cores are not fully resolved at the coarse-grained scale of N-MPCD. Thus,n c presents no singularities and its spatial derivatives are well-defined.

Boundary and Anchoring Conditions
Confinement is simulated by imposing bounce-back boundary conditions, which completely reverse the velocity of particles encountering solid walls at x 3 = 0 and x 3 = L 3 , while periodic boundary conditions are implemented alongê 1 andê 2 . As it is customary in MPCD methods, the simulation box suffers a random displacement along a vector with components in the range [0, a], where prior collisions take place. This is intended to preserve Galilean invariance in the system [56]. Due to bounce-back boundary conditions, empty spaces appear after this random displacement and virtual particles can be used to fill them. Virtual particles have the same mass and degrees of freedom as solvent particles and are allowed to participate during the momentum exchange event. To this aim, they are inserted randomly within two additional layers at the top and bottom of the simulation box with the same numerical density as the nematic fluid and velocities sampled from the Maxwell distribution at temperature T. The centre of mass velocity in the additional layers is cancelled. These conditions guarantee that momentum flux at the boundary will be correctly simulated during the multiparticle velocity exchange step [57]. Furthermore, virtual particles can be used to control anchoring conditions on the confining surfaces. Orientations of all virtual particles are chosen alongê 3 . This favours the alignment ofn c along the x 3 axis in collision cells containing virtual particles, since they are also allowed to contribute to the orientation exchange operator. Thus, their inclusion promotes effective homeotropic anchoring becauseû i ê 3 for particles close to x 3 = 0 and x 3 = L 3 . Similarly, the solvent experiences homeotropic anchoring conditions over the surface of the magnetic capsule. These are implemented with the aid of two mechanisms. First, nematic particles that satisfy the condition i ≤ 2 1/24 σ during the velocity-Verlet step-see Equation (1)-are reoriented along the perpendicular direction to the surface of the capsule at r i (t). Second, an additional set of massless virtual particles is attached at the surface of the CP, all of which point along the normal direction to its surface. They are randomly distributed over this surface with a superficial numerical density σ S and reorient solvent particles in their vicinity analogously as virtual particles used for the confining walls. It has been shown that these two mechanisms promote strong homeotropic anchoring on the colloidal surface and permit the formation of topological defects around the capsule for equilibrium and nonequilibrium flow conditions [36].

Analytical Description of Capsule Rotation
The analysis in this paper is supplemented by a simplified dynamical model that permits validation of the simulation method. In this model, the magnetic CP is assumed to perform rotations on theê 1 −ê 3 plane while it follows B(t). The magnetic potential is F mag = −µB cos ϑ, where ϑ is the angle betweenê 3 and B(t). The nematic environment imposes an effective potential on the CP, which can be approximated as F nem = F ⊥ + F − F ⊥ cos 2 θ, where F and F ⊥ are configuration energies achieved at parallel and perpendicular orientations of the CP with respect ton 0 , θ = 0 and θ = π/2, respectively. This specific form for F nem was originally proposed by Burylov and Raikher in Reference [58] and explains the tendency of elongated particles to achieve perpendicular orientation with respect ton 0 for F − F ⊥ = ∆ F > 0. The capsule's orientation is considered to obey the equation of motion where η R is the rotational friction coefficient and τ(t) is a stochastic torque. The latter is assumed to be a zero average Gaussian-Markov process with where brackets indicate stochastic average and δ(t − t) denotes Dirac delta function.
As it is shown in Appendix A, in the limit of large magnetic torques, the average stationary rotation of the capsule is retarded from B(t) and slightly deviated from pure sinusoidal dependence according to where the phase retardation is φ θ = η R ω 0 /(µB) and φ ϑ = tan −1 −2ω 0 η R / µB − 4ω 2 0 J 1 .
Thermal noise around the main orientation can be described in terms of ϑ, representing the angle betweenê 3 and the instantaneous magnetic field. In Appendix A, it is shown that the normalised self-correlation function of ϑ in Fourier domain can be approximated as where ω B and β are defined in Appendix A and β 1 = J −1 1 (µB) 2 − η 2 R /4. The model summarised by Equations (21) and (22) show a good agreement with simulations based on the MD-N-MPCD algorithm in the correct limit.

Simulation Parameters
In this paper, results are presented in simulation units using a, m, and k B T as units of length, mass, and energy, respectively (this corresponds to fixing a = 1, k B T = 1, and m = 1 in the simulation code) [49]. All other mechanical quantities have units derived from the previous ones. Table 1 summarises these derived units as well as the parameters used during the simulation stage. In Table 1, t 0 and η 0 indicate units of time and viscosity, respectively;N c is the average numerical density of the solvent (number of N-MPCD particles per collision cell); ρ col is the mass density of the magnetic capsule, which is set five times larger than the solvent's density. Concerning the magnetic field, four different strengths were used. In terms of the magnetic energy, they are expressed as µB = 0, 10, 20, and 50 k B T. Two rotation frequencies were considered, namely, ω 0 = π × 10 −4 t −1 0 and 2π × 10 −4 t −1 0 . Table 1. Units and parameters used for MD-N-MPCD simulations. All these parameters must be given independently to specify a particular numerical execution.

Derived Units Fluid Parameters Colloid Parameters
The size of the simulated capsule could be estimated in terms of the defect core radius, which, for N-MPCD, is found to be r c a [29]; thus, L is scaled to L 100 nm. The shear viscosity of the simulated solvents η can be calculated from the fluid parameters listed in Table 1 and the well known formulae of MPCD analytical descriptions. Specifically [48], which yields η = 12.11η 0 . It is worth noticing that due to their anisotropic character, actual nematic solvents have more than one viscosity coefficient, whereas N-MPCD approximates the momentum transport through the solvent as isotropic [55]. However, this approximation has shown consistent results for simulations of CPs driven by nematic flow [29,36] and can be considered as a first attempt to reproduce the hydrodynamic behaviour on the mesoscopic scale. The selection of tumbling parameter λ = 1.2 will produce solvents with shear alignment behaviour. Notice that in Table 1, two nematicities are given U = 6.5 k B T and 10 k B T, which yield simulations for two distinct solvents with average order parametersS c = 0.60 and 0.74, respectively. The elastic coefficients for these solvents are also different. They can be estimated from order parameter fluctuations, as described in [36,51]. This yields K 151 k B T/a and K 232 k B T/a, for U = 6.5 k B T and U = 10 k B T, respectively.
All simulations were executed for a total time interval 2.5 × 10 4 t 0 , and a thermalisation stage lasting 0.5 × 10 4 t 0 is considered.

Capsule's Rotatory Motion
The simulated CP performs translational and rotational Brownian motion when immersed in the N-MPCD fluid. A topological defect in the form of an elongated Saturn ring appears in the vicinity of the CP, as illustrated in Figure 2a-d, which are four snapshots of the nematic-capsule system taken at regular time intervals after the thermalisation stage of simulations with B = 0. Disclination rings are bent in opposite directions close to the extremes of the CP and are essentially flat when the latter has perpendicular orientation relative ton 0 , θ = π/2. Supplementary Video S1 provides an animation of the motion of the capsule at B = 0. In the absence of magnetic field, configurations close to perpendicular are visited more frequently than others. This indicates thatê 3 ⊥n 0 is the most stable state for the capsule and ∆ F > 0. To estimate ∆ F precisely, the numerical probability distribution of θ is obtained and compared with the canonical distribution P(θ) sin θ dθ = exp −∆ F cos 2 θ /k B T sin θ dθ/Z, where Z is the partition function. By considering ∆ F as an adjustable parameter, values ∆ F = 1.75 k B T and 2.5 k B T are obtained for U = 6.5 k B T and 10 k B T, respectively. As expected, the capsule's orientations is closer to theê 1 −ê 2 plane at larger nematicity. Figure 2e shows the comparison between the numerical and canonical orientation distributions. The capsule's rotations are induced when the magnetic field is applied. These take place on theê 1 −ê 3 plane, though deviations outside this plane can be observed and become appreciable mainly in the case µB = 10 k B T. The dynamic orientation of the nanocapsule is represented through the function cos θ(t) in Figure 3, where the noisy trajectory indicates the actual time series observed in simulations with ω 0 = 2π × 10 −4 t −1 0 and the continuous curve is the approximation cos(θ stat (t)), obtained from the model in Section 2.4. For comparative purposes, the function B 3 (t)/B = cos(ω 0 t) has been included in Figure 3 as a dashed curve representing the rotation of B. To calculate the stationary solution θ stat (t) in Equation (21), it is considered that the rotational friction coefficient of the capsule is similar to that presented by a long rod of length L and radius σ [59]: where b is a factor that can be adjusted to correct specific differences in slip boundary conditions and shape. It has been estimated to be b 0.6 [36,44]. Figure 3 exhibits a good agreement between simulations and theory. For large magnetic fields (µB = 50 k B T, in Figure 3a), the model predicts no appreciable difference between θ stat (t) and ω 0 t. Correspondingly, the capsule follows the magnetic field very closely, performing small thermal fluctuations around the main trajectory. The phase retardation between cos(θ stat (t)) and B 3 (t)/B predicted by the model becomes noticeable for moderate fields (µB = 20 k B T, Figure 3a). This retardation is present in simulations and can be appreciated since, in Figure 3b, the noisy trajectory is closer to the continuous curve than to the dashed curve. Notice also that fluctuations in the case µB = 20 k B T are stronger than those in the case µB = 50 k B T. These effects are more pronounced for µB = 10 k B T (Figure 3c), where retardation in simulations is larger than that predicted analytically. This can be explained by noticing that stronger thermal fluctuation reduce the average projection ofê 3 on theê 1 −ê 3 plane. Thus, the effective magnetic torque alongê 2 decreases from the value assumed in the two-dimensional linear model, µB, which, in turn, produces an increment of φ θ in Equation (21). Fluctuations around the average trajectory are analysed in terms of the normalised correlation χ(ω) defined by Equation (22). This correlation is calculated in simulations using the resulting time-series ϑ(t) assuming stationarity and applying a simple filtering process based on a moving average. Results for the specific case of ω 0 = π × 10 −4 t −1 0 are compared with predictions of the model in Figure 4. Numerical correlations are wellapproximated by the analytical curves for the three considered magnetic fields, which cover underdamped (µB = 50 k B T) and overdamped dynamics (µB = 10 k B T and µB = 20 k B T). The agreement is better in the case µB = 50 k B T (Figure 4c), while in cases µB = 10 k B T and µB = 20 k B T (Figure 4a,b, respectively), the analytical model slightly underestimates the actual width of the autocorrelation function. In addition, systematic oscillations can be identified at the tails of χ(ω) in the simulation case. Both of these effects are due to the fact that the stochastic process ϑ(t) is nonstationary but its autocorrelation has been estimated under the stationary approximation. Nevertheless, the good correspondence presented in Figure 4 provides confidence in using the proposed MD-N-MPCD method to simulate the nonequilibrium stochastic dynamics of the magnetic capsule.

Elastic Energy Changes and Topological Defect Deformations
The rotational Reynolds number and the Ericksen number are two dimensionless quantities that can be used to characterise the rotation of the capsule in the nematic liquid crystal. The former is defined as the ratio of the inertial torque to the viscous torque on the CP [60] and can be estimated as where ρ is the solvent mass density. Then, it can be noticed that simulations were conducted with Re (r) = 0.05 and Re (r) = 0.1 for ω 0 = π × 10 −4 t −1 0 and ω 0 = 2π × 10 −4 t −1 0 , respectively. Er is the ratio of viscous-to-elastic effects and is customarily defined as where γ 1 is a rotational viscosity and v char (l char ) is a characteristic velocity (length). The ratio K/γ 1 can be estimated from the spectrum of director fluctuations as discussed in References [55,61]. Using the same set of MPCD parameters as here but with U = 20 k B T, it is obtained K/γ 1 0.62a 2 t −1 0 [61]. However, this estimation is not expected to vary noticeably over a wide range of U values [55]. On the other hand, l char and v char can be taken from mapping the angular velocity of the CP into translational velocity, i.e., l char ∼ L and v char ∼ Lω 0 /2π. Thus, Er ∼ 0.032 and Er ∼ 0.016 are estimated for simulations with ω 0 = π × 10 −4 t −1 0 and ω 0 = 2π × 10 −4 t −1 0 , respectively. This analysis permits situating simulations in the regime of small Reynolds and Ericksen numbers, where viscous effects dominate over inertial effects and elastic forces dominate over viscous forces.
It is observed thatn is periodically distorted as the capsule rotates due to the resulting moving anchoring imposed on its surface. The rotatory motion of the magnetic CP is accompanied by a corresponding rotation of the disclination ring. These effects can be studied in terms of time-dependent elastic energy F elas (t), which will hereafter be simplified to the normalised form where F eq elas is the equilibrium average elastic energy at B = 0. The behaviour of f elas (t) is illustrated in Figures 5 and 6 for solvents with U = 6.5 k B T and U = 10 k B T, respectively. For each nematicity, four results are presented for the combination of parameters B = 10 k B T, B = 50 k B T, ω 0 = π × 10 −4 t −1 0 , and ω 0 = 2π × 10 −4 t −1 0 . Energy exhibits periodic increments. They appear with the frequency of the rotational motion of the capsule. These increments are, in general, stronger for larger magnetic fields and at some particular events, they have the shape of a remarkable peak that could be as large as 4% the elastic energy of the solvent in the simulation box. A comparison between Figures 5 and 6 suggests that these large energy increments are more common for U = 6.5 k B T. By inspecting the time evolution of the capsule, it can be concluded that f elas increments appear whenê 3 gets close to ±n 0 . At these states, noticeable transformations of the topological defect can take place. This process is illustrated in Figure 7 for the specific case of parameters U = 10 k B T, ω 0 = 2π × 10 −4 t −1 0 , and µB = 50 k B T. Figure 7a-f are snapshots of the capsule-defect couple as it traverses the energy peak located in the time interval t ∈ (14, 250 t 0 , 15, 750 t 0 ) in Figure 5. The CP rotates, and when its long symmetry axis coincides withê 3 , the form of the companion defect changes from a disclination ring to a broader structure that almost completely surrounds the cylindrical portion of the CP (Figure 7d,e). This transformation is reverted as capsule's orientation departs from the director field at far distances. Supplementary Video S2 provides an animation of the capsule's rotation and the topological defect changes for the same parameters used in Figure 7. Similar changes have been observed recently in defects around elongated particles driven by an external shear flow [36]. Here, their occurrence can be explained in terms of the orientation angle θ as follows. When θ = 0, the elongated, bent Saturn ring defect has been proven to be the structure with the smallest elastic energy [11,36,62]. However, for θ 0, the most stable structure is a small ring located close to one of the ends of nanocapsule, while the disclination ring is a metastable configuration of larger elastic energy [11,36]. Due to the rotating magnetic field, the capsule-defect pair is driven to the state θ 0, possessing large elastic energy. Then, the system attempts to reduce this energy by changing the shape of the defect from the unstable large ring to the small ring at the top or bottom of the CP. During this process, the region of small order spreads around the capsule. However, the transition to the small ring is frustrated because the continuously changing field removes the anisotropic CP from θ 0, and the stability of the elongated Saturn ring is recovered. Interestingly, Figure 6c suggests that morphological transformations of the disclination rings can be eluded since elastic energy increments are kept moderate during the whole simulation. Indeed, no such transformations were observed for that case, corresponding to U = 10 k B T, µB = 10 k B T, and ω 0 = 2π × 10 −4 t −1 0 . The reason is that, for small magnetic fields and large nematicity, it is expected that the torque exerted by the solvent on the CP will be able to act against the magnetic torque. Then, as the capsule-defect pair approaches the state θ = 0 and the nematic torque tilts the CP out of theê 1 −ê 3 plane. This permits preserving the shape of the topological defect around the capsule. Figure 8a,b illustrate, respectively, the trajectories of the orientation vectorê 3 over the unit sphere obtained from simulations at (U = 10 k B T, µB = 10 k B T) and (U = 10 k B T, µB = 50 k B T). In the first case, the inclination about theê 1 −ê 3 plane is significant, whereas in the second one, the strength of the magnetic field maintains rotations of the capsule very close to that plane. By comparing trajectories in Figure 8a,b with the corresponding energy changes, Figure 7c,d, respectively, it can be concluded that the tilt induced by elastic torques prevents large energy peaks associated with defect shape transformations. This also explains why such peaks appear more frequently in the case U = 6.5 k B T ( Figure 5) than in the case U = 10 k B T ( Figure 6).

Discussion and Conclusions
Techniques from MD and N-MPCD were combined into a single algorithm with multiscale characteristicsthat recreated the dynamics of individualmagnetic nanocapsules immersed in nematic liquid crystals and subjected to external fields. The method shows that time-varying magnetic fields can be used to effectively control the rotation of the nanocapsules in the liquid crystal along with the topological defects around them. Capsules rotate at the pace of the magnetic field when the latter rotates at small frequencies. A phase retardation between magnetic field and capsule rotations could be observed, which was in agreement with predictions of a simplified dynamical model. The power spectrum of capsule orientation fluctuations along the main rotating trajectory was quantified and found to be in agreement with theoretical approximations. This encourages us to propose that MD-N-MPCD can be used in complex nonequilibrium situations involving the interaction of liquid crystals with moving objects. Consequently, it could be particularly well-suited for simulating particles in anisotropic fluids where external forces cause translational and rotational motion of the immersed particles. It could also be used to study the self-assembly process of anisotropic CPs or the dynamics of active CPs in liquid crystals.
During the rotation of the anisotropic CP, the associated topological defect suffers morphological changes and a consequential elastic energy gain is observed in the solvent. These changes appear as a response of the solvent to the action of the magnetic field that drives the capsule-defect pair to an orientational state of high elastic energy. Then, a reorientation process of the director field around the CP occurs that intends to transform the disclination ring into a defect of smaller elastic energy. However, such transformation is prevented by the magnetic field itself that drives the CP back to weak nonequilibrium states where the disclination defect is more convenient again. This nonequilibrium effect is similar to that observed for anisotropic particles rotating in a sheared nematic [36]. However, it should be noticed that energy gains sustained by magnetically induced colloidal rotations avoid global nematic flow and, therefore, they could be easier to produce in practice. Furthermore, the aforementioned morphological changes are obtained for rotational Reynolds numbers Re (r) ∼ 0.05 and Ericksen numbers Er ∼ 10 −2 , which are smaller than those considered in the nonequilibrium situation induced by shear. This suggests that they could be experimentally accessible.
When magnetic torques were moderate compared with torques exerted by the nematic solvent, morphological transformations of defects were not observed. This is because the elastic torque on the capsule counteracts the magnetic torque and prevents the system from being driven far from equilibrium. This result is in very good qualitative agreement with recent experimental observations in which peanut-shaped dipolar CPs were rotated by magnetic fields and tilted out of the rotation plane by the nematic liquid crystal [16]. All this reinforces the reliability of MD-N-MPCD as a promising method to deal with nonequilibrium dynamics in liquid crystal colloids.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

CP
Colloidal particle MD Molecular Dynamics N-MPCD Nematic Multiparticle Collision Dynamics