Hydrodynamic correlations in isotropic fluids and liquid crystals simulated by multi-particle collision dynamics
H. Híjar


Optimal Policies for Convex Symmetric Stochastic Dynamic Teams and their Mean-field Limit
This research was supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada. A summary of some of the results here is included in <cit.> presented at the 2019 IEEE Conference
on Decision and Control (CDC). The authors are with the Department of Mathematics and Statistics,
     Queen's University, Kingston, ON, Canada.
     Email: {16sss3,yuksel@queensu.ca}
    Sina Sanjari and Serdar Yüksel
    
===========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





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.


ulti-scale simulation techniques, particle-based simulations of 
          fluids, thermal fluctuations, nematic liquid crystals, structure 
	  of liquids and liquid crystals

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





§ 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 <cit.>. 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 <cit.>, 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 <cit.>.

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 <cit.> and the lattice Boltzmann method <cit.> 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 <cit.>, 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 <cit.>.

A considerable number of soft condensed matter systems have been successfully 
simulated under MPCD schemes. Recent examples include sediment-water interface 
flow <cit.>, bistable biochemical systems <cit.>, and two-dimensional 
one-component plasma <cit.>, to mention but a few. Furthermore, modified 
MPCD rules have been used to extend simulations towards complex solvents, e.g., 
viscoelastic fluids <cit.> and binary mixtures <cit.>. Very recently,
algorithms for simulating nematic liquid crystals (NLCs) using the principles of 
MPCD have been also proposed <cit.>. 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 <cit.>, 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↔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. 




§ 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,α and v_i,α, respectively, where indices i and α 
indicate, respectively, particle number and Cartesian coordinate. Thus,
i = 1, 2, 3, …, N; and α = x,y,z. Quantities r_i,α and 
v_i,α are considered to be continuous functions of the simulation time, 
t. Multi-particle collision dynamics promotes the evolution of r_i,α 
and v_i,α through a succession of streaming and collision events taking 
place at regular time intervals of size Δ t. 
The precise steps of the MPCD algorithm are described below.




 §.§ Streaming




For simulations of homogeneous systems, streaming consists of a simple ballistic 
displacement of the particles during the time-step Δ t. More precisely, 
given the current state variables, r_i,α(t) and v_i,α(t), 
positions at next time-step are given by

    r_i,α(t+Δ t) = r_i,α(t) + v_i,α(t) Δ t.


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 <cit.>. 




 §.§ 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_xa, 
L_y = n_y a, and L_z = n_za, 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

    v_i,α (t+Δ t) = v̅^ c_α(t)
    			  + ξ_i,α - ξ̅^ c_α ,

where v̅^ c_α(t) is the centre of mass velocity of the 
cell containing the i-th particle, ξ_i,α is the assigned stochastic 
velocity, and 
ξ̅^ c_α = ∑_j∈ cξ_j,α/N^c. In the 
last definition, summation extends over all the N^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 (<ref>) preserves linear momentum after 
collision. Nevertheless, it generates the change in angular momentum,

    Δ L^c_α 
    = ∑_j ∈ cε_αα^'β^'(r_j,α^'-r̅^ c_α^')
     (v_j,β^'-ξ_j,β^'),

where r̅^ c_α is the centre of mass position in the cell, 
ε_αα^'β^' 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 (<ref>) the amount of 
angular momentum arising from stochastic velocity sampling, namely <cit.>

    v_i,α(t + Δ t) = v̅^ c_α(t) + ξ_i,α 
                               - ξ̅^ c_α 
    			   -ε_αα^'β^'
    			    J^c_α^'α^'' (t) 
                                Δ L^ c_α^''(t)
    			    ( r_i,β^' - r̅^ c_β^'),

where the last term on the right-hand side is responsible for angular 
momentum conservation. There, J^c_αβ 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.




 §.§ 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_m = Δ t(k_BT/m)^1/2, 
where k_B is the Boltzmann constant. The reason is that for small 
l_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 <cit.>. It is usual to conduct the grid shift independently along each 
Cartesian axis, using uniform distributions within the range 
[-a/2,a/2]. 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.





§ HYDRODYNAMIC FLUCTUATIONS IN MPCD







 §.§ Linearized hydrodynamics




Steps in MPCD are designed to satisfy the mass and momentum conservation. Malevanets 
and Kapral <cit.> 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 <cit.>. 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 <cit.>. 
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, ρ(r⃗,t) and V_α(r⃗,t),
respectively, where r⃗ is the position vector. These fields obey the 
continuity and momentum conservation equations,

    ∂ρ/∂ t 
    = -∂/∂ x_α( ρ V_α) ,


    ρ( 
         ∂/∂ t 
         + V_β∂/∂ x_β) 
         V_α 
    = - ∂/∂ x_βσ_αβ ,

respectively, where σ_αβ = σ_αβ (r⃗,t) 
is the stress tensor. 
By introducing the thermodynamic pressure, p(r⃗,t), and assuming a linear 
dependence of the viscous stress on the velocity gradient 
∂ V_α/∂ x_β, σ_αβ can 
be written in the form

    σ_αβ = -p δ_αβ 
                          + η_αβα^'β^'∂ V_β^'/∂ x_α^'  ,

where η_αβα^'β^' is the viscous 
tensor, whose particular form depends on the employed MPCD algorithm. Explicit 
expressions for viscous contributions in MPCD-AT±a have been derived, e.g., 
in reference <cit.>. In case of no angular momentum conservation,  
σ_αβ, reads as <cit.>

    σ_αβ = -p δ_αβ 
    		      + η̅^ kin( ∂ V_β/∂ x_α
                                +∂ V_α/∂ x_β
                                -2/3δ_αβ∂ V_α^'/∂ x_α^')
    		      +η̅^ col∂ V_α/∂ x_β ,

where viscosity coefficients η̅^ kin and η̅^ 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, 
η̅^ kin and η̅^ col explicitly read 
as <cit.>

    η̅^ kin = N̅^c k_B T Δ t/a^3( N̅^c/N̅^c-1 + ^-N̅^c
    		          -1/2),


    η̅^ col = m/12 a Δ t(N̅^c-1 + ^-N̅^c) ,

where N̅^c represents the average numerical density at collision cells.

For MPCD-AT+a, the stress tensor can be written in the form <cit.> 

    σ_αβ = -p δ_αβ 
                          +η( ∂ V_β/∂ x_α
                                +∂ V_α/∂ x_β
                                -2/3δ_αβ∂ V_α^'/∂ x_α^')
    		      +η^Vδ_αβ∂ V_α^'/∂ x_α^' ,

where, η = η^kin + η^col with

    η^kin = k_B T N̅^cΔ t/a^3(N̅^c/N̅^c - 5/4 -1/2),


    η^col = m N̅^c/24 a Δ t( 1 - 7/5 N̅^c),

while η^V, i.e., the bulk viscosity coefficient, reduces to 
η^V = η̅^ col/3.

The system of hydrodynamic equations for MPCD-AT±a can be closed through the 
thermal equation of state, p(r,t) = (T,ρ(r⃗,t)), a relation that for 
MPCD fluids has been proven to be of the ideal gas type <cit.>. 
Density and velocity fluctuations, δρ(r⃗,t) and 
δ V_α(r⃗,t), are defined through

    ρ(r⃗,t) = ρ^eq + δρ(r⃗,t),

and 

    V_α(r⃗,t) = V^eq_α + δ V_α (r⃗,t),

respectively, where ρ^eq and V^eq_α = 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 σ_αβ, Σ_αβ, 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

    ⟨Σ_αβ(r⃗^ ',t^') 
    	     Σ_α^'β^'(r⃗,t) 
    ⟩ 
    = 2 k_B T η_αβα^'β^'δ(r⃗^ '- r⃗)
    δ( t^' - t ),

where δ(r⃗^ '- r⃗) and 
δ( t^'-t ) 
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

    ⟨Σ̃_αβ(k⃗^ ',ω^') 
    Σ̃_α^'β^'(k⃗,ω) 
    ⟩ 
    = 2 (2)^4 k_B T η_αβα^'β^'δ(k⃗^ '+k⃗)
    δ( ω^'+ω).


By introducing unit vectors ê_1 = k⃗/k, 
ê_2 = k [ ê_z - (ê_1·ê_z) 
ê_1]/(k_x^2+k_y^2)^1/2, and 
ê_3 = ê_1×ê_2, and projections 
δṼ_α = ê_α·δṼ⃗̃,
hydrodynamic fluctuations can be split into independent variables
δṼ_2, Ṽ_3, and 
{δρ̃, δṼ_1}. 
Dynamic equations for these quantities are 
obtained after replacing equations (<ref>) and 
(<ref>) into equations (<ref>) 
and (<ref>), retaining only linear terms in fluctuations, and  
incorporating a stochastic stress. Such equations explicitly read as

    (-ω  + ν k^2) δṼ_α  
    =  k/ρ_0Σ̃_1α

for the transverse velocity components, α = 2,3; and

    (
    	-ω    - k ρ_0
    
    	- k c_T^2ρ_0^-1    -ω  + D_l k^2) ·(δρ̃
    δṼ_1) =
             k/ρ_0(
          0 
    Σ̃_11)

for the so-called longitudinal fluctuations. In equations (<ref>) and 
(<ref>), ν = (η^kin+η^col)/ρ_0 
is the kinematic viscosity of the fluid; 
c_T = [(∂ p/∂ρ)_T]^1/2 is the isothermal sound speed; 
and D_l is the so-called longitudinal viscosity. For MPCD-AT-a 
D_l = (4η^kin/3 + η^col)/ρ_0, while for 
MPCD-AT+a one has D_l = (4η/3 + η^V)/ρ_0.

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

    C_ρ(k⃗,ω) = ⟨δρ̃(k⃗,ω) 
                               δρ̃^*(k⃗,ω) 
                               ⟩,

and

    C_Vα(k⃗,ω) = ⟨δṼ_α(k⃗,ω) 
                                  δṼ_α^*(k⃗,ω) 
                                  ⟩,

where no summation over repeated indices is implied. Explicit expressions for 
C_ρ(k⃗,ω) and C_Vα(k⃗,ω) can be obtained by 
solving equations (<ref>)–(<ref>) 
and using the FDT in Fourier space, equation (<ref>). This 
procedure yields the classical expressions

    C_ρ(k⃗,ω) 
    = A ρ_0 D_l k^4/q^2(k) + D_l^2k^4{1-ω/[2q(k)]/[ω-q(k)]^2 + D_l^2k^4
    +1+ω/[2q(k)]/[ω+q(k)]^2 + D_l^2k^4},


    C_V1(k⃗,ω) 
    = ω^2/ρ_0^2k^2 
      C_ρ(k⃗,ω),


    C_V2(k⃗,ω) 
    = C_V3(k⃗,ω) 
    = A 1/ρ_0ν k^2/ω^2 + ν^2 k^4.


In equations (<ref>)–(<ref>), q(k) has been 
defined as q(k) = c_T k [1-(D_lk/c)^2]^1/2, and 
A is a constant defined by 
A = 2 (2 )^4 k_B T δ(k⃗ = 0) δ(ω = 0). 

It can be noticed that C_ρ consist of two symmetric peaks, commonly 
referred to as the Brillouin peaks, centred at ± q(k). In MPCD-AT±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 <cit.>. 
Temperature fluctuations produce the previously mentioned Rayleigh peak and have 
a measurable effect on the propagating modes <cit.>. In the presence of 
temperature fluctuations, Brillouin peaks are centered at the isentropic sound 
speed instead of being located at c_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 <cit.>, 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 <cit.>.

It can be further observed that function C_V1 is also double 
peaked. However, due to the factor ω^2 in front of it, C_V1 
vanishes at ω = 0. Correlation functions C_V2 and 
C_V3 consist of a single Lorentzian peak at ω = 0. 




 §.§ Measurements




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

    δρ̅(k⃗ , t_l) 
    = ρ^eq∑_j=1^Nexp[-k⃗·r⃗_j(t_l) ],

and

    δV̅_α(k⃗ , t_l) 
    = ∑_j=1^N v_j,αexp[-k⃗·r⃗_j(t_l) ],

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.,

    C_ρ(k⃗, t_l+m - t_l) 
    = 1/N_s - m +1∑_l = 0^N_s-mδρ̅(k⃗, t_l+m)   
                                    δρ̅(k⃗, t_l),

where N_s denotes the size of the recorded time-series. Finally, 
discrete Fourier transformation over time domain can be applied to get ACFs in
{k⃗,ω} 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  / L. In this way, 
summations in equations (<ref>) and 
(<ref>) automatically take care of minimum images of 
simulated particles.





Numerical experiments return ACF with the mathematical features summarised by 
equations (<ref>)–(<ref>), which can be 
supplemented with an expression for viscosity coefficients, 
equations (<ref>)–(<ref>), to verify
the correctness of the numerical implementation. 
Fixed units of time (Δ t=1.0), energy (k_B T = 1.0), and 
mass (m = 1.0) were used during simulations. Units of length were derived 
from a = Δ t √(k_BT/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 N̅^c = 20. 
Scattering geometry was defined by the wave vector k⃗ = (2/L,0,0).
In correlations obtained from MPCD-AT±a implementations, constant A 
in equations (<ref>)–(<ref>) was considered an 
adjustable parameter. Reduced functions defined by
G_ρ = C_ρ/A, 
and 
G_Vα
=C_Vα/A, 
are compared against their theoretical counterparts obtained from 
equations (<ref>)–(<ref>) and 
(<ref>)–(<ref>) in 
figures <ref>–<ref>. 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.






§ 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 <cit.>. 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 <cit.>. 
Modern LC technologies have promising applications in colloid manipulation, 
detection of biological agents, and medical diagnosing <cit.>. 
Numerical simulations could play an important role in helping developers to test 
models of these emerging technologies <cit.>.

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,α = u_i,α(t), which is a unit
vector. Furthermore, MPCD rules described in
section <ref> are augmented to cover three
main aspects, namely: multi-particle collision events cause changes in 
u_i,α; velocity gradients produce torques on u_i,α; and 
reorientation induces flow (backflow effect). Implementation of these effects is 
described below.




 §.§ 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

    Q^c_αβ      
    		       = 1/2 N^c∑_j ∈ c( 
                             3 u_j,α u_j,β - δ_αβ),

where N_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^c_αβ, called the scalar 
order parameter, S^c. The unit eigenvector associated to 
S^c, n^c_α, is called the director and 
represents the average molecular orientation. If nematogen i is found in a 
cell with the order parameter S^c and director 
n^c_α, it is assumed that it interacts with nematogens 
within the same cell according to the (mean-field) Maier-Saupe potential 
energy <cit.>,

    U_mf = -3/2 U S^c(
    		 u_i,αu_i,β
                     n^c_αn^c_β - 1
    		 ),

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,α^', from the canonical distribution
P = Z^-1exp[ -U_mf/(k_BT) ]sinθ, where 
Z^-1 is the normalisation constant and θ is the angle between 
u_i,α and n^c_α.




 §.§ Reorientation by flow




Flow effects on u_i,α 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 <cit.>

    u_i,α(t + Δ t) 
    = u_i,α^' 
    + χ_HI[
      Ω^c_αβu_i,β^'
    +λ(D^c_αβ u_i,β^' 
    -D^c_α^'β^'u_i, β^'^' 
                                          u_i,α^'^' 
                                          u_i,α^')
      ],

where Ω^c_αβ and 
D^c_αβ stand for the antisymmetric and symmetric
parts of ∂ v_α/∂ x_β, 
respectively.  Equation (<ref>) introduces quantities 
λ and χ_HI. The former is the tumbling parameter. 
It determines the tendency of u_i,α to be aligned by shear at a 
stationary angle (λ>1), or to continuously rotate in a direction 
consistent with the vorticity of the flow (λ < 1). During simulations, 
λ can be fixed to produce the desired dynamic effect <cit.>. 
Parameter χ_HI is an auxiliary one. It is initialised within the 
range [0,1] to control the overall flow-orientation coupling. By setting 
χ_HI = 0, such a coupling vanishes, while the maximum reorientation by 
flow is simulated for χ_HI=1.




 §.§ Backflow




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

    w^c_α
    =-1/Δ t∑_j ∈ cε_αα^'β^' 
    u_j,α^'(t+Δ t)
    u_j,β^'(t).

Under the assumption of over-damped rotational dynamics with viscosity 
γ_R, this implies a net torque, 
Γ^c_α = -γ_R w^c_α, and the 
production of an angular momentum in the cell 
Δℒ^c_α=Γ^c_αΔ t.

Backflow is taken into account by adding this angular momentum contribution to
Δ L^c_α in the MPCD-AT+a collision rule, 
equation (<ref>). 




§ FLUCTUATING NEMATODYNAMICS




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

    n_α(r⃗,t) = n^eq_α 
                                     + δ n_α(r⃗,t),

where n^eq_α 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 <ref>. 
Specifically, following general hydrodynamic models of NLCs <cit.>, the 
relaxation equation for the director field is proposed to have the form

    ( ∂/∂ t 
          + V_β∂/∂ x_β) 
    n_α = Y_α + Υ_α  ,

where Y_α is sometimes termed a quasi-current due to the fact 
that the director is not a conserved quantity and, consequently, the integral of 
Y_α over a surface cannot be a flux <cit.>. In 
equation (<ref>), Υ_α is a stochastic 
contribution to Y_α, which can be assumed to be a Gaussian-Markov 
process with zero mean, strength γ, and FDT

    ⟨Υ_α(r⃗^ ',t^')  
                 Υ_β(r⃗,t)  
    ⟩ = 2 k_B T γ(
                    δ_αβ - n_α n_β) 
                    δ(r⃗^ ' -r⃗)
                    δ(t^'-t),

where matrix δ_αβ - n_α n_β projects vectors 
on a direction perpendicular to n_α. This permits to satisfy the 
condition n_αδ n_α = 0, imposed by the normalisation
of n_α.

In MPCD-N, time evolution of orientations is promoted by shear, 
equation (<ref>). Furthermore, streaming and collision 
dynamics could be anticipated to give rise to orientation diffusion, then 
Y_α is written as

    Y_α    =  -D_n( 
                              δ_αβ - n_α n_β) 
                              ∂/∂ x_β^'∂/∂ x_β^' 
                              n_β 
                 + χ_HI[ 
    	       Ω_αβ n_β
                 + λ(
    		     D_αβ n_β
    		    -D_α^'β^' 
                         n_β^'n_α^'n_α)
    	     ],

where the orientation diffusion coefficient is D_n and 
Ω_αβ and D_αβ are the antisymmetric and symmetric
parts of ∂ V_α/∂ x_β. 

Dynamic equations for director fluctuations can be obtained by replacing 
equations (<ref>) and (<ref>) 
into equations (<ref>) and (<ref>), and
retaining only linear contributions. 
In a first approximation, mass and momentum transport can be assumed to be still 
given by equations (<ref>) and (<ref>), thus 
neglecting the effects of backflow <cit.>. In 
section <ref>, the validity of this assumption will be 
explored. It will be shown to be a very satisfactory approximation
for those values γ_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 ê_1 = k⃗/k, 
ê_2 = (n̂^eq×ê_1)/sinθ^eq, and 
ê_3 = ê_1×ê_2, where θ^eq is the angle 
between n̂^eq and k⃗.
For concreteness, the scattering geometry defined by k⃗ = (k,0,0), and 
n̂^eq = (0,0,1) will be considered hereafter. Then, the following system
of independent equations is obtained for fluctuating fields in MPCD-N,

    ( 
           
           -ω + D_n k^2    0 
    
           0     -ω + ν k^2) 
    ·(
           δñ_2
    δṼ_2) 
    =
    (
           Υ̃_2
     kρ_0^-1Σ̃_12) ,



    (
            -ω     - k ρ_0   
    
    	- k c_T^2ρ_0^-1    -ω  + D_l k^2) ·(δρ̃
    δṼ_1) =
             k/ρ_0(
          0 
    Σ̃_11) ,

and


    ( 
           
           -ω + D_n k^2    k χ_HI (λ -1)/2 
    
           0     -ω + ν k^2) 
    ·(
           δñ_1
    δṼ_3) 
    =
    (
           Υ̃_1
     kρ_0^-1Σ̃_13) .


It can be observed from 
equations (<ref>)–(<ref>) that in the present
nematodynamic model, density and velocity fluctuations are not perturbed
by orientation fluctuations though fluctuations δṼ_3 affect 
the dynamics of δñ_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 (<ref>) and (<ref>). Such a solution, 
together with the FDT in Fourier space,

    ⟨Υ̃_α( k⃗,ω) 
    	     Υ̃^*_β( k⃗,ω)⟩
    = 2 (2)^4γ k_B T (δ_αβ
    	                           -n_α n_β) 
                                  δ( k⃗^' +k⃗) 
    	                      δ(ω^'+ω),

can be used to obtain the following orientation ACFs

    C_n2(k⃗,ω)
           = ⟨δñ_2( k⃗,ω) 
             δñ_2^*( k⃗,ω) 
             ⟩
           = A γ/ω^2 + D_n^2k^4 ,


    C_n1(k⃗,ω)
         =  ⟨δñ_1( k⃗,ω) 
             δñ_1^*( k⃗,ω) 
             ⟩
         = A[
    	    γ/ω^2 + D_n^2k^4 
    	    +χ_HI^2(λ-1)^2/4 ρ_0ν k^4/(ω^2 + D_n^2k^4)
                                  (ω^2 +   ν^2k^4)] .


Consequently, C_n2 is not expected to be 
modified by velocity-orientation coupling, while the intensity of 
C_n1 is anticipated to increase with a term
proportional to χ_HI^2 (λ-1)^2.

Another correlation function that could bring information about the correct
coupling between velocity and orientation fluctuations is the cross correlation
function, C_n1,V3, given by

    C_n1,V3(k⃗,ω) =
    ⟨δñ_1(k⃗,ω)
             δṼ_3(k⃗,ω)
    ⟩ 
    =
    A χ_HI (λ - 1)/2 ρ_0ν k^3 (ω + D_n k^2)/(ω^2 + D_n^2k^4)
            (ω^2 +   ν^2k^4).





 §.§ Measurements




Measurements of correlations C_n1, C_n2, and 
C_n1,V3 can be carried out by introducing the time-series of 
spatial Fourier transform of orientation fluctuations

    δn̅_α(k⃗, t_l) 
    = ∑_j=1^N(
    		 u_j,α - n^eq_α)
    		 exp[-k⃗·r⃗_j(t_l) ].


Then, the same procedure used in section <ref> 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 (<ref>)–(<ref>), 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 n̂^eq 
towards the z <cit.>. This direction is 
considered as the ê_3 direction, while x and y are considered to 
coincide with ê_1 and ê_2, respectively.

Numerical experiments were conducted using the same simulation parameters 
Δ t, k_B T, m, N̅^c, a, and L, given 
in section <ref>. Nematic order was simulated using 
U=20k_BT. For this mean-field strength, it has been shown that the average 
scalar order parameter produced in collision cells, 
S̅^c=0.947, is in agreement with predictions of 
self-consistent models of NLCs <cit.>. Two
main features of hydrodynamic correlation functions in MPCD-N will be explored.
Namely, the effects of backflow controlled by the parameter 
γ_R; and the effects of velocity-orientation coupling 
controlled by parameters χ_HI and λ.




  §.§.§ Backflow effects




With the purpose of analysing the extent at which the backflow perturbs 
hydrodynamic fluctuations, γ_R was varied over the range 
[0.01,1.0], which extends by one order of magnitude the values 
considered before the present contribution. 
Eighteen experiments were conducted using combinations between parameters 
λ = 0.5, 2.0; χ_HI = 0.0, 0.5, 1.0, and 
γ_R = 0.01, 0.1, 1.0. Correlation functions 
C_ρ and C_Vα were measured and compared with those
obtained from simple MPCD-AT+a. No systematic differences were observed
between correlations C_ρ, C_V1 and C_V2 in 
MPCD-N and those measured in MPCD-AT+a, in good agreement with 
equations (<ref>) and (<ref>). A similar 
conclusion applies for function C_V3 for small and moderate values of the
product γ_Rχ_HI, namely, for 
γ_Rχ_HI < 0.5. Nevertheless, backflow effects
become noticeable, increasing the intensity of C_V3 for λ =0.5,
and decreasing it for λ = 2.0. This is illustrated in 
figure <ref>, where the results for C_V3 obtained for the whole
experimental setup are presented. 
Therefore, it can be concluded that if γ_R⩽ 0.1, as it has 
been used in references <cit.>, hydrodynamic fluctuations of flow and 
density in MPCD-N can be considered to be independent of orientation 
fluctuations.






  §.§.§ Velocity-orientation coupling




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

    G_n2(k⃗,ω)
    = C_n2(k⃗,ω)/A 
    = γ/ω^2 + D_n^2 k^4 ,


    G_n1(k⃗,ω) 
    = C_n1(k⃗,ω)/A 
    = γ/ω^2 + D_n^2 k^4 
    	+χ_HI^2(λ-1)^2/4 ρ^eqν k^4/(ω^2 + D_n^2k^4)
                                  (ω^2 +   ν^2k^4) ,

and

    G_n1,V3(k⃗,ω) 
    	  = C_n1,V3(k⃗,ω)/A 
    	  = χ_HI (λ - 1)/2 ρ^eqν k^3 ( -ω +  D_n k^2)/(ω^2 + D_n^2k^4)
                (ω^2 +   ν^2k^4)  .


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 γ and D_n 
could be inferred in terms of simulation parameters, N̅^c, 
Δ t, m, etc. Consequently, the subsequent discussion is conducted in a 
semiempirical form, using values for γ and D_n derived from 
the numerical experiments.
Correlations G_n2 and G_n1 showed the Lorentzian shape 
predicted from equations (<ref>) and 
(<ref>) with estimates, in s.u., 
γ_R=5.5× 10^-4 and D_n = 0.615. 
G_n2 was found to be independent of 
simulation parameters, as expected from equation (<ref>),
except for simulations conducted at γ_R = 1.0, λ = 2.0, 
and χ_HI= 1.0. In this case, G_n2 is more intense 
and exhibits two small lateral peaks located at ± q(k), indicating an 
apparent coupling with propagating modes. These features are illustrated in 
figure <ref>. 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 <cit.>. At large values of 
λ and χ_HI, this effect could be large enough to produce
a reduction of global order in the system. Specifically, for λ = 2.0 and
χ_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 <ref> where 
snapshots of systems simulated at 
{λ = 0.5,γ_R=1.0,χ_HI = 1.0} and 
{λ = 2.0,γ_R=1.0,χ_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 
n̂^eq, from which separation of longitudinal and transverse 
variables is constructed, becomes weaker. This implies that in the
case {λ = 2.0,γ_R=1.0,χ_HI = 1.0}, 
fluctuation δñ_y that is used to estimate G_n2,
could be coupled with the propagating velocity components, a situation not expected 
for the transverse variable δñ_2.





Correlation G_n1 was found to exhibit a
systematic dependence on the product (λ -1)^2χ_HI^2.
Results are in good agreement with equation (<ref>), as it
can be observed in figure <ref>. Analogous results were found in
the case of the cross correlation G_n1,V3,
whose real part obtained from equation (<ref>) is illustrated 
in figure <ref> 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 γ_R= 1.0, and 
χ_HI = 1.0, seem to deviate from this rule. In the case 
λ = 0.5, G_n1 and G_n1,V3 had a larger height 
than the one predicted by the linearized model. On the contrary, for 
λ = 2.0, both G_n1 and G_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_V3 already presented in figure <ref>. To be specific,
the increment (decrement) of C_V3 for λ=0.5 (λ=2.0)
causes a corresponding increment (decrement) of G_n1 and 
G_n1,V3 through the hydrodynamic coupling in 
equation (<ref>).








§ 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±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±a.

Afterwards, an extension of MPCD towards nematic phases has been presented. This 
implementation was based on the original one introduced in reference <cit.>. 
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 <cit.>. 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 <cit.>, 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, 
γ_R and χ_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 <cit.>, 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.


99


001 Malevanets A., Kapral R., 
              
              J. Chem. Phys., 1999, 110, 8605–8613, 
              10.1063/1.478857.
002 Malevanets A., Kapral R.,
              
              J. Chem. Phys., 2000, 112, 7260–7269, 
              10.1063/1.481289.
003 Frenkel D., Smit D., 
              Understanding Molecular Simulation: From Algorithms to Applications, 
              Academic Press, San Diego, 2002. 
003aIgnatyuk V.V., Mryglod I.M., Bryk T.,
              
              Condens. Matter Phys., 2018, 21, 13001,
              10.5488/CMP.21.13001.
006 Brady J.F., Bossis G., 
              
              Annu. Rev. Fluid Mech., 1988, 20, 111–157, 
              10.1146/annurev.fl.20.010188.000551.
007 Succi S., 
              The Lattice Boltzmann Equation for Fluid Dynamics and Beyond,
              Oxford University Press, Oxford, 2001.
008 Híjar H., Sutmann G., 
              
              Phys. Rev. E, 2011, 83, 046708,
              10.1103/PhysRevE.83.046708.
009 Híjar H., 
              
              Phys. Rev. E, 2015, 91, 022139,
              10.1103/PhysRevE.91.022139.
010 Matyka M., 
              
              Transp. Porous Media, 2017, 120, 633–641,
              10.1007/s11242-017-0944-7.
011 Chen J.X., Zhu J.X., Zhao Y.H., Sun W.G., Xu J.R., Ying H.P., 
              
              Commun. Nonlinear Sci. Numer. Simul., 2014, 19, 2505–2512,
              10.1016/j.cnsns.2013.11.013.
012 Di Cintio P., Livi R., Lepri S., Ciraolo G., 
              
              Phys. Rev. E, 2017, 95, 043203,
              10.1103/PhysRevE.95.043203.
013 Kowalik B., Winkler R.G.,
              
              J. Chem. Phys., 2013, 138, 104903,
              10.1063/1.4792196.
014 Eisenstecken T., Hornung R., Winkler R.G., Gompper G., 
              
              EPL, 2018, 121, 24003,

              10.1209/0295-5075/121/24003.
015 Lee K.W., Mazza M.G., 
              
              J. Chem. Phys., 2015, 142, 164110,
              10.1063/1.4919310.
016 Shendruk T.N., Yeomans J.M., 
              
              Soft Matter, 2015, 11, 5101–5110,
              10.1039/C5SM00839E.
021 Gompper G., Ihle T., Kroll D.M., Winkler R.G., 
              In: Advanced Computer Simulation Approaches for Soft Matter Sciences III. Advances in Polymer Science, Vol. 221,
              Holm C., Kremer K. (Eds.), Springer, Berlin, Heidelberg, 2009, 1–87. 
022 Ihle T., Kroll D.M., 
              
              Phys. Rev. E, 2003, 67, 066705,
              10.1103/PhysRevE.67.066705.
023 Bedeaux D., Mazur P., 
              
              Physica, 1974, 76, 247–258, 
              10.1016/0031-8914(74)90198-0. 
023aNoguchi H., Gompper G., 
              
              Phys. Rev. E, 2008, 78, 016706,
              10.1103/PhysRevE.78.016706.
023bTheers M., Winkler R.G., 
              
              Phys. Rev. E, 2015, 91, 033309,
              10.1103/PhysRevE.91.033309.
021aBryk T., Mryglod I., Scopigno T., Ruocco G., Gorelli F., Santoro M.,
	      
	      
              J. Chem. Phys., 2010, 133, 024502,
              10.1063/1.3442412.
030 De Gennes P.G., Prost J., 
              The Physics of Liquid Crystals,
	      Clarendon Press, Oxford, 1995. 
031 Guzmán O., Abbott N.L., de Pablo J.J.,
              
              J. Chem. Phys., 2005, 122, 184711,
	      10.1063/1.1896354.
031aKurochkin O., Atkuri H., Buchnev O., Glushchenko A., Grabar O., Karapinar R., Reshetnyak V., West J., Reznikov Yu.,
              
              Condens. Matter Phys., 2010, 13, 33701,
              10.5488/CMP.13.33701.
032 Ohzono T., Fukuda J., 
              
              Nat. Commun., 2012, 3, 701–708,
              10.1038/ncomms1709.
033 Zhou S., 
              Lyotropic Chromonic Liquid Crystals. From Viscoelastic Properties to Living Liquid Crystals,
	      Springer, Cham, 2017. 
034 Genkin M.M., Sokolov A., Lavrentovich O.D., Aranson I.S., 
              
              Phys. Rev. X, 2017, 7, 011029,

              10.1103/PhysRevX.7.011029.
035 Lee K.W., Pöschel T., 
              
              RSC Adv., 2017, 7, 42218–42224,
	      10.1039/C7RA06757G.
035aHíjar H., de Hoyos D.M., Santamaría-Holek I.,
              
              J. Chem. Phys., 2012, 136, 114109,
	      10.1063/1.3692600.
038 Doi M., Edwards S.F., 
              The Theory of Polymer Dynamics, 
	      Oxford University Press, New York, 1998.
038aPleiner H., Brand H.R.,
              
              In: Pattern Formation in Liquid Crystals,
	      Partially Ordered Systems, Buka A., Kramers L. (Eds.),
              Springer-Verlag, New York, 1996, 15–67.  
039 Híjar H., Halver R., Sutmann G., Fluctuation Noise Lett., 2019, 18, 1950011, 10.1142/S0219477519500111. 
040 Cleaver D.J., Allen M.P., 
              
              Phys. Rev. A, 1991, 43, 1918–1931,
	      10.1103/PhysRevA.43.1918.








 ,    ,  06140, , 




    —          .             .    ,               ,       .     ,     .                 . 
        . ,             .      . ,          ,          ,      .   ,        .


    ,  ,   ,     

