text:
# Static and Dynamic Monte Carlo Simulations of Br Electrodeposition on Ag(100)

## 1 Introduction

A number of experimental techniques have emerged in recent years which enable studies of the structure and dynamics of electrode-electrolyte interfaces in situ, i.e., without removing the electrode from the liquid . This new atomic scale information allows greater understanding of the deposition process where, in contrast to high-vacuum surface-science methods, the electrochemical environment provides direct experimental control over the electrochemical potentials of the adsorbates, rather than over their coverages or fluxes. These new in situ experimental techniques, together with new methods to manipulate adsorbates at the interface on the nanometer scale , promise to make electrodeposition an even more important method for low-cost, high-quality materials synthesis than it is today . These developments emphasize the importance of developing theoretical and computational methods to study the dynamical aspects of electrodeposition in microscopic detail.

Motivated by these developments in electrochemical surface science, we here present a static and dynamic Monte Carlo (MC) simulation study of a simple system which can serve as a prototype for more complicated ones of applied interest. This system is a submonolayer of Br electrodeposited onto a single-crystal Ag(100) electrode. Recently it was investigated at room temperature by in-situ surface X-ray scattering (SXS) and chronocoulometry by Ocko, Wang, and Wandlowski . That study detected a second-order phase transition between a disordered low-coverage phase at more negative electrode potentials and a doubly degenerate c$`(2\times 2)`$ ordered phase with coverage $`1/2`$ at more positive potentials. (This phase has also been observed by low-energy electron scattering (LEED) for Br adsorbed on Ag(100) in vacuum .) The transition was shown to belong to the two-dimensional Ising universality class . The scattering data  also indicate that Br adsorbs at the four-fold hollow sites between the surface Ag atoms, consistent with recent in situ X-ray absorption fine structure (XAFS) results . The fluctuations around the adsorption sites are small, suggesting that a lattice-gas approximation for the adlayer structure and dynamics is appropriate. Here we investigate in detail such a lattice-gas model with a nearest-neighbor excluded-volume interaction and lateral dipole-dipole repulsive interactions . The dipole interaction strength and the electrosorption valency, which are parameters of the model, are estimated by fitting simulated equilibrium adsorption isotherms to experimental chronocoulometric data .

Cyclic voltammograms (CVs) at room temperature show two main features : a broad shoulder at more negative potentials which ends in a sharp peak at more positive potentials. By SXS it was confirmed that the sharp peak corresponds to a second-order phase transition in the layer of adsorbed Br . However, the cause of the broad shoulder has been less clear. It was previously suggested  that interactions with surface water are responsible for the shoulder, but equilibrium MC simulations by Koper , as well as our own work, show the same features in simulations without water. We show here that the broad shoulder in the simulated CV is due to fluctuations within the disordered adlayer. Locally, these fluctuations resemble other ordered phases with coverage 1/4, which exist in the model at much lower temperatures.

Going beyond these equilibrium simulations, we further present a dynamic MC study of the ordering and disordering processes between the disordered phase and the c$`(2\times 2)`$ phase, assuming a perfectly stirred solution so that bulk diffusion effects can be ignored. For ordering, we expect similar behavior to the random sequential adsorption process of the hard-square model with lateral diffusion . At late stages of the ordering, this dynamic model is expected to be in the same dynamic universality class as the ferromagnetic  or antiferromagnetic  Ising model in zero magnetic field following sudden cooling from infinite temperature to a temperature below the critical point. These phase ordering dynamics with non-conserved order parameter, often called Model A in the literature on dynamics of phase ordering and phase separation , are driven by the reduction of interface curvature . We further present dynamic results for the disordering process and discuss how simulation results for this process can be compared with experiment.

We also apply the dynamic MC algorithm to study sweep-rate dependence in simulated CV experiments. Results are presented for the response of the system to a wide range of sweep rates, and we suggest some comparisons between simulation results and experiment. Some preliminary results of this study were included in Refs. .

The remainder of this paper is organized as follows. In Sec. 2 we introduce the lattice-gas model and the relation between the electrochemical potential and the experimentally measured electrolyte concentration and electrode potential. In Sec. 3 we present our equilibrium results, and in Sec. 4 we introduce our dynamic MC algorithm. In Sec. 5 we present dynamic simulations of sudden potential-step experiments, and in Sec. 6 we present dynamic simulations of CV experiments. In Sec. 7 we give a summary and our conclusions. The Appendix describes the nonlinear fitting procedure used to determine the lattice-gas parameters. Animations and additional figures are available on the World Wide Web .

## 2 Lattice-Gas Model

The lattice-gas model consists of a two-dimensional $`L\times L`$ square array of Br adsorption sites, corresponding to the four-fold hollow sites on the Ag(100) surface. Periodic boundary conditions were used to eliminate edge effects in the relatively small systems it is possible to simulate. The grand-canonical lattice-gas Hamiltonian is

$$=\underset{i<j}{}\varphi _{ij}c_ic_j\overline{\mu }\underset{i=1}{\overset{L^2}{}}c_i,$$
(1)
where $`i`$ and $`j`$ denote sites on the lattice, $`c_i`$ denotes the occupation state of site $`i`$, which is either 0 (empty) or 1 (occupied), $`_{i<j}`$ denotes a sum over all pairs of sites on the lattice, $`\varphi _{ij}`$ denotes the lateral interaction energy of the pair $`(i,j)`$, and $`\overline{\mu }`$ is the electrochemical potential. The sign convention used here is such that $`\varphi _{ij}<0`$ denotes a repulsive interaction, and $`\overline{\mu }>0`$ favors adsorption. We measure $`\varphi _{ij}`$ and $`\overline{\mu }`$ in units of meV/pair and meV/particle, respectively. \[For brevity we write both simply as meV.\]

The electrochemical potential $`\overline{\mu }`$ is the intensive variable which controls the average number of Br in the adlayer. In the weak-solution approximation it is related to the electrode potential by

$$\overline{\mu }=\overline{\mu }_0+k_\mathrm{B}T\mathrm{ln}\frac{[C]}{[C_0]}e\gamma E,$$
(2)
where $`\overline{\mu }_0`$ is an arbitrary reference level, $`k_\mathrm{B}`$ is Boltzmann’s constant, $`T`$ is the absolute temperature, $`e`$ is the elementary charge unit, $`[C]`$ is the concentration of Br- in solution, $`[C_0]`$ is an arbitrary reference concentration (here we use $`C_0=1`$ mM), $`\gamma `$ is the electrosorption valency ($`\gamma e`$ is the charge transfered through the external circuit for each Br ion adsorbed on the electrode ), and $`E`$ is the electrode potential in mV.

A variety of lateral interactions have been previously explored for this model, including finite and infinite nearest-neighbor repulsion and screened dipole-dipole interactions . From these previous studies and our own exploratory simulations, we find that the interactions are quite adequately described with a nearest-neighbor excluded-volume interaction (infinite nearest-neighbor repulsion) related to the large ionic radius of Br- (1.96 Å ), plus a long-range dipole-dipole repulsion. The pairwise lateral interactions used here are

$$\varphi (r)=\{\begin{array}{cc}\mathrm{}\hfill & r=1\hfill \\ \frac{2^{3/2}\varphi _{\mathrm{nnn}}}{r^3}\hfill & \sqrt{2}r5\hfill \\ 0\hfill & r>5\hfill \end{array},$$
(3)
where $`r`$ is the separation of an interacting Br pair within the adlayer, measured in units of the Ag(100) lattice spacing ($`a=2.889`$ Å ), and $`\varphi _{\mathrm{nnn}}`$ (which is negative) is the lateral dipole-dipole repulsion between next-nearest neighbors. In other work it has been noticed that reducing the infinite nearest-neighbor repulsion to finite but large values has very little effect on the isotherms . All interactions are truncated beyond five lattice spacings to simplify computation. This truncation can be seen as an approximate cut-off due to screening. With this cut-off, the total interaction energy of a fully ordered c$`(2\times 2)`$ phase is approximately $`87\%`$ of the interaction energy without a cut-off.

The Ag(100) substrate provides a square adsorption lattice of hollow sites for the commensurate Br adlayer. The nearest-neighbor excluded-volume interactions give rise to two different sublattices which have a geometry corresponding to the two colors on a checkerboard. Each of the two degenerate c$`(2\times 2)`$ phases can occupy only one of the two sublattices, which we label $`A`$ and $`B`$. The $`A`$ sublattice coverage is the fraction of occupied sites on sublattice $`A`$, defined as $`\mathrm{\Theta }_A=N_A^1_{iA}^{N_A}c_i`$, where $`N_A`$ is the number of sites on sublattice $`A`$ and $`_{iA}^{N_A}`$ denotes a sum over all sites on that sublattice. The $`B`$ sublattice coverage, $`\mathrm{\Theta }_B`$, is analogously defined. For this system, there are two linearly independent observables of interest: the total Br coverage, $`\mathrm{\Theta }=(\mathrm{\Theta }_A+\mathrm{\Theta }_B)/2`$, and the staggered coverage, $`\mathrm{\Theta }_\mathrm{S}=\mathrm{\Theta }_A\mathrm{\Theta }_B`$, which is the order parameter for the c$`(2\times 2)`$ phase.

## 3 Equilibrium Monte Carlo Simulations

Monte Carlo simulations were performed at room temperature ($`k_\mathrm{B}T=25`$ meV, corresponding to $`T290`$ K) to determine the equilibrium properties of the Br adlayer for different values of $`\overline{\mu }`$. For these equilibrium simulations, we used the Metropolis algorithm  with single-site updates at randomly selected sites.

The main purpose of these equilibrium simulations was to obtain adsorption isotherms for parameter estimation by fitting to experimental data. The most important region of comparison is in the disordered phase, where the fluctuations are small. It was therefore sufficient to use relatively small systems with $`L=`$32. No significant differences were noticed between data for $`L=32`$ and $`L=128`$ (see below). No systematic finite-size scaling was performed near the phase transition. Thus, the results obtained in the critical region are only approximate.

For a given temperature and chemical potential, we generated a set of samples from which we computed the averages $`\mathrm{\Theta }(T,\overline{\mu })`$ and $`|\mathrm{\Theta }_\mathrm{S}(T,\overline{\mu })|`$, where $``$ denotes an average over a set of statistically independent equilibrium samples. For convenience, we drop the $``$ notation for the equilibrium isotherms.

The critical value of the electrochemical potential and the shape of the adsorption isotherm in the disordered phase are quite sensitive to the strength of the long-range repulsive interactions. To determine the dipole interaction strength, $`\varphi _{\mathrm{nnn}}`$, and the electrosorption valency, $`\gamma `$, we therefore simultaneously compared the simulated room-temperature adsorption isotherm, $`\mathrm{\Theta }(\overline{\mu })`$, with experimental isotherms, $`\mathrm{\Theta }^{\mathrm{exp}}(E)`$, obtained from chronocoulometry . The nonlinear fitting method is described in the Appendix. The resulting parameter values are $`\varphi _{\mathrm{nnn}}=26\pm 2`$ meV, $`\gamma =0.73\pm 0.03`$, and $`\overline{\mu }_0=0.67\pm 0.03`$ meV. These results are consistent with those found previously by others

The fitted MC isotherm is shown together with the chronocoulometric coverage isotherms in Fig. 1. For the two lowest concentrations, the agreement is very good over the whole range of electrode potentials, while the agreement in the ordered phase is somewhat less satisfactory for the highest concentration. This disagreement at higher concentrations suggests a break-down in the weak-solution approximation and the ultimate need for a more detailed treatment of the solution.

Figure 2 shows the simulated $`\mathrm{\Theta }(\overline{\mu })`$ and $`|\mathrm{\Theta }_\mathrm{S}(\overline{\mu })|`$ at $`T290`$ K, with $`\varphi _{\mathrm{nnn}}=26`$ meV as obtained from the fits to the experimental isotherms. The order parameter, $`\mathrm{\Theta }_\mathrm{S}`$, is zero in the disordered phase, while it is near $`\pm 1`$ deep in either of the two degenerate ordered phases. A second-order phase transition between the low-coverage disordered phase and the c$`(2\times 2)`$ ordered phase is observed at a critical electrochemical potential, $`\overline{\mu }_c180\pm 5`$ meV. As $`\overline{\mu }\overline{\mu }_c^+`$, $`|\mathrm{\Theta }_\mathrm{S}|(\overline{\mu }\overline{\mu }_c)^{1/8}`$ , as expected for the Ising universality class.

The coverage isotherm shows a long, gradual slope in the disordered phase, which is not present in a model with only nearest-neighbor exclusion (called the hard-square (HS) model) . As a result, $`\overline{\mu }_c`$ is considerably larger than that of the pure hard-square model, $`\overline{\mu }_c^{\mathrm{HS}}33`$ meV 222 The value $`\overline{\mu }_c^{\mathrm{HS}}33`$ meV is obtained as $`\overline{\mu }_c^{\mathrm{HS}}=k_\mathrm{B}T\mathrm{ln}z_c`$, using precise numerical values for the critical activity of the hard-square model, $`z_c=3.7962(1)`$  and $`z_c=3.7959(1)`$ , obtained by series expansions and transfer-matrix finite-size scaling, respectively.. The simulated and experimental isotherms show coverages at $`\overline{\mu }_c`$ of about $`0.37\pm 0.01`$ and $`0.36\pm 0.02`$, respectively. Within statistical errors, this is the same result as the pure hard-square model, $`\mathrm{\Theta }_c^{\mathrm{HS}}=0.368`$ 333 Precise numerical results have been obtained for the pure hard-square model by transfer-matrix calculations, including $`\mathrm{\Theta }_c^{\mathrm{HS}}=0.36776(1)`$  and $`\mathrm{\Theta }_c^{\mathrm{HS}}=0.368(1)`$ .. The adsorption isotherm, $`\mathrm{\Theta }(\overline{\mu })`$, has its maximum slope at $`\overline{\mu }_c`$. Corresponding maxima in the slopes are also seen in the experimental isotherms in Fig. 1. Simulation data are included in Fig. 2 for both $`L=32`$ and $`L=128`$. On the scale of this figure, differences between the two sizes are visible only in the data for $`|\mathrm{\Theta }_\mathrm{S}|`$ very close to $`\overline{\mu }_c`$, where $`\mathrm{d}|\mathrm{\Theta }_\mathrm{S}|/\mathrm{d}\overline{\mu }`$ is expected to diverge with increasing $`L`$ as $`L^{7/4}`$ 444 At $`\overline{\mu }_c`$ in a finite system of size $`L`$, $`\mathrm{d}|\mathrm{\Theta }_\mathrm{S}|/\mathrm{d}\overline{\mu }L^{\gamma /\nu }`$, where $`\gamma `$ and $`\nu `$ are the critical exponents for the ordering susceptibility and the correlation length, respectively. For the two-dimensional Ising model $`\gamma =7/4`$ and $`\nu =1`$ ..

A quasi-equilibrium CV (the CV in the limit of vanishing potential sweep rate) can be obtained from the simulated equilibrium isotherm using the simple relation,

$$𝒥=\gamma eM\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t}=\gamma eM\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}\overline{\mu }}\frac{\mathrm{d}\overline{\mu }}{\mathrm{d}E}\frac{\mathrm{d}E}{\mathrm{d}t}=\gamma ^2e^2M\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}\overline{\mu }}\frac{\mathrm{d}E}{\mathrm{d}t}\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}\overline{\mu }},$$
(4)
where $`𝒥`$ is the voltammetric current density (oxidation currents positive), $`M`$ is the total number of adsorption sites per unit area ($`a^2=1.198\times 10^{15}`$ sites/cm2 ), and $`\mathrm{d}E/\mathrm{d}t`$ is the sweep rate. Figure 3 shows $`\mathrm{d}\mathrm{\Theta }/\mathrm{d}\overline{\mu }`$, which exhibits the same broad shoulder and sharp peak seen in experimental CVs . The shoulder has previously been attributed to interactions with surface water . However, the fact that the lattice-gas model does not contain explicit water molecules indicates that the shoulder is most likely caused by configurational fluctuations in the Br adlayer, as was also pointed out by Koper . Differences between the simulation data for $`L=32`$ and 128 are only visible near the peak, which is expected to diverge logarithmically with $`L`$ 555 At $`\overline{\mu }_c`$ $`\mathrm{d}\mathrm{\Theta }/\mathrm{d}\overline{\mu }`$ is expected to be proportional to the heat capacity of the adlayer. The heat capacity of the two-dimensional Ising model diverges logarithmically with $`L`$ .. This quasi-equilibrium CV is compared with dynamically simulated CVs for different potential sweep rates in Sec. 6.

Having obtained an estimate for the lateral interaction constant by equilibrium simulations at room temperature, we briefly explore the full phase diagram of the model with this interaction strength, shown in Fig. 4. We expect the phase diagram for this system to be topologically similar to that for the hard-square model with next-nearest neighbor repulsion  (even though the low-coverage ordered phases in the two systems are different). In addition to the disordered phase and the two equilibrium ordered phases \[the low-temperature c$`(4\times 2)`$ phase with $`\mathrm{\Theta }=1/4`$ and the c$`(2\times 2)`$ phase with $`\mathrm{\Theta }=1/2`$ which is observed at all temperatures for sufficiently large $`\overline{\mu }`$\], there is also a metastable phase, p$`(2\times 2)`$ with $`\mathrm{\Theta }=1/4`$, which lies closely above the c$`(4\times 2)`$ phase in energy. To facilitate equilibration at low temperatures, in this simulation we used a heat-bath algorithm which updates clusters  consisting of randomly chosen nearest-neighbor pairs of sites. The presence of the p$`(2\times 2)`$ phase made low-temperature equilibration very difficult, even with the pair-update algorithm. Other phases with $`\mathrm{\Theta }<1/4`$ probably also exist as a consequence of the long-range repulsion; however, we expect these phases to be present only at extremely low temperatures, and their effects should be negligible at room temperature. The phase boundaries have been determined by ground-state calculations for $`T=0`$, by the location of the sharp peak in the quasi-equilibrium CV for $`T>100`$ K, and by inspection of the pair-update MC sample configurations for $`0<T<100`$ K.

The broad shoulder in the CV is caused by configurational fluctuations in the disordered phase, which locally resemble the $`\mathrm{\Theta }=1/4`$ phases . At the experimentally unrealistic low temperatures where the c$`(4\times 2)`$ phase is observed in the model, the shoulder breaks up into two transitions between the disordered phase and the c$`(4\times 2)`$ phase. Although the c$`(4\times 2)`$ phase is not experimentally observable in this system, local fluctuations which resemble the $`\mathrm{\Theta }=1/4`$ phases can be seen in the model in equilibrium at room temperature and should be observable in the diffuse intensity of SXS experiments. We therefore next present simulated equilibrium SXS intensities in the room-temperature disordered phase at $`\mathrm{\Theta }1/4`$.

We define the structure factor as $`S(\stackrel{}{k})=|\widehat{c}(\stackrel{}{k})|^2`$, where $`|\widehat{c}(\stackrel{}{k})|^2`$ denotes the squared absolute value of the complex Fourier transform of the real-space adlayer configuration $`\{c_i\}`$, and $`\stackrel{}{k}=(k_x,k_y)`$ is the wave vector in the plane of the adlayer (given in units of the inverse lattice constant, $`a^1`$). Equivalently, $`S(\stackrel{}{k})`$ can be expressed as the Fourier transform of the adsorbate correlation function, $`C(\stackrel{}{r_i}\stackrel{}{r_j})=c_ic_j`$. The Fourier transform is normalized such that $`S(\stackrel{}{k}=(0,0))=(L\mathrm{\Theta })^2`$. In the Born approximation, the structure factor is proportional to the SXS intensity . In Fig. 5(a) and (b) we show the simulated structure factor for the disordered room-temperature phase at $`\overline{\mu }=+100`$ meV, where $`\mathrm{\Theta }1/4`$. Small domains of c$`(2\times 2)`$, c$`(4\times 2)`$, and p$`(2\times 2)`$ symmetry are noticeable in snapshots of the adsorbate configurations \[Fig. 5(c)\], and they also affect the simulated diffuse scattering intensity. Referring to Fig. 5(b), Bragg peaks would occur at the points marked $`\times `$ for a fully ordered $`c(2\times 2)`$ phase, at the points marked $``$ for $`p(2\times 2)`$, and at the points marked $`\mathrm{}`$ for $`c(4\times 2)`$. Due to the disorder, all of these peaks \[except for the one at $`\stackrel{}{k}=(0,0)`$\] are very broad and combine to give the anisotropic diffuse scattering intensity shown in the figure. The individual peaks would become narrower and more distinct at temperatures less than room temperature (but still above the critical temperature for the low-temperature ordered phase) or with stronger dipole-dipole interactions.

## 4 Dynamic Monte Carlo Simulations

Next we investigate the dynamics of the Br adlayer using a dynamic MC algorithm. The microscopic processes considered are adsorption/desorption, in which a particle enters or leaves the adlayer, as well as nearest-neighbor and next-nearest neighbor lateral diffusion. The thermally activated, stochastic barrier-hopping picture of the microscopic dynamics used here is described in further detail elsewhere .

We assume a simple linear form for the local energy landscape, which is equivalent to the Butler-Volmer approximation with symmetry parameter 1/2 . See Fig. 6. The initial and final states, $`I`$ and $`F`$, respectively, are represented by lattice-gas configurations and differ by a single microscopic move. We label each type of microscopic move (adsorption, desorption, nearest-neighbor and next-nearest neighbor diffusion) by the index $`\lambda `$. The state $`T_\lambda `$ represents an intermediate transition state of higher energy between $`I`$ and $`F`$, which is associated with the move $`\lambda `$. This intermediate state cannot be represented by a lattice-gas configuration, and we associate with it a “bare” barrier $`\mathrm{\Delta }_\lambda `$. The energy of the transition state is then approximated as

$$U_{T_\lambda }=\frac{U_I+U_F}{2}+\mathrm{\Delta }_\lambda ,$$
(5)
where $`U_{T_\lambda }`$, $`U_I`$, and $`U_F`$ are the energies of the states $`T_\lambda `$, $`I`$, and $`F`$, respectively.

Although other choices of the transition probability are used in the literature as well , we here approximate the probability per MC step per site (MCSS) of making a transition from $`I`$ to $`F`$ by the simple one-step Arrhenius rate

$$R(F|I)=\nu \mathrm{exp}\left(\frac{U_{T_\lambda }U_I}{k_\mathrm{B}T}\right)=\nu \mathrm{exp}\left(\frac{\mathrm{\Delta }_\lambda }{k_\mathrm{B}T}\right)\mathrm{exp}\left(\frac{U_FU_I}{2k_\mathrm{B}T}\right),$$
(6)
where the rate constant $`\nu `$ (in units of (MCSS)-1) sets the overall time scale for the simulation and is assumed to be the same for all microscopic processes.

Dynamic Monte Carlo steps proceed in the following way. First, a lattice site $`j`$ is randomly selected. If $`j`$ is empty, only adsorption may be attempted. If $`j`$ is occupied, both desorption and lateral diffusion may be attempted. In general, some moves will be excluded by the presence of other adatoms (corresponding to $`U_F=+\mathrm{}`$). Next, a weighted list is constructed of new potential configurations $`F`$ which would result from successfully attempting each of the allowed moves for site $`j`$. The probability to undergo each transition is then calculated from Eq. (6). Blocked moves have transition probability zero. To ensure microscopic reversibility (detailed balance), the possibility that no change in configuration occurs ($`F=I`$) must also be included in the list of new configurations. The probability to remain in the initial state is $`R(I|I)=1_{FI}R(F|I)`$.

In general, there is a possibility that $`R(I,I)`$ could be negative. However, this essentially means that the time corresponding to one MCSS is chosen too long for a consistent stochastic description. Thus, $`\nu `$ must be chosen sufficiently small that $`R(I|I)`$ remains positive. As long as this condition is satisfied, changing $`\nu `$ only corresponds to changing the overall time scale for the simulation. In our simulations we used $`\nu =1`$ (MCSS)-1, which ensured that $`R(I|I)`$ remained positive with the barriers used here.

The main price paid for the simplicity of coding our algorithm is that it is rather slow in parameter regimes where the acceptance probabilities are very low. However, except for extremely fast processes, in which the acceptance probabilities are large and change significantly within a single MCSS, it is equivalent to more sophisticated continuum-time algorithms .

The adsorption/desorption barrier is $`\mathrm{\Delta }_\mathrm{a}`$, the nearest-neighbor diffusion barrier is $`\mathrm{\Delta }_{\mathrm{nn}}`$, and the next-nearest neighbor diffusion barrier is $`\mathrm{\Delta }_{\mathrm{nnn}}`$. In the simulations, the diffusion barriers are $`\mathrm{\Delta }_{\mathrm{nn}}=100`$ meV and $`\mathrm{\Delta }_{\mathrm{nnn}}=200`$ meV, based on ab-initio calculations  of binding energies for a single Br ion to a Ag(100) substrate in vacuum. The barrier $`\mathrm{\Delta }_{\mathrm{nn}}`$ is approximated by the difference in binding energy between the four-fold hollow site and the bridge site, while $`\mathrm{\Delta }_{\mathrm{nnn}}`$ is approximated by the difference in binding energy between the four-fold hollow site and the on-top site.

Theoretical estimates of $`\mathrm{\Delta }_\mathrm{a}`$ are extremely sensitive to the details of both the ion-surface and ion-water interactions . Based on data from calculations of the potential of mean force for halide ions in water near a Cu(100) surface , our very rough estimates indicate that values between 200 and 500 meV are not unreasonable. Since we also expect $`\mathrm{\Delta }_\mathrm{a}`$ to be significantly larger than $`\mathrm{\Delta }_{\mathrm{nnn}}`$, we examined some aspects of the behavior for $`\mathrm{\Delta }_\mathrm{a}=300`$ meV and $`\mathrm{\Delta }_\mathrm{a}=400`$ meV. Finding few qualitative differences between simulation results for these two values, we chose the lower value, $`\mathrm{\Delta }_\mathrm{a}=300`$ meV, to increase the simulation speed.

## 5 Simulated Potential Steps

One way to probe the dynamics of an adsorbed layer is to subject it to a sudden change in electrochemical potential. Dynamic simulations of such potential-step experiments are started by equilibrating the room-temperature system at a potential $`\overline{\mu }_1`$, using the simulation technique described in Sec. 3. Then, at time $`t=0`$, the potential is instantaneously stepped to a new value $`\overline{\mu }_2`$, and the simulation continues by the dynamic MC algorithm described in Sec. 4. To study the dynamics of the phase ordering and disordering processes, we examined four different potential steps: two from the disordered phase into the ordered phase, and two from the ordered phase into the disordered phase.

Before discussing each of these potential steps in detail, we give the initial adsorption, desorption, and diffusion rates at $`\overline{\mu }_2`$, immediately after the sudden potential change. We define the probability during one MCSS of accepting (when possible) an adsorption move as $`R_\mathrm{a}`$, a desorption move as $`R_\mathrm{d}`$, a nearest-neighbor diffusion move as $`R_{\mathrm{nn}}`$, and a next-nearest neighbor diffusion move as $`R_{\mathrm{nnn}}`$. As an example, we present the rate calculations for a deep disorder-to-order step with $`\overline{\mu }_1=200`$ meV and $`\overline{\mu }_2=+600`$ meV. Each of the four rates is calculated from Eq. (6). The initial and final energies, $`U_I`$ and $`U_F`$, respectively, are calculated from Eq. (1). At $`\overline{\mu }_1`$, $`\mathrm{\Theta }0`$; therefore, at early times the surface is nearly empty. For adsorption, $`I`$ is approximated as having $`\mathrm{\Theta }=0`$, and $`F`$ contains only one occupied site. Thus, $`U_I=0`$, $`U_F=\overline{\mu }_2`$, and $`R_\mathrm{a}=\nu \mathrm{exp}[(\mathrm{\Delta }_\mathrm{a}(U_FU_I)/2)/k_\mathrm{B}T]=\nu \mathrm{exp}[0.0]`$. For desorption, $`I`$ must contain at least one occupied site for the process to be possible; therefore, $`I`$ is approximated as having only one occupied site, and $`F`$ is approximated as an empty surface. Thus, $`U_I=\overline{\mu }_2`$, $`U_F=0`$, and $`R_\mathrm{d}=\nu \mathrm{exp}[(\mathrm{\Delta }_\mathrm{a}(U_FU_I)/2)/k_\mathrm{B}T]=\nu \mathrm{exp}[24.0]`$. For nearest-neighbor diffusion, $`I`$ is also approximated as having only one occupied site, and $`F`$ is one of the four possible configurations which would occur from a nearest-neighbor diffusion move. Thus, $`U_I=U_F=\overline{\mu }_2`$, and $`R_{\mathrm{nn}}=4\nu \mathrm{exp}[(\mathrm{\Delta }_{\mathrm{nn}}(U_FU_I)/2)/k_\mathrm{B}T]=4\nu \mathrm{exp}[4.0]`$, where the factor of 4 comes from adding the contributions from each of the four final states. Similarly, for next-nearest neighbor diffusion, $`U_I=U_F=\overline{\mu }_2`$, and $`R_{\mathrm{nnn}}=4\nu \mathrm{exp}[(\mathrm{\Delta }_{\mathrm{nnn}}(U_FU_I)/2)/k_\mathrm{B}T]=4\nu \mathrm{exp}[8.0]`$. The initial rates for all four potential steps are summarized in Table 1. Note that the rates change significantly during the course of the time evolution following each step.

We now discuss each of the four potential steps. First we consider a disorder-to-order step with $`\overline{\mu }_1=200`$ meV and $`\overline{\mu }_2=+600`$ meV. At $`\overline{\mu }_1`$, $`\mathrm{\Theta }0`$, and the final potential $`\overline{\mu }_2`$ is deep into the ordered phase, approximately $`420`$ meV on the positive side of the transition at $`\overline{\mu }_c180`$ meV. Figure 7(a) shows time series for both $`\mathrm{\Theta }`$ and $`|\mathrm{\Theta }_\mathrm{S}|`$. For deep steps like these, the desorption rate is negligible (see Table 1). Thus, the adsorption dynamics are essentially described by the process of random sequential adsorption with diffusion (RSAD) . The coverage quickly reaches the jamming coverage for random sequential adsorption of hard squares without diffusion, $`\mathrm{\Theta }_\mathrm{J}0.364`$ 666 Precise values for $`\mathrm{\Theta }_\mathrm{J}`$ in the hard-square model, 0.364 13(1) and 0.364 133 0(5), have been obtained by MC simulations  and series expansions , respectively.. This is only slightly less than the critical coverage, $`\mathrm{\Theta }_c0.368`$, which is reached in approximately 37 MCSS. During this time, small, uncorrelated domains of the ordered phase are formed on both sublattices. As the lattice fills up, the interfaces between these domains contain empty sites, but the nearest-neighbor exclusion generally prevents the adsorption of additional Br ions. This is followed by the growth of adsorbate domains on the two sublattices which is well described by dynamical scaling (see below). At this stage, the coverage can increase only where domain walls move together, opening a gap large enough to fit an additional Br. In the absence of single-particle diffusion, the domain walls can move only by a series of very improbable desorption/adsorption events. However, nearest-neighbor diffusion allows some Br ions to diffuse between adjacent domains, thus providing a much faster mechanism for interface motion. For deep steps with the barriers used here, interface motion is almost exclusively due to nearest-neighbor diffusion. Excluding very small fluctuations on short time scales, the coverage is always driven toward the equilibrium value, $`\mathrm{\Theta }_{\mathrm{eq}}=1/2`$. Thus, the total length of interfaces is steadily decreasing. Finally, selection of one or the other of the two degenerate sublattices takes place at late times. The phase-selection time is expected to depend on $`L`$.

In the dynamical scaling regime, the correlation length $`\xi `$ of the system should grow as a power law with time $`t`$ ,

$$\xi t^n.$$
(7)
In models with non-conserved order parameter such as the one studied here, the domain growth is driven solely by the tendency to reduce the amount of interface. For such systems (known as Model A in the nomenclature of critical dynamics ) $`n=1/2`$ . We have used two different measures for $`\xi `$. One common measure uses the width of the diffuse structure-factor peak corresponding to the ordered phase, here the $`(1/2,1/2)`$ or $`\stackrel{}{k}=(\pi ,\pi )`$ peak. We define the correlation length given by the structure factor as $`\xi _S=2\pi /k`$, where

$$k=\frac{_{|\stackrel{}{k}\stackrel{}{k}_{\pi ,\pi }|>0}^{}|\stackrel{}{k}\stackrel{}{k}_{\pi ,\pi }|S(|\stackrel{}{k}|)}{_{|\stackrel{}{k}\stackrel{}{k}_{\pi ,\pi }|>0}^{}S(|\stackrel{}{k}|)}.$$
(8)
Here $`S(|\stackrel{}{k}|)`$ is the circularly averaged $`S(\stackrel{}{k})`$ and $`_{|\stackrel{}{k}\stackrel{}{k}_{\pi ,\pi }|>0}^{}`$ runs over all allowed values of $`|\stackrel{}{k}||\stackrel{}{k}_{\pi ,\pi }|`$, such that $`|\stackrel{}{k}\stackrel{}{k}_{\pi ,\pi }|`$ is less than $`\pi `$ to avoid picking up contributions from the small diffuse scattering surrounding the $`\stackrel{}{k}=(0,0)`$ Bragg peak. Figure 8 shows three-dimensional views of the structure factor and typical system configurations at three different times following this deep step into the c$`(2\times 2)`$ ordered phase. The diagonal anisotropy in $`S(\stackrel{}{k})`$ at late times corresponds to domain walls oriented at 45 to the lattice directions, which are stable under diffusion. Figure 9 shows the circularly averaged structure factor at various times. At late times the structure factor is proportional to $`|\stackrel{}{k}\stackrel{}{k}_{\pi ,\pi }|^3`$. This asymptotic behavior, known as Porod’s law , indicates that the domain walls are microscopically thin.

For this system, we are also able to exploit the area covered by the empty domain walls as a second measure of $`\xi `$ , since a decrease in interfacial length is equivalent to an increase in the correlation length . We define this second measure as $`\xi _W=(\mathrm{\Theta }_{\mathrm{eq}}\mathrm{\Theta })^1`$, where $`\mathrm{\Theta }_{\mathrm{eq}}`$ is the equilibrium value of the coverage at $`\overline{\mu }_2`$ (here, $`\mathrm{\Theta }_{\mathrm{eq}}=1/2`$). The inset in Fig. 7(a) shows $`\xi _S`$ and $`\xi _W`$ as functions of $`t^{1/2}`$ for $`L=256`$. The agreement with the expected $`t^{1/2}`$ behavior is excellent. This behavior has also been found in previous studies of RSAD . For an infinite system, $`\xi `$ would continue to grow without bound; however, when $`\xi `$ reaches the same order of magnitude as $`L`$, the system enters the phase-selection regime in which one of the two degenerate ordered phases is randomly selected.

For the other disorder-to-order step simulation, we take $`\overline{\mu }_1=200`$ meV and $`\overline{\mu }_2=+250`$ meV. This step ends approximately $`70`$ meV on the positive side of the transition. Figure 7(b) shows the dynamic behavior following this shallow step into the ordered phase. As one would expect, the overall response is much slower than for the deep step. Although the behavior at early times is more complicated than for the deep step, $`\xi _S`$ does show $`t^{1/2}`$ domain growth at late times (see inset). The length scale evaluated from the coverage, $`\xi _W`$, is not an accurate measure for shallow steps in the time range we have simulated, since the interfaces are not as clearly defined as in the deep step.

As expected, the ordering dynamics are well described by the dynamical scaling prediction of Model A. Although dynamical scaling has yet to be experimentally measured in this system, it has been measured in other electrochemical adsorption systems , and we believe that $`t^{1/2}`$ dynamical scaling will be experimentally observed in the future for this system as well.

Potential steps that induce disordering are also of interest. For deep negative steps, $`\overline{\mu }_1`$=$`+600`$ meV and $`\overline{\mu }_2`$=$`200`$ meV, the disordering process is relatively uneventful, as shown in Fig. 10(a). Particles simply desorb at a roughly constant rate, leading to an essentially exponential relaxation to equilibrium.

Much more interesting behavior is seen after a shallow negative step, such as $`\overline{\mu }_1`$=$`+600`$ meV and $`\overline{\mu }_2`$=$`+100`$ meV, shown in Fig. 10(b). The simulation begins with all particles on sublattice $`A`$ and relaxes to the disordered phase at $`\overline{\mu }_2`$, where $`\mathrm{\Theta }1/4`$. Four different dynamical regimes are seen. In the first regime, particles simply desorb from the surface, so that $`\mathrm{d}\mathrm{\Theta }_\mathrm{S}/\mathrm{d}t2\mathrm{d}\mathrm{\Theta }/\mathrm{d}t`$. In the second regime, both desorption and diffusion contribute significantly to the disordering process. As more sites become vacant by desorption, small domains of the $`B`$ phase are created, both by lateral diffusion from sublattice $`A`$ and by re-adsorption from the solution. As a result, the diffuse scattering peak narrows. In the third regime, the net desorption slows down, and diffusion and to some extent adsorption/desorption are the only significant contributions to the disordering process. Here, $`\mathrm{d}\mathrm{\Theta }_\mathrm{S}/\mathrm{d}t`$ is largest. In the final regime, the two sublattices become approximately equally populated, yielding $`\mathrm{\Theta }_\mathrm{S}0`$. Similar behavior is seen in a dynamic mean-field theory .

The inset in Fig. 10(b) shows the correlation length obtained from the diffuse (1/2,1/2) scattering peak as a function of time for the step to +100 meV. This length can be interpreted as the typical distance over which the order-parameter fluctuations are correlated. Thus, at early times, when the Br are predominantly on sublattice $`A`$, the fluctuations onto sublattice $`B`$ are very small. We believe that the temporary increase in the correlation length obtained from the diffuse scattering is caused by the system passing near the critical point of the phase transition. The length scale of the fluctuations decreases again as the coverage falls below its critical value and the macroscopic order is lost.

For shallow disordering steps, time series for $`\mathrm{\Theta }`$ and $`\mathrm{\Theta }_\mathrm{S}`$ are quite sensitive to the choice of the adsorption/desorption barrier. Here we discuss how this dependence can be exploited to estimate both $`\mathrm{\Delta }_\mathrm{a}`$ and the rate constant $`\nu `$ from future experimental results. For simplicity, we assume that our choices of diffusion barriers are realistic. Neglecting this assumption does not affect our qualitative conclusions since the time scale for the disordering process is dominated by $`\mathrm{\Delta }_\mathrm{a}\mathrm{\Delta }_{\mathrm{nn}}`$. We first discuss $`\mathrm{\Theta }(t^{})`$, which depends on time only through the dimensionless rescaled time $`t^{}=t\nu \mathrm{exp}(\mathrm{\Delta }_\mathrm{a}/k_\mathrm{B}T)`$. See Fig. 11. By comparing experimental measurements of $`\mathrm{\Theta }(t)`$ with the rescaled simulation results, the product $`\nu \mathrm{exp}(\mathrm{\Delta }_\mathrm{a}/k_\mathrm{B}T)`$ can be estimated. However, knowledge of $`\mathrm{\Theta }(t)`$ alone is not enough to independently estimate both $`\nu `$ and $`\mathrm{\Delta }_\mathrm{a}`$. To estimate both parameters, knowledge of $`|\mathrm{\Theta }_\mathrm{S}(t)|`$ is also required. The dependence of $`|\mathrm{\Theta }_\mathrm{S}|`$ on $`\mathrm{\Delta }_\mathrm{a}`$ can be seen in Fig. 11, but the dependence is more readily seen in a parametric plot of $`|\mathrm{\Theta }_\mathrm{S}(\mathrm{\Theta })|`$, which is independent of $`\nu `$. Such a plot is shown in Fig. 12. In this form, the shape dependence on $`\mathrm{\Delta }_\mathrm{a}`$ is more obvious. Also, the result that $`\mathrm{d}|\mathrm{\Theta }_\mathrm{S}|/\mathrm{d}t=2\mathrm{d}\mathrm{\Theta }/\mathrm{d}t`$ at short times is easily seen in the parametric plot. By plotting $`|\mathrm{\Theta }_\mathrm{S}(\mathrm{\Theta })|`$ from experiments and simulations together, it should be possible to estimate $`\mathrm{\Delta }_\mathrm{a}`$, assuming that the diffusion barriers used are realistic. Once $`\mathrm{\Delta }_\mathrm{a}`$ has been estimated, the product $`\nu \mathrm{exp}(\mathrm{\Delta }_\mathrm{a}/k_\mathrm{B}T)`$, known from comparing the time scales of the experimental and simulated processes, can be used to estimate $`\nu `$. It is our hope that experimental measurement of both $`\mathrm{\Theta }(t)`$ and $`|\mathrm{\Theta }_\mathrm{S}(t)|`$ will be possible in the future for this shallow disordering process.

## 6 Simulated Cyclic Voltammograms

Next we use the dynamic Monte Carlo method described in Sec. 4 to study the dynamic response of the adlayer to a continuously varying potential, as realized in a CV experiment. The system displays a phase lag between the ordering and disordering processes and the time dependent $`\overline{\mu }`$. For these simulated CVs, $`\overline{\mu }`$ is linearly ramped from $`\overline{\mu }_1`$ to $`\overline{\mu }_2`$ with a sweep rate $`\rho `$ and then ramped back to $`\overline{\mu }_1`$ at the same rate. A wide range of sweep rates were simulated, from 10-1 to $`3\times 10^4`$ meV/MCSS. Figure 13 shows $`\mathrm{\Theta }`$ as a function of $`\overline{\mu }`$ and is analogous to the familiar hysteresis plots for magnetic systems, showing the magnetization as a function of the applied magnetic field . The net flux of Br- into the adlayer is proportional to the voltammetric current. Figure 14 shows $`\rho ^1\mathrm{d}\mathrm{\Theta }/\mathrm{d}t`$ as a function of $`\overline{\mu }`$ for the three slowest sweep rates from Fig. 13, along with the quasi-equilibrium CV from Fig. 3. For the faster sweep rates, we used $`L=256`$, while for the slowest sweep rates, $`L=128`$ was used so that simulations could be performed in a reasonable amount of time (about 30 days of CPU time per cycle for $`\rho =3\times 10^4`$ meV/MCSS on a Pentium III 500 MHz processor). To reduce the effects of statistical noise, the derivatives were calculated by the Savitzky-Golay smoothing method using a quadratic fit to the $`\mathrm{\Theta }(t)`$ data over a $`\overline{\mu }`$ range of about 60 meV . This accounts for the reduced peak height for $`\rho =0`$, compared to Fig. 3.

As expected, the CVs approach the quasi-equilibrium result as $`\rho `$ approaches zero, and the phase lag between the response and the potential depends on the sweep rate. For non-zero sweep rates, there is asymmetry between the positive-going (ordering) and negative-going (disordering) scans caused by the differences between the phase-ordering and disordering processes discussed in Sec. 5. We denote the position of the sharp peak in $`\mathrm{d}\mathrm{\Theta }/\mathrm{d}t`$ for the positive-going scan as $`\overline{\mu }_\mathrm{p}`$ and for the negative-going scan as $`\overline{\mu }_\mathrm{n}`$. Figure. 15(a) shows $`\overline{\mu }_\mathrm{p}\overline{\mu }_c`$ and $`|\overline{\mu }_\mathrm{n}\overline{\mu }_c|`$ as functions of $`\rho `$. The peak can be distinguished in the positive-going scans for all of the simulated sweep rates but is only visible in negative-going scans for the three slowest rates. We therefore report the separation between the positive-going and negative-going peak positions, $`\overline{\mu }_\mathrm{p}\overline{\mu }_\mathrm{n}`$, for only these slowest rates, shown in Fig. 15(b). No simple behavior is seen, but by examining experimental peak separations , it is obvious that our slowest sweep rate is much faster than typical experimental sweep rates in the range of 1-10 mV/s. At the slowest simulated sweep rate, our peak separation is about 60 mV (in experimental electrode-potential units). From the shape of the curve it is tempting to suggest that by decreasing the sweep rate by another order of magnitude, we might obtain peak separations in the experimental range of about 10 mV.

## 7 Summary and Conclusion

We have studied the static and dynamic behavior of Br electrosorption on the Ag(100) single-crystal surface by Monte Carlo simulation. We have made comparison with equilibrium experiments to determine the strength of the lateral dipole-dipole repulsion and the electrosorption valency. The simulations reproduce the experimentally observed second-order phase transition between a disordered adlayer and a doubly degenerate c$`(2\times 2)`$ ordered phase at room temperature. As was also found previously by Koper , the simulations reproduce the broad shoulder seen in experimental cyclic voltammograms. We explain the shoulder by small configurational fluctuations in the disordered phase that locally resemble phases with $`\mathrm{\Theta }=1/4`$, one of which is present for this model in equilibrium at very low temperatures. Stronger effects of this low-temperature equilibrium phase might be observed in the behavior of other halides on Ag(100).

We also studied the dynamics of phase ordering and disordering processes by stepping the potential between the disordered and ordered phases. The ordering process is well described by $`t^{1/2}`$ dynamical scaling, consistent with growth of domains of the two degenerate c($`2\times 2`$) ordered phases. For disordering, no such competition between degenerate domains exists, and the process is well described by random desorption with diffusion. Study of such disordering processes may provide a way to estimate the adsorption/desorption barrier and obtain a relation between MC time and physical time.

In simulations of CVs for a wide range of potential sweep rates we found a pronounced asymmetry between the positive-going and negative-going scans. This asymmetry is caused by the differences between the ordering and disordering kinetics. No simple behavior is seen for the peak separation as a function of the sweep rate, but comparison with experiments suggests that our slowest sweep rate is considerably faster than experimental sweep rates which are in the range of 1-10 mV/s.

We leave comparison of simulated and experimental dynamics for future work when it hopefully becomes possible to experimentally study the dynamics of this system in full detail.

## Acknowledgments

We thank J. X. Wang, Th. Wandlowski, and B. M. Ocko for providing access to their experimental data and for useful discussions and correspondence. We also thank O. M. Magnussen for suggesting the importance of studying order-to-disorder potential steps and V. Privman for correspondence on random sequential adsorption. Comments on the manuscript by M. A. Novotny and Th. Wandlowski are also gratefully acknowledged. Supported in part by NSF grants No. DMR-9634873 and DMR-9981815, by Florida State University through the Center for Materials Research and Technology, the Supercomputer Computations Research Institute (U.S. Department of Energy contract No. DE-FC05-85FR2500), and the School of Computational Science and Information Technology.

## Appendix

In this Appendix we discuss the nonlinear fitting of simulated adsorption isotherms to the experimental data from Ref. .

Fixing a value of $`\varphi _{\mathrm{nnn}}`$, we generated a room-temperature isotherm by MC simulation, in this Appendix denoted by $`\mathrm{\Theta }^{\mathrm{sim}}(\overline{\mu };\varphi _{\mathrm{nnn}})`$. We compared this isotherm to three experimental coverage isotherms \[$`\mathrm{\Theta }^{\mathrm{exp}}(E;C_1)`$, $`\mathrm{\Theta }^{\mathrm{exp}}(E;C_2)`$, and $`\mathrm{\Theta }^{\mathrm{exp}}(E;C_3)`$\] at three different Br- concentrations \[$`C_1`$, $`C_2`$, and $`C_3`$\]. For convenience, we set the reference concentration $`C_0=1`$ mM, corresponding to the units used for the experimental KBr concentrations. Equation (2) was inverted to obtain $`E(\overline{\mu };\overline{\mu }_0,\gamma ;C)`$, assuming trial values of the parameters $`\overline{\mu }_0`$ and $`\gamma `$, and the $`\overline{\mu }`$-axis was rescaled into an $`E`$-axis for each concentration. Thus, we obtained three simulated approximations to the individual experimental isotherms \[$`\mathrm{\Theta }^{\mathrm{sim}}(E;\varphi _{\mathrm{nnn}},\gamma ,\overline{\mu }_0;C_1)`$, $`\mathrm{\Theta }^{\mathrm{sim}}(E;\varphi _{\mathrm{nnn}},\gamma ,\overline{\mu }_0;C_2)`$, and $`\mathrm{\Theta }^{\mathrm{sim}}(E;\varphi _{\mathrm{nnn}},\gamma ,\overline{\mu }_0;C_3)`$\].

In order to obtain the “best” values of the parameters $`\varphi _{\mathrm{nnn}}`$, $`\gamma `$, and $`\overline{\mu }_0`$, we defined the simultaneous squared difference between the simulation at $`\varphi _{\mathrm{nnn}}`$ and the experimental isotherms,

$$\chi ^2(\varphi _{\mathrm{nnn}},\gamma ,\overline{\mu }_0)=\underset{k=1}{\overset{3}{}}\underset{l=1}{\overset{l_{\mathrm{max}}(k)}{}}[\mathrm{\Theta }^{\mathrm{sim}}(E_l^k;\overline{\mu }_0,\gamma ,\varphi _{\mathrm{nnn}};C_k)\mathrm{\Theta }^{\mathrm{exp}}(E_l^k;C_k)]^2.$$
(9)
Here $`C_k`$ is the $`k^{th}`$ experimental concentration value, $`E_l^k`$ is the $`l^{th}`$ value of the electrode potential for the $`k^{th}`$ experimental isotherm, $`_l`$ is the sum over all experimental data points in an isotherm, and $`_k`$ is the sum over the three experimental concentrations. The $`\chi ^2`$ calculation requires that the simulation and experimental isotherms be evaluated at the same $`E`$ values. This was accomplished by linear interpolation of the simulation data.

The parameters $`\varphi _{\mathrm{nnn}}`$, $`\gamma `$, and $`\overline{\mu }_0`$ can, in principle, be found by minimizing $`\chi ^2`$ with respect to all three parameters. However, generating a simulated isotherm for a particular value of $`\varphi _{\mathrm{nnn}}`$ is too computationally intensive to allow simultaneous minimization of $`\chi ^2`$ with respect to all three parameters. We therefore minimized $`\chi ^2`$ in two steps, which saves time and allows better control over the fitting process. First, we generated a family of simulation isotherms for $`\varphi _{\mathrm{nnn}}`$ between $``$98 and $``$10 meV at intervals of 2 meV with a $`\overline{\mu }`$ resolution of 50 meV over the electrochemical potential range $``$200 to $`+1000`$ meV. For each of these simulated isotherms, we minimized $`\chi ^2`$ with respect to $`\gamma `$ and $`\overline{\mu }_0`$ and denoted the minimum error obtained by $`\chi _{\mathrm{min}}^2(\varphi _{\mathrm{nnn}},\gamma ^{\mathrm{min}},\overline{\mu }_0^{\mathrm{min}})`$, where $`\gamma ^{\mathrm{min}}(\varphi _{\mathrm{nnn}})`$ and $`\overline{\mu }_0^{\mathrm{min}}(\varphi _{\mathrm{nnn}})`$ are the values of $`\gamma `$ and $`\overline{\mu }_0`$ which minimize $`\chi ^2`$ for a particular value of $`\varphi _{\mathrm{nnn}}`$. In the second step, we minimized $`\chi _{\mathrm{min}}^2(\varphi _{\mathrm{nnn}},\gamma ^{\mathrm{min}},\overline{\mu }_0^{\mathrm{min}})`$ with respect to $`\varphi _{\mathrm{nnn}}`$. The final parameter results are $`\varphi _{\mathrm{nnn}}=26\pm 2`$ meV, $`\gamma =0.73\pm 0.03`$, and $`\overline{\mu }_0=0.67\pm 0.03`$. These errors are estimated from the $`\varphi _{\mathrm{nnn}}`$ resolution and the differences in best-fit parameter values around $`\varphi _{\mathrm{nnn}}=26`$ meV. The results of the fit are shown in Fig. 1.


format: text

text:
# Statistical Mechanics of an Urban Cross : A Solvable Model

\[

## Abstract

We propose a model for the intersection of two urban streets. The traffic status of the crossroads is controlled by a set of traffic lights which periodically switch to red and green with a total period of$`T`$. Two different types of crossroads are discussed. The first one describes the intersection of two one-way streets, while the second type models the intersection of a two-way street with an one-way street. We assume that the vehicles approach the crossroads with constant rates in time which are taken as the model parameters. We optimize the traffic flow at the crossroads by minimizing the total waiting time of the vehicles per cycle of the traffic light. This leads to the determination of the optimum green-time allocated to each phase.

\]

Over the past decade, the majority of the studies in the inter-disciplinary field ” physics of traffic  ” have been allocated to the highway traffic  ( for a review see Ref.  ) and less attention has been paid to the urban traffic. The simulation of traffic flow in a large-sized city is a formidable task and many degrees of freedom have to be involved (see e.g. ). In practice, streets of a city form a network of junctions that are linked togetether. Each crossroads receives demands (vehicles attempting to pass the cross ) and at each crossroads there exists a traffic light which, with some certain programming, controls the transportation. The first model introduced by statistical physicists for the description of the city network which is known as the BML model in the literature, uses a deterministic cellular automata framework  and predicts a sort of phase transition from free-flow to a jammed state, but the model assumptions are too simple to be taken seriously for practical purposes. Later, a more serious model was introduced by Schadschneider and Chowdhury . This model combines the BML mode together with the Nagel-Schreckenberg model of highway traffic . Despite its nice formulation, there are still a lot of simplifications which prevent it from being an effective  and applicable model to a city network. The empirical mechanism by which the traffic lights are controlled is generally divided into two distinct methods : fixed time and real-time. In the fixed-time method, a fixed value of time is allocated to the traffic light as well as its sub-phase times. In the real-time method, which is becoming increasingly popular in great cities, the ensemble of crossroads are intelligently controlled by a central controller. The control mechanism is usually based on the concept of producing a kind of green waves  between the crossroads. These waves interact with each other, and, if the passing demands are of high values, the green waves may have destructive effects on each other, and, hence, the concept of green wave may fail to be the final solution for optimizing the overall flow.
In this paper we aim to analyze a single crossroads in detail in order to find a better insight to the problem of optimizing the total flow in cities. Our single crossroads is the result of the intersection of two perpendicular streets. In their simplest structure, these streets can each direct a one-way traffic flow. With no loss of generality, we take them as one-way South to North (S-N) and West to East (W-E) streets. Cars arrive at the south and the west entrances of the crossroads. Here we assume the arrival rates of the cars , i.e. the number of cars reaching the crossroads per second, are constant in time. Although from everyday driving experiences in cities, we know that these rates have fluctuations in the course of time; yet in definite time intervals , the assumption of the constant arrival rates could be justified at least on a average level. We take the arrival rates to be $`\alpha _1`$ (for S-N cars ) and $`\alpha _2`$ (for the W-E cars) respectively. Also we denote the passing-rate of cars ( number of cars passing the crossroads in the unit of time during the green-phase ) by $`\beta _1`$ and $`\beta _2`$. The period of the traffic lights is taken to be a definite value $`T`$ which is assumed to remain constant in the time interval under consideration. The starting time of each cycle of the traffic light is at the moment which the light turns green for the S-N street. The S-N light remains green until $`T_1`$. At $`T_1`$ the traffic lights turn red for the S-N street and simultaneously changes to green for the W-E street. This is the beginning of the second phase which continues from $`T_1`$ to $`T`$ (end of one cycle). During Phase I ( $`0tT_1`$ ), the S-N cars can pass the crossroads northward and W-E cars are stopped for the red light. Over Phase II ( $`T_1tT`$ ), the S-N cars must wait behind the red light whilst the W-E cars are passing through the crossroads.

Now the basic question is how should traffic engineers adjust the value of $`T_1`$ in order to optimize the traffic flow throughout the intersection? Regarding the fact that drivers are identical human beings and with the assumption that the number of passengers of each car takes an equal average value for each direction, the optimization task is realized by minimizing the total waiting time of the cars per cycle of the traffic lights. For this purpose, we introduce two quantities $`N_1`$ and $`N_2`$ which represent the number of cars stopping at the red lights in the red phases of the S-N and W-E streets respectively. Clearly $`N_1`$ and $`N_2`$ are functions of time and, in general, are divided into different lanes on the streets. For the first cycle of the steady time-interval ( in which the arrival rate are taken to be constant ) the following equations are easily written down:

$$N_1(T_1)=[N_1(0)+(\alpha _1\beta _1)T_1]\theta $$
(1)
$$N_2(T_1)=N_2(0)\theta +\alpha _2T_1$$
(2)
$$N_1(T)=[N_1(0)+(\alpha _1\beta _1)T_1]\theta +\alpha _1(TT_1)$$
(3)
$$N_2(T)=[N_2(0)\theta +\alpha _2T\beta _2(TT_1)]\theta $$
(4)
The $`\theta `$ symbols ensures the positiveness of the quantities in the brackets, i.e. the value of $`\theta `$ is one if the quantity in the bracket is positive and zero elsewhere. This limitation is dictated to us since, by definition, the quantities $`N_1`$ and $`N_2`$ can only take positive values.

We now define the Total Waiting Time (TWT) of the vehicles per cycle of the traffic lights. It is the total time wasted by the vehicles during their stop in the red phases. The TWT is sum of the sub-waiting times of each direction. Denoting the TWT and the sub waiting times by $`T^W`$, $`T^{(W,1)}`$ and $`T^{(W,2)}`$ respectively, the following equations could be written for the $`n`$-th cycle

$$T_{nn+1}^{(W,1)}=\frac{1}{2}\alpha _1(TT_1)^2+N_1(nT+T_1)(TT_1)$$
(5)
$$T_{nn+1}^{(W,2)}=N_2(nT)T_1+\frac{1}{2}\alpha _2T_1^2$$
(6)
In the above equations, one must take into account the fact that during the red phase of a street, not all the cars reaching the crossroads wait equally. The cars reaching the intersection just after the light goes red ( unlucky drivers) have the most waiting time, while those reaching the crossroads just before the lights turns green (lucky drivers) spend no time waiting. This leads to an average factor $`\frac{1}{2}`$ in the above equations. The value of $`N_2(0)`$ refers to the number of the stopped cars in the W-E street just after the instant when the light turns red for the W-E cars (in the first cycle). Here two distinct situations are identified. In the first case, $`N_2(0)`$ is nearly zero which describes the situation when the W-E street is not crowded. The other case is $`N_2(0)0`$ which corresponds to the situation of heavy traffic in the W-E street. This means that even in the final stages of the green phase for the W-E direction, still there is a considerable number of W-E cars desiring to pass the cross eastwards.The same scenario goes with the value of $`N_1(T_1)`$, and, as we will see in what follows, these different cases lead to different prescriptions of the optimization. Before investigating these cases in details, it should be mentioned that for the arbitrary values of the rates we do not have a periodical structure in the values of the $`N_1`$ and $`N_2`$. This can be easily verified by comparing $`N_1(T)`$ and $`N_2(T)`$ with their values at $`t=0`$ which implies that, in general, one has $`N_i(0)N_i(T)`$ for $`i=1,2`$. Now let us discuss the situation where $`N_1(T_1)=N_2(0)=0.`$ This refers to a light traffic flow in both directions during the first cycle of the interval. One can easily check that if the following stability conditions $`\alpha _2T\beta _2(TT_1)0`$ and $`\alpha _1T\beta _1T_10`$ hold, Then we always have the similar conditions in the next cycles $`N_1(nT+T_1)=N_2(nT)=0`$ where $`n`$ denotes the number of the traffic light cycle. In this particular case which here after will be called State I, the TWT is independent of the cycle-number, and minimizing the TWT with respect to $`T_1`$ leads to the following equation:

$$T_1=\frac{\alpha _1}{\alpha _1+\alpha _2}T$$
(7)
this answer must be consistent with the stability conditions which leads to the constraints $`\beta _2\alpha _1+\alpha _2`$ and $`\beta _1\alpha _1+\alpha _2`$ among the rates. We now investigate a totally different situation i.e. a crowded crossroads. Supposing that the first cycle is characterized by the conditions $`N_2(0),N_1(T_1)>0`$. one could easily verify that provided the relations $`\alpha _1>\beta _1`$ and $`\alpha _2T>\beta _2(TT_1)`$ hold, we have a quasi-stationary condition in the next cycles i.e. $`N_1(nT+T_1),N_2(nT)>0`$. In sharp contrast to the State I, (this state which will be referred to as the State II), the values of $`N_1`$ and $`N_2`$ would be functions of the cycle number. It can be shown that in the large cycle-number limit one has $`T_{nn+1}^{(W,1)}n(TT_1)(\alpha _1T\beta _1T_1)`$ and $`T_{nn+1}^{(W,2)}nT_1(\alpha _2T\beta _2(TT_1))`$. We now minimize the TWT with respect to $`T_1`$ which leads to the following equation:

$$T_1=\frac{\beta _1+\beta _2+\alpha _1\alpha _2}{2(\beta _1+\beta _2)}T$$
(8)
the consistency of this solution with stability conditions yields the following constraints: $`\alpha _1>\beta _1`$ and $`\alpha _1\beta _2+2\alpha _2\beta _1+\alpha _2\beta _2\beta _1\beta _2+\beta _2^2`$. Also positiveness of $`T_1`$ itself imposes the extra restriction $`\alpha _1\alpha _2<\beta _1+\beta _2`$. Next we consider the situation (State III) where in the first cycle of the traffic light, one has $`N_1(T_1)=0`$ but $`N_2(0)>0`$ this correspond to the situation in which the S-N street has a light traffic flow while the W-E street has a heavy one. The conditions for having these situations repeated in the next cycles are $`\alpha _1T\beta _1T_10`$ and $`\alpha _2T\beta _2(TT_1)0`$. In the large cycle-number limit, the minimizing of the TWT leads to the value

$$T_1=\frac{\beta _2\alpha _2}{2\beta _2}T$$
(9)
It could be easily checked that the above solution is inconsistent with the stability conditions, and, hence, is not acceptable as an optimization of traffic light. This inconsistency will remain in the final state ( State IV ) which is characterized by $`N_1(nT+T_1)>0`$ and $`N_2(nT)=0`$. In this case, the solution obtained by minimizing the TWT is again inconsistent with the stability conditions. In these last two states, one must numerically find the consistent minima. We next consider another more usual configuration of a crossroads. Here we let the vehicles move in the both directions on the S-N street but still the vehicles in the W-E street are restricted to move eastwards. This situation describes a one-way to two-way urban intersection. Here each cycle consists of three phases. In the first phase which lasts for $`0tT_1`$, the traffic light is green for the S-N cars and red for the other two directions. During the second phase which starts at $`T_1`$ and finishes at $`T_2`$, the traffic light is green for the N-S cars and red for the other two directions. In the final phase, which last for $`T_2tT`$, the traffic light remains green for the W-E cars and red for the N-S as well as S-N directions. The entrance rates are taken to be $`\alpha _1,\alpha _2`$ and $`\alpha _3`$ and we denote the passing rate by $`\beta _1,\beta _2`$ and $`\beta _3`$ for each direction respectively.The starting time of the cycles is chosen to be the moment at which the traffic light turns green for the S-N direction. Similar equations for the Vehicle-Numbers $`N_1,N_2`$ and $`N_3`$ could be written down and one in principle can evaluate the TWT in terms of these quantities. The exact form of the TWT strongly depends on the positiveness of the vehicle-numbers just after the traffic lights goes red for the respective direction. In the case under consideration, eight different possibilities can occur. Let us only discuss the most probable one which corresponds to the the following situation: $`N_1(nT+T_1)=N_2(nT+T_2)=N_3(nT)=0`$ In this case the traffic status is light in each of the three directions. Following the same discussion made for the previous crossroads, it could be easily verified that the triple stability condition for the validity of the above assumptions are as follows: $`\alpha _1T\beta _1T_10`$, $`\alpha _2T\beta _2(T_2T_1)0`$ and $`\alpha _3T\beta _3(TT_2)`$. It can be shown that minimizing the TWT with respect to $`T_1`$ and $`T_2`$

leads to the following fixation of $`T_1,T_2`$

$$T_1=T\frac{\alpha _1(\alpha _2+\alpha _3)\alpha _2\alpha _3}{\alpha _1(\alpha _2+\alpha _3)+\alpha _2\alpha _3}$$
(10)
$$T_2=T\frac{2\alpha _1\alpha _2}{\alpha _1(\alpha _2+\alpha _3)+\alpha _2\alpha _3}$$
(11)
One directly observes that in the symmetric case of equal arrival rates, $`T_1`$ and $`T_2`$ take the expected values $`\frac{T}{3}`$ and $`\frac{2T}{3}`$ respectively. The other point which must be mentioned is that in this traffic status, $`T_1`$ and $`T_2`$ do not depend on the passing rates, and are solely determined by the arrival rates. The remaining seven status will be discussed elsewhere.
For comparison of our model to the empirical data, a time-series analysis on two of Tehran intersections were carried out. The central part of Tehran is under control of the SCATS (Sydney Coordinated Adaptive Traffic System) that is an intelligent traffic controller system. The strategy followed by most of urban traffic control systems is based on establishing green-wave along the major streets of cities. We considered two different intersections. The first one which is located in Tehran downtown connects Valiasr St. to Takhtejamshid St. Both of these are one-way and major streets. The other crossroads connects Abbasabad St.(one-way) to Mahnaz St.(one-way) .In contrast to the previous case, here the first street is a major while the second street is a minor one. The data set is provided by magnetic counting loops which are installed just before the pedestrian-line of each crossroads. The data were collected on second of July, 2000. The three following figures show a time-series of data for our major-to-minor crossroads.

Fig. 1: The cycle-ratio of green-time of Abbasabad street to the total period of the traffic light cycle. The numbers on the horizontal axes are the cycle number.

Fig. 2: The cycle-ratio of the passed vehicles (during the green phases) of Abassabad street to Mahnaz street of each traffic light cycle. The numbers on the horizontal axis denote the cycle number.

Fig. 3: the above graph denotes the total number of passed vehicles (from both streets) during fifteen-minute time intervals. The lower graph shows the number of passed vehicles form Mahnaz street during fifteen-minute intervals.
As seen from the graphs, the cycle-ratio of the green time allocated to Abbasabad str. strongly fluctuates due to demand fluctuations received by the cross. In order to have a rough estimation of our model parameters, we considered a two-hour time interval between 12:30-14:30 . We approximated the ratio of $`\frac{\alpha _1}{\alpha _2}`$ by the ration of the total number of vehicles which have passed through Abbasabad Str. to that number for Mahnaz Street. This yields the value $`\frac{\alpha _1}{\alpha _2}=0.35`$. Putting this value into equation (7) yields $`\frac{T_1}{T}=0.74`$. On the other hand the averaged value of the cycle-ratio of the green time of Abbasabad street over the two-hour period leads to result $`\frac{T_1}{T}=0.64`$. The following figures belong to our major-to-major crossroads.

Fig. 4: the cycle-ratio of green time allocated to Takhtejamshid Street. The numbers on the horizontal axis denote the cycle number.

Fig. 5: cycle-ratio of the passed vehicles from Takhtejamshid street to Valiasr street.

Fig. 6: the above graph denotes the total number of passed vehicles (from both streets) during fifteen-minute time intervals. The lower graph shows the number of passed vehicles form Takhtejamshid str. during fifteen-minute intervals.
Here we focused on the interval 13:30 - 15:30. The averaged value of $`\frac{T_1}{T}`$ is 0.52 while the same amount evaluated from eq.(7) is 0.56.

In conclusion, we have developed a prescription for the traffic-light programming at a single urban crossroads. The method is based on minimizing the TWT of cars stopping in the red phases of the traffic light. The value of $`T`$ should be so tuned that during the green phases, time-headway is less than a certain value. In fact, in the model the values of the passing rate refer to the maximum capacity of cross ( maximum number of cars passing the cross in the unit of time ), which is plausible if the $`T`$ is appropriately adjusted with the congestion of the crossroads. Optimizing the traffic at each crossroads is the stating point of the more comprehensive problem of city network. Nevertheless, our model best suits those marginal intersections of cities where the effect of the other crossroads is suppressed.

Acknowledgments: M.E.F is grateful to Tehran Traffic Control Company  for the data support. We would like to express our gratitude to R. Sorfleet, M.Arab salmani, M. Salmani and M. Khoshechin for useful helps.


text:
# STRANGE QUARKS, MAGNETIC MOMENTS AND SU(3)

ERRATUM: The preprint listed as hep-ph/0007271 had a sign error which rendered false most of the conclusions drawn.


text:
# Tutorial on quasi-sparse eigenvector diagonalization

## 1 Introduction

Quasi-sparse eigenvector (QSE) diagonalization is a new computational method which finds the low-lying eigenvalues and eigenvectors for a general quantum field Hamiltonian . It is able to handle the exponential increase in the size of Fock space for large systems by exploiting the sparsity of the Hamiltonian. QSE diagonalization can even be applied directly to infinite-dimensional systems. The method is most effective when the splitting between low-lying eigenvalues is not too small compared to the size of the off-diagonal Hamiltonian matrix entries. In such cases the low-lying eigenvectors are quasi-sparse, meaning that the vector is dominated by a small fraction of its largest components. The QSE algorithm can then be applied to the Hamiltonian $`H`$ using the following steps:

1. Select a subset of orthonormal basis vectors $`\{e_{i_1},\mathrm{},e_{i_n}\}`$ and call the corresponding subspace $`S`$.
2. Diagonalize $`H`$ restricted to $`S`$ and find one eigenvector $`v.`$
3. Sort the basis components $`e_{i_j}|v`$ according to their magnitude and remove the least important basis vectors.
4. Replace the discarded basis vectors by new basis vectors. These are selected at random from a pool of candidate basis vectors which are connected to the old basis vectors through non-vanishing matrix elements of $`H`$.
5. Redefine $`S`$ as the subspace spanned by the updated set of basis vectors and repeat steps 2 through 5.

If the subset of basis vectors is sufficiently large, the exact eigenvectors will be stable fixed points of the update process.

The purpose of this short tutorial is to provide additional information for those interested in writing or using codes that implement QSE methods. We discuss the calculation of Hamiltonian matrix elements in the momentum Fock space representation and the generation and selection of new basis states. We also include a simple program which applies the QSE method to find the ground state of an input matrix. Readers interested in an introductory overview of QSE diagonalization are encouraged to consult the references in .

## 2 Fock states and the Hamiltonian matrix

It is most effective to describe the details of a QSE calculation in the context of a specific example. In this discussion we consider $`\varphi ^4`$ theory in $`1+1`$ dimensions, the same example as in . Some of the techniques described here have been specifically optimized for scalar field theories. Other systems require a somewhat different set of tools and will be discussed in the future.

The Hamiltonian density for $`\varphi ^4`$ theory in $`1+1`$ dimensions has the form

$$=\frac{1}{2}\left(\frac{\varphi }{t}\right)^2+\frac{1}{2}\left(\frac{\varphi }{x}\right)^2+\frac{\mu ^2}{2}\varphi ^2+\frac{\lambda }{4!}\text{:}\varphi ^4\text{:},$$

where the normal ordering is with respect to the mass $`\mu `$. We consider the system in a periodic box of length $`2L`$. We then expand in momentum modes and reinterpret the problem as an equivalent Schrödinger equation . We will use a momentum cutoff scheme where the modes $`q_n`$ are restricted in momentum so that $`\left|n\right|N_{max}`$. The resulting Hamiltonian is

$`H`$ $`=\frac{1}{2}{\displaystyle \underset{n}{}}\frac{}{q_n}\frac{}{q_n}+\frac{1}{2}{\displaystyle \underset{n}{}}\left(\omega _n^2(\mu )\frac{\lambda b(\mu )}{4L}\right)q_nq_n`$ (1)
$`+\frac{\lambda }{4!2L}{\displaystyle \underset{n_1+n_2+n_3+n_4=0}{}}q_{n_1}q_{n_2}q_{n_3}q_{n_4},`$
where

$$\omega _n(\mu )=\sqrt{\frac{n^2\pi ^2}{L^2}+\mu ^2}$$
(2)
and $`b(\mu )`$ is the coefficient for the mass counterterm

$$b(\mu )=\underset{n}{}\frac{1}{2\omega _n(\mu )}.$$
(3)
It is convenient to split the Hamiltonian into free and interacting parts with respect to an arbitrary mass $`\mu ^{}`$:

$$H_{free}=\frac{1}{2}\underset{n}{}\frac{}{q_n}\frac{}{q_n}+\frac{1}{2}\underset{n}{}\omega _n^2(\mu ^{})q_nq_n,$$
(4)
$`H`$ $`=H_{free}+\frac{1}{2}{\displaystyle \underset{n}{}}\left(\mu ^2\mu ^2\frac{\lambda b(\mu )}{4L}\right)q_nq_n`$ (5)
$`+\frac{\lambda }{4!2L}{\displaystyle \underset{n_1+n_2+n_3+n_4=0}{}}q_{n_1}q_{n_2}q_{n_3}q_{n_4}.`$
$`\mu ^{}`$ will be used to define the basis states of our Fock space. Since $`H`$ is independent $`\mu ^{}`$, it is useful to perform calculations for different $`\mu ^{}`$ to obtain an estimate of the error. It is also helpful to find the range of values for $`\mu ^{}`$ which maximizes the quasi-sparsity of the eigenvectors. This can be used to improve the accuracy of the QSE calculation.

Let us define ladder operators with respect to $`\mu ^{}`$,

$`a_n(\mu ^{})`$ $`=\frac{1}{\sqrt{2\omega _n(\mu ^{})}}\left[q_n\omega _n(\mu ^{})+\frac{}{q_n}\right]`$ (6)
$`a_n^{}(\mu ^{})`$ $`=\frac{1}{\sqrt{2\omega _n(\mu ^{})}}\left[q_n\omega _n(\mu ^{})\frac{}{q_n}\right].`$ (7)
These satisfy the usual commutation relations,

$`[a_n,a_m^{}]`$ $`=\delta _{nm}`$ (8)
$`[a_n,a_m]`$ $`=[a_n^{},a_m^{}]=0.`$ (9)
The Hamiltonian can now be rewritten as

$`H`$ $`={\displaystyle \underset{n}{}}\omega _n(\mu ^{})a_n^{}a_n+\frac{1}{4}(\mu ^2\mu ^2\frac{\lambda b(\mu )}{4L}){\displaystyle \underset{n}{}}\frac{\left(a_n+a_n^{}\right)\left(a_n+a_n^{}\right)}{\omega _n(\mu ^{})}`$ (10)
$`+\frac{\lambda }{192L}{\displaystyle \underset{n_1+n_2+n_3+n_4=0}{}}\left[\frac{\left(a_{n_1}+a_{n_1}^{}\right)}{\sqrt{\omega _{n_1}(\mu ^{})}}\frac{\left(a_{n_2}+a_{n_2}^{}\right)}{\sqrt{\omega _{n_2}(\mu ^{})}}\frac{\left(a_{n_3}+a_{n_3}^{}\right)}{\sqrt{\omega _{n_3}(\mu ^{})}}\frac{\left(a_{n_4}+a_{n_4}^{}\right)}{\sqrt{\omega _{n_4}(\mu ^{})}}\right].`$
In (10) we have omitted constants contributing only to the vacuum energy.

We can represent any momentum-space Fock state as a string of occupation numbers,

$$|o_{N_{\mathrm{max}}},\mathrm{},o_{N_{\mathrm{max}}},$$
(11)
where

$$a_n^{}a_n|o_{N_{\mathrm{max}}},\mathrm{},o_{N_{\mathrm{max}}}=o_n|o_{N_{\mathrm{max}}},\mathrm{},o_{N_{\mathrm{max}}}.$$
(12)
From the usual ladder operator relations, it is straightforward to calculate the matrix element of $`H`$ between two arbitrary Fock states. In  auxiliary cutoffs were used to render the dimension of Fock space finite. This is in fact an unnecessary restriction and QSE diagonalization encounters no difficulties dealing with the full infinite-dimensional space.

## 3 New basis vectors

Aside from calculating matrix elements, the only other fundamental operation involving basis vectors is the generation of new basis vectors. As mentioned before, the new states should be connected to an old basis vector through non-vanishing matrix elements of $`H`$. Let us refer to the old basis vector as $`|e`$. For this example there are two types of terms in our interaction Hamiltonian, a quartic interaction

$$\underset{n_1,n_2,n_3}{}\left(a_{n_1}+a_{n_1}^{}\right)\left(a_{n_2}+a_{n_2}^{}\right)\left(a_{n_3}+a_{n_3}^{}\right)\left(a_{n_1n_2n_3}+a_{n_1+n_2+n_3}^{}\right),$$
(13)
and a quadratic interaction

$$\underset{n}{}\left(a_n+a_n^{}\right)\left(a_n+a_n^{}\right).$$
(14)
To produce a new vector from $`|e`$ we simply choose one of the possible operator monomials

$`a_{n_1}a_{n_2}a_{n_3}a_{n_1n_2n_3},a_{n_1}^{}a_{n_2}a_{n_3}a_{n_1n_2n_3},\mathrm{},`$ (15)
$`a_na_n,a_n^{}a_n^{},\mathrm{}`$
and let it act upon $`|e`$. Our experience is that the interactions involving the small momentum modes are generally more important than those for the large momentum modes, a signal that the ultraviolet divergences have been properly renormalized. For this reason it is best to arrange the selection probabilities such that the smaller values of $`\left|n_1\right|`$, $`\left|n_2\right|`$, $`\left|n_3\right|`$ and $`\left|n\right|`$ are chosen more often.

We note that the new basis vector selection will occasionally fail. Either the vector is zero due to an annihilation operator acting on an unoccupied state or the vector is already in our subset of basis vectors and therefore not new. In either case we simply select a new basis vector again. If the selection process fails repeatedly then a different old basis vector $`|e^{}`$ is used.

## 4 Sample code

We have considered two basic operations, the calculation of matrix elements and the generation of new basis states. The subroutines which perform these tasks will depend on the form of the quantum field Hamiltonian. The main structure of the QSE algorithm, however, is independent of the details of the Hamiltonian. We demonstrate the essential features of the algorithm with the following short MATLAB program. In this example we find the ground state of a finite matrix. To keep the example as simple as possible, we avoid calculating matrix elements and generating new basis states by assuming that the entire Hamiltonian matrix can be stored in memory.

niter = 30;

nactive = 100;

nretain = 80;

H = loadsparse(’samplematrix’);

vecs = \[1:nactive\];

for iter = 1:niter

\[v,d\] = eigs(H(vecs,vecs),’SR’,1);

weights1 = v.^2;

\[sortelements, sortorder\] = sort(-weights1);

vecs = vecs(sortorder(1:nretain));

wtsum1 = cumsum(weights1)/sum(weights1);

while (length(vecs) $`<`$nactive)

over1 = find(rand\*wtsum1(nretain) $`<`$wtsum1);

seedvec = vecs(over1(1));

candidates = find(H(:,seedvec));

candidates = candidates(candidates ~= seedvec);

weights2 = abs(H(candidates,seedvec));

wtsum2 = cumsum(weights2)/sum(weights2);

over2 = find(rand $`<`$wtsum2);

newvec = candidates(over2(1));

vecs = unique(\[vecs newvec\]);

end

end

save results.mat vecs v d

niter is the number of iterations of the algorithm. nactive is the number of basis states used to form the subspace $`S`$, and nretain is the number of basis states kept after removing the tail of the eigenvector distribution. loadsparse is a program which reads in a datafile called samplematrix. This data is stored in the matrix variable H. The remaining code can be broken down in terms of the five basic steps of the QSE algorithm.

Select a subset of orthonormal basis vectors $`\{e_{i_1},\mathrm{},e_{i_n}\}`$ and call the corresponding subspace $`S`$.

vecs = \[1:nactive\];

The variable vecs is a row vector which labels the basis vectors in the subspace $`S`$. In our example vecs is initialized to be the first nactive basis vectors. We are assuming that the ground state of the system has a non-zero overlap with at least one of the first nactive basis vectors. A more general initialization would be to choose a random subset of basis vectors.

Diagonalize $`H`$ restricted to $`S`$ and find one eigenvector $`v.`$

\[v,d\] = eigs(H(vecs,vecs),’SR’,1);

The function eigs is a sparse matrix Arnoldi diagonalization routine. It is being called with the parameters ’SR’,1 which tells eigs to find the eigenvalue and eigenvector with the smallest real part. The Hamiltonian submatrix corresponding with basis vectors in vecs is diagonalized. The eigenvector is stored in the column vector v and the eigenvalue is stored as the variable d.

Sort the basis components $`e_{i_j}|v`$ according to their magnitude and remove the least important basis vectors.

weights1 = v.^2;

\[sortelements, sortorder\] = sort(-weights1);

vecs = vecs(sortorder(1:nretain));

weights1 is a column vector which contains the squares of the components of the eigenvector v. sortelements is a list of the elements of -weights1 from smallest to greatest. sortorder is the corresponding list of locations for these elements. If for example weights1 = \[.2 .5 .3\], then sortorder = \[2 3 1\]. The list sortelements is not used. Only the nretain most important basis vectors are kept in vecs.

Replace the discarded basis vectors by new basis vectors. These are selected at random from a pool of candidate basis vectors which are connected to the old basis vectors through non-vanishing matrix elements of $`H`$.

wtsum1 = cumsum(weights1)/sum(weights1);

while (length(vecs) $`<`$nactive)

over1 = find(rand\*wtsum1(nretain) $`<`$wtsum1);

seedvec = vecs(over1(1));

candidates = find(H(:,seedvec));

candidates = candidates(candidates ~= seedvec);

weights2 = abs(H(candidates,seedvec));

wtsum2 = cumsum(weights2)/sum(weights2);

over2 = find(rand $`<`$wtsum2);

newvec = candidates(over2(1));

vecs = unique(\[vecs newvec\]);

end

cumsum is a function which takes the input list and creates another list consisting of partial sums. wtsum1 is therefore an ascending string of numbers between $`0`$ and 1 such that the difference between an entry and its preceding entry is proportional to the corresponding value of weights1. over1 is a location list of elements of wtsum1 which are greater than the product of a generated random number between 0 and 1 and wtsum1(nretain). over1(1$`)`$ is the first entry of over1. We note that over1(1$`)`$ must be less than or equal to nretain. seedvec is the basis vector which corresponds with over1(1). candidates is a column vector consisting of all basis vectors which have a non-zero matrix element with seedvec. Since the task is to find new basis vectors, seedvec is removed from candidates.

weights2 is a column vector containing the absolute values of the matrix transition elements from seedvec. wtsum2 is an ascending string between $`0`$ and 1 such that the difference between an entry and its preceding entry is proportional to the value of weights2. over2 is a location list of elements of wtsum2 greater than a generated random number. newvec is the basis vector which corresponds with over2(1). newvec is appended to the set vecs. The function unique sorts the input list and removes any repetitions. This loop is repeated until vecs has nactive elements.

Redefine $`S`$ as the subspace spanned by the updated set of basis vectors and repeat steps 2 through 5.

for iter = 1:niter

...

end

save results.mat vecs v d

The algorithm is iterated niter times and the final results are saved in the file results.mat. More examples of QSE diagonalization codes for various field theory systems written in MATLAB, Fortran, and/or C++ will be made available in the future.

## 5 Summary

In this short tutorial we have listed some supplemental material for those interested in writing or using codes that implement QSE methods. We have discussed the calculation of Hamiltonian matrix elements and the selection of new basis states. We have also presented a simple MATLAB program which applies the QSE method to find the ground state of a matrix.

The discussion here was limited to the basic QSE method. As presented in the Guangzhou workshop there have also been interesting new developments in QSE diagonalization regarding the technique of stochastic error correction. This technique will be presented in a forthcoming work .

## Acknowledgments

I thank my collaborators on the papers cited here and the organizers and participants of the International Workshop on Nonperturbative Methods and Lattice QCD in Guangzhou. Financial support was provided by the National Science Foundation.


text:
# Global Characteristics of Nucleus-Nucleus Collisions in an Ultrarelativistic Domain

## 1 Introduction

Global characteristics of nucleus-nucleus collisions available in the future experiments at high energy (the LHC energy scale) are differential distributions of total transverse energy flows, $`dE_T/d\eta `$, as well as electromagnetic energy and charged multiplicity ones, $`dE_\gamma /d\eta `$ and $`dN_{ch}/d\eta `$. Measurements of the global distributions allows one to establish quite general rules of nucleus-nucleus collision dynamics up to the highest collision energy frontier over a widest rapidity range and verify some critical predictions of quark-gluon plasma formation models in a sufficiently simple way.

In addition to the general physics interest the global event response provides an estimate of a collision impact parameter for correlated studies in specific reaction channels like jet and $`W`$, $`Z^0`$ production.

## 2 Jet Quenching Manifestation in Global Energy Flows

The dense matter formation (or quark-gluon plasma) in relativistic nuclear collisions is predicted by many models . One of the discussed features of such a state of a nuclear matter is energy losses of scattered partons in final state interactions with a dense nuclear matter called jet quenching . Among other effects originated by a jet quenching one may expect a significant modification in the differential distributions mentioned above. Indeed an indication is found for ultrarelativistic energy domain on the appearance of a wide bump in the interval $`2<\eta <2`$ over a pseudorapidity plateau of such distributions due to jet quenching . The bump existence problem rises a new question of principle – whether the asymptotic behaviour of the distributions established already at the lower energy scale will be broken in a new energy domain, or not. A calorimeter with a wide enough acceptance ($`5<\eta <5`$) gives a chance to obtain a definite answer to the question.

Our calculations has been done in the framework of the HIJING model  for the energy range from 1 up to 5 TeV/nucleon. In accordance to the HIJING already $``$80% of the total transverse energy flow are calculated by a perturbative QCD application. A minijet production (i. e. jets with $`p_t`$ 2 GeV) becomes a dominant feature of multiple processes in the considered energy region. Soft hadron-hadron processes are simulated as classical strings with kinks and valence quark ends following the FRITIOF model . The number of jets in inelastic nucleon-nucleon collisions is calculated in the eikonal approximation. The model provides the dependence of the parton structure function on a collision impact parameter.

The most important feature of the HIJING model is the energy losses $`dE/dx`$ of partons traversing a dense nuclear matter. High energy quark and gluon losses in a hot chromodynamic matter are estimated in . It was shown that the dominant mechanism of energy losses is a radiative one due to a gluon emission inside a narrow jet cone (Bethe-Heitler limit) in soft final state interactions (bremsstrahlung) in accordance with the relation used in the HIJING:

$$\frac{dE}{dx}\frac{C_F\alpha _s}{\pi }\mu ^2ln\frac{3ET}{2\mu ^2}(ln\frac{9E}{\pi ^3T}+\frac{3\pi ^2\alpha _s}{2\mu ^2}T^2)$$
(1)
where $`E`$ is the propagating parton energy, $`\mu _{D}^{}{}_{}{}^{1}1/gT`$ the screening constant (Debye cutoff scale), $`\alpha _s`$ the running coupling constant. In the HIJING model this mechanism is fulfilled as consecutive transmissions of parton energy fractions from one string configuration to another.

It is supposed that energy losses are proportional to a distance $`l`$ passed by a jet after a last interaction and occur only in the transverse direction within the nucleus radius. An interaction proceeds until a parton jet stays inside the considered volume and the jet energy exceeds a jet production energy limit. The space-time picture of the evolution of a parton shower is not considered in the HIJING and the dense matter formation is introduced by a phenomenological way.

In fig.1 differential transverse energy distributions $`dE_T/d\eta `$ are presented for minimum bias Pb-Pb collision at the CMS energy $`\sqrt{s_{NN}}`$ equal to 5, 3, 1, 0.5, 0.2, 0.1 TeV/nucleon. A distinctive feature of the distributions is the appearance of a central pseudorapidity bump (-2$`<\eta <`$2).

One may conclude that the bump becomes distinguishable over a plateau starting from an energy value larger than 3 TeV/nucleon. For lower energy values this effect is not so profound due to smaller value of parton energy losses in the nuclei. This fact induced by decreasing the Debye screening constant of nuclear matter at a lower energy.

We followed a jet quenching sensitivity to a mass number for lighter colliding nuclei. Energy losses of hard parton jets depend on a parton path in a dense matter like $`\mathrm{\Delta }E=ldE/dx`$. Therefore reduction of the colliding nucleus radii might lead to a bump reduction. This provides an additional verification of the jet quenching effect . Fig.2 shows that lead-lead collisions demonstrate maximum quenching dependence while in lighter ion cases the bump becomes less and less profound. It’s interesting to note that the shadowing induced modification of the parton distribution functions can be estimated by means of absolute value measurements of the transverse energy flow density for various ion species.

## 3 Impact Parameter Estimation

A secondary interaction effect as a jet quenching modifies only the central rapidity part leaving both fragmentation regions ($`3|\eta |5`$) practically unchanged. This circumstance can be used for estimation of a collision impact parameter. The average number of minijets produced in nucleus-nucleus collisions in the pseudorapidity interval $`\mathrm{\Delta }\eta `$ is related with a collision impact parameter $`b`$:

$$\overline{N}_{AA}(b,\sqrt{s},p_0)_{\mathrm{\Delta }\eta }=T_{AA}(b)\sigma _{jet}(\sqrt{s},p_0)_{\mathrm{\Delta }\eta }$$
(2)
where $`T_{AA}(b)`$ is the nuclear density overlap function of two colliding nuclei calculated with the assumption of the Wood-Saxon nuclear density distribution $`\rho (r)`$, $`\sigma _{jet}`$ the minijet production cross-section.

An average transverse energy flow in $`3|\eta |5`$ intervals is related with a collision impact parameter as

$$\overline{E}_T(b,\sqrt{s},p_0,3|\eta |5)=T_{AA}(b)\sigma _{jet}(\sqrt{s},p_0,3|\eta |5)<E_t>_{\mathrm{\Delta }\eta }$$
(3)
where $`E_t`$ is the mean transverse energy per a minijet. The experimental data shows that parton structure functions using to calculated $`\sigma _{jet}(\sqrt{s},p_0,3|\eta |5)<E_t>_{\mathrm{\Delta }\eta }`$ stay relevant until transverse energy per minijet $`E_t`$300 GeV. In the considered pseudorapidity range the minijet energy $`E_t`$ is of the order of 10 GeV. Thus, in this way the nuclear collision geometry (i.e. impact parameter) can well be defined by means of the wide acceptance calorimeter with minimal dependence on collision dynamics details in the central region.

Fig.3 shows correlations between $`E_T`$ and $`b`$ for CaCa, NbNb, PbPb collisions allowing to conclude that an impact parameter can be estimated with 1.5 fm precision.