format: text

text:
# Robust creation of entanglement between ions in spatially separate cavities

## Abstract

We present a protocol that allows the generation of a maximally entangled state between individual atoms held in spatially separate cavities. Assuming perfect detectors and neglecting spontaneous emission from the atoms, the resulting idealized scheme is deterministic. Under more realistic conditions, when the the atom-cavity interaction departs from the strong coupling regime, and considering imperfect detectors, we show that the scheme is robust against experimental inefficiencies and yields probabilistic entanglement of very high fidelity.

The ability to reliably create entanglement between spatially separate parties is of paramount importance for the actual implementation of any quantum communication protocol Plenio V 98  and is also a pre-requisite for distributed quantum computation Cirac EHM 00 . Atoms or ions trapped inside optical resonators provide a promising set up for demonstrating the feasibility of quantum networking. Proposed ion trap quantum gates Gates  allow the coherent processing of quantum information stored in long lived electronic states. Indeed, sequential gate operation allowed the first quantum algorithm to be implemented in a linear ion trap Blatt 03 . Transferring quantum information between distant sites could be achieved by mapping the electronic degrees of freedom of the ions onto the photonic degrees of freedom of the cavity, which can then be used to transmit the quantum information to a distant site. First experimental progress towards this direction has been recently reported, demonstrating that individual ions can be positioned inside an optical resonator achieving sub-wavelength position control ionsincavities . The next step would be the controlled transfer of quantum information between electronic and photonic qubits, which should then be mapped out of the cavity.

However, once the photon has left the resonator through one of its highly reflecting mirrors it is not a straightforward task to feed it into another cavity. Ingenious schemes using careful pulse shaping have been devised to achieve this goal Repeater  but their experimental implementation remains challenging. A conceptually different approach consists of relaxing the condition that the quantum information is transferred via a photon leaving cavity A and entering cavity B.

Several schemes have been proposed for the generation entanglement between atoms, by detecting photons, in such a way that it is impossible to distinguish from which site they were emitted Cabrillo ; Bose KPV 99 ; Protsenko ; feng ; simon . For example, one could imagine a setting as in Fig. 1 where photons are allowed to leave both cavities and are then mixed on a beam-splitter prior to ordinary photon detection. Which-path information is destroyed in the beam-splitter and subsequent detection of a photon can then lead to a projection of the electronic degrees of freedom of the atom (which are entangled with the photonic degrees of freedom) onto some maximally entangled state Bose KPV 99 .

In its original formulation, protocol Bose KPV 99  employed sudden excitation of the ions, which, assuming otherwise perfect experimental conditions, limited its efficiency to $`50\%`$. Besides that constraint, two further problems would be difficult to overcome in practice. The first and most serious one is that the mapping between an ion and an optical cavity usually takes place within the weak coupling regime, defined by the relationship $`g^2/\kappa \gamma 1`$. Here $`g`$ is the ion-cavity coupling for a relevant set of atomic levels, $`\kappa `$ is the decay rate of the optical cavity and $`\gamma `$ denotes the spontaneous decay rate on the transition driven by the cavity mode. Within weak coupling, it is very likely that the atom will suffer an incoherent spontaneous emission, resulting in a photon leaving the cavity undetected to the sides before the electronic degree of freedom has been mapped onto the photonic degree of freedom. As it is very difficult to detect a photon that is emitted to the sides of the cavities, this event severely damages the quantum entanglement that one intends to create.

Additionally, most protocols assume perfect detectors while in practice they are generally not available. This problem is compounded by the fact that in a number of setups for optical cavities the mirrors possess considerable absorption which can be as high as $`50\%`$ of the photons that are not reflected from the cavity Hinds 03 . Therefore, any proposed scheme aimed to be demonstrated with current technology needs to be highly insensitive to detector inefficiencies. Many of the above problems, such as weak coupling, poor detector efficiencies or absorption in the mirrors, also occur, if one wishes to entangle two ions in a single optical cavity by detecting photons as they leak out of the mirrors. For this setting a number of schemes have been put forward recently, see e.g. Plenio HBK 99 ; Sorensen M 02 . In Plenio HBK 99  an entangled state between the ions is prepared conditional on the failure to detect a photon leaking from the cavity. In practice the fidelity of the state decreases very rapidly when one enters the weak coupling limit or when one has imperfect detectors or absorption in the mirrors. The second scheme Sorensen M 02  is more robust within the weak coupling regime, but requires single photons pulses and suffers strong loss of fidelity when faced with imperfect detectors or absorption in the mirrors Duanaswell .

In the following we present a scheme to entangle ions trapped individually in spatially separated cavities which (i) succeeds with $`100\%`$ probability under ideal conditions, (ii) allows the achievement of high fidelity entanglement outside the strong coupling regime upon the detection of a photon, (iii) is robust against detector inefficiencies and absorption losses in the cavity mirrors and (iv) can be adapted, with the same efficiency, to entangle ions trapped in a single optical cavity.

The method proposed here has its roots in the scheme presented in Bose KPV 99 , where a teleportation protocol between two cavities that employs the leakage of photons through the cavity mirrors was discussed. The same method can also be used to establish entanglement between the ions trapped in separate cavities. We briefly describe this approach here to illustrate its limitations and to motivate how to overcome them. Consider the set up depicted in Fig. 1, where each cavity contains a single trapped ion with an internal level structure which is given in Fig. 2. Light that may leak through the cavity mirrors is mixed on a $`50/50`$ beam splitter and subsequently observed by photo-detectors. The qubit is represented by the lower two energy levels which are coupled via a far detuned Raman-like transition. In Bose KPV 99  it was envisaged that the ions are both initially prepared in state $`|2`$. Then, identical far-detuned classical light pulses are applied to both ions such that, under ideal conditions, the state of the global system is given by

$`|\psi _{tot}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(|2_A,2_B|v_A,v_B|1_A,1_B|p_A,p_B`$ (1)
$`+`$ $`i(|2_A,1_B|v_A,p_B+|1_A,2_B|p_A,v_B)),`$
where $`|v_A`$ represents the vacuum state in cavity A and $`|p_B`$ denotes the one-photon Fock state in cavity B. Following this pulse, one waits to allow photons to leak through the cavity mirrors, mix at the beam splitter and reach the detectors. If a single click occurs, then the system is projected onto one of the two entangled states $`\{(|2_A,1_B+|1_A,2_B)/\sqrt{2},(|2_A,1_B|1_A,2_B)/\sqrt{2}\}`$. If no photon is detected or two photons are detected, then the ions are projected onto a product state and the procedure has failed.

Apart from the sensitivity of this scheme to losses due to spontaneous emission and detector inefficiencies, the procedure fails even under ideal conditions in $`50\%`$ of the cases. The reason for the $`50\%`$ failure rate of the scheme is that we excite both ions suddenly, which leads to a very high probability for the two photon detection event, leaving the ions in a product state. Furthermore, the scheme is not robust to spontaneous decay of the ions (this is particularly relevant in the weak coupling regime) or detector inefficiencies. If spontaneous emission occurs in one of the ions, the photon emitted will escape undetected, the detection of one photon in the photo-detectors after this event will then lead to the generation of the product state $`|1_A,1_B`$, not the desired entangled state. Secondly, if the detector is inefficient, then only one of the photons of a two photon event, i.e. a failure of the protocol, might be detected. Both these cases will therefore lead to a potentially severe reduction in fidelity of the states produced by the protocol.

However, an additional ingredient can make this scheme highly robust against any of these error sources. This is achieved by relaxing the condition of sudden excitation of the two ions and replacing it by more gentle driving. In the following we will show that weak driving, under ideal conditions, allows the scheme to succeed with arbitrarily high probability (and unit fidelity).

The Hamiltonian of the combined ion-cavity system, with ion internal level structure as given in Fig. 2, in a suitable interaction picture and setting $`\mathrm{}=1`$, a convention we will use throughout this paper, is given by

$`H`$ $`=`$ $`{\displaystyle \underset{i=A,B}{}}(\mathrm{\Delta }|3_i{}_{i}{}^{}3|+g|3_{i}^{}{}_{i}{}^{}1|c_i+g|1_{i}^{}{}_{i}{}^{}3|c_i^{}`$ (2)
$`+\mathrm{\Omega }|3_i_i2|+\mathrm{\Omega }|2_i_i3|),`$
where we have assumed that the two ions are subjected to identical laser fields on the $`|2|3`$ transition and that this laser field has the same detuning as the cavity field that couples to the $`|1|3`$ transition. The annihilation and creation operators for the cavity photons are denoted by $`c_i`$ and $`c_i^{}`$. The upper level $`|3`$ of both ions can decay to levels $`|1`$ and $`|2`$ with a rate of $`2\gamma _{31}`$ and $`2\gamma _{32}`$, respectively. Each of the cavities has a decay rate $`2\kappa `$. The full master equation for the density operator $`\rho `$ then takes the form

$`\dot{\rho }`$ $`=`$ $`i(H_{\text{eff}}\rho \rho H_{\text{eff}}^{})+2\kappa {\displaystyle \underset{i=A,B}{}}c_i\rho c_i^{}`$ (3)
$`+{\displaystyle \underset{i=A,B}{}}2\gamma _{31}|1_i_i3|\rho |3_i_i1|`$
$`+{\displaystyle \underset{i=A,B}{}}2\gamma _{32}|2_i_i3|\rho |3_i_i2|,`$
where we have defined the effective non-Hermitian Hamiltonian

$`H_{\text{eff}}=Hi\kappa {\displaystyle \underset{i=A,B}{}}c_i^{}c_ii(\gamma _{31}+\gamma _{32}){\displaystyle \underset{i=A,B}{}}|3_i_i3|.`$ (4)
This effective Hamiltonian will be used in the quantum jump approach to describe the system dynamics under the condition that neither a spontaneous emission nor a cavity photon have been detected Plenio K 98 .

For the sake of simplicity we now consider the system in the strong coupling limit, by setting $`\gamma _{31}=\gamma _{32}=0`$. We will later relax this condition to show that our method also works in the weak coupling limit. The requirement of weak driving means that the condition $`\frac{g\mathrm{\Omega }}{\mathrm{\Delta }}\kappa `$ is satisfied. Intuitively this implies that the rate of transitions between levels $`|1`$ and $`|2`$ of the ions will be weaker than the cavity decay. This in turn implies that the population in level $`|1`$ of the atoms will be small, unless a photon is detected. Indeed, after adiabatic elimination of the upper level $`|3`$ we obtain that the weak driving dynamics is governed by the new master equation pellizzari ,

$`\dot{\rho }`$ $`=`$ $`i(H_{\text{ad}}\rho \rho H_{\text{ad}}^{})+2\kappa {\displaystyle \underset{i=A,B}{}}c_i\rho c_i^{},`$ (5)
where we have defined

$$\begin{array}{cc}\hfill H_{\text{ad}}& =\underset{i=A,B}{}[\frac{g\mathrm{\Omega }}{\mathrm{\Delta }}(|2_i{}_{i}{}^{}1|c_i+h.c.)+\frac{g^2}{\mathrm{\Delta }}|1_i{}_{i}{}^{}1|\hfill \\ & +\frac{\mathrm{\Omega }^2}{\mathrm{\Delta }}|2_i_i2|i\kappa c_i^{}c_i].\hfill \end{array}$$
(6)
Under the condition that no detection has been registered, the time evolution is governed by $`H_{ad}`$. Given an initial state $`|2_A,2_B|v_A,v_B`$, the state of the systems will quickly approach the form

$`|\psi `$ $`=`$ $`|2_A,2_B|v_A,v_B+x(|2_A,1_B|v_A,p_B`$
$`+|1_A,2_B|p_A,v_B)+O(x^2),`$
where $`xi\frac{g\mathrm{\Omega }}{\mathrm{\Delta }\kappa }`$. Therefore, the rate $`R`$ at which one observes photons in one of the detectors is proportional to $`R4\kappa \left(\frac{g\mathrm{\Omega }}{\mathrm{\Delta }\kappa }\right)^2`$. Note that in this regime this rate is essentially unaffected by the level shifts in $`H_{\text{ad}}`$, thus it is not necessary for these to be compensated for. If one of the photo-detectors clicks, then a maximally entangled state has been prepared to a high precision and one switches off the lasers so that the entangled state is then preserved as the ions decouple from the cavity. The mean time before the first detection event will be

$$T_{\text{av}}\frac{\mathrm{\Delta }^2\kappa }{4(g\mathrm{\Omega })^2}.$$
(8)
In such a time interval, there is a small probability that two photons are detected, however this probability scales as $`|x|^42\kappa T_{\text{av}}\frac{1}{2}\left(\frac{g\mathrm{\Omega }}{\mathrm{\Delta }\kappa }\right)^2`$ and can therefore be made arbitrarily small in the limit of large detuning. We then observe that we can prepare a perfectly entangled state with arbitrarily high fidelity if we choose a sufficiently high detuning or sufficiently weak coupling strengths $`g`$ and $`\mathrm{\Omega }`$. As a consequence, by choosing a detuning that is very large, i.e. driving the $`|1|2`$ transition very slowly, we can ensure that any detection event is linked to a single photon and that therefore the fidelity of the prepared state will be very close to unity. This demonstrates our first claim that the scheme can achieve perfect fidelity and unit success probability.

However, this result is still only valid in the strong coupling limit as we have so far neglected the effect of finite $`\gamma _{31}`$ and $`\gamma _{32}`$. The strong coupling limit is not easy to achieve experimentally and it would be desirable to have a procedure that generates a very high fidelity entangled state even away from this limit, ideally with reasonable success probability. In the following we will demonstrate that our scheme can also successfully produce high fidelity maximally entangled states outside the strong coupling limit where we allow $`\frac{g^2}{\kappa \gamma _{31}}1`$ or even $`\frac{g^2}{\kappa \gamma _{31}}1`$ and $`\gamma _{32}`$ is allowed to be non-vanishing. In fact, in this case we can still achieve very high fidelities at good success rates. This observation is confirmed by a numerical simulation for the following choice of parameters $`\mathrm{\Omega }=g,\kappa =10g,\gamma _{31}=\gamma _{32}=0.1g`$ and $`\mathrm{\Delta }=20g`$ and a waiting time of $`T=100/g`$ we obtain numerically a fidelity of $`F=0.98`$ with a success probability of the scheme of $`p0.1`$. Indeed, analytically the success probability is approximately given by

$$p_{suc}\left(\frac{g\mathrm{\Omega }}{\mathrm{\Delta }\kappa }\right)^24\kappa T$$
(9)
where $`T`$ is the time one is willing to wait for the first detection. If no click is observed after this time, the experiment is deemed a failure and the systems will be e-prepared with both atoms in state $`|2`$.

The previous considerations still assume that the detection efficiency for photons that leak out of the cavity is unity. However, there are important sources of losses in experiments that make this assumption unrealistic. Firstly, there may be absorption in the mirrors themselves Hinds 03  and secondly the detectors may only have a finite efficiency. A scheme that can work in a practical environment should therefore also be robust against detector inefficiencies. Fortunately, the present method exhibits exactly such a robustness. In terms of the detector efficiency $`\eta `$ we find that the success probability simply scales linearly as

$$p_{\text{suc}}\left(\frac{g\mathrm{\Omega }}{\mathrm{\Delta }\kappa }\right)^24\kappa \eta T.$$
(10)
With falling detector efficiency, the fidelity of the resulting state will decrease, because it will now contain an admixture from events where two photons have been emitted from the cavities, but only a single one has been detected. Indeed, from this argument one expects a weak linear reduction of the fidelity.

In Fig. 3 we have plotted both the success probability and the achieved error (1-fidelity) for fixed $`\mathrm{\Omega }=g,\kappa =10g,\gamma _{31}=\gamma _{32}=0.1g`$ and $`\mathrm{\Delta }=20g`$ and a waiting time of $`T=100/g`$ against the detector efficiency. This figure confirms the approximate analytical formulas presented above and underlines that our scheme is robust against variations in the detector efficiencies.

A further experimental imperfection which must be considered is the presence of “dark counts”, i.e. when the detector fires although no light is incident upon it. Clearly, this will degrade to some extent the performance of all schemes which rely on a single detector click to generate an entangled state, and will lead, in general, to a loss in fidelity of the state produced. However, in the present scheme, the time-window in which a click due to a photon should occur is far shorter than the mean time between dark counts. For example, in Buttler  a dark count rate of approximately 1400 s<sup>-1</sup> is reported, thus the mean time between dark counts is on the order of ms. In the optical regime, the atom-cavity coupling $`g`$, detuning $`\mathrm{\Delta }`$, cavity decay rate $`\kappa `$ and the coupling with the classical field $`\mathrm{\Omega }`$ will all be at least on the order of MHz Kimble . Thus, using Eq. (8) one can estimate that $`T_{\text{av}}`$, the mean time before a proper click occurs in this scheme, is on the order of $`\mu `$s. Since the time-window for detection in this scheme can thus be made much smaller than the mean time between dark counts, their effect on this scheme can be made very small.

So far we have considered the case where we were faced with the task of entangling two spatially separated ion-cavity systems. The key ingredient in the detection scheme was the beam-splitter that erased the which-path information from the system, so that a photo-detection event could lead to entanglement between the cavities. However, the above method could also be used to entangle two ions trapped in a single cavity whose decay is monitored by a single photo-detector if the system is set-up such that the detection of a photon does not provide any information about which ion the photon was emitted from.

In summary, we have presented here an approach that, under ideal conditions, allows for the deterministic generation of perfect entanglement between individual ions in distant cavities. In the ideal scenario, the unit success probability also allows for the generalization of this scheme to the direct implementation of quantum gates. The scheme can be adapted to entangling multiple ions in a single optical cavity. Most importantly, the scheme is robust to realistic experimental imperfections, and in particular it allows for the probabilistic generation of high fidelity entanglement when operated within the weak coupling limit and monitored by inefficient detectors.

###### Acknowledgements.

The authors would like to thank P. L. Knight for helpful advice. This work was supported by the US-Army through grant DAAD 19-02-0161, the Nuffield Foundation and Hewlett Packard Ltd. via an EPSRC CASE award.


text:
# Transport coefficients of a mesoscopic fluid dynamics model

## I Introduction

Investigating the hydrodynamics of complex fluids is a challenging task because of the large number of degrees of freedom involved, the interplay between the flow and thermodynamic properties and the existence of widely varying length and time scales. Examples include the dynamics of polymers and biopolymers in solution, polymer melts, amphiphilic systems, particulate suspensions and chemically reacting fluids.

The complexity of these systems points to numerical approaches and it has become apparent that mesoscale algorithms such as dissipative particle dynamics, smoothed particle hydrodynamics and lattice Boltzmann algorithms are increasingly providing useful numerical tools to investigate complex fluids and, in particular, the role of hydrodynamics in their behaviour. These approaches are essentially novel ways of solving the Navier-Stokes equations designed to incorporate the most relevant microscopic physics of a given complex fluid.

Recently Malevanets and Kapral described one such mesoscopic approach which we shall term stochastic rotation dynamics. This method solves the thermohydrodynamic equations of motion by following the paths of particles moving in discrete time but continuous space. The algorithm has two major advantages. Firstly, because it is a particulate method, it is easy to couple it to a solute, say, polymers or colloids, which can be treated using a molecular dynamics algorithm. Secondly, the algorithm is rather simple and therefore analytical work is possible which can be aimed at increasing our understanding of mesoscopic approaches to the hydrodynamics of complex fluids.

Malevanets and Kapral mapped the stochastic rotation algorithm to a Boltzmann equation and showed that the model possesses an H-theorem. Ihle and Kroll extended the algorithm to correct for spatial correlation effects and gave approximate results for the viscosity and diffusion coefficient in two dimensions. Allahyarov and Gompper used a Poiseuille flow geometry to mesure the viscosity in three dimensions. Stochastic rotation dynamics has been applied to colloids in solution, dilute polymer solutions, amphiphilic systems, flow in porous media, binary fluid mixtures, cluster dynamics and flow around solid objects.

To understand as fully as possible how stochastic rotation dynamics can be used to simulate hydrodynamic behaviour it is helpful to increase our understanding of its properties. Therefore this paper presents an investigation of the viscosity and diffusion coefficient of the model in two and three-dimensions.

The algorithm is described in section $`2`$ following Malevanets and Kapral and Ihle and Kroll. We then show how it is possible to model a shear flow by imposing Lees-Edwards boundary conditions. The ability to simulate shear is important as it provides a computationally feasible way of measuring viscosities over a wide parameter range. It will also extend the range of applications of the approach to study viscoelastic effects such as the shear-thinning of a polymer solution.

In section $`3`$ we discuss the shear viscosity, $`\eta `$, of stochastic rotation dynamics. We argue that it is useful to divide the viscosity into two contributions, a kinetic and a collisional term. We first provide analytical expressions for the kinetic and collisional parts of the viscosity. We then compare these to numerical results and discuss their relative magnitudes.

Section $`4`$ deals with the coupling of a test particle to the solvent model. A massive particle obeys Langevin dynamics and we obtain an expression for the friction coefficient describing the coupling to the solvent. We stress that the friction coefficient is related to the collisional part of the viscosity, not to the total viscosity, a distinction unimportant in a real fluid, but pertinent here.

Finally a conclusion summarises the paper and suggests possibilities for future work.

## II THE MODEL

### A Stochastic rotation dynamics

The solvent is modelled by a large number $`N`$, typically $`10^5`$, of point-like particles of mass $`m`$ which move in continuous space with a continuous distribution of velocities, but discretely in time . The algorithm proceeds in two steps. In the first of these, a free streaming step, the positions of the solvent particles at time $`t`$, $`𝒓_i(t)`$, are updated simultaneously according to

$`𝒓_i(t+\delta t)=𝒓_i(t)+𝒗_i(t)\delta t`$ (1)
where $`𝒗_i(t)`$ is the velocity of a particle and $`\delta t`$ is the value of the discretized time step of the solvent which we take to be unity.

The second part of the algorithm is the collision step. The system is coarse-grained into $`L^d/a^d`$ cells of a $`d`$-dimensional lattice with side $`a`$. There is no restriction on the total number of solvent particles in each cell, (although the total number of particles in the system is conserved). Stochastic multi-particle collisions are performed within each individual cell, by rotating the velocity of each particle relative to the centre of mass velocity $`𝒗_{cm}(t)`$ of all the particles within that cell

$`𝒗_i(t+\delta t)=𝒗_{cm}(t)+𝑹\left(𝒗_i(t)𝒗_{cm}(t)\right).`$ (2)
$`𝑹`$ is a rotation matrix which rotates velocities by a fixed angle $`\alpha `$ around an axis generated randomly for each cell and at each time step.

The aim of the collision step is to transfer momentum between the particles while conserving the total momentum and energy of each cell. Because mass, momentum and energy are conserved locally the thermohydrodynamic equations of motion are captured in the continuum limit. Hence hydrodynamic interactions can be propagated by the solvent. Note, however, that any molecular details of the solvent are excluded: this allows the hydrodynamic interactions to be modelled with minimal computational expense.

The volume in phase space is invariant under both the free streaming and collision steps. Hence the system is described by a microcanonical distribution at equilibrium. Also it has been proved that the algorithm obeys an H-theorem.

A particularly useful feature of stochastic rotation dynamics is the ease with which hydrodynamic interactions can be “turned-off” thus replacing the hydrodynamic heat bath by a Brownian (random) heat bath. This is achieved by randomly interchanging the velocities of all the solvent particles after each collision step, thus relaxing the constraint of momentum conservation from a local to a global one. Accordingly the velocity correlations which result in hydrodynamic interactions disappear from the fluid although the equilibrium Maxwell-Boltzmann distribution is maintained. This greatly facilitates pinpointing the effects of the hydrodynamic interactions in a given situation.

### B Simulation parameters

Unless otherwise stated $`m=1`$, $`a=1`$ and the number of cells was $`L^3=32^3`$. The initial solvent distribution was generated by assigning positions randomly within the system with an average number of particles per unit cell $`\gamma =5`$. Thus the total number of particles was $`N=163,840`$. The velocities were assigned from a uniform distribution $`\left(v_{max}v_\alpha v_{max}\right)`$, $`\alpha =x,y,z`$ where $`v_{max}=\sqrt{\frac{3k_BT}{m}}`$. The distribution relaxed rapidly ($`t100`$) to the equilibrium Maxwell-Boltzmann form corresponding to the temperature $`T`$. Because of the long time-scales studied, any small non-zero net momentum in the system can eventually cause a large shift of a velocity distribution along the $`y`$ axis. To overcome this we perform a Galilean transformation to remove any net momentum from the system and then rescale the velocities to restore the correct temperature.

### C Lees-Edwards shear boundary condition

To achieve a steady shear flow we use Lees-Edwards boundary conditions, which are rather naturally incorporated into stochastic rotation dynamics. Our aim is to introduce a flow in the $`x`$-direction, with a constant velocity gradient along $`y`$. Periodic boundary conditions are applied along $`x`$ and $`z`$. Along $`y`$ the boundary conditions are modified so that if a particle crosses the lower boundary at ($`x,y=0,z`$) it returns to the upper boundary at the position ($`x=\left(x+ut\right)_{mod.L},y=L_,z`$) with velocity ($`v_x=v_x+u,v_y,v_z`$) where $`mod.L`$ denotes modulus L. Similarly if a particle leaves the upper boundary at ($`x,y=L,z`$), it reappears at the lower boundary at ($`x=\left(xut\right)_{mod.L},y=0,z`$) with velocity ($`v_x=v_xu,v_y,v_z`$). This algorithm gives a steady shear flow, and the system mimics one of infinite size.

When an external force is applied to the system, energy is constantly produced. Moreover a small net momentum can be introduced by the Lees-Edwards boundary condition for any system with a finite number of particles. To sustain a zero net momentum and temperature the appropriate shear velocity is subtracted from each particle velocity, a Galilean transformation is applied, the solvent velocities are rescaled to give a temperature $`T`$ and the shear velocity is restored. This procedure is carried out after each collision step.

### D Grid shift for periodic and shear boundary conditions

Ihle and Kroll  pointed out that at low temperatures the transport coefficients of stochastic rotation dynamics showed anomalies. This happens because at low temperatures particles in a given cell can remain in that cell over several time steps and participate in a series of collisions. Accordingly the molecular chaos assumption breaks down. However Ihle and Kroll  showed that it is possible to improve the behaviour of the algorithm by placing the cubic grid in a random position at each collision step.

In practice the easiest way to implement such a grid shift is to move all the particles with the same random vector with components in the interval $`[a/2,a/2]`$ before each collision step. After the collision step the particles are returned to their original positions. If the displacement takes the particles across the edge of the simulation box Lees-Edwards or periodic boundary conditions are applied as appropriate.

We note that the grid shift accelerates momentum transport between cells. As a result the system relaxes more quickly. It must also be taken into account when calculating transport coefficients such as the viscosity.

## III Viscosity

One of the most important parameters characterising a fluid is its viscosity. In this section we investigate the shear viscosity associated with the stochastic rotation model. The viscosity is usefully divided into two contributions. The kinetic viscosity results from the momentum transferred during the free streaming step, and the collisional viscosity results from the momentum transferred during the velocity rotations.

### A Calculation of the kinetic viscosity

To calculate the kinetic viscosity we consider a system undergoing shear with rate $`\dot{\gamma }=\frac{u_x\left(y\right)}{y}`$ in the $`x`$ direction and use a kinetic theory approach. Consider first two dimensions. On average the velocity profile is given by $`𝐯=(\dot{\gamma }y,0)`$. In the steady state there is a frictional force acting on any plane in the fluid perpendicular to $`y`$ which is given by the element of the stress tensor

$`\sigma _{xy}=\eta {\displaystyle \frac{u_x\left(y\right)}{y}}`$ (3)
where $`\eta `$ is the shear viscosity. This is the experimental definition of the viscosity. Our aim is to calculate the stress tensor during the streaming step: $`\sigma _{xy}=`$(flux of $`x`$-momentum crossing a plane of constant $`y`$).

For simplicity we take $`y=y_0=0`$. Consider a particle at position ($`x,y`$) with velocity ($`v_x,v_y`$). Particles will only cross the plane during a time step if they are moving towards the plane and have a velocity such that $`|v_y\delta t|`$ is greater than the distance to the plane. Therefore, the stress tensor can be written

$`\sigma _{xy}=`$ $`{\displaystyle \frac{\rho }{\delta t}}`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑v_x{\displaystyle _{\mathrm{}}^0}𝑑y{\displaystyle _{y/\delta t}^{\mathrm{}}}𝑑v_yv_xP(v_x\dot{\gamma }y,v_y)`$ (4)
$`+`$ $`{\displaystyle \frac{\rho }{\delta t}}`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑v_x{\displaystyle _0^{\mathrm{}}}𝑑y{\displaystyle _{\mathrm{}}^{y/\delta t}}𝑑v_yv_xP(v_x\dot{\gamma }y,v_y)`$ (5)
where $`P(v_x,v_y)`$ is the velocity probability distribution of particles in the rest frame of the fluid and $`\rho =\frac{m\gamma }{a^d}`$ is the mass density. $`\gamma `$ denotes the average number of particles per cell and $`d=2,3`$ is the dimensionality. By making the change of variable $`v_x^{}=v_x\dot{\gamma }y`$, and changing the order of integration this reduces to

$`\sigma _{xy}={\displaystyle \frac{\dot{\gamma }\rho \delta t}{2}}v_y^2\rho v_xv_y`$ (6)
where the averaging is over the probability distribution $`P`$. The first and second terms result from the shear profile and from induced correlations between $`v_x`$ and $`v_y`$ respectively. In the steady state the velocity distribution $`P`$ deviates from the Maxwell Boltzmann form because of these correlations.

To find the behaviour of $`v_xv_y`$ we calculate separately the effect of the streaming and collision operations. Consider the velocity distribution at the plane $`y=y_0`$. At time $`t`$, at $`y_0`$, particles with positive velocity $`v_y`$ come from $`y_0v_y\delta t`$ due to the streaming step, and correspondingly have a smaller $`x`$ velocity. Conversly particles from $`y>y_0`$ tend to have a higher $`x`$ velocity. This means that the velocity probability distribution is sheared by the streaming. That is $`P^{after}(v_x,v_y)=P^{before}(v_x+\dot{\gamma }v_y\delta t,v_y)`$. Averaging $`v_xv_y`$ over this new distribution and changing variables gives

$`v_xv_y^{after}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑v_x{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑v_yv_xv_yP(v_x+\dot{\gamma }v_y\delta t,v_y)`$ (7)
$`=`$ $`v_xv_y\dot{\gamma }\delta tv_y^2.`$ (8)
Thus the streaming operation changes the value of $`v_xv_y`$ by $`\dot{\gamma }\delta tv_y^2`$, and tends to make $`v_x`$ and $`v_y`$ increasingly anticorrelated.

Next we consider the effect of the collision step on the velocity distribution. Collisions tend to reduce correlations. To see this we consider a cell of $`n`$ particles and the two dimensional collision equation

$`𝒗(t+\delta t)=\left(\begin{array}{cc}cos\alpha & \pm sin\alpha \\ sin\alpha & cos\alpha \end{array}\right)\left(𝒗(t)𝒗_{𝒄𝒎}\right)+𝒗_{𝒄𝒎}`$ (11)
where the sign corresponds to the direction of the rotation axis. The centre of mass velocity $`𝒗_{𝒄𝒎}`$ can be divided into two contributions, $`𝒗_{𝒄𝒎}=\left(𝒗+\widehat{𝒗}\right)/n`$, the first from the test particle of velocity $`𝒗`$ and the second $`\widehat{𝒗}=_{i=1}^{n1}𝒗_𝒊`$ from the $`n1`$ particles sharing the same cell. We assume molecular chaos, that is the velocities of different particles are uncorrelated. Then, $`v_x\widehat{v}_y=\widehat{v}_xv_y=0`$ and $`\widehat{v}_x\widehat{v}_y=(n1)v_xv_y`$ and

$`v_x(t+\delta t)v_y(t+\delta t)`$ (12)
$`=`$ $`\left[1\left({\displaystyle \frac{n1}{n}}\right)\left(1\mathrm{cos}2\alpha \right)\right]v_x(t)v_y(t).`$ (13)
where we have averaged over the direction of the rotation axis. In the simulations $`n`$ is not constant. Therefore, we must consider density fluctuations. The probability of $`n`$ particles being in a given cell is given by the Poisson distribution $`P\left(n\right)=\frac{e^\gamma \gamma ^n}{n!}`$. The probability of a given particle being in a cell with a total of $`n`$ particles is $`\frac{nP(n)}{\gamma }`$, so taking an average over this distribution

$`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{nP(n)}{\gamma }}v_x(t+\delta t)v_y(t+\delta t)`$ (14)
$`=`$ $`\left[1\left({\displaystyle \frac{\gamma 1+e^\gamma }{\gamma }}\right)\left(1\mathrm{cos}2\alpha \right)\right]v_x(t)v_y(t)`$ (15)
$`=`$ $`f(\alpha ,\gamma )v_x(t)v_y(t).`$ (16)
During one time step we find that $`v_xv_y`$ is first reduced by the streaming operation and then multiplied by a constant factor in the collision operation. In the steady state it will therefore oscillate between two values and we have the self-consistency condition $`\left(v_xv_y\dot{\gamma }\delta tv_y^2\right)f=v_xv_y`$. Thus

$`v_xv_y={\displaystyle \frac{\dot{\gamma }\delta tf}{1f}}v_y^2.`$ (17)
Substitution into (6) gives

$`\sigma _{xy}=\rho \dot{\gamma }\delta tv_y^2\left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{f}{1f}}\right).`$ (18)
Using equipartition of energy, eqn. (16) and the definition of viscosity (3) yields

$`\eta _{kin}^{2D}={\displaystyle \frac{\gamma k_BT\delta t}{a^2}}\left[{\displaystyle \frac{\gamma }{(\gamma 1+e^\gamma )(1\mathrm{cos}2\alpha )}}{\displaystyle \frac{1}{2}}\right].`$ (19)
Applying the same proceedure in three dimensions gives

$`\eta _{kin}^{3D}=`$ (20)
$`{\displaystyle \frac{\gamma k_BT\delta t}{a^3}}\left[{\displaystyle \frac{5\gamma }{(\gamma 1+e^\gamma )(42\mathrm{cos}\alpha 2\mathrm{cos}2\alpha )}}{\displaystyle \frac{1}{2}}\right].`$ (21)
### B Calculation of the collisional viscosity

The collisional viscosity $`\eta _{col}`$ can also be calculated using a kinetic theory approach. Consider a cubic cell with side $`a`$. As in the case of the kinetic viscosity a shear rate $`\dot{\gamma }=\frac{u_x(y)}{y}`$ is applied along the $`x`$-axis. The system is divided into two subcells by a plane at $`y=y_0`$ where $`0y_0a`$. Let the subcell at $`y_0<ya`$ contain $`n_1`$ particles with an average collective flow velocity $`𝒖_\mathrm{𝟏}`$ and that at $`0<yy_0`$ contain $`n_2`$ particles with an average velocity $`𝒖_\mathrm{𝟐}`$. We have $`u_{1x}u_{2x}=\frac{n_1+n_2}{n_2}\left(u_{1x}v_{cmx}\right)`$ where $`v_{cmx}`$ denotes the $`x`$-component of the centre of mass velocity of the cell. As the average distance between the subcell velocities is $`\delta y=\frac{a}{2}`$ the shear rate is

$`\dot{\gamma }={\displaystyle \frac{u_x}{y}}`$ $`=`$ $`{\displaystyle \frac{u_{1x}u_{2x}}{\delta y}}`$ (22)
$`=`$ $`{\displaystyle \frac{2n}{a\left(nn_1\right)}}\left(u_{1x}v_{cmx}\right)`$ (23)
where $`n=n_1+n_2`$ is the number of particles in the cell.

Our aim is to calculate the momentum crossing the plane at $`y_0`$. We consider the $`n_1`$ particles at $`y_0<ya`$. The momentum transfer between the two subcells

$`\sigma _{xy}=\left({\displaystyle \underset{i}{}}p_{ix}(t+\delta t){\displaystyle \underset{i}{}}p_{ix}(t)\right)/\left(a^{d1}\delta t\right)`$ (24)
where $`i`$ runs over the $`n_1`$-particles, is calculated using the collision operation (2). Averaging over an isotropic distribution of the rotation axis gives

$`\sigma _{xy}={\displaystyle \frac{m}{a^{d1}\delta t}}\left[{\displaystyle \frac{2}{d}}n_1\left(1\mathrm{cos}\alpha \right)\left(u_{1x}v_{cmx}\right)\right]`$ (25)
where $`m`$ is the mass of the solvent particle. Using eqns. (23), (25) and the definition of the viscosity (3) we obtain

$`\eta _{col}={\displaystyle \frac{mn_1\left(nn_1\right)}{a^{d2}dn\delta t}}\left(1\mathrm{cos}\alpha \right).`$ (26)
As $`n_1+n_2`$ is in general small we must again consider fluctuations in the particle density. The numbers in the subcells, $`n_1`$ and $`n_2`$ are binomially distributed and averaging over them gives

$`\eta _{col}={\displaystyle \frac{m\left(1\mathrm{cos}\alpha \right)}{a^{d2}d\delta t}}(n1)({\displaystyle \frac{y_0}{a}})\left(1{\displaystyle \frac{y_0}{a}}\right).`$ (27)
Here we used $`n_1(nn_1)P(n_1)_{bin}=n^2pnpqn^2p^2`$ where $`p=1\frac{y_0}{a}`$ and $`q=\frac{y_0}{a}`$. Using the Poisson distribution to average over $`n`$

$`\eta _{col}`$ $`=`$ $`{\displaystyle \frac{m\left(1\mathrm{cos}\alpha \right)}{a^{d2}d\delta t}}({\displaystyle \frac{y_0}{a}})(1{\displaystyle \frac{y_0}{a}}){\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}\left(n1\right){\displaystyle \frac{e^\gamma \gamma ^n}{n!}}`$ (28)
$`=`$ $`{\displaystyle \frac{m\left(1\mathrm{cos}\alpha \right)}{a^{d2}d\delta t}}(\gamma 1+e^\gamma )({\displaystyle \frac{y_0}{a}})\left(1{\displaystyle \frac{y_0}{a}}\right)`$ (29)
where the summation runs from $`2`$ to $`\mathrm{}`$ because there is no momentum transfer unless $`n_11`$ and $`n_21`$. Finally, averaging over all planes $`0y_0a`$ gives

$`\eta _{col}={\displaystyle \frac{m\left(1\mathrm{cos}\alpha \right)}{6a^{d2}d\delta t}}\left(\gamma 1+e^\gamma \right).`$ (30)
In the large density limit the effect of the density fluctuations vanish and the theory agrees with the expression obtained by Ihle and Kroll for the two dimensional case ($`d=2`$). Also for case $`\alpha =\frac{\pi }{2}`$ and $`d=2`$ the result presented here is identical to that given in by Malevanets and Kapral. When $`d=3`$ we find that our expression differs from that given in  for $`\alpha =\frac{\pi }{2}`$.

### C Measuring the shear viscosity numerically

The computationally most efficient approach to measure the shear viscosity directly exploits the experimental definition of viscosity (3). A shear is applied along $`x`$ setting up a velocity gradient along $`y`$. The stress tensor $`\sigma _{xy}`$, equal to the flux of $`x`$-momentum across a plane perpendicular to the $`y`$ axis can easily be measured numerically.

Again we consider the momentum transfer in two parts, that due to the streaming step and that due to the collision step. A subtlety that needs to be taken into account occurs because of the finite size of the grid imposed in the collision step. This means that the streaming and collisional contributions to the momentum transfer depend on the position of the plane chosen for the measurement. For example, if the plane lies at the edge of a cell, there is no momentum transfer during the collision step whereas if the plane is taken at the centre of a cell the momentum transfer in the collision is maximised. The streaming contribution also varies, in such a way that the total momentum transfer in one simulation time step is independent of the position of the measuring plane. This is illustrated in figure 1. To separate out contributions to the kinetic and collisional viscosities it is necessary to average over all plane positions within a cell: this is another advantage of imposing the grid shift which performs the averaging over the cell automatically.

To measure the shear viscosity we used a simulation of size $`32\times 32\times 32`$. Four thousand time steps were required to reach the steady state and the stress tensor was averaged over ten thousand time steps. It is also necessary to average over many measurement planes to minimise thermal noise. We choose $`25`$ planes from the $`32`$ cubed box. Planes near the edge of the simulation box are omitted to avoid edge effects resulting from the Lees-Edwards boundary conditions. We also note that the pressure, $`P=\sigma _{yy}`$, can be measured in the same way, and takes the expected value for an ideal gas.

To check that the algorithm corresponds to a Newtonian fluid we plot in figure 2 the dependence of the shear viscosity on the shear rate. The viscosity is independent of the shear rate, as expected, except at very low and very high shear rates. The deviations at low shear are due to statistical errors because the shear represents only a small perturbation to the Maxwellian velocity distribution. Those at high shear result from finite-size effects. They occur when the distance moved by the walls in each time step approaches the size of the system. The value of $`\dot{\gamma }`$ chosen to measure the shear viscosity was $`0.0625`$ at $`k_BT=0.8`$.

The shear viscosity $`\eta _{GK}`$ can also be measured using an equilibrium approach based on the Green-Kubo formula

$`\eta _{GK}=\underset{t\mathrm{}}{lim}{\displaystyle \frac{V}{k_BT}}{\displaystyle _0^t}𝑑S\sigma _{xy}(S)\sigma _{xy}(0).`$ (31)
For discrete time steps this can be approximated by

$`\eta _{GK}=`$ (32)
$`{\displaystyle \frac{V\delta t}{k_BT}}\left({\displaystyle \frac{1}{2}}\sigma _{xy}(0)\sigma _{xy}(0)+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\sigma _{xy}(n\delta t)\sigma _{xy}(0)\right)`$ (33)
where $`\delta t`$ is the time step of the solvent dynamics.

Further details of the equilibrium approach is described in . However thermal fluctuations mean that very long runs are needed to obtain statistically meaningful results – using shear is faster by a factor $`30`$. Therefore this approach was used only as a check on the validity of the non-equilibrium method.

### D Results for the viscosity

Consider first the kinetic viscosity. Figures 3 and 4 present the dependence of the three dimensional kinetic viscosity on the rotation angle $`\alpha `$ at two different temperatures. Agreement between the theoretical result eqn (21) and the simulations is excellent. As $`\alpha `$ goes to zero the kinetic viscosity diverges. This corresponds to particles moving with an infinite mean free path, unimpeded by collisions. From eqn (16) and (17), as $`\alpha 0`$ the correlation $`v_xv_y\mathrm{}`$. In two dimensions the kinetic viscosity is symmetrical about $`\alpha =\frac{\pi }{2}`$ and hence also diverges at $`\alpha =\pi `$. This differs from the three dimensional case where the viscosity reaches a local maximum at $`\alpha =\pi `$.

A small deviation between the theory and the numerical results is observed at low temperatures. We believe that this is due to one of two effects. At such low temperatures, even with a gridshift, the assumption of molecular chaos is an approximation. Also at these temperatures the equilibrium distribution no longer approximates to a Maxwell-Boltzmann distribution due to the effect of the shear. Figure 5 compares the density dependence of the kinetic viscosity with theory, again showing excellent agreement.

We now turn to the collisional viscosity. This is as predicted, independent of temperature (eqn (30)). Figures 6 and 7 show its dependence on rotation angle $`\alpha `$ and particle density respectively. Agreement with the theory is again excellent except at extremely low temperatures.

The kinetic viscosity is proportional to the temperature and we recall that the collisional viscosity does not depend on the temperature. Hence the former dominates at high temperatures, and the latter at low temperatures. This has important indications for the balance of viscous and diffusive transport as discussed below.

Finally we check for any finite size effects in the simulations. Figure 8 shows that for the systems studied there is only a very small effect on the viscosity at the smallest system size.

## IV Coupling to a solute particle

To place solute particles such as polymers or colloids in the stochastic rotation dynamics solvent it is important to understand the coupling between them and the solvent. Therefore we calculate here the friction coefficient acting on a particle of mass $`M`$ and velocity $`𝒗`$ in a Brownian solvent, and then show that the fluctuation-dissipation theorem holds as expected.

### A Friction coefficient

The discretized Langevin equation for a monomer is

$`M{\displaystyle \frac{𝒗(t+\delta t)𝒗(t)}{\delta t}}=\xi 𝒗(t)+𝜼(t)`$ (34)
where $`\xi `$ denotes the friction and $`\delta t`$ is the time step. $`𝜼(t)`$ is stochastic noise with time average $`𝜼=0`$. The evolution of the velocity is given by eqn (2) which may be rewritten

$`𝒗(t+\delta t)𝒗(t)=\left(𝑹𝑰\right)\left(𝒗(t)𝒗_{𝒄𝒎}\right)`$ (35)
where $`𝑹`$ is the usual rotation matrix and $`𝑰`$ is the identity matrix.

We consider a particle of mass $`M`$ in a cell with $`n`$ solvent particles of mass $`m`$. The centre of mass velocity of the cell is given by

$`𝒗_{𝒄𝒎}={\displaystyle \frac{M𝒗(t)+_{i=1}^nm𝒖_𝒊}{M+mn}}.`$ (36)
Averaging eqn (35) over an isotropic distribution of the rotation axis and over $`𝒖_𝒊`$, assuming $`𝒖_𝒊=0`$, gives

$`𝒗(t+\delta t)𝒗(t)={\displaystyle \frac{2}{3}}{\displaystyle \frac{mn}{M+mn}}\left(1\mathrm{cos}\alpha \right)𝒗(t).`$ (37)
In time the number of particles within a given cell, $`n`$, changes. The probability is given by the Poisson distribution, $`P(n)=\frac{e^\gamma \gamma ^n}{n!}`$, where $`\gamma `$ is the average number per cell. By comparing (37) with the time averaged Langevin equation (34), we identify the friction coefficient as

$`\xi ^{3D}={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{e^\gamma \gamma ^n}{n!}}{\displaystyle \frac{mn}{M+mn}}{\displaystyle \frac{2M}{3\delta t}}\left(1\mathrm{cos}\alpha \right).`$ (38)
In two dimensions

$`\xi ^{2D}={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{e^\gamma \gamma ^n}{n!}}{\displaystyle \frac{mn}{M+mn}}{\displaystyle \frac{M}{\delta t}}\left(1\mathrm{cos}\alpha \right).`$ (39)
### B Numerical results for the friction coefficient

The friction coefficient is measured numerically using eqn (34). Multiplying (34) by $`𝒗(t)`$ and taking a time average

$`𝒗(t+\delta t)𝒗(t)𝒗^2(t)={\displaystyle \frac{\delta t}{M}}\xi 𝒗^2(t)`$ (40)
where we use $`𝜼(t)𝒗(t)=𝜼(t)𝒗(t)=0`$. Therefore, using the equipartition theorem

$`\xi ={\displaystyle \frac{M}{\delta t}}\left(1{\displaystyle \frac{M}{dk_BT}}𝒗(t+\delta t)𝒗(t)\right)`$ (41)
Theory and numerical results are in good agreement for different $`M`$ as shown in figure 9.

### C Fluctuation-dissipation theorem

Next we verify fluctuation-dissipation theorem. Substituting (35) into (34) and rearranging gives

$`𝜼(t)={\displaystyle \frac{M}{\delta t}}\left(𝑹𝑰\right)\left(𝒗(t)𝒗_{𝒄𝒎}\right)+\xi 𝒗(t).`$ (42)
To calculate the equal time correlation function we must assume equipartition of energy and molecular chaos, therefore $`v_\alpha u_{i\beta }=v_\alpha u_{i\beta }=0`$ and $`u_{i\alpha }u_{j\beta }=\delta _{ij}\delta _{\alpha \beta }k_BT/m`$. Using eqn (36) and averaging over an isotropic distribution of the rotation axis and the solvent velocity distribution we find

$`𝜼(t)𝜼(t)`$ $`=`$ $`{\displaystyle \frac{2d\xi k_BT}{\delta t}}\left(1{\displaystyle \frac{mn\left(1\mathrm{cos}\alpha \right)}{d\left(M+mn\right)}}\right)`$ (43)
for dimension $`d`$. The first term in the bracket is that obtained from the fluctuation-dissipation theorem. The second is a small correction, which goes to zero in the limit of large $`M`$.

### D Diffusion coefficient

The diffusion coefficient can be calculated using the definition

$`D`$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{\left(𝒓\left(t\right)𝒓\left(0\right)\right)^2}{2dt}}`$ (44)
$`=`$ $`{\displaystyle \frac{1}{2d}}𝒗\left(0\right)𝒗\left(0\right)\delta t+{\displaystyle \frac{1}{d}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}𝒗\left(n\delta t\right)𝒗\left(0\right)\delta t.`$ (45)
We note that the resulting equation (45) is the same as the Green-Kubo formula, which is derived from linear response theory.

For a Brownian particle the velocity-velocity autocorrelation is easily calculated using eqn (40)

$`𝒗(t+n\delta t)𝒗(t)=\left(1{\displaystyle \frac{\delta t}{M}}\xi \right)^n𝒗^2(t).`$ (46)
Substitution into (45) gives,

$`D={\displaystyle \frac{k_BT}{\xi }}\left(1{\displaystyle \frac{\delta t\xi }{2M}}\right).`$ (47)
This is the normal Einstein relation with an error term which becomes vanishingly small in the large $`M`$ limit.

## V Conclusion

We have provided analytic expressions for the shear viscosity of a two- and three-dimensional mesoscale model of hydrodynamics, stochastic rotation dynamics. The derivation assumes molecular chaos. The viscosity has two contribution $`\eta _{kin}`$ and $`\eta _{col}`$ from the streaming and collision steps of the algorithm respectively. The former dominates at high temperatures and the latter at low temperatures.

We also obtained an expression for the viscous drag $`\xi `$ acting on a test particle. Note that $`\xi `$ and $`\eta _{col}`$ given by eqns (38) and (30) have the same dependence on the rotation angle $`\alpha `$ and are both independent of temperature. This is to be expected as they both result from momentum transfer in the collision step of the algorithm. It is useful to realise that in the stochastic rotation method for a fluid with a given total viscosity the coupling to a test particle will depend on temperature and will be very weak at high temperatures where $`\eta _{kin}\eta _{col}`$.

The viscosity in three dimensions was measured numerically by using Lees-Edwards boundary conditions to apply a shear to the system. The analytic and numerical results were found to be in excellent agreement thus providing a check on both calculations. It is important to have a robust way of measuring transport coefficients for this model as it can easily be coupled to a solute such as polymers or colloids where the non-Newtonian flow properties of the resultant solution are of interest.