Effect of applied magnetic fields on the morphology of nematic nanobridges in slit pores

Molecular Physics An International Journal at the Interface Between Chemistry and Physics ISSN: (Print) (Online) Journal homepage: www.tandfonline.com/journals/tmph20 Effect of applied magnetic fields on the morphology of nematic nanobridges in slit pores Pablo Romero-Llorente & José Manuel Romero-Enrique To cite this article: Pablo Romero-Llorente & José Manuel Romero-Enrique (19 Nov 2024): Effect of applied magnetic fields on the morphology of nematic nanobridges in slit pores, Molecular Physics, DOI: 10.1080/00268976.2024.2428818 To link to this article: https://doi.org/10.1080/00268976.2024.2428818 Published online: 19 Nov 2024. Submit your article to this journal Article views: 16 View related articles View Crossmark data Full Terms & Conditions of access and use can be found at https://www.tandfonline.com/action/journalInformation?journalCode=tmph20 MOLECULAR PHYSICS e2428818 https://doi.org/10.1080/00268976.2024.2428818 RULL-ABASCAL SPECIAL ISSUE FOR STATISTICAL PHYSICS IN SPAIN (BY INVITATION ONLY) Effect of applied magnetic fields on the morphology of nematic nanobridges in slit pores Pablo Romero-Llorenteaand José Manuel Romero-Enrique b,c aDepartamento de Sistemas Físicos, Químicos y Naturales, Universidad Pablo de Olavide, Seville, Spain; bDepartamento de Física Atómica, Molecular y Nuclear, Área de Física Teórica, Universidad de Sevilla, Seville, Spain; cInstituto Carlos I de Física Teórica y Computacional, Campus Universitario Fuentenueva, Granada, Spain ABSTRACT In this paper, we report a molecular dynamics study of the effect of the application of magnetic fields on the morphology of nematic nanobridges of 32,000 oblate Gay-Berne particles in slit pores favouring homeotropic anchoring. In absence of magnetic fields, previous studies show that there are different conformations of the nanobridge, depending on the slit pore width Dand the wettability of the walls. Under the application of uniform magnetic field cycles, in which the intensity of the field is increased stepwise until a maximum value and then decreased at the same rate, switching between different nanobridge conformations can be observed if the magnetic field is applied in a perpendicular direction to the global nematic director of the nanobridge. However, there are situations in which the initial bridge conformation is recovered after the magnetic field application, indicating that the switching can only be observed if there are different locally stable nanobridge conformations under the same thermodynamic conditions. Moreover, magnetic fields can destabilise the nanobridge, leading to its breakdown into isolated nanodroplets attached to a single wall. ARTICLE HISTORY Received 1 August 2024 Accepted 5 November 2024 KEYWORDS Nematic liquid crystals; confinement; Gay–Berne model; magnetic field; molecular dynamic simulations 1. Introduction Nematic liquid crystals are fluids in which long-range orientational order is present [1]. However, interplay between elastic distortions, surface tension and anchoring can easily frustrate the orientational order of nematic liquid crystals in the presence of boundary surfaces, such as walls or interfaces, leading in general to the emergence of topological defects [2], which are regions where the orientational and/or positional order abruptly vanishes and alter the long-range order around them, generating structures which stable against perturbations. In nematics, the most frequent topological defects are point defects or hedgehogs, and disclination lines, where the defect cores are either points or lines, respectively. Both types of defects are characterised by their winding number or topological charge, which accounts for the number CONTACT José Manuel Romero-Enrique [email protected] Departamento de Física Atómica, Molecular y Nuclear, Área de Física Teórica, Universidad de Sevilla, Avenida de Reina Mercedes s/n, Seville 41012, Spain; Instituto Carlos I de Física Teórica y Computacional, Campus Universitario Fuentenueva, Calle Dr. Severo Ochoa, Granada 18071, Spain of times the nematic director field rotates in a counterclockwise way if we follow a closed continuous path which encloses the topological defect core. For hedgehogs, it is usually either +1or−1, where the minus sign means that the nematic director rotates in a clockwise way. On the other hand, for disclination lines it is either +1/2 or −1/2, which means that the nematic director rotates πradians, either counterclockwise or clockwise, respectively. The manipulation of the orientational ordering by external electric or magnetic fields is known since the seminal work by Fréedericksz and Repiewa [3], which is the basis of the nowadays liquid crystal displays technology. In many cases, there are locally stable conformations of the bounded or confined nematics which differ by the orientational ordering of the liquid crystal, such as zenithal bistable nematic devices © 2024 Informa UK Limited, trading as Taylor & Francis Group 2P. ROMERO-LLORENTE AND J. M. ROMERO-ENRIQUE with different surface relief structures [4–7]. The application of electric or magnetic fields can switch between these conformations which, in many cases, have different optical properties. Due to this fact, these systems have been proposed as key ingredients for practical applications such as privacy windows and other electro-optic devices, in which liquid crystal droplets dispersed in a polymer matrix change their optical properties under the application of external fields [8–14]. We will address the study of capillary bridges confined in slit capillaries. For this purpose, the Gay–Berne model for nematogens [15] is considered. This model has been used to model realistic molecular systems. This potential is parametrised by κ, which accounts for the anisotropy in the repulsive part of the potential, κfor the anisotropy intheattractivepartofthepotential,andµand ν,which determine the orientational dependence of the potential depth. A more precise definition of these parameters will begiveninthenextsection.Forexample,theGay–Berne model with (κ,κ,μ,ν) =(4.4, 20, 1, 1)has been used to parametrise p-terphenyl molecular interactions [16, 17] (which do not show experimentally liquid crystalline behaviour), and (0.345, 0.2, 1, 2)for the nematogen triphenylene core discotic model [18,19]. The bulk phase diagram of this model has been extensively studied for both calamitic [20–28] and discotic [18,29–32]cases.Herewe consider an oblate Gay-Berne (0.5, 1, 2, 1),whichshows homeotropic anchoring on the vapour-nematic interface, i.e. the direction associated to the shortest molecular axis orients preferentially parallel to the normal to the interface [33,34]. In addition, the wall-particle interactions are chosen to promote homeotropic anchoring on the wall. This choice is motivated by the experimental work reported in [35] about 5CB capillary bridges on airinaslitporeformedbyparallelglasssurfaceswith homeotropic anchoring conditions,which reports different nematic configurations in the capillary bridge by tuning its diameter-to-height aspect ratio and shape, which can be either barrel-like or hourglass-like. So, for large aspect-ratio bridges a ring defect located on the equatorial plane is observed, while a point-defect structure is observed for small aspect-ratio bridges. On the other hand, the defects are either hyperbolic for hourglass-like bridges or radial for barrel-like bridges. However, the results obtained from Monte Carlo simulations for our oblate Gay–Berne model [36] disagree for small aspect ratios with experiments. So, instead of the point-defect structure, simulations show biaxial bridge textures in which an annular disclination is arranged in a vertical plane. This can be qualitatively rationalised by the different length scales involving simulations and experiments [36]. In this paper, we will study how the application of magnetic fields can affect the nanobridge conformations. Forthispurpose,wewillstartfromthecasesreported in Ref. [36]and,ineachcase,wewillexplorehowthe application of an uniform magnetic field with intensity modulated gradually in time can affect the morphology of the nanobridge. We will restrict ourselves to situations in which the magnetic field is applied to a direction perpendicular to the global nematic director. The reason forthisisthat,forlargeaspect-ratiobridgestheglobal nematic director is perpendicular to the walls, i.e. the z axis, but for small aspect ratios is on the transversal xy plane. Thus, switching between configurations may be expected if the magnetic field is applied in the xdirection for large aspect-ratio bridges and in the zdirection for small aspect-ratio bridges. The paper is organised as follows. The potential model and computer simulation methods will be presented in Section 2. Section 3will describe the principal findings of our simulations. We end up the paper with a discussion of our results and main conclusions. 2. Methodology 2.1. Potential model Particles are subject to a pairwise potential energy, the socalled Gay–Berne potential [15], which is a modification of the Lennard–Jones potential that has been extensively considered in the study of nematogen fluids. It takes the form Uij(rij,ui,uj)=4ε(ˆ rij,ui,uj)ρ−12 ij −ρ−6 ij ,(1) where ρij =rij −σ(ˆ rij,ui,uj)+σ0 σ0 .(2) Here uiistheunitvectoralongthesymmetryaxisofparticle i,rij =|ri−rj|is the distance along the intermolecular vector rij joining the centres of mass of particles i and j,andˆ rij =rij/rij. The anisotropic contact distance, σ(ˆ rij,ui,uj), and the depth of the interaction energy, ε(ˆ rij,ui,uj), depend on the orientational unit vectors, the length-to-breadth ratio of the particle, κ=σee/σss, and the energy depth anisotropy, κ=ee/ss which are both defined as the ratio of the size and energy interaction parameters in the end-to-end (ee) and side-by-side (ss) configurations, respectively. With these definitions, κ>1 corresponds to prolate particles while κ<1cor- responds to oblate particles. Their expressions are given MOLECULAR PHYSICS 3 by σ(ˆ rij,ui,uj) σ0=1−χ 2(ˆ rij ·ui+ˆ rij ·uj)2 1+χ(ui·uj) +(ˆ rij ·ui−ˆ rij ·uj)2 1−χ(ui·uj)−1/2 (3) and ε(ˆ rij,ui,uj) 0=[1(ui,uj)]ν×[2(ˆ rij,ui,uj)]μ,(4) where 1(ui,uj)=[1 −χ2(ui·uj)2]−1/2,(5) 2(ˆ rij,ui,uj)=1−χ 2(ˆ rij ·ui+ˆ rij ·uj)2 1+χ(ui·uj) +(ˆ rij ·ui−ˆ rij ·uj)2 1−χ(ui·uj),(6) χ=(κ2−1)/(κ2+1)and χ=[(κ)1/μ −1]/[(κ)1/μ +1]. Here σ0is the side-by-side intermolecular collision diameter, and 0is [2κ/(κ2+1)]νtimes the minimum intermolecular potential energy between two molecules in the side-by-side configuration. We will consider the original Gay and Berne choice for μ=2andν=1[15]. The length-to-breadth geometrical ratio κand κplays an important role in the formation of ordered phases, and also determine the nematic phase anchoring on the nematic-vapour phase [33,37]. Thus, the nematic anchors homeotropically to the nematic-vapour interface if κ≤κ, being random-planar otherwise. In addition to the fluid-fluid interactions, particles are subject to two external potentials. The first one is the interaction with the slit pore of width Dalong the zaxis Uw,i(zi)=UGB,w(zi)+UGB,w(D−zi),(7) where ziis the height of the particle iabove the wall. The single-wall-particle potential Uw,iisa9−3potential: UGB,w(zi)=a2π 3ρwσ3 ww2 15 σw zi9 −σw zi3. (8) Here ρwσ3=0.988, σw/σ0=1.096, w/0=1.277, which are the parameters for the wall-particle interaction in the Ar–CO2system [38]. Finally, ais an adimensional wall interaction strength which can be varied. Although his potential does not depend on the orientation of the GB particles (unlike in Ref. [39]), this wall promotes homeotropic anchoring to the nematic phase [36]. Finally,theinteractionofeachparticlewithanuniform magnetic field B=BuBis given by UB i=λ 21−3(ui·uB)2,(9) where λis the magnetic interaction strength parameter, which has the expression [40] λ=χB2 3μ0 , (10) with uibeing the symmetry axis orientation, χ the magnetic polarisability and μ0the vacuum magnetic permeability. χ the anisotropy in the molecular magnetic polarisability, which is usually positive for calamitic nematogens and negative for discotic nematogens [41], and μ0thevacuummagneticpermeability.Despitein this paper a model for discotic liquid crystals is considered, we assume that χ is positive. This potential favours particle to orient along the magnetic field or opposite to it by inducing a torque TB igiven by TB i=3λ(ui·uB)(ui×uB). (11) 2.2. Simulation procedure We carry our molecular dynamics simulations of a system of N=32,000 GB oblate particles with κ=0.5 and κ= 1 by using the LAMMPS open source package [42–44] in an in-house-modified version to include the interaction between particles and the magnetic field. We used LAMMPS implementation of the GB potential [45,46], which is unshifted, and we choose a cutoff radius of 4σ0. Energies, lengths and times are given in units of 0,σ0 and σ0√m/0,respectively,wheremis the mass of the particle. Accordingly, the reduced temperature T∗and number density ρ∗are defined as T∗=kBT 0 ,ρ∗=ρσ3 0. (12) Hereafter we will drop the asterisk, so quantities are given in reduced units. Previous Monte Carlo simulations [33,34] showed that nematic-vapour coexistence appears in a narrow range of temperatures T=0.4−0.5. We choose T=0.5, in which the density of the nematic phase at coexistence was estimated as ρ=1.86(2),with a nematic order parameter S=0.50 [33,34]. In addition, the nematic-vapour interface shows homeotropic anchoring [34], as expected for κ<κ . The simulation box is a square cuboid of dimensions L×L×D,withL=45 and different values of D.Periodic boundary conditions were imposed in the xand y directions,whiletwowallsareplacedatz=0andz=D. 4P. ROMERO-LLORENTE AND J. M. ROMERO-ENRIQUE In this paper, we consider slit pores with D=20, 30 and 40, similarly to the cases considered in Ref. [36], and adimensional wall interactions strengths a=1 and 1.25. Simulations are performed in the NV T ensemble with the Nose-Hover thermostat [47,48] for a fixed temperature T=0.5. The time step is set to t=0.002 and the thermostat relaxation time to τ=0.01. The magnetic field application cycle is composed by two stages. In the first stage, we perform simulations in which the parameter λ,Equation(10),isstepwise increased by λ =0.025 each 2 ×105steps from a null value up to a maximum value λmax. Afterwards, during the second stage the values of λare stepwise decreased again by λ =0.025 each 2 ×105steps until the zeromagnetic field case is reached. This gradual modulation in time of the magnetic field intensity is chosen to allow the nanobridge to equilibrate in most stages of the magnetic field cycle. For each value of D,thestartingconfiguration is the bridge configuration obtained from the Monte Carlo simulations for a=1, with velocities and angular velocities chosen randomly following the corresponding equilibrium Maxwell–Boltzmann distribution. As we wish to induce a switch between different orientational textures by the magnetic field application, magnetic fields favour different orientations from the initial ones. Thus, for D=20, in which the initial condition is axisymmetric about the zaxis, magnetic fields are applied along the xaxis with a maximum value λmax =0.1. For this value, the bridge-averaged nematic order parameter Sis slightly above the bulk value, and most particles arealignedtothemagneticfieldthroughoutthebridge. On the other hand, for D=30 and D=40, in which a disclination line on a vertical plane is observed in the initial configuration, magnetic fields are applied along the zaxis with a maximum value λmax =0.15. This maximum value was selected in a similar way to the case of a magnetic field applied along the xaxis. For each value of λ, different physical quantities, such as the mean potential energy and the global orientational ordering, are monitored during the simulation. The instantaneous global nematic ordering is characterised by the tensor order parameter in each step: Q=1 N N  i=1 3ui⊗ui−I 2(13) where uiistheunitvectorparalleltothesymmetryaxis of the particle i,⊗represents the tensorial product and Iis the identity matrix. Its largest eigenvalue, S,andits associated eigenvector N, characterise the instantaneous global nematic order. These properties help us to characterise the reorientation process when the magnetic field is applied. In addition, the shape and nematic texture within the bridge are characterised by calculating density and orientational profiles. The density profile, ρ(x,y,z),isobtained as: ρ(x,y,z)≡N  i=1 δ(xi−x)δ(yi−y)δ(zi−z), (14) where (xi,yi,zi)are the Cartesian coordinates of the particle i,δdenotes Dirac’s delta and ...stands for the time average, typically over the last 105steps of each simulation. On the other hand, the orientational order profile is given by the local tensor order parameter Q(x,y,z)≡1 ρ(x,y,z)N  i=1 3ui⊗ui−I 2δ(xi−x) δ(yi−y)δ(zi−z), (15) The scalar nematic order parameter profile S(x,y,z)and the local director field n(x,y,z)are obtained as the largesteigenvalueandthecorrespondingeigenvectorof Q(x,y,z), respectively. It is also possible to get the biaxiality profile as the difference between the two smallest eigenvalues, but we see that this quantity does not providefurtherinformation,sowewillskipitinouranalysis. For the numerical evaluation of the density profile and local tensor order parameter, the box is divided in a grid of cubic voxels of unit side length (in reduced units), so thedensityandlocaltensororderprofilesareobtained by averaging in each voxel. This voxel size is chosen as a compromise between coarse-graining and statistical precision, since smaller sizes are subject to larger statistical fluctuations, while larger sizes lead to a blurred picture of the density and order parameter profiles. 3. Results 3.1. a =1 We first considered slit pores in which the adimensional wall interaction strength a=1. Note that corresponds tothelargestvalueconsideredinRef.[36]associatedto hourglass-like bridges. Surprisingly, the first simulations at zero magnetic field for each Dchange the shapes of the bridges, characterised by bigger contact angles but still hourglass-like, and in some cases lead to different orientational profiles. The explanation for these discrepancies mayberelatedtothefactthat,intheMonteCarlosimulations, the GB potential was shifted so takes zero value at the cutoff radius, while in the GB potential implemented inLAMMPSisunshifted.Theshiftdoesnotaffectthe intermolecular forces, but it will alter the torques with MOLECULAR PHYSICS 5 Figure 1. Snapshots of xy,yz and xz cross sections of a bridge in a slit pore of two horizontal walls separated by a distance D=20 and a=1 at the end of the first λ=0 simulation, at λ=0.1 and at the end of the last λ=0 simulation of the magnetic field application cycle.Thecolourcodeassociated totheGB particle orientationsis the following:redifthemainsymmetry axesoftheparticlesarealigned with the zaxis, blue if they are aligned with the yaxis and green if they are aligned with the xaxis. respecttotheunshiftedversionduetoitsorientation dependence. However, as we will see, we still observe a transition from a nearly axisymmetrical configuration where most the particles are oriented along the zaxis for D=20toaconfigurationinwhichmostparticlesareori- ented along a transversal axis in the xy plane for D=30 and D=40. So, we will explore if transitions between different bridge states can be obtained from the magnetic field application. 3.1.1. D =20 We start with the case D=20. We consider in this case an x−direction magnetic field application cycle. Figure 1 shows different cross sections of a bridge at the final configurations of the first λ=0 simulation of the cycle, at the maximum value λ=0.1 and the last λ=0simula- tion of the cycle. In the first zero-field configuration, we observe that the bridge is slightly sheared (not perfectly axisymmetric). Most of the particles orient parallel to the zaxis, except in the interfacial region, in which they preferentially orient homeotropically to the interface. In this region, the nematic order parameter Ssmoothly decays from the bulk value inside the bridge core to zero outside the bridge, as well as the density (not shown by simplicity). This is a common feature of all the systems consideredinthiswork.Thesefindingsqualitativelyagree with the results reported for Monte Carlo simulation in Ref. [36]. When the magnetic field is turned on, particles orient quickly along the xaxis (in fact, at the end of the λ=0.025 simulation most of the particles are already oriented along the xaxis). But, additionaly, the bridge is stretched along the yaxis (the larger the value of λis, the larger the stretching). This effect may be explained by the fact that this arrangement maximises the surface perpendicular to the magnetic field. Such conformation optimises the alignment of particles along the xaxis, reducing the magnetic contribution to the energy. Finally, the configurations corresponding to the final λ=0sim- ulation is quite similar to the one from the initial λ=0 simulation, indicating that, under these conditions, the magnetic field cycle does not induce a transition to a different nanobridge conformation. Figure 2plots the order parameter Sprofiles and nematic director fields corresponding to selected cross sections (equatorial plane z=10, as well as intermediate heights z=5andz=15) of the bridge for the initial and final λ=0 simulations, as well as the simulation for λmax =0.1. For the first λ=0 simulation, the nematic director field confirms our previous observation that most particles orient along the zaxis in the central part of the bridge but mostly homeotropically on the interfacial region. However, there is no indication of a disclination ring located on the equatorial plane as it was reported in Ref. [36], but instead a smooth crossover compatible with an escaped radial texture [49,50]. In the equatorial plane 6P. ROMERO-LLORENTE AND J. M. ROMERO-ENRIQUE Figure 2. Plots of the nematic order parameter Sprofile (colour map) and nematic director field (nail representation, where the nail is parallel to the projection of tilted nematic director in the figure plane, and a perpendicular segment is plotted on the extreme of the nail with the lowest component of the director in the direction normal to the figure plane) of a cross section of the bridge for a=1 and D=20 at heights z=5 (upper row), z=10 (middle row) and z=15 (lower row). Left column corresponds to the initial λ=0 simulation, central column to the λ=0.1 simulation and right column to the final λ=0 simulation. there is indication of two +1/2 vertical disclinations lines located in nearly opposite positions close to the interface, where the nematic order parameter Sis significantly reduced with respect to the bulk value and the nematic director rotates +πradians in a counterclock closed loop around the disclination core. The z=5andz=15 cross sections indicate that both disclinations emerge from the interface. For λ=0.1, the nematic order parameter Sis uniform up to the interfacial region, with no indications of disclination lines. Note that the Svalues are larger than thesecorrespodingtotheλ=0cases.Ontheotherhand, it is confirmed that the bridge shape is elongated along the ydirection while the nematic director is oriented mainly on the xaxis throughout the bridge. Finally, the last λ=0 simulation of the cycle shows a nematic order parameter profile and nematic director field very similar to the corresponding one in the initial λ=0 simulation. This strongly suggests that, after the cycle, the bridge returns to the initial state. Thus, the magnetic field cycle has not induced any transition in the bridge configuration for this case. 3.1.2. D =30 and D =40 We turn now to the D=30 and D=40 cases. Similarly to the D=20 case, the initial λ=0 simulation from an initial condition obtained from Monte Carlo simulations deform the bridges, but they have similar features to those reported in Ref. [36], as it can be seen in Figure 3.Under the application of a magnetic field along the zdirection, particles reorient vertically, but in addition it amplifies any top-bottom asymmetry that the bridge may present, so now the contact area becomes smaller in one wall and MOLECULAR PHYSICS 7 Figure 3. Final snapshots of the xy,yz and xz cross sections of a bridge in a slit pore of two horizontal walls with a=1and separated by a distance D=30 (left column) and D=40 (right column) after a λ=0 simulation with 2 ×105steps. The colour code associated to particle orientations is the same as in Figure 1. larger in the other. If the magnetic field is strong enough, the bridge detaches from the wall where its contact surface is smaller and a single droplet is formed. Figure 4 shows this detachment for the D=30 case and λ= 0.125 when the magnetic field intensity starts to decrease. A similar instability is observed for D=40. This process is irreversible and cannot be reverted by switching off the magnetic field. Figure 4. Time evolution of the bridge collapse for D=30, a=1andλ=0.125. The xz and yz central sections of the bridge snapshots are shown for the time step t=105(left column), t=1.2 ×105(middle column) and t=1.4 ×105(right column). 3.2. a =1.25 Now we increased the wall-particles interaction parameter to a=1.25. In this case, bridges have lower contact angles with the walls. As in the previous value of a,we analyse the effect of the magnetic field application on bridges with D=20, D=30 and D=40. 3.2.1. D =20 Figure 5shows the evolution of the instantaneous nematic order parameter Sin a magnetic field application cycle. In a first stage, we followed a similar procedure to the a=1 case, in which the magnetic field is applied along the xaxis. The initial λ=0 simulation leads to a similar configuration to the obtained for the a=1case (see left columns of Figures 6and 7). As in the previous case, the bridge is slightly sheared and most of the particles are parallel to the zaxis but homeotropicallyoriented in the interfacial region. Again, no indication of a disclination ring is observed, but instead a smooth crossovercompatiblewithanescapedradialtexture. In addition, two vertical disclinations emerge from the interface. As the magnetic field is turned on, the system evolvesinasimilarwayasthea=1case:particlesorient preferentially along the xaxis and the bridge stretches along the ydirection (not shown by simplicity). It is worth noticing that Schanges abruptly at the beginning of the simulations for λ=0.025 and λ=0.05, indicating major particle reorientations in the bridge as the region of particles oriented along the xaxis grows from the interface neighbourhoodtothecentreofthebridge.However,as the value of λisreducedto0,weobservethatthefinal stage differs from the initial one, as it can be seen at the middle columns of Figures 6and 7.After4×105 8P. ROMERO-LLORENTE AND J. M. ROMERO-ENRIQUE Figure 5. Time evolution of the instantaneous global nematic order parameter Salong the magnetic field application cycle for D=20 and a=1.25. For simplicity, the values are represented each 102steps. The values of λin blue correspond to the application of the magnetic field along the xaxis. On the other hand, red values of λrepresent magnetic fields along the zaxis. Figure 6. Snapshots of the xy,yz and xz cross sections of a bridge in a slit pore of two horizontal walls separated by a distance D=20 and a=1.25. The left column corresponds to the final configuration of the initial λ=0 simulation, middle column to the final configuration of the intermediate λ=0 simulation after the application of magnetic fields along the xaxis, and the right column to the final configuration of the last λ=0 simulation after the application of magnetic fields along the zaxis. The colour code associated to particle orientations is the same as in Figure 1. MOLECULAR PHYSICS 15 [43] W.M. Brown, M.K. Petersen, S.J. Plimpton and G.S. Grest, J. Chem. Phys. 130 (4), 044901 (2009). doi:10.1063/1. 3058435 [44] A.P. Thompson, H.M. Aktulga, R. Berger, D.S. Bolint i n e anu , W. M. B row n , P. S . C roz i er, P. J. i n ’ t Vel d, A . Kohlmeyer, S.G. Moore, T.D. Nguyen, R. Shan, M.J. Stevens, J. Tranchida, C. Trott and S.J. Plimpton, Comp. Phys. Comm. 271, 108171 (2022). doi:10.1016/j.cpc.2021. 108171 [45] R. Berardi, C. Fava and C. Zannoni, Chem. Phys. Lett. 297 (1–2), 8–14 (1998). doi:10.1016/S0009-2614(98)010 90-2 [46]J.W.Perram,J.Rasmussen,E.PræstgaardandJ.L. Lebowitz, Phys. Rev. E 54 (6), 6565–6572 (1996). doi:10.1103/PhysRevE.54.6565 [47] S.Nosé,J.Chem.Phys.81 (1), 511–519 (1984). doi:10. 1063/1.447334 [48] W.G. Hoover, Phys. Rev. A 31 (3), 1695–1697 (1985). doi:10.1103/PhysRevA.31.1695 [49] P.E. Cladis and M. Kléman, J. Phys. France 33 (5–6), 591–598 (1972). doi:10.1051/jphys:01972003305-605 9100 [50] R.B. Meyer, Phil. Mag. 27 (2), 405–424 (1973). doi:10. 1080/14786437308227417