Out-of-equilibrium Monte Carlo simulations of a classical gas with Bose-Einstein statistics

If transition probabilities between different states of a system are known, then the kinetic Monte Carlo algorithm can be used to numerically reproduce a correct description of the transient or non-equilibrium behavior Bortz et al. (1975); Fichthorn and Weinberg (1991); Voter (2005); Jansen (2012). Nevertheless, it is not uncommon that the information available is restricted to the state’s energy, and transition probabilities have to be estimated using, for example, Glauber or Metropolis algorithms. Such algorithms guarantee a correct description of the equilibrium state, but not of the out-of-equilibrium transient. Convergence towards equilibrium is ensured by imposing the detailed balance condition on the transition probabilities, see for example Binder (1997); Binder and Heermann (2010).

Let us consider a system of particles, at temperature $T$, divided into cells; transition probabilities describe jumps of one particle between neighboring cells. In Refs. Suárez et al. (2015); Martínez and Hoyuelos (2018, 2019a), it has been shown that the detailed balance condition can be used to derive a class of transition probabilities characterized by an interpolation parameter $\theta$. If $\theta=-1$, the transition probability depends on the one-particle potential in the origin cell; if $\theta=1$, it depends on the potential in the target cell; and if $\theta=0$, it depends on the potential energy change (a frequent choice in Monte Carlo simulations). The interpolation algorithm consists in the implementation of the kinetic Monte Carlo method using transition probabilities that are obtained from the potential and the interpolation parameter Suárez et al. (2015); Martínez and Hoyuelos (2018, 2019a). The purpose of this paper is to verify that this algorithm correctly describes an out-of-equilibrium system. To do this we select a specific example. The chosen system is a gas of classical particles with an effective attractive interaction which reproduces Bose-Einstein statistics. The main reason for this choice is that this system exhibits clear differences among the results generated by the mentioned numerical algorithms, in contrast, as shown below, with purely repulsive interaction of hard core.

So, we consider the following stationary non-equilibrium state. A constant force $F$ is applied along direction $x$; the external potential is $U_{i}=-Fx_{i}$, where $x_{i}$ is the position of cell $i$. The system has periodic boundary conditions along the $x$ axis and, after some time, a stationary current is established. Particles have a mean velocity proportional to $F$ and the proportionality constant is the mobility $B$. We consider $B$ (for small values of $F$) as the parameter which characterizes the non-equilibrium state, and analyze its dependence on the density for the different algorithms.

This article is organized as follows. In Sec. II, the three algorithms used for transition probabilities: Glauber, Metropolis and interpolation, are described. In Sec. III, the hard-core interaction is briefly analyzed. In Sec. IV, the basic formulae for the Bose-Einstein gas are derived. In Sec. V, we present numerical and analytical results of the mobility using Glauber, Metropolis and interpolation algorithms for the Bose-Einstein gas. The Monte Carlo simulations are performed in a one-dimensional lattice gas with periodic boundary conditions. For completeness and as a verification, in Sec. VI, we compare the previous results with the ones obtained through Molecular Dynamics. This corresponds to a more detailed kinetic description, where the position and velocity of each particle inside the cells are considered sup Finally, conclusions are presented in Sec. VII.

According to the Glauber algorithm Glauber (1963), in a Monte Carlo step the transition from state $i$ to state $j$, with energies $E_{i}$ and $E_{j}$ respectively, has a probability

where $\beta=(k_{B}T)^{-1}$.

In the Metropolis or Metropolis-Hastings algorithm Hastings (1970), the probability is

Note that Metropolis algorithm is faster: if $E_{i}=E_{j}$, we have $p_{i,j}^{G}=1/2$ and $p_{i,j}^{M}=1$. We have to take this difference into account to have the same time scale in both cases. If $P$ is the number of jump attempts per unit time, then the transition probabilities are

where a factor 2 is included in $W_{i,j}^{G}$ to compensate the speed difference between algorithms.

The energy of a given configuration $\{n_{i}\}$ is

where $\phi_{n_{i}}$ represents the interaction energy (or configuration energy) among the $n_{i}$ particles in the cell, in local equilibrium, and $U_{i}$ is an external potential. Only processes where one particle jumps to a neighboring cell are allowed.

A different perspective is adopted for the interpolation algorithm. Instead of using the state energy, the mean field potential $V_{i}$ for one particle, produced by all the other particles in the cell, is considered. We call $\theta_{i}$ the interpolation parameter in cell $i$. The transition probability from cell $i$ to cell $i+1$ is given by

where $\Delta U=U_{i+1}-U_{i}$ and $\Delta V=V_{i+1}-V_{i}$. See Ref. Vattulainen et al. (1998) for a related approach in which the transition probability depends on the sum of energies in origin and target sites. Note that $\theta_{i}$ and $V_{i}$ can be functions of the number of particles, $n_{i}$, in a cell $i$.

A relationship between $V_{i}$ and $\phi_{n_{i}}$ has to be established. This can be done using the equilibrium distribution of particles $\bar{n}_{i}$. From the transition probabilities (6), a Fokker-Planck equation can be derived, whose equilibrium solution is

where $\mu$ is the chemical potential, see Ref. Suárez et al. (2015). On the other hand, $\bar{n}_{i}$ can also be obtained from the grand partition function for cell $i$ with energy $\phi_{n_{i}}+n_{i}U_{i}$. Combining both results, it can be shown Suárez et al. (2015) that

where $V_{i}$ is evaluated at the average number $\bar{n}_{i}$ ($V_{i}$ is a continuous function of $\bar{n}_{i}$ but, when implementing the algorithm for specific realizations, $V_{i}$ has to be evaluated only at integer values of $n_{i}$).

The interpolation parameter $\theta_{i}$ can be simplified in the definition of transition probabilities (6) using the following relationship that holds in the absence of a phase transition:

where $V^{\prime}_{i}$ is the the derivative of $V_{i}$ with respect to the number of particles. The demonstration of (9) can be found in Ref. Martínez and Hoyuelos (2019b); it is obtained from the fact that the energy difference between final and initial states is equal to the difference between the energy barriers (the exponents in $W_{i,i+1}$) of the forward and backward processes. The transition probability form cell $i$ to $i+1$ becomes

In summary, transition probabilities in the interpolation algorithm are not obtained directly from the state energy, but are functions of the potential $V_{i}$ that is obtained from (8). The method is more involved than Glauber or Metropolis, and it only applies to system of particles with local interactions, where the energy can be written as in (5). The advantage is that it holds out of equilibrium.

Here, models where $\phi_{n_{i}}$ is a simple function of $n_{i}$ are considered. Long range interactions, such as the Coulomb potential, can not be represented by these methods, where the energy of one cell, given by $\phi_{n_{i}}+n_{i}U_{i}$, does not depend on the number of particles in nearby cells. Interaction between particles in different cells is neglected. For a specific (short range) interaction potential, the configuration energy $\phi_{n_{i}}$ in general has to be obtained numerically. If the interaction potential is $\mathcal{U}_{i}(\mathbf{q}_{1},\dots,\mathbf{q}_{n_{i}})$, with $\mathbf{q}_{j}$ the position of particles in cell $i$, the configuration energy $\phi_{n_{i}}$ is given by

where average $\langle\ \rangle^{0}$ is computed with the probability distribution of non-interacting particles (that is, particles randomly distributed).

The main purpose of this paper is to check the different algorithms with an effective Bose-Einstein interaction. But before that, we wish to briefly analyze other simple potential: hard core. These potentials represent two extreme situations: purely attractive or purely repulsive interactions. It is shown in this section that the mobility of hard spheres against concentration is the same independently of the simulation algorithm. Any of the three options, Glauber, Metropolis or interpolation, can be used for the description of the out-of-equilibrium behavior of hard spheres in a discrete lattice.

The hard-core interaction is given by

Hard core imposes that, in any cell, $n_{i}$ takes values 0 or 1 only. Let us consider a jump from cell $i$ to cell $i+1$. The energy of the initial configuration $\{\cdots,n_{i},n_{i+1},\cdots\}$ is $E_{\text{ini}}$, and the energy of the final configuration $\{\cdots,n_{i}-1,n_{i+1}+1,\cdots\}$ is $E_{\text{fin}}$. The external energy change is $\Delta U=-Fa$, with $a$ the cell size. Then, the energy change is

where it was assumed that $n_{i}=1$ so that the jump process from $i$ to $i+1$ is possible (there must be a particle to make the jump).

Let us call $v$ the mean velocity of particles when a force $F$ is applied. Analytical expressions for the mobility, defined as $B=v/F$, can be derived for the three algorithms to compare with numerical results. Mobility is obtained from the mean current:

since $J=v\,\bar{n}/a=BF\bar{n}/a$, so that $B=J\,a/(F\bar{n})$. We call the mobilities for Glauber and Metropolis algorithms $B^{G}$ and $B^{M}$ respectively. In the ideal system, where effects of interactions are neglected, the mobility is given by the Einstein relation: $B_{0}=\beta D_{0}$, where $D_{0}=Pa^{2}$ is the diffusivity.

First, for the Metropolis algorithm, the transition probability for a jump from $i$ to $i+1$, assuming that $F>0$, is

expression that can be written as

And, for the jump from $i+1$ to $i$,

where it is assumed that the force and the cell size satisfy $Fa\ll k_{B}T$ so that $e^{-\beta Fa}\simeq 1-\beta Fa$; this approximation is necessary to get a linear relation between mean velocity $v$ and force $F$, from which the mobility is obtained; the approximation corresponds to the linear regime of non equilibrium thermodynamics. Using (18) and (19) to evaluate the current (16), we get

and the mobility is

It is convenient to present the mobility relative to $B_{0}$, so that the result is independent of the jump rate $P$ or the cell size $a$. In the last expressions, subindex $i$ is not used for the mean number of particles, $\bar{n}$, since a homogeneous system is considered.

For the Glauber algorithm, using the definition (1), it can be shown that

Using (16), the result for the current $J^{G}$ turns out to be equal to $J^{M}$, and, of course, also the mobility:

For the interpolation algorithm, the first step is to calculate $V_{i}$. We have that

then, from Eq. (8), we have

and

When replacing this expressions in Eq. (10) for the transition probability, as mentioned before, $V_{i}$ and $V_{i}^{\prime}$ have to be evaluated at integer values of $n_{i}$. Then, using also that

it can be shown that the transition probabilities for the interpolation algorithm are equal to those of the Glauber algorithm:

Therefore, the mobility is also the same,

Figure 1 shows numerical results of the mobility using Metropolis and Glauber (or interpolation) algorithms for the hard-core interaction. The results match the expression obtained for $B/B_{0}$.

In summary, hard spheres are not useful to discriminate among algorithms. Any of them seems to work. This is not the case for the attractive interaction of a classical Bose-Einstein gas, and this is the main reason why we focus our attention on it.

In the lattice gas description, the system is divided into cells of size $a$ much smaller than the characteristic length of particle density variations. Temperature $T$ and chemical potential $\mu$ are homogeneous. Interaction energy between particles in different cells is neglected. Since spatial variations are smooth, cells are point-like in the continuous limit, and local thermal equilibrium holds. This means that each cell can be considered as an equilibrium system, although the whole system is out of equilibrium. We can write the classical grand partition function of cell $i$ as

where the energy in cell $i$ is $\phi_{n_{i}}+n_{i}U_{i}$.

We consider a configuration energy

to reproduce Bose-Einstein statistics, since in this case we have

and the average number of particles is

If we know that the system has a total number of particles $N$, the chemical potential is obtained from the relationship $N=\sum_{i}\bar{n}_{i}$, where the sum is performed in all cells. From Ec. (38), the chemical potential is

Let us consider a cell in which the potential $U_{i}$ takes a value $\mu^{\circ}$, i.e. a reference chemical potential. Then, in the previous expression, we can recognize the ideal and residual parts of the chemical potential:

This simplified description in terms of jumps between neighboring cells is intended to correctly reproduce the behavior of a system of particles which move with given velocities and interact with a space dependent potential. The effective potential, also known as statistical potential, between two bosons at distance $r$ is

where $\lambda$ determines the range of the interaction (it is equal to the de Broglie wavelength in the quantum case) (Pathria and Beale, 2011, p. 138). Notice that here we use the term “boson” to informally refer to classical particles with Bose-Einstein statistics. The statistical potential holds for small concentration, so we can expect descriptions to agree only in this limit. The concentration of a cell (in local equilibrium) is $\rho_{i}=\bar{n}_{i}/a^{3}$. Also, in the limit of small concentration, we can use the cluster expansion to obtain the chemical potential:

where $b_{2}$ is the coefficient of two-particle clusters given by

see for example Eq. (5.32) in Ref. Kardar (2007). The cluster expansion is based on the quantity $f(\mathbf{r})=e^{-\beta v(r)}-1$ as a convenient expansion parameter, with $v(r)$ the interaction potential. For short-range hard-core interactions, it is equal to $-1$ for $r\rightarrow 0$ and decays to zero for increasing $r$. In our case, it takes the value $1$ for $r\rightarrow 0$ and vanishes exponentially as $r$ increases.

From Eqs. (42) and (43) we have that the chemical potential is

where, as in Eq. (40), we can identify the ideal and residual parts.

By imposing the condition in which the chemical potential in the lattice gas description, Eq. (40), is equal to the chemical potential in the kinetic description, Eq. (44), more specifically, by matching the residual parts of the chemical potential, we have that $\bar{n}_{i}=\rho_{i}\lambda^{3}/2^{3/2}$, and therefore,

But setting this condition poses a problem in the calculation of the grand partition function, since in Eq. (35) interactions with particles in neighboring cells are neglected. This approximation is valid as long as $a$ is much larger than the interaction range, i.e. $a\gg\lambda$, a condition that is not fulfilled in (45). Still, the qualitative behavior of the chemical potential is the same in both descriptions. We expect equivalent qualitative behaviors also for the mobility although, due to the inconsistency of approximations, we can not expect quantitative agreement.

The energy change for a jump from cell $i$ to $i+1$ with effective Bose-Einstein interaction is

where the condition $n_{i}\geq 1$ has to be fulfilled in order to have the possibility of a jump from cell $i$ (for cells with $n_{i}=0$ there is no jump process and no need to calculate $\Delta E$).

With fixed boundary conditions and an external force, particles accumulate on one side of the system and, in equilibrium, they have the Bose-Einstein distribution given by Eq. (38). We have verified that Metropolis or Glauber algorithms converge to the Bose Einstein distribution in this case, as shown in Fig. 2.

Now we turn to the non-equilibrium situation. Instead of fixed boundary conditions, let us consider periodic boundary conditions. After some time, the system is in a stationary non-equilibrium state. The mean velocity of a particle in terms of the transition probabilities between cells of size $a$ is

For small concentration (ideal gas), interactions can be neglected and $\Delta E\simeq-Fa$. It can be shown that for both algorithms, Glauber and Metropolis, the mean velocity for small concentration is $v=Pa^{2}\beta F$. For Metropolis at small concentration and a positive force $F$, transition probabilities are

Replacing in (47), we obtain

For the Glauber algorithm we have

And, again, replacing in (47), we obtain

The mobility, $B$, is defined as $v/F$. The previous results for the velocity, for both algorithms, are consistent with the Einstein relation; in both cases, $B_{0}=\beta Pa^{2}=\beta D_{0}$ is obtained. Subindex 0 identifies the small concentration regime.

As mentioned before, mobility is obtained from the mean current, Eq. (16). It can be shown that, for the Glauber algorithm,

and that the corresponding expression for the Metropolis algorithm can be simplified to

see the Appendix for details.

In Fig. 3, we show the numerical results of the mobility obtained with both algorithms, compared with the analytical expressions (V) and (49). The agreement between numerical and analytical results supports their validity. Both algorithms coincide in the limit of small concentration, where $B=B_{0}$ and interactions can be neglected. But as the mean number of particles is increased, discrepancies grow. There is a quantitative difference of about 30% between both predictions of $B$ for larger $\bar{n}$. It is well known that we can not expect a correct description of a non-equilibrium state with these algorithms. Nevertheless, it is interesting to determine how far from the correct result they are. This is the purpose of the next paragraphs.

For the interpolation algorithm, using Eq. (36) for $\phi_{n_{i}}$ in (8), we have

Notice that $V_{i}$ is equal to the residual part of the chemical potential, see Eq. (40), and Eq. (8) corresponds to the Widom insertion method, where $\phi_{n_{i}+1}-\phi_{n_{i}}$ is the insertion energy, i.e. the interaction energy needed to insert one particle. Using this result for $V_{i}$, and its derivative, in (10), the transition probability is

With this information, we can calculate the mobility $B^{I}$ for the interpolation algorithm (see the Appendix). The result is

As in the previous results, a homogeneous system, in which $\bar{n}_{i}=\bar{n}_{i+1}=\bar{n}$, is considered for the calculation of mobility. Figure 4 shows a good agreement between this theoretical result and kinetic Monte Carlo simulations.

The result for the mobility obtained with the interpolation algorithm is qualitatively different from the results of Glauber and Metropolis algorithms. While for the interpolation algorithm we have a mobility which increases with $\bar{n}$, for the other algorithms it decreases (see Fig. 3). As mentioned before, we know that Glauber and Metropolis algorithms are designed to give the correct equilibrium state; they should not be applied to non-equilibrium situations, but it is interesting to evaluate the error. According to the interpolation algorithm, which is designed for non-equilibrium states (with limitations that are summarized in the conclusions), the error increases with concentration. Glauber and Metropolis algorithms give the correct result for the mobility only in the limit of small $\bar{n}$, where interactions can be neglected.

The objective of this section is to obtain the mobility of a boson gas in the context of a classical kinetic description which includes velocity of particles, and compare with the results of the previous sections.

The method is to numerically obtain the self-diffusivity $D$ and use the Einstein relation to calculate the mobility $B=\beta D$. The Green-Kubo formula (see for example (Kubo et al., 1998, Sec. 4.6.2) or (Reichl, 1998, Sec. S10.G)) is used to obtain the self-diffusivity from the velocity autocorrelation function:

We wish to emphasize that with this method the mobility is obtained from simulations of the system in equilibrium, without a force applied. The main purpose of this paper is to check algorithms for Monte Carlo simulations out of equilibrium. The motivation of this section is to obtain mobility using molecular dynamics, and in this context we can choose the method based on simplicity.

To compare with results of the previous sections, we have to calculate the mobility $B$ relative to its value for the ideal gas, $B_{0}$. From the kinetic theory of transport in dilute gases (see for example (McQuarrie, 2000, Sec. 16-1)), we know that the self-diffusivity in the ideal case, and therefore $B_{0}$, behaves as $1/\rho$, where $\rho$ is the particle density. The proportionality constant, between $B_{0}$ and $1/\rho$, is numerically set in our results so as to have $B/B_{0}\rightarrow 1$ for $\rho\rightarrow 0$.

Then, we perform molecular dynamic simulations of a system of particles with a given density $\rho$, in equilibrium, which interact among them with the statistical potential of Eq. (41). We do this for different densities and obtain the mobility through the velocity autocorrelation function. We focus our attention on the slope of the mobility for small concentrations, since this is the limit for which the statistical potential (41) holds.

In Fig. 5, we show the data obtained for $B$ as a function of $\bar{n}=\rho\lambda^{3}/2^{3/2}$. This result supports the validity of the interpolation algorithm since the same qualitative behavior is obtained: an approximate linear increase of mobility with concentration with slopes of the same order, see Fig. 5.

We study mobility in a system of classical particles with interactions. Two different approaches are possible: a lattice gas with transition probabilities among cells, or a kinetic description in a continuous space. A fundamental problem of the lattice gas description is to determine transition probabilities when only state energies are known. Monte Carlo simulations with Glauber or Metropolis algorithms correctly describe the equilibrium state but they are not supposed to hold out of equilibrium. Instead, the interpolation algorithm was recently introduced Suárez et al. (2015); Martínez and Hoyuelos (2018, 2019a, 2019b) in order to obtain transition probabilities that hold out of equilibrium; it is summarized in Eqs. (8) and (10). Limitations of the method are: the system should be in local thermal equilibrium (i.e. deviations from equilibrium have a characteristic length much larger than the cell size), no phase transition occur and the interaction is local [more precisely, the interaction energy can be written as in (5)]. Also, the information provided by the algorithm is incomplete if the jump rate $P$ has a non-trivial dependence on concentration, as it happens, for example, for diffusion in a solid Martínez and Hoyuelos (2019b).

We calculate the mobility in a non-equilibrium stationary state: a force is applied to a one-dimensional array of cells with periodic boundary conditions and homogeneous density. For hard core interaction, the same results are obtained for Glauber, Metropolis or interpolation algorithms. Instead, for an attractive potential such that the Bose-Einstein statistics is reproduced, important differences are obtained among algorithms. There is a difference of about 30% between the Glauber and Metropolis predictions for large concentration; in both cases, mobility decreases with concentration towards an asymptotic value. The mobility obtained from the interpolation algorithm qualitatively differs from the ones of the other methods. Instead of decreasing with concentration, it increases. This means an unbounded increasing error for the Glauber and Metropolis algorithms. They can be used to calculate mobility only in the limit of small concentration, where $B\rightarrow B_{0}$.

The lattice gas description should be consistent with the kinetic description. So, we also consider particles moving in a continuous space and interacting with a statistical potential (41) which corresponds to the Bose-Einstein distribution. The mobility obtained from molecular dynamics simulations, Fig. 5, is in qualitative agreement with the prediction of the interpolation algorithm.

In this appendix, the expressions (V) and (49) for the mobility for Glauber and Metropolis algorithms are derived.

Using (46) for the energy change, the transition probabilities for the Glauber algorithm are

Replacing these expressions in the current $J=\langle\mathcal{J}\rangle$, with

and making approximations for $\beta Fa\ll 1$, we have

The first term of the right side cancels due to the symmetry of the homogeneous stationary state. Knowing that $B=J\,a/(F\bar{n})$ and that $B_{0}=\beta Pa^{2}$, we have

The probability of having $n_{i}$ particles knowing that the average value is $\bar{n}$, for the Bose-Einstein distribution, is given by

see for example (Pathria and Beale, 2011, p. 152). Using this probability, we obtain Eq. (V) for the mobility in the Glauber algorithm. Actually, a correction of order $\beta Fa$ should be added to $p_{\text{BE}}(n_{i})$ to use the probability which corresponds to the stationary non-equilibrium state, since Eq. (59) holds in equilibrium. But this correction cancels when only terms up to order $\beta Fa$ are kept in the equation for the current.

A similar process is applied to obtain the mobility for the Metropolis algorithm. Replacing the expression (46) for $\Delta E$ in the Metropolis transition probabilities, we have

Assuming that $\beta Fa\ll 1$ and considering all possible combinations of $n_{i}$ and $n_{i+1}$, we get

The average of this expression can be written as

These sums can be simplified. Replacing $p_{\text{BE}}(n_{i})=q^{n_{i}}/(\bar{n}+1)$ with $q=\bar{n}/(\bar{n}+1)$, and using a symbolic manipulator such as Maxima Maxima, a Computer Algebra System , we obtain

Using the relationship between current and mobility, $B=J\,a/(F\bar{n})$, Eq. (49) is immediately obtained.

Finally, the calculation of the mobility for the interpolation algorithm is simpler. We have the transition probabilities in (51), and the mean current is

where it was considered that $\langle n_{i}\rangle=\langle n_{i+1}\rangle=\bar{n}$ and that fluctuations in different cells are independent, so $\langle n_{i}n_{i+1}\rangle=\bar{n}^{2}$. From this equation for $J$, the result for the mobility (52) is obtained.