format: text

text:
# Tunneling and Electric-Field Effects on Electron-Hole Localization in Artificial Molecules

## Abstract

I theoretically investigate the Stark shift of the exciton goundstate in two vertically coupled quantum dots as a function of the interdot distance. The coupling is shown to enhance the tuneability of the linear optical properties, including energy and oscillator strength, as well as the exciton polarizability. The coupling regime that maximizes these properties results from the detailed balance between the effects of the single-particle tunneling, of the electric field and of the carrier-carrier interaction. I discuss the relevance of these results to the possible implementation of quantum-information processing based on semiconductor quantum dots: in particular, I suggest the identification of the qubits with the exciton levels in coupled- rather than single-dots.

Semiconductor quantum dots (QDs) are quasi-zero-dimensional heterostructures, where the strong spatial confinement of electrons and holes gives rise to discrete optical spectra hawrylak:98 . As compared to natural ones, such “artificial atoms” show a wide tuneability of the structural, electronic and optical properties, of the coupling with each other and with external fields. This makes them an ideal laboratory for the many-body physics, where the relative importance of the different kinds of interaction can be modified in order to produce regimes which are precluded to atomic and molecular physics (see, e.g., Refs. troiani:02 ; rontani:01 , and references therein). In particular, the application of an electric field on single and coupled QDs and the resulting quantum confined Stark effect (QCSE) has been theoretically barker:01 ; sheng:02 ; janssens:02  and experimentally fry:00 ; findeys:01 ; patane:00 ; passaseo:01  investigated, in order to gain insight into the spatial distribution of the carriers and into the effects of possible strain-induced or built-in electric fields.

In recent years an increasing attention has been devoted to QDs as candidates for a solid-state implementation of quantum-information processing (QIP). Many possible schemes have been proposed, involving the use of either spin loss:98  or orbital barenco:95  degrees freedom as qubits. The latter include those that identify the computational degrees of freedom with specific exciton levels and propose to manipulate them on the ps timescale by means of coherent-carrier control techniques troiani:00 ; biolatti:00 ; chen:01 . In each of these QD-based schemes the coupling between dots is a crucial issue and is generally supposed to be tuneable by means of external electric or magnetic fields. In particular, the use of the electric field has been suggested in order to induce relevant dipoles in the exciton states, which would result in a dipole-dipole interaction between excitons (qubits) localized in neighbouring dots biolatti:00 . The more critical aspects of similar proposals are related to the need of spectrally resolving each single qubit within a large array and to the efficient decohering phenomena, such as the radiative recombination of the electron-hole pairs decay , that result in decoherence times between tens and hundreds of ps. In these two respects an increased tuneability of the exciton energy and of the oscillator strength might represent a relevant improvement: one of the main goals of the present work is to demonstrate that this can be achieved by switching from single to coupled QDs. In particular, I will focus on the range of electric-field values that is compatible with finite exciton oscillator strength, as required in order to allow optical manipulation of the exciton states.

As a first step, I compute the single-particle eigenfunctions and eigenvalues for electrons and holes within the envelope-function and effective-mass approximations hawrylak:98 . The external confinement of the carriers in the double dot is described by a prototypical confinement potential, which is double-box-like along $`z`$ and parabolic in the $`(x,y)`$-plane; an additional term in the single-particle Hamiltonian $`H_{SP}`$ accounts for the electric field $`F`$, directed along $`z`$:

$`H_{SP}={\displaystyle \underset{\alpha =e,h}{}}{\displaystyle }\widehat{\psi }_\alpha ^{}(𝒓)[{\displaystyle \frac{\mathrm{}^2}{2m_\alpha ^{}}}^2+V_\alpha ^{DW}(z)+`$
$`{\displaystyle \frac{1}{2}}m_\alpha ^{}\omega _\alpha ^2(x^2+y^2)q_\alpha Fz]\widehat{\psi }_\alpha (𝒓)d𝒓,`$ (1)
where $`\widehat{\psi }_e(𝒓)`$ and $`q_e`$ ($`\widehat{\psi }_h(𝒓)`$ and $`q_h`$) are the electron (hole) field operators and electric charge respectively, $`V_{e,h}^{DW}(z)=V_{e,h}^0[\theta (|z|d/2l)+\theta (d/2|z|)]`$, $`l=10`$ nm and $`d`$ are the width of the wells and the interdot barrier respectively, $`V_{e,h}^0=400,215`$ meV are the well depths for electrons and holes.

Starting from these ingredients, I compute the interacting electron-hole states by means of a full configuration interaction (CI) calculation banyai:93 . For each number of electron and holes, I generally truncate the Hilbert space by keeping only the $`N_H`$ configurations of lowest single-particle energy, $`N_H`$ beeing determined on the grounds of a convergence criterion. Within such subspace, I diagonalize the Hamiltonian matrix $`H_{SP}+H_C`$, where

$$H_C=\frac{1}{2}\underset{\alpha \beta =e,h}{}\frac{d𝒓d𝒓^{}}{(\kappa |𝒓𝒓^{}|)}\widehat{\psi }_\alpha ^{}(𝒓^{})\widehat{\psi }_\beta ^{}(𝒓)\widehat{\psi }_\beta (𝒓)\widehat{\psi }_\alpha (𝒓^{})$$
(2)
accounts for the carrier-carrier Coulomb interactions and for the resulting correlation effects ($`\kappa `$ is the semiconductor dielectric constant, $`𝒓𝒓^{}`$). In the investigation of the single exciton ground state 20 SP states for both electrons and holes are kept and no further truncation of the Hilbert space is performed (therefore $`N_H=400`$ configurations).

As a preliminary step, let me consider the effects of the interplay between tunneling and electric field on the SP states, which is described by the $`z`$-dependent part of $`H_{SP}`$. In Fig. 1 I plot the distributions along $`z`$ corresponding to the electron and hole SP groundstates as a function of $`F`$ and of the interdot distance $`d`$. These eigenstates result from a competition between the coherent coupling, that tends to delocalize the wavefunctions over the double dot, and the field, that instead tends to localize them in one dot or in the other. While the tunneling-induced energy splitting between the bonding and antibonding states is conventionally denoted by $`2t`$, the one that the field induces between states localized in two different dots is approximately given by $`2\mathrm{\Delta }e(d+l)F`$. The ratio $`x\mathrm{\Delta }/t`$ thus provides a good estimate of the relative strength of such two effects and of the resulting degree of localization of each carrier: this is found to increase from electrons to holes (typically $`m_e^{}<m_h^{}`$ and thus $`t_e>t_h`$) and with increasing $`d`$ (see Fig. 1(a-d)). In the same range of values of $`F`$ the SP eigenstates and eigenfunctions corresponding to a single dot (not reported here) are hardly affected by the field. Within the present picture, such a difference can be understood by thinking of the single QD as a double dot in the in the limit $`d0`$: the values of $`x\mathrm{\Delta }/t`$ corresponding to an electric field $`F8`$ kV/cm turn out to be very small both for electrons ($`x_e0.038`$) and for holes ($`x_h0.18`$).

The structural parameters which determine the carrier confinement in each dot, namely $`l`$ and $`\mathrm{}\omega _{e,h}`$, are kept constant throughout the following calculations: while a detailed understanding of their role would require a full exploration of the parametric space, a few comments can be done on the grounds of the present results. In particular, the possible reduction of the well width $`l`$ would mainly result in an enhancement of the carrier penetration in the barrier and therefore of the interdot coupling and of $`t`$, whose role and interplay with the electric field are discussed above. The in-plane confinement, instead, does not directly affect the interdot coupling (and therefore $`t`$ or $`x`$). A possible increase of $`\mathrm{}\omega _{e,h}`$ would essentially produce an enhancement of the intradot Coulomb matrix elements and therefore of the Coulomb-induced effects (see the discussion below).

The properties of the exciton states result from the interplay between the abovementioned features of the SP states and the electron-hole Coulomb interaction. In Fig. 2 I plot the Stark shifts and the oscillator strengths of the exciton groundstate $`|X_0`$ as a function of $`F`$ and $`d`$ (panel (a)) and the contributions to the total energy arising from $`H_{SP}`$ and $`H_C`$ (Fig. 2(b-c)); the results for the single dot are also reported for a comparison. Quite generally, the magnitude of the QCSE is greatly enhanced and the wavefunctions are much more affected by the field in coupled dots as compared to single ones. While in single QDs the Stark shift $`\mathrm{\Delta }EE(F)E(0)`$ typically shows a quadratic dependence on $`F`$, due to the perturbative nature of the contribution from the field, a nonparabolic behaviour is expected in the case of “artificial molecules” sheng:02 ; janssens:02 . The field induces a monotonic decrease of the single-particle energy $`E_{SP}`$, which determines the negative sign of the Stark shift, and a progressive spatial separation of the carriers, as demonstrated by the increasing Coulomb energy $`E_C`$ and oscillator strengths (Fig. 2(a-c)).

A closer inspection of the plots reveals the occurence of two regimes. In the strong-field region ($`F6`$ kV/cm) the overlap between the electron and hole distributions rapidly decraeses and $`\mathrm{\Delta }E`$ depends quasi linearly on $`F`$ (as in the case of a field-rigid dipole interaction): both trends are enhanced at larger values of the interdot distance $`d`$. In the weak-field region ($`F6`$ kV/cm) the interplay between tunneling, Coulomb interaction and electric field is more balanced: rather counterintuitively, the energies $`E`$, $`E_{SP}`$ and $`E_C`$ get closer to the single-dot value as $`d`$ is increased. As discussed in the following, this can be attributed to the fact that the strength of the interdot coupling determines to which extent the carrier-carrier Coulomb interaction can oppose the effects of the electric field.

As discussed in Ref. sheng:02 , the possible presence of strain effects is expected to meaningfully affect the valence-band states and thus the optical properties, e.g., in self-assembled InAs/GaAs QDs. More specifically, the authors show how the biaxial components of the strain field decrease the height of the interdot barrier, increase the energy levels of the confined holes, and induce a stronger confinement of the holes in the upper dot as compared to the lower one. The first two features result in an effective enhancement of the interdot coupling for the holes, the latter implies an asymmetry in the QCSE with respect to the direction of the applied electric field. According to the present analysis, neither of these two effects is required for the anomalous Stark shift to occur in coupled QDs: I thus expect it to be observed also in case of softened or negligible strain fields in and around the dots.

The spatial distribution of the carriers provides some more physical insight in the exciton groundstate $`|X_0`$. In Fig. 3(a)\[(b)\] I plot the probability of finding the electron (hole) in the upper dot for $`d=`$ 1, 2 and 3 nm (red, green and blue lines), and in the upper part of the single dot (black dotted line). As expected from the previous discussion on the SP states, the hole is highly localized in the upper dot already at small values of the field and especially for large values of $`d`$. On the other hand, the behaviour of the electron again reflects the existence of the abovementioned two regimes. In the strong-field region the electron is seen to approach the lower dot, under the direct influence of the field; as observed in Fig. 2 for the energies, the derivatives of the curves in this region ($`|N_{ud}^{e,h}/F|`$) increase with increasing $`d`$. In the weak-field region the Coulomb interaction with the hole dominates on the interaction with the field and the electron tends to be localized in the upper dot. This spatial arrangement is also reflected by the small changes in this region of the oscillator strength and Coulomb interaction (Fig. 2(b-c)). The curves plotted in Fig. 3(a) cross at $`F7`$ kV/cm: both in the weak- and in the strong-field regimes the degree of localization of the electron, in the upper or in the lower dot respectively, is opposed by the tunneling. For $`d=1`$ nm the electron is basically frozen on a delocalized bonding orbital and $`N_{ud}^e`$ does not deviate significantly from $`1/2`$ in the considered range of $`F`$ values; at large interdot distances ($`d=3`$ nm) $`t_e`$ is small enough for the antibonding orbital to be also occupied in order to produce the localization of the electron in either dot, depending on the relative strength of the external- and of the hole-induced electric field.

The polarizability of the exciton groundstate reflects this behaviour. While the value of $`D`$ is on average greatly enhanced in coupled QDs as compared to the single one (black dotted vs. coloured lines in Fig. 3(c)), the sign of $`D/d`$ actually depends on the range of electric-field values of interest. For $`F6`$ kV/cm, where the electron tends to follow the hole in the upper dot, the larger degree of localization resulting from the larger values of $`d`$ tends to suppress the dipole; for $`F6`$ kV/cm, where the electron interaction with the field dominates on that with the hole, a larger degree of localization results in a larger induced dipole. Depending on the detailed balance between SP and Coulomb terms in the Hamiltonian, either the strongly or the weakly coupled QDs show the highest polarizability.

Let me finally give a rough estimate of the dipole-dipole interaction that a field of 8 kV/cm could induce between two pairs of vertically coupled dots, given that $`E_{int}[\stackrel{}{D}_1\stackrel{}{D}_23(\stackrel{}{D}_1\widehat{r}_{12})(\stackrel{}{D}_2\widehat{r}_{12})]/\kappa r_{12}^3`$, and with $`\stackrel{}{D}_1\stackrel{}{D}_2\stackrel{}{r}_{12}`$, $`|\stackrel{}{D}_{1,2}|=10e\times `$nm and $`r_{12}=20`$ nm. The energy shift produced by the exciton-exciton interaction is $`E_{int}1.4`$ meV, i.e. the same order of magnitude of the exciton-exciton interactions in a single QD of the same sizes.

In the perspective of the QD-based QIP schemes, the large enhancements of the exciton polarizability and of the optical tuneability suggest the identification of the computational degrees of freedom with exciton levels in vertically-coupled dots rather than in single ones, as proposed in previous works troiani:00 ; biolatti:00 . The tuneability of the exciton levels by means of local electric fields can be exploited in order to favour the spectral resolution of specific transitions (qubits) within an array of neighbouring dots of similar dimensions and to bring specific transitions in and out of resonance with the laser fields used for the gate implementation. The polarizability of the exciton states, instead, is required in order to induce a relevant dipole-dipole interaction $`E_{int}`$ between excitons that are localized in different and (quantum-mechanically) decoupled dots. In fact, the value of $`E_{int}`$ determines a lower bound to the duration $`\tau _g`$ of the required laser pulses ($`\tau _g\mathrm{}/E_{int}`$): an exciton-exciton interaction of the order of one meV, as the one estimated above, allows in principle to perform the quantum gates on the ps timescale.

The suggested geometry can be integrated in a planar array of vertically coupled dots: each double dot would couple to its nearest neighbours only by means of the field-induced dipole-dipole interaction, while the in-plane SP tunneling should be suppressed in order for the required tensor-product structure of the Hilbert space to be preserved jozsa:97 . As compared to the ones that consider an in-plane electric field, the present geometry has a twofold advantage: (i) it avoids the breaking of the optical selection rules related to the cylinder symmetry of the dot and the resulting thickening of the optica spectra banyai:93 , which interferes with the identification and selective addressing of the single absorbtion lines; (ii) it does not require the etching of a mesa structure and thus a strong limitation on the number of qubits that compose the quantum hardware.

To summarize, I have investigated the effects of a vertical electric field on the exciton groundstate and on the related optical properties of two vertically coupled QDs. The present results suggest the identification of each qubit with the exciton groundstate in a pair of coupled dots rather than in a single one, in view of the higher degree of tuneability of the linear optical properties, including the energy levels and the corresponding oscillator strength. Besides, the polarizability of the lowest exciton states in artificial molecules is enhanced by one order of magnitude or more as compared to the one in artificial atoms, resulting in relevant dipole-dipole interactions between excitons localized in different pairs of coupled dots. More specifically, due to the interplay between tunneling, carrier-carrier Coulomb interactions and electric field, both the tuneability of the optical properties and the polarizability are maximized either in the strong- or in the weak-coupling regime, depending on the range of $`F`$ values.

The author is grateful to E. Molinari and U. Hohenester for stimulating and encouraging discussion. This work has been supported in part by INFM through PRA-99-SSQI and by the EU under the IST programme “SQID”.


text:
# A schematic model for QCD, II: finite temperature regime.

## 1 Introduction

In Ref.  (hereon referred to as (I)) a simple model, representative of QCD, was introduced and applied to the calculation of the spectrum of mesons. It is a Lipkin type model  for the quark sector, coupled to a boson level which is occupied by gluon pairs with spin zero. The four parameters of the model were adjusted in order to reproduce 13 known meson states with spin zero or one. The calculated spectra, for mesons with spin different from the ones used in the fit, were found to be in qualitative agreement with data. As reported in (I), the calculated meson states contain many quarks, antiquarks and gluons. The gluon contributions were found to be of the order of 30$`\%`$. The model predictions (I) are free of the so-called multiplicity problem, i.e. that a given state can be described in many ways, which is removed due to the action of particle mixing interaction. The model itself resembles the one of Ref.  which treats nucleons coupled to pions. Also, it is related to the work of Ref. , which describes quarks and uses particle conserving interactions. Generally speaking, the model of (I) belongs to the class of models described in Refs. . The gluon part in (I) is fixed  and does not contain any new parameters. The validity of the basic theoretical assumptions, and the applications to low and high temperature regimes, has been studied for mesons with flavor (0,0) and spin 0 . The aim of these studies was to formulate a manageable, schematic, albeit realistic, model to describe qualitatively QCD at low and high energies. Since the model is algebraic, i.e. all matrix elements are analytic, and exactly solvable, it can provide a non-perturbative description based on QCD relevant degrees of freedom, like quarks, antiquarks, and gluons. This, in turn, allows to test other microscopic many body techniques previously applied to the non-perturbative treatment of real QCD . Although the proposed model (I) is probably too simple to describe real QCD, it contains all basic ingredients of real QCD. These are the correct number of degrees of freedom associated to color, flavor and spin, and the orbital degree of freedom , which is contained in the degeneracy $`2\mathrm{\Omega }`$ of each of the quark levels.

In this work we investigate the behavior of the model, in the finite temperature regime. By starting from the model predictions of the meson spectrum, we calculate the partition function and different thermodynamical quantities, like the energy density and the heat capacity as a function of temperature. Next, we focus on the calculation of meson production rates. As we shall show, these production rates are in qualitative agreement with the experiments. Also, we calculate absolute production rates for pions and kaons. Finally, we concentrate on the transition from the Quark-Gluon-Plasma (QGP)  to the hadron gas. The results support the notion that the present model may be taken as a toy model for QCD.

The paper is organized as follows: In section 2 the model is shortly outlined, since the details have been presented in (I). In section 3 we calculate the partition function and give the expressions for the relevant observables. In section 4 the model is applied to the description of the QGP. There, we present and discuss the results corresponding to some branching ratios and absolute production rates. Finally, conclusions are drawn in section 5.

## 2 The model

As described in (I), the fermion (quarks and antiquarks) sector of the model consists of two levels at energies $`+\omega _f`$ and $`\omega _f`$, each level with degeneracy $`2\mathrm{\Omega }`$ $`=`$ $`n_cn_fn_s`$, where $`n_c=3`$, $`n_f=3`$ and $`n_s=2`$ are the color, flavor and spin degrees of freedom, respectively (see Fig. 1 of (I)). Each level can be occupied by quarks. Antiquarks are described by holes in the lower level. Equivalently, one can use only the positive energy level and fill it with quarks and antiquarks with positive energy. The Dirac picture is useful because it gives the connection to the Lipkin model as used in nuclear physics. The quarks and antiquarks are coupled to gluon pairs with spin zero. The energy of the gluon level is 1.6 GeV , and the energy $`\omega _f`$ is fixed at the value $`\omega _f=0.33`$ GeV , which is the effective mass of the constituent quarks.

The basic dynamical constituent blocks of the model are quark-antiquark pairs $`𝑩_{\lambda f,SM}^{}`$ which are obtained by the coupling of a quark and an antiquark to flavor $`\lambda `$ ($`\lambda =0,1`$) and spin $`S`$ ($`S=0,1`$). The index $`f`$ is a short hand notation for hypercharge $`Y`$, isospin $`I`$ and its third component $`I_z`$. Under complex conjugation the operator obeys the phase rule defined in .

The states of the Hilbert space can be classified according to the group chain

$`[1^N]`$ $`[h]=[h_1h_2h_3]`$ $`[h^\mathrm{t}]`$
$`U(4\mathrm{\Omega })`$ $`U({\displaystyle \frac{\mathrm{\Omega }}{3}})`$ $`U(12)`$ (1)
$``$ $``$
$`(\lambda _C,\mu _C)SU_C(3)(\lambda _f,\mu _f)`$ $`SU_f(3)SU_S(2)S,M,`$
where the irreducible representation (irrep) of the different unitary groups are attached to the symbols of the groups. The irrep of $`U(4\mathrm{\Omega })`$ is completely antisymmetric (fermions) and the one of $`U(\frac{\mathrm{\Omega }}{3})`$, the color group $`U(3)`$ for $`\mathrm{\Omega }=9`$, and $`U(12)`$ are complementary . The color irrep $`(\lambda _C,\mu _C)`$ of the color group $`SU_C(3)`$ is related to the $`h_k`$ via $`\lambda _C=h_1h_2`$ and $`\mu _C=h_2h_3`$. The reduction of the $`U(12)`$ group to the flavor ($`SU_f(3)`$) and spin group ($`SU_S(2)`$) is obtained by using the procedure described in . In (1) no multiplicity labels are indicated (see (I)).

The classification appearing in (1) is useful to determine the dimension and content of the Hilbert space. Instead of working in the fermion space we have introduced a boson mapping . The quark-antiquark boson operators are mapped to

$`𝑩_{\lambda f,SM}^{}`$ $``$ $`𝒃_{\lambda f,SM}^{}`$
$`𝑩_{\lambda f,SM}`$ $``$ $`𝒃_{\lambda f,SM}.`$ (2)
where the operators on the right hand side satisfy exact boson commutation relations.

The model Hamiltonian is defined completely in the boson space and it is given by

$`𝑯`$ $`=`$ $`2\omega _f𝒏_𝒇+\omega _b𝐧_𝐛+`$ (3)
$`{\displaystyle \underset{\lambda S}{}}V_{\lambda S}\{[(𝒃_{\lambda S}^{})^2+2𝒃_{\lambda S}^{}𝒃_{\lambda S}+(𝒃_{\lambda S})^2](1{\displaystyle \frac{𝒏_𝒇}{2\mathrm{\Omega }}})𝒃+`$
$`𝒃^{}(1{\displaystyle \frac{𝒏_𝒇}{2\mathrm{\Omega }}})[(𝒃_{\lambda S}^{})^2+2𝒃_{\lambda S}^{}𝒃_{\lambda S}+(𝒃_{\lambda S})^2]\},`$
where $`(𝒃_{\lambda S}^{})^2`$ $`=`$ $`(𝒃_{\lambda S}^{}𝒃_{\lambda S}^{})`$ is a short hand notation for the scalar product. Similarly for $`(𝒃_{\lambda S})^2`$ and $`(𝒃_{\lambda S}^{}𝒃_{\lambda S})`$. The factors $`(1\frac{𝒏_𝒇}{2\mathrm{\Omega }})`$ simulate the terms which would appear in the exact boson mapping of the quark-antiquark pairs. The $`𝒃^{}`$ and $`𝒃`$ are boson creation and annihilation operators of the gluon pairs with spin $`S=0`$ and color $`\lambda =0`$. The interaction describes scattering and vacuum fluctuation terms of fermion and gluon pairs. The strength $`V_{\lambda S}`$ is the same for each allowed value of $`\lambda `$ and $`S`$, due to symmetry reasons, as shown in (I). The matrix elements of the Hamilton operator are calculated in a seniority basis. The interaction does not contain terms which distinguish between states of different hypercharge and isospin. It does not contain flavor mixing terms, either. The procedure used to adjust the four parameters (values of $`V_{\lambda S}`$), was discussed in detail in (I).

The disadvantage posed by working in the boson space is the appearance of un-physical states. In (I) we have presented a method which is very efficient to eliminate spurious states, as we shall show in this paper. The suitability of the Hamiltonian (3) to describe the gluon pair and quark-antiquark pair contents of mesons has been discussed in details in (I).

As a next step, in this paper, we have introduced temperature and discussed the transition to and from the QGP. As it can be expected, because of the schematic nature of the model, we may attempt to describe only the general trends of the observables. To achieve this goal, further assumptions have to be made respect to the volume of the system because the model, as it has been proposed in (I), has no a priory information about the volume of the particle.

## 3 The partition function, some state variables and observables

The group classification of the basis (1), allows for a complete book-keeping of all possible states belonging to the Hilbert space of (3). The corresponding partition function, which contains the contribution of the quark-antiquark and gluon pairs configurations introduced in the previous section, is given by

$`Z_{qg_0^2}`$ $`=`$ $`{\displaystyle \underset{[h]}{}}dim(\lambda _C,\mu _C){\displaystyle \underset{(\lambda _f,\mu _f)}{}}{\displaystyle \underset{J}{}}(2J+1){\displaystyle \underset{P=\pm }{}}{\displaystyle \underset{i}{}}`$ (4)
$`mult(E_i)e^{\beta (E_i\mu _BB\mu _ss\mu _II_z)},`$
where $`\mu _B`$, $`\mu _s`$ and $`\mu _I`$ are the baryon, strange and isospin chemical potentials, respectively. The sum over $`[h]=[h_1h_2h_3]`$ denotes all color irreps of $`U(3)`$ with $`_kh_k=N`$, where $`N`$ is the total number of quarks in the two levels (Dirac’s picture). The transposed Young diagram $`[h]^\mathrm{t}`$, obtained by interchanging rows and columns, denotes the $`U(12)`$ irrep. The index ”i” refers to all states with the same color, flavor, spin and parity ($`P`$). These states are obtained after the diagonalization of the Hamiltonian (3). For mesons belonging to the $`\pi `$-$`\eta `$ and $`\omega `$-$`\rho `$ octet, the mass values entering in (4) do not take into account flavor mixing. The eigenvalues $`E_i`$ are denoted by the eigenvalue index $`i`$, and they are functions of all the numbers needed to specify the allowed configurations, namely: $`s`$, $`(\lambda _f,\mu _f)`$, $`P`$, $`S`$, $`[h]`$ and of the cutoff for the different boson species $`[\lambda ,S]`$ (see (I)). The quantities $`B`$ and $`I_z`$ are the baryon number and the third component of the isospin. According to the experimental evidences the value $`\mu _I=0`$ is a reasonable approximation, and we have consistently adopted it in our calculations. The dimension corresponding to color configurations is given by $`mult(\lambda _f,\mu _f)`$ $`=`$ $`\frac{1}{2}(\lambda _f+1)(\mu _f+1)(\lambda _f+\mu _f+2)`$ and of the spin by $`(2J+1)`$.

Since the eigenstates of the Hamiltonian (3) have been calculated after performing a boson mapping, as described in (I), we have consistently fixed the corresponding cutoff-values at the values $`2\mathrm{\Omega }`$, $`\mathrm{\Omega }`$, $`\frac{2\mathrm{\Omega }}{3}`$, and $`\frac{\mathrm{\Omega }}{3}`$ for the boson pair species , ,  and  respectively (see (I)). These values are adequate when the fermion (quark-antiquark) configurations entering in the boson states correspond to a full occupation of the fermion lower state ($`\omega _f`$). These numbers may be modified when the fermion configurations correspond to states where the upper level is partially occupied and the lower level is partially unoccupied. The distribution of the occupation in the upper and lower fermion levels is fixed for a lowest weight state $`|lw>`$ of a given $`U(12)`$ irrep, defined by $`B_{\lambda f,SM}|lw>=0`$. The irrep of $`U(12)`$ is given by a Young diagram  with $`m_k`$ boxes in the $`k`$’th row. The lowest weight state is given by $`_{k=1}^6m_k`$ quarks in the lower level and $`_{k=7}^{12}m_k`$ quarks in the upper level. The highest weight state is obtained by interchanging the occupation. The difference of the number of quarks in the upper level, appearing in the highest and lowest weight states, gives the maximal number of quarks we can excite for a given U(12) irrep. This number is given by

$$2J=\underset{k=1}{\overset{6}{}}m_k\underset{k=7}{\overset{12}{}}m_k.$$
(5)
For the case $`[3^60^6]`$, used in (I), we found $`J=\mathrm{\Omega }`$. Therefore, $`2J`$ is the maximal number of quarks we can shift to the higher level, i.e. it is equal to the maximal number of quark-antiquark pairs which can be put on top of the lowest weight state of a given $`U(12)`$ irrep, which is also the state with the lowest energy in absence of interactions.

The total partition function is given by

$$Z=Z_{qg_0^2}Z_g,$$
(6)
where $`Z_g`$ is the contribution of all gluon states  which do not include contributions of gluon pairs with spin zero. It is written $`Z_g=_\alpha exp(\beta E_\alpha )`$. The values $`E_\alpha `$ can be deduced by using Eq. (40) of Ref. . Except for the gluon pairs with spin zero, all other gluon states are treated as spectators because the Hamiltonian in (3) does not contain interactions with these other states. Note, that the interaction between the gluons is taken into account explicitly in the model of Ref. . As a short hand notation we will abbreviate the partition function by $`Z=_ie^{\beta _i}`$, taking into account that $`_i`$ contains the information about the chemical potential, the contributions of the quarks and gluons.

The observables $`<𝑶>`$ are calculated via

$`<𝑶>={\displaystyle \frac{_i𝑶e^{\beta _i}}{Z_a}},`$ (7)
where the index $`a`$ denotes the color configurations, i.e; $`a`$=(0,0) when only color zero states are considered and $`a=c`$ when also states with definite color are allowed. This distinction is needed to investigate the phase where color confinement is effective and the phase where color is allowed over a wide area of space. The quantities to calculate are the internal energy ($`<E>`$), heat capacity ($`<C>`$), average baryon number ($`<B>`$), strangeness ($`<s>`$) and the expectation value of different particle species ($`<n_k>`$), where $`k`$ refers to the quantum numbers of a particular particle and $`_k<n_k>=1`$, with the sum over all possible quantum numbers.

The particle expectation values have a simple expression because they select one of the eigenvalues at the time, thus if the state of a given particle is denoted by ”$`i`$” and $`E_i`$ is its energy, the particle expectation value is given by

$`<n_i>`$ $`=`$ $`{\displaystyle \frac{e^{\beta (E_i\mu _BB\mu _ss)}}{Z_a}},`$ (8)
(where we have used the value $`\mu _I=0`$, for the isospin chemical potential ).

At this point we have to make an assumption upon the volume considered. The whole reaction volume can be divided in elementary volumes, and we assume that the elementary volume ($`V_{el}`$) is of the size of a hadron, corresponding to a sphere with a radius of the order of 1 fm. Later on we shall show that this choice is reasonable, as seen from the calculated thermodynamic properties of the whole system. Another assumption is related to the interaction, which does not take into account confinement. We shall discuss two scenarios, namely: a) no additional interaction related to color is taken into account for temperatures above a critical (de-confinement) value, and b) confinement is operative for temperatures below the critical de-confinement temperature. In the regime (a) the lowest state with color (1,0) lies at the energy $`2\omega _f`$, and it corresponds to put one quark in the upper fermion level. Although it is a possible configuration, the Hamiltonian (3) cannot act upon it. This is consistent with the fact that above a certain temperature, $`T_c`$ (de-confinement temperature), only color non-singlet states are allowed. The actual value of $`T_c`$ will then give us an idea about the regime where hadronization is operative. In the real world hadronization, i.e. confinement, should set in below a critical temperature, as a true phase transition. In our model this will be signaled by a sharp transition from a state where color non-singlet states are still allowed ($`T>T_c`$) and a state where confinement is effective ($`TT_c`$).

## 4 Description of the high temperature regime: the QGP

.

We shall first discuss the case where no additional color interaction is taken into account. The states with energy $`E_i`$ for a given flavor, spin and parity are obtained from the diagonalization of the model Hamiltonian, now including flavor mixing and the corrections due to the Gel’man-Okubo mass formula for the two lowest meson nonets (one with spin zero and the other with spin 1). In Figure 1 we show the internal energy as a function of the temperature $`T`$, with and without interactions. The results shown by a dashed line have been obtained by calculating the internal energy in the fermion space without interactions. In this case the Hamiltonian has a simple image in the fermion space and the calculation of the partition function can be performed exactly. The dotted line shows the results obtained by working with the boson mapping, and by enforcing the corresponding cut-offs in the maximal number of bosons as explained above and in (I). The internal energy is a good indicator of the number of active states in the Hilbert space. Note that the results shown by both curves, the doted and dashed lines, practically coincide. This suggests that the number of active states is nearly the same in the boson and fermion spaces for a wide range of temperatures. This does not imply that all un-physical states have disappeared but rather that the approximate method of cutting un-physical states works reasonable well. The curve shown by a solid line, in the same Figure 1, gives the internal energy obtained from the calculation performed in the boson space and in presence of interactions. Although the curve does not show a clear phase transition of first order (e.g: a sharp increase of the energy in a narrow interval around $`T_c`$ ) the behavior around $`T_c=0.170`$ GeV is pretty suggestive of it. We have interpreted the observed smearing-out of the curve as follows: for $`T=0`$ the vacuum state is dominated by pairs of the type  and in this channel a quantum phase transition  does indeed take place. By this we mean that the pairs  are effectively blocked at high temperature. The other channels contribute less significantly to the ground state (see (I)) and interact weakly than the () channels, therefore, they remain in a perturbative regime. Thus, as the temperature increases, an approximate first order phase transition takes place in the channel , a mechanism similar to the one shown in Ref.  for the  channel, while the other configurations remain un-affected. The superposition of these two mechanism leads to the smearing out of the curve around $`T_c`$, as shown in Figure 1.

In Figure 2 the expectation value of the Casimir operator ($`C_2=<𝑪_\mathrm{𝟐}>`$) of color and its variation ($`\mathrm{\Delta }C_2=\sqrt{<𝑪_{\mathrm{𝟐}}^{}{}_{}{}^{2}><𝑪_\mathrm{𝟐}>^2}`$) are shown. The eigenvalue of the Casimir operator, for an irrep with color numbers $`(\lambda _C,\mu _C)`$, is given by $`C_2(\lambda _C,\mu _C)=\lambda _C^2`$ $`+`$ $`\lambda _C\mu _C`$ $`+`$ $`\mu _C^2`$ $`+`$ $`3(\lambda _C+\mu _C)`$. As a reference, for a color (1,0) irrep $`C_2=3`$ while the irrep (1,1) has $`C_2=9`$. We assume that $`\mu _B`$ and $`\mu _s`$ are zero. It is interesting to observe that, in the present model, the variation of the color is approximately symmetric around $`T=0.170`$ GeV. A possible interpretation is the following: at high energy the probability to have a color non-singlet state is large (the variation is not large enough to allow color singlet states) and a QGP is formed where color is effective over a wide range in space. From $`T=0.170`$ GeV on the probability to find a state in color (0,0) is significantly increased, since the variation is large enough to allow color singlet states. In lowering the temperature the variation is much larger than the average color and the whole QGP dissolves in droplets of color zero. Within the present model, these results, of the average color and its variation, are signals of the transition to the hadronic phase. Accordingly, we assume that it takes place near $`T_c=0.170`$ GeV, for $`\mu _B=\mu _s=0`$. We may now calculate the bag pressure and construct the $`T\mu _B`$ diagram. At $`T=0.170`$ GeV the pressure is determined via the expression $`p=\frac{Tln(Z)}{V_{el}}`$, where $`\mathrm{\Phi }=Tln(Z)`$ is the grand canonical partition function  and $`V_{el}`$ is the elementary volume $`V_{el}=\frac{4\pi }{3}r_{el}`$. For $`r_{el}`$=1 fm we obtain a bag pressure $`p^{\frac{1}{4}}`$ of about 0.17 GeV, which is in reasonable agreement with standard values. When the chemical potentials $`\mu _B`$ and $`\mu _s`$ are different from zero, the temperature dependence of the internal energy changes and also changes the value of the temperature for which the pressure is equal to the bag pressure $`p^{\frac{1}{4}}`$ $`=`$ $`0.18`$ GeV. The results are shown in Figure 3. Assuming that the local strangeness is $`<s>=0`$, we arrived at a functional relation between $`T`$, $`\mu _B`$ and $`\mu _s`$, i.e. $`f(T,\mu _B,\mu _s)=0`$, which fixes $`\mu _s`$ as a function of $`\mu _B`$. The results of this functional relation are displayed in Figure 4.

Once the chemical potential $`\mu _s`$ is adjusted, by using the results shown in Figures 3 and 4, the chemical potential $`\mu _B`$ and the transition temperature $`T_c`$ can be consistently determined.

Up to now we did not take into account an interaction which generates confinement. This has to be done by hand. One possibility is to assume that the transition from the QGP to the hadronic phase takes place within a very small range of temperatures around the critical temperature $`T_c`$. We require that the partition function above $`T_c`$ allows any color while for $`T<T_c`$ it contains only color zero states. Finally, chemical equilibrium connecting both phases, the QGP and the hadron gas, is understood. Figure 5 shows the results of the internal energy, without confinement (upper curve) and with confinement (lower curve). The solid line connecting both curves indicates the values for which confinement vanishes above $`T_c`$. As seen from the results, the transition is now of first order. Figure 6 shows the heat capacity calculated for the case without confinement (dotted line) and with confinement below $`T_c`$ (dashed line). The solid line interpolates between them, as in the case of the internal energy (see Figure 5).

The model can first be tested in the energy region below the transition temperature $`T_c`$, where the hadron gas should prevail. The confinement is effective and therefore we have to use the partition function $`Z_{a=(0,0)}`$ in the equations (7) and (8). We take, as an example, the measured total production rates of $`\pi ^+`$ at 10 GeV/A, as reported in the SIS-GSI experiment . The system considered was Au+Au and we assume that all particle participate, i.e. $`N_{part}`$=394. In Figure 2.3 of Ref.  a temperature of about $`T=0.13`$ GeV is reported. Assuming local strangeness conservation , $`<s>=0`$, we obtain a relation of $`\mu _s`$ versus $`T`$, which is depicted in Figure 7. From there we obtain, for the reported temperature, $`\mu _s0.128`$ GeV and via Figure 4 a value $`\mu _B0.55`$ GeV. In Figure 8 we show for a fixed value of $`\mu _B`$, physically acceptable at temperatures near $`T=0.13`$ GeV, the resulting total production rate of $`\pi ^+`$ ($`N_\pi `$) as a function in the temperature $`T`$. For $`T`$ = 0.13 GeV the production rate is approximately 180 pions $`\pi ^+`$, which is close to the value 160, which we have obtained by using Figure 2.3 of Ref. . The good qualitative agreement with the experiment demonstrates that the present model is able, indeed, to describe, approximately, observed QCD features.

We have also determined ratios of particle production and some absolute production rates. The particle production is calculated for temperatures just below $`T_c`$, where only color zero states are allowed. This implies that the partition function to use is $`Z_{(0,0)}`$. We then apply Eq. (8). Note that in the expression of the particle-production ratios the partition function cancels out and only the dependence on the mass of the particles and the chemical potential remains. Figure 9 shows results for some particle-production ratios for beam energies $`\sqrt{s}=130`$A GeV. The experimental values are taken from Ref. , based on the experiment described in Ref.  (see also ). For baryons, only the ratios of particle and anti-particle production are shown because these expressions are independent of the mass of the baryon. As noted in (I) the masses of the baryons are not well reproduced because they are considered as consisting of three idealized fermions on top of the meson sea. The interaction to the meson sea is not taken into account yet, but indications about how to do it are given in (I).

The central value of the $`K^{}/K^+`$ production ratio, shown in Figure 9, was reproduced with the values $`\mu _s=0.012`$ GeV and $`\mu _B=0.044`$ GeV. The other ratios are predicted by the model. Considering the simplicity of the model, the ratios are found to be in a reasonable agreement with data.

In order to obtain the total yields for kaons and for the $`\pi ^+`$ pions it is necessary to introduce further assumptions about the size of the QGP. The baryon density is given by $`\frac{<B>}{V_{el}}`$, where $`V_{el}`$ is the size of the representative volume, as explained before. In order to conserve, on the average, the baryon number we multiply the baryon density by the total volume and require that it must be equal to the total baryon number, given by the number of participants $`N_{part}`$. This leads to the total volume

$$V_{tot}=\frac{N_{part}V_{el}}{<B>_a},$$
(9)
where the index $`a`$ refers to, as in the partition function, color ($`a=c`$) when the average value is calculated in the QGP and $`a=(0,0)`$ when it is calculated in the hadron gas. In the QGP the average value $`<B>_c`$ before the transition is smaller than the average value $`<B>_{(0,0)}`$ in the hadron gas after the transition. This is due to the small value of $`Z_{(0,0)}`$ at $`T_c`$, since many other possible color states are excluded from it which do contribute to $`Z_c`$. As a consequence, for $`T=T_c`$ the volume of the QGP phase, as a function of $`N_{part}`$, is much smaller than the one in the hadron gas phase. Assuming a sphere, the radius of the QGP phase is about 8 fm, and it changes to about 20 fm after the transition. This implies volumes of the order of approximately 2 $`10^3`$ fm<sup>3</sup> and 3.4 $`10^4`$ fm<sup>3</sup>, respectively.

This transition is, as pointed out earlier, assumed to take place suddenly at $`T_c`$ (probably it should be smeared out but we cannot describe it with the present model). This implies that within the scenario assumed there is a rapid expansion of the volume caused by the transition from the QGP to the hadron phase, which should be observed as a large outward motion. One possible interpretation is that most of the pions are produced during the transition, liberating energy and provoking a rapid expansion of the system. The energy gained is represented by the jump between the lower and upper curves of Figure 5, but its origin cannot be explained by the present model, where confinement was shifted by hand.

Figure 10 shows the total pion yield as a function of the temperature $`T`$, corresponding to the Au+Au collision. The upper curve describes the total yield when all nucleons participate, while for the lower one we have taken $`N_{part}=250`$. This value agrees better with the experiment, as seen from the results, and it means that in the collision about 250 nucleons participate in the QGP. In Figure 11 the total kaon production rate is displayed. The upper curve corresponds to $`K^+`$ and the lower one to $`K^{}`$ absolute production rates, respectively. In both cases $`N_{part}=250`$ was used, the same value used previously in the calculation of the pion yield. The ratio of the curves was already adjusted at the point corresponding to the $`K^+/K^{}`$ ratio. The absolute production rate and the shape of the curve, however, is a prediction of the model (as far as we can talk about ”predictions” within this toy model). Considering the simplicity of the model it is surprising that the absolute production rate is well reproduced. This feature is common to other thermodynamical descriptions of the transition from the QGP to the hadron-gas .

Finally, in Figure 12, we show the calculated expectation values of the number of quark and gluon pairs as a function of the temperature T. At $`T=0`$ GeV the results correspond to the fractions of gluon pairs and fermion (quark-antiquark) pairs in the physical vacuum state. At high temperatures the gluon part increases and takes over the fermion part, which shows saturation. However, at the temperatures of interest, i.e. around the point of the phase transition $`T_c0.16`$ GeV, the gluon number is still suppressed with respect to the fermion pair number. This might be in favor of the ALCOR model  which supposes a suppression of gluons in the QGP and takes only constituent quarks and antiquarks into account. Note, that at $`T=0.170`$ GeV still a sensible amount of gluon pairs are present.

## 5 Conclusions

We have presented a toy model of QCD. The model is described in (I) and in this paper we have focused on the thermodynamic properties, at equilibrium, emerging from the model. We have calculated the partition function with and without color, and studied the temperature dependence of some observables, like the internal energy, the heat capacity, and the production rates of particles. The parameters of the model were determined in (I), adjusting the meson spectrum at low energy. Without further parameters the internal energy, the heat capacity and some particle ratios were determined, as explained in the text.

We have applied the model to the case of the Au+Au collision at 10 GeV/A  and shown that it can reproduce qualitatively the absolute production rate of $`\pi ^+`$. At this energy the QGP has not yet formed and, therefore, the results show that the model can be applied to study schematically the thermodynamics of a hadron gas.

Next, we have applied the model to energies where one assumes that the QGP has been already formed. The absolute production rate of $`\pi ^+`$ and kaons were calculated, just below the transition temperature, by taking the number of participant nucleons ($`N_{part}`$) as an input. The agreement between calculated and experimental values was found to be satisfactory. Also, the resulting production rate was described reasonable well, once the $`\mu _s`$ chemical potential was fixed to yield the correct (observed central value) $`\frac{K^+}{K^{}}`$ ratio. Some mass-independent baryon-antibaryon ratios, were qualitatively reproduced by the model predictions.

This demonstrates that the model is able to describe the general trend of QCD, in the finite temperature domain, and the transition to and from the quark gluon plasma.

## 6 Acknowledgment

We acknowledge financial support through the CONACyT-CONICET agreement under the project name Algebraic Methods in Nuclear and Subnuclear Physics and from CONACyT project number 32729-E. (S.J.) acknowledges financial support from the Deutscher Akademischer Austauschdienst (DAAD) and SRE, (S.L) acknowledges financial support from DGEP-UNAM. Financial help from DGAPA, project number IN119002, is also acknowledged.


text:
# Exact computation of the mean velocity, molecular diffusivity and dispersivity of a particle moving on a periodic lattice

## I Introduction

Lattice Monte Carlo simulation constitutes a powerful method for computing the average transport rates of a particle moving through a fixed network of obstacles by modeling the continuous transport process as a discrete random walk, biased or unbiased, across a lattice. At the heart of this technique lies a master equation van Kampen (1981) which quantifies the probability of jumping from one point on the lattice to an adjacent point in a given period of time. Owing to the stochastic (probabilistic) nature of the master equation, computer simulations based upon lattice Monte Carlo techniques require numerous realizations of the system in order to yield statistically significant results. Moreover, since average transport rates are valid only for long times, each simulation must be sufficiently long in order to satisfy this criteria.

Gel electrophoresis is one particularly relevant application of lattice Monte Carlo techniques, where the fibers of the gel are modeled as obstacles in the network which “reject” solute jumps onto them as a result of hard-core repulsion forces. In a recent series of publications Slater and Guo (1996a, b); Slater and Treurniet (1997); Mercier and Slater (1998); Labrie et al. (2000); Boileau and Slater (2001); Mercier and Slater (2001); Mercier et al. (2001); Gauthier and Slater (2002), Slater and coworkers invoked lattice Monte Carlo-type schemes to compute the effective solute mobility, $`\overline{𝐌}^{}`$, for transport in spatially periodic model gels in the “Ogston-sieving” regime of gel electrophoresis. The latter works are unique in that they do not rely upon computer simulations at all. Rather, analytical solutions are obtained for the solute mobility, based upon the fact that, for periodic lattices, unit-cell probabilities can be computed analytically in the long-time (steady-state) limit. As a consequence, Slater and coworkers have been able to consider a vast array of periodic model gels, both 2-D and 3-D, thereby critically assessing the validity of existing gel electrophoresis theories and investigating the effects of gel structure on the resultant average solute mobility. Results obtained from these types of calculations are presumably exact, at least in an asymptotic sense, for perfectly periodic structures, such as those which may be realized in microfluidic environments Chou et al. (1999); Cabodi et al. (2002); Doyle et al. (2002). Results for pseudo-random structures, which may be more representative of real gels, can be obtained by repeating the periodic calculation scheme numerous times in conjunction with random samplings of possible obstacle configurations (using the spatial periodicity to mimic macroscopic-sized gels), and then averaging the ensemble of results Slater and Treurniet (1997); Koplik (1982). The efficiency of analytical methods renders such calculations very feasible.

Although the aforementioned series of papers Slater and Guo (1996a, b); Slater and Treurniet (1997); Mercier and Slater (1998); Labrie et al. (2000); Boileau and Slater (2001); Mercier and Slater (2001); Mercier et al. (2001); Gauthier and Slater (2002) has made significant progress in addressing the issue of the average solute mobility (or, equivalently, the mean solute velocity, $`\overline{𝐔}^{}`$), much less progress has been made towards calculating solute dispersion arising as a result of the stochastic nature of the transport on the lattice. In the context of separation sciences Giddings (1991), it is essential to quantify both the dispersion and the average mobility, since the relative mobility of the particles being separated determines the degree of the separation, while the dispersion governs the separation sharpness. Moreover, while the local solute mobility dyadic, $`𝐌`$, and diffusivity dyadic, $`𝐃`$, are related at each point in space by the Stokes-Einstein equation 111 In free solution, these dyadics are isotropic and the solute mobility reduces to its Stokes mobility. ,

$$𝐃=𝐌kT,$$
(1)
with $`kT`$ the Boltzmann factor, the comparable Nernst-Einstein relationship between the dispersivity dyadic, $`\overline{𝐃}^{}`$, and the average mobility dyadic, $`\overline{𝐌}^{}`$,

$$\overline{𝐃}^{}=\overline{𝐌}^{}kT,$$
(2)
is only valid in the limit of vanishingly small applied fields.

Mercier et al. Mercier et al. (1999) exploited this property of eq. (2) to successively compute the components of $`\overline{𝐃}^{}`$ in the zero-field limit 222The zero-field limit dispersion coefficient is often referred to as the diffusion coefficient for the system, which should not be confused with the molecular diffusivity of the solute. For clarity, we will only use the term diffusion in the context of local motion and the term dispersion (or effective diffusion) in the context of averaged global motion. . Explicitly, an analytical scheme Slater and Guo (1996a) was invoked to compute the different mean velocity vectors, say, $`\overline{𝐔}_x^{}`$ and $`\overline{𝐔}_y^{}`$, corresponding to infinitesimally small forces being applied in the subscripted spatial directions. The average mobility was then obtained by using the linear relationship,

$$\overline{𝐔}^{}=\overline{𝐌}^{}𝐅,$$
(3)
in conjunction with each of the mean velocities. The resulting average mobility tensor proves to be independent of field strength, whereupon eq. (2) immediately furnishes the dispersivity dyadic, $`\overline{𝐃}^{}`$. The resulting dispersivity was shown to agree exactly with that obtained by first-passage time analysis Mercier et al. (1999). Given the relative ease of the mobility calculation compared to its first-passage time counterpart, Mercier et al. advocate the use of the former.

While calculating the dispersion (effective diffusion) coefficient for an unbiased random walk on the lattice is clearly a non-trivial endeavor, this is the limiting case of the much broader problem of dispersion occurring during biased random walks, the bias being induced by the presence of finite electric fields (or some other external force) 333 By an unbiased walk, we mean that all jumps are equally likely at a given lattice site, but where some of the jumps may be rejected by the presence of the obstacles. In contrast, for a biased random walk, the likelihood of making a jump in a particular direction is biased by the presence of the imposed force, resulting in a preferred direction of motion. Jumps onto sites occupied by the obstacles are still rejected in the biased random walk.. Keller et al. Keller et al. (2002) recently proposed a method for computing the dispersivity dyadic from biased lattice walks in the low-field limit, based upon what is essentially a volume-averaging scheme. In §VI.2, we will demonstrate that, in contrast to the present theory, the latter volume-averaged scheme does not furnish the correct effective diffusivity for the anistropic array discussed above Mercier et al. (1999). While this does not address the underlying validity (or lack thereof) of the volume-averaged approach, it does bring into question its applicability to the present class of problems. Consequently, no simple method exists in the literature for computing solute dispersion from a generic lattice model, and it has been noted Gauthier and Slater (2002) that such a calculation scheme would be extremely useful for analyzing trapping-based separations.

The present contribution addresses this need by developing an asymptotically exact, analytical scheme for computing the mean velocity vector, $`\overline{𝐔}^{}`$, and dispersivity dyadic, $`\overline{𝐃}^{}`$, from lattice Monte Carlo models, by way of discrete generalized Taylor-Aris dispersion theory for spatially periodic networks Adler and Brenner (1984); Dorfman and Brenner (2002). The latter theory was initially developed as a means to compute $`\overline{𝐔}^{}`$ and $`\overline{𝐃}^{}`$ from lumped-parameter models of transport in model porous media Adler and Brenner (1984) and microfluidic networks Dorfman and Brenner (2002). We will show that the lattice random walk has much in common with the existing lumped-parameter transport problem, whereupon the moment scheme invoked previously Adler and Brenner (1984); Dorfman and Brenner (2002) may be readily employed here to homogenize the master equation governing transport on the periodic lattice. Consideration of the obstacle-free case furnishes generic relationships between the jump time, transition probabilities and the lattice spacing. These generic relationships between the parameters are used in an illustrative example of a solute moving on a simple periodic lattice in the presence of a finite external field. The theory simplifies substantially in the absence of a field, leading to a compact scheme for computing the effective diffusion coefficient. The simplified theory is invoked to compute the zero-field dispersion coefficients for the asymmetric lattice problem originally considered by Mercier et al. Mercier et al. (1999), demonstrating that the present theory successfully reproduces their results.

## II Biased Random Walk on a Periodic Lattice

### II.1 Geometry

Consider a periodic lattice of the type depicted in Fig. 1. The available sites of the lattice are indicated by white boxes and the unavailable sites, corresponding to the obstacles in the array, are indicated by black boxes. Transport between lattice sites occurs by discrete “jumps” between them. The spatial periodicity of the lattice is reflected by the presence of a repetitive unit cell, which reproduces the entire unbounded medium upon translation by its base lattice vectors $`(𝐥_1,𝐥_2,𝐥_3)`$. Each cell is characterized by its respective discrete location vector, $`𝐈=(I_1,I_2,I_3)`$, with the origin arbitrarily chosen to be located at the cell $`𝐈_\mathrm{𝟎}=\mathrm{𝟎}`$. The position of the locator point, say, the centroid, of cell $`𝐈`$ is given by the discrete vector $`𝐑_𝐈=I_1𝐥_\mathrm{𝟏}+I_2𝐥_\mathrm{𝟐}+I_3𝐥_\mathrm{𝟑}`$.

The unbounded, periodic lattice may be represented by a series of three graphs (although only the unit-cell based local graph, $`\mathrm{\Gamma }_l`$, will be necessary for computing $`\overline{𝐔}^{}`$ and $`\overline{𝐃}^{}`$), where the nodes of the graph correspond to the lattice sites and the edges of the graph correspond to the possible jumps between sites. The specifics of the graphical construction are discussed at length in Ref. Dorfman and Brenner, 2002, whereupon the following exposition will be brief. Each lattice site (node) is assigned a discrete location $`(𝐈,i)`$, where $`i`$ represents the intracellular position (i.e. $`a`$-$`g`$ in Fig. 1). As a consequence of the spatial periodicity of the network, it only proves necessary to consider edges of the type indicated in Fig. 1. Let $`\mathrm{\Omega }^+(i)`$ be the subset of the edges $`j`$ entering node $`(𝐈,i)`$ from node $`(𝐈^{},i^{})`$, and $`\mathrm{\Omega }^{}(i)`$ be those edges exiting node $`(𝐈,i)`$ and entering node $`(𝐈^{},i^{})`$ 444 For those edges contained entirely in unit cell $`𝐈`$, $`𝐈=𝐈^{}`$. . The global graph, $`\mathrm{\Gamma }_g`$, contains all of the infinitely many nodes and all of the edges connecting them, equipollent with the unbounded lattice. As an intermediate step, we construct the basic graph, $`\mathrm{\Gamma }_b`$, depicted in Fig. 1, which includes: (i) all of the nodes in the unit cell; (ii) all those edges contained wholly within the unit cell; (iii) those edges entering the unit cell; and (iv) their associated homologous vertices, the latter indicated with the prime affix in Fig. 1.

A macroscopic jump vector,

$$𝐑(j)=𝐑_𝐈𝐑_𝐈^{},$$
(4)
is assigned to each of the edges, this vector being identically zero for those edges which do not cross the boundary of the unit cell. This permits one to distinguish between edges connecting, say, $`e^{}b`$ and $`eb`$ in the basic graph, since their respective macroscopic jump vectors are $`𝐑=𝐥_y`$ and $`𝐑=\mathrm{𝟎}`$. Use of the macroscopic jump vectors prevents some (potentially significant) geometrical simplifications achieved in the lattice models proposed by Slater and coworkers Slater and Guo (1996a, b), since “identical” sites in their model now typically differ with respect to the macroscopic jump vectors in $`\mathrm{\Omega }^+(i)`$. Since the velocity may be computed without use of the these vectors \[cf. eq. (17)\], their scheme Slater and Guo (1996b) for computing $`\overline{𝐔}^{}`$ is more efficient for certain geometries. However, the primary contribution here is a a scheme for computing $`\overline{𝐃}^{}`$ \[cf. eq. (18)\], which requires using the macroscopic jump vectors to compute the moments of the probability \[cf. eq. (9)\].

The local graph, $`\mathrm{\Gamma }_l`$, depicted in Fig. 2, is constructed by combining the homologous vertices on $`\mathrm{\Gamma }_b`$ and contracting those edges between them. The dashed lines in Fig. 2 correspond to edges which cross the boundary of the unit cell in the basic graph (and thus have non-zero macroscopic jump vectors). Reference to the discrete position vector $`𝐈`$ is no longer necessary in the local graph, since the periodic structure of the medium has been embedded in the macroscopic jump vectors.

### II.2 Master Equation

Consider the probability $`P(𝐈,i,N|𝐈_0,i_0,0)`$ that the particle will be located at site $`(𝐈,i)`$ on $`\mathrm{\Gamma }_g`$ at step $`N`$, given its initial location at site $`(𝐈_0,i_0)`$. In the usual manner Dorfman and Brenner (2002), the final results only depend upon the displacement, $`𝐈𝐈_0`$, from the arbitrarily positioned origin. In light of the prior choice $`𝐈_\mathrm{𝟎}=\mathrm{𝟎}`$, the probability adopts the canonical form $`PP(𝐈,i,N|i_0)`$.

During each time step of duration $`\tau `$, a particle located at a node adjacent to node $`i`$, say, node $`i^{}`$, may enter node $`i`$ via a jump through an edge in $`\mathrm{\Omega }^+(i)`$. If the particle is located at node $`i`$, it may either remain there, say, due to a rejected jump caused by the presence of an obstacle, or exit through one of the edges contained in $`\mathrm{\Omega }^{}(i)`$. For each edge $`j`$, denote the transition probability $`w(j)`$ that a particle located at the initial vertex of edge $`j`$ will chose to move through that edge during time $`\tau `$, and for each node $`i`$, let $`w(i)`$ denote the probability that the particle will choose to remain at node $`i`$ during time $`\tau `$. Since a particle located at node $`i`$ must either remain there or exit therefrom at each time step, there exists the conservation relationship 555 In the context of the original discrete Taylor-Aris dispersion theory Dorfman and Brenner (2002), it is possible to account for the rejected jumps by adding loops to the graph which return the particle to its original position. This alters the form of the master equation \[cf. eq. (8)\], requiring an additional sum over the outlet edges in $`\mathrm{\Omega }^{}(i)`$. Since the underlying logic of the lattice random walk is that the particle must make a jump at each time step (albeit with the possibility that this jump will be rejected), the resulting calculation scheme is greatly simplified by requiring a unitary probability \[cf. eq. (5)\] of jumping and accounting for rejected jumps by the factor $`w(i)`$. :

$$w(i)+\underset{j\mathrm{\Omega }^{}(i)}{}w(j)=1(i).$$
(5)
The evolution of $`P`$ is governed by the following master equation quantifying the change in probability at each node on the global graph, $`\mathrm{\Gamma }_g`$, between time step $`N`$ and $`N+1`$:

$`P(𝐈,i,N+1)P(𝐈,i,N)=`$ $`{\displaystyle \underset{j\mathrm{\Omega }^+(i)}{}}w(j)P(𝐈^{},i^{},N)`$ (6)
$`[1w(i)]P(𝐈,i,N)`$
$`+\delta (𝐈,\mathrm{𝟎})\delta (i,i_0)\times `$
$`\times \delta (N,0),`$
where the dependence upon the initial condition $`i_0`$ has been suppressed, and the $`\delta (i,j)`$ are Kronecker delta functions; $`\delta (i,j)=1`$ if $`i=j`$, $`\delta (i,j)=0`$ otherwise.

A rate equation may be obtained by dividing eq. (6) by the jump time, $`\tau `$. Interest here is focused upon the asymptotic, long-time behavior of $`P`$ after many jumps ($`N1`$), whence, to a good approximation 666 It is expected that moving from a finite-difference approximation to a continuous derivative reflects the physical process being modeled, since the lattice model merely represents a discrete approximation to a continuous process.,

$$\frac{P(𝐈,i,N+1)P(𝐈,i,N)}{\tau }\frac{dP(𝐈,i,N)}{dt}.$$
(7)
The time, $`tN\tau `$, is viewed as an essentially continuous variable for the case $`N1`$, which is equivalent to the long-time limit in conventional Taylor-Aris dispersion theory Dorfman and Brenner (2002). Use of the latter in eq. (6) furnishes the following differential algebraic equation

$`{\displaystyle \frac{dP(𝐈,i,N)}{dt}}=`$ $`{\displaystyle \underset{j\mathrm{\Omega }^+(i)}{}}\left[{\displaystyle \frac{w(j)}{\tau }}\right]P(𝐈^{},i^{},N)`$ (8)
$`\left[{\displaystyle \frac{1w(i)}{\tau }}\right]P(𝐈,i,N)`$
$`+\tau ^1\delta (𝐈)\delta (i,i_0)\delta (N).`$
## III Moment Scheme

The master equation (8) possesses a mathematical structure similar to the governing equation appearing in Ref. Dorfman and Brenner, 2002, albeit with a significantly different physical interpretation. Consequently, the moment-matching scheme employed previously may be invoked here to compute the mean velocity vector, $`\overline{𝐔}^{}`$, and dispersivity dyadic, $`\overline{𝐃}^{}`$, for the lattice transport problem defined by eq. (8). The scheme is outlined in what follows, and the reader is referred to Ref. Dorfman and Brenner, 2002 for further details.

Define the local moment of the probability as the $`m`$-adic,

$$𝐏_m(i,N|i_0)\stackrel{\text{def.}}{=}\underset{𝐈}{}(𝐑_𝐈)^mP(𝐈,i,N|i_0),$$
(9)
where $`(𝐑_𝐈)^m𝐑_𝐈𝐑_𝐈\mathrm{}𝐑_𝐈`$ ($`m`$ times) and the sum over $`𝐈`$ represent the sum over all cells. (The spatial location of cell $`𝐈_\mathrm{𝟎}=\mathrm{𝟎}`$ is $`𝐑_{𝐈_0}=\mathrm{𝟎}`$.) The total moments of the probability are the $`m`$-adics,

$$𝐌_m(N|i_0)\stackrel{\text{def.}}{=}\underset{iV\mathrm{\Gamma }_l}{}𝐏_m(i,N|i_0),$$
(10)
where $`V\mathrm{\Gamma }_l`$ are the vertices of the local graph. The zeroth moment is trivial,

$$M_0=1,$$
(11)
which reflects the conservation of the probability on the lattice. The next two higher-order total moments are expected to grow linearly with time, and possess the respective forms:

$`{\displaystyle \frac{d𝐌_1}{dt}}\overline{𝐔}^{},`$ (12)
$`{\displaystyle \frac{d}{dt}}\left(𝐌_2𝐌_1𝐌_1\right)2\overline{𝐃}^{}.`$ (13)
We briefly outline here the moment-matching scheme Dorfman and Brenner (2002) for computing $`\overline{𝐔}^{}`$ and $`\overline{𝐃}^{}`$: (i) form the equations governing the local moments from eqs. (8) and (9); (ii) compute the asymptotic solutions to the first two local moments; (iii) form asymptotic total moments from the latter result and eq. (10); and (iv) match the asymptotic moments of the exact solution to eqs. (12)-(13). Only the results of the moment scheme will be presented, since the mathematical manipulations required to arrive at these results are lengthy and follow a procedure identical to that appearing elsewhere Dorfman and Brenner (2002).

### III.1 Mean Velocity Vector

The mean velocity vector is computed by forming the sum

$$\overline{𝐔}^{}=\underset{jE\mathrm{\Gamma }_l,j\mathrm{\Omega }^+}{}\left[\frac{w(j)}{\tau }\right]𝐑(j)P_0^{\mathrm{}}(i^{}),$$
(14)
where $`E\mathrm{\Gamma }_l`$ are the edges of the local graph, and the index $`i^{}`$ refers to the initial vertex of edge $`j`$. The vertex field $`P_0^{\mathrm{}}(i)`$, corresponding to the asymptotic, steady-state probability of locating the particle at node $`i`$ in the unit cell, is the solution of the linear set of equations,

$$\underset{j\mathrm{\Omega }^+(i)}{}w(j)P_0^{\mathrm{}}(i^{})[1w(i)]P_0^{\mathrm{}}(i)=0.$$
(15)
The latter are not linearly independent, and must be supplemented by the normalization condition

$$\underset{iV\mathrm{\Gamma }_l}{}P_0^{\mathrm{}}(i)=1.$$
(16)
The present method of computing $`\overline{𝐔}^{}`$ differs somewhat from that employed by Slater and coworkers (see, for example, Ref. Slater and Guo, 1996a). In the latter, the mean velocity is computed by first solving eqs. (15)-(16) for $`P_0^{\mathrm{}}(i)`$. The velocity vector for each lattice site, say, $`𝐯(i)`$, is then defined as the difference in transition probability between taking a forward step and a backwards step at site $`i`$, each step being weighted by its respective displacement (i.e. zero displacement for rejected jumps). The mean velocity is then computed by the sum

$$\overline{𝐔}^{}=\underset{iV\mathrm{\Gamma }_l}{}𝐯(i)P_0^{\mathrm{}}(i).$$
(17)
Equations (14) and (17) furnish identical values for $`\overline{𝐔}^{}`$, since the difference between these two methods is simply a manifestation of Gauss’ divergence theorem in a discrete sense; eq. (14) is the “surface” term and eq. (17) is the “volume” term.

### III.2 Dispersivity Dyadic

The dispersivity dyadic is computed by the edge-sum

$$\overline{𝐃}^{}=\frac{1}{2}\underset{jE\mathrm{\Gamma }_l,j\mathrm{\Omega }^+}{}\left[\frac{w(j)}{\tau }\right]P_0^{\mathrm{}}(i^{})\stackrel{~}{𝐛}(j)\stackrel{~}{𝐛}(j),$$
(18)
wherein

$$\stackrel{~}{𝐛}(j)\stackrel{\text{def.}}{=}𝐑(j)𝐁(i)+𝐁(i^{}),j=\{i^{},i\}$$
(19)
for an edge directed from $`i^{}`$ to $`i`$. The node-based vectors $`𝐁(i)`$ are the solutions of

$`{\displaystyle \underset{j\mathrm{\Omega }^+(i)}{}}\left[{\displaystyle \frac{w(j)}{\tau }}\right]P_0^{\mathrm{}}(i^{})[𝐁(i^{})+𝐑(j)]`$
$`\left[{\displaystyle \frac{1w(i)}{\tau }}\right]P_0^{\mathrm{}}(i)𝐁(i)`$
$`=P_0^{\mathrm{}}(i)\overline{𝐔}^{}.`$ (20)
Equations (20) only define the $`𝐁`$ vectors to within an arbitrary, additive constant vector. Consequently, this degree of freedom may be employed to set the $`𝐁`$ vector at one of the nodes, say $`i^{}`$, equal to zero, $`𝐁(i^{})=\mathrm{𝟎}`$.

## IV Obstacle-Free Case

Up to this point, no restrictions were placed upon the transition probabilities, $`w(j)`$, the lattice spacing, $`l`$, or the jump time, $`\tau `$, aside from the conservation statement (5), which ensures that the particle makes a decision (either to jump to a new site or remain at the present site) during each time step. However, by considering the trivial case of transport on an obstacle-free lattice, additional relationships may be derived between the latter parameters, at least for the case of square lattices.

Attention is restricted here to a two dimensional lattice, whose basic graph is depicted in Fig. 3, since generalizations to other dimensions are trivial. In the present circumstances, it only proves necessary to consider the case where the force is directed in a single direction on the lattice, say, the $`x`$-direction, since the spatial orientation of the lattice is aribtrary. The unit cell consists of a single node with four entering edges, each of the latter possessing the transition probability $`w_{+x}`$, $`w_x`$ or $`w_y`$, depending on the edge orientation.

The governing equations (15)-(16) and (20) are satisfied here by the trivial solutions $`P_0^{\mathrm{}}=1`$ and $`𝐁=\mathrm{𝟎}`$ \[whereupon $`\stackrel{~}{𝐛}(j)=𝐑(j)`$\]. From eq. (14), the mean velocity adopts the form

$$\overline{𝐔}^{}=𝐢_x\left(w_{+x}w_x\right)\frac{l}{\tau },$$
(21)
and, from eq. (18), the dispersivity dyadic is

$$\overline{𝐃}^{}=𝐢_x𝐢_x\left[\left(w_{+x}+w_x\right)\frac{l^2}{2\tau }\right]+𝐢_y𝐢_y\left(w_y\frac{l^2}{\tau }\right).$$
(22)
For this trivial problem, though, the proper results are known a priori, namely $`\overline{𝐔}^{}=M𝐅`$, where $`M`$ is the free solution mobility and $`𝐅=𝐢_xF`$ is the applied force, and $`\overline{𝐃}=𝐈D`$, where $`D`$ is the molecular diffusivity and $`𝐈`$ is the idemfactor. As a consequence, the usual velocity restriction Slater et al. (1986) is recovered,

$$\left(w_{+x}w_x\right)\frac{l}{\tau }=MF,$$
(23)
along with two restrictions governing the diffusive behavior,

$$\left(w_{+x}+w_x\right)\frac{l^2}{2\tau }=\frac{w_yl^2}{\tau }=D.$$
(24)
It is our contention that the latter restrictions (23)-(24) must be satisfied by any Monte Carlo algorith of this type, in order for the moments of the master equation to properly reflect the microscale physics. For example, Gauthier and Slater Gauthier and Slater (2002) noted that their high-field Monte Carlo algorithm, which furnishes the correct velocity (23), predicts that the obstacle-free diffusivity depends upon the field strength. Indeed, it is readily confirmed that their choices of transition probabilities (and $`\tau `$) fail to satisfy eq. (24). Consequently, the criteria of eq. (24) serves as a litmus test for newly proposed algorithms.

If the following normalization condition is enforced,

$$w_{+x}+w_x+2w_y=1,$$
(25)
then eqs. (23)-(25) constitute 4 equations for 5 unknowns. If the lattice spacing, $`l`$, is left as an ajustable parameter, then one recovers the usual “small-bias” algorithm (see, for example, Ref. Slater and Guo, 1996a)

$$w_{\pm x}=\frac{1\pm ϵ}{4},w_y=\frac{1}{4},\tau =\frac{l^2}{4D},$$
(26)
where $`ϵ`$ is the reduced field,

$$ϵ\stackrel{\text{def.}}{=}\frac{Fl}{2kT}.$$
(27)
Remarkably, the latter results were derived solely from considering the obstacle-free velocity and diffusion coefficient, along with the normalization condition (25), without ever explicitly stating that $`ϵ1`$. However, the transition probabilities are only sensible in the latter limit, or at least for the inequality $`ϵ<1`$, since the physical interpretation of a negative transition probability is suspect.

## V Transport on an Isotropic Lattice in a Finite External Field

In the following illustrative example, we analyze transport occurring in a simple isotropic system, depicted in Fig. 4, under the influence of a finite field. The lattice, consisting of square obstacles of size $`l`$ separated by a distance $`l`$, constitutes the simplest configuration for demonstrating the present scheme.

In the parlance of the present theory, and with the edge numbering depicted in Fig. 4, the transition probabilities are $`w(j)=w_{+x}`$ for $`(j=1,2)`$, $`w(j)=w_x`$ for $`(j=3,4)`$, and $`w(j)=w_y`$ otherwise. For now, the specific forms of the transition probabilities, $`w(j)`$, and the jump time, $`\tau `$, will be left unspecified, although subject to the generic physical restrictions (23)-(24). As a consequence of rejected jumps, the probabilities of remaining at the different lattice sites are $`w(a)=w_{+x}+w_x`$, $`w(b)=2w_y`$ and $`w(c)=0`$. The problem specification is completed by the macroscopic jump vectors,

$$𝐑(j)=2l\{\begin{array}{cc}\hfill 𝐢_x,& j=1,\hfill \\ \hfill 𝐢_x,& j=4,\hfill \\ \hfill 𝐢_y,& j=5,\hfill \\ \hfill 𝐢_y,& j=8,\hfill \\ \hfill \mathrm{𝟎},& \text{otherwise}.\hfill \end{array}$$
(28)
It is readily apparent that the solution $`P_0^{\mathrm{}}=1/3`$ satisfies eqs. (15)-(16), a consequence of the isotropy of the lattice. Substituting the relevant parameters into eq. (14) furnishes the mean velocity,

$$\overline{𝐔}^{}=\frac{2}{3}\left(w_{+x}w_x\right)\frac{l}{\tau }𝐢_x,$$
(29)
which, with use of eq. (23), adopts the form

$$\overline{𝐔}^{}=\frac{2}{3}M𝐅.$$
(30)
The latter result is invariant to the (arbitrary) choice of axis labels, whereupon it is recognized that isotropy of the lattice geometry renders the mean mobility tensor isotropic,

$$\overline{𝐌}^{}=𝐈\overline{M}^{},\overline{M}^{}=\frac{2}{3}M.$$
(31)
The effective mobility, $`\overline{M}^{}`$, is less than the free solution value, $`M`$, owing to the inability to make forward jumps on lattice site $`a`$. The latter observation, as well as the equivalence of eqs. (14) and (17), is rendered transparent by this considering the local velocity Slater and Guo (1996b) at each node,

$$𝐯(i)=\{\begin{array}{cc}M𝐅,\hfill & i=b,c,\hfill \\ \mathrm{𝟎},\hfill & i=a.\hfill \end{array}$$
(32)
Use of the latter and eq. (17) furnishes the result (30).

Since nodes $`a`$ and $`b`$ communicate via node $`c`$, the remaining calculations are greatly simplified by choosing $`𝐁(c)=\mathrm{𝟎}`$, whereupon eq. (20) directly furnishes the values of the remaining $`𝐁`$ vectors,

$`{\displaystyle \frac{𝐁(a)}{l}}`$ $`=`$ $`{\displaystyle \frac{2}{3}}\left({\displaystyle \frac{w_xw_{+x}}{2w_y}}\right)𝐢_x𝐢_y,`$ (33)
$`{\displaystyle \frac{𝐁(b)}{l}}`$ $`=`$ $`{\displaystyle \frac{2}{3}}\left({\displaystyle \frac{2w_{+x}+w_x}{w_x+w_x}}\right)𝐢_x,`$ (34)
which have been simplified with use of eq. (25).

Compute $`\stackrel{~}{𝐛}(j)`$ from eq. (19), with use of eqs. (47) and (74), and use the result in eq. (45) to ultimately arrive at the dispersivity dyadic

$$\frac{\overline{𝐃}^{}}{D}=\overline{D}_{xx}^{}𝐢_x𝐢_x+\overline{D}_{yy}^{}𝐢_y𝐢_y,$$
(35)
where the components of the dispersivity,

$`{\displaystyle \frac{\overline{D}_{xx}^{}}{D}}`$ $`=`$ $`{\displaystyle \frac{4}{27}}\left[\begin{array}{c}\frac{\left(2w_{+x}+w_x\right)^2+\left(2w_x+w_{+x}\right)^2}{\left(w_{+x}+w_x\right)^2}+\hfill \\ +\frac{\left(w_{+x}w_x\right)^2}{2w_y^2}\hfill \end{array}\right],`$ (38)
$`{\displaystyle \frac{\overline{D}_{yy}^{}}{D}}`$ $`=`$ $`{\displaystyle \frac{2}{3}},`$ (39)
have been rendered in terms of the molecular diffusivity, $`D`$, with use of eq. (24).

Clearly, the lateral dispersivity, $`\overline{D}_{yy}^{}`$, obeys the Nernst-Einstein relationship (2) for any choice of the transition parameters, as must be the case since there exists no bias in the lateral direction. In contrast, the axial dispersivity, $`\overline{D}_{xx}^{}`$, only obeys the Nernst-Einstein relationship in the limit where $`w=1/4(j)`$, i.e. for pure molecular diffusion. Importantly, the axial dispersion is also invariant to the arbitrary orientation of the $`x`$-axis. If one adopts the small-bias parameters (26), the axial dispersion reduces to the form

$$\frac{\overline{D}_{xx}^{}}{D}=\frac{2}{3}+\frac{10}{27}ϵ^2.$$
(40)
For finite field strengths, there exists a so-called “Taylor” contribution (or convective dispersivity) in the direction of the field, the latter depending quadratically upon $`ϵ`$, or, equivalently, the square of the Peclet number. The latter result agrees with conventional continuum analyses of the convective dispersivity.

## VI Unbiased Random Walk

### VI.1 Simplification of the General Theory

Significant simplifications are realized when $`w(j)=w=\text{const.}(j)`$, which is equivalent in its consequences to an unbiased random walk (or molecular diffusion) on the lattice, albeit with the possibility of rejected steps due to the presence of the obstacles. For a unit cell with $`\eta `$ available sites, the asymptotic probability adopts the form

$$P_0^{\mathrm{}}(i)=\eta ^1(i),$$
(41)
reflecting the fact that it is equally likely for the particle to be located at any of the available sites 777 When using the Nernst-Einstein relationship (2) to compute the effective diffusivity in the zero-field limit, the presence of an infinitesimally small force induces a bias in $`P_0^{\mathrm{}}`$. (As an illustrative example, consider the effect of curved field lines Mercier et al. (2001).) The latter bias appears here only in the $`𝐁`$-field, since the force is identically zero. . With use of eq. (5), it is easily shown that eq. (41) satisfies eqs. (15)-(16). Substitution of eq. (41) into (14) reveals that

$$\overline{𝐔}^{}=\left[\frac{w}{\eta \tau }\right]\underset{jE\mathrm{\Gamma }_l}{}𝐑(j).$$
(42)
Since $`E\mathrm{\Gamma }_l`$ only contains edges entering the unit cell, the lattice structure guarantees that each edge possessing the macroscopic jump vector $`𝐑(j)`$ can be paired with another edge whose macroscopic jump vector is $`𝐑(j)`$ (see, for example, Fig. 2). Consequently, the summation appearing in eq. (42) is identically zero, whereupon

$$\overline{𝐔}^{}=\mathrm{𝟎},$$
(43)
as must be the case for a pure diffusion problem.

The $`𝐁`$-equation (20) undergoes a similar simplification in the unbiased limit,

$$\underset{j\mathrm{\Omega }^+(i)}{}[𝐁(i^{})+𝐑(j)]\left[\frac{1w(i)}{w}\right]𝐁(i)=\mathrm{𝟎}.$$
(44)
The $`\stackrel{~}{𝐛}(j)`$ are still computed by eq. (19), whereupon $`\overline{𝐃}^{}`$ adopts the form

$$\overline{𝐃}^{}=\left[\frac{w}{2\eta \tau }\right]\underset{jE\mathrm{\Gamma }_l}{}\stackrel{~}{𝐛}(j)\stackrel{~}{𝐛}(j),$$
(45)
For certain geometries, the present scheme for computing the effective diffusivity is computationally more efficient than the alternative procedure Mercier et al. (1999) of evaluating $`\overline{𝐌}^{}`$ in the low-force limit and then using eq. (2) to compute $`\overline{𝐃}^{}`$. First, the coefficient matrix appearing in eq. (44) is invariant to spatial direction. Consequently, it is only necessary to specify the new solution vector in order to calculate a second (or third) component of $`𝐁`$. In contrast, the method of Ref. Mercier et al., 1999 requires computing a separate set of probabilities corresponding to the external force being oriented in each direction, each calculation (generally) involving a new coefficient matrix. While both methods require the same number of matrix inversions, the $`𝐁`$ matrices are one dimension smaller. However, we expect this slight edge to be offset by the more involved calculations required to compute $`\overline{𝐃}^{}`$ by eq. (45). Importantly, though, the present method corresponds exactly to the case of no applied force, which eliminates the need to obtain the perturbation to the uniform probability distribution in the low-force limit, although this can be done fairly efficiently Mercier and Slater (1998). The Nernst-Einstein equation scheme is preferable, however, when significant geometrical simplifications can be effected by eliminating the macroscopic jump vectors and making use of “identical” nodes. That being said, both the present technique and the Nernst-Einstein limit technique are markedly less complicated than first-passage time analyses, and either method serves as an efficient analytical alternative to conventional simulation techniques, as well as a check on complex numerical codes.

### VI.2 Pure Diffusion in an Asymmetric Array

As a final example, we compute the effective diffusivity of a particle moving on the lattice shown in Fig. 1. The effective diffusivity for this particular geometry was previously computed by Mercier et al. Mercier et al. (1999) using the techniques mentioned §I, and we demonstrate here that the present theory reproduces their results.

There are seven available locations in the unit cell, $`\eta =7`$, and the jump probability is $`w=1/4`$. The jump time, $`\tau `$, is given by the Brownian time scale of eq. (26). The probability of remaining on node $`i`$ is equal to the number of adjacent unavailable sites divided by four,

$$w(i)=\{\begin{array}{cc}0,\hfill & i=b,c,d,f,\hfill \\ 1/4,\hfill & i=e,g,\hfill \\ 1/2,\hfill & i=a.\hfill \end{array}$$
(46)
With the edge numbers depicted in Fig. 2, the macroscopic jump vectors are

$$𝐑(j)=2l\{\begin{array}{cc}\hfill 2𝐢_x,& j=1,\hfill \\ \hfill 2𝐢_x,& j=2,\hfill \\ \hfill 𝐢_y,& j=35,\hfill \\ \hfill 𝐢_y,& j=68,\hfill \\ \hfill \mathrm{𝟎},& \text{otherwise}.\hfill \end{array}$$
(47)
Substituting the relevant parameters into eq. (44) and choosing $`i^{}=d`$ (i.e. $`𝐁(d)=\mathrm{𝟎}`$),

$`\left[\begin{array}{cccccc}2& 1& 0& 0& 0& 0\\ 1& 4& 1& 2& 0& 0\\ 0& 1& 4& 0& 2& 0\\ 0& 2& 0& 3& 1& 0\\ 0& 0& 2& 1& 4& 1\\ 0& 0& 0& 0& 1& 3\end{array}\right]\left[\begin{array}{c}𝐁(a)\\ 𝐁(b)\\ 𝐁(c)\\ 𝐁(e)\\ 𝐁(f)\\ 𝐁(g)\end{array}\right]`$ (60)
$`=4l\left[\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]𝐢_x+2l\left[\begin{array}{c}0\\ 1\\ 1\\ 1\\ 1\\ 1\end{array}\right]𝐢_y.`$ (73)
The solution of this system of equations is

$$\left[\begin{array}{c}𝐁(a)\\ 𝐁(b)\\ 𝐁(c)\\ 𝐁(e)\\ 𝐁(f)\\ 𝐁(g)\end{array}\right]=\frac{l}{4}\left[\begin{array}{c}11\\ 6\\ 3\\ 5\\ 3\\ 1\end{array}\right]𝐢_x+l\left[\begin{array}{c}0\\ 0\\ 0\\ 1\\ 1\\ 1\end{array}\right]𝐢_y.$$
(74)
Compute $`\stackrel{~}{𝐛}(j)`$ from eq. (19), with use of eqs. (47) and (74), and use the result in eq. (45) to furnish the dispersivity dyadic

$$\overline{𝐃}^{}=\left(\frac{5}{28}𝐢_x𝐢_x+\frac{3}{14}𝐢_y𝐢_y\right)\frac{l^2}{\tau },$$
(75)
or, in terms of $`D`$,

$$\overline{𝐃}^{}=\left(\frac{5}{7}𝐢_x𝐢_x+\frac{6}{7}𝐢_y𝐢_y\right)D,$$
(76)
The latter result agrees exactly with Ref. Mercier et al., 1999. Although not shown here, the present calculation scheme also reproduces the low-field mobilities obtained previously Mercier et al. (1999) during the course of the Nernst-Einstein calculation, although knowledge of the latter mobilities proves unnecessary to compute $`\overline{𝐃}^{}`$.

Having furnished a third independent verification of the value of $`\overline{𝐃}^{}`$ for this array, we turn our attention to the volume-averaging theory Keller et al. (2002) and investigate its comparable result. The pure diffusion problem includes no bias in the transition probabilities and results in zero mean velocity, whereupon their generic result for the $`yy`$ component of the dispersivity adopts the form

$$\overline{D}_{yy}^{}=\mathrm{\Gamma }\left(1\underset{W}{}\widehat{n}_yg^y\right)v_yg^y,$$
(77)
where $`\mathrm{\Gamma }`$ is their dimensionless diffusivity, $`_W`$ is the summation over all sites adjacent to the obstacle, $`\widehat{n}_y`$ are outward pointing normal vectors from the obstacles, $``$ is the sum over all available sites, $`v_y`$ is the microscopic drift velocity, and $`g^y`$ is the volume-averaging quantity which plays a similar role to our $`𝐁`$ field. Refering back to Fig. 1, the quantity $`v_y=0`$, since the transport process only includes diffusion in the $`y`$-direction. Consequently, the second term in eq. (77) makes no contribution to the effective diffusivity. Moreover, only site $`a`$ has a non-zero value of $`\widehat{n}_y`$. Indeed, the latter site contains both a positive and negative unit vector, corresponding to the lower and upper obstacle sites. As such, the summation $`_W`$ makes no contribution to $`\overline{D}_{yy}^{}`$, and one recovers the physically incorrect result $`\overline{D}_{yy}^{}=\mathrm{\Gamma }=D`$.

## VII Concluding Remarks

The present contribution has developed a generic scheme for computing the mean velocity vector, $`\overline{𝐔}^{}`$, and dispersivity dyadic, $`\overline{𝐃}^{}`$, from lattice models of transport in spatially periodic arrays of obstacles, the latter representing simple models of gel electrophoresis. To the best of our knowledge, this constitutes the only analytical method available at the present time for correctly computing the dispersivity dyadic from lattice Monte Carlo models for finite fields. The power and simplicity of the present calculation scheme renders it appropriate for revisiting the Ogston sieving problems considered previously Slater and Guo (1996a, b); Slater and Treurniet (1997); Mercier and Slater (1998); Labrie et al. (2000); Boileau and Slater (2001); Mercier and Slater (2001); Mercier et al. (2001); Gauthier and Slater (2002). in order to investigate the non-trivial solute dispersion arising during such transport processes.

Furthermore, it may be possible to use the general framework developed here to analyze DNA separations in microfluidic devices Chou et al. (1999); Cabodi et al. (2002); Doyle et al. (2002), provided that rational rules can be established for the transition probabilities prevailing therein. In the case of so-called rectified Brownian motion Chou et al. (1999); Cabodi et al. (2002); Duke and Austin (1998) (or vector chromatography Dorfman and Brenner (2001)), a directional separation is achieved by applying the force at an angle with respect to the symmetry axes of the array. One should be able to analyze these devices in a manner similar to Ref. Keller et al., 2002, although it is recommended that the present approach be employed to compute the dispersivity. It would be interesting to consider a mesh which is sufficiently fine in order to account for local variations in field strength and direction arising from the fact that the electric field does not penetrate the obstacles. The latter variations were recently demonstrated to have a strong effect on the experimental realization of directional separation schemes Huang et al. (2002).

In the case of magnetosensitive arrays Doyle et al. (2002), the separation is achieved by a size-specific “hold-up” of the different DNA strands on the posts of the array. Consequently, application of the present scheme would require augmenting the hard-core repulsion of the obstacles to also account for the “attractive” residence time of the molecules when they become hooked. This problem is analogous to the inclusion of attractive forces exerted by the fibers of the gel Labrie et al. (2000), although exactly defining what is meant by the “attraction” to the posts is somewhat more ambiguous, since the attraction is not caused by a physical potential well.

###### Acknowledgements.

The author would like to thank Gary Slater of University of Ottowa for an enlightening discussion on Monte Carlo simulations. The encouragement of Jean-Louis Viovy of Institut Curie is greatly appreciated. This work was supported by a Postdoctoral Fellowship from Institut Curie.