Hydrodynamic correlations in isotropic fluids and liquid crystals simulated by multi-particle collision dynamics

TL;DR

This paper analyzes multi-particle collision dynamics for simulating isotropic fluids and liquid crystals, focusing on hydrodynamic fluctuations, relaxation, and the extension of the method to nematic liquid crystals.

Contribution

It extends the multi-particle collision dynamics method to include nematic liquid crystals and analyzes hydrodynamic and orientational fluctuations.

Findings

Hydrodynamic fluctuations relax towards equilibrium as described by a linearized theory.

Flow fluctuations can induce orientation fluctuations in liquid crystals.

Observable effects on velocity correlations occur only at certain simulation parameters.

Abstract

Multi-particle collision dynamics is an appealing numerical technique aiming at simulating fluids at the mesoscopic scale. It considers molecular details in a coarse-grained fashion and reproduces hydrodynamic phenomena. Here, the implementation of multi-particle collision dynamics for isotropic fluids is analysed under the so-called Andersen-thermostatted scheme, a particular algorithm for systems in the canonical ensemble. This method gives rise to hydrodynamic fluctuations that spontaneously relax towards equilibrium. This relaxation process can be described by a linearized theory and used to calculate transport coefficients of the system. The extension of the algorithm for nematic liquid crystals is also considered. It is shown that thermal fluctuations in the average molecular orientation can be described by an extended linearized scheme. Flow fluctuations induce orientation fluctuations. However, orientational changes produce observable effects on velocity correlation functions only when simulation parameters exceed their values from those used in previous applications of the method. Otherwise, the flow can be considered to be independent of the orientation field.

Equations

$$r_{i,\alpha }(t+\Delta t)=r_{i,\alpha }(t)+v_{i,\alpha }(t)\Delta t.$$

$$v_{i,\alpha }(t+\Delta t)={\bar{v}}_{\alpha }^{\text{c}}(t)+{\xi }_{i,\alpha }-{\bar{\xi }}_{\alpha }^{\text{c}},$$

$$\Delta L_{\alpha }^{\text{c}}=\sum _{j\in c}{\epsilon }_{\alpha {\alpha }^{\prime }{\beta }^{\prime }}\left(r_{j,{\alpha }^{\prime }}-{\bar{r}}_{{\alpha }^{\prime }}^{\text{c}}\right)\left(v_{j,{\beta }^{\prime }}-{\xi }_{j,{\beta }^{\prime }}\right),$$

$$v_{i,\alpha }(t+\Delta t)={\bar{v}}_{\alpha }^{\text{c}}(t)+{\xi }_{i,\alpha }-{\bar{\xi }}_{\alpha }^{\text{c}}-{\epsilon }_{\alpha {\alpha }^{\prime }{\beta }^{\prime }}{J}_{{\alpha }^{\prime }{\alpha }^{\prime \prime }}^{\text{c}}(t)\Delta L_{{\alpha }^{\prime \prime }}^{\text{c}}(t)\left(r_{i,{\beta }^{\prime }}-{\bar{r}}_{{\beta }^{\prime }}^{\text{c}}\right),$$

$$\frac{\partial \rho }{\partial t}=-\frac{\partial }{\partial x_{\alpha }}\left(\rho V_{\alpha }\right),$$

$$\rho \left(\frac{\partial }{\partial t}+V_{\beta }\frac{\partial }{\partial x_{\beta }}\right)V_{\alpha }=-\frac{\partial }{\partial x_{\beta }}{\sigma }_{\alpha \beta },$$

$${\sigma }_{\alpha \beta }=-p{\delta }_{\alpha \beta }+{\eta }_{\alpha \beta {\alpha }^{\prime }{\beta }^{\prime }}\frac{\partial V_{{\beta }^{\prime }}}{\partial x_{{\alpha }^{\prime }}},$$

$${\sigma }_{\alpha \beta }=-p{\delta }_{\alpha \beta }+{\bar{\eta }}^{\text{kin}}\left(\frac{\partial V_{\beta }}{\partial x_{\alpha }}+\frac{\partial V_{\alpha }}{\partial x_{\beta }}-\frac{2}{3}{\delta }_{\alpha \beta }\frac{\partial V_{{\alpha }^{\prime }}}{\partial x_{{\alpha }^{\prime }}}\right)+{\bar{\eta }}^{\text{col}}\frac{\partial V_{\alpha }}{\partial x_{\beta }},$$

$${\bar{\eta }}^{\text{kin}}=\frac{{\bar{N}}^{\text{c}}{k}_{\text{B}}T\Delta t}{a^3}\left(\frac{{\bar{N}}^{\text{c}}}{{\bar{N}}^{\text{c}}-1+{\mathrm{e}}^{-{\bar{N}}^{\text{c}}}}-\frac{1}{2}\right),$$

$${\bar{\eta }}^{\text{col}}=\frac{m}{12a\Delta t}\left({\bar{N}}^{\text{c}}-1+{\mathrm{e}}^{-{\bar{N}}^{\text{c}}}\right),$$

$${\sigma }_{\alpha \beta }=-p{\delta }_{\alpha \beta }+\eta \left(\frac{\partial V_{\beta }}{\partial x_{\alpha }}+\frac{\partial V_{\alpha }}{\partial x_{\beta }}-\frac{2}{3}{\delta }_{\alpha \beta }\frac{\partial V_{{\alpha }^{\prime }}}{\partial x_{{\alpha }^{\prime }}}\right)+{\eta }^{\text{V}}{\delta }_{\alpha \beta }\frac{\partial V_{{\alpha }^{\prime }}}{\partial x_{{\alpha }^{\prime }}},$$

$${\eta }^{\text{kin}}=\frac{k_{\text{B}}T{\bar{N}}^{\text{c}}\Delta t}{a^3}\left(\frac{{\bar{N}}^{\text{c}}}{{\bar{N}}^{\text{c}}-\frac{5}{4}}-\frac{1}{2}\right),$$

$${\eta }^{\text{col}}=\frac{m{\bar{N}}^{\text{c}}}{24a\Delta t}\left(1-\frac{7}{5{\bar{N}}^{\text{c}}}\right),$$

$$\rho (\vec{r},t)={\rho }^{\text{eq}}+\delta \rho (\vec{r},t),$$

$$V_{\alpha }(\vec{r},t)=V_{\alpha }^{\text{eq}}+\delta V_{\alpha }(\vec{r},t),$$

$$⟨{\Sigma }_{\alpha \beta }\left({\vec{r}}^{\prime },t^{\prime }\right){\Sigma }_{{\alpha }^{\prime }{\beta }^{\prime }}\left(\vec{r},t\right)⟩=2k_{\text{B}}T{\eta }_{\alpha \beta {\alpha }^{\prime }{\beta }^{\prime }}\delta \left({\vec{r}}^{\prime }-\vec{r}\right)\delta \left(t^{\prime }-t\right),$$

$$\big{\langle}\tilde{\Sigma}_{\alpha\beta}\big{(}\vec{k}^{\,\prime},\omega^{\prime}\big{)}\tilde{\Sigma}_{\alpha^{\prime}\beta^{\prime}}\big{(}\vec{k},\omega\big{)}\big{\rangle}=2(2\piup)^{4}k_{\text{B}}T\eta_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\delta\big{(}\vec{k}^{\,\prime}+\vec{k}\big{)}\delta\left(\omega^{\prime}+\omega\right).$$

$$\left(-\mathrm{i}\omega +\nu k^2\right)\delta {\tilde{V}}_{\alpha }=\frac{\mathrm{i}k}{{\rho }_0}{\tilde{\Sigma }}_{1\alpha }$$

$$\left(\begin{array}{cc}-\mathrm{i}\omega & -\mathrm{i}k{\rho }_0\\ -\mathrm{i}k{c}_{\text{T}}^2{\rho }_0^{-1}& -\mathrm{i}\omega +D_{\text{l}}{k}^2\end{array}\right)\cdot \left(\begin{array}{c}\delta \tilde{\rho }\\ \delta {\tilde{V}}_1\end{array}\right)=\frac{\mathrm{i}k}{{\rho }_0}\left(\begin{array}{c}0\\ {\tilde{\Sigma }}_{11}\end{array}\right)$$

$$C_{\rho}\big{(}\vec{k},\omega\big{)}=\big{\langle}\delta\tilde{\rho}\big{(}\vec{k},\omega\big{)}\delta\tilde{\rho}^{*}\big{(}\vec{k},\omega\big{)}\big{\rangle},$$

$$C_{\textrm{V}\alpha}\big{(}\vec{k},\omega\big{)}=\big{\langle}\delta\tilde{V}_{\alpha}\big{(}\vec{k},\omega\big{)}\delta\tilde{V}_{\alpha}^{*}\big{(}\vec{k},\omega\big{)}\big{\rangle},$$

$$C_{\rho}\big{(}\vec{k},\omega\big{)}=A\frac{\rho_{0}D_{\textrm{l}}k^{4}}{q^{2}(k)+D_{\textrm{l}}^{2}k^{4}}\left\{\frac{1-\omega/[2q(k)]}{[\omega-q(k)]^{2}+D_{\textrm{l}}^{2}k^{4}}+\frac{1+\omega/[2q(k)]}{[\omega+q(k)]^{2}+D_{\textrm{l}}^{2}k^{4}}\right\},$$

$$C_{\textrm{V1}}\big{(}\vec{k},\omega\big{)}=\frac{\omega^{2}}{\rho_{0}^{2}k^{2}}C_{\rho}\big{(}\vec{k},\omega\big{)},$$

$$C_{\textrm{V2}}\big{(}\vec{k},\omega\big{)}=C_{\textrm{V3}}\big{(}\vec{k},\omega\big{)}=A\frac{1}{\rho_{0}}\frac{\nu k^{2}}{\omega^{2}+\nu^{2}k^{4}}.$$

$$\delta\bar{\rho}\big{(}\vec{k},t_{l}\big{)}=\rho^{\textrm{eq}}\sum_{j=1}^{N}\exp\big{[}-\mathrm{i}\vec{k}\cdot\vec{r}_{j}(t_{l})\big{]},$$

$$\delta\bar{V}_{\alpha}\big{(}\vec{k},t_{l}\big{)}=\sum_{j=1}^{N}v_{j,\alpha}\exp\big{[}-\mathrm{i}\vec{k}\cdot\vec{r}_{j}(t_{l})\big{]},$$

$$C_{\rho}\big{(}\vec{k},t_{l+m}-t_{l}\big{)}=\frac{1}{N_{\text{s}}-m+1}\sum_{l=0}^{N_{\text{s}}-m}\delta\bar{\rho}\big{(}\vec{k},t_{l+m}\big{)}\delta\bar{\rho}\big{(}\vec{k},t_{l}\big{)},$$

$$Q_{\alpha \beta }^{\text{c}}=\frac{1}{2N^{\text{c}}}\sum _{j\in c}\left(3u_{j,\alpha }{u}_{j,\beta }-{\delta }_{\alpha \beta }\right),$$

$$U_{\text{mf}}=-\frac{3}{2}U{S}^{\text{c}}\left(u_{i,\alpha }{u}_{i,\beta }{n}_{\alpha }^{\text{c}}{n}_{\beta }^{\text{c}}-1\right),$$

$$u_{i,\alpha }\left(t+\Delta t\right)=u_{i,\alpha }^{\prime }+{\chi }_{\text{HI}}\left[{\Omega }_{\alpha \beta }^{\text{c}}{u}_{i,\beta }^{\prime }+\lambda \left(D_{\alpha \beta }^{\text{c}}{u}_{i,\beta }^{\prime }-D_{{\alpha }^{\prime }{\beta }^{\prime }}^{\text{c}}{u}_{i,{\beta }^{\prime }}^{\prime }{u}_{i,{\alpha }^{\prime }}^{\prime }{u}_{i,\alpha }^{\prime }\right)\right],$$

Full text

Key words: multi-scale simulation techniques, particle-based simulations of fluids, thermal fluctuations, nematic liquid crystals, structure of liquids and liquid crystals

PACS: 61.20.Ja, 05.20.Jj, 05.40.-a, 83.80.Xz

Ключовi слова: ìóëüòèìàñøòàáíèé ìåòîä ìîäåëþâàííÿ, чàñòèíêîâî áàçîâàíå ìîäåëþâàííÿ ïëèíiâ, òåïëîâi ôëóêòóàöi¿, íåìàòèчíi ðiäêi êðèñòàëè, ñòðóêòóðà ðiäèí i ðiäêèõ êðèñòàëiâ

1 Introduction

Simulation of complex fluids, e.g., chemically reacting systems, colloids, and liquid crystals (LCs), represents a major task in computational physics of condensed matter [1, 2]. Such fluids are characterised by phenomena occurring at widely separated length and time-scales, all of which are relevant in determining the observed behaviour. The interest in complex systems usually focuses on the dynamics of some microscopic degrees of freedom interacting with a solvent. Although it is essential to reproduce the correct behaviour of the solvent over long distances and large times, its molecular details are irrelevant. Traditional simulation techniques, e.g., molecular dynamics (MD), are very successful in describing equilibrium states of atomistic ensembles [3], though it is not feasible to use them to simulate the solvent due to the enormous number of degrees of freedom and the extremely large time ranges that must be covered [4].

To incorporate the collective effects, e.g., thermal fluctuations and flow, in simulations of liquids is a challenge that has motivated the development of novel computational approaches which take into account essential dynamical properties and allow for coupling with the interesting microscopic degrees of freedom, and yet are simple enough to be simulated for long times and distances. Stokesian dynamics [5] and the lattice Boltzmann method [6] are examples of such approaches. Multi-particle collision dynamics (MPCD) is probably the most successful of such mesoscopic algorithms. It was introduced by Malevantes and Kapral [1, 2], and simulates fluids by means of particles whose positions and velocities are considered as continuous variables. The microscopic details of these particles are not specified. Instead, their time evolution is treated in a simplified form through collective stochastic collisions which preserve momentum and energy. This property allows MPCD to recover, over long simulation times, the hydrodynamic equations of mass, momentum, and heat propagation. In addition, the stochastic character of MPCD produces fluctuations and Brownian forces [7, 8].

A considerable number of soft condensed matter systems have been successfully simulated under MPCD schemes. Recent examples include sediment-water interface flow [9], bistable biochemical systems [10], and two-dimensional one-component plasma [11], to mention but a few. Furthermore, modified MPCD rules have been used to extend simulations towards complex solvents, e.g., viscoelastic fluids [12] and binary mixtures [13]. Very recently, algorithms for simulating nematic liquid crystals (NLCs) using the principles of MPCD have been also proposed [14, 15]. They allow one to reproduce hydrodynamic and elastic characteristics of such phases, still being potentially able to be coupled with microscopic degrees of freedom. These algorithms simulate isotropic-nematic phase transitions, consistent dynamics for annihilation of defects, and molecular reorientation under flow.

Here, basic implementations of MPCD for simple liquids and NLCs are considered. Simulations for equilibrium states are presented, where systems are in contact with the thermal bath based on an Andersen thermostat that keeps them at a fixed temperature. In order to exhibit the ability of MPCD to reproduce hydrodynamic behaviour and to emphasise its stochastic character, attention is focused on the analysis of the spectra of hydrodynamic fluctuations produced by the algorithm in both simple liquids and nematics. There is discussed a very good agreement that exists between correlation functions obtained from the numerical implementation with those derived from linearized hydrodynamic models of liquids and LCs. Following the MPCD model for nematics (MPCD-N) introduced by Shendruk and Yeomans in reference [15], simulated particles are assumed to be slender rods and velocity gradients generate torques on them. On the other hand, the effect of reorientation on flow, commonly referred to as backflow, is incorporated in the velocity update stage through a term that converts the angular momentum generated by reorientation into orbital angular momentum. Numerical experiments are conducted in favour of exploring the resulting effects of this orientation$\leftrightarrow$velocity coupling scheme on the relaxation of hydrodynamic fluctuations. It is found that the influence of flow on orientation is very well accounted for by the linearized form of Jeffery’s shear alignment rule. By contrast, backflow does not produce appreciable effects on the spectra of hydrodynamic correlations as long as its associated simulation parameter is maintained within the range proposed in the original version of the method. Through these results it is shown that MPCD has the potential for simulating hydrodynamic phenomena at the mesoscopic scale when extended to mesomorphic phases.

2 MPCD algorithm for isotropic fluids

Multi-particle collision dynamics simulates fluids composed of point particles. In the simplest case, all these particles have the same mass, $m$. They move within a simulation box with sizes $L_{x}$, $L_{y}$, and $L_{z}$ in the $x$, $y$, and $z$ directions of a Cartesian coordinate system, respectively. In what follows, the total number of simulated particles will be denoted by $N$, whereas positions and velocities of particles will be grouped in the matrices $r_{i,\alpha}$ and $v_{i,\alpha}$, respectively, where indices $i$ and $\alpha$ indicate, respectively, particle number and Cartesian coordinate. Thus, $i=1,2,3,\ldots,N$; and $\alpha=x,y,z$. Quantities $r_{i,\alpha}$ and $v_{i,\alpha}$ are considered to be continuous functions of the simulation time, $t$. Multi-particle collision dynamics promotes the evolution of $r_{i,\alpha}$ and $v_{i,\alpha}$ through a succession of streaming and collision events taking place at regular time intervals of size $\Delta t$. The precise steps of the MPCD algorithm are described below.

2.1 Streaming

For simulations of homogeneous systems, streaming consists of a simple ballistic displacement of the particles during the time-step $\Delta t$. More precisely, given the current state variables, $r_{i,\alpha}(t)$ and $v_{i,\alpha}(t)$, positions at next time-step are given by

[TABLE]

It is worth noting that streaming must be followed by a rule for handling the escape of particles through the boundaries of the simulation box. To approximate the behaviour of macroscopic systems in thermodynamic equilibrium, typical periodic boundary conditions (PBCs) are preferred [3].

2.2 Collision

The physical aim of the collision step is to promote the momentum exchange between particles. Collisions are performed by imagining the simulation box as composed of equal-sized cubic cells. For this purpose, it is considered that $L_{x}$, $L_{y}$, and $L_{z}$, are multiples of a unit length, $a$, i.e., $L_{x}=n_{x}a$, $L_{y}=n_{y}a$, and $L_{z}=n_{z}a$, with $n_{x}$, $n_{y}$, and $n_{z}$ integers. Particles located within the same cell collide. Therefore, the imaginary cells are referred to as collision cells. The number of particles within collision cells could be different and vary as function of time because particles might come in and out from them.

Here, for concreteness, discussion will be focused on the MPCD method based on the application of an Andersen thermostat (MPCD-AT), where relative velocities within collision cells are replaced by new velocities sampled from the Maxwell-Boltzmann distribution at temperature $T$. Two possible versions of MPCD-AT exist differing from one another by their capacity to preserve the angular momentum. For simulations without angular momentum conservation (MPCD-AT$-$a), velocities are updated according to

[TABLE]

where $\bar{v}^{\,\text{c}}_{\alpha}(t)$ is the centre of mass velocity of the cell containing the $i$-th particle, $\xi_{i,\alpha}$ is the assigned stochastic velocity, and $\bar{\xi}^{\,\text{c}}_{\alpha}=\sum_{j\in c}\xi_{j,\alpha}/N^{\text{c}}$. In the last definition, summation extends over all the $N^{\text{c}}$ particles containing the cell where $i$ is located. Notice that hereafter a superscript “c” will be used to indicate properties measured at MPCD collision cells, and should not be confused with a variable that can take numerical values. It can be easily verified that equation (2.2) preserves linear momentum after collision. Nevertheless, it generates the change in angular momentum,

[TABLE]

where $\bar{r}^{\,\text{c}}_{\alpha}$ is the centre of mass position in the cell, $\varepsilon_{\alpha\alpha^{\prime}\beta^{\prime}}$ is the Levi-Civita symbol, and the summation over Greek repeated indices is implied, a convention that will be followed from now on, unless the contrary is explicitly indicated. Simulations with conservation of local angular momentum, MPCD-AT$+$a, are accomplished by subtracting from equation (2.2) the amount of angular momentum arising from stochastic velocity sampling, namely [16]

[TABLE]

where the last term on the right-hand side is responsible for angular momentum conservation. There, $J^{\text{c}}_{\alpha\beta}$ is the inverse of the local moment of inertia tensor.

In both MPCD-AT$-$a and MPCD-AT$+$a, temperature is automatically controlled by the stochastic velocity sampling, so simulations proceed in the canonical ensemble.

2.3 Grid shift

Implementations of MPCD based solely on streaming and collision are not Galilean invariant, an effect that becomes more notorious in simulations conducted at small mean free paths, $l_{\textrm{m}}=\Delta t(k_{\textrm{B}}T/m)^{1/2}$, where $k_{\textrm{B}}$ is the Boltzmann constant. The reason is that for small $l_{\textrm{m}}$, particles collide repeatedly with those in their neighbourhood and become correlated over long time periods. In order to restore Galilean invariance, Ihle and Kroll proposed to perform a stochastic displacement of the grid of collision cells before the momentum exchange event [17]. It is usual to conduct the grid shift independently along each Cartesian axis, using uniform distributions within the range $\left[-a/2,a/2\right]$. The stochastic grid shift helps particles to exchange momentum with a broader set of neighbours and allows the system to recover molecular chaos. In practice, it is easier to program this step by displacing all the particles in the system by the same random vector, applying the corresponding boundary conditions for those particles leaving the simulation box. Once particles have been accommodated, collision takes place. Then, particles are moved back by reversing the stochastic displacement and applying the proper boundary conditions.

3 Hydrodynamic fluctuations in MPCD

3.1 Linearized hydrodynamics

Steps in MPCD are designed to satisfy the mass and momentum conservation. Malevanets and Kapral [1, 2] showed that when these steps are applied over long simulation times, they permit the system to solve the Navier-Stokes equations. Moreover, MPCD is intrinsically stochastic, therefore density and velocity fields suffer from small fluctuations at the local level. Analysing the behaviour of fluctuations produced by the implementation is fundamental for diverse reasons. First, it is important to verify that spontaneous changes in hydrodynamic fields are consistent with a hydrodynamic description, since such changes will be responsible for producing Brownian forces on solute particles [18]. Second, by analysing the spatial and temporal decay of fluctuations it is possible to measure the properties of the simulated systems, namely, viscosity, diffusion, and sound attenuation coefficients [7]. The results of such measurements can be compared against predictions of theoretical treatments of the algorithm to reinforce the understanding of the method and give reliance in using it in more complex situations.

Since MPCD-AT$+$a and MPCD-AT$-$a are intrinsically thermostatted, fluids simulated by them require the specification of two fields only. They are the density and velocity fields, $\rho(\vec{r},t)$ and $V_{\alpha}(\vec{r},t)$, respectively, where $\vec{r}$ is the position vector. These fields obey the continuity and momentum conservation equations,

[TABLE]

[TABLE]

respectively, where $\sigma_{\alpha\beta}=\sigma_{\alpha\beta}(\vec{r},t)$ is the stress tensor. By introducing the thermodynamic pressure, $p(\vec{r},t)$, and assuming a linear dependence of the viscous stress on the velocity gradient $\partial V_{\alpha}/\partial x_{\beta}$, $\sigma_{\alpha\beta}$ can be written in the form

[TABLE]

where $\eta_{\alpha\beta\alpha^{\prime}\beta^{\prime}}$ is the viscous tensor, whose particular form depends on the employed MPCD algorithm. Explicit expressions for viscous contributions in MPCD-AT$\pm$a have been derived, e.g., in reference [19]. In case of no angular momentum conservation, $\sigma_{\alpha\beta}$, reads as [20]

[TABLE]

where viscosity coefficients $\bar{\eta}^{\,\text{kin}}$ and $\bar{\eta}^{\,\text{col}}$ represent contributions due to kinetic and collision processes, respectively. The former comes from the transverse momentum transport produced by motion of particles and is the dominant contribution to the viscosity of gas phases. Collisional viscosity is related with the redistribution of momentum caused by multi-particle collision events. For MPCD-AT$-$a, $\bar{\eta}^{\,\text{kin}}$ and $\bar{\eta}^{\,\text{col}}$ explicitly read as [20]

[TABLE]

[TABLE]

where $\bar{N}^{\textrm{c}}$ represents the average numerical density at collision cells.

For MPCD-AT$+$a, the stress tensor can be written in the form [20]

[TABLE]

where, $\eta=\eta^{\text{kin}}+\eta^{\text{col}}$ with

[TABLE]

[TABLE]

while $\eta^{\text{V}}$, i.e., the bulk viscosity coefficient, reduces to $\eta^{\text{V}}=\bar{\eta}^{\,\text{col}}/3$.

The system of hydrodynamic equations for MPCD-AT$\pm$a can be closed through the thermal equation of state, $p(r,t)=(T,\rho(\vec{r},t))$, a relation that for MPCD fluids has been proven to be of the ideal gas type [18]. Density and velocity fluctuations, $\delta\rho(\vec{r},t)$ and $\delta V_{\alpha}(\vec{r},t)$, are defined through

[TABLE]

and

[TABLE]

respectively, where $\rho^{\textrm{eq}}$ and $V^{\textrm{eq}}_{\alpha}=0$, are the equilibrium density and flow fields.

In addition, it must be considered that due to the stochastic motion of MPCD particles, spontaneous stress may occur at the local level which is not originated by velocity gradients. Stochastic stress is considered to be an additive term to $\sigma_{\alpha\beta}$, $\Sigma_{\alpha\beta}$, which as customary is approximated as a Gaussian-Markov process (white noise). The fluctuation-dissipation theorem (FDT) relates two stochastic stress components at two different times and positions as

[TABLE]

where $\delta\left(\vec{r}^{\,\prime}-\vec{r}\right)$ and $\delta\left(t^{\prime}-t\right)$ are Dirac delta functions in space and time domains, respectively.

Fluctuating fields are more conveniently studied in Fourier space. Hereafter, Fourier transforms are to be indicated by a tilde over the corresponding function. In Fourier domain, the FDT for stochastic stress reads as

[TABLE]

By introducing unit vectors $\hat{e}_{1}=\vec{k}/k$, $\hat{e}_{2}=k\left[\hat{e}_{z}-\left(\hat{e}_{1}\cdot\hat{e}_{z}\right)\hat{e}_{1}\right]/(k_{x}^{2}+k_{y}^{2})^{1/2}$, and $\hat{e}_{3}=\hat{e}_{1}\times\hat{e}_{2}$, and projections $\delta\tilde{V}_{\alpha}=\hat{e}_{\alpha}\cdot\delta\tilde{\vec{V}}$, hydrodynamic fluctuations can be split into independent variables $\delta\tilde{V}_{2}$, $\tilde{V}_{3}$, and $\left\{\delta\tilde{\rho},\delta\tilde{V}_{1}\right\}$. Dynamic equations for these quantities are obtained after replacing equations (3.10) and (3.11) into equations (3.1) and (3.2), retaining only linear terms in fluctuations, and incorporating a stochastic stress. Such equations explicitly read as

[TABLE]

for the transverse velocity components, $\alpha=2,3$; and

[TABLE]

for the so-called longitudinal fluctuations. In equations (3.14) and (3.15), $\nu=(\eta^{\text{kin}}+\eta^{\text{col}})/\rho_{0}$ is the kinematic viscosity of the fluid; $c_{\textrm{T}}=[(\partial p/\partial\rho)_{T}]^{1/2}$ is the isothermal sound speed; and $D_{\textrm{l}}$ is the so-called longitudinal viscosity. For MPCD-AT$-$a $D_{\textrm{l}}=(4\eta^{\text{kin}}/3+\eta^{\text{col}})/\rho_{0}$, while for MPCD-AT$+$a one has $D_{\textrm{l}}=(4\eta/3+\eta^{\text{V}})/\rho_{0}$.

Autocorrelation functions (ACFs) between fluctuations in Fourier coordinates are defined as

[TABLE]

and

[TABLE]

where no summation over repeated indices is implied. Explicit expressions for $C_{\rho}(\vec{k},\omega)$ and $C_{\text{V}\alpha}(\vec{k},\omega)$ can be obtained by solving equations (3.14)–(3.15) and using the FDT in Fourier space, equation (3.13). This procedure yields the classical expressions

[TABLE]

[TABLE]

[TABLE]

In equations (3.18)–(3.20), $q(k)$ has been defined as $q(k)=c_{\textrm{T}}k[1-(D_{\textrm{l}}k/c)^{2}]^{1/2}$, and $A$ is a constant defined by $A=2(2\piup)^{4}k_{\text{B}}T\delta(\vec{k}=0)\delta(\omega=0)$.

It can be noticed that $C_{\rho}$ consist of two symmetric peaks, commonly referred to as the Brillouin peaks, centred at $\pm q(k)$. In MPCD-AT$\pm$a, the classical scattering spectrum lacks the central (Rayleigh) peak because the algorithm is thermostatted and no temperature fluctuations are sustained. In more general cases, temperature fluctuations in fluids cannot be neglected. They relax towards equilibrium following a diffusion type equation with an extra term that couples them with longitudinal velocity components [7]. Temperature fluctuations produce the previously mentioned Rayleigh peak and have a measurable effect on the propagating modes [21]. In the presence of temperature fluctuations, Brillouin peaks are centered at the isentropic sound speed instead of being located at $c_{\textrm{T}}$. It is worth noting that some MPCD variations could be used to incorporate temperature fluctuations. With this aim, a natural alternative to use is stochastic rotation dynamics [1, 2], the first MPCD algorithm, because it incorporates mesoscopic heat flow. Stochastic rotation dynamics could be used at the momentum exchange step instead of the collision rule based on the Andersen thermostat and it can be even modulated to give conditions ranging from adiabatic to isothermal [7].

It can be further observed that function $C_{\textrm{V1}}$ is also double peaked. However, due to the factor $\omega^{2}$ in front of it, $C_{\textrm{V1}}$ vanishes at $\omega=0$. Correlation functions $C_{\textrm{V2}}$ and $C_{\textrm{V3}}$ consist of a single Lorentzian peak at $\omega=0$.

3.2 Measurements

In practice, measurements of ACFs are conducted by recording time-series of spatial Fourier transforms of density and velocity fluctuations, namely

[TABLE]

and

[TABLE]

where $t_{l}$ indicates the simulation time and the over-bar transformation over spatial domain only. Then, ACFs at two different times are calculated using the time-origin moving scheme, e.g.,

[TABLE]

where $N_{\text{s}}$ denotes the size of the recorded time-series. Finally, discrete Fourier transformation over time domain can be applied to get ACFs in $\{\vec{k},\omega\}$ coordinates.

It should be stressed that due to the use of PBC, wave vector components are restricted to be multiples of the inverse system size, $2\piup/L$. In this way, summations in equations (3.21) and (3.22) automatically take care of minimum images of simulated particles.

Numerical experiments return ACF with the mathematical features summarised by equations (3.18)–(3.20), which can be supplemented with an expression for viscosity coefficients, equations (3.5)–(3.9), to verify the correctness of the numerical implementation. Fixed units of time ($\Delta t=1.0$), energy ($k_{\textrm{B}}T=1.0$), and mass ($m=1.0$) were used during simulations. Units of length were derived from $a=\Delta t\sqrt{k_{\textrm{B}}T/m}$. Throughout the following, quantities will be reported in these simulation units (s.u.). Simulated systems had $L_{x}=L_{y}=L_{z}=L=20\,a$, and $\bar{N}^{\textrm{c}}=20$. Scattering geometry was defined by the wave vector $\vec{k}=(2\piup/L,0,0)$. In correlations obtained from MPCD-AT$\pm$a implementations, constant $A$ in equations (3.18)–(3.20) was considered an adjustable parameter. Reduced functions defined by $G_{\rho}=C_{\rho}/A$, and $G_{\textrm{V}\alpha}=C_{\textrm{V}\alpha}/A$, are compared against their theoretical counterparts obtained from equations (3.18)–(3.20) and (3.5)–(3.9) in figures 1–3. They show the ability of MPCD to support hydrodynamic fluctuations that relax towards equilibrium as fluctuations do in a fluid kept in contact with a thermal bath. The good agreement found between numerical and theoretical results exhibits a very satisfactory understanding about MPCD that makes it a reliable method for simulations of soft condensed matter systems.

4 MPCD for nematic phases

Liquid crystals are matter phases characterised by possessing a structural order lower than the one present in crystals, but larger than the corresponding to ordinary isotropic liquids [22]. They have been subject of continuous interest through decades since they exhibit anisotropic optical properties that can be manipulated by small energies of electromagnetic or surface-liquid interaction origin [23]. Modern LC technologies have promising applications in colloid manipulation, detection of biological agents, and medical diagnosing [24, 25, 26, 27]. Numerical simulations could play an important role in helping developers to test models of these emerging technologies [28].

In nature, NLCs are constituted by anisotropic molecules of elongated or discotic types, commonly referred to as nematogens. Accordingly, in order to generalise MPCD to nematics, simulated particles are supplied with an orientation degree of freedom $u_{i,\alpha}=u_{i,\alpha}(t)$, which is a unit vector. Furthermore, MPCD rules described in section 2 are augmented to cover three main aspects, namely: multi-particle collision events cause changes in $u_{i,\alpha}$; velocity gradients produce torques on $u_{i,\alpha}$; and reorientation induces flow (backflow effect). Implementation of these effects is described below.

4.1 Orientation exchange

Orientation exchange takes into account the interaction of nematogens with their surroundings in a coarse-grained fashion. More precisely, order in NLCs is measured through the quadrupolar moment of the orientation distribution, also referred to as the order parameter tensor. This quantity is estimated in MPCD-N at the cell level through

[TABLE]

where $N_{\textrm{c}}$ is the number of particles in the cell after streaming and grid shift. For small-size nematic systems far away from the isotropic-nematic transition, the amount of the order in the cell can be quantified by the largest eigenvalue of $Q^{\textrm{c}}_{\alpha\beta}$, called the scalar order parameter, $S^{\textrm{c}}$. The unit eigenvector associated to $S^{\textrm{c}}$, $n^{\textrm{c}}_{\alpha}$, is called the director and represents the average molecular orientation. If nematogen $i$ is found in a cell with the order parameter $S^{\textrm{c}}$ and director $n^{\textrm{c}}_{\alpha}$, it is assumed that it interacts with nematogens within the same cell according to the (mean-field) Maier-Saupe potential energy [29],

[TABLE]

where $U$ represents the strength of the mean-field interaction. Accordingly, an orientation collision operator can be proposed that updates orientations by sampling new vectors, $u_{i,\alpha}^{\prime}$, from the canonical distribution $P=Z^{-1}\exp\left[-U_{\text{mf}}/(k_{\text{B}}T)\right]\sin\theta$, where $Z^{-1}$ is the normalisation constant and $\theta$ is the angle between $u_{i,\alpha}$ and $n^{\textrm{c}}_{\alpha}$.

4.2 Reorientation by flow

Flow effects on $u_{i,\alpha}$ are modelled as if nematogens were slender rods that experience torques induced by local velocity gradients. The velocity gradient tensor at a given collision cell is measured by considering the difference between the centre of mass velocities of its neighbouring planes. Then, nematogens are reoriented following the rule [30]

[TABLE]

where $\Omega^{\textrm{c}}_{\alpha\beta}$ and $D^{\textrm{c}}_{\alpha\beta}$ stand for the antisymmetric and symmetric parts of $\partial v_{\alpha}/\partial x_{\beta}$, respectively. Equation (4.3) introduces quantities $\lambda$ and $\chi_{\text{HI}}$. The former is the tumbling parameter. It determines the tendency of $u_{i,\alpha}$ to be aligned by shear at a stationary angle ($\lambda>1$), or to continuously rotate in a direction consistent with the vorticity of the flow ($\lambda<1$). During simulations, $\lambda$ can be fixed to produce the desired dynamic effect [15]. Parameter $\chi_{\text{HI}}$ is an auxiliary one. It is initialised within the range $[0,1]$ to control the overall flow-orientation coupling. By setting $\chi_{\text{HI}}=0$, such a coupling vanishes, while the maximum reorientation by flow is simulated for $\chi_{\text{HI}}=1$.

4.3 Backflow

During the orientation update stage, cells receive an effective angular velocity

[TABLE]

Under the assumption of over-damped rotational dynamics with viscosity $\gamma_{\text{R}}$, this implies a net torque, $\Gamma^{\textrm{c}}_{\alpha}=-\gamma_{\text{R}}w^{\textrm{c}}_{\alpha}$, and the production of an angular momentum in the cell $\Delta{\mathcal{L}}^{\textrm{c}}_{\alpha}=\Gamma^{\textrm{c}}_{\alpha}\Delta t$.

Backflow is taken into account by adding this angular momentum contribution to $\Delta L^{\textrm{c}}_{\alpha}$ in the MPCD-AT$+$a collision rule, equation (2.4).

5 Fluctuating nematodynamics

In MPCD-N, in addition to density and velocity fluctuations discussed in section 3, there occur spontaneous changes in the orientation field. Orientation fluctuations are defined through

[TABLE]

where $n^{\textrm{eq}}_{\alpha}$ is the global director field, which is assumed to be fixed in space. Thermal director fluctuations can be analysed by extending the linearized scheme of section 3. Specifically, following general hydrodynamic models of NLCs [31], the relaxation equation for the director field is proposed to have the form

[TABLE]

where $Y_{\alpha}$ is sometimes termed a quasi-current due to the fact that the director is not a conserved quantity and, consequently, the integral of $Y_{\alpha}$ over a surface cannot be a flux [31]. In equation (5.2), $\Upsilon_{\alpha}$ is a stochastic contribution to $Y_{\alpha}$, which can be assumed to be a Gaussian-Markov process with zero mean, strength $\gamma$, and FDT

[TABLE]

where matrix $\delta_{\alpha\beta}-n_{\alpha}n_{\beta}$ projects vectors on a direction perpendicular to $n_{\alpha}$. This permits to satisfy the condition $n_{\alpha}\delta n_{\alpha}=0$, imposed by the normalisation of $n_{\alpha}$.

In MPCD-N, time evolution of orientations is promoted by shear, equation (4.3). Furthermore, streaming and collision dynamics could be anticipated to give rise to orientation diffusion, then $Y_{\alpha}$ is written as

[TABLE]

where the orientation diffusion coefficient is $D_{n}$ and $\Omega_{\alpha\beta}$ and $D_{\alpha\beta}$ are the antisymmetric and symmetric parts of $\partial V_{\alpha}/\partial x_{\beta}$.

Dynamic equations for director fluctuations can be obtained by replacing equations (3.11) and (5.1) into equations (5.2) and (5.4), and retaining only linear contributions. In a first approximation, mass and momentum transport can be assumed to be still given by equations (3.1) and (3.2), thus neglecting the effects of backflow [32]. In section 5.1, the validity of this assumption will be explored. It will be shown to be a very satisfactory approximation for those values $\gamma_{\textrm{R}}$ used up to now in previous implementations of MPCD-N. As in the case of simple fluids, the resulting system is more conveniently studied in Fourier space. In addition, separation of transverse and longitudinal variables is worth recommending. Taking into account nematic symmetry, basis vectors are defined by $\hat{e}_{1}=\vec{k}/k$, $\hat{e}_{2}=(\hat{n}^{\textrm{eq}}\times\hat{e}_{1})/\sin\theta^{\textrm{eq}}$, and $\hat{e}_{3}=\hat{e}_{1}\times\hat{e}_{2}$, where $\theta^{\textrm{eq}}$ is the angle between $\hat{n}^{\textrm{eq}}$ and $\vec{k}$. For concreteness, the scattering geometry defined by $\vec{k}=(k,0,0)$, and $\hat{n}^{\textrm{eq}}=(0,0,1)$ will be considered hereafter. Then, the following system of independent equations is obtained for fluctuating fields in MPCD-N,

[TABLE]

[TABLE]

and

[TABLE]

It can be observed from equations (5.5)–(5.7) that in the present nematodynamic model, density and velocity fluctuations are not perturbed by orientation fluctuations though fluctuations $\delta\tilde{V}_{3}$ affect the dynamics of $\delta\tilde{n}_{1}$. Consequently, it is expected that density and velocity ACFs in MPCD and MPCD-N will be the same.

Orientation fluctuations can be easily solved from equations (5.5) and (5.7). Such a solution, together with the FDT in Fourier space,

[TABLE]

can be used to obtain the following orientation ACFs

[TABLE]

[TABLE]

Consequently, $C_{\textrm{n2}}$ is not expected to be modified by velocity-orientation coupling, while the intensity of $C_{\textrm{n1}}$ is anticipated to increase with a term proportional to $\chi_{\text{HI}}^{2}(\lambda-1)^{2}$.

Another correlation function that could bring information about the correct coupling between velocity and orientation fluctuations is the cross correlation function, $C_{\textrm{n1,V3}}$, given by

[TABLE]

5.1 Measurements

Measurements of correlations $C_{\textrm{n1}}$, $C_{\textrm{n2}}$, and $C_{\textrm{n1,V3}}$ can be carried out by introducing the time-series of spatial Fourier transform of orientation fluctuations

[TABLE]

Then, the same procedure used in section 3.2 to measure density and velocity ACFs can be followed. It is worth stressing that in order to make these measurements compatible with the use of PBC and with the model yielding equations (5.9)–(5.11), it is necessary to ensure that global director will remain fixed during the computation stage. This can be done by introducing an additional step in the MPCD-N algorithm where all nematogens are rotated by the same angle in order to align $\hat{n}^{\textrm{eq}}$ towards the $z$ [33]. This direction is considered as the $\hat{e}_{3}$ direction, while $x$ and $y$ are considered to coincide with $\hat{e}_{1}$ and $\hat{e}_{2}$, respectively.

Numerical experiments were conducted using the same simulation parameters $\Delta t$, $k_{\textrm{B}}T$, $m$, $\bar{N}^{\textrm{c}}$, $a$, and $L$, given in section 3.2. Nematic order was simulated using $U=20k_{\text{B}}T$. For this mean-field strength, it has been shown that the average scalar order parameter produced in collision cells, $\bar{S}^{\textrm{c}}=0.947$, is in agreement with predictions of self-consistent models of NLCs [32]. Two main features of hydrodynamic correlation functions in MPCD-N will be explored. Namely, the effects of backflow controlled by the parameter $\gamma_{\textrm{R}}$; and the effects of velocity-orientation coupling controlled by parameters $\chi_{\text{HI}}$ and $\lambda$.

5.1.1 Backflow effects

With the purpose of analysing the extent at which the backflow perturbs hydrodynamic fluctuations, $\gamma_{\textrm{R}}$ was varied over the range $\left[0.01,1.0\right]$, which extends by one order of magnitude the values considered before the present contribution. Eighteen experiments were conducted using combinations between parameters $\lambda=0.5,2.0$; $\chi_{\textrm{HI}}=0.0,0.5,1.0$, and $\gamma_{\text{R}}=0.01,0.1,1.0$. Correlation functions $C_{\rho}$ and $C_{\textrm{V}\alpha}$ were measured and compared with those obtained from simple MPCD-AT$+$a. No systematic differences were observed between correlations $C_{\rho}$, $C_{\textrm{V1}}$ and $C_{\textrm{V2}}$ in MPCD-N and those measured in MPCD-AT$+$a, in good agreement with equations (5.5) and (5.6). A similar conclusion applies for function $C_{\text{V3}}$ for small and moderate values of the product $\gamma_{\textrm{R}}\chi_{\textrm{HI}}$, namely, for $\gamma_{\textrm{R}}\chi_{\textrm{HI}}<0.5$. Nevertheless, backflow effects become noticeable, increasing the intensity of $C_{\textrm{V3}}$ for $\lambda=0.5$, and decreasing it for $\lambda=2.0$. This is illustrated in figure 4, where the results for $C_{\text{V3}}$ obtained for the whole experimental setup are presented. Therefore, it can be concluded that if $\gamma_{\textrm{R}}\leqslant 0.1$, as it has been used in references [15, 32], hydrodynamic fluctuations of flow and density in MPCD-N can be considered to be independent of orientation fluctuations.

5.1.2 Velocity-orientation coupling

In order to analyse the dynamics of director fluctuations, reduced correlations are introduced as follows:

[TABLE]

[TABLE]

and

[TABLE]

Numerical results are normalised using the value of $A$ obtained in simulations of isotropic phases.

It is worth stressing that, up to the present, there are no analytical treatments of MPCD-N from which values of $\gamma$ and $D_{\textrm{n}}$ could be inferred in terms of simulation parameters, $\bar{N}^{\textrm{c}}$, $\Delta t$, $m$, etc. Consequently, the subsequent discussion is conducted in a semiempirical form, using values for $\gamma$ and $D_{\textrm{n}}$ derived from the numerical experiments. Correlations $G_{\textrm{n2}}$ and $G_{\textrm{n1}}$ showed the Lorentzian shape predicted from equations (5.13) and (5.14) with estimates, in s.u., $\gamma_{\text{R}}=5.5\times 10^{-4}$ and $D_{\text{n}}=0.615$. $G_{\textrm{n2}}$ was found to be independent of simulation parameters, as expected from equation (5.13), except for simulations conducted at $\gamma_{\text{R}}=1.0$, $\lambda=2.0$, and $\chi_{\textrm{HI}}=1.0$. In this case, $G_{\textrm{n2}}$ is more intense and exhibits two small lateral peaks located at $\pm q(k)$, indicating an apparent coupling with propagating modes. These features are illustrated in figure 5. This behaviour can be explained by noting that velocity fluctuations play a destructive role on the orientation order, since they increase the strength of director fluctuations [15]. At large values of $\lambda$ and $\chi_{\textrm{HI}}$, this effect could be large enough to produce a reduction of global order in the system. Specifically, for $\lambda=2.0$ and $\chi_{\text{HI}}=1.0$, the average global order in the sample was found to be $S=0.74$, while in all remaining cases it was found to be close to $S=0.94$. This is schematically illustrated in figure 6 where snapshots of systems simulated at $\{\lambda=0.5,\gamma_{\textrm{R}}=1.0,\chi_{\text{HI}}=1.0\}$ and $\{\lambda=2.0,\gamma_{\textrm{R}}=1.0,\chi_{\text{HI}}=1.0\}$ are shown. Large director gradients can be observed in the latter case, for which it could be expected that significant changes of the global director will be produced between simulation steps. Then, the approximation of fixed $\hat{n}^{\textrm{eq}}$, from which separation of longitudinal and transverse variables is constructed, becomes weaker. This implies that in the case $\{\lambda=2.0,\gamma_{\textrm{R}}=1.0,\chi_{\text{HI}}=1.0\}$, fluctuation $\delta\tilde{n}_{y}$ that is used to estimate $G_{\textrm{n2}}$, could be coupled with the propagating velocity components, a situation not expected for the transverse variable $\delta\tilde{n}_{2}$.

Correlation $G_{\textrm{n1}}$ was found to exhibit a systematic dependence on the product $(\lambda-1)^{2}\chi_{\text{HI}}^{2}$. Results are in good agreement with equation (5.14), as it can be observed in figure 7. Analogous results were found in the case of the cross correlation $G_{\textrm{n1,V3}}$, whose real part obtained from equation (5.15) is illustrated in figure 8 and is compared with that obtained from numerical experiments. Thus, it is shown that the coupling between orientation and velocity fields predicted from the linearized model of nematic phases is well reproduced by MPCD-N. Again, systems simulated at $\gamma_{\textrm{R}}=1.0$, and $\chi_{\text{HI}}=1.0$, seem to deviate from this rule. In the case $\lambda=0.5$, $G_{\textrm{n1}}$ and $G_{\textrm{n1,V3}}$ had a larger height than the one predicted by the linearized model. On the contrary, for $\lambda=2.0$, both $G_{\textrm{n1}}$ and $G_{\textrm{n1,V3}}$ decreased their maximum height with respect to the value anticipated by the theoretical model. This discrepancy is related with the behaviour of the velocity ACF $C_{\textrm{V3}}$ already presented in figure 4. To be specific, the increment (decrement) of $C_{\textrm{V3}}$ for $\lambda=0.5$ ($\lambda=2.0$) causes a corresponding increment (decrement) of $G_{\textrm{n1}}$ and $G_{\textrm{n1,V3}}$ through the hydrodynamic coupling in equation (5.7).

6 Conclusions

The capacity of MPCD-AT for reproducing collective phenomena of the hydrodynamic type has been analysed by exploring its spectra of thermal fluctuations in implementations for simple fluids and LCs. First, it was illustrated that relaxation of spontaneous fluctuations in simple fluids simulated under MPCD-AT$\pm$a rules can be described by linearized hydrodynamics. Analytical expressions of viscosity coefficients have been used for comparison purposes, in order to exhibit a good agreement existing between theoretical models and numerical implementations of MPCD-AT$\pm$a.

Afterwards, an extension of MPCD towards nematic phases has been presented. This implementation was based on the original one introduced in reference [15]. Considering the fundamental physical and mathematical aspects involved in the study of LCs, as well as the technological relevance of these materials, such an extension could be very fruitful in the near future [28]. The analysis of correlation functions in MPCD-N was performed in a semiempirical form, since there are no explicit expressions for this method relating transport coefficients and simulation parameters. Numerical experiments showed that hydrodynamic fluctuations in MPCD-N can be described by a simplified linearized model where viscosity and elastic effects are isotropic. In addition, if simulation parameters are kept within the ranges of the original MPCD-N proposal [15], velocity and orientation fluctuations are only one-way coupled. Namely, while velocity fluctuations affect the orientation fluctuations, backflow can be neglected in the linearized regime. Nonetheless, for large values of the parameters that couple the flow and orientation fields, $\gamma_{\textrm{R}}$ and $\chi_{\textrm{HI}}$, it was demonstrated that the velocity-orientation interaction can be strong enough to cause changes in the global orientation state. Though this effect has been already discussed in reference [15], in the context of isotropic-nematic phase transition in MPCD-N, its implications on the measurements of hydrodynamic correlations were exhibited here for the first time. In this respect, it was concluded that large shear-induced reorientations could make the linearized model incompatible with correlation measurements.

Numerical experiments were carried out in small-sized systems. Ongoing tests are in progress for systems with larger sizes. Their purpose is to assess the importance of size-dependence on the previously mentioned effects. Similarly, activities are currently under way that have the purpose of bringing MPCD-N properties closer to those observed in real nematics. They include the introduction of director-dependent collision rules to simulate systems with anisotropic viscosity and elasticity, as well as the use of adapted mean field potentials in order to simulate smectic phases.

Acknowledgements

Author thanks La Salle University Mexico for final support under grant NEC-$08/18$.