text:
# Possible 𝑍₂ phase and spin-charge separation in electron-doped cuprate superconductors

## I Introduction

A complete understanding of the behavior displayed by high temperature superconducting cuprates is still lacking. In these materials, d-wave superconducting (dSC) samples are obtained upon doping the parent compounds. The latter are Heisenberg antiferromagnets. The presence of low energy properties consistent with standard BCS theory in dSC samples is also generally undisputed. However, the phase-space region intervening between the two aforementioned regimes hosts an unusual phenomenology – and a quite debated one as well. In this paper, we will adopt the point of view that the unconventional behavior in underdoped cuprates is evidence of doped Mott insulator physics. To capture the important properties of these materials, which are then controlled by the large Mott gap and the two-dimensional spin interactions in the copper-oxide planes A8796 , we will use the slave-boson approach developed in Refs. SLAVE\_BOSON,  and WENLEENN, .

In underdoped electron-doped (ED) materials the antiferromagnetic (AF) phase is very robust. However that is not the case in their hole-doped (HD) counterparts. Instead, a pseudogap metallic normal state appears in underdoped samples. This is a paramagnetic state with a $`d_{x^2y^2}`$ gaped spectral function. Reconciling the gap with the strong local AF fluctuations is not trivial though. RW\_SPIN  One way to understand it is offered by the slave-boson approach to the $`tJ`$ model. In particular, in Ref. WENLEENN,  a translation symmetric pseudogap metallic state (the staggered-flux or s-flux state) was proposed. It contains spinons and holons interacting via a $`U(1)`$ gauge field. This state was shown to bear enhanced AF fluctuations compatible with experiments.RW\_SPIN  The $`U(1)`$ gauge fluctuations play a key role in the AF instability at very low doping as well as in the destruction of quasiparticle peaks.RW\_SPECTRAL ; FT0103  However, as discussed in Ref. RW\_SPIN, , the appearance of the spin-pseudogap below a certain energy and the existence of well defined electronlike quasiparticles in the nodal direction above $`T_c`$ suggest that spin and charge recombine before the s-flux state reaches the dSC state.

There is a competing scenario to describe this crossover from the s-flux to the dSC phase which involves the emergence of an intervening state with a $`Z_2`$ gauge structure.SF0050  This state would result from the condensation of spinon bilinears.RS9173 ; WENsrvb ; WEN  The dSC state would then follow from holon condensation. $`Z_2`$ spin liquids are very interesting states with true spin-charge separation. They were predicted theoretically over ten years ago.RS9173 ; WENsrvb  Yet they lack experimental realization. The present work proposes to study the conditions favoring the emergence of these spin liquid states as well as to discuss their relevance for the cuprate superconductors.

Despite the effort to check the applicability of $`Z_2`$ spin liquids in the context of high $`T_c`$ superconductivity, so far all the experimental tests gave negative results.Z2CHECK\_EX ; Z2CHECK\_TH  However, this does not imply that $`Z_2`$ states do not exist in high $`T_c`$ samples. In this paper we find that $`Z_2`$ states are indeed unlikely to exist in HD compounds. But a $`Z_2`$-linear state (i.e. with linear gapless spinon excitations) is likely to appear as the pseudogap metallic state for ED materials. In superconducting (SC) samples without AF order spin liquid pseudogap signatures may emerge once SC coherence is destroyed by thermal fluctuations or an external magnetic field. In the first case such a behavior is still to be reported. Experimental evidence in the magnetic field driven normal state, however, seems to be consistent with $`Z_2`$ pseudogap phenomenology. KLEEFISCH\_BISWAS\_ALFF ; HILL

The $`tJ`$ model has a particle-hole symmetry. In order to break it we introduce the next-nearest-neighbor (NNN) hopping $`t^{}`$ term. Taking both signs of $`t^{}`$ covers both ED ($`t^{}>0`$) and HD ($`t^{}<0`$) cases. In Sec. II we introduce the $`SU(2)`$ MF theory for the $`tt^{}J`$ model. In Sec. III we present the resulting MF phase diagrams. For $`t^{}<0`$ our calculations are consistent with the previous studies for $`t^{}=0`$ where only the s-flux state was obtained. However we find that NNN hopping can stabilize a $`Z_2`$-linear state in the $`tt^{}J`$ model for values of $`t^{}0.5J`$. These results show that unfrustrated hopping favors fractionalized phases (Sec. IV). In Sec. V we establish the link between our results and some properties of ED cuprates. We draw our main conclusions in Sec. VI.

## II  $`SU(2)`$ MF theory for the $`tt^{}J`$ model

The $`tt^{}J`$ model is given by the Hamiltonian

$$H_{tt^{}J}=H_{hf}+H_{hop}+H^{}$$
(1)
with

$`H_{hf}`$ $`=`$ $`{\displaystyle \underset{<ij>}{}}J_{ij}(\stackrel{}{S_i}.\stackrel{}{S_j}{\displaystyle \frac{1}{4}}),`$
$`H_{hop}`$ $`=`$ $`{\displaystyle \underset{<ij>}{}}t_{ij}𝑷(c_i^{}c_j+c_j^{}c_i)𝑷,`$
$`H^{}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{<ij>}{}}J_{ij}\left(1n_in_j\right),`$ (2)
where sums are taken over pairs of sites $`<ij>`$. The exchange coupling $`J_{ij}`$ is $`J`$ for $`<ij>`$ nearest-neighbor (NN) sites while the hopping parameter $`t_{ij}`$ equals $`t`$ for $`<ij>`$ NN sites, $`t^{}`$ for $`<ij>`$ NNN sites. $`\stackrel{}{S_i}`$ is the spin on site $`i`$, $`c_i^{}`$ and $`n_i=c_i^{}c_i`$ are electron creation and occupation number operators and $`𝑷`$ projects out doubly occupied sites.

To deal with the no double occupancy constraint the electron operator is decoupled according to $`c_i=\frac{1}{\sqrt{2}}(\psi _{i1}b_{i1}^{}+\psi _{i2}b_{i2}^{})`$ and $`c_i=\frac{1}{\sqrt{2}}(\psi _{i2}^{}b_{i1}^{}\psi _{i1}^{}b_{i2}^{})`$ and the spin operators are expressed as $`\stackrel{}{S_i}=\frac{1}{2}f_{i\alpha }^{}\stackrel{}{\tau }_{\alpha \beta }f_{i\beta }`$ where the $`SU(2)`$ doublets $`\mathrm{\Psi }_i^{}=[\psi _{i1}^{}\psi _{i2}^{}]=[f_{i1}^{}f_{i2}]`$ and $`h_i^{}=[b_{i1}^{}b_{i2}^{}]`$ represent spinons (fermionic particles carrying spin $`\frac{1}{2}`$) and holons (bosonic particles carrying unit charge) respectively and $`\stackrel{\mathbf{}}{𝝉}`$ are the Pauli matrices. These operators were introduced in Ref. WENLEENN,  to keep the $`SU(2)`$ gauge structure of $`H_{tt^{}J}`$ even away from half-filling. Parton operators introduce unphysical degrees of freedom which are projected out by implementing the constraint $`\stackrel{}{Q}_i\mathrm{\Psi }_i^{}\stackrel{}{\tau }\mathrm{\Psi }_i+h_i^{}\stackrel{}{\tau }h_i=0`$ (i.e. physical states are $`SU(2)`$ singlets on every site). Introducing a Lagrange multiplier per lattice site ($`\stackrel{}{a}_{0i}`$) the problem is reduced to one of lattice gauge theory.WENLEENN  Performing a Hartree-Fock type of decoupling of $`H_{hf}`$ and $`H_{hop}`$ (neglecting $`H^{}`$ for the time being) and implementing the projection constraint at MF level only, the following MF Hamiltonian is obtained:

$`H_{MF}`$ $`=`$ $`{\displaystyle \frac{3}{16}}{\displaystyle \underset{<ij>}{}}J_{ij}Tr\left(U_{ij}U_{ji}\right)+{\displaystyle \frac{1}{2}}{\displaystyle \underset{<ij>}{}}t_{ij}Tr\left(U_{ij}B_{ji}\right)`$ (3)
$`+`$ $`{\displaystyle \underset{i}{}}\stackrel{}{a}_{0i}.\left(\mathrm{\Psi }_i^{}\stackrel{}{\tau }\mathrm{\Psi }_i+h_i^{}\stackrel{}{\tau }h_i\right)\mu _b{\displaystyle \underset{i}{}}h_i^{}h_i`$
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{<ij>}{}}t_{ij}(h_i^{}U_{ij}h_j+h.c.)`$
$``$ $`{\displaystyle \underset{<ij>}{}}[\mathrm{\Psi }_i^{}({\displaystyle \frac{3}{8}}J_{ij}U_{ij}+{\displaystyle \frac{1}{2}}t_{ij}B_{ij})\mathrm{\Psi }_j+h.c.]`$
where $`\mu _b`$ is the chemical potential of holons.

Despite ignoring gauge fluctuations, we expect MF theory to account for some qualitatively right features. Indeed, we have the following:
i) The energetics of the model is correctly captured at MF level. This fact is suggested by the $`SU(2)`$ slave-boson MF theory phase diagram for the $`tJ`$ model – it includes the pseudogap regime (s-flux state) that turns into the dSC state as holons become coherent, as well as the strange-metal and Fermi-liquid regimes.WENLEENN  In particular, the pseudogap metallic state was predicted by the slave-boson approach prior to experimental observation;KL
ii) The nature of fluctuations in each phase is also consistent with experiments. For instance, the $`U(1)`$ gauge structure in the s-flux state accounts for the spin excitation spectrum in underdoped HD materials.RW\_SPIN  It also leads to the lack of well defined quasiparticles.RW\_SPECTRAL  In the strange-metal regime $`SU(2)`$ gauge fluctuations are responsible for the observed non-Fermi liquid behavior;
iii) MF states have properties that survive projection even in the presence of long range $`U(1)`$ gauge fields. Take the s-flux state, which includes both staggered current and d-wave pairing fluctuations, as an example. Remarkably, in Ref. IVANOVLW,  staggered vorticity correlations were reported to emerge from Gutzwiller projecting the dSC wave function. There is also numerical evidence that $`Z_2`$ MF states lead to fractionalized phases after performing Gutzwiller projection.IVANOVSENTHIL

## III  MF phase diagrams

To implement the self-consistent MF equations for the $`tt^{}J`$ model, the uniform ansatz considered in Ref. WENLEENN,  must be extended to include the NNN MF parameters $`U_{i,i+\widehat{x}\pm \widehat{y}}`$. Here we only consider translation invariant ansatsë which do not break any symmetry. In our MF calculation we use the following ansatz which properly comprises both the s-flux and dSC states:

$`U_{i,i+\widehat{x}}`$ $`=`$ $`\chi \tau ^3+\eta \tau ^1,`$
$`U_{i,i+\widehat{y}}`$ $`=`$ $`\chi \tau ^3\eta \tau ^1,`$
$`U_{i,i+\widehat{x}\pm \widehat{y}}`$ $`=`$ $`\gamma \tau ^3,`$
$`a_0^3`$ $``$ $`0.`$ (4)
Depending on the values of $`\chi `$, $`\eta `$, and $`\gamma `$, the above ansatz can describe phases with the same symmetry but different quantum orders.WEN  These phases may contain $`SU(2)`$, $`U(1)`$, or $`Z_2`$ gauge fluctuations. Some phases display a large Fermi surface while others have only Fermi points. These different MF states and their labels are summarized in Table 1.

The MF phase diagrams are determined for $`t=3J`$ and $`t^{}`$ between $`1.5J`$ and $`1.5J`$. These values are representative of the cuprates (note that $`t^{}`$ is close to $`J`$ for ED samples). Holons may condense at nonzero temperatures due to interactions or to the projection constraint (as proposed in Ref. WENLEENN, ). However, at MF level holons are free and conventional holon condensation only occurs at $`T=0`$. In our calculations we take $`T>0`$ and, hence, only consider the unconventional mechanism.

The resulting $`(t^{},x)`$ phase diagram for $`T=0^+`$ is shown in Fig.1(a).

For $`t^{}0.5J`$ the s-flux state survives all the way to $`T=0^+`$ – it becomes the dSC state for $`T=0`$. However, for $`t^{}0.5J`$ new MF states are obtained. In the underdoped regime, where the $`SU(2)`$ theory was proposed to be relevant, the Z2 phase emerges. When $`T=0`$ holons condense in the band bottom and the Z2 state changes into the dSC phase. For high enough doping we find the dSC state even at nonzero temperature – resulting from projection constraint driven condensation. As the pseudogap closes ($`\eta 0`$) a state with $`U(1)`$ gauge structure and gapless spinons arises (U1). Upon holon condensation it becomes a Fermi liquid (FL).

For $`t^{}0.5J`$ we obtain the $`(x,T)`$ phase diagram reported previously.WENLEENN  The MF finite temperature phase diagrams are illustrated for $`t^{}=0.8J`$ and $`t^{}=1.2J`$ in Figs. 2(a) and 2(b) respectively. As, by construction, the AF phase is absent, the $`Z_2`$-linear state appears as the dominant phase at low doping. ED samples can have three different normal metallic states: (a) a Z2 state with a $`d`$-wave pseudogap, true spin-charge separation and a nontrivial topological order;RS9173 ; WENsrvb ; SF0050  (b) a strange metal (the U1 phase) with a large Fermi surface and a $`U(1)`$ gauge interaction, and (c) a FL phase.

The above discussion of the MF phase diagrams is qualitatively correct even after including the dropped $`H^{}`$ term. In terms of spinons and holons, $`H^{}`$ has the form

$`H^{}`$ $`=`$ $`{\displaystyle \frac{J}{2}}{\displaystyle \underset{i}{}}(\mathrm{\Psi }_i^{}\stackrel{}{\tau }\mathrm{\Psi }_i)(h_i^{}\stackrel{}{\tau }h_i)`$ (5)
$`{\displaystyle \frac{J}{16}}{\displaystyle \underset{i}{}}{\displaystyle \underset{\widehat{u}=\widehat{x},\widehat{y}}{}}h_i^{}h_ih_{i+\widehat{u}}^{}h_{i+\widehat{u}}+O(x^3),`$
(which is $`x^2`$ at the MF level). While the second term correctly describes the attraction between holons resulting from the background with AF fluctuations, the first one should have the opposite sign to correctly account for the attraction between spinons and holons caused by gauge fluctuations. Hence, $`H^{}`$ should be considered together with an effective contribution arising from the gauge interaction. At MF level the constraint is implemented only globally, i.e. $`<_i\stackrel{}{Q}_i>=0`$. Including the following effective interaction, $`H^{\prime \prime }=\frac{3J}{16}_{i,\widehat{u}=\widehat{x},\widehat{y}}(\stackrel{}{Q}_i\stackrel{}{Q}_{i+\widehat{u}})^2`$, the fluctuations of $`(\stackrel{}{Q}_i)^2`$ away from zero are reduced. The prefactor in $`H^{\prime \prime }`$ is such that the exchange constant is renormalized to $`J_{eff}=1.5J`$ (note that Gutzwiller projection is known to increase the effective exchange constant ZHANG ). $`H^{\prime \prime }`$ also makes holons more massive. Including $`H^{}`$ and $`H^{\prime \prime }`$ at the MF level extends the region where the Z2 state emerges (Fig.1(b)). Hence gauge fluctuations stabilize the Z2 phase. As expected, projection has a smaller effect in this phase.

The above results support that the $`Z_2`$-linear spin liquid introduced in Sec. II is a competitive state for higher values of $`t^{}`$ and low doping levels. In particular, the lower bound $`t^{}0.5J`$ makes this state relevant for ED cuprates.

## IV  Spin-charge separation: the role of $`t^{}`$

The $`SU(2)`$ slave-boson approach is designed to address the low doping regime. As it is clear from Fig.1, for low doping there are two dominant phases, namely, the s-flux state and the Z2 state. The former becomes the latter as $`t^{}`$ is increased above a certain critical value. This parameter driven transition links two states with the same symmetries (they both describe totally symmetric spin liquids). However they have different quantum orders.WEN  This statement can be made more explicit by emphasizing the quite different nature of their low-energy effective behavior.

In the s-flux state massless Dirac fermions (spinons) and charged bosons (holons) interact with a long-range $`U(1)`$ gauge field. The gapless nature of these interacting excitations is protected by the quantum order.WEN  Still, the gauge field has a drastic effect on the nature of the excitations – this state is in fact a manifestation of a two-dimensional Luttinger liquid RW\_SPIN ; RW\_SPECTRAL ; FT0103  which has no well defined quasiparticles.

On the other hand, in the Z2 state the gauge interaction becomes short-ranged. Hence, the (linearly dispersing) spinons and holons are well defined quasiparticles. For that reason there is true spin-charge separation. The remnant discrete $`Z_2`$ gauge structure reflects the topological order of the state.WENsrvb  This topological order is related to the ground-state degeneracy of the system when it is embedded in a manifold with nontrivial topology. Such degenerate ground-states are locally similar but have different global (topological) properties.

The different gauge structure in the two states ($`U(1)`$ and $`Z_2`$) naturally gives rise to qualitatively distinct spectra for the collective modes. Indeed, the transition is accompanied by the opening of a gap in the fluctuations around the MF saddle-point.RIB\_WEN

Physically, the transition from the s-flux state to the Z2 state corresponds to the emergence of coherent diagonal (intrasublattice) charge carrier hopping – i.e. $`<c_i^{}c_{i+\widehat{u}}+c_{i+\widehat{u}}^{}c_i>`$, with $`\widehat{u}=\pm \widehat{x}\pm \widehat{y}`$, becomes nonzero. It results from a combination of two factors: (a) increasing $`t^{}`$ unfrustrates diagonal hopping and (b) electron/hole intrasublattice hopping is not frustrated and does not frustrate the background AF correlations (unlike intersublattice hopping). These play a crucial role in the dynamics of charge carriers (and of spin degrees of freedom as well).

In the s-flux state intrasublattice hopping is depleted by the negative $`t^{}`$. As charge carriers hop between different sublattices they interact strongly with the background AF fluctuations. Hence charge carrier motion and AF correlations frustrate each other. That is the motive underlying the staggered current and staggered chiral spin correlations.S\_CHIRAL  The presence of such correlations in the context of the $`tJ`$ model was reported in exact diagonalization studies.LEUNG  Moreover, charge carrier density and vorticity correlations are proportional to each other.IVANOVLW ; LEUNG  The emerging picture consists of charge carriers of opposite vorticity attracting each other in a background of staggered chiral spin fluctuations. This attraction eventually induces dSC coherence.

In the Z2 state charge carriers hop coherently along the diagonals and the staggered current correlations decrease. In particular, the energy difference between the dSC and staggered current states is seen to increase with $`t^{}`$. Therefore, staggered current fluctuations are reduced in the Z2 state. Exact diagonalization studies in doped $`tt^{}J`$ ladders show the same trend – i.e. staggered currents are stabilized by decreasing $`t^{}`$.TSUTSUI  According to our calculations, spinon pairing ($`\eta `$) is also reduced by the emergence of coherent diagonal hopping (therefore, charge carrier pairing may be expected to decrease as well). This is further consistent with the fact that the $`Z_2`$-vortex gap is physically connected to the reduction of staggered chiral spin fluctuations.RIB\_WEN ; SMALL\_CLUSTERS

Topological order in the Z2 state is associated with the intrasublattice hopping parameter $`\gamma `$. Studies of Gutzwiller-projected wave functions in Ref. IVANOVSENTHIL,  further corroborate the importance of the NNN hopping term (i.e. a non-vanishing $`\gamma `$) in realizing true fractionalization. Consequently, for underdoped samples, the $`Z_2`$-vortex gap in this fractionalized metallic pseudogap phase is set by $`\gamma `$ which is $`x`$. This is qualitatively different from the suggestion of Senthil and Fisher where the $`Z_2`$-vortex gap is set by the pseudogap scale ($`\eta `$).Z2CHECK\_TH

The arguments above support and delineate how the parameter $`t^{}`$ induces spin-charge separation. However, to fit experimental band structure, an extra hopping term ($`t^{\prime \prime }`$) must be included.KIM\_TSUTSUI  This requires an extension of the ansatsë considered in our analysis. Nevertheless, in a tight-binding model with $`t>0`$, $`t^{}>0`$ and $`t^{\prime \prime }<0`$ (as it is the case for ED cuprates), $`t^{\prime \prime }`$ is seen to frustrate hopping along the lattice bonds while unfrustrating hopping along the diagonals – in this sense enhancing the effect of $`t^{}`$.

## V Application to ED cuprates

As mentioned before, in ED cuprates the AF state is very stable. Such stability can be addressed within the $`tt^{}J`$ model framework. Indeed, an increasing $`t^{}`$ unfrustrates intrasublattice hopping and renders the charge carriers less effective in frustrating the background AF correlations.TOHYAMA\_GOODING  However, our MF results suggest that, if AF order can be destroyed by tuning parameters in the Hamiltonian, the resulting spin liquid may contain a pseudogap and correspond to the Z2 state.

A distinctive experimental feature of the Z2 phase is the presence of a $`d`$-wave pseudogap.AF\_GAP  Well, in the different Nd<sub>2-x</sub>Ce<sub>x</sub>CuO<sub>4-y</sub>, Pr<sub>2-x</sub>Ce<sub>x</sub>CuO<sub>4-y</sub> (PCCO), and La<sub>2-x</sub>Ce<sub>x</sub>CuO<sub>4-y</sub> (LCCO) SC thin films, the SC gap ($`510meV`$) is reported to survive in the magnetic field driven normal state.KLEEFISCH\_BISWAS\_ALFF  The value of the normal state gap is comparable to the SC gap and decreases with doping. Moreover, the low temperature normal state in PCCO SC crystals was reported to violate the Wiedemann-Franz law.HILL  This provides evidence for a non-Fermi liquid state compatible with spin-charge separation. According to our MF results, the aforementioned normal state may realize the Z2 state. On the other hand, if the Z2 state is to exist above $`T_c`$ it probably only covers a small fraction of the phase diagram due to the robust AF phase (Figs.2(c) and 2(d)).FLIQUID

Despite still disputed, there are a number of experiments reporting dSC order in ED samples.TSUEI  In our picture, the $`d_{x^2y^2}`$ gap in the SC phase is inherited from spinon pairing in the Z2 state. Recently, however, a SC pairing symmetry transition across optimal doping, from d-wave to extended s-wave, was reported in both LCCO and PCCO.PAIRING  We suggest the extended s-wave pairing may result in the FL state (Fig.2(d)) from exchange of collective modes – as the interaction between the Fermi sea and the massive gauge-fields raises an instability in the s-wave pairing channel.PATRICKLEE  We remark that this mechanism is unrelated to the one underlying dSC.

## VI Conclusions

In this paper we use a self-consistent $`SU(2)`$ slave-boson approach to study the role of NNN hopping in the context of the $`tt^{}J`$ model. We find that $`t^{}`$ induces a transition between the s-flux state (an algebraic spin liquid RW\_SPIN ; RW\_SPECTRAL ) and the Z2 state (a spin liquid with true spin-charge separation). This transition is interpreted as the emergence of coherent intrasublattice hopping in the Z2 state. The parameter $`t^{}`$ is also found to decrease staggered current and staggered chiral spin fluctuations, as well as to reduce pairing between charge carriers.

The range of parameters involved leads us to propose that a pseudogap metallic state in ED cuprates is likely to be described by a $`Z_2`$-linear phase. In contrast, the pseudogap state for HD samples is found to be described only by the s-flux phase with $`U(1)`$ gauge interaction.

Realizing fractionalized states in real systems is a major challenge. The present work suggests that the search for a pseudogap metallic state in ED cuprate samples may achieve such a goal.

We would like to acknowledge many discussions with W. Rantner. TCR was partially supported by the Praxis XXI Grant No. BD/19612/99 (Portugal). TCR and XGW were also supported by NSF Grant No. DMR–01–23156 and by NSF-MRSEC Grant No. DMR–02–13282.


format: text

text:
# How would you integrate the equations of motion in dissipative particle dynamics simulations?

## I Introduction

The mesoscopic phenomena of so-called “soft matter” physics Dao99 ; Cat00 ; deGennes96 , embracing a diverse range of systems including liquid crystals, colloids, and biomembranes, generally involve some form of coupling between different characteristic time- and length-scales. Computational modeling of such multi-scale effects requires new methodology applicable beyond the realm of traditional techniques such as ab initio and classical molecular dynamics Fre02 ; Ber98  (the methods of choice in the microscopic regime), and phase field modeling Eld01  or the lattice-Boltzmann method Rot97  (usually concerned with the macroscopic regime).

Dissipative particle dynamics (DPD) Hoo92 ; Esp95 ; War98 ; Gro97  is a particularly appealing technique in this regard. The “particles” of DPD correspond to coarse-grained entities, representing a collection of molecules or molecular groups rather than individual atoms. Coarse-graining leads to soft pair potentials allowing the particles to overlap (Forrest and Suter For95 ).

Although coarse-graining might also be considered implicit in Brownian and Langevin dynamics simulation, DPD offers the explicit advantage of a proper description of hydrodynamic modes significant in the physical approach towards a system’s equilibrium. This is achieved in DPD by implementing a thermostat in terms of pairwise random and dissipative forces such that the total momentum of the system is conserved. Due to these reasons, DPD has been used in studies covering a wide range of aspects in soft matter systems, including the structure of lipid bilayers Ven99 ; Gro01 , self-assembly Yam02 , and the formation of polymer-surfactant complexes Gro00 .

In practice, the pairwise coupling of particles through random and dissipative forces makes the integration of the equations of motion a nontrivial task. The main difficulty arises from the dissipative force, which depends explicitly on the relative velocities of the particles, while the velocities in turn depend on the dissipative forces. An accurate description of the dynamics requires a self-consistent solution.

The considerable computational load associated with this task has motivated the development of schemes Gro97 ; Nov98 ; Gib99 ; Pag98 ; Bes00 ; Vat02 ; Ott01 ; Sha01  which account for the velocity dependence of dissipative forces in some approximate manner, allowing the integration to be carried out to a sufficient degree of computational efficiency. The search for a satisfactory such integration scheme is ongoing, since many of the recent proposal have been found to exhibit non-physical behavior, such as systematic drift in temperature, and artificial structures in the radial distribution function Pag98 ; Bes00 ; Vat02 .

In order to overcome these problems, a number of new integration schemes for DPD simulations have been developed in the past few years. Self-consistent determination schemes exist on the one hand Pag98 ; Bes00 ; Vat02 , but these are rather elaborate.

Alternative proposals include (i) a parameterization of the integrator based on the specific application being modeled by den Otter and Clarke Ott01 , (ii) operator splitting by Shardlow Sha01 , and (iii) an elegant Monte Carlo-based approach due to Lowe Low99  which completely avoids the problems arising from random and dissipative forces as it does not use random or dissipative forces at all.

In this article, we apply these schemes respectively to specific model systems, with the objective of assessing their relative performance. To the authors’ knowledge, the latter three have yet to be tested and compared independently. The self-consistent approach has been tested, although not extensively, and the previously tested so-called DPD-VV (DPD version of the “velocity-Verlet” scheme) will be used here as a benchmark. By monitoring a number of physical observables including temperature, radial distribution function, radius of gyration for polymers, and tracer diffusion, we find that the methods by Lowe Low99  and Shardlow Sha01  give the best overall performance and are superior also to the integrators tested in previous studies Bes00 ; Vat02 . As will be discussed in the last section, a direct comparison between these two is not straightforward, since they are based on essentially different conceptual views of dissipative particle dynamics.

The paper is organized as follows. First, we briefly summarize the dissipative particle dynamics method and introduce the three model systems used here. In Sect. III we describe the integration algorithms and the update schemes in detail, and in Sect. IV we present the results from the tests. Finally, in Sect. V the findings and their relevance are discussed.

## II Methods and Models

Below we give a short summary of the dissipative particle dynamics method and describe the three model systems used in this work. For more thorough accounts on DPD, see e.g. Refs. War98 ; Gro97 .

### II.1 Dissipative Particle Dynamics

Dissipative particle dynamics describes a system in terms of $`N`$ particles having masses $`m_i`$, positions $`\stackrel{}{r}_i`$, and velocities $`\stackrel{}{v}_i`$. Interactions are composed of pairwise conservative, dissipative, and random forces exerted on particle $`i`$ by particle $`j`$, respectively, and are given by

$$\begin{array}{ccc}\stackrel{}{F}_{ij}^C\hfill & =& F_{ij}^{(c)}(r_{ij})\stackrel{}{e}_{ij},\hfill \\ \stackrel{}{F}_{ij}^D\hfill & =& \gamma \omega ^D(r_{ij})(\stackrel{}{v}_{ij}\stackrel{}{e}_{ij})\stackrel{}{e}_{ij},\hfill \\ \stackrel{}{F}_{ij}^R\hfill & =& \sigma \omega ^R(r_{ij})\xi _{ij}\stackrel{}{e}_{ij},\hfill \end{array}$$
(1)
where $`\stackrel{}{r}_{ij}\stackrel{}{r}_i\stackrel{}{r}_j`$, $`r_{ij}|\stackrel{}{r}_{ij}|`$, $`\stackrel{}{e}_{ij}\stackrel{}{r}_{ij}/r_{ij}`$, and $`\stackrel{}{v}_{ij}\stackrel{}{v}_i\stackrel{}{v}_j`$. The variables $`\gamma `$ and $`\sigma `$ are the strengths of the dissipative and random forces, respectively. The $`\xi _{ij}`$ are symmetric Gaussian random variables with zero mean and unit variance, and are independent for different pairs of particles and different times. The condition $`\xi _{ij}=\xi _{ji}`$ is employed to ensure momentum conservation.

The pairwise conservative force $`F_{ij}^{(c)}`$ is not specified by the DPD formulation and it can be chosen to include any forces that are appropriate for a given system, such as van der Waals and electrostatic interactions. In addition, it is important to notice that it is completely independent of the random and dissipative forces. Since one of the main motivations for using DPD is to be able to simulate systems at coarse-grained, or mesoscopic, scales, the conservative force is often chosen to be soft repulsive. Here, we use the “classical” DPD choice, i.e., soft repulsive conservative forces in the cases of model A and model B, and hard Lennard-Jones interactions combined with harmonic spring forces in the model polymer system.

In contrast to the conservative force, the random and dissipative forces are not independent, but are coupled through a fluctuation-dissipation relation. This coupling is due to the requirement that in thermodynamic equilibrium the system must have canonical distribution. The necessary conditions were first derived by Español and Warren in 1995 using a Fokker-Planck equation Esp95 .

The requirement of canonical distribution sets two conditions linking the random and dissipative forces in Eq. (1). The first one couples the weight functions through $`\omega ^D(r_{ij})=[\omega ^R(r_{ij})]^2`$, and the second one the strengths of the random and dissipative forces via $`\sigma ^2=2\gamma k_BT^{}`$. The latter condition fixes the temperature of the system $`T^{}`$ ($`k_B`$ being the Boltzmann constant) and relates it to the two DPD parameters $`\gamma `$ and $`\sigma `$.

Like classical molecular dynamics (MD) simulations, DPD allows studies of dynamical properties since the time evolution of particles can be described by the Newton’s equations of motion

$$\begin{array}{ccc}d\stackrel{}{r}_i\hfill & =& \stackrel{}{v}_idt,\hfill \\ d\stackrel{}{v}_i\hfill & =& \frac{1}{m_i}(\stackrel{}{F}_i^Cdt+\stackrel{}{F}_i^Ddt+\stackrel{}{F}_i^R\sqrt{dt}).\hfill \end{array}$$
(2)
Here $`\stackrel{}{F}_i^C=_{ij}\stackrel{}{F}_{ij}^C`$ is the total conservative force acting on particle $`i`$, and $`\stackrel{}{F}_i^D`$ and $`\stackrel{}{F}_i^R`$ are defined correspondingly. It is important to notice that the velocity increment due to the random force in Eq. (2) has the factor $`\sqrt{dt}`$ instead of $`dt`$. It can be justified by a Wiener process as in stochastic differential equations. Here, it suffices to notice that physically the Wiener process models intrinsic (thermal) noise in the system and provides the simplest approach to modeling Brownian motion using stochastic processes (see Refs. Esp95 ; Gro97  for a detailed discussion). The above continuous-time version of DPD satisfies detailed balance and describes the canonical NVT ensemble. In practice, however, the time increments in Eq. (2) are finite and the equations of motion must be solved by some integration procedure. We will return to this issue in Sect. III.

### II.2 An alternative approach to DPD

Dissipative particle dynamics described above can be thought of as a momentum conserving thermostat that allows one to study a system within the NVT ensemble with full hydrodynamics. The key features are therefore momentum and temperature conservation. As discussed above, momentum conservation arises naturally from pairwise forces. Temperature conservation, in turn, arises from the random and dissipative forces that are chosen to satisfy the fluctuation-dissipation theorem.

An alternative approach was formulated by Lowe in 1999 Low99 . It does not use dissipative or random forces at all, yet provides the same conservation laws and is similar in spirit to DPD as it is aimed for studies of coarse-grained models in terms of soft interactions. In Lowe’s method, one first integrates Newton’s equations of motion with a time step $`\mathrm{\Delta }t`$, and then thermalizes the system using the Andersen thermostat And80  for pairs of particles. We will discuss this method in detail in Sect. III.5.

Lowe’s approach is appealing for a number of reasons. First of all, since there are no dissipative forces we can assume that this method does not suffer from the same drawback as DPD: While DPD requires a self-consistent solution of the equations of motion, Lowe’s approach is easier to use and most likely performs well even with integration schemes that are commonly used in classical MD simulations. Secondly, the rate of how often the particle velocities are thermalized may be varied over a wide range, which implies that the dynamical properties of the system can be tuned in a controlled fashion. Lowe has demonstrated this idea by showing how some dimensionless variables (such as the Schmidt number $`Sc`$) can be tuned to match values found in actual fluids Low99 .

### II.3 Model systems

In this study, we evaluate the performance of a number of novel integration schemes that have been recently suggested for large-scale DPD simulations (see Sect. III). We test these integrators using three different model systems. The first two are based on a 3D model fluid system with a fixed number of identical particles. The first of them is aimed to clarify the performance of integration schemes in weakly interacting systems dominated by the random and dissipative forces, while the second model is relevant for systems in which the conservative interactions are of major importance. Finally, to gain insight into problems associated with hybrid models in which both soft and hard interactions are included, we consider a model of an individual polymer chain in a hydrodynamic solvent.

#### II.3.1 Model A

Model A describes the case characterized by the absence of conservative forces ($`F_{ij}^{(c)}=0`$). This choice corresponds to an ideal gas and it is customarily called “ideal DPD fluid”. The reason for using this model is that it provides us with some exact theoretical results that can be compared to results from model simulations. Here, the dynamics of the system arises only from thermal noise and dissipative coupling between pairs of particles. In DPD simulations, the random force strength is chosen to be $`\sigma =3`$ in units of $`k_BT^{}`$, and the strength of the dissipative force $`\gamma `$ is then determined by the fluctuation-dissipation relation $`\sigma ^2=2\gamma k_BT^{}`$.

The random and dissipative forces are chosen to be soft-repulsive,

$$\omega (r_{ij})=\{\begin{array}{cc}1r_{ij}/r_c\hfill & \text{for }r_{ij}<r_c;\hfill \\ 0\hfill & \text{for }r_{ij}>r_c,\hfill \end{array}$$
(3)
with a cut-off distance $`r_c`$ Gro97  and $`\omega ^R(r_{ij})=\omega (r_{ij})`$. This is the most common choice in DPD simulations and it has been used in recent investigations of integration schemes Bes00 ; Vat02 , thus allowing for a comparison of the present results with those of previous works. Although Eq. (3) has been used in virtually all published studies using DPD, it should be noted that the fluctuation-dissipation theorem does not specify the functional form of the weight function. The simple form of Eq. (3) simply provides a convenient choice.

In our simulations, a 3D simulation box of size $`10\times 10\times 10`$ with periodic boundary conditions is used. The length scale is defined by setting $`r_c=1`$, and a particle number density $`\rho =4`$. The energy scale is defined by setting the desired thermal energy to unity via $`k_BT^{}=1`$. All particles are identical, and thus $`m_i=m`$ for all $`i`$.

#### II.3.2 Model B

Model B is a simple interacting DPD fluid. Its main difference to model A is the presence of a conservative force, which we choose to be of the form $`F_{ij}^{(c)}(r_{ij})=𝒜\omega (r_{ij})`$. The amplitude of the force was chosen to be $`𝒜=25`$. This functional form for the conservative force is by far the most common choice in DPD simulations.

#### II.3.3 Model polymer system

The last model system considered in this work describes an individual polymer chain in an explicit hydrodynamic solvent. Our interest in a system of this kind originates from the fact that various soft matter systems such as liquid crystals and lipid bilayers are composed of particles which are essentially chain-like molecules. DPD serves well for studies of these systems due to the hydrodynamic nature of the solvent which plays an important role in various soft matter systems. However, while it is often desirable to describe chain-like molecules on a molecular level by hard interactions, complemented with bending and torsional potentials to account for the most relevant microscopic degrees of freedom, the solvent can often be described on a simpler level in terms of soft pair potentials. Thus, one possibility for efficiently modeling polymeric systems is a hybrid approach of chain-like molecules in a coarse-grained solvent.

The gain of using a hybrid approach in complex macromolecular systems is evident. It allows one to reduce the computational burden of dealing with an explicit solvent, while the molecular description of macromolecules is still accounted for in detail. However, the practical implications of including both hard and soft interactions in a model are not well understood. In a previous study for an ensemble of spherical particles described by hard conservative and soft dissipative forces, we found certain features which differentiated integration schemes from each other Vat02 . However, a full study of the performance of integration schemes within a true hybrid approach of a macromolecular system has been lacking up till now.

This model system aims to quantify the effects of integration schemes under conditions that combine both soft and hard interactions for a model polymer system. Here, the idea is to optimize the efficiency of the model by using a minimal approach. We thus describe the solvent as an ensemble of identical particles which interact via soft pairwise forces and satisfy momentum conservation, while the polymer chain is described on a more microscopic level in terms of (hard) Lennard-Jones interactions and harmonic springs.

The linear polymer chain is described as a chain of $`M`$ monomers connected by harmonic bonds whose potential follows the form

$$U_{\mathrm{harm}}=\frac{k}{2}|\stackrel{}{r}_i\stackrel{}{r}_{i1}|^2,i=2,3,\mathrm{},M,$$
(4)
with a spring constant $`k=7`$. Within the chain, the conservative monomer-monomer interactions are given by the truncated and shifted Lennard-Jones potential

$$U_{\mathrm{LJ}}(r_{ij})=\{\begin{array}{cc}4ϵ\left[\left(\frac{\mathrm{}}{r_{ij}}\right)^{12}\left(\frac{\mathrm{}}{r_{ij}}\right)^6+\frac{1}{4}\right]\hfill & ,r_{ij}r_c\hfill \\ 0\hfill & ,r_{ij}>r_c\hfill \end{array}$$
(5)
such that the potential is purely repulsive and decays smoothly to zero at $`r_c`$. We choose $`\mathrm{}=2^{1/6}`$ and $`ϵ=k_BT^{}`$, and therefore $`r_c=\mathrm{}\mathrm{\hspace{0.17em}2}^{1/6}=1`$. The pairwise conservative force acting on a monomer due to other monomers in a chain therefore follows directly from $`\stackrel{}{F}^C=(U_{\mathrm{harm}}+U_{\mathrm{LJ}})`$. The dissipative and random forces acting on the monomers are chosen to follow Eqs. (1) and (3) with $`\sigma =3`$.

The monomer-solvent and solvent-solvent interactions are described as in model A with $`\sigma =3`$ (and $`𝒜=0`$), i.e., the random and dissipative parts are used as a momentum conserving thermostat.

The justification for this choice of interactions lies in a wish to clarify the size of artifacts due to integration schemes in a case, where the polymer chain is described on a molecular level in terms of hard interactions, while the solvent is coarse grained as much as possible and is thus described by an ideal gas. The coupling between the solvent and the polymer chain comes from the dissipative forces which give rise to hydrodynamic modes. This eventually results in a minimal model of a polymer chain with full hydrodynamics under good solvent conditions. This was confirmed by studying the scaling behavior of the radius of gyration.

We consider polymers of size $`M=20`$ and use a 3D box of size $`10\times 10\times 10`$ with periodic boundary conditions, where the length scale is defined by $`r_c=1`$. The particle number density is chosen to be $`\rho =4`$ while the energy scale is defined by setting the desired thermal energy to unity via $`k_BT^{}=1`$. All particles are identical, and thus $`m_i=m`$ for both solvent and monomer particles.

#### II.3.4 Choice of random numbers

In the present work, uniformly distributed random numbers $`uU(0,1)`$ are used such that $`\xi _{ij}=\sqrt{3}(2u1)`$. This approach is highly efficient and yields results that are indistinguishable from those generated by Gaussian random numbers Gro97 . However, in the case of Lowe’s approach, the $`\xi _{ij}^{(g)}`$ used are true Gaussian random numbers.

## III Integrators

The integration schemes tested in this work have been chosen from the most recent ones that have been suggested in the literature but not tested and compared to other methods. They complement each other in the sense that the velocity dependence of the dissipative forces is accounted for in all cases, but the approaches differ substantially. We feel that all of the integrators considered here are promising candidates for large-scale simulations of soft matter systems. However, due to the lack of comparative studies in which all of these schemes would have been tested on equal footing, their relative performance has remained an open question. Here, we clarify this situation.

### III.1 Velocity-Verlet based integration scheme DPD-VV

In Table 1 we summarize the simplest integrator tested in this study. DPD–VV Bes00 ; Vat02  is based on the standard molecular dynamics velocity-Verlet algorithm Fre02 ; Ver67 ; All93  which is a time-reversible and symplectic second-order integration scheme. These properties can be proven by a straightforward application of the Trotter expansion Fre02 . The standard velocity-Verlet has been shown to be relatively accurate in typical MD simulations especially at large time steps All93 . The simplicity and good overall performance of the velocity-Verlet algorithm thus makes it a good starting point for further development.

We use the acronym DPD–VV for the modified velocity-Verlet. DPD–VV differs from the standard velocity-Verlet scheme in one important respect. As discussed above, the dissipative forces in DPD depend on the velocities, which in turn are governed by the dissipative forces \[see Eqs. (1) and (2)\]. This matter is not accounted for by the standard velocity-Verlet scheme. The DPD–VV, however, accounts for this complication in an approximate fashion by updating the dissipative forces \[step (5) in Table 1\] for a second time at the end of each integration step. This improves its performance considerably yet keeping it computationally efficient since the additional update of dissipative forces is not particularly time-consuming. In previous studies, the DPD–VV scheme has shown good overall performance Bes00 ; Vat02  for which reason we have chosen it as the “minimal standard” to which other integrators are compared.

### III.2 Self-consistent velocity-Verlet integrator

The update scheme of a self-consistent variant of DPD–VV is presented in Table 1. This SC–VV algorithm Bes00 ; Vat02  determines the velocities and dissipative forces self-consistently through functional iteration, and the convergence of the iteration process is monitored by the instantaneous temperature $`k_BT`$. This approach is similar in spirit to the self-consistent leap-frog scheme introduced recently by Pagonabarraga et al. Pag98 , which is the only other published self-consistent DPD integration scheme in addition to SC–VV (to the authors’ knowledge). The SC–VV scheme has the advantage of being very easy to implement as seen from Table 1. A recent study of the SC–VV scheme confirmed that it is a promising approach for interacting DPD systems Vat02 . That is particularly the case for the structural properties using long time steps in dense systems. As a drawback, the SC–VV has no advantage in temperature control as compared to other methods. For that, there exists a variant of the SC–VV integrator with a Nosé-Hoover type additional thermostat Bes00 ; Vat02 .

### III.3 Integrator by den Otter and Clarke

In 2001 den Otter and Clarke Ott01  proposed an approach which uses a leap-frog algorithm with predefined variables $`\alpha `$ and $`\beta `$. The idea in the OC-integrator, as it is called here, is to try to determine these parameters such that the effects due to the velocity dependence of dissipative forces are reduced as much as possible. The OC algorithm implementation\_of\_oc  is described in Table 2, in which the parameters $`\alpha `$ and $`\beta `$ describe the relative weight of the random forces with respect to the dissipative and conservative ones. They are determined prior to the actual DPD simulation by calculating the averages $`\stackrel{}{F}_i^D\stackrel{}{v}_i`$, $`\stackrel{}{F}_i^D\stackrel{}{F}_i^D`$, and $`\stackrel{}{F}_i^R\stackrel{}{F}_i^R`$ from an ensemble in which both the kinetic and configurational temperature equal the desired temperature $`k_BT^{}`$ Ott01 . Once they have been calculated, one obtains $`\alpha `$ and $`\beta `$ with a desired time step $`\mathrm{\Delta }t`$ from the equations

$$\alpha =\frac{1}{G\mathrm{\Delta }t}(1e^{G\mathrm{\Delta }t})$$
(6)
with $`G=\stackrel{}{F}_i^D\stackrel{}{v}_i/\stackrel{}{v}_i\stackrel{}{v}_i`$, and

$$\beta ^2=\frac{2m\alpha \stackrel{}{F}_i^D\stackrel{}{v}_i\alpha ^2\mathrm{\Delta }t\stackrel{}{F}_i^D\stackrel{}{F}_i^D}{\stackrel{}{F}_i^R\stackrel{}{F}_i^R}.$$
(7)
Note that both $`\alpha `$ and $`\beta `$ depend explicitly on $`\mathrm{\Delta }t`$. Since the averages $`\stackrel{}{F}_i^D\stackrel{}{v}_i`$, $`\stackrel{}{F}_i^D\stackrel{}{F}_i^D`$, and $`\stackrel{}{F}_i^R\stackrel{}{F}_i^R`$ can be derived analytically only for a limited number of systems, one usually has to calculate them from simulation (with a very small time step). This might be a problem in cases where properties such as density and temperature are varied over a wide range, since the parameters $`\alpha `$ and $`\beta `$ should (at least in principle) be calculated separately for all different conditions. Nevertheless, den Otter and Clarke have shown Ott01  that the OC algorithm performs well in both the ideal gas and a softly interacting DPD fluid.

For the three models considered in this work, we determined the parameters $`\stackrel{}{F}_i^D\stackrel{}{v}_i`$, $`\stackrel{}{F}_i^D\stackrel{}{F}_i^D`$, and $`\stackrel{}{F}_i^R\stackrel{}{F}_i^R`$ separately for all cases. Their values are shown in Table 3.

Note that the parameters for the model polymer system in Table 3 have been determined by averaging over all particles in the system. An alternate way would be to find $`(M+1)`$ different values for $`\alpha `$ and $`\beta `$ by averaging over the solvent and monomer particles separately. However, this would require a major computational study prior to actual simulations and is therefore not feasible. Besides, it might be against the usual spirit as integration schemes are typically based on prefactors that are identical for all particles in a system.

### III.4 Shardlow’s splitting method

The most recent addition to DPD integrators has been introduced by Shardlow Sha01 . He applied ideas commonly used in solving differential equations to the case of integrating the equations of motion in DPD. The key idea is to factorize the integration process such that the conservative forces are calculated separately from the dissipative and random terms. After this splitting the conservative part can be solved using traditional molecular dynamics methods, while the fluctuation-dissipation part is solved separately as a stochastic differential (Langevin) equation. To this end, Shardlow suggested two integrators, called S1 and S2, based on splitting the equations of motion up to first and second order, respectively.

The formal approach involves a first order splitting using the Trotter expansion Fre02 ; tuckerman92  (integrator S1) and a second order splitting using the Strang expansion strang  (integrator S2). The mathematical details of the derivation can be found in Sharlow’s original article Sha01 . It is important to notice that the power of the Trotter (Strang) expansion lies in the fact that it provides a general method for deriving symplectic algorithms. Importantly, the method works for both Hamiltonian and non-Hamiltonian systems; see Ref. tuckerman92  for a detailed discussion of the Trotter expansion and its applications.

Based on our simulation studies and the results presented in Ref. Sha01  for a system related to the model B in the present work, the performance of both of the two splitting methods was found to be excellent with S2 displaying slightly better overall characteristics. Since the first order method (S1) is more efficient, we have chosen it to be on the spotlight. Here we use the same naming convention and call it S1. The algorithm is presented in Table 4. To the best of our knowledge, further tests of S1 have not been reported yet.

The curious fact that the algorithm (Table 4) appears asymmetric for particles $`i`$ and $`j`$ is a result of the use of the fluctuation-dissipation theorem and Newton’s third law. The form in which the algorithm is presented keeps the notation otherwise symmetric.

### III.5 Lowe’s approach – Lowe-Andersen method

The approach introduced by Lowe in 1999 Low99  is presented in Table 5, which shows how one first integrates the Newton’s equations of motion with a time step $`\mathrm{\Delta }t`$, and then thermalizes the system as follows. For all pairs of particles for which $`r_{ij}<r_c`$, one decides with a probability $`\mathrm{\Gamma }\mathrm{\Delta }t`$ whether to take a new relative velocity from a Maxwell distribution. For each pair of particles whose velocities are to be thermalized, one works on the component of the velocity parallel to the line of centers and generates a relative velocity $`\stackrel{}{v}_{ij}^{\mathrm{\hspace{0.17em}0}}\stackrel{}{e}_{ij}`$ from a distribution $`\xi _{ij}^{(g)}\sqrt{2k_BT^{}/m}`$. Here $`\xi _{ij}^{(g)}`$ are Gaussian random numbers with zero mean and unit variance. This approach has its origin in the Andersen thermostat And80 , hence the name Lowe-Andersen method.

The key factor in Lowe’s method is the parameter $`1/\mathrm{\Gamma }`$ which describes the decay time for relative velocities. Since the condition $`0<\mathrm{\Gamma }\mathrm{\Delta }t1`$ is obvious, one finds that for $`\mathrm{\Gamma }\mathrm{\Delta }t=1`$ the particle velocities are thermalized at every time step, while for $`\mathrm{\Gamma }\mathrm{\Delta }t0`$ the model system is only weakly coupled to the thermostat. Thus the dynamical properties of the system can be tuned by the choice of $`\mathrm{\Gamma }`$ as shown by Lowe Low99 .

Although the present version of the algorithm follows the original one Low99  and is based on the velocity-Verlet scheme, it is clear that other approaches such as the leap-frog are equally useful, if desired. Further, based on the work by Lowe Low99 , this approach seems very promising although it has received only little attention by far Low99b .

In the present work for the three model systems considered here, we set $`\mathrm{\Gamma }`$ such that the tracer diffusion properties of the fluid are similar with those of DPD systems with chosen $`\sigma `$ in the limit $`\mathrm{\Delta }t0`$. In this fashion, we end up to a value of $`\mathrm{\Gamma }=0.745`$ for model A and to a value of $`\mathrm{\Gamma }=0.44`$ for model B. In the model polymer system we used the same value as in model A since the solvent is described in a similar fashion in both cases.

## IV Performance of integrators

### IV.1 Physical quantities studied

We characterize the integrators by studying a number of physical observables. After equilibrating the system, we first calculate the average kinetic temperature

$$k_BT=\frac{m}{3N3}\underset{i=1}{\overset{N}{}}\stackrel{}{v}_i^{\mathrm{\hspace{0.17em}2}},$$
(8)
whose conservation is one of the main conditions for reliable simulations in the canonical ensemble.

Next, in the cases of models A and B, we examine the radial distribution function $`g(r)`$ Boo80  which is one of the most central observables in studies of structural properties of liquids and solids. For the ideal gas (model A), the radial distribution function provides an excellent test for the integrators since then $`g(r)1`$ at the continuum limit. Therefore, any deviation from unity has to be interpreted as an artifact due to the integration scheme employed.

In model B, in which conservative interactions are present, there are no exact theoretical predictions for $`g(r)`$ that would allow a straightforward comparison of different integration schemes. A comparison is possible, though, in terms of physical observables such as the compressibility and the coordination number that are based on integrating $`g(r)`$. In the present study, we have chosen to consider the coordination number defined as

$$N_c=4\pi \rho _0^{r_1}\mathrm{dr}g(r)r^2,$$
(9)
where $`\rho `$ is the particle number density of the system and $`r_1`$ is the radial distance at which $`g(r)`$ has its first minimum after the leading (first) peak.

The radial distribution function reflects equilibrium (time-independent) properties of the system. To complement the comparison of different integrators, we also consider the tracer diffusion coefficient

$$D_T=\underset{t\mathrm{}}{lim}\frac{1}{6t}[\stackrel{}{r}_i(t)\stackrel{}{r}_i(0)]^2,$$
(10)
which can provide us with information of possible problems on the dynamics of the system. Here $`\stackrel{}{r}_i(t)`$ is the position of a tagged particle in models A and B, and the mean-square displacement is then averaged over all particles in a system to get better statistics for $`D_T`$. In the model polymer system, $`\stackrel{}{r}_i(t)`$ describes the center-of-mass position of the polymer chain via

$$\stackrel{}{r}_{\mathrm{cm}}(t)=\frac{1}{M}\underset{i=1}{\overset{M}{}}\stackrel{}{r}_i(t),$$
(11)
where the index runs over monomers in a chain.

For the model polymer chain, we further calculate the radius of gyration $`R_g\sqrt{R_g^2}`$ defined as

$$R_g^2\frac{1}{M}\underset{i=1}{\overset{M}{}}[\stackrel{}{r}_i\stackrel{}{r}_{\mathrm{cm}}]^2,$$
(12)
which shows that $`R_g`$ is a measure of polymer size. It is actually one of the most central quantities in polymer science and therefore serves as an excellent measure for our purposes.

In the following, the errors are stated in the figure captions and given as the magnitude of standard deviation.

### IV.2 Results for model A

As discussed above, model A is characterized by the absence of conservative forces, and thus any artifacts arising from the velocity-dependent forces are expected to be pronounced in this model. To study this possibility, we first discuss the deviations of the observed kinetic temperature $`k_BT`$ from the desired temperature $`k_BT^{}`$. The results for $`k_BT`$ shown in Fig. 1 indicate that DPD–VV and SC–VV are reasonably good at small time increments, but larger time steps lead to major deviations from the desired temperature. For OC, $`k_BT`$ decreases monotonically with $`\mathrm{\Delta }t`$ for $`\mathrm{\Delta }t0.15`$, after which the temperature increases rapidly. Nevertheless, the deviation is greater than in the case of DPD–VV but considerably smaller than for SC–VV. The Shardlow S1 integrator, however, has very good temperature control and the deviations remain less than 0.5 % up to $`\mathrm{\Delta }t=0.2`$. The best temperature control is found for the method by Lowe, however, yielding $`k_BT=k_BT^{}`$ for all time steps $`\mathrm{\Delta }t`$. In this case, we have extended the studies further and tested the behavior of $`k_BT`$ with various values of $`\mathrm{\Gamma }`$ between 0.1 and 10, but the conclusions remain the same.

Results for $`g(r)`$ are shown in Fig. 2. We find that the deviations from the ideal gas limit $`g(r)=1`$ are pronounced for OC, indicating that even for small time steps this integration scheme gives rise to unphysical correlations. The performance of DPD–VV is considerably better, although artificial structures are yet rather pronounced, while SC–VV and S1 lead to a radial distribution function that is close to the theoretically predicted one. Completely structureless $`g(r)`$ is found only for the integrator by Lowe, however. Again, in this case, we have tested the behavior of $`g(r)`$ with various values of $`\mathrm{\Gamma }`$, but the results remain the same. This confirms the expectation that $`\mathrm{\Gamma }`$ does not influence the equilibrium properties of the system.

The results for the diffusion coefficient $`D_T`$ in Fig. 3 are essentially consistent with the conclusions above. The integrator OC is not very useful in a system of the present kind, since it seems to lead to substantial deviations from the expected behavior. The SC–VV and the integrator by Lowe perform much better, while the integrators S1 and DPD–VV are most stable in this case.

### IV.3 Results for model B

Results shown in Fig. 4 for the observed kinetic temperature $`k_BT`$ reveal that the differences between the integration schemes deviations are weaker in model B than in model A. As it turns out below, this conclusion is generic and holds for all quantities studied here.

The deviations of $`k_BT`$ from the desired temperature $`k_BT^{}`$ are very minor for all integrators at small time steps $`\mathrm{\Delta }t0.01`$. Differences between the integrators become evident only at larger time steps. We first find how the temperature in SC–VV first decreases, then has a pronounced minimum around $`\mathrm{\Delta }t0.05`$, after which $`k_BT`$ increases very rapidly and the integrator eventually becomes unstable. The temperature conservation of OC is also relatively poor at time steps beyond $`\mathrm{\Delta }t=0.025`$. Best performance in this respect is found for the remaining integration schemes DPD–VV, S1, and Lowe, whose behavior is quite similar and comparable to each other.

Demonstrative results for the radial distribution functions in model B are shown in Fig. 5. It is clear that large time steps lead to major problems with regard to pair correlations. This is particularly clear at small distances ($`r<0.2`$). To quantify these changes, we calculated the coordination number \[see Eq.(9)\]. The results shown in Table 6 highlight the fact that all integration schemes converge to the same result at small time steps $`\mathrm{\Delta }t0.05`$, while for larger time steps there are significant deviations from the correct behavior found in the limit $`\mathrm{\Delta }t0`$. However, it is somewhat surprising that the coordination numbers obtained by different integration schemes are essentially similar within error limits, while based on $`g(r)`$’s there are noticeable differences between the pair correlation properties of different integrators. This is a consequence of the fact that due to the integration in Eq. (9), the deviations in $`g(r)`$ toward too small and too large values compensate each other. A similar effect has been observed recently in compressibility Vat02 , which is also defined as an integral over $`g(r)`$.

Tracer diffusion data shown in Fig. 6 is consistent with the conclusions above. For small time steps the results of all integration schemes are comparable, while for large time steps we can find how the differences become more and more pronounced. The scatter in the data does not allow us to make conclusive statements of the relative merits of the integrators, however.

### IV.4 Results for the model polymer system

Before we consider the results for the model polymer system, we would like to emphasize certain similarities it has with model A. Namely, in both cases the solvent is an ideal gas governed by dissipative and random forces only. Further, the model polymer system is dilute (more than 99.5% of the particles in a system are solvent particles) and there are no conservative interactions between the monomers and the solvent particles. This suggests that the artifacts due to integration schemes in the model polymer system would be essentially similar to those found in model A. It turns out below, however, that this is not the case.

Results shown in Fig. 7(a) indicate that the observed kinetic temperature $`k_BT`$ only rather weakly depends on the integration scheme. Results for S1, Lowe, and DPD–VV are all within one percent up to $`\mathrm{\Delta }t0.02`$ above which the integration schemes become unstable. The OC integrator is somewhat less reliable in this case, as it leads to a monotonous increase of $`k_BT`$, thus differing rather clearly from the results of other integration schemes.

A comparison of Figs. 1 and 7(a) provides one with an intriguing view of effects that arise from the hybrid approach. We first note that Lowe’s approach as well as S1 are equally good, in agreement with the conclusions made in model A. However, in model A the integration schemes were found to be stable up to very large time steps on the order of $`\mathrm{\Delta }t0.4`$, while in the model polymer system the largest time steps possible are about 0.02. This is due to the hard conservative monomer-monomer interactions used in describing the polymer chain, for which reason the size of the time step has to be reduced considerably as compared to model A. This suggests that, for practical purposes in hybrid models, one should seriously consider integration schemes with two different time steps, one for the solvent and another for the polymer degrees of freedom. Another interesting feature concerns the behavior of $`k_BT`$ in the case of OC. In model A, the observed kinetic temperature decreased monotonously down to timesteps $`\mathrm{\Delta }t0.15`$, while in the model polymer system the trend is the opposite. This is consistent with our results for model B, and suggests that the artifacts due to integration schemes depend significantly on the interactions chosen for the system.

To gain further insight into the performance of the integration schemes in a hybrid approach, we studied a number of physical quantities that specifically characterize the properties of the polymer chain. First, we studied the observed kinetic temperature of the polymer, defined as

$$k_BT_{\mathrm{chain}}=\frac{m}{3M}\underset{i=1}{\overset{M}{}}\stackrel{}{v}_i^{\mathrm{\hspace{0.17em}2}}.$$
(13)
This quantity characterizes the thermal fluctuations of the polymer chain, and therefore if there are any serious problems due to the integration schemes, then we expect that they are manifested in the behavior of $`k_BT_{\mathrm{chain}}`$.

We find that the results in Fig. 7(b) are wholly consistent with those presented in Fig. 7(a). Essentially, this implies that the temperature deviations in the whole system arise from the hard interparticle interactions used to describe the solute even though the system is dilute. We think that this finding is of generic nature and applies to both hybrid models and other DPD simulations in which all interactions are approximately of equal magnitude. In particular, it allows us to suggest that the artifacts due to integration schemes in DPD are predominated by the interactions that dictate the size of the time step.

In systems with hard conservative interactions the time step must be small. Otherwise, gradients in forces become too large and the system becomes unstable. Consequently, if the conservative force is stronger than the dissipative and random forces the artifacts due to integration schemes are driven by the conservative forces just like in classical molecular dynamics simulations. Naturally, the velocity-dependent dissipative forces are still playing a role but their effect is not as important as the influence of conservative interactions. From a practical point of view, this means that there is no particular reason to use an integrator which accounts for the velocity dependence of dissipative forces.

On the other hand, if all interactions are soft, or if the conservative forces are weak compared to the dissipative forces, then it is plausible that the velocity-dependence of dissipative forces is the underlying reason for artifacts due to the integration procedure. This is the case when the time step is determined by the dissipative force rather than the hard-core of the conservative potential. In this situation the quality of the integration scheme is very important, and the velocity dependence of the dissipative forces has to be accounted for by the integration scheme. It is clear that this matter warrants attention and should be accounted for in all subsequent studies by DPD.

Further studies of the radius of gyration and the diffusion coefficient of the polymer chain revealed the expected and undesired fact that computational studies of a dilute polymer system in an explicit solvent are very time consuming, and therefore the error bars remained rather large despite major computational efforts. Consequently, we found that the results for $`R_g`$ and $`D_T`$ of different integrators were essentially equal within error limits (data not shown). It is likely that more extensive calculations would have expressed deviations between different integration schemes, but we concluded that such studies were not worthwhile.

### IV.5 Computational efficiency

Besides the strength of the artifacts, we have paid attention to the computational efficiency of the integration schemes by calculating the cpu time needed for a single time step $`\mathrm{\Delta }t`$. Although this depends on various practical matters such as the computer architecture, the implementation of the algorithms, and the size of Verlet neighbor tables, we think this approach can provide one with the essential information of the relative speed of the integrators tested in this work.

The tests for efficiency have been carried out on a Compaq Alpha Station XP1000 with a 667 MHz processor. To this end, we used model B (with 4 000 particles at a density of $`\rho =4`$) with standard Verlet neighbor tables All93  and a time step of $`\mathrm{\Delta }t=0.05`$.

We calculated the cpu time needed for the integration of a single time step (averaged over 1 000 consecutive steps). The times needed to update the Verlet neighbor tables (or to calculate any physical quantities) are not included in these results. Thus the DPD–VV and SC–VV schemes were considered over steps (1) – (5) in Table 1, the OC approach over steps (1) – (3) in Table 2, Shardlow’s approach over steps (1) – (5) in Table 4, and Lowe’s integration scheme over steps (1) – (5) in Table 5. In the case of SC–VV, steps (4b) and (5) were repeated six times which guaranteed self-consistency in this case. Note that Lowe’s approach is based on normally distributed random numbers while other integration schemes use uniformly distributed ones. Since there are various methods available for the generation of normally distributed random numbers, this may affect the efficiency of Lowe’s approach to some extent randomnumbers .

The results shown in Table 7 indicate that Lowe’s method is substantially faster than OC and S1, which in turn are clearly faster than DPD–VV. Finally, the result that SC–VV is considerably slower than DPD–VV is not surprising due to the iteration process for the velocities and dissipative forces.

When these times are compared to each other, one should also keep in mind that DPD–VV, S1, and Lowe’s method are significantly easier to deal with compared to SC–VV and OC. In SC–VV, one needs to find the sufficient number of iterations prior to actual simulations to guarantee self-consistency. In OC, the preliminary work required prior to simulations is even more extensive, as one has to determine the parameters $`\alpha `$ and $`\beta `$ for a given system under desired thermodynamic conditions. This task may indeed take some time.

## V Discussion and Conclusions

In this work, we have tested several novel schemes on an equal footing through DPD simulations of three different model systems. The first of the models corresponds to a case where conservative interactions play no role, while the second model describes fluid-like systems with relatively strong but soft conservative potentials. Finally, the third model aims to characterize the quality of integration schemes in a hybrid approach for a dilute polymer system.

Of the integration schemes considered here, DPD–VV and SC–VV have recently been examined in Refs. Bes00 ; Vat02 . The results of the present study are consistent with previous findings: DPD–VV exhibits good overall performance, indicating that it presents a relatively accurate means to integrate the equations of motion at a reasonable computational cost.

Of the previously untested methods the Otter-Clarke (OC) method Ott01  is fast, performing especially well in interacting systems in which conservative forces are important. A drawback is that the parameters $`\alpha `$ and $`\beta `$ need to be determined prior to actual simulations through time-consuming precursory simulations with a very small time step. We note, however, that the properties of the OC scheme are in fact relatively insensitive to slight changes in $`\alpha `$ and $`\beta `$. Thus, for example, studies of the model polymer system using specifically determined $`\alpha `$ and $`\beta `$ values yielded results almost identical with studies based on the parameters of model A.

The Shardlow S1 integrator Sha01  is possibly the brightest star in this work. It performed very well in all models, and it is fast and rather easy to implement. We feel that it presents the best choice of integration schemes within the “usual” conceptual framework of DPD.

Interestingly, however, we have also found that the elegant and conceptually distinct method of Lowe Low99  performed excellently and is easy to implement. Furthermore, and what is important when Lowe’s method is compared to Shardlow’s integration scheme, it provides an alternative and a very attractive description of dissipative particle dynamics. Thus, a direct comparison of S1 and Lowe’s method is not meaningful. Instead, we discuss the pros and cons of these two approaches.

The usual DPD description is based on the idea that soft matter systems can be described in terms of softly interacting particles with some of the degrees of freedom coarse grained out and replaced with random noise coupled to dissipation. Temperature conservation is achieved through the fluctuation-dissipation theorem and the correct hydrodynamic behavior is guaranteed by momentum conservation Esp95 . Various studies have extended these ideas further. For example, Flekkøy et al. developed a DPD framework starting from a microscopic description Fle99 ; Fle00 . Español and coworkers, in turn, studied the dependence of transport properties of DPD fluids on the length and time scales Rip01  and a generalization of DPD to energy conserving systems Esp97 . DPD has recently been used together with molecular dynamics to coarse grain aqueous salt solutions Lyu02  in which the effective interactions used in DPD simulations were obtained from MD simulations by the inverse Monte Carlo procedure Lyu95 . The last ten years have been very successful on both the analytical and the computational fronts–the theoretical basis of DPD is now well established, and the number of applications has increased at a steady pace.

Lowe’s Low99  approach is a very recent inception and thus far has received limited attention. Although the theoretical foundations of Lowe’s method have yet to be fully worked out, it offers promising aspects that are not obvious in the traditional DPD description. To clarify these aspects, let us first remind ourselves that Lowe’s method does not include dissipation in the usual sense. Rather, it is based on a thermostat that thermalizes the velocities of pairs of particles at a rate which depends on the dynamical parameter $`\mathrm{\Gamma }`$. This parameter tunes the dynamical properties of the system. Lowe pointed out that the soft interactions used in DPD lead to a situation where the ratio of the kinematic viscosity and the diffusion coefficient of solvent particles (known as the Schmidt number Sc) is of the order of one. This value corresponds to a situation often found in gases, while in fluids $`\mathrm{𝑆𝑐}10^3`$ or even larger. To get closer to more realistic values for Sc, one can reduce the diffusion rate by using harder interparticle interactions, but this is against the philosophy of DPD and would reduce some of the benefits of the DPD approach.

Lowe’s approach is very different in this respect. It allows one to adjust the viscosity of the system to a desired value by varying the dynamical parameter $`\mathrm{\Gamma }`$ while the diffusive properties are not considerably affected since the conservative interparticle interactions remain soft. As a result, the Schmidt number can obtain values as large as $`10^7`$ Low99 . When compared to the usual DPD description, this implies that Lowe’s approach may be more feasible for describing hydrodynamic systems in which one needs to worry about the time scales of momentum diffusion and mass transfer with respect to the size of the colloidal particle.

There still remains the issue of the practical viability of Lowe’s approach, since we are not aware of any applications where the method by Lowe has been used. However, we are positive that this approach is a promising technique. For example, we have recently applied Lowe’s method to microphase separation of block copolymers in the spirit of Groot and Madden Gro98 , and it turned out that Lowe’s method was able to reproduce their results. Finally, in Figs. 8 and 9 we show how the velocity autocorrelation function depends on the choice of $`\mathrm{\Gamma }`$ for models A and B (Fig. 8), and how it is affected when $`\mathrm{\Gamma }`$ is varied but keeping the product $`\mathrm{\Gamma }\mathrm{\Delta }t`$ constant in the case of model B (Fig. 9). As discussed above, it is clear that large values of $`\mathrm{\Gamma }`$ lead to faster decay. However, the qualitative behavior of the velocity autocorrelation function is not seemingly affected, as illustrated by Fig. 9, and the effect of $`\mathrm{\Gamma }`$ on the diffusion coefficient was found not to be important for the studied combinations of $`\mathrm{\Gamma }`$ and $`\mathrm{\Delta }t`$.

To conclude, we have studied the performance of various novel integration schemes that have been designed specifically for DPD simulations. We have tested these integration schemes in three different model systems by varying the nature of interactions and found that the artifacts due to the integration scheme are essentially driven by the interactions that dictate the size of the time step. Thus, the artifacts and the performance of integrators are model dependent. Overall, we have found that there are two approaches whose performance is above the others. Of these, Shardlow’s integration scheme is based on splitting the equations of motion and can be applied to the usual DPD picture, while the approach by Lowe is distinctly different in nature and is related to the classical work by Andersen.

###### Acknowledgements.

We would like to thank Gerhard Besold, Wouter den Otter, Alex Bunker, and James Polson for useful discussions at the early stages of this work, Tony Shardlow for sharing his results prior to publication, and Nick Braun for a critical reading of the manuscript. This work has, in part, been supported by the Academy of Finland (I.V.) and by the Academy of Finland Grant No. 54113 (M.K.).


text:
# Replicator dynamical systems and their gradient and Hamiltonian properties

## 1 Introduction

Replicator equations were introduced by Fisher to capture Darwin’s notion of the survival of the fittest  and replicator dynamics is one of the most important dynamic models arising in biology and ecology , evolutionary game theory  and economics , traffic simulation systems and distributed computing  etc. It is derived from the strategies that fare better than average and thriving is on the cost of others at the expense of others (see e.g. ). This leads to the fundamental problem of how complex multi-agent systems widely met in nature can adapt to changes in the environment when there is no centralized control in the system. For complex “living” system such problem has been considered in . By term complex multi-agent systems under consideration we mean one that comes into being, provides for itself, and develops pursuing its own goals . In the same time replicator dynamics arises if the agents have to deal with conflicting goals and the behavior of such systems is quite different from the adaptation problem considered in . In this article we develop a new dynamic model to describe the multi-agent interaction in complex systems of interacting agents sharing common but limited resources.

Based on this model we consider a population composed of $`nZ_+`$ distinct competing “varieties” with associated fitnesses $`f_i(v)`$, $`i=\overline{1,n},`$ where $`v[0,1]`$ is the vector of relative frequencies of the varieties $`(v_1,v_2,\mathrm{},v_n)`$. The evolution of relative frequencies is described by the following equations:

$$dv_i/dt=v_i(f_i(v)f(v)),$$
(1)
where $`i=\overline{1,n},`$

$$f(v)=\underset{i=1}{\overset{n}{}}v_if_i(v).$$
(2)
The essence of (1) and (2) is simple: varieties with an above-average fitness will expand, those with a below-average fitness will contract.

Since $`v_i[0,1],i=\overline{1,n}`$ have to be nonnegative for all time, the system (1), (2) is defined on the nonnegative orthand

$$_+^n=\{v^n:v_i0\}.$$
(3)
The replicator equation describes the relative share dynamics and thus holds on the unit a n-1 dimensional space simplex

$$S_n=\{v_+^n:\underset{i=1}{\overset{n}{}}v_i=1\}.$$
(4)
Let us write the Fisher’s model in the following form

$$dv_i/dt=v_i\left(\underset{j=1}{\overset{𝑛}{}}a_{ij}v_j\underset{j,k=1}{\overset{𝑛}{}}a_{jk}v_jv_k\right),$$
(5)
where $`a_{ij},`$ $`v_i[0,1],`$ $`i,j=\overline{1,n},`$ $`n_+`$.

One can check the system (5) can be written in such a matrix commutative form:

$$dP/dt=[[D;P],P].$$
(6)
Here by definition,

$$P=\{\left(v_iv_j\right)^{\frac{1}{2}}:i,j=\overline{1,n}\},$$
(7)
$$D=\frac{1}{2}diag\{\underset{k=1}{\overset{𝑛}{}}a_{jk}v_k:j=\overline{1,n}\}.$$
(8)
It can be checked easily that the matrix $`P`$End$`E^n`$ is a projector in $`E^n`$, that is $`P^2=P`$ for all $`t`$, what appears to be very important for our further studying the structure of vector field (6) on the corresponding projector matrix manifold $`𝒫`$ . In particular, the expression (6) ensures that there exists the time invariant vector subspace Im$`PE^n`$ in the Euclidian phase vector space $`E^n,`$ necessary for the replicating system under regard to be the information data conserved.

## 2 Vector field analysis

In order to study the structure of the flow (6) on the projector matrix manifold $`𝒫P`$ let us consider a functional $`\mathrm{\Psi }:𝒫,`$ where by definition the usual variation

$$\delta \mathrm{\Psi }(P):=\text{Sp}(D(P)\delta P),$$
(9)
with $`D\mathrm{E}ndE^n,`$ $`\mathrm{S}p:\mathrm{E}ndE^n^1`$ being the standard matrix trace. Taking into account the natural metrics on $`𝒫`$, we consider the projection of the usual gradient vector field $`\mathrm{\Psi }`$ to the tangent space $`T(𝒫)`$ under the following conditions:

$$\phi (X;P):=\mathrm{S}p(P^2P,X)=0,\mathrm{S}p(\phi ,_\phi \mathrm{\Psi })_𝒫=0,$$
(10)
holding on $`𝒫`$ for all $`X\mathrm{E}ndE^n.`$The first condition is evidently equivalent to $`P^2P=0,`$ that is $`P𝒫`$ . Thereby we can formulate such a lemma.

Lemma 1. The functional gradient $`_\phi \mathrm{\Psi }(P),`$ $`P𝒫`$, at condition (10) has the following commutator representation:

$$_\phi \mathrm{\Psi }(P)=[[D,P],P].$$
(11)
Proof. Consider the projection of the usual gradient $`\mathrm{\Psi }(P)`$ to the tangent space $`T(𝒫)`$ of the manifold $`𝒫`$having assumed that $`P\mathrm{E}ndE^n`$:

$$_\phi \mathrm{\Psi }(P)=\mathrm{\Psi }(P)_\phi (\mathrm{\Lambda },P),$$
(12)
where $`\mathrm{\Lambda }\mathrm{E}ndE^n`$ is some unknown matrix. Taking into account the conditions (10), we find

$`_\phi \mathrm{\Psi }(P)`$ $`=`$ $`D\mathrm{\Lambda }P(D\mathrm{\Lambda })(D\mathrm{\Lambda })P+PD+DP`$ (13)
$`=`$ $`PD+DP+2P\mathrm{\Lambda }P.`$ (14)
where we made use the conditions

$$_\phi \mathrm{\Psi }(P)=D\mathrm{\Lambda }+P\mathrm{\Lambda }+\mathrm{\Lambda }P)$$

and

$$P(D\mathrm{\Lambda })+(D\mathrm{\Lambda })P+2P\mathrm{\Lambda }P=D\mathrm{\Lambda }.$$

Now one can see from (14) and the second condition in (12), that

$$P\mathrm{\Lambda }P=PDP$$
(15)
for all $`P𝒫`$, giving rise to the final result

$$_\phi \mathrm{\Psi }(P)=PD+DP2PDP,$$
(16)
coinciding exactly with commutator (11).

It should be noted that the manifold $`𝒫`$ is also a symplectic Grassmann manifold (), whose canonical symplectic structure is given by the expression:

$$\omega ^{(2)}(P):=\mathrm{S}p(PdPdPP),$$
(17)
where $`d\omega ^{(2)}(P)=0`$ for all $`P𝒫`$, and the differential form (17) is non-degenerate  upon the tangent space $`T(𝒫).`$

Let us assume that $`\xi :𝒫`$ is an arbitrary smooth function on $`𝒫`$ . Than Hamiltonian vector field $`X_\xi :𝒫T(𝒫)`$ on $`𝒫`$ generated by this function relative to the symplectic structure (17) is given as follows:

$$X_\xi =[[D_\xi ,P],P],$$
(18)
where $`D_\xi \mathrm{E}ndE^n`$ is a certain matrix. The vector field $`X_\xi :𝒫T(𝒫)`$ generates on $`𝒫`$ the flow

$$dP/dt=X_\xi (P)$$
(19)
being defined globally for all $`t`$. This flow by construction is evidently compatible with the condition $`P^2=P`$. This means in particular that

$$X_\xi +PX_\xi +X_\xi P=0.$$
(20)
Thus, we stated that dynamical system (6)being considered on the Grassmann manifold $`𝒫`$is Hamiltonian what makes it possible to formulate the following statement.

Statement 1. A gradient vector field of the form (18) on Grassmann manifold $`𝒫`$ is Hamiltonian with respect to the canonical symplectic structure (17) and certain Hamiltonian function $`\xi :𝒫`$, satisfying

$$\xi (P)=[D_\xi ,P],D_\xi =D,$$

where simultaneously $`_\phi \mathrm{\Psi }(P)=X_\xi (P)`$ for all $`P𝒫`$.

Consider now the (n-1)-dimensional Riemannian space

$$M_g^{n1}=\{v_i_+:i=\overline{1,n},\text{ }\underset{i=1}{\overset{𝑛}{}}v_i=1\}$$

with the metrics

$$ds^2(v):=d^2\mathrm{\Psi }_𝒫(v)=\underset{i,j=1}{\overset{𝑛}{}}g_{ij}(v)dv_idv_j_𝒫,$$

where

$$g_{ij}(v)=\frac{^2\mathrm{\Psi }(v)}{v_iv_j},i,j=\overline{1,n},\underset{i=1}{\overset{𝑛}{}}v_i=1.$$

Subject to the metrics on $`M_g^{n1}`$ we can calculate the gradient $`_\phi \mathrm{\Psi }`$ of the function $`\mathrm{\Psi }:𝒫`$ and set on $`M_\phi ^{n1}`$ the gradient vector field

$$dv/dt=_\psi \mathrm{\Psi }(v),$$
(21)
where $`vM_\psi ^{n1}`$, or $`\underset{i=1}{\overset{𝑛}{}}v_i=1`$ is satisfied. Having calculated (21), we can formulate the following statement.

Statement 2. The gradient vector fields $`_\phi \mathrm{\Psi }`$ on $`𝒫`$ and $`_g\mathrm{\Psi }`$ on $`M_g^{n1}`$ are equivalent or in another words vector fields

$$dv/dt=_g\mathrm{\Psi }(v)$$
(22)
and

$$dP(v)/dt=[[D(v),P(v)],P(v)]$$
(23)
generates the same flow on $`M_g^{n1}`$.

As a result from the Hamiltonian property of the vector field $`_\phi \mathrm{\Psi }`$ on the Grassmann manifold $`𝒫`$ we get a new statement.

Statement 3. The gradient vector field $`_g\mathrm{\Psi }`$ (21) on the metric space $`M_g^{n1}`$ is Hamiltonian subject to the non-degenerate symplectic structure

$$\omega _g^{(2)}(v):=\omega ^{(2)}(P)_{M_g^{n1}}$$
(24)
for all $`vM_g^{2m}`$ with the Hamiltonian function $`\xi _\psi :M_g^{2m}`$, where $`\xi _\psi :=\xi _{M_g^{n1}},`$ $`\xi :𝒫`$ is the Hamiltonian function of the vector field $`X_\xi `$ (18) on $`𝒫`$ . Otherwise if $`nZ_+`$ is arbitrary our two flows (22) and (23) are on $`𝒫`$ only Poissonian.

## 3 Spectral properties

Consider the eigenvalue problem for a matrix $`P𝒫`$, depending on evolution parameter $`t`$:

$$P(t)f=\lambda f,$$
(25)
where $`f^n`$ is an eigenfunction, $`\lambda `$ is a real eigenvalue $`P^{}=P`$, i.e. matrix $`P𝒫`$ is symmetric. It is seen from expression (22) that $`\mathrm{s}pecP(t)=\{0,1\}`$for all $`t`$. Moreover, taking into account the invariance of $`\mathrm{S}pP=1`$ we can conclude that only one eigenvalue of the matrix $`P(t),`$ $`t`$, is equal to 1, all others being equal to zero. So, we can formulate the next lemma.

Lemma 2. The image $`\mathrm{I}mPE^n`$ of the matrix $`P(t)𝒫`$ for all $`t`$ is k-dimensional $`k=rankP,`$and the kernel $`\mathrm{ker}PE^n`$ is $`(nk)`$-dimensional, where $`k_+`$ is constant, not depending on $`t`$ . As a consequence of the lemma we establish that at $`k=1`$ there exists a unique vector $`f_0E^n/(\mathrm{ker}P)`$ for which

$$Pf_0=f_0,f_0E^n/(\mathrm{ker}P).$$
(26)
Due to the statement above for projector $`P:E^nE^n`$ we can write down the following invariant in time expansion in the direct sum of mutually orthogonal subspaces:

$`E^n=\mathrm{ker}P\mathrm{I}mP.`$

Take now $`f_0E^n`$ satisfying the condition (26). Than in accordance with (22) the next lemma holds.

Lemma 3. The vector $`f_0E^n`$ satisfies the following evolution equation:

$$df_0/dt=[D(v),P(v)]f_0+C_0(t)f_0,$$
(27)
where $`C_0:`$ is a certain function depending on the choice of the vector $`f_0\mathrm{I}mP`$. At some value of the vector $`f_0\mathrm{I}mP`$ we can evidently ensure the condition $`C_00`$ for all $`t^n`$. Moreover one easily observes that for the matrix $`P(t)𝒫`$ one has  the representation $`P(t)=f_0f_0`$, $`f_0,f_0=1`$ , giving rise to the system (5) if $`f_0:=\{\pm \sqrt{v_j}_+:j=\overline{1,n}\}E^n`$.

## 4 Lyapunov Function

Let us consider gradient vector fields $`_\phi \mathrm{\Psi }`$ on $`𝒫`$ and $`_g\mathrm{\Psi }`$ on $`M_g^{n1}`$. It is easy to state that the function $`\mathrm{\Psi }:𝒫`$ given by (9) and equal on $`M_g^{n1}`$following expression

$$\mathrm{\Psi }(v)=\underset{i,j=1}{\overset{𝑛}{\frac{1}{4}}}a_{ij}v_iv_jc\underset{i=1}{\overset{𝑛}{}}v_i$$
(28)
under the condition $`c_+,`$ $`\underset{i=1}{\overset{𝑛}{}}v_i=1`$ being at the same time a Lyapunov function for the vector fields $`_\phi \mathrm{\Psi }`$ on $`𝒫`$ and $`_g\mathrm{\Psi }`$ on $`M_g^{n1}`$. Indeed:

$`{\displaystyle \frac{d\mathrm{\Psi }}{dt}}`$ $`=`$ $`_g\mathrm{\Psi },{\displaystyle \frac{dv}{dt}}_{T(M_g^{n1})}`$
$`=`$ $`_g\mathrm{\Psi },_g\mathrm{\Psi }_{T(M_g^{n1})}0,`$
$`{\displaystyle \frac{d\mathrm{\Psi }}{dt}}`$ $`=`$ $`\mathrm{S}p(DdP)/dt=\mathrm{S}p(D{\displaystyle \frac{dP}{dt}})`$
$`=`$ $`\mathrm{S}p(D,[[D,P],P])=\mathrm{S}p([P,D],[P,D])`$
$`=`$ $`\mathrm{S}p([P,D],[P,D]^{})0`$
for all $`t`$ , where $`.,._{T(M_g^{n1})}`$ is the scalar product on $`T(M_g^{n1})`$ obtained via the reduction of the scalar product $`.,.`$ on $`^n`$upon $`T(M_g^{n1})`$ under the constraint $`\underset{i=1}{\overset{𝑛}{}}v_i=1`$.

## 5 Multi-agent replicator system

In the above we have considered the case when $`rankP(t)=1,`$ $`t`$. It is naturally to study now the case when $`nrankP(t)=k>1`$ being evidently constant for all $`t`$ too. This means therefore that there exists some orthonormal vectors $`f_\alpha \mathrm{I}mP,`$ $`\alpha =\overline{1,k},`$ that

$$P(t)=\underset{\alpha =1}{\overset{k}{}}f_\alpha (t)f_\alpha (t)$$
(31)
for all $`t`$. Put now $`f_\alpha :=\{\pm \sqrt{v_i^{(\alpha )}}:i=\overline{1,n}\}ImP,`$ $`\underset{i=1}{\overset{n}{}}v_i^{_{(\alpha )}}=1`$ for all $`\alpha =\overline{1,k}.`$ The necessary orthonormality condition $`<f_\alpha ,f_\beta >=\delta _{\alpha \beta },`$ $`\alpha ,\beta =\overline{1,k},`$ can be automatically satisfied if one put $`f_{\alpha ,i}:=(\mathrm{exp}h)_{i,\alpha },i=\overline{1,n},`$ for a general skew-symmetric matrix $`h=h^{}`$ in the Euclidian vector space $`E^n.`$ Then, evidently, all signs at vector components $`\pm \sqrt{v_i^\alpha }=f_{\alpha ,i},i=\overline{1,n}`$ will be defined exactly for each $`\alpha =\overline{1,k}.`$ As a simple consequence of the above representation for stable vector $`f_\alpha E^n,\alpha =\overline{1,k}`$ one derives that in the multiagent case one can not write down the resulting replicator dynamics equations in the terms of positive concentration frequencies. The corresponding gradient flow on $`𝒫`$ then takes modified commutator form with the Lyapunov type function variation $`\delta \mathrm{\Psi }:𝒫`$ as follows:

$$\delta \mathrm{\Psi }=\underset{\alpha =1}{\overset{k}{}}Sp(D_\alpha \delta P^{(\alpha )}),$$
(32)
where by definition $`D_\alpha ^{}=D_\alpha endE^n,`$ and $`P^{(\alpha )}=f_\alpha (t)f_\alpha (t),`$ $`P^{(\alpha )}P^{(\beta )}=P^{(\alpha )}\delta _{\alpha \beta },`$ $`\alpha ,\beta =\overline{1,k},`$ $`\underset{\alpha =1}{\overset{k}{}}P^{(\alpha )}=P.`$ The resulting flow on $`𝒫`$ is given as

$$dP/dt=_\phi \mathrm{\Psi }(P),$$
(33)
where the gradient $`_\phi :𝒟(𝒫)T(𝒫)`$ is calculated taking into account the set of natural constraints:

$$\phi _{_{\alpha ,\beta }}(X;P):=Sp\left[X_{\alpha ,\beta }(P_\alpha P_\beta \delta _{\alpha ,\beta }P_\beta )\right]=0$$
(34)
for every $`\alpha \beta =\overline{1,k.}`$ Assume now that $`k=rankP(t)`$ for all $`t`$ as was stated before. Then the gradient flow (33) brings about

$$dP_\alpha /dt=[[D_\alpha ,P_\alpha ],P_\alpha ]\underset{\beta =1,\beta \alpha }{\overset{k}{}}(P_\alpha D_\alpha P_\beta +P_\beta D_\alpha P_\alpha ).$$
(35)
The flow (35) one can also obtain as a Hamiltonian one with respect to the following canonical symplectic structure on $`𝒫`$:

$$\omega ^{(2)}(P)=\underset{\alpha =1}{\overset{k}{}}Sp(P_\alpha dP_\alpha dP_\alpha P_\alpha )$$
(36)
with an element $`P=\underset{\alpha =1}{\overset{𝑘}{}}P_\alpha 𝒫`$, natural constraints $`P_\alpha P_\beta =P_\alpha \delta _{\alpha \beta },`$ $`\alpha ,\beta =\overline{1,k},`$ and some Hamiltonian function $`HD(P)`$ which has to be found making use of the expression

$$i_{_\phi \mathrm{\Psi }}\omega ^{(2)}=dH.$$
(37)
Straightforward but tedious calculations of (37) give rise to the same expression (35).

The system of equations (35) is a natural generalization of the replicator dynamics for description of a multi-agent interaction. They can describe for example economic communities trying to adapt to a changing environment. In this approach agents update their behavior in order to get maximum payoff under the given matrices of strategies $`\{a_{ik}^{(\alpha )}\}`$ in response to the information received from other agents. The structure of the system (35) is quite different from the system of equations (5). The first term on the right hand side describes the individual evolution of each economic agent $`\alpha `$in accordance to its own independent replicator dynamics and the second term describes average payoffs of all others $`k`$ agents except $`\alpha `$.