format: text

text:
# Regular binary thermal lattice-gases

## I Introduction

Lattice-gas automata (LGA), as introduced by Frisch, Hasslacher, and Pomeau Frisch et al. (1986), are ideal testing systems for kinetic theory. Although they have a simple structure, which makes them extremely efficient simulation tools, they still address the full many-body particle problem. Much of the efficiency originates from the discretization of the positions and velocities of the point-like particles onto a lattice. The simplified dynamics is a cyclic process: a streaming step, where all particles propagate to neighboring lattice sites, followed by a local collision step. The collisions typically conserve mass and momentum, and in addition energy if the model under consideration is thermal.

Lattice-gases are capable of simulating macroscopic fluid flow Frisch et al. (1986); Rivet and Boon (2001), and can be used for studying flow through porous media Chen et al. (1991), immiscible multicomponent fluids Gunstensen and Rothman (1991), reaction-diffusion Chopard et al. (1994), von Karman streets Malevanets and Kapral (1998), Rayleigh-Bénard convection Baran et al. (1997), interfaces and phase transitions Rothman (1994). However, there are practical problems when using LGA’s: the models are not Galilean-invariant, temperature is not well-defined, the transport coefficient have unexpected behavior as a function of temperature and density, and many models contain spurious invariants Brito et al. (1991); Das and Ernst (1992); Ernst and Das (1992). It is therefore difficult to make connections with realistic systems. Although in the early years the LGA’s were thought to seriously compete with other fluid flow solvers, the main line of research has shifted to a testing ground for concepts in kinetic theory. The main motivation for this paper lies in the exploration of the limits of thermal lattice-gas capabilities concerning diffusion-like phenomena.

Diffusion can be incorporated in LGA’s in several ways. The computational most efficient method is the *color mixture* Rivet and Boon (2001); Ernst and Dufty (1990); Hanon and Boon (1997), which we have analyzed in detail for a thermal model Blaak and Dubbeldam (2001a). The otherwise identical particles are in such a mixture painted with a probability corresponding to the concentration and a color-blind observer would not notice any difference, i.e. all transport properties are the same as for the single-component fluid. But in addition there is an extra hydrodynamical and color-dependent diffusion mode that does not interfere with the other modes.

In this paper we consider a different way to include diffusion, namely the regular binary mixture Ernst and Dufty (1990). In the regular binary mixture one has two distinct species of particles, i.e. different mass. For particles of the same species the exclusion principle holds and hence there can be at most one particle of a given species in any velocity channel. There is, however, no mutual exclusion for particles off different species. Consequently, a given velocity channel can be occupied by two particles, provided they both belong to a different species.

One can interpret this as each species living on its own but identical lattice with exclusion. Since each specie is restricted to its own lattice there can be no mass exchange between the two lattices and there is local mass conservation for both species. Upon interaction, however, particles of both lattices corresponding to the same node are able to exchange momentum and energy, provided this does no violate local mass conservation and exclusion for each species. Consequently, the dynamics of the regular binary mixture is much more involved than that of the color mixture, but it is closer to the dynamics of real fluids. For simplicity we here will assume the special case of the two species having the same physical properties, i.e. the same mass, and distinguish both species by a different color label.

We analyze the model at the microscopic level and focus mainly on the behavior of the power spectrum. The coupling between energy transport and diffusion that arises in these systems, results in a more complicated structure of the power spectrum Mountain and Deutch (1969); Boon and Yip (1980); Blaak and Dubbeldam (2001b), in which the central peak now contains combined effects of both entropy fluctuations and concentration fluctuations. Macroscopically this coupling manifests itself as the Dufour effect (a concentration gradient induces a heat-flow) and the Soret effect (a temperature gradient induces a diffusion flux).

The remainder of this paper is organized as follows. We start by introducing the regular binary thermal lattice gas model in section II and use the molecular chaos assumption to obtain the linearized collision operator. Section III is concerned with a mode analysis of the linearized system, revealing the appearance of an extra diffusion mode. The modes related to thermal diffusivity and mass-diffusion are coupled and form two non-propagating, diffusive modes. Furthermore, in section IV we derive a Landau-Placzek analogue, a formula which describes the power spectrum in the hydrodynamic limit (long wavelength, small frequency domain), which we check with simulation results performed on this system in section V. In the final section we give a brief overview of the main results and some concluding remarks.

## II The regular binary GBL-model

The LGA model we like to consider here is a thermal model consisting of two different interacting species of particles with identical mass. They can be thought of to live on separate two-dimensional hexagonal lattices on which they propagate independently, but red and blue particles corresponding to the same node interact during the collision. The velocities are discretized and have a spatial layout as shown in figure 1. They are distributed over four rings: one ring contains a single null-velocity, and three rings each contain six velocities with magnitude $`1`$, $`\sqrt{3}`$, and $`2`$. The velocity set is of size $`19`$, and equivalent to the GBL-model, proposed by Grosfils, Boon, and Lallemand Grosfils et al. (1992). Since the model under consideration is basically a combination of two GBL models, particles of different type (red and blue) can have the same velocity. However, there can not be more than a single particle of each type in a given velocity state. The multiple rings correspond each to a different energy level and are a necessary requirement in order to introduce thermal properties. This particular velocity set guarantees the absence of spurious invariants and results in macroscopically isotropic behavior Grosfils et al. (1993).

The state of a node can be specified by a set of boolean occupation numbers $`n_{i\mu }`$, denoting the presence or absence of a particle of type $`\mu =\{r,b\}`$ in velocity channel $`𝒄_i`$, where $`i`$ is a label running over all 19 velocities. Due to the boolean nature of the LGA, the ensemble average of the occupation numbers $`f_{i\mu }`$ in equilibrium, is described by a Fermi-Dirac distribution Rivet and Boon (2001)

$$f_{i\mu }n_{i\mu }=\frac{1}{1+e^{\alpha _\mu +\frac{1}{2}\beta 𝒄_i^2𝜸𝒄_i}},$$
(1)
where $`\alpha _r`$, $`\alpha _b`$, $`\beta `$, and $`𝜸`$ are Lagrange multipliers and fixed by setting the value of the average red density $`\rho _r=_{ir}f_{ir}`$, the average blue density $`\rho _b=_{ib}f_{ib}`$ (or alternatively the total density $`\rho =\rho _r+\rho _b`$ and the fraction $`P_r`$ of red particles), the average velocity $`\rho 𝒖=_{i\mu }f_{i\mu }𝒄_i`$, and the energy density $`\rho e=\frac{1}{2}_{i\mu }f_{i\mu }𝒄_i^2`$. Here, $`\beta `$ is the inverse temperature, $`\alpha _r`$ and $`\alpha _b`$ fulfill a chemical potential role, and $`𝜸`$ is a parameter conjugate to the flow velocity. In the remainder of this paper we will restrict ourselves to the global zero-momentum case by putting $`𝜸=0`$.

The lattice-gas Boltzmann equation, describing the time-evolution of the average occupation numbers, is given by Rivet and Boon (2001)

$$f_{i\mu }(𝒓+𝒄_i,t+1)=f_{i\mu }(𝒓,t)+\mathrm{\Delta }_{i\mu }(f),$$
(2)
where the nonlinear collision term $`\mathrm{\Delta }_{i\mu }(f)`$ is a summation over all pre-and post-collision states $`s`$ and $`s^{}`$

$$\mathrm{\Delta }_{i\mu }(f)=\underset{s,s^{}}{}P(s)A(ss^{})(s_{i\mu }^{}s_{i\mu }).$$
(3)
The collision rules that are used are taken in to account by the transition matrix $`A`$, which contains the probability that on collision a state is transformed into an “equivalent” state. We denote the collection of states that can be transformed into each other by an equivalence class $`\mathrm{\Gamma }=(M_r,M_b,𝑷,E)`$, i.e. a class having the same red mass $`M_r=_is_{ir}`$, the same blue mass $`M_b=_is_{ib}`$, the same momentum $`𝑷=_{i\alpha }s_{i\alpha }𝒄_i`$, and the same energy $`E=_{i\alpha }\frac{1}{2}s_{i\alpha }𝒄_i^2`$. We adopt maximal collision rules, i.e. under collision a state can transform in any other state within the same equivalence class, including itself, thus fulfilling all conservation laws. Hence $`A(\mathrm{\Gamma })=\frac{1}{|\mathrm{\Gamma }|}`$, where we use $`|\mathrm{\Gamma }|`$ to denote the number of elements in the class $`\mathrm{\Gamma }`$. The density and energy density enter the collision operator through $`P(s)`$, the probability of occurrence of state $`s`$ in equilibrium. In the Boltzmann approximation, the velocity channels are assumed to be independent, hence $`P(s)`$ can be written as

$$P(s)=\underset{i\mu }{}f_{i\mu }^{s_{i\mu }}(1f_{i\mu })^{1s_{i\mu }}.$$
(4)
If the velocity fluctuations are sufficiently small a Taylor expansion of the collision term in the neighborhood of the equilibrium distribution is justified Rivet and Boon (2001), yielding the linearized collision operator $`\mathrm{\Omega }`$

$$(\mathrm{\Omega }\kappa )_{i\mu ,j\nu }=\underset{s,s^{}}{}P(s)A(ss^{})(s_{i\mu }^{}s_{i\mu })s_{j\nu },$$
(5)
where the diagonal matrix $`\kappa _{i\mu ,j\nu }`$ is determined by $`\kappa _{i\mu ,i\mu }=f_{i\mu }(1f_{i\mu })`$, which is the variance in the occupation number of a channel with the corresponding labels. It can easily be checked from Eq. (5) that $`\mathrm{\Omega }\kappa `$ is a symmetric $`38\times 38`$ matrix. In contrast with the colored LGA Blaak and Dubbeldam (2001a), $`\kappa `$ is a 38-dimensional diagonal matrix due to the fact that there is no mutual exclusion of a red and blue particle with the same velocity. It also enables us to introduce the colored thermal scalar product Hanon and Boon (1997); Blaak and Dubbeldam (2001a); Ernst and Das (1992),

$`A|B`$ $`={\displaystyle \underset{i\mu }{}}A(𝒄_{i\mu })B(𝐜_{i\mu })\kappa _{i\mu },`$ (6)
$`A|\mathrm{\Omega }|B`$ $`={\displaystyle \underset{i\mu j\nu }{}}A(𝒄_{i\mu })(\mathrm{\Omega }\kappa )_{i\mu ,j\nu }B(𝒄_{j\nu }),`$ (7)
where we adopt the convention that the matrix $`\kappa `$ is attached to the right vector.

## III Perturbation theory

In first approximation the behavior of the system can be obtained by considering the fluctuations $`\delta f_{i\nu }(𝒓,t)=f_{i\nu }(𝒓,t)f_{i\nu }`$ and analyze these deviations from the uniform equilibrium state in terms of eigenmodes. The solutions of the Boltzmann equation (2) are then of the form

$$\delta f_{i\nu }=\psi _\mu (𝒌\mathbf{,}𝒄_𝒊)e^{ı𝒌𝒓+z_\mu (𝒌)t},$$
(8)
where $`z_\mu (𝒌)`$ represents the eigenvalue of the mode $`\psi _\mu `$ at wavevector $`𝒌`$. The hydrodynamic modes are related to the collisional invariants and satisfy $`z_\mu 0`$ for $`𝒌0`$. There are five independent collision invariants $`a_n`$ corresponding to the conservation of the red and blue mass, the conservation of momentum, and the conservation of energy

$$a_n|\mathrm{\Omega }=0\mathrm{\Omega }|a_n=0,$$
(9)
$$|a_n=\{|R,|B,|c_x,|c_y,|\frac{1}{2}c^2\},$$
(10)
where $`|R_{i\alpha }=\delta _{\alpha r}`$, $`|B_{i\alpha }=\delta _{\alpha b}`$, and can be combined to form $`|\rho =|R+|B`$. It turns out that an equivalent and more convenient set of invariants is given by the following five combinations which are mutual orthogonal with respect to the colored thermal scalar product (6)

$$|a_n=\{|c_x,|c_y,|p=|\frac{1}{2}c^2,|s=|pc_s^2|\rho ,|d\},$$
(11)
where $`|s`$ is the analogue to the microscopic entropy for the GBL model Grosfils et al. (1993) and $`c_s^2=\frac{p|p}{p|\rho }`$ the speed of sound following from the orthogonality requirement $`p|s=0`$. The remaining invariant $`|d`$ is, contrary to the one in the colored GBL model Blaak and Dubbeldam (2001a), not simply the weighted difference between the red and blue densities $`\frac{1}{R|R}|R\frac{1}{B|B}|B`$, but is given by

$$|d=\left(\frac{|R}{s|R}\frac{|B}{s|B}\right)\left(\frac{p|R}{s|R}\frac{p|B}{s|B}\right)\frac{|p}{p|p}.$$
(12)
This form is determined by the requirement that it is orthogonal to the four other invariants. In the case of equal red and blue density the last term on the righthand side vanishes, but in general this is not the case. The origin of this term is in fact the absence of equipartition in this model due to the exclusion principle, which causes the ratio $`f_{ir}/f_{ib}`$ to be dependent on the velocity of the particles as can be seen from Eq. (1). Note that this last invariant is also perpendicular to the density, i.e. $`\rho |d=0`$. Other useful relations between these invariants are $`\rho |p=\rho |c_x^2=\rho |c_y^2`$ and $`p|p=p|c_x^2=p|c_y^2`$.

The choice of $`|s`$ and $`|d`$ made here is based on the correspondence of the definition of $`|s`$ with the one made for the GBL model and its colored counterpart Blaak and Dubbeldam (2001a); Grosfils et al. (1993). Since $`|s`$ nor $`|d`$, with the exception of some special limits, is going to be the zeroth order of an eigenmode the choice is somewhat arbitrary and any two linear combinations that are mutual orthogonal could be used as well. The current choice, however, will facilitate us to make a connection with continuous theory and the proper transport coefficients.

Note that this is similar to what is done in the case of continuous theory and generalized versions of hydrodynamics (See Ref. Boon and Yip (1980) and references therein), where one also needs to obtain a set of variables that are statistically independent based on thermodynamic fluctuation theory. This does not uniquely fix these variables, and the freedom that remains can be used in order to select an appropriate, orthogonal set for the specific problem.

Following the method of Résibois and Leener Résibois and de Leener (1977) we need to find the $`𝒌`$-dependent eigenfunctions and eigenvalues of the single-time step Boltzmann propagator

$$e^{ı𝒌𝒄}(\mathrm{𝟏}+\mathrm{\Omega })|\psi (𝒌)=e^{z(𝒌)}|\psi (𝒌),$$
(13)
$$\varphi (𝒌)|e^{ı𝒌𝒄}(\mathrm{𝟏}+\mathrm{\Omega })=e^{z(𝒌)}\varphi (𝒌)|,$$
(14)
where $`e^{ı𝒌𝒄}`$ has to be interpreted as a diagonal matrix, and $`\mathrm{𝟏}`$ is the identity matrix. The symmetries of the matrices cause the left and right eigenvectors to be related by $`\varphi _\mu (𝐤)=e^{ı𝒌𝒄}\psi _\mu (𝐤)/_\mu `$, and form a complete biorthonormal set

$$\underset{\mu }{}|\psi _\mu \varphi _\mu |=\mathrm{𝟏}\varphi _\mu |\psi _\nu =\delta _{\mu \nu },$$
(15)
where we used $`\mu `$ and $`\nu `$ to label the different eigenfunctions and introduced the normalization constants $`_\mu `$.

To obtain the hydrodynamic modes characterized by $`z(𝒌)0`$ in the limit $`𝒌0`$, we make a Taylor expansion of the eigenfunctions and eigenvalues

$$\psi _\mu (𝒌)=\psi _\mu ^{(0)}+(ık)\psi _\mu ^{(1)}+(ık)^2\psi _\mu ^{(2)}+\mathrm{},$$
(16)
$$z_\mu (𝒌)=(ık)z_\mu ^{(1)}+(ık)^2z_\mu ^{(2)}+\mathrm{},$$
(17)
where we already used that $`z_\mu ^{(0)}=0`$. The functions $`\psi _\mu (𝒌)`$ will be determined up to a normalization factor. As this normalization is not allowed to be observable, and in fact will not be observable, we can choose it in a convenient way by

$$\psi _\mu ^{(0)}|\psi _\mu (𝒌)=\psi _\mu ^{(0)}|\psi _\mu ^{(0)},$$
(18)
which leads to

$$\psi _\mu ^{(0)}|\psi _\mu ^{(n)}=\delta _{n0}\psi _\mu ^{(0)}|\psi _\mu ^{(0)}.$$
(19)
Substitution of the expansions in Eq. (13) and grouping according to the same order in $`k`$ gives

$$\mathrm{\Omega }|\psi _\mu ^{(0)}=0,$$
(20)
$$\mathrm{\Omega }|\psi _\mu ^{(1)}=(c_{\mathrm{}}+z_\mu ^{(1)})|\psi _\mu ^{(0)},$$
(21)
$$\mathrm{\Omega }|\psi _\mu ^{(2)}=(c_{\mathrm{}}+z_\mu ^{(1)})|\psi _\mu ^{(1)}+\left[z_\mu ^{(2)}+\frac{1}{2}(c_{\mathrm{}}+z_\mu ^{(1)})^2\right]|\psi _\mu ^{(0)},$$
(22)
with $`c_{\mathrm{}}=\widehat{𝒌}𝒄`$ and its orthogonal counterpart $`c_{}=\widehat{𝒌}_{}𝒄`$.

The solution of the zeroth order equation is straightforward and gives

$$|\psi _\mu ^{(0)}=\underset{n}{}A_{\mu n}|a_n,$$
(23)
with some unknown coefficients $`A_{\mu n}`$. Substitution in the first order equation and multiplying with $`a_m|`$ at the left gives for each $`m`$ a linear equation in the $`A_{\mu n}`$

$$\underset{n}{}A_{\mu n}a_m|c_{\mathrm{}}+z_\mu ^{(1)}|a_n=0,$$
(24)
where we used that $`a_m|\mathrm{\Omega }=0`$. This set of equations has a solution provided the determinant is zero, leading to the following eigenvalue equation for $`z_\mu ^{(1)}`$

$$\left(z_\mu ^{(1)}\right)^3\left(p|pc_{\mathrm{}}|c_{\mathrm{}}\left(z_\mu ^{(1)}\right)^2p|c_{\mathrm{}}^2\right)=0,$$
(25)
where it is used that combinations containing an odd power of $`c_{\mathrm{}}`$ or $`c_{}`$ are necessarily zero. This determines two of the five zeroth order eigenfunctions

$$|\psi _\sigma ^{(0)}=|p+\sigma c_s|c_{\mathrm{}}z_\pm ^{(1)}=\sigma c_s,$$
(26)
where $`\sigma =\pm `$ denotes opposite directions parallel to $`𝒌`$ and

$$\left(z_\sigma ^{(1)}\right)^2=\frac{p|c_{\mathrm{}}^2^2}{p|pp|\rho }=c_s^2.$$
(27)
For the remaining three eigenfunctions we can only conclude at this level that they satisfy $`z_\mu ^{(1)}=0`$ and are formed by combinations of $`|d`$, $`|s`$, and $`|c_{}`$ only. On physical grounds it can be argued that $`|c_{}`$ will be a mode on itself, however, it also will follow in a natural way later on.

The combination $`|j_\mu =c_{\mathrm{}}+z_\mu ^{(1)}|\psi _\mu ^{(0)}`$ is called the current. We will extend the definitions of the current $`j_\mu `$ and $`\psi _\mu ^{(0)}`$ to include $`s`$ and $`d`$ as a possible value for $`\mu `$. This is for convenience, since strictly speaking $`|s`$ and $`|d`$ will in general not be the source of an eigenfunction, but both will be a linear combination of two “true” modes of the system.

With the introduction of these currents it follows immediately from Eqns. (9) and (21) that invariants and currents are orthogonal

$$\psi _\nu ^{(0)}|j_\mu =0,$$
(28)
hence the currents lie in the complement of the null-space. Since this is not affected by applying $`\mathrm{\Omega }`$ we obtain

$$\psi _\nu ^{(0)}|\mathrm{\Omega }^m|j_\mu =0$$
(29)
for any integer $`m`$, including negative values. Note that for negative values of $`m`$ the expression $`\mathrm{\Omega }^m|j_\mu `$ is uniquely determined by the requirement that the result has no component in the null-space. Therefore the solution of the first order equation (21) can formally be written as

$$|\psi _\mu ^{(1)}=\frac{1}{\mathrm{\Omega }}|j_\mu +\underset{\nu }{}B_{\mu \nu }|\psi _\nu ^{(0)},$$
(30)
with yet to be determined coefficients $`B_{\mu \nu }`$ and also three still unknown currents. The coefficients $`B_{\mu \mu }=0`$ as follows from the chosen normalization (19) of the eigenfunctions.

Multiplying the second order equation (22) with $`\psi _\lambda ^{(0)}|`$ at the left side we find

$`0=`$ $`\psi _\lambda ^{(0)}|(c_{\mathrm{}}+z_\lambda ^{(1)})|\psi _\mu ^{(1)}+(z_\mu ^{(1)}z_\lambda ^{(1)})\psi _\lambda ^{(0)}|\psi _\mu ^{(1)}+`$ (31)
$`z_\mu ^{(2)}\psi _\lambda ^{(0)}|\psi _\mu ^{(0)}+{\displaystyle \frac{1}{2}}\psi _\lambda ^{(0)}|(c_{\mathrm{}}+z_\lambda ^{(1)})(c_{\mathrm{}}+z_\mu ^{(1)})|\psi _\mu ^{(0)}+`$
$`{\displaystyle \frac{1}{2}}(z_\mu ^{(1)}z_\lambda ^{(1)})\psi _\lambda ^{(0)}|(c_{\mathrm{}}+z_\mu ^{(1)})|\psi _\mu ^{(0)}.`$
Using the definition of the currents, substitution of the formal solution of the first order equation (30), and using the orthogonality relations (28) this leads to

$$z_\mu ^{(2)}\psi _\lambda ^{(0)}|\psi _\mu ^{(0)}+(z_\mu ^{(1)}z_\lambda ^{(1)})B_{\mu \lambda }\psi _\lambda ^{(0)}|\psi _\lambda ^{(0)}=j_\lambda |\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_\mu .$$
(32)
For $`\lambda =\mu `$ this gives the five transport coefficients

$$z_\mu ^{(2)}=\frac{j_\mu |\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_\mu }{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}},$$
(33)
and for $`\lambda \mu `$

$$(z_\mu ^{(1)}z_\lambda ^{(1)})B_{\mu \lambda }\psi _\lambda ^{(0)}|\psi _\lambda ^{(0)}=j_\lambda |\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_\mu ,$$
(34)
from which some of the values of $`B_{\mu \lambda }`$ can be obtained, provided that $`(z_\mu ^{(1)}z_\lambda ^{(1)})`$ is nonzero, and of which the resulting expressions can be found in Appendix A. Note that at this stage three of the zeroth order eigenfunctions and their corresponding currents are still undetermined.

In the space of $`|d`$, $`|s`$, and $`|c_{}`$ one can rewrite Eq. (32) in the form of a new eigenvalue problem in the still undetermined coefficients $`A_{\mu n}`$

$$\underset{n}{}\left(j_m|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_n+z_\mu ^{(2)}a_m|a_n\right)A_{\mu n}=0.$$
(35)
Equating the determinant to zero we find that there are three different eigenvalues and we can obtain the form of the remaining three zeroth order eigenfunctions. Using the fact that there is no degeneracy we know that we can diagonalize the matrix in terms of the proper functions. Hence the off-diagonal matrix elements $`j_\mu |\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_\nu `$ need to vanish, which could already be seen from Eq. (34) since for these modes we have $`z_\mu ^{(1)}=0`$. It can also easily be checked from symmetry considerations that $`|c_{}`$ is one of the modes as was suggested earlier. The solution of the remaining two zeroth order eigenfunctions is straightforward and leads to

$`|\psi _{}^{(0)}`$ $`=|c_{},`$ (36)
$`|\psi _{s_\pm }^{(0)}`$ $`=|s_\pm ={\displaystyle \frac{|s}{\sqrt{s|s}}}+\left(\xi \pm \sqrt{1+\xi ^2}\right){\displaystyle \frac{|d}{\sqrt{d|d}}},`$ (37)
where the prefactors are determined by the requirements

$`s_\pm |s_{}`$ $`=0,`$ (38)
$`j_{s_\pm }|{\displaystyle \frac{1}{\mathrm{\Omega }}}+{\displaystyle \frac{1}{2}}|j_s_{}`$ $`=0,`$ (39)
and we have introduced the short hand notation

$$\xi \frac{\chi 𝒟}{2Q},$$
(40)
with

$`\chi `$ $`{\displaystyle \frac{j_s|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_s}{s|s}},`$ (41)
$`𝒟`$ $`{\displaystyle \frac{j_d|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_d}{d|d}},`$ (42)
$`Q`$ $`{\displaystyle \frac{j_s|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_d}{\sqrt{s|sd|d}}}.`$ (43)
As the form of the last three quantities suggests by comparison with Eq. (33), these quantities are related to transport coefficients. In the low density limit and in the limit of only one specie, one can identify $`\chi `$ as being the thermal diffusivity of the GBL model. In the low density limit $`𝒟`$ will correspond to the self-diffusion coefficient. Their interpretation as transport coefficients is justified in Appendix B. The remaining value $`Q`$ can not be interpreted as a transport value, but rather is some measure of the coupling between different modes.

Using Eq. (33) we can evaluate the five transport coefficients corresponding to the hydrodynamic modes. In the case of the two soundmodes (26) this leads to the sounddamping $`\mathrm{\Gamma }`$

$$z_\pm ^{(2)}=\mathrm{\Gamma }=\frac{j_\pm |\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_\pm }{\psi _\pm ^{(0)}|\psi _\pm ^{(0)}},$$
(44)
while the perpendicular mode (36) gives rise to the viscosity $`\nu `$

$$z_{}^{(2)}=\nu =\frac{j_{}|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_{}}{c_{}|c_{}}.$$
(45)
By writing out the definitions of $`|s_\pm `$ and $`|j_{s_\pm }`$ and using the introduced quantities (40) - (43) where needed, we can rewrite the second order eigenvalues of the two non-propagating modes $`|s_\pm `$ as

$$s_\pm ^{}=\frac{j_{s_\pm }|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_{s_\pm }}{s_\pm |s_\pm }=\frac{1}{2}(\chi +𝒟)\pm Q\sqrt{1+\frac{(\chi 𝒟)^2}{4Q^2}}.$$
(46)
In the low density limit one finds that the two eigenmodes (37) up to a normalization factor will converge to $`|s`$ and $`|d`$ as defined in Eqns. (11) and (12). Consequently one obtains $`s_\pm ^{}\chi ,𝒟`$, which is an illustration of the decoupling of entropy and concentration fluctuations in this limit. In general, however, these two “true” transport coefficients do not seem to correspond with a conventional transport coefficient, but rather they always appear in combination with each other, and it is only in the appropriate limits that they reduce to the thermal diffusivity and diffusion coefficient. This is completely in agreement with the results known for the continuous theory Boon and Yip (1980), where the coupling of fluctuations in concentration and entropy results in the same effect.

In the present model, however, there is some ambiguity in the choice made for the basic invariants. To arrive at this last formula it was only necessary to assume that the invariants $`|s`$ and $`|d`$, are mutually orthogonal, and span a 2-dimensional subspace in the null-space which is perpendicular to $`|c_{\mathrm{}}`$, $`|c_{}`$, and $`|p`$. A good choice necessarily leads to the proper modes in the fully known limit of one specie. Although this puts some limitations on the possible choices it will not uniquely fix the basis. The present choice for $`|s`$ in (11) is however the most natural extension of the conventional formulation Grosfils et al. (1993) and we will refer to $`\chi `$ as being the generalized thermal diffusivity.

This problem is more easily identified in the diffusion-like mode $`|d`$. In the present formulation (12), it contains in general a contribution proportional to $`|p`$, rather than being a weighted difference of the red and blue densities only as found in the colored GBL model Blaak and Dubbeldam (2001a). This already suggests that this vector is not the proper generalization of the self-diffusion mode. This is consistent with the continuous theory Boon and Yip (1980), where the analogue of the $`s_\pm ^{}`$ also depend on more than the thermal diffusivity and diffusion only. We will address this subject again in a later section.

Similar to the case of the GBL model we can introduce $`\gamma `$, the ratio of specific heats by

$$\gamma =1+\frac{s|s}{p|p}=\frac{p|p\rho |\rho }{p|\rho ^2}=\frac{c_s^2}{c_T^2},$$
(47)
with $`c_s^2=p|p/p|\rho `$ the adiabatic speed of sound, $`c_T^2=p|\rho /\rho |\rho `$ the isothermal speed of sound. Two additional useful definitions are $`\tau _{xy}=c_{\mathrm{}}c_{}`$ and $`\tau _{xx}=\frac{1}{2}(c_{\mathrm{}}^2c_{}^2)`$. These allow us to rewrite the currents related to the sound modes as $`|j_\pm =|j_s\pm c_s|\tau _{xx}`$. Consequently this leads to

$$\mathrm{\Gamma }=\frac{j_s|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_s+c_s^2\tau _{xx}|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|\tau _{xx}\pm 2c_s^2j_s|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|\tau _{xx}}{p|p+c_s^2c_{\mathrm{}}|c_{\mathrm{}}\pm 2c_sp|c_{\mathrm{}}}.$$
(48)
From symmetry considerations it follows that $`p|c_{\mathrm{}}=0`$ and $`j_s|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|\tau _{xx}=0`$, and since $`\tau _{xx}|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|\tau _{xx}=\tau _{xy}|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|\tau _{xy}`$ because of the isotropy of the lattice, we obtain the following relation for the main transport coefficients

$$\mathrm{\Gamma }=\frac{1}{2}\left(\nu +(\gamma 1)\chi \right).$$
(49)
An example of the full wave-vector dependent real eigenvalue spectrum is shown in Fig. 2, which also confirms the absence of spurious invariants. Of all eigenvalues only the five related to the hydrodynamic modes go to zero in the limit $`𝒌0`$. It is with these modes that the binary-mixture responds to deviations from thermal equilibrium

$`z_\pm (𝒌)`$ $`=\pm ıc_s𝒌\mathrm{\Gamma }𝒌^2,`$ (50)
$`z_{}(𝒌)`$ $`=\nu 𝒌^2,`$ (51)
$`z_{s_\pm }(𝒌)`$ $`=s_\pm ^{}𝒌^2.`$ (52)
The first two eigenvalues describe sound propagation in the two opposite directions parallel to $`𝒌`$ with $`c_s`$ the adiabatic sound speed, the third eigenvalue describes the shear mode, and the last two eigenvalues represent purely diffusive, non-propagating processes. In this hydrodynamic regime characterized by $`k\lambda 1`$, where $`\lambda `$ is the mean free path length, one can exploit the fact that the real component of the eigenvalues of the hydrodynamic modes is much smaller than that of the kinetic modes. We will use this in the next section in order to obtain the Landau-Placzek formula (59) for the power spectrum.

For larger wavevectors (here roughly $`0.5k1.5`$) the relations (50)-(52) start to deviate from the true values. This is the generalized hydrodynamic regime ($`k\lambda 1`$), and the transport coefficients become $`k`$ dependent. The hydrodynamic modes are, however, still smaller than the kinetic modes yielding a reasonable accurate Landau-Placzek formula. In the kinetic regime ($`k\lambda 1`$) the hydrodynamic modes and kinetic modes become of the same order of magnitude. The distinction between fast and slow modes can not be made and all modes contribute to the power spectrum.

In Fig. 3 the diffusive transport properties are shown at fixed reduced temperature $`\theta =\mathrm{exp}(\frac{1}{2}\beta )`$. The modes $`s_\pm ^{}`$ are the modes observed in the eigenvalue spectrum and are combinations of $`\chi `$, $`𝒟`$, and $`Q`$. In the low density limit these two modes converge to the thermal diffusivity and diffusion, i.e. $`s_+^{}\chi `$, $`s_{}^{}𝒟`$, and is caused by the decoupling of the fluctuations in the concentrations and entropy. In our model this means that the ratio $`Q/(\chi 𝒟)`$ vanishes. The value of $`Q`$, however, will in general remain small but finite due to the divergencies of the transport coefficients in the low density limit of LGA. The same is observed in the high density limit, which is merely an illustration of the duality of the LGA-model if one interchanges particles and holes. In both cases this is a direct consequence of the fact that the fluctuations in the occupation numbers become linear in the average occupation densities, thus leading to equipartition.

As an additional remark we like to mention that the value of $`Q`$ can be both positive and negative. Therefore the correct curves of $`s_\pm ^{}`$ need not be continuous as a function of the density, but can contain discontinuities located at points where $`Q`$ changes sign. Consequently, depending on the system parameters, the role of $`s_+^{}`$ and $`s_{}^{}`$ is interchanged with respect to $`\chi `$ and $`𝒟`$. This is merely due to the choice in convention we have used in Eq. (37) and of no physical importance.

For some intermediate values of the density one also finds that the values of $`s_\pm ^{}`$ and $`\chi ,𝒟`$ coincide. This is, however, not caused by decoupling of the fluctuations. In these cases there is no effective equipartition, but one obtains $`Q=0`$ as a consequence of the cancellation of terms. Moreover, the location of these points depends in a non-trivial way on the system parameters.

The diffusive transport properties versus the relative concentration of both species at fixed density and reduced temperature $`\theta `$ is shown in Fig. 4. In the low density limit we have $`s_\pm ^{}\chi ,𝒟`$ for all relative concentrations. This is not generally true as can be seen from the higher density figure. It does, however, reveal that the decoupling also occurs in general for the single specie limits $`P_r0`$ and $`P_r1`$. This is not surprising, as in those limits the model effectively reduces to a normal GBL model with only a single diffusive mode related to the thermal diffusivity. This is actually also the case in the high density limit of Fig. 3, because there one has a situation in which the red lattice is almost completely filled and therefore an effective blue system that remains.

Finally, decoupling can also be observed in the case $`P_r=P_b`$. In fact this is a rather special limit and could be analyzed completely in a manner similar to the one used for the colored GBL model Blaak and Dubbeldam (2001a), because based on the symmetry of the problem one can decompose the linearized Boltzmann operator in two type of contributions, i.e. $`|R+|B`$ and $`|R|B`$.

For $`P_r=0`$ the diffusion mode vanishes, since we then recover the original GBL model, but for $`P_r0`$ the diffusion mode remains finite. This case is illustrated in Fig. 5. Here, in the low density limit the diffusion mode becomes equal to the self-diffusion of the GBL model (and hence also equal to the diffusion mode in the GBL color mixture), but for higher density they deviate. A behavior that finds its origin in the fact that the color mixture does not allow more than one particle in a single velocity channel and, although the number of particles of the second specie gets very small, they still have a large impact on the diffusion.

## IV Landau-Placzek theory

### IV.1 Full spectrum

In the hydrodynamic regime of small wavevectors ($`k0`$) and the frequency $`\omega `$ being linear in $`k`$, the spectral density $`S(𝒌,\omega )`$ can be expanded in powers of $`ık`$. In this long-wavelength, small frequency limit the hydrodynamic modes are well separated from the kinetic modes, which can be neglected due to their exponential decay. By keeping only terms up to $`𝑶(k^2)`$ ($`\omega /k`$ is kept constant for consistency), one obtains the Landau-Placzek approximation.

It has been shown that the dynamic structure factor can be evaluated by Grosfils et al. (1993)

$$\frac{S(𝒌,\omega )}{S(𝒌)}=2\text{Re}\frac{\rho |(e^{ı\omega +ı𝒌𝒄}1\mathrm{\Omega })^1+\frac{1}{2}|\rho }{\rho |\rho }.$$
(53)
The closely related spectral function $`F(𝒌,\omega )`$ can be written as

$`F(𝒌,\omega )\rho |(e^{ı\omega +ı𝒌𝒄}1\mathrm{\Omega })^1+{\displaystyle \frac{1}{2}}|\rho =2{\displaystyle \underset{\mu }{}}\text{Re}𝒩_\mu 𝒟_\mu ,`$ (54)
where

$`𝒟_\mu `$ $`={\displaystyle \frac{1}{e^{ıwz_\mu }1}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{z_\mu ^{(2)}ı(wz_\mu ^{(1)}k)}{(z_\mu ^{(2)}k^2)^2+(wz_\mu ^{(1)}k)^2}},`$ (55)
$`𝒩_\mu `$ $`=\rho |\psi _\mu \varphi _\mu |\rho .`$ (56)
The coefficients $`𝒩_\mu `$ are evaluated in Appendix A for small $`𝒌`$ and yield

$`𝒩_\sigma `$ $`={\displaystyle \frac{\rho |\rho }{2\gamma }}\left(1+{\displaystyle \frac{ı\sigma k}{c_s}}\left[\mathrm{\Gamma }+(\gamma 1)\chi \right]\right),`$ (57)
$`𝒩_{s_\pm }`$ $`={\displaystyle \frac{(\gamma 1)\rho |\rho }{2\gamma }}\left(1\pm {\displaystyle \frac{\chi 𝒟}{s_+^{}s_{}^{}}}\right).`$ (58)
From symmetry considerations one can conclude that the shear mode will not contribute to the spectrum and one finds $`𝒩_{}=0`$. Combining these results with the expressions (54) and (55), we finally obtain the Landau-Placzek formula, describing the power-spectrum in the hydrodynamical domain

$`{\displaystyle \frac{S(𝒌,\omega )}{S(𝒌)}}={\displaystyle \frac{\gamma 1}{\gamma }}\left[\left(1+{\displaystyle \frac{\chi 𝒟}{s_+^{}s_{}^{}}}\right){\displaystyle \frac{s_+^{}k^2}{\omega ^2+(s_+^{}k^2)^2}}+\left(1{\displaystyle \frac{\chi 𝒟}{s_+^{}s_{}^{}}}\right){\displaystyle \frac{s_{}^{}k^2}{\omega ^2+(s_{}^{}k^2)^2}}\right]`$
$`+{\displaystyle \frac{1}{\gamma }}\left[{\displaystyle \frac{\mathrm{\Gamma }k^2}{(\omega +c_sk)^2+(\mathrm{\Gamma }k^2)^2}}+{\displaystyle \frac{\mathrm{\Gamma }k^2}{(\omega c_sk)^2+(\mathrm{\Gamma }k^2)^2}}\right]`$
$`+{\displaystyle \frac{1}{\gamma }}\left[\mathrm{\Gamma }+(\gamma 1)\chi \right]{\displaystyle \frac{k}{c_s}}\left[{\displaystyle \frac{\omega +c_sk}{(\omega +c_sk)^2+(\mathrm{\Gamma }k^2)^2}}{\displaystyle \frac{\omega c_sk}{(\omega c_sk)^2+(\mathrm{\Gamma }k^2)^2}}\right].`$ (59)
The spectrum contains an unshifted central peak that is formed by two Lorentzians due to the two processes related to the non-propagating modes $`s_\pm ^{}`$. The two propagating modes lead to the presence of the two shifted Brillouin lines. Their width at half-height can be used as a measurement of $`\mathrm{\Gamma }𝒌^2`$, the position of the peaks can be used a measurement of $`\pm c_s𝒌`$. The last two terms in Eq. (59) give an asymmetric correction to the Brillouin peaks and induce a slight pulling of the peaks toward the central peak.

The symmetry of the different contributions is such that the ratio of the integrated contributions of the central peak and the Brillouin components is constant and given by

$$\frac{_\sigma 𝑑\omega Re(𝒩_{s_\sigma }𝒟_{s_\sigma })}{_\sigma 𝑑\omega Re(𝒩_\sigma 𝒟_\sigma )}=\gamma 1.$$
(60)
In Fig. 6 the Landau-Placzek formula is compared with the full Boltzmann spectrum for two different wavevectors. For the smallest wavevector inside the hydrodynamic regime they coincide, while for the larger they differ considerably and indicates that this wavevector is outside the hydrodynamic, and in fact inside the kinetic regime.

In general the contributions forming the central peak cannot easily be separated. Even in the case they differ sufficiently in order to fit the central peak with two Lorentzians, one only obtains information on the values of $`s_\pm ^{}`$, which is not enough to determine the more interesting values of the thermal diffusivity and diffusion. An example of this is shown in Fig. 7 where we decomposed the central peak in the two Lorentzians using the theoretical expressions. For comparison we included the least square fit according to Eq.(59) to the full range of the spectrum, including the Brillouin peaks, where we only used a single Lorentz for the central peak.

For this reason it is important to consider the limits where decoupling occurs, since in those cases the the relevant transport values can be obtained from the spectra. Here we have however an interesting difference with respect to light scattering experiments Boon and Yip (1980). Whereas in that case the sensitivity of the dielectric fluctuations with respect to density is much larger than with respect to temperature, the amplitude factors will differ orders of magnitude. This results in the observation that the diffusion component will usually dominate the spectrum. For LGA this does not apply, and by using the decoupling limit $`s_\pm ^{}\chi ,𝒟`$ on the first term off the Landau-Placzek formula (59), one immediately obtains that the central component of the spectrum is given by a single peak characterized by the thermal diffusivity

$$\frac{S^{\text{cen}}(𝒌,\omega )}{S(𝒌)}=\frac{\gamma 1}{\gamma }\frac{2\chi k^2}{\omega ^2+(\chi k^2)^2}.$$
(61)
One could in principle use this also in some of the intermediate cases seen in Fig. 3. The location of these points in terms of density and relative fractions, however, depends in a non-trivial way on the system parameters. Moreover, the identification with the thermal diffusivity and diffusion can only be made if one interpret $`\chi `$ and $`𝒟`$ as generalizations of these quantities.

### IV.2 Red spectrum

The red dynamic structure factor can be evaluated by Hanon and Boon (1997); Blaak and Dubbeldam (2001a)

$$\frac{S^{\text{red}}(𝒌,\omega )}{S(𝒌)}=2\text{Re}\frac{R|(e^{ı\omega +ı𝒌𝒄}1\mathrm{\Omega })^1+\frac{1}{2}|R}{R|R}.$$
(62)
It is possible to follow the same route for the red dynamic structure factor as for the full spectrum. However, since $`\rho |d=0`$, but $`R|d0`$, it is more convenient to use a different set of basic invariants

$$|a_n=\{|c_x,|c_y,|p=|\frac{1}{2}c^2,|s_r=|p\frac{p|p}{p|R}|R,|d_r\},$$
(63)
where $`|d_r`$ is constructed to be perpendicular to the other conserved quantities and defined by

$$|d_r=|p\left(\frac{p|R}{R|R}|R\frac{p|R^2p|pR|R}{R|Rp|B}|B\right).$$
(64)
Note that by construction it is perpendicular to the red density, i.e. $`R|d=0`$. Analogous to the case of the normal density we introduce some transport-like coefficients related to this basis

$`\chi _\text{r}`$ $`={\displaystyle \frac{s_r|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|s_r}{s_r|s_r}},`$ (65)
$`𝒟_\text{r}`$ $`={\displaystyle \frac{d_r|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|d_r}{d_r|d_r}},`$ (66)
$`\gamma _\text{r}`$ $`=1+{\displaystyle \frac{s_r|s_r}{p|p}}.`$ (67)
Obviously, the sound damping $`\mathrm{\Gamma }`$, the modes $`s_\pm ^{}`$, and the speed of sound $`c_s`$ are all unchanged, since these are true transport values and independent of a chosen set of basis functions. Also in this case we can derive a Landau-Placzek formula, describing the red power-spectrum in the hydrodynamical domain

$`{\displaystyle \frac{S^{\text{red}}(𝒌,\omega )}{S(𝒌)}}={\displaystyle \frac{\gamma _\text{r}1}{\gamma _\text{r}}}\left[\left(1+{\displaystyle \frac{\chi _\text{r}𝒟_\text{r}}{s_+^{}s_{}^{}}}\right){\displaystyle \frac{s_+^{}k^2}{\omega ^2+(s_+^{}k^2)^2}}+\left(1{\displaystyle \frac{\chi _\text{r}𝒟_\text{r}}{s_+^{}s_{}^{}}}\right){\displaystyle \frac{s_{}^{}k^2}{\omega ^2+(s_{}^{}k^2)^2}}\right]`$
$`+{\displaystyle \frac{1}{\gamma _\text{r}}}\left[{\displaystyle \frac{\mathrm{\Gamma }k^2}{(\omega +c_sk)^2+(\mathrm{\Gamma }k^2)^2}}+{\displaystyle \frac{\mathrm{\Gamma }k^2}{(\omega c_sk)^2+(\mathrm{\Gamma }k^2)^2}}\right]`$
$`+{\displaystyle \frac{1}{\gamma _\text{r}}}\left[(\nu \mathrm{\Gamma })+2(\gamma _\text{r}1)\chi _\text{r}\right]{\displaystyle \frac{k}{c_s}}\left[{\displaystyle \frac{\omega +c_sk}{(\omega +c_sk)^2+(\mathrm{\Gamma }k^2)^2}}{\displaystyle \frac{\omega c_sk}{(\omega c_sk)^2+(\mathrm{\Gamma }k^2)^2}}\right].`$ (68)
In principle the quantities $`\chi _\text{r}`$, $`𝒟_\text{r}`$, and $`\gamma _\text{r}`$ can be expressed in the normal transport values. However, these relations lead to more complex expressions and are for that reason omitted here. In Fig. 6 this formula is compared with the full Boltzmann red spectrum for a wave vector in the hydrodynamic regime, yielding a satisfactory approximation, and one in the kinetic regime with a large discrepancy.

### IV.3 Diffusion spectrum

Finally we like to consider the spectra of the diffusive processes only. Due to the nature of LGA, however, it is not obvious how to define diffusion properly. On the one hand we have the requirement that a diffusive process is not propagating, on the other hand one expects the diffusion only to be related to the densities of both components in the mixture. As we will show in this section, it turns out that these two views are not completely compatible.

In the colored GBL model Blaak and Dubbeldam (2001a) and in the continuous case McLennan (1988) a diffusion spectrum can be obtained by considering fluctuations in the normalized difference of the red and blue density

$$\rho _{\text{diff}}=\frac{\rho _{\text{red}}}{P_r}\frac{\rho _{\text{blue}}}{P_b}.$$
(69)
The proper translation in terms of invariants is given by

$$|\text{diff}=\frac{|R}{P_r}\frac{|B}{P_b}.$$
(70)
However, a spectrum based on this vector, does not lead to purely diffusive peaks only, but also includes parts of the propagating modes. This can easily be checked since the vector is in general not orthogonal to the two soundmodes due to the absence of equipartition.

A simple solution is to subtract the propagating part by adding the appropriate term proportional to $`|p`$

$$|\text{diff}=\frac{|R}{P_r}\frac{|B}{P_b}\left(\frac{p|R}{P_r}\frac{p|B}{P_b}\right)\frac{|p}{p|p}.$$
(71)
In a binary athermal mixture Ernst and Dufty (1990), however, it was suggested to use

$$|\text{diff}=\frac{|R}{R|R}\frac{|B}{B|B}.$$
(72)
Unfortunately this choice also leads to a propagating mode in the thermal case. One could again subtract the propagating part, but the result would be different from (71). The origin of this problem is that we have two different diffusive modes and any linear combination would lead to a diffusive spectrum. However, due to the lack of equipartition the different choices of fluctuations one wishes to consider do not coincide, not even after the propagating part is eliminated.

Naively one would expect a combination of the red and blue component only and in addition it should be perpendicular to the propagating modes. This leads to the following generalization of the athermal result (72), which in the case of an athermal model would coincide

$$|\text{diff}=\frac{|R}{p|R}\frac{|B}{p|B}.$$
(73)
Obviously there is some freedom here in order to choose the generalization of the diffusion. The natural extension would be one that that satisfies the appropriate limits. However, in the low density limit and in the limit $`P_r1`$, they all converge to the same value equal to the one found for the single component GBL model. The special case of an equal density for red and blue particles is not helpful either. This limit can be completely analyzed in a manner analogous to what is done for colored GBL model Blaak and Dubbeldam (2001a). For reasons of symmetry this will cause none of the definitions to have a propagating character and the resulting diffusions would coincide.

In Appendix B we derive the proper definition (100) for the transport values related to these quantities. Since none of them are proportional to a single eigenmode the usual formulation (33) is no longer valid. The results can be found in Fig. 8, where the diffusions obtained for the various diff’s are shown as function of the density. In general the five different diff’s described above all lead to different values for the corresponding transport value, although depending on the choice of system parameters this difference might be marginal. Also compare with Fig. 3 to see the difference with respect to $`𝒟`$.

Although the choice (73) is the most natural generalization, the diffusion is only properly identified in the limiting cases. In general some arbitrariness remains in a binary thermal lattice gas.

The corresponding diffusion spectrum is obtained by

$$\frac{S(𝒌,\omega )}{S(𝒌)}=2\text{Re}\frac{\text{diff}|(e^{ı\omega +ı𝒌𝒄}1\mathrm{\Omega })^1+\frac{1}{2}|\text{diff}}{\text{diff}|\text{diff}}.$$
(74)
and shown in Fig. 9 for two diff’s: the $`|\text{diff}`$ defined in Eq. (72) based on an athermal model and the correct adjustment for the thermal model Eq. (73). For comparison we also included the curve for which the propagating part in Eq. (72) is subtracted as was described above. The first diffusion spectrum has a Brillouin-like pair of peaks which is a manifestation of the propagating part of $`|\text{diff}`$ according to (72), while the other two are different superpositions of two Lorentzians characterized by $`s_{\pm }^{}{}_{}{}^{(0)}`$.

## V Simulation results

We have verified our results with simulation results. Some notes on the implementation are in order. The model we consider is a true $`38`$-bits model, i.e. the collision operator can act on $`2^{38}`$ different states. In a previous article we reported the simulation results of a colored GBL model Blaak and Dubbeldam (2001a). The number of states in that model was $`3^{19}`$, since no red and blue particle could exist with the same velocity. The simulations could, however, be performed by realizing that the collisions could be separated in a normal GBL collision followed by a redistribution of the colors over the different particles present. These two separated processes can easily be performed by any computer.

In the present model this approach can not be used. One could of course make use of the symmetries present in the model. The hexagonal lattice leads to six rotations and two reflections. In addition we can make use of the interchange of red and blue particles, and the interchange of particles and holes. At best this would lead to a reduction with a factor of 48 on the total number of different states, but in practice this factor is less because a large fraction of the states is invariant with respect to some of these symmetries. This is however still too large in order to be applied in any but the largest supercomputers available at present. An alternative, but rather inefficient scheme, would be to store only part of the collision table and generate other collisions on the fly.

Binary mixtures, however, do allow for a more convenient solution, which is less efficient than storing the complete collision table, but still has a relative good performance (about a factor 2-5 slower than the colored GBL model). Rather than storing a collision table based on the states, we make one based on the different classes. In the GBL model there are $`29926`$ different classes $`𝒞=(M,𝑷,E)`$, characterized by the total mass, momentum, and energy. Since the the binary mixture is the combination of two GBL models there are $`29926^2`$ different combinations $`(𝒞_R,𝒞_B)`$ of a red and blue class which form a total of $`4478568`$ classes $`\mathrm{\Gamma }=(M_r,M_b,𝑷,E)`$. If we combine this with the 48 symmetry operations we get a working algorithm that already can be used on a computer with 256 MB of memory.

An arbitrary input state is now analyzed in order to determine to which class $`\mathrm{\Gamma }`$ it belongs. From its characterization one can easily find which transformation is needed in order put it in a form with $`M_RM_B`$, $`M_R+M_B19`$, and $`0P_yP_x/\sqrt{3}`$. The last inequality confines the total momentum to an angle of $`\pi /6`$. For each $`\mathrm{\Gamma }`$ in this limited set of classes only, all the classes $`𝒞_R`$, that with the appropriate class $`𝒞_B`$ can give rise to it, have been stored, including the number of states in each $`(𝒞_R,𝒞_B)`$ combination. From this we generate with the proper weight an “outgoing” class combination to which the inverse transformation is applied. Finally we only need to determine a random state in $`𝒞_R`$ and $`𝒞_B`$, which is just a GBL-like process.

Simulation results along with the Boltzmann approximation and the Landau-Placzek formula are shown in Fig. 6 for both the normal and red spectrum. In the Figs. 6a and 6b the spectra at low $`𝒌`$ value completely overlap and the Landau-Placzek formula describes the power spectrum very well. If we increase the wavevector $`𝒌`$, and/or lower the density, the hydrodynamical regime is left and deviations start to appear. In that generalized regime, the transport properties are $`𝒌`$-dependent. Figs. 6c and 6d show a spectrum even further away from the hydrodynamic regime, in the kinetic regime. Several kinetic modes invade the spectrum, and a parameterization into Lorentzians has lost all physical meaning. The full spectrum based on the Boltzmann approximation, however, still leads to a very good description, supporting the molecular chaos assumption.

## VI Discussion

We have constructed a thermal binary lattice gas mixture. The model is characterized by cross-effects between energy transport and diffusion. The Landau-Placzek formula derived in the the low wavevector, low frequency domain gives an excellent description. For larger wavevectors it will fail, but the spectra can still accurately be described by the Boltzmann approximation.

The Landau-Placzek formula for a regular binary thermal mixture is quite similar in structure as the one for the continuous case. The main feature is a central peak formed by two Lorentzians due to the coupled entropy-concentration fluctuations. In the limits of low density, single specie, their dual interpretations, and equal red and blue density it is possible for the central peak of the spectrum to be decomposed into two Lorentzians with linewidths given by $`\chi 𝒌^2`$ and $`𝒟𝒌^2`$, otherwise the linewidths depend on both these values.

In the low density limit the true modes $`s_\pm ^{}`$ converge to $`\chi `$ and $`𝒟`$. Moreover, they can be identified with the thermal diffusivity and mass-diffusion coefficient of the GBL model and coincide with the low density limit of the continuous binary mixture. In contrast with continuous theory, however, it is the thermal diffusivity that dominates the central part of the spectrum, rather than the diffusion. The interpretation of $`\chi `$ as the thermal diffusivity can be extended to the general situation. For $`𝒟`$ this not true, since the related spectrum will contain Brillouin peaks, which is a not a purely diffusive property. Although it is possible to correct for this, there remains some ambiguity in the choice for the generalized diffusion coefficient.

The analysis and results presented here are in general true for binary thermal lattice gasses and not restricted to this particular model only. It also reduces automatically to the proper formulation for an athermal model, in which case some of the ambiguities are removed.

## Acknowledgments

We would like to thank H. Bussemaker and D. Frenkel for helpful discussions. R.B. acknowledges the financial support of the EU through the Marie Curie Individual Fellowship Program (contract no. HPMF-CT-1999-00100).

## Appendix A

To calculate terms of the type $`𝒩_\mu `$ given in Eq. (56) we express the left eigenvectors in terms of right eigenvectors, and expand $`𝒩_\mu `$ to linear order in $`k`$

$`𝒩_\mu =`$ $`{\displaystyle \frac{\rho |\psi _\mu ^{(0)}^2}{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}}}(1+2ık[{\displaystyle \frac{\rho |\psi _\mu ^{(1)}}{\rho |\psi _\mu ^{(0)}}}{\displaystyle \frac{\psi _\mu ^{(0)}|\psi _\mu ^{(1)}}{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}}}]+`$ (75)
$`ık[{\displaystyle \frac{\rho |c_{\mathrm{}}\psi _\mu ^{(0)}}{\rho \psi _\mu ^{(0)}}}{\displaystyle \frac{\psi _\mu ^{(0)}|c_{\mathrm{}}\psi _\mu ^{(0)}}{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}}}]+𝒪(k^2)).`$
As can be seen from the second term on the righthand side, partial knowledge of the $`|\psi _\mu ^{(1)}`$ is required. From the formal solution (30) of the first order equation and the orthogonality relation (29) it follows that we only need to evaluate the coefficients $`B_{\mu \nu }`$. From the normalization (19) we already know that $`B_{\mu \mu }=0`$. In addition we like to mention that from symmetry considerations one can conclude that the mode related to the viscosity will not contribute. Hence the values of $`B_\mu `$ and $`B_\mu `$ need not be determined.

In order to facilitate the calculations we first list some simple relations

$$\psi _\sigma ^{(0)}|\psi _\sigma ^{(0)}=2p|p,$$
(76)
$$\rho |\psi _\sigma ^{(0)}=\rho |p,$$
(77)
$$\rho |s_\pm =\frac{(1\gamma )}{\sqrt{s|s}}\rho |p,$$
(78)
$$s_\pm |s_\pm =2\sqrt{1+\xi ^2}\left(\sqrt{1+\xi ^2}\xi \right).$$
(79)
Most coefficients can be evaluated from the relation (34). In the case of the soundmodes this leads to

$$B_{\sigma ,\sigma }=\frac{j_\sigma |\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_\sigma }{(z_\sigma ^{(1)}z_\sigma ^{(1)})\psi _\sigma ^{(0)}|\psi _\sigma ^{(0)}}$$
(80)
and with the use of $`|j_\pm =|j_s\pm c_s|\tau _{xx}`$ this can be rewritten in terms of the basic transport coefficients

$$B_{\sigma ,\sigma }=\sigma \frac{\nu (\gamma 1)\chi }{4c_s}=\frac{\sigma (\nu \mathrm{\Gamma })}{2c_s}.$$
(81)
The next coefficient we need to evaluate is

$$B_{s_\pm ,\sigma }=\frac{j_\sigma |\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|j_{s_\pm }}{(z_{s_\pm }^{(1)}z_\sigma ^{(1)})\psi _\sigma ^{(0)}|\psi _\sigma ^{(0)}},$$
(82)
and here we can use the same relation to rewrite the current of the soundmodes. For the other current the relation (37) can be used. Realizing that the terms containing $`|\tau _{xx}`$ will vanish, the remaining terms can be manipulated to yield

$$B_{s_\pm ,\sigma }=(\gamma 1)\frac{\sigma s_\pm ^{}}{2c_s\sqrt{s|s}}.$$
(83)
From Eq. (34) we can now also see that

$$B_{\sigma ,s_\pm }=\frac{\psi _\sigma ^{(0)}|\psi _\sigma ^{(0)}}{s_\pm |s_\pm }B_{s_\pm ,\sigma }=\left(1\pm \frac{\xi }{\sqrt{1+\xi ^2}}\right)\frac{\sigma s_\pm ^{}\sqrt{s|s}}{2c_s}.$$
(84)
The last two coefficients which we will need are $`B_{s_\pm ,s_{}}`$. However, since $`z_{s_\pm }^{(1)}=0`$ they can not be obtained from (34). In order for them to be determined we need to make use of the third order equation of the eigenvalue problem (13)

$`\mathrm{\Omega }|\psi _\mu ^{(3)}`$ $`=`$ $`(c_{\mathrm{}}+z_\mu ^{(1)})|\psi _\mu ^{(2)}+\left[z_\mu ^{(2)}+{\displaystyle \frac{1}{2}}(c_{\mathrm{}}+z_\mu ^{(1)})^2\right]|\psi _\mu ^{(1)}+`$ (85)
$`\left[z_\mu ^{(3)}+z_\mu ^{(2)}(c_{\mathrm{}}+z_\mu ^{(1)})+{\displaystyle \frac{1}{6}}(c_{\mathrm{}}+z_\mu ^{(1)})^3\right]|\psi _\mu ^{(0)}.`$
As we only want to determine the values of $`B_{s_\pm ,s_{}}`$, we do not attempt to solve the complete third order equation but restrict ourselves to the two equations from which they can be obtained. Note that we have $`z_{s_\pm }^{(1)}=0`$ and hence $`|j_{s_\pm }=|c_{\mathrm{}}s_\pm `$.

Substituting these results in the third order equation, replacing $`\mu `$ by $`s_\pm `$, and multiplying on the left with the appropriated term $`s_{}|`$ we find

$$0=s_{}|c_{\mathrm{}}|\psi _{s_\pm }^{(2)}+s_{}|s_\pm ^{}+\frac{1}{2}c_{\mathrm{}}^2|\psi _{s_\pm }^{(1)}+s_{}|z_{s_\pm }^{(3)}+s_\pm ^{}c_{\mathrm{}}+\frac{1}{6}c_{\mathrm{}}^3|s_\pm .$$
(86)
The last term will disappear due to the odd power of $`c_{\mathrm{}}`$ and because $`s_{}|s_\pm =0`$. In order to proceed we not only need the solution of the first order equation (30) but also solution of the second order equation. Fortunately the later does not have to be computed completely but the formal solution will suffice

$`|\psi _\mu ^{(2)}={\displaystyle \frac{1}{\mathrm{\Omega }}}(c_{\mathrm{}}+z_\mu ^{(1)})|\psi _\mu ^{(1)}+{\displaystyle \frac{1}{\mathrm{\Omega }}}\left(z_\mu ^{(2)}+{\displaystyle \frac{1}{2}}(c_{\mathrm{}}+z_\mu ^{(1)})^2\right)|\psi _\mu ^{(0)}+`$
$`{\displaystyle \underset{\nu }{}}C_{\mu \nu }|\psi _\nu ^{(0)},`$ (87)
where the $`C_{\mu \nu }`$ are some unknown coefficients. Substitution of both solutions in Eq. (86) gives

$`0=s_{}|c_{\mathrm{}}{\displaystyle \frac{1}{\mathrm{\Omega }}}c_{\mathrm{}}|\psi _{s_\pm }^{(1)}+s_{}|c_{\mathrm{}}{\displaystyle \frac{1}{\mathrm{\Omega }}}\left(s_\pm ^{}+{\displaystyle \frac{1}{2}}c_{\mathrm{}}^2\right)|s_\pm +{\displaystyle \underset{\nu }{}}C_{s_\pm \nu }j_s_{}|\psi _\nu ^{(0)}+`$
$`s_{}|\left(s_\pm ^{}+{\displaystyle \frac{1}{2}}c_{\mathrm{}}^2\right){\displaystyle \frac{1}{\mathrm{\Omega }}}|c_{\mathrm{}}s_\pm +{\displaystyle \underset{\nu }{}}B_{s_\pm \nu }s_{}|\left(s_\pm ^{}+{\displaystyle \frac{1}{2}}c_{\mathrm{}}^2\right)|\psi _\nu ^{(0)}.`$ (88)
The second and fourth term on the right are zero because of the odd power in $`c_{\mathrm{}}`$, the third term is zero because of the orthogonality relations (28). Another substitution of the first order solution leads to

$`0=`$ $`s_{}|c_{\mathrm{}}{\displaystyle \frac{1}{\mathrm{\Omega }}}c_{\mathrm{}}{\displaystyle \frac{1}{\mathrm{\Omega }}}c_{\mathrm{}}|s_\pm +{\displaystyle \underset{\nu }{}}B_{s_\pm \nu }s_{}|c_{\mathrm{}}{\displaystyle \frac{1}{\mathrm{\Omega }}}c_{\mathrm{}}|\psi _\nu ^{(0)}+`$ (89)
$`{\displaystyle \underset{\nu }{}}B_{s_\pm \nu }s_{}|\left(s_\pm ^{}+{\displaystyle \frac{1}{2}}c_{\mathrm{}}^2\right)|\psi _\nu ^{(0)}.`$
The first term is again zero because of the odd power in $`c_{\mathrm{}}`$. Writing out the sum, using the orthogonality relations, and realizing that the shear mode dependent term cancels due to symmetries this leads to

$`B_{s_\pm ,s_{}}s_\pm ^{}s_{}|s_{}+B_{s_\pm ,s_{}}j_s_{}|{\displaystyle \frac{1}{\mathrm{\Omega }}}+{\displaystyle \frac{1}{2}}|j_s_{}+`$
$`B_{s_\pm ,+}j_s_{}|{\displaystyle \frac{1}{\mathrm{\Omega }}}+{\displaystyle \frac{1}{2}}|c_{\mathrm{}}\psi _+^{(0)}+B_{s_\pm ,}j_s_{}|{\displaystyle \frac{1}{\mathrm{\Omega }}}+{\displaystyle \frac{1}{2}}|c_{\mathrm{}}\psi _{}^{(0)}=0.`$ (90)
From symmetry considerations it follows that the terms with $`|\psi _\pm ^{(0)}`$ only contribute via $`|p`$ and by using the definition of the transport values on the second term we obtain

$$B_{s_\pm ,s_{}}(s_\pm ^{}s_{}^{})s_{}|s_{}+(B_{s_\pm ,+}+B_{s_\pm ,})j_s_{}|\frac{1}{\mathrm{\Omega }}+\frac{1}{2}|c_{\mathrm{}}p=0.$$
(91)
Using the expressions for $`B_{s_\pm ,\sigma }`$ and the fact that $`s_+^{}`$ and $`s_{}^{}`$ are different, it follows that $`B_{s_\pm ,s_{}}=0`$.

The evaluation of the $`𝒩`$ now becomes straightforward. The last term at the righthand site of Eq. (75) will vanish because of symmetry reasons. From the second term it can be observed that even if we would not have chosen $`B_{\mu \mu }=0`$ as a normalization, this coefficient would not contribute to the spectrum. For the viscosity we immediately obtain $`𝒩_{}=0`$. In the case of $`\mu =\sigma `$ this results in

$`𝒩_\sigma `$ $`=`$ $`{\displaystyle \frac{\rho |p^2}{2p|p}}\left(1+2ık\left[B_{\sigma ,\sigma }+B_{\sigma ,s_+}{\displaystyle \frac{\rho |s_+}{\rho |\psi _\sigma ^{(0)}}}+B_{\sigma ,s_{}}{\displaystyle \frac{\rho |s_{}}{\rho |\psi _\sigma ^{(0)}}}\right]\right),`$ (92)
$`=`$ $`{\displaystyle \frac{\rho |\rho }{2\gamma }}\left(1+{\displaystyle \frac{ı\sigma k}{c_s}}\left[\mathrm{\Gamma }+(\gamma 1)\chi \right]\right),`$
where in the second line we have eliminated the viscosity. The $`|s_\pm `$ lead to

$`𝒩_{s_\pm }`$ $`=`$ $`{\displaystyle \frac{(\gamma 1)^2\rho |p^2}{2s_\pm |s_\pm }}\left(1+2ık\left[B_{s_\pm ,+}{\displaystyle \frac{\rho |\psi _+^{(0)}}{\rho |s_\pm ^{(0)}}}+B_{s_\pm ,}{\displaystyle \frac{\rho |\psi _{}^{(0)}}{\rho |s_\pm }}\right]\right),`$ (93)
$`=`$ $`{\displaystyle \frac{(\gamma 1)\rho |\rho }{2\gamma }}\left(1\pm {\displaystyle \frac{\xi }{\sqrt{1+\xi ^2}}}\right),`$
where the imaginary part cancels completely.

## Appendix B

In general one defines a transport coefficient via the decay of small fluctuations with respect to the equilibrium distribution. If we take $`\delta f(𝒌)`$ to be such a fluctuation, we obtain for a single timestep

$$e^{z(𝒌)}=\frac{\delta f|e^{ı𝒌𝒄}(\mathrm{𝟏}+\mathrm{\Omega })|\delta f}{\delta f|\delta f},$$
(94)
where the second order term in $`k`$ of $`z(𝒌)`$ will be the transport coefficient.

Naturally these fluctuations can always be written in terms of the eigenfunctions of the linearized Boltzmann operator

$$|\delta f(𝒌)=\underset{\mu }{}\frac{\psi _\mu |\delta f}{\psi _\mu |\psi _\mu }|\psi _\mu (𝒌),$$
(95)
where the coefficients on the right-hand site are independent of $`𝒌`$ and therefore can be obtained from the $`k=0`$ limit $`\psi _\mu |\delta f/\psi _\mu |\psi _\mu =\psi _\mu ^{(0)}|\delta f^{(0)}/\psi _\mu ^{(0)}|\psi _\mu ^{(0)}`$. Combining this formal expansion with Eq. (94) and using the results of the eigenvalue problem (13) we obtain

$$e^{z(𝒌)}=\frac{1}{\delta f|\delta f}\underset{\mu }{}\frac{\psi _\mu |\delta f}{\psi _\mu |\psi _\mu }\delta f|e^{z_\mu }|\psi _\mu .$$
(96)
If we restrict ourselves to fluctuations proportional to the hydrodynamic modes, and thus assume the exponential fast decay of the kinetic modes, the eigenvalues $`e^{z_\mu }`$ can be expanded in terms of $`𝒌`$. This results in

$$e^{z(𝒌)}=1+\underset{\mu }{}\frac{\psi _\mu |\delta f\delta f|\psi _\mu }{\psi _\mu |\psi _\mu \delta f|\delta f}\left(ıkz_\mu ^{(1)}k^2\left[z_\mu ^{(2)}+\frac{1}{2}(z_\mu ^{(1)})^2\right]\right)+𝒪(k^3).$$
(97)
Note that the prefactors also depend on $`𝒌`$, but can be rewritten as

$`{\displaystyle \frac{\psi _\mu |\delta f\delta f|\psi _\mu }{\psi _\mu |\psi _\mu \delta f|\delta f}}=\left({\displaystyle \frac{\psi _\mu ^{(0)}|\delta f^{(0)}}{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}}}\right)^2{\displaystyle \frac{\psi _\mu |\psi _\mu }{\delta f|\delta f}}=`$
$`{\displaystyle \frac{\psi _\mu ^{(0)}|\delta f^{(0)}\delta f^{(0)}|\psi _\mu ^{(0)}}{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}\delta f^{(0)}|\delta f^{(0)}}}+𝒪(k^2).`$ (98)
This allows us to calculate the lowest order terms of $`z(𝒌)`$ via

$$z^{(1)}=\underset{\mu }{}\frac{\psi _\mu ^{(0)}|\delta f^{(0)}\delta f^{(0)}|\psi _\mu ^{(0)}}{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}\delta f^{(0)}|\delta f^{(0)}}z_\mu ^{(1)}$$
(99)
$$z^{(2)}+\frac{1}{2}\left(z^{(1)}\right)^2=\underset{\mu }{}\frac{\psi _\mu ^{(0)}|\delta f^{(0)}\delta f^{(0)}|\psi _\mu ^{(0)}}{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}\delta f^{(0)}|\delta f^{(0)}}\left(z_\mu ^{(2)}+\frac{1}{2}(z_\mu ^{(1)})^2\right)$$
(100)
It can easily be checked that in the case of a single eigenmode for the fluctuation $`\delta f`$ these equations reduce to the normal results, in particular Eq. (33). In general, however, such a formulation in terms of currents is not valid. The fact that it nevertheless is consistent for Eqns. (41) and (42) is a consequence of the orthogonality relation (39) between the currents of the two diffusive eigenmodes.

In the case of the different diffusions calculated in Sec. IV.3 we need to make use of these equations, because the fluctuations under consideration are not eigenmodes. From symmetry arguments it follows that in all those cases $`z^{(1)}=0`$, and strictly speaking none of these fluctuations will therefore be propagating. In the case of Eqns. (70) and (72), however, the soundmodes will contribute to the transport value (100) and spectra based on these fluctuations will contain Brillouin peaks.

We finally like to mention that these results are only valid in the limit $`k0`$ and for small times, because if one considers the decay for larger time intervals one would obtain

$$e^{z(𝒌)t}=\frac{\delta f|\left[e^{ı𝒌𝒄}(\mathrm{𝟏}+\mathrm{\Omega })\right]^t|\delta f}{\delta f|\delta f}=\underset{\mu }{}\left(\frac{\psi _\mu ^{(0)}|\delta f^{(0)}}{\psi _\mu ^{(0)}|\psi _\mu ^{(0)}}\right)^2\frac{\psi _\mu |\psi _\mu }{\delta f|\delta f}e^{z_\mu t}.$$
(101)
This results in a different behavior for short and long times in the case one considers fluctuations that are not proportional to a single eigenmode (See also Ref. Mountain and Deutch (1969)).


text:
# Universal measurement apparatus controlled by quantum software

## Abstract

We propose a quantum device that can approximate any projective measurement on a qubit. The desired measurement basis is selected by the quantum state of a “program register”. The device is optimized with respect to maximal average fidelity (assuming uniform distribution of measurement bases). An interesting result is that if one uses two qubits in the same state as a program the average fidelity is higher than if he/she takes the second program qubit in the orthogonal state (with respect to the first one). The average information obtainable by the proposed measurements is also calculated and it is shown that it can get different values even if the average fidelity stays constant. Possible experimental realization of the simplest proposed device is presented.

Programmable quantum “multimeters” are devices that can realize any desired generalized quantum measurement from a chosen set (either exactly or approximately) DuBu . Their main feature is that the particular positive operator valued measure (POVM) is selected by the quantum state of a “program register” (quantum software). In this sense they are analogous to universal quantum processors Nielsen97 ; Vidal00 ; Hillery02 . Quantum multimeters could play an important role in quantum state estimation and quantum information processing.

In this paper, we will describe a programable quantum device that can approximately accomplish any projective von Neumann measurement on a single qubit. Since it is impossible to encode an arbitrary unitary operation (acting on a finite-dimensional Hilbert space) into a state of a finite-dimensional quantum system Nielsen97  it is also impossible to encode arbitrary projective measurement on a qubit into such a state DuBu . However, it is still possible to encode POVM’s that represent, in a certain sense, the best approximation of the required projective measurements.

Suppose we would like to measure a qubit in the basis represented by two orthogonal vectors $`|\psi `$ and $`|\psi _{}`$. We want this measurement basis be controlled by the quantum state of a program register, $`|\varphi _p(\psi )`$. An *ideal* multimeter would map the composite state of the measured system and the program register to two fixed orthogonal pure states $`|0`$ and $`|1`$ according to

$$|\psi |\varphi _\mathrm{p}(\psi )|0,|\psi _{}|\varphi _\mathrm{p}(\psi )|1.$$
(1)
As mentioned above, such transformation cannot be implemented exactly. Thus, our task is to find a realistic linear trace-preserving completely positive (CP) map that represents the closest approximation to this non-realistic map. We focus on the scenario when we always obtain one of the two measurement results $`|0`$ or $`|1`$ but errors, i.e. deviations from the ideal map (1), may appear. Our aim is to minimize the probability of error, i.e., we will maximize the probability of the correct discrimination between states $`|\psi `$ and $`|\psi _{}`$. In general we could optimize both the program and the fixed transformation so as to optimally approximate map (1) for a given dimension of the program register. However, this is an extremely hard problem that we will not attempt to solve in its generality. Instead, we shall optimize the fixed transformation for two natural choices of the program.

First we shall assume that the program register contains $`N`$ copies of the state $`|\psi `$, $`|\varphi _\mathrm{p}(\psi )=|\psi ^N`$. Our second choice of the program – the two-qubit state $`|\psi |\psi _{}`$ – is motivated by recent results on optimum quantum state estimation. Gisin and Popescu showed that the state of two orthogonal qubits $`|\psi |\psi _{}`$ encodes the information on the state $`|\psi `$ better than state of two identical qubits $`|\psi |\psi `$ Gisin99 . If we possess one copy of the state $`|\psi |\psi _{}`$ than we can estimate $`|\psi `$ with fidelity $`_{}=(1+1/\sqrt{3})/20.7886`$ which is slightly higher than the fidelity of optimum estimation on one copy of two identical qubits, $`_{||}=3/4`$ Gisin99 ; Massar00 . One would thus expect that the state $`|\psi |\psi _{}`$ should also give an advantage when used as a program of the multimeter. Rather surprisingly, this is not the case and we shall see that it is better to use two identical qubits $`|\psi |\psi `$.

In what follows we shall benefit from the isomorphism between CP maps and bipartite positive semidefinite operators Jamiolkowski72 ; Fiurasek01 . Let $``$ and $`𝒦`$ denote the Hilbert spaces of input and output states, respectively. Choose basis $`|i`$ in $``$, define a maximally entangled state $`_i|i|i`$ on $`^2`$ and apply the CP map to one part of this state. The density matrix $`\chi `$ of the resulting state on Hilbert space $`𝒦`$ represents the CP map and the relation between input and output density matrices reads Fiurasek01

$$\rho _{\mathrm{out}}=\mathrm{Tr}_{\mathrm{in}}[\chi \rho _{\mathrm{in}}^T𝟙_{\mathrm{out}}]$$
(2)
where $`T`$ stands for the transposition in the basis $`|i`$ and $`𝟙`$ denotes an identity operator. The CP map is trace preserving if the positive semidefinite operator $`\chi `$ satisfies the condition $`\mathrm{Tr}_{\mathrm{out}}[\chi ]=𝟙_{\mathrm{in}}`$.

Let us define the fidelity $`F(\psi )`$ of our multimeter projecting onto states $`|\psi `$ and $`|\psi _{}`$ as the probability that a correct measurement result will be obtained when we send the states $`|\psi `$ or $`|\psi _{}`$ to the input randomly each with probability one half. This probability can be interpreted as a success rate of the discrimination between two orthogonal states $`|\psi `$ and $`|\psi _{}`$. Assuming the program state to be $`|\psi ^N`$, the two relevant input states of the multimeter read

$$|\mathrm{\Psi }=|\psi |\psi ^N,|\mathrm{\Psi }_{}=|\psi _{}|\psi ^N.$$
(3)
The input Hilbert space of the multimeter is a tensor product of the Hilbert space of signal qubit $`_\mathrm{s}`$ and symmetric (bosonic) subspace $`_+^N`$ of the Hilbert space of $`N`$ qubits, $`=_\mathrm{s}_+^N`$ and $`\mathrm{dim}=2(N+1)`$. A trace preserving completely positive map $`\chi `$ transforms the input states onto the states of a single output qubit that is subsequently measured in the computational basis. The outcome $`|0`$ corresponds to the projection onto $`|\psi `$ while $`|1`$ is associated with the projection onto $`|\psi _{}`$. Making use of the input-output relation (2) we have,

$`F(\psi )`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}\left[\chi (|\mathrm{\Psi }\mathrm{\Psi }|)^T|00|\right]`$ (4)
$`+`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}\left[\chi (|\mathrm{\Psi }_{}\mathrm{\Psi }_{}|)^T|11|\right].`$
The figure of merit that we would like to maximize is the mean fidelity obtained on averaging $`F(\psi )`$ over all pure qubit states $`|\psi `$, i.e., over the surface of the Bloch sphere,

$$F=_\psi 𝑑\psi F(\psi )=\mathrm{Tr}[\mathrm{R}\chi ].$$
(5)
The positive semidefinite operator $`R`$ reads

$$R=R_+|00|+R_{}|11|,$$
(6)
where the operators $`R_+`$ and $`R_{}`$ acting on the input Hilbert space $``$ are given by integrals

$$R_+^T=\frac{1}{2}_\psi 𝑑\psi |\mathrm{\Psi }\mathrm{\Psi }|,R_{}^T=\frac{1}{2}_\psi 𝑑\psi |\mathrm{\Psi }_{}\mathrm{\Psi }_{}|.$$
(7)
A straightforward calculation reveals that $`R_+`$ is proportional to the projector onto symmetric subspace $`_+^{(N+1)}`$ of the Hilbert space of $`N+1`$ qubits,

$$R_+=\frac{1}{2(N+2)}\mathrm{\Pi }_+^{(N+1)}.$$
(8)
Furthermore, the sum of operators $`R_+`$ and $`R_{}`$ is proportional to the identity operator on the input Hilbert space, $`R_++R_{}=𝟙_{\mathrm{in}}/[\mathrm{𝟚}(+\mathrm{𝟙})]`$. Thus we immediately have

$$R_{}=\frac{1}{2(N+1)}𝟙_{\mathrm{in}}\frac{\mathrm{𝟙}}{\mathrm{𝟚}(+\mathrm{𝟚})}\mathbb{\Pi }_+^{(+\mathrm{𝟙})}.$$
(9)
The determination of the optimum CP map amounts to maximization of the linear function (5) under the constraints $`\chi 0`$ and $`\mathrm{Tr}_{\mathrm{out}}[\chi ]=𝟙_{\mathrm{in}}`$. The optimum CP map that maximizes the mean fidelity (5) must satisfy the extremal equations Fiurasek01 ; Audenaert01

$`(\lambda 𝟙_{\mathrm{out}})\chi =\mathrm{𝟘},`$ (10)
$`\lambda 𝟙_{\mathrm{out}}\mathrm{𝟘},`$ (11)
where $`\lambda `$ is positive definite operator on the input Hilbert space. Notice that the extremal equations (10) and (11) resemble the Helstrom equations for optimum POVM that maximizes the success rate in ambiguous quantum state discrimination helst . One can prove that if both Eqs. (10) and (11) are satisfied, then $`\chi `$ is indeed optimum CP map and $`F`$ attains its global maximum on the convex set of trace-preserving CP maps Audenaert01 ; Fiurasek01b .

The extremal equations can be very efficiently solved numerically by means of repeated iterations Fiurasek01 . From the numerical results, we were able to conjecture the optimum CP map,

$$\chi =\mathrm{\Pi }_+^{(N+1)}|00|+\mathrm{\Pi }_{}^{(N+1)}|11|.$$
(12)
where $`\mathrm{\Pi }_{}^{(N+1)}=𝟙_{\mathrm{in}}\mathbb{\Pi }_+^{(+\mathrm{𝟙})}`$ is a projector onto subspace orthogonal to the symmetric subspace of $`N+1`$ qubits. On inserting the expressions (6) and (12) into Eq. (5), we obtain the mean fidelity

$$F=\frac{2N+1}{2N+2}.$$
(13)
Prior to analyzing the properties of the CP map (12) in detail, let us prove its optimality. Taking into account the trace-preservation condition $`\mathrm{Tr}_{\mathrm{out}}[\chi ]=𝟙_{\mathrm{in}}`$, we find from Eq. (10) that $`\lambda =\mathrm{Tr}_{\mathrm{out}}[R\chi ]`$. After some algebra we arrive at

$$\lambda =\frac{1}{2(N+1)}𝟙_{\mathrm{in}}\frac{\mathrm{𝟙}}{\mathrm{𝟚}(+\mathrm{𝟙})(+\mathrm{𝟚})}\mathbb{\Pi }_+^{(+\mathrm{𝟙})}.$$
(14)
It is easy to check that the first extremal equation (10) is satisfied. The inequality (11) splits into two independent inequalities for operators acting on input Hilbert space,

$$A_1=\lambda R_{}0,A_2=\lambda R_+0.$$

A straightforward calculation yields explicit formulas for $`A_1`$ and $`A_2`$,

$`A_1`$ $`=`$ $`{\displaystyle \frac{N}{2(N+1)(N+2)}}\mathrm{\Pi }_+^{(N+1)},`$
$`A_2`$ $`=`$ $`{\displaystyle \frac{1}{2(N+1)}}\mathrm{\Pi }_{}^{(N+1)}.`$ (15)
Since both these operators are positive semidefinite, the condition (11) is satisfied. This concludes our proof of the optimality of the CP map (12).

The structure of the optimum CP map (12) indicates that this map is a joint generalized quantum measurement on the signal qubit and $`N`$ program qubits and the corresponding POVM has two elements $`\mathrm{\Pi }_+^{(N+1)}`$ and $`\mathrm{\Pi }_{}^{(N+1)}`$. We can now determine the effective POVM carried out on the signal qubit,

$`\mathrm{\Pi }_{||}`$ $`=`$ $`\mathrm{Tr}_\mathrm{p}\left[𝟙_\mathrm{s}(|\psi \psi |)^{}\mathbb{\Pi }_+^{(+\mathrm{𝟙})}\right],`$
$`\mathrm{\Pi }_{}`$ $`=`$ $`\mathrm{Tr}_\mathrm{p}\left[𝟙_\mathrm{s}(|\psi \psi |)^{}\mathbb{\Pi }_{}^{(+\mathrm{𝟙})}\right],`$ (16)
where $`\mathrm{Tr}_\mathrm{p}`$ denotes trace over the program qubits. The outcome $`\mathrm{\Pi }_{}`$ cannot occur if the input state is $`|\psi `$ because the input state $`|\mathrm{\Psi }`$ belongs to the symmetric subspace of $`N+1`$ qubits and $`\mathrm{\Pi }_{||}`$ clicks with certainty. Hence the POVM element $`\mathrm{\Pi }_{}`$ must be proportional to the projector $`|\psi _{}\psi _{}|`$. Since the sum of POVM elements (16) is an identity operator, we have the following ansatz,

$`\mathrm{\Pi }_{||}`$ $`=`$ $`|\psi \psi |+(1p)|\psi _{}\psi _{}|,`$
$`\mathrm{\Pi }_{}`$ $`=`$ $`p|\psi _{}\psi _{}|.`$ (17)
The probability $`p`$ that $`\mathrm{\Pi }_{}`$ clicks when the input state is $`|\psi _{}`$ is given by $`p=\mathrm{\Psi }_{}|\mathrm{\Pi }_{}^{(N+1)}|\mathrm{\Psi }_{}`$. After some algebra we get $`p=N/(N+1)`$ and the effective POVM representing our universal multimeter reads

$`\mathrm{\Pi }_{||}`$ $`=`$ $`{\displaystyle \frac{1}{N+1}}𝟙+{\displaystyle \frac{}{+\mathrm{𝟙}}}|\psi \psi |,`$
$`\mathrm{\Pi }_{}`$ $`=`$ $`{\displaystyle \frac{N}{N+1}}|\psi _{}\psi _{}|.`$ (18)
Notice that the POVM (18) is asymmetric, which reflects the asymmetry of the program register. Furthermore, the fidelity $`F(\psi )`$ is independent of $`\psi `$ and equal to the mean fidelity (13). In the limit of infinitely large program register $`(N\mathrm{})`$, the POVM (18) approaches the ideal projective measurement.

We now turn our attention to the program $`|\psi |\psi _{}`$. The optimum CP map for this program can be found following the same procedure as described above for the program $`|\psi ^N`$. Briefly, one has to calculate the operator $`R`$ and solve extremal equations (10) and (11). We will not give the details of calculations here and only present the results. Similarly as before, the optimum CP map is a generalized measurement on the signal qubit and two program qubits. The two elements of this three-qubit POVM read,

$`\mathrm{\Pi }_+`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Pi }_+^{(3)}+|\varphi _1\varphi _1|+|\varphi _2\varphi _2|,`$
$`\mathrm{\Pi }_{}`$ $`=`$ $`𝟙_\mathrm{𝟛}\mathbb{\Pi }_+,`$ (19)
where $`𝟙_\mathrm{𝟛}`$ is an identity operator on Hilbert space of three qubits and

$`|\varphi _1`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{3}}}[(\sqrt{3}+1)|0_\mathrm{s}|01_\mathrm{p}(\sqrt{3}1)|0_\mathrm{s}|10_\mathrm{p}`$
$`2|1_\mathrm{s}|00_\mathrm{p}],`$
$`|\varphi _2`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{3}}}[(\sqrt{3}+1)|1_\mathrm{s}|10_\mathrm{p}(\sqrt{3}1)|1_\mathrm{s}|01_\mathrm{p}`$ (20)
$`2|0_\mathrm{s}|11_\mathrm{p}].`$
Here the subscripts s and p label the states of signal and program qubits, respectively. After some algebra, we find the effective POVM carried out on the signal qubit,

$`\mathrm{\Pi }_{||}^{}`$ $`=`$ $`{\displaystyle \frac{3\sqrt{3}}{6}}𝟙+{\displaystyle \frac{\sqrt{\mathrm{𝟛}}}{\mathrm{𝟛}}}|\psi \psi |,`$
$`\mathrm{\Pi }_{}^{}`$ $`=`$ $`{\displaystyle \frac{3\sqrt{3}}{6}}𝟙+{\displaystyle \frac{\sqrt{\mathrm{𝟛}}}{\mathrm{𝟛}}}|\psi _{}\psi _{}|.`$ (21)
This POVM is symmetric (reflecting the symmetry of the program $`|\psi |\psi _{}`$). The fidelity $`F(\psi )`$ is state independent and equal to the mean fidelity

$$F^{}=\frac{1}{2}\left(1+\frac{1}{\sqrt{3}}\right).$$
(22)
Notice that $`F^{}=_{}`$. This is not a mere coincidence, the optimum strategy for program $`|\psi |\psi _{}`$ is to carry out an optimal estimation of $`|\psi `$ and then measure the signal qubit in the basis formed by estimated state $`|\psi _{\mathrm{est}}`$ and its orthogonal counterpart. The POVM (19) is an explicit implementation of this procedure. We emphasize here that $`F^{}`$ is a maximum fidelity attainable with program $`|\psi |\psi _{}`$, because the corresponding CP map solves the extremal equations (10) and (11). With the program $`|\psi |\psi `$ we achieve the fidelity $`5/60.8333`$ which is higher than $`F^{}`$, hence the program $`|\psi |\psi `$ exhibits better performance than $`|\psi |\psi _{}`$.

The aim of the measurement is to obtain some information on the system. It can be illustrated on the following example: Alice encodes bits of a message into quantum states and Bob tries to read the message making measurements on these states. The average Bob’s information (per bit) on the message can be quantified as

$$I=\underset{i}{}r_i\left[1+p_i\mathrm{log}_2p_i+(1p_i)\mathrm{log}_2(1p_i)\right],$$
(23)
where $`r_i`$ is a fraction of bits that Bob knows with probability $`p_i`$ ($`_ir_i=1`$). It is reasonable to look for measurements that maximize average information per . However, as the average information is a non-linear function of probabilities this is usually a very difficult problem. We will not try to search for a CP map (or POVM) that would maximize information instead of fidelity. Nevertheless, we will calculate the average information obtainable by the measurements described above that are optimal with respect to the fidelity. Let us define fidelities $`F_{||}=\psi |\mathrm{\Pi }_{||}|\psi `$ and $`F_{}=\psi _{}|\mathrm{\Pi }_{}|\psi _{}`$. If the average rates of both the states are the same and equal to $`1/2`$ (generalization to an arbitrary rate is straightforward) then

$`r_{||}=\left[F_{||}+(1F_{})\right]/2,p_{||}=F_{||}/(2r_{||}),`$
$`r_{}=\left[F_{}+(1F_{||})\right]/2,p_{}=F_{}/(2r_{}).`$ (24)
For example, if $`F_{||}=1`$ and $`F_{}=1/2`$ Bob gets result “0” always when Alice sends $`|\psi `$ and in half cases when she sends $`|\psi _{}`$. So, he knows the original bit with probability $`p_{||}=2/3`$ and the average relative number of these results is $`r_{||}=3/4`$. For result “1” the probability $`p_{}=1`$ (this result appears only if Alice sends $`|\psi _{}`$) and $`r_{}=1/4`$. This leads to average information $`I0.311`$. The described situation corresponds to the above mentioned measurement with the program state $`|\psi `$. Note that a symmetrized measurement with the *same* average fidelity ($`F_{||}=F_{}=3/4`$) would lead to *lower* information ($`I0.189`$).

For the case when the program state is $`|\psi |\psi `$ and $`F_{||}=1`$, $`F_{}=2/3`$ the average information gets value $`I0.459`$. For the program state $`|\psi |\psi _{}`$, when $`F_{||}=F_{}=\frac{1}{2}(1+\frac{1}{\sqrt{3}})`$, the information is $`I0.256`$. Notice that this information is *lower* than the one obtained with the device programmed by the single-qubit program state $`|\psi `$ ($`0.256<0.311`$).

The proposed universal measurement devices can be, at least in principle, realized experimentally. As an example, we describe now a possible realization of the simplest one that uses a single-qubit program register. Such a device can be built up from one Fredkin (controlled swap) and two Hadamard gates as shown in Fig. 1 Ekert01 . Signal and program enter the Fredkin gate as “controlled” qubits; the result can be read out from the ancilla serving as a “control” qubit. There are several ways how to implement the presented scheme on real physical systems radim . For example, we can employ an experimental setup suggested previously in cavity-QED experiments Raimond97  that is sketched in Fig. 2. The Fredkin gate is realized by two cavities C0 and C1 coupled by a waveguide and it is controlled by an atomic “probe” affecting cavity field through a nonlinear atom-field interaction in a large-detuning limit Brune92 . This gate is inserted between two Hadamard gates implemented as Ramsey atomic interferometers. The state of cavity C1 represents the signal, the state of cavity C0 the program. The program state can be set employing an atom-field interaction too. An atom is prepared in the desired program state $`a|\mathrm{g}+b|\mathrm{e}`$. This is done in zones B and R1 (atoms come from an oven O and only those with an appropriate velocity are selected in zone V). An atom transition frequency must coincide with the cavity frequency. Passing through the cavity C0 the atom changes its photonic state to $`a|0+b|1`$. Once the program and the signal are encoded into the states of the cavities, 50:50 coupling between C0 and C1 is established. Then the second atom, initially prepared in the excited state, is sent through the Ramsey interferometer. This atom is exposed to the nonlinear interaction with the cavity field in a large-detuning limit (in contrast to the first one). Finally, the result of our programmed measurement process encoded in the state of the atom can be revealed by the field-ionization counters De end Dg.

###### Acknowledgements.

This research was supported under the project LN00A015 of the Ministry of Education of the Czech Republic.