1. Introduction
Anisotropic interacting particle systems exist widely in nature, i.e. spin systems [1], molecules [2], patchy colloidal particles [3, 4], liquid crystals [5], etc. Some of these systems, such as liquid crystals, already have essential applications to everyday life. The study of the dynamics of these systems has both high theoretical and practical significance.
The diffusion motion of anisotropic particles is important for many processes, such as the diffusion and transport of biological macromolecules in cells [6], colloidal self-assembly [7–9], diffusion-limited aggregation [10], and the design of artificial molecular machines [11–13]. Diffusion also plays an important role in some phase transitions. For example, in some solid–solid phase transitions, the transition between two solid structures is caused by the diffusion of particles in a mediated liquid state [14]; the anisotropic diffusion of colloids can induce a solid–liquid transition [15]. The diffusion of anisotropic particles differs from that of isotropic particles and depends on the particle shape or interaction potential. Previous studies have shown that the Brownian motion of a single nonspherical particle in solution [16–20] is different from that of an isotropic particle and is significantly impacted by rotational-translational coupling [21–24]. Liquids composed of particles with different degrees of anisotropy also have different collective diffusion or self-diffusion behaviours [25].
The Janus particle, a patchy particle with a single patch, is a typical anisotropic interacting particle that has been extensively studied over the past two decades [3, 26–28]. This particle is a useful model for understanding the mechanism of phase transitions and the dynamics of anisotropic particle systems [29, 30]. The anisotropic properties of some Janus particles derive from the surface modification of spherical colloidal particles [4, 31–33]. The interactions between such Janus particles are clearly anisotropic, but the particle properties are significantly different from those of rod-like [34–36] and ellipsoid-like particles [18, 37].
In this study, we focus on two-dimensional (2D) anisotropic liquids composed of Janus particles. These particles are usually spherical in shape, so the hydrodynamic effects exerted on them by the background liquid can be considered isotropic. On the other hand, the interaction between the particles is anisotropic. Hereinafter we refer to such particles as nearly spherical particles. Molecular dynamics simulations are performed to study the collective diffusion of a 2D liquid of Janus particles. The static spatial structure and dynamic time correlation of this anisotropic liquid are analysed to identify similarities to and differences from isotropic Lennard-Jones (LJ) liquids [38–41]. We show that anisotropic liquids composed of nearly spherical particles exhibit anisotropy only on short timescales. Long-term averaging of the system results in similar diffusion behaviour to that of an isotropic liquid, similar to the high temperature behavior of the 2D crystal composed of Janus particles [29, 30]. Finally, we discuss the origin of the differences in the diffusion properties of liquids consisting of nonspherical (rod-like) [34, 36] and nearly spherical particles.
2. Model and methods
We study patchy particles with attractive interactions similar to antiferromagnetic interactions [42–45]. Those so-called inverse patchy particles behave like a ferromagnetic system [46–48]. The modulated LJ pair potential proposed in [49, 50] is used to describe the attractive interaction between patchy particles. In this model, an anisotropic factor dependent on the orientations of the particles is introduced to vary the depth of the LJ potential well
$\begin{eqnarray}{U}_{{ij}}({{\boldsymbol{r}}}_{{ij}},{\alpha }_{i},{\alpha }_{j})=f({\hat{{\boldsymbol{r}}}}_{{ij}},{\alpha }_{i},{\alpha }_{j}){U}_{{\rm{LJ}}}({r}_{{ij}}).\end{eqnarray}$
The LJ potential is $\begin{eqnarray}{U}_{\mathrm{LJ}}({r}_{{ij}})=\left\{\begin{array}{ll}4\varepsilon \left[{\left(\displaystyle \frac{\sigma }{{r}_{{ij}}}\right)}^{12}-{\left(\displaystyle \frac{\sigma }{{r}_{{ij}}}\right)}^{6}\right], & {r}_{{ij}}\lt {r}_{c};\\ 0, & {r}_{{ij}}\geqslant {r}_{c},\end{array}\right.\end{eqnarray}$
where rij is the distance between the centres of mass of particles i and j, σ characterizes the particle diameter, ε is the depth of the potential well, and rc is the cutoff distance. The anisotropic factor is $\begin{eqnarray}f({\hat{{\boldsymbol{r}}}}_{{ij}},{\alpha }_{i},{\alpha }_{j})=\left\{\begin{array}{ll}1, & {r}_{{ij}}\leqslant \sigma ;\\ \exp \left(-\displaystyle \frac{{\alpha }_{i}^{2}+{\alpha }_{j}^{2}}{2{\sigma }_{p}^{2}}\right), & {r}_{{ij}}\gt \sigma .\end{array}\right.\end{eqnarray}$
αi is the angle between the patch vector pi and the interparticle vector rij (figure 1(a)). σp represents the angular range of the attractive patch and regulates the sensitivity of f to αi,j. When σp ≫ π, f ∼ 1, and the interaction between particles tends to an isotropic LJ potential. For simplicity, when σp < π, we introduce the patch coverage $\chi \,=[1-\cos ({\sigma }_{p})]/2$ to replace σp; in the simulation we set σp = 30 ≫ π to make the system exhibit isotropic LJ interactions approximately. The contour plot of the pair potential is schematised in figure 1(b), where $U(r={r}_{{ij}},\alpha \,=\sqrt{{\alpha }_{i}^{2}+{\alpha }_{j}^{2}})=\exp (-{\alpha }^{2}){U}_{\mathrm{LJ}}({r}_{{ij}})$. Two particles with facing patches have the strongest attraction at αi = αj = 0, and the potential energy reaches the minimum at a distance of $\sqrt[6]{2}\sigma $. When αi or αj is large relative to σp, f(rij > σ) → 0, only the repulsive part in the potential takes effect, and the two spheres behave as hard spheres with a diameter σ. The interaction force is set to zero beyond the cutoff distance rc = 3.5. rc = 5 is tested in some simulations and no difference is found in the results. From the perspective of interaction potential, patchy particles are still nearly spherical (figure 1(b)). The rotations of Janus particles are not hindered by space occupation as liquid crystal particles do. We use the reduced dimensionless units: the mass m = 1, energy ε = 1, and diameter σ = 1. Thus, the time unit is ${t}^{* }=\sqrt{\epsilon /(m{\sigma }^{2})}=1$.
Figure 1. (a) Parameters of two interacting patchy particles in equations ( |
A Langevin thermostat is applied to simulate the NVT ensemble of the liquid state of the patchy sphere system. The Langevin thermostat [51] simultaneously defines the dynamics of the system in this study:
$\begin{eqnarray}{m}_{i}^{\mu }{\ddot{X}}_{i}^{\mu }=-{\partial }_{{X}_{i}^{\mu }}U(\{{X}_{i}\})-{\gamma }_{{ij}}^{\mu }{\dot{X}}_{j}+{\eta }_{i}^{\mu }(t),\end{eqnarray}$
where μ = T, R and denotes the translational/rotational motion of a particle, respectively. For μ = T(, R), mμ = m(, I) is the mass (rotational inertia) of the particle, Xμ is the position (orientation), γμ is the translational (rotational) friction coefficient, and ημ is the random force (torque) applied to the particle. γ is the hydrodynamic damping coefficient matrix of the environment; here, no hydrodynamic coupling between particles is considered, so the diagonal form is used, i.e. γij = δijγii. The random force and torque ${\eta }_{i}^{\mu }$ $\begin{eqnarray}\left\langle {\eta }_{i}^{\mu }(t)\right\rangle =0,\end{eqnarray}$
$\begin{eqnarray}\left\langle {\eta }_{i}^{\mu }(t){\eta }_{j}^{\nu }(t^{\prime} )\right\rangle =2{\gamma }^{\mu }{k}_{{\rm{B}}}T{\delta }_{i,j}{\delta }_{\mu ,\nu }\delta (t-t^{\prime} ).\end{eqnarray}$
The temperature is measured in ε, i.e. kB = 1.The system consists of 40 000 patchy spheres in a 2D box with periodic boundary conditions in the NVT ensemble. All simulations are carried out using LAMMPS [52], where we insert a piece of code to implement the potential defined in equations (1 )–(3 ). At every time integral step, the z-components of all the forces are cancelled to confine the system in 2D, whereas the torques in all directions are taken into account to model the 3D rotations of spheres. The time integration is performed using the LAMMPS' default velocity-Verlet algorithm with a time-step dt = 0.002. A simulation is performed for each set of parameters (σp, ρ, T), where the spheres are placed on a simple square lattice in the initial state and then relaxed for a sufficiently long time (t ≥ 2 × 104) to reach an equilibrium state. All statistics and analyses are performed in the equilibrium state. We simulate patch coverages of χ = 0.3, 0.5 and LJ cases. T is varied from 0.15 to 0.9, and ρ is varied from 0.4 to 0.9 for all temperatures. For some parameters, the system is in a solid or gas state. However, all the systems discussed in this paper are in a liquid state.
3. Static structure of the anisotropic liquid
Anisotropic liquids exhibit many different properties from isotropic liquids. These differences depend on the factors that cause anisotropy, such as geometric asphericity. Liquid crystals are typical complex liquids with a large anisotropy due to asphericity. Compared to liquid crystal particles, Janus particles are more spherical (as can be seen from the contour plot of the potential presented in figure 1(b)), necessitating a more refined analysis of the anisotropic characteristics of liquid Janus particles.
We investigate the static structural properties of patchy-sphere liquids by introducing some spatial correlation functions. The first is the radial distribution function
$\begin{eqnarray}{g}_{r}(r)=\displaystyle \frac{1}{2\pi {rN}\rho }\left\langle \sum _{i}^{N}\sum _{j\ne i}^{N}\delta (r-{r}_{{ij}})\right\rangle ,\end{eqnarray}$
where N is the number of particles, ρ is the number density of the system, and ⟨ · ⟩ is the ensemble averaging. gr is an important measure of the static structure of liquids. Because we study the properties of systems with different number densities, for convenience we use the local number density ρ(r) = ρgr(r) rather than the radial distribution function gr(r) to describe the spatial correlation of particles. Figure 2 shows how ρr(r) depends on the patch coverage χ and temperature T. In figure 2(a), the first few peaks of g(r) of the anisotropic liquid are lower in height than those of the LJ liquid. This result can be attributed to the anisotropic interactions, which produce more dispersion in the distance between neighbouring particles than that for the LJ liquid. The peak position of g(r) of a patchy particle liquid is slightly smaller than that of the LJ case. This difference is due to the statistical contribution of neighbouring particles with opposite orientations, which behave as hard spheres with a diameter σ (comparing to the equilibra distance $\sqrt[6]{2}\sigma $ of the isotropic LJ particles). The heights of the peaks in g(r) in figure 2(b) decrease with increasing temperature. This result is consistent with that for LJ liquids, in which a higher temperature results in a larger fluctuation of the particle distance. The system shown in figure 2(a) (ρ = 0.8) has more peaks than that shown in figure 2(b) (ρ = 0.6). This result means that the system has larger positional correlation lengths at higher densities.
Figure 2. The local number density ρ(r) for the system against (a) patch coverage degree χ and (b) temperature T. The system density is ρ = 0.8 and 0.6 for (a) and (b), respectively. |
We then define the spatial correlation function of the particle orientation as
$\begin{eqnarray}{g}_{\phi }(r)=\langle \cos (2[{\phi }_{i}(0)-{\phi }_{j}(r)])\rangle ,\end{eqnarray}$
which is plotted in figure 3. This correlation depends on the orientations of the particles, so it is only valid for anisotropic liquids, not for isotropic liquids. Particle orientations are negatively correlated in the first nearest neighbour range (r < 2) due to attractive interactions between patches. The particle orientation is positively correlated in the second nearest neighbour range (2 < r < 3) and so forth. This result is similar to the low-temperature behaviour of anisotropic crystals. For the system with χ = 0.3 and T = 0.15 (figure 3(a)), an orientation correlation is still observed for the fourth nearest neighbour. As the patch coverage increases (e.g. χ = 0.5 in figure 3(a)), the interaction tends to be more isotropic, making the orientation correlation decrease. The correlation length of the orientation decreases accordingly, thus the orientation correlation for the system with χ = 0.5 can only be observed within the range of the first two nearest neighbours. Increasing the temperature also leads to the shortening of the orientation correlation length, and the orientation of the system becomes more disordered (e.g. T = 0.7 in figure 3(b)).
Figure 3. The spatial correlation functions of the orientation gφ for the system with ρ = 0.8 against the (a) patch coverage degree χ and (b) temperature T, respectively. |
4. Diffusion behaviour of the anisotropic liquid
The collective diffusion of liquid particles is usually characterized by the mean squared displacement (MSD):1 ). For the system with χ = 0.3, the translational (rotational) harmonic period is
$\begin{eqnarray}{\rm{MSD}}(t)={ \langle {\left({\boldsymbol{r}}({t}_{0}+t)-{\boldsymbol{r}}({t}_{0})\right)}^{2} \rangle }_{{t}_{0}}.\end{eqnarray}$
The long-time behaviour of the MSD (t → ∞) can be used to characterise the type of diffusion and thereby infer the diffusion mechanism. For normal diffusion, the MSD is proportional to the time interval t, and obeys the Einstein relation. In the 2D case, $\begin{eqnarray}\mathrm{MSD}(t)=4{Dt}.\end{eqnarray}$
If the MSD(t) is significantly nonlinear, the diffusion is anomalous, where $\begin{eqnarray}\mathrm{MSD}(t)\sim {t}^{\alpha }.\end{eqnarray}$
Figure 4 shows the long- and short-time MSD against the patch coverage χ and temperature T. The system density is 0.8. To obtain an intuitive understanding of the time interval shown in the figure, it is useful to consider the period of the harmonic motion of two particles with a pair interaction at the minimum potential energy, i.e. ${r}_{{ij}}=\sqrt[6]{2}\sigma $ and αi = αj = 0 in equation ( ${t}^{T}=2\pi \sqrt{m/| {\partial }^{2}U/\partial {r}_{{ij}}^{2}| }=0.83$
${(t}^{R}=2\pi \sqrt{I/| {\partial }^{2}U/\partial {\alpha }_{i}^{2}| }=1.63$).
For the system with χ = 0.5, the rotational period increases to tR = 2.21, while the translational period remains unchanged. The long-time MSD curve of the anisotropic liquid satisfies normal diffusion (MSD ∼ t0.99), which is consistent with the behaviour of the LJ liquid (figures 4(a) and (c)). The diffusion coefficient decreases with increasing patch coverage, which is evident from a comparison of the two curves for χ = 0.3, 0.5 and the LJ curve in figure 4(a). This result is obtained because as the patch coverage decreases, the attractive interaction among particles becomes weaker, the liquid viscosity decreases, and particle diffusion is faster. For both anisotropic and LJ liquids, the short-time (t < 0.1) MSD follows a power-law relationship t1.9, which is close to ballistic diffusion. The diffusion behaviour in this range originates from the rarity of interparticle collision events.
Figure 4. Left column: the MSD for the system with ρ = 0.8 against the (a) patch coverage degree χ and (c) temperature T. Right column: (b) and (d) show the MSD⊥,∥ in the short-time range corresponding to (a) and (c), respectively. The black dashed line in the figure is tα obtained by fitting equation ( |
For anisotropic liquids, the displacement vector of the particles can be correlated with the particle orientation. Therefore, the particle displacement over a time interval can be decomposed into two vectors that are parallel and perpendicular to the orientation p of the particles (figure 4(e)). We introduce the parallel (perpendicular)-MSD [19] as
$\begin{eqnarray}{\mathrm{MSD}}_{\parallel }(t)={\left\langle {\left({\rm{\Delta }}{\boldsymbol{r}}(t)\cdot {\boldsymbol{p}}({t}_{0})\right)}^{2}\right\rangle }_{{t}_{0}};\end{eqnarray}$
$\begin{eqnarray}{\mathrm{MSD}}_{\perp }(t)={\left\langle {\left({\rm{\Delta }}{\boldsymbol{r}}(t)\cdot {\boldsymbol{n}}({t}_{0})\right)}^{2}\right\rangle }_{{t}_{0}},\end{eqnarray}$
in which Δr(t) = r(t0 + t) − r(t0) and n is a unit vector perpendicular to p in the plane of the particles' translational motion.The difference between MSD∥ and MSD⊥ in anisotropic liquids is only observed in systems with low temperatures and small patch coverages (e.g. see the curves for T = 0.15, χ = 0.3 in figures 4(b) and (d)) for mid-time intervals (0.5 < t < 5). This result shows the particularity of anisotropic liquids composed of Janus particles. Microscopically, Janus particles always tend to form small clusters in which the attractive patches are near to each other (figure 4(e)). When the time interval t exceeds the short-time range, these particle clusters are more frequently destroyed and reorganised due to the effect of thermal noise. The particles undergo different perpendicular and parallel motions relative to their orientation vectors during this process. The neighbours of an anisotropic particle i in a small cluster have a high probability to locate ahead in the direction of the particle i's orientation vector pi, and the potential well for the parallel motion of particle i is deeper, i.e. it is easier for particle i to leave the cluster by perpendicular motion than by parallel motion. Therefore, the diffusion coefficient associated with MSD∥ is smaller than that associated with MSD⊥ (see the blue curves in figures 4(b) and (d)). When χ increases, the differences between MSD∥ and MSD⊥ disappear (see the red curves in figure 4(b)). We recognize this feature to be unique to the diffusion of nearly spherical particles.
For anisotropic liquids, diffusion is also associated with the rotational degrees of freedom of particles. The MSDφ for the patch vector's azimuth φ is employed to investigate the rotational diffusion:
$\begin{eqnarray}{\mathrm{MSD}}_{\varphi }(t)={\left\langle {\left(\varphi ({t}_{0}+t)-\varphi ({t}_{0})\right)}^{2}\right\rangle }_{{t}_{0}}.\end{eqnarray}$
MSDφ is shown in figure 5. The speed of rotational diffusion can be increased either by increasing the patch coverage χ (figure 5(a)) or the temperature T (figure 5(b)). More importantly, the rotational diffusion is normal diffusion in the short time range; at t ∼ 10, MSDφ reaches π/2 and does not increase further at larger t. ${\mathrm{MSD}}_{\varphi }\to \tfrac{\pi }{2}$ does not result from diffusion stagnation but from the uniform distribution of φ ∈ [0, π]. As the patch coverage or the temperature increases, the characteristic time that the rotational diffusion used to produce a uniform distribution of φ decreases. After this characteristic time, rotational diffusion leads to completely disordered orientations for the particles. Therefore, the particle orientations disappear statistically under time-averaging. This result also explains why the long-time behaviour of the MSD of the anisotropic liquid is consistent with the isotropic result (figures 4(a) and (c)) and the difference between MSD∥ and MSD⊥ disappears over long time intervals (figures 4(b) and (d)).
Figure 5. The rotational MSDφ for the system with ρ = 0.8 against the (a) patch coverage degree χ and (b) temperature T, respectively. |
To further investigate the inconsistency between MSD∥ and MSD⊥, figure 6 shows the joint probability distributions of the particle orientations φ(t0) and φ(t0 + t) before and after the time interval t, where φ is the clockwise angle between the patch vector and the displacement vector. In this system, χ = 0.3. At a higher temperature T = 0.5 (see the upper row of figure 6), φ(t0) and φ(t0 + t) are distributed uniformly in the range (0, 2π) for all t. However, at a lower temperature T = 0.15 (see the lower row of figure 6), the distribution is not uniform at mid-time intervals (i.e. t = 1 in figure 6(e)). The probability distribution exhibits two peaks at approximately π/2 and 3π/2, showing that the perpendicular motion of the particles has a higher probability than parallel motion. As t becomes larger (figure 6(f)), the nonuniformity of this distribution weakens or even disappears, and the correlation between φ(t0) and φ(t0 + t) becomes weaker. Going from figures 6(a) to (c) (or from figures 6(d) to (f)), the bright stripe widens with increasing t. φ(t0) and φ(t0 + t) eventually become statistically irrelevant. This result is consistent with the long-time behaviour of MSDφ shown in figure 5.
Figure 6. Upper row: (a), (b) and (c) joint probability distributions of the particle orientations φ(t0) and φ(t0 + t) for the system at T = 0.5 before and after the time intervals t = 0.1, 1, and 10, respectively. Lower row: (d)–(f) show the corresponding distributions for the system at T = 0.15. |
5. Discussion
In this study, we investigate the structure and collective diffusion of a 2D anisotropic liquid composed of Janus particles, and compare the results against those of an isotropic LJ liquid. Janus particles characterized by the model given in equation (1 ) are more spherical than typical liquid crystal [5], ellipsoidal [18, 37], or rod-like [35] particles. This result can be observed simply by examining the equipotential surfaces of the corresponding potential functions (figure 1(b)).
For nonspherical particles, as long as the particle distance is smaller than the long ellipsoidal axis, the rotation of one particle hinders another particle from passing through all orientations. The translational properties of the particles along the directions of the long and short axes are quite different, and the positional order of the system is destroyed before the orientational order vanishes. The special structure of the liquid crystal system, e.g. the nematic phase, derives from this mechanism.
The anisotropy of the Janus particles studied here arises from the inhomogeneous nature of the spherical particle surface [32, 33], not the geometric shape. Therefore, for various models at low temperatures, the orientational order of the system can be destroyed before the positional order vanishes [29]. Such systems are more similar to spin systems than liquid crystals. As Janus particles can rotate freely in the liquid state, the long-time average behaviour of the liquid phase is close to that of an isotropic liquid. Anisotropic properties manifest only during short time intervals. This conclusion can be extended to anisotropic liquids composed of nearly spherical particles.
To summarize, differences between isotropic liquids and anisotropic liquids composed of nearly spherical and nonspherical particles have been identified in this study. The results of this study can be used to deepen our understanding of the diffusion of anisotropic liquid and subjects such as molecular machine design and colloidal self-assembly.