text:
# Electron-Electron Interactions in Isolated and Realistic Quantum Dots: A Density Functional Theory Study

## I Introduction

The interplay of quantum mechanical interference and electron-electron interaction in semiconductor quantum dots (QDs) has attracted a lot of experimental and theoretical interest in recent years.Kouwenhoven et al. (1997); Alhassid (2000); Aleiner et al. (2002); Reimann and Manninen (2002) Experimentally the ground state properties of a QD can be probed by Coulomb blockade (CB) measurements at near-zero temperature, in which the conductance through the dot is blocked except at particular gate voltages, causing the conductance as a function of gate voltage to be a series of sharp peaks. In particular, the spacing between neighboring peaks is related to the second difference of the ground state energy with respect to electron number $`N`$, $`\mathrm{\Delta }_2E(N)E_{\mathrm{gs}}(N+1)+E_{\mathrm{gs}}(N1)2E_{\mathrm{gs}}(N)`$, which is often called the addition energy. The ground state spin of the dot can be inferred from the motion of the peak position in a magnetic field.Folk et al. (2001); Lindemann et al. (2002) Except in small dots with high geometrical symmetry,Tarucha et al. (1996) whose ground state properties show shell structure following Hund’s rules, most experimental QDs are asymmetric due to either the irregular shape of the confining gate or impurity disorder, and the number of electrons is usually several hundred.Sivan et al. (1996); Patel et al. (1998); Simmel et al. (1997); Lüscher et al. (2001) The properties of QDs show apparently random but reproducible fluctuations as some external parameters such as the chemical potential, the shape of the confining potential, external magnetic field and so on, vary. In this case, only statistical properties are of physical significance.

While the statistical behavior of peak height fluctuations can be well described by a non-interacting random matrix theory (RMT) model,Alhassid (2000) a similar model, the constant interaction (CI) plus RMT, fails dramatically in describing the peak spacing distribution. The importance of e-e interaction was well recognized, but different theoretical approaches differ quantitatively or even qualitatively in many aspects of e-e interactions effects. The numerical studies for small dotsSivan et al. (1996); Prus et al. (1996); Berkovits (1998); Cohen et al. (1999); Walker et al. (1999); Levit and Orgad (1999); Ahn et al. (1999); Bonci and Berkovits (1999); Hirose and Wingreen (2002) yield a Gaussian-like peak spacing distribution without any even/odd pairing effect, which agrees with experimental results on much larger dots. On the other hand, the statistical approaches for large dotsUllmo and Baranger (2001); Usaj and Baranger (2001, 2002) predict a much more pronounced even/odd effect at zero temperature that is absent in experimental results.rem (a)

Besides conductance peak spacings, another quantity that has attracted a lot of interest is the ground state spin in the CB regime. The occurrence of non-minimal ground state spin has been suggested as an explanation for the absence of a bi-modal distribution in conductance peak spacings,Prus et al. (1996) and for the kinks in the parametric motion of the peaks.Baranger et al. (2000) A recent experimental study on the change in CB peak position as a function of magnetic filed by Folk, et al.Folk et al. (2001) found that three out of five even-$`N`$ states had $`S=1`$. As in the case of peak spacing distributions, previous theoretical studies on ground state spin distributions in QDs gave very different results: Simulations on small tight-binding lattice models with random site energies predict a dominant fraction of high spin ground states,Kurland et al. (2000) but results from KS-SDFT studies of small disordered parabolic quantum dotsHirose and Wingreen (2002) showed that the probability of high-spin ground states never exceeds $`50\%`$ even in the strong interaction regime.

Despite intensive efforts to understand the role of e-e interaction in QD systems, there are still many questions that require further investigation. In this paper, we focus on the following issues: (1) The first one is concerned with the evolution of statistical properties of QDs as the interaction strength increases. The e-e interaction strength is usually characterized by the Wigner-Seitz radius $`r_s`$, which is defined in 2D as $`r_s1/\sqrt{\pi n}a_0`$, where $`n`$ is the electron density and $`a_0`$ is Bohr radius.rem (b) Previous theoretical studies gave different answers to this question, indicating the necessity of further theoretical exploration. In the meantime, this issue has attracted the experimental interest, and some relevant experiments are underway. (2) The second issue to be explored is the effect of confining potential. Previous theoretical studies were based on abstract statistical assumptions, chaotic model external potentials, or Gaussian impurities. None of them have considered the realistic confining potential imposed in actual experiments. Though it is usually assumed that the statistical properties of QDs are of universal significance, the validity of that is still to be verified. (3) The third issue that is worthy of further elaboration is how statistical properties of QDs evolve as the size of a QD system, mainly the electron number, increases, which is important not only in solving the puzzling discrepancy between numerical small-dot calculations and semi-analytic large-dot statistical analysis, but also in understanding the fundamental question in mesoscopic physics of how macroscopic properties emerge from microscopic physics via the mesoscopic regime.

In our previous study,Jiang et al. (2003a, b) we reported a Kohn-Sham spin density functional theory (KS-SDFT) study of spin and conductance peak spacing distributions in 2D isolated quantum dots in a classically chaotic quartic oscillator potential with electron number $`N`$ up to 200 at $`r_s1.5`$. We found that for large electron number and asymmetric external potential, the peak spacing distribution is mainly Gaussian-like and there is no observable even/odd effects. We also found enhanced probability of high spin ground state implying stronger interaction effects than expected from RPA approaches. In this paper, we extend the previous study in two directions: (1) Similar calculations are done for weaker ($`r_s0.2`$) interactions to investigate how statistical properties of quantum dots evolve as the interaction strength varies. (2) We extended previous studies in isolated QDs, in which electrons are confined by an analytic chaotic potential without explicit consideration of the coupling with external gate, to QD systems in which the confining potential is calculated by solving the electrostatic Poisson equation with realistic boundary conditions. In the latter case, with more controllable parameters available, we are able to study statistical properties of QDs as a function of both e-e interaction strength ($`r_s`$) and electron number ($`N`$) in a more systematic way.

The outline of the paper is as follows: In the next section, the KS-SDFT method and the numerical techniques used in our study are concisely described. In Section III, we report the statistical properties of isolated quantum dots. In Section IV, the results for realistic quantum dots are presented. Section V concludes the paper with discussion about the significance of these results.

## II Method

### II.1 Kohn-Sham spin density functional theory

In KS-SDFT, the ground state energy of an interacting system with electron number $`N`$ and total spin $`S`$ in the local external potential $`V_{\mathrm{ext}}(𝐫)`$ is written as a functional of spin densities $`n^\sigma `$ with $`\sigma =\alpha ,\beta `$ denoting spin-up and spin-down, respectively,

$`E[n^\alpha ,n^\beta ]=T_s[n^\alpha ,n^\beta ]+{\displaystyle n(𝐫)V_{\mathrm{ext}}(𝐫)𝑑𝐫}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \frac{n(𝐫)n(𝐫^{})}{\left|𝐫𝐫^{}\right|}𝑑𝐫𝑑𝐫^{}}+E_{\mathrm{xc}}[n^\alpha ,n^\beta ]`$ . (1)
(Effective atomic units are used through the paper: for GaAs QDs with the dielectric constant $`\epsilon =12.9`$ and effective electron mass $`m^{}=0.067m_e`$, the values are 10.08 meV for energy and 10.95 nm for length.) $`T_s[n^\alpha ,n^\beta ]`$ is the kinetic energy of the KS non-interacting reference system which has the same ground state spin density as the interacting one, and $`E_{\mathrm{xc}}[n^\alpha ,n^\beta ]`$ is the exchange-correlation energy functional. The spin densities $`n^\sigma `$ satisfy the constraint $`n^\sigma (𝐫)𝑑𝐫=N^\sigma `$ with $`N^\alpha =(N+2S)/2`$ and $`N^\beta =(N2S)/2`$.

Assuming that the ground state of the non-interacting reference system is non-degenerate, the non-interacting kinetic energy is given by $`T_s[n^\alpha ,n^\beta ]=_{i,\sigma }\psi _i^\sigma \left|\frac{1}{2}^2\right|\psi _i^\sigma `$, and the ground state spin density is uniquely expressed as

$$n^\sigma (𝐫)=\underset{i}{\overset{N^\sigma }{}}\left|\psi _i^\sigma (𝐫)\right|^2,\sigma =\alpha ,\beta .$$
(2)
Here $`\psi _i^\sigma `$ are the lowest single-particle orbitals which are obtained from

$$\left\{\frac{1}{2}^2+V_{\mathrm{ext}}(𝐫)+V_H[n;𝐫]+V_{\mathrm{xc}}^\sigma [n^\alpha ,n^\beta ;𝐫]\right\}\psi _i^\sigma (𝐫)=\epsilon _i^\sigma \psi _i^\sigma (𝐫),$$
(3)
where $`V_H[n;𝐫]`$ and $`V_{\mathrm{xc}}^\sigma [n^\alpha ,n^\beta ;𝐫]`$ are the Hartree and exchange-correlation potentials, respectively,

$`V_H[n;𝐫]`$ $``$ $`{\displaystyle \frac{n(𝐫^{})}{|𝐫𝐫^{}|}d^3𝐫^{}},`$ (4)
$`V_{\mathrm{xc}}^\sigma [n^\alpha ,n^\beta ;𝐫]`$ $``$ $`{\displaystyle \frac{\delta E_{\mathrm{xc}}[n^\alpha ,n^\beta ]}{\delta n^\sigma (𝐫)}}.`$ (5)
We have used the local spin-density approximation (LSDA)Parr and Yang (1989); Dreizler and Gross (1990) for the exchange-correlation energy functional,

$$E_{\text{xc}}[n^\alpha ,n^\beta ]n(𝐫)ϵ_{\text{xc}}(n(𝐫),\zeta (𝐫))𝑑𝐫,$$
(6)
where $`\zeta =(n^\alpha n^\beta )/n`$ is the spin polarization, and $`ϵ_{\mathrm{xc}}`$ is the exchange-correlation energy per electron in the homogeneous electron gas. We use Tanatar-Ceperley’s parameterized form of $`ϵ_{\mathrm{xc}}`$ for the spin-compensated case ($`\zeta =0`$) and fully spin polarized case ($`\zeta =1`$).Tanatar and Ceperley (1989) $`ϵ_{\mathrm{xc}}`$ for arbitrary spin polarization is obtained by the following interpolation,Reimann and Manninen (2002)

$$ϵ_{\text{xc}}(n,\zeta )=ϵ_{\text{xc}}(n,0)+f(\zeta )\left[ϵ_{\text{xc}}(n,1)ϵ_{\text{xc}}(n,0)\right],$$
(7)
with

$$f(\zeta )=\frac{(1+\zeta )^{3/2}+(1\zeta )^{3/2}2}{2^{3/2}2},$$
(8)
which is exact for the exchange energy density, but only approximate for the correlation energy density. Recently, Attaccalite, et al.Attaccalite et al. (2002) parameterized a new LSDA exchange-correlation functional form based on more accurate quantum Monte Carlo modeling data that includes spin polarization explicitly. The new functional form was shown to be more accurate than Tanatar-Ceperley’s in both zero and finite magnetic field when compared to quantum Monte Carlo data.Saarikoski et al. (2003) Numerical tests undertaken by us showed, however, that in the regime with which this paper is concerned, the difference between the two forms is minor. In this paper, we will use Tanatar-Ceperley’s form. Comparisons with quantum Monte-Carlo calculations for $`N12`$ and interaction strength parameter $`r_s8.0`$ have shown that LSDA works well for both the ground state spin and energy in 2D parabolic quantum dots,Egger et al. (1999); Pederiva et al. (2000); Umr  as well as for many spin excitation energies.Umr

### II.2 Numerical techniques

To investigate large electron number regime, we developed efficient techniques for the implementation of KS-SDFT method which have been reported elsewhere.Jiang et al. (2003b) Essentially our algorithm includes the following components: (1) We use a particle-in-the-box basis for the representation of KS orbitals and for the action of the kinetic energy operator on the orbital wave functions; (2) The Kohn-Sham total energy is directly minimized by a conjugate gradient method that is based on Teter, Payne and Allan’s methodTeter et al. (1989) but with important modifications, which include a more efficient line minimizaiton scheme and a delayed update of the effective potential;Jiang et al. (2003b) (3) We use a Fourier convolution method to calculate Hartree potential efficiently;Martyna and Tuckerman (1999) (4) A simple one-way multi-grid technique is used to accelerate the convergence.Lee et al. (2000) The high efficiency of our method makes the numerical study of statistical properties of large quantum dots with several hundred electrons become computationally feasible.Jiang et al. (2003a, b)

## III Spin and peak spacing distributions in isolated quantum dots

We first investigate statistical properties of quantum dots in a 2D isolated model system. Since experimental dots usually have irregular shapes that result in chaotic classical dynamics, we use the following quartic oscillator potential (QOP) as the external confining potential,

$$V_{ext}(𝐫)=a\left[\frac{x^4}{b}+by^42\lambda x^2y^2+\gamma (x^2yxy^2)r\right].$$
(9)
Both the classical dynamics and the single-particle quantum mechanics at $`\gamma =0`$ have been extensively studied:Bohigas et al. (1993) the system evolves continuously from integrable to fully chaotic as $`\lambda `$ changes from 0 to 1. The parameter $`\gamma `$ is introduced to break the four-fold symmetry. The prefactor $`a`$ is used to adjust the flatness of the confining potential and hence the electron density.

Fig. 1 shows the contours of a typical $`V_{\mathrm{ext}}`$ used in this study (upper panel), and the corresponding electron density at $`N=200`$ (lower panel). A noticeable feature is that the electron density in isolated QD systems is not uniform: Electrons form peaks near the boundary region because of classical Coulomb repulsion between electrons. In real experimental systems, because of the coupling with the top gate, we expect to see an essentially uniform electron density distribution. Whether such a difference in electron density will result in significant deviation in statistical properties will be discussed in the next section.

For a given $`V_{\mathrm{ext}}`$, the ground state energy $`E_{\mathrm{gs}}`$ and spin $`S_{\mathrm{gs}}`$ as a function of $`N`$ are determined by calculating several spin configurations for each $`N`$ and selecting the one with minimal energy. The addition energy is then calculated as the second difference of $`E_{\mathrm{gs}}(N)`$. As shown in Fig. 2(a), $`\mathrm{\Delta }_2E(N)`$ decreases slowly as $`N`$ increases. This is mainly a classical effect – the increasing capacitance as the dot becomes large. Since we are mainly interested in the quantum mechanical effects, the classical smooth trend, denoted $`\mathrm{\Delta }_2E(N)`$, is fit by a 4-th order polynomial of $`N`$ and then subtracted from $`\mathrm{\Delta }_2E(N)`$. To compare with experiments, the addition energy is further scaled by the mean level spacing, $`\mathrm{\Delta }`$, which is calculated by $`\mathrm{\Delta }=2\pi \mathrm{}^2/m^{}A_{\mathrm{eff}}`$ with $`A_{\mathrm{eff}}`$ estimated as the classical allowed area; the resulting dimensionless spacing is denoted $`\xi `$,

$$\xi (N)\frac{\mathrm{\Delta }_2E(N)\mathrm{\Delta }_2E(N)}{\mathrm{\Delta }}.$$
(10)
To obtain good statistics, we choose five sets of ($`\lambda ,\gamma `$) – (0.53, 0.2), (0.565, 0.2), (0.6,0.1), (0.635, 0.15), and (0.67,0.1) – and for each set of parameters, we calculate a range of electron numbers $`N`$ from 1 to 200.

Fig. 3 plots peak spacing distributions in isolated QDs in different interaction regimes; $`r_s`$ ranges from about $`0.2`$ to $`1.3`$, and $`N=120`$-$`200`$. In the large density or weak interaction regime, $`r_s0.2`$ \[Fig. 3(a)\], the peak spacing distributions for even and odd $`N`$ are sharply different: the distribution for even $`N`$ is asymmetric and close to the Wigner surmise while that for odd $`N`$ is a narrow sharp symmetric peak, which agrees very well with the CI-RMT model. As the interaction strength increases, such an even/odd effect is gradually reduced: the odd $`N`$ distribution becomes increasingly wider, and the even $`N`$ distribution becomes more symmetric. At $`r_s0.9`$ \[Fig. 3(c)\], the only difference between the two is a weak tail in the even $`N`$ distribution at positive spacing. As $`r_s`$ is further increased to about $`1.3`$ \[Fig. 3(d)\], which is the usual experimental regime, the even/odd effect is hardly perceptible, and the shape of the distribution is mainly Gaussian-like.

The insets in Fig. 3 show ground state spin distributions in the corresponding interaction regimes. At $`r_s0.2`$, the minimal spin configuration dominates in both even and odd $`N`$ cases, as one expects according to CI-RMT model. But even in this weak interaction regime, the portion of $`S=1`$ for even $`N`$ is already significant, $`P_{S=1}=0.1`$. The probability of high-spin ground state increases quickly as the e-e interaction becomes stronger. At $`r_s1.3`$, the high spin ground state is more probable than the minimal spin state.Jiang et al. (2003a) Such a result is remarkable considering that $`r_s1.3`$ is only a moderate interaction strength.

## IV Spin and peak spacing distributions in 2D realistic quantum dots

### IV.1 Model

The boundary condition used in isolated quantum dots is very different from that in the experimental systems. Fig. 4 illustrates a typical experimental quantum dot structure, which is fabricated based on a GaAs-AlGaAs heterostructure. Numerical modeling of realistic quantum dotsKumar et al. (1990); Jovanovic and Leburton (1994); Stopa (1996); Nagaraja et al. (1997); Matagne and Leburton (2002) is still a great challenge in terms of computational efforts. The difficulty mainly comes from the fact that although the volume where electrons are confined is small, the boundary conditions that couple the dot with the external confining gate are defined on a much larger scale. Previous numerical studies of realistic quantum dots considered only small dots with less than a few tens of electrons, except in Ref. Stopa, 1996 where by reducing Schrödinger equations from 3D to multi-component 2D the author was able to study realistically a lateral quantum dot with $`N`$ up to 100.

We note that the most important factor in experimental conditions that is absent in the isolated quantum dot model is the coupling between the dot and the top gate. To take that into account, we introduce semi-infinite boundary conditions.Davies et al. (1995); Yakimenko et al. (2001) Using $`𝐑=(𝐫,z)`$ to denote the 3D spatial coordinate, $`𝐑=(𝐫,z)`$ where $`z`$ is the height in the growth direction and $`𝐫`$ the coordinates perpendicular to the growth direction, we take

$`\varphi (𝐫,z=0)`$ $`=`$ $`V_g(𝐫)+V_s,`$
$`{\displaystyle \frac{\varphi }{𝐑}}|_R\mathrm{}`$ $`=`$ $`0,`$ (11)
where $`V_s`$ is the Schottky barrier due to surface states,Davies (1988); Davies et al. (1995) and $`V_g(𝐫)`$ is the gate voltage imposed on the top surface. In 2D cases, the inclusion of semi-infinite realistic boundary conditions is quite straightforward. For $`V_{\mathrm{ext}}`$ in Eq. (3), instead of an analytic form, it is calculated from

$$V_{\mathrm{ext}}^{\mathrm{R2D}}(𝐫;z)=\frac{1}{2\pi }𝑑𝐫^{}V_g(𝐫^{})\frac{|z|}{(|𝐫𝐫^{}|^2+z^2)^{3/2}}+V_{\mathrm{QW}}(z),$$
(12)
where $`V_{\mathrm{QW}}(z)`$ is the confining potential in $`z`$-direction due to the quantum well or heterojunction structure from which the quantum dot is fabricated. In the case of 2D quantum dots, $`z`$ is taken as the distance between the top surface and the mean electron height near the middle of the QD, $`z_0`$. Because of the coupling with the gate, the Hartree potential has an additional image term,Davies et al. (1995); Yakimenko et al. (2001)

$$V_\mathrm{H}(𝐫)=𝑑𝐫^{}n(𝐫^{})\left[\frac{1}{|𝐫𝐫^{}|}\frac{1}{(|𝐫𝐫^{}|^2+4z_0^2)^{1/2}}\right].$$
(13)
This 2D quantum dot model will be denoted R2D to simplify the notation.

Besides the more realistic boundary conditions included in the R2D model, its other advantage is that we can obtain statistics by changing the shape of the confining gate.Chan et al. (1995) In this study, we use a confining gate with the shape illustrated in Fig. 5. By taking $`x_1`$, $`x_2`$, $`y_1`$, and $`y_2`$ randomly, we get a series of irregular confining potentials that are statistically independent as indicated by correlation analysis. Fig. 6 shows the contour of a typical confining potential and a sample electron density for $`N=200`$. We see that indeed the electron density in R2D is much more uniform than that in the isolated case. The availability of a lot of statistically independent irregular potentials enables us to investigate how spin and peak spacing distributions evolve as a function of $`N`$ and $`r_s`$ in a more systematic way.

### IV.2 Comparison of distributions from isolated and realistic quantum dots

Fig. 7 shows spin and peak spacing distributions from isolated and realistic QD systems at two $`r_s`$ values ($`r_s1.0`$ and $`1.5`$) within comparable electron number ranges. Qualitatively, they have similar features: At $`r_s1.0`$, distributions from both QOP and R2D have a weak tail on the positive side for even $`N`$ compared to odd $`N`$, and at $`r_s1.5`$, the tail vanishes. The spin distributions from two types of model systems are also the same within statistical uncertainty. Quantitatively there are indeed some observable differences in both spin and peak spacing distributions, but with the limited statistics available in this study, it is difficult to attribute physical significance to these subtle features. Our study therefore verifies the validity of using model systems to study statistical properties of quantum dots, and provides a direct demonstration of the universal feature of mesoscopic fluctuations.

### IV.3 Spin distributions as a function $`r_s`$

Fig. 8 shows the spin distribution as a function of the interaction strength in the electron number range $`N=11`$-$`40`$. This study extends previous investigations for weak and moderate interaction strength in isolated dots to the stronger interaction regime. In agreement with previous studies, the fraction of high-spin ground states increases as the electron-electron interaction becomes stronger. There are several features that are worth comment. (i) Crossing: For even $`N`$, the probability of $`S=0`$ crosses with that of $`S=1`$ at $`r_s3.0`$. (ii) Saturation: Probabilities of ground state spins for even $`N`$ saturate at about $`r_s3.5`$. The spin distribution for odd $`N`$ is not saturated in the interaction strength range studied here, but the trend indicates that saturation is likely. (iii) Higher spin states: The probabilities of $`S=2`$ and $`S=5/2`$ remain small for all interaction strengths, though $`P(S=2)`$ becomes significant at $`r_s>3.0`$.

### IV.4 Size Effects

To investigate how statistical properties of quantum dots evolve as a function of electron number $`N`$, we calculated peak spacing and spin distributions in different electron number ranges for $`N`$ up to 200 at $`r_s1.0`$

Fig. 9 shows peak spacing distributions in different $`N`$ ranges. For $`N=11`$-$`40`$, the distributions for even and odd $`N`$ are obviously different. The difference decreases as $`N`$ increases, and becomes barely observable for $`N=131`$-$`160`$ and $`N=171`$-$`200`$. Our results show clearly that the conductance peak spacing distribution depends on electron number, and is significantly different between small ($`N<50`$) and large ($`N>100`$) dots.

Fig. 10 shows the trend in ground state spin as $`N`$ increases. For this quantity, the dependence on $`N`$ is even more significant: the probability of high spin ground states increases as $`N`$ increases. In particular, $`P(S=1)`$ increases from 0.27 at $`N=11`$-$`40`$ to 0.45 at $`N=171`$-$`200`$. Another notable feature is that the spin distribution saturates at large $`N`$.

### IV.5 Two-level model analysis

KS-SDFT calculations reported above reveal some interesting features in conductance peak spacing and ground state spin distribution. To obtain a clearer picture of the underlying physics, and, in particular, to understand the main factors that determine the spin distribution in different regimes, we use a simple two-level model (TLM) to investigate the probability of singlet and triplet ground state spin for even $`N`$ cases. The same model was used as an interpretive tool in Hirose and Wingreen’sHirose and Wingreen (2002) work on spin distributions in small disordered parabolic quantum dots, and in doing that, they attribute saturation of the spin distribution to the effects of off-diagonal interaction matrix elements. Our analysis using single-particle wave functions in truly ballistic chaotic systems, however, leads to very different conclusions.

The basic idea of the two-level model is to consider a two-electron system with only two single-particle levels. The interacting matrix elements (IME)

$$V_{ijkl}𝑑𝐫𝑑𝐫^{}\psi _i^{}(𝐫)\psi _j^{}(𝐫^{})v_{\mathrm{sc}}(𝐫,𝐫^{})\psi _k(𝐫^{})\psi _l(𝐫).$$
(14)
are calculated using a screened interaction kernel $`v_{\mathrm{sc}}(𝐫,𝐫^{})`$ to implicitly account for the effects of all other electrons. In particular, we use the screened interaction obtained from the 2D random-phase approximation (RPA),

$$v_{\mathrm{sc}}(q)=\frac{v_\mathrm{c}(q)}{1v_\mathrm{c}(q)\chi (q)}$$
(15)
where $`v_\mathrm{c}(q)=2\pi /q`$ is the Fourier transform of the bare Coulomb interaction, and $`\chi (q)`$ is the susceptibility of the 2D electron gas. It is straightforward to include the screening due to the top gate by using $`v_\mathrm{c}(q)=(2\pi /q)(1\mathrm{e}^{2z_0q})`$, but in the regime considered in this paper, the effects of the gate screening on the interacting matrix elements are negligible.

In this TLM system, one can obtain singlet and triplet state energies either by exact diagonalization (ED) or in the framework of Hartree-Fock (HF) theory (ie., first-order perturbation theory). In TLM-HF, the difference of energy between singlet and triplet states has the following simple form,

$$\delta EE_{S=1}E_{S=0}=\delta \epsilon [J_{11}(J_{12}K_{12})],$$
(16)
where $`\delta \epsilon `$ is the single particle level spacing, $`J_{ij}V_{ijji}`$ are the direct Coulomb interaction terms, and $`K_{12}V_{1212}`$ is the exchange term.

The main difference between our TLM analysis and Hirose and Wingreen’s is in the way to choose single-particle energies and wave functions. In our analysis, we use single-particle eigen-energies and wavefunctions in the extended Thomas-Fermi effective potentialBrack and Bhaduri (1997) that is obtained within the same realistic confining gate structure as in KS-SDFT calculations. Wave functions obtained in this way are guaranteed to be those of ballistic systems, and possible complications related with the disorder strength (ie., localization) can therefore be avoided. We note that the two-level model with this way of calculating interacting matrix elements can be regarded as a simplified version of the model proposed by Ullmo, et al. in Ref. Ullmo and Baranger, 2001.

Fig. 11(a) shows the probability of finding triplet ground states, $`P(S=1)`$, as a function of $`r_s`$ obtained from TLM-ED, TLM-HF, as well as KS-SDFT calculations. We see a striking agreement between TLM-ED and KS-SDFT results, which is indeed remarkable considering the extreme simplicity of our two-level model compared to KS-SDFT. Both TLM-HF and ED reproduce correctly the saturation of spin distribution at $`r_s3.5`$. To explore possible mechanisms for the saturation, Fig. 11(b) plots the mean values of $`J_{11}`$, $`J_{12}`$, $`K_{12}`$, and the off-diagonal term $`V_{1121}`$, (the other off-diagonal term, $`V_{2221}`$, is close to $`V_{1121}`$, and is therefore not plotted here) as a function of $`r_s`$. It is clearly shown that off-diagonal terms are an order of magnitude smaller than the Coulomb and exchange terms, indicating that the off-diagonal terms are negligible in terms of their contribution to the spin distribution. This is verified by the general agreement between TLM-HF and TLM-ED. On the other hand, Fig. 11(b) shows that though $`J_{12}`$ and $`K_{12}`$ do saturate at high $`r_s`$, $`J_{11}`$ does not. This indicates that considering mean values of $`J_{11}`$, $`J_{12}`$, and $`K_{12}`$ alone can not explain the saturation in the spin distribution.

We also use the TLM to interpret the effect of finite size on the spin distributions, and the results are plotted in Fig. 12. The simple two-level model also grasps most of the physics: the probability of the triplet ground state calculated from both TLM-ED and HF increases as $`N`$ increases, but saturates at large $`N`$, in agreement with the KS results. We note that the difference between different $`N`$ ranges is mainly due to wave functions rather than eigenvalues (which are always random-matrix like for our system); for a $`N`$-electron system, the wave functions involved in the TLM are $`\psi _{N/2}`$ and $`\psi _{N/2+1}`$. Our TLM results show clearly that the $`N`$-dependence of the spin distributions comes mainly from the changing statistics of the wave functions. This is more directly demonstrated in Fig. 12(b), where mean values of $`J_{11}`$, $`J_{12}`$, $`K_{12}`$ and $`V_{1121}`$ are plotted. We note that the $`N`$-dependence of wave-function statistics is usually neglected in most statistical theories of quantum dots, where the large $`N`$ limit is usually taken.Alhassid (2000); Ullmo and Baranger (2001); Aleiner et al. (2002); Usaj and Baranger (2002) Since many experimental data were obtained with $`N`$ less than 200, our studies show clearly the necessity of considering finite $`N`$ effects in interpreting those data.

## V Summary

In this paper, we use Kohn-Sham spin density functional theory to study electron-electron interaction effects in quantum dots. In particular, we calculated conductance peak spacing and ground state spin distributions in an isolated quantum dot with classically chaotic quartic external potential and in a realistic quantum dot with explicit consideration of the coupling between the dot and the external confining gate. In the former case, as an extension of our previous study, we investigated how the statistical properties evolve when the e-e interaction strength increases from the weak ($`r_s0.2`$) to moderate regimes ($`r_s1.5`$) for $`N=120`$-$`200`$. In the case of the realistic quantum dots, we first made a direct comparison between isolated and realistic quantum dots, and found that in spite of the difference in the boundary conditions, their statistical properties are qualitatively similar. We further investigated spin and peak spacing distributions in strong interaction regimes with $`r_s`$ up to 4.0 and their dependence on the electron number $`N`$. The even/odd effect in peak spacing distributions vanishes as $`r_s`$ increases and becomes weaker as $`N`$ increases.

We close with two comments on our results. First, our study, as well as other theoretical investigations, have made it quite clear that the lack of even/odd pairing in experimental peak spacing distributions is due to the e-e interaction, and the CI-RMT model is therefore not valid for the description of mesoscopic fluctuations in Coulomb blockade peak spacings. In theoretical condensed matter physics, it is well accepted that the e-e interaction at $`r_s1.5`$ can be accurately described by many-body perturbation theory, most notably the RPA model. However, the results obtained from the RPA model combined with a RMT description of single particle properties disagree significantly from those obtained with numerical approaches. Understanding the sources for such discrepancies will deepen our general physical picture of quantum dots.

Second, our results were obtained under the local spin density approximation within KS-SDFT, whose validity, though well established in 3D bulk material systems through decades of experience, is not necessarily well justified for 2D large quantum dots with irregular confining potential and/or in the strong interaction regime. In the small dots for $`N`$ up to about 30, it is still computationally feasible to check the validity of LSDA by comparing with more accurate quantum Monte Carlo modeling, but not for the large-dot regime ($`N>50`$). On the other hand, the statistical quantities investigated above are experimentally accessible. Comparing future experimental data with LSDA predictions will provide a unique chance to check LSDA in real systems.

###### Acknowledgements.

We appreciate discussions with G. Usaj and C. J. Umrigar. This work was supported in part by NSF Grant No. DMR-0103003.


format: text

text:
# The Opinion Game: Stock price evolution from microscopic market modelling

## 1. Introduction

The financial markets constitute an intriguing and complex system that has not failed to attract mathematicians and scientists from other fields for a long time. Only rather recently, however, has mathematical finance, and more specifically the theory of derivatives on the stock markets become a major field of mathematics and one of the major sources of inspiration for probability theory in general and stochastic analysis in particular. The reason for this development is simple: It is based on the apparent success of the so-called Black-Scholes formula for the fair price of an option as a tool for the actual trader on the market. Indeed, the very existence of this mathematical theory appears to be largely responsible for the recent growth and diversification of the derivative market itself, which in its present form would have been impossible without an underlying mathematical theory. On the other side, the great success of this same theory in the mathematical sciences is due, beyond the obvious advantages it provides for the careers of students trained in this field, largely to the fact that there is a very clear mathematical setting for this theory with clearly spelled out axioms and assumptions which allows for the mathematician to bring his traditional weapons to bear in a familiar terrain.

One of the crucial issues in the financial mathematics is the modeling of prices of commodities (stocks, currencies etc.) with help of stochastic processes. The main approaches that have been used in this context are the following:

* generalized Black-Scholes (BS) theory. Originally the BS theory  (for good textbook exposition see, e.g., ) emerged as a theoretical foundation of pricing of derivatives (options) of underlying financial instruments (stocks, currencies etc.). Initially the price of the underlying was taken as geometric Brownian motion. Later this theory was generalized and put on the axiomatic background based on the assumption of the so-called non-arbitrage condition. This led to the conclusion that prices are described by semi-martingale measures, and thus are essentially given by solutions of certain stochastic differential equations. This framework has been the main driving force of the rapid growth of financial mathematics in the last decades. Note that the main purpose of this theory is not to derive specific models for the underlying but rather to deduce consequences from generally accepted principles. In this respect this approach can be compared to classical thermodynamics.
* statistical approach. Empirical studies of share price data try to model the data by certain stochastic process. Popular classes of models are ARCH, GARCH, ARMA etc. There seem to be no totally conclusive results, but certain interesting phenomena have been observed, such as universal exponents in certain correlation data (see ). Quite frequently, analogies to phenomena like turbulence are drawn. This approach does not usually intend to derive the model from any underlying economic theory.
* agent based models. “As prices are generated by the demand of agents who are active on the financial market for the given asset, a \[…\] model \[…\] should be explained in terms of the interaction of these agents.” Based on this observation stated by Föllmer in , a large number of agent based models for price evolution have been developed. Föllmer in  suggested a model of diffusion in random environment by viewing the price process as a sequence of temporary equilibria in a market with agents. Among the most popular microscopic models for financial markets are so-called agent based models that are more or less sophisticated versions of the “minority game” (MG); for reviews, see Jeffries and Johnson , Bouchaud and Giardina , and references therein. In all these models there exists a collection of traders, each endowed with the possibility to make a decision (typically of the type “buy”, “sell”, or “hold”) concerning a given investment. To reach such a decision, the trader disposes of a certain number of “strategies”, a strategy typically being a function of the price history into the set of decisions. The game consists of the traders choosing their strategies in a way as to reach certain objectives (in the minority game to be “in the minority”, in reality to make a maximal profit from the transactions undertaken). At each time step, the asset price is updated according to an empirical rule as a function of the number of “buyers” and “sellers”. There are currently a large number of versions of such models around, including models with additional stochastic components. These models exhibit a rather rich dynamical structure. However, they are rather heavy handed both analytically (where little or nothing is known on a mathematical level) and numerically. Moreover, the large number of assumptions and parameters entering the models makes their predictive power somewhat limited. Purely deterministic models of this type have been criticized before (Bouchaud and Giardina ) and stochastic models have been proposed that should allow to take into account irrational behaviour of agents.

In this paper we want to propose agent based market models that are much simpler and that, at the same time, are build reasonably close to what is actually happening in financial markets. The main distinction between our model and the MG type models is that it focuses on collective effects of a market while not attempting to model the actual reasoning process of an individual trader. The basic paradigm of our modeling approach is the notion of a price. Prices of a share of stock or other commodity arise from trading. There are various developed theories in economics concerning prices based on some equilibrium assumptions, but, fortunately, in the stock market in particular, the price of a share is obtained by a well defined procedure which is easily implementable in an algorithm provided a sufficient amount of information about the state of the market at any given time is available. The basic principle of our approach is the modeling of a set of interacting agents in a way that allows to extract the price from the current trading state in a rigorous way.

The remainder of this article is organized as follows. In Section 2 we explain the basic principles of our modeling approach in detail and illustrate them in some special cases. In Section 3 we analyze some of the phenomena that can already be observed in this simplest setup and relate them to certain features of real world data. In Section 4 we introduce a number of additional features that should be implemented to obtain more realistic models, and discuss some of the additional effects they produce. In Section 5 we present our conclusions.

## 2. Our approach

The basic procedure of price determination in real stock markets is done by trading. That is to say, participating traders propose prices at which they are willing to sell respectively buy a certain number of shares of a given stock. The list of these offers at a stock exchange is called the orderbook. On the basis of this orderbook the market maker is matching buy and sell offers according to certain rules until no further transactions are possible. That is to say, after the matching procedure (clearing) the highest price proposed for buying a share is smaller than the lowest price at which a share is offered. These two values are then quoted as the bid price and the ask price. The dynamics of the real stock market thus has two components:

* the changing of buy and sell offers made by traders and
* the matching of these offers by the market maker which fixes the quoted price at a given time.

Our purpose is to develop a class of models which reflects this mechanism of pricing and allows for diverse modeling on different levels of complexity of the behaviour of the agents, while maintaining the pricing mechanism by the market maker. To do so, a minimal requirement for the description of the state space of the trading agents is that it must allow us to recover the state of the orderbook at any given moment in time.

Our idea is thus to consider the time evolution of a virtual orderbook or a “trading state” containing the opinions of each participating agent about the “value”111We will distinguish the notion of the value from that of the price. The value is what agents have an opinion about, while the price is determined by the market. The opinion on the value can be driven by fundamental considerations (e.g. earning or dividend expectations, typically coming from outside information), or speculative considerations (e.g. predictions based on partial knowledge on the current state of the opinions of other traders), or both. of the stock. The evolution is driven by the change in opinion of the agents and the action of the market maker.

A minimal model in which this idea can be implemented can be described as follows. We consider trading in one particular stock. Assume that there are $`N`$ “traders” and $`M<N`$ shares of the stock. We make the simplifying assumption that each trader can own at most one share. The state of each  trader $`i`$ is given by its opinion $`p_i`$ about the logarithm of the value of the stock, and by the number of shares he owns $`n_i\{0,1\}`$. This is to say, the trader $`i`$ would be willing to sell his share at the price $`e^{p_i}`$, if he owns one ($`n_i=1`$), respectively buy a share at this price, if he does not own one ($`n_i=0`$). We say that a trading state is *stable*, if the $`M`$ traders having the $`M`$ highest opinions $`p_i`$ all own a share. This means in particular that in a stable state one can infer the set of owners of shares from the knowledge of the state of opinions $`𝐩=(p_1,\mathrm{},p_N)`$. Thus a stable trading state is completely determined by the set of $`N`$ values $`p_i`$, and we will in the sequel identify stable trading states with the vector $`𝐩`$. As we will normally only work with stable trading state, we suppress the qualifier stable when no confusion can arise.

Given a stable trading state $`𝐩`$ we denote by $`\widehat{𝐩}=(\widehat{p}_1,\mathrm{},\widehat{p}_N)`$ its order statistics, that is $`\widehat{p}_i=p_{\pi _i}`$ for a permutation $`\pi \pi (𝐩)`$ of the set of $`N`$ elements such that $`\widehat{p}_1\mathrm{}\widehat{p}_N`$. Then the number of shares owned by traders, $`𝐧(𝐩)=(n_1(𝐩),\mathrm{},n_N(𝐩))`$ satisfies

$$n_i(𝐩)=\{\begin{array}{cc}0,\hfill & \text{ if }\pi _i(𝐩)NM\text{,}\hfill \\ 1,\hfill & \text{ if }\pi _i(𝐩)NM+1\text{.}\hfill \end{array}$$

To a trading state $`𝐩`$ we associate the ask price

$$p^ap^a(𝐩)=\widehat{p}_{NM+1},$$

and the bid price

$$p^bp^b(𝐩)=\widehat{p}_{NM}.$$

Obviously, $`p^a(𝐩)`$ is the lowest price asked by traders owning a share and $`p^b(𝐩)`$ is the highest price offered by traders wanting to buy a share. For convenience we will refer to $`\frac{1}{2}(p^a+p^b)`$ as the current price in the sequel.

Any dynamics $`𝐩(t)`$ defined on the trading state induces the dynamics of $`\widehat{𝐩}`$ and in particular of the pair $`(p^a,p^b)`$.

Our next simplifying assumption is that the above trading state $`𝐩`$ evolves in time as a (usually time-inhomogeneous) Markov chain222We will discuss the Markovian assumption later. $`𝐩(t)`$ with state space $`^N`$. We will further assume that time is discrete (this is inessential but more convenient for computer simulations) and that the updating proceeds asynchronously, i.e. typically at a given instant only a single opinion changes. This dynamics can be considered as an interacting particle system, however, some special features should be incorporated that reflect the peculiarities of a market.

The first and most obvious one is that the ask and bid prices $`p^a(t),p^b(t)`$ are likely to be important quantities which will influence the updating probabilities.

Moreover, it is reasonable to distinguish between transitions that leave $`𝐧`$ unchanged, and those for which $`𝐧\left(𝐩(t+1)\right)𝐧\left(𝐩(t)\right)`$. In the latter case we say that a transaction has occurred.

Before we discuss some more specific implementations of this general setup a few remarks concerning some features that may appear offending are in place.

The first is doubtlessly the assumption that the process is Markovian. This appear unnatural because the most commonly available information on a stock is the history of its price, the “chart”, and most serious traders will take this information at least partly into account when evaluating a stock, with some making it the main basis for any decision.333A frequently heard remark being that all information on a stock is in its price. Certainly one could retain such information and formulate a non-Markovian model, as e.g. the model of Bouchaud and Giardina . However, if one starts to think about this, one soon finds that it is very difficult to formulate reasonable transition rules on the basis of the price history. On a more fundamental level, one will also come to the conclusion that the analysis of the history of the share price is in fact performed in order to obtain information of the current opinion of the traders concerning the value of the stock with the hope of inferring information on the future development of the price. For instance, if one knew that there are many people willing to sell shares at a price not much higher than $`p^a`$, one knows that it will be difficult for the price to break through this level (this is know as a “resistance” by chart-analysts and usually inferred from past failures to break through such a level). Therefore, instead of devising rules based on past price history, we may simply assume that the market participants have some access to the prevailing current opinions, obtained through various sources (chart analysis, rumours, newspaper articles, etc.) and take this into account when changing their own estimates.

The second irritating point is that money does not appear in our model except in the form of opinions about values. In particular, we do not keep track of the cashflow of a given investor (that is to say we do not care whether a given investor wins or looses money). There are various reasons to justify this. First we consider that the market participants do not invest a substantial fraction of their assets in this one stock, so that shortage of cash will not prevent anyone to buy if she deems opportune to do so (in the worst case money can be obtained through credits). Then, money is not conserved, but the total value of the stock can inflate as long as there is enough confidence. Also, we do not keep track of the objective success of a trader, because we do not know how this will eventually influence her decisions. While a given trader may follow her personal strategy with the hope of making profits, we cannot be sure that these strategies will succeed. What is important and what is built into our model, however, is the fact that any trader444We could enlarge the model to incorporate a small fraction of traders which do not act according to common sense, or against their own convictions (e.g. traders that have bought their stock on credit and are executed by their creditors on falling prices. It may be interesting to consider the effect of that in the context of market crashes). will have the subjective impression to make a profit at any transaction.555Which also implies that this subjective opinion must be wrong at least for one of the traders involved. But this seems to reflect reality. Thus we feel that opinions about values are the correct variables to describe such a market rather than the actual flow of capital, at least at the level of a simple model.

The above setting suggests a rather general and flexible class of models of a stock market. Its main feature is that it describes the time evolution of a share price as the result of an interacting random process that reflects the change of the opinions of individual traders concerning the value of the stock. Even when this last process is modelled as a Markov process, the resulting price process $`(p^a(t),p^b(t))`$ will in general not be a Markov process.

## 3. Examples

### 3.1. Ideal gas approximation

Obviously the simplest model for the dynamics of the trading state $`𝐩(t)`$ is to choose it as a collection of independent identically distributed one-dimensional Markov processes (“random walks”) $`p_i(t)`$. This corresponds to the ideal gas approximation in statistical mechanics. In this case, the price process is simply obtained from the order statistics of independent processes and asymptotic results for $`M`$ and/or $`N`$ large (recall that $`M`$ denotes the number of traded shares and $`N`$ the total number of traders) can be obtained rather easily. While this model is somewhat simplistic, some rather interesting phenomena can already be modelled in this context, as we will explain now.

We may be interested in a situation where some macro-economic model may predict several stable (respectively metastable) values of the stock price, realized as the minimum of some utility function $`V`$.

In such a situation it seems not unreasonable to model the process of a single trader as a one-dimensional diffusion process with drift obtained from a potential function $`V`$, i.e. we can take $`p_i(t)`$ to be a solution of the stochastic differential equation

$$dp_i(t)=V^{}\left(p_i(t)\right)dt+\sqrt{\epsilon }dW_i(t)$$

with $`W_i(t)`$ i.i.d. standard Brownian motions, and $`\epsilon >0`$ a parameter measuring the diffusivity. Alternatively, we can take discrete approximations of this process, as will always be done in numerical simulations.

Let us consider the situation when there are two (meta) stable values of the price, $`q_1`$ and $`q_2`$, i.e. the situation where the potential $`V`$ has two minima (wells) at $`q_1`$ and $`q_2`$. If the potential is strong, resp. $`\epsilon `$ is small, an individual trader would typically spent long periods of time near one of the favoured values $`q_1`$ or $`q_2`$.

Let $`w_1`$ be the escape rate from the well $`q_1`$ to the well $`q_2`$ and let $`w_2`$ be the escape rate in the opposite direction. It is well known that these escape rates are exponentially small, $`w_i\mathrm{exp}\left\{2(V^{}V_{q_i})/\epsilon \right\}`$, where $`V^{}`$ is the value of the maximum of $`V`$ on the interval $`[q_1,q_2]`$, if $`\epsilon `$ is small. Denote by $`A_t`$ the number of traders in the right well $`q_2`$ at time $`t0`$. Suppose that $`A_0`$ is much larger than $`M`$ implying that the actual price at the initial moment is situated near $`q_2`$. We are interested in describing the moment of the “crash”, i.e. when the price moves from the right well $`q_2`$ to the left well $`q_1`$.

Then we can approximate the individual processes $`p_i(t)`$ by a two-state Markov chains with state space $`\{q_1,q_2\}`$ and transition rates $`w_1`$ and $`w_2`$. In this approximation we can compute the normalized expected number $`a_t𝐄A_t/N`$ of traders in state $`q_2`$ at time $`t`$ as the solution to the ordinary differential equation

$$\frac{d}{dt}a_t=w_2a_t+w_1(1a_t).$$

We get

$$a_t=\frac{w_1}{w_2+w_1}+\left(a_0\frac{w_1}{w_2+w_1}\right)\mathrm{exp}\left\{(w_2+w_1)t\right\}$$

and the crash time $`T_c`$ can be defined as $`t`$ such that $`a_t=M/N`$,

$$T_c=\frac{1}{w_2+w_1}\mathrm{log}\frac{a_0(w_2+w_1)w_1}{\frac{M}{N}(w_2+w_1)w_1}$$

If the energy barrier $`\mathrm{\Delta }VV^{}V_{q_2}`$ is large enough, the time for each single buyer to escape from the initial well is much larger than the relaxation time for the system of $`A_t`$ particles in the right well, and thus it is natural to expect that the system will pass through the sequence of local equilibrium states corresponding to $`A_t`$ independent random walkers. Using this observation, the evolution of the price can be described in terms of the $`B_tA_tM`$’s order statistic in a system of $`A_t`$ random variables whose distribution is approximately Gaussian with parameters $`q_2`$ and $`\left(V^{\prime \prime }(q_2)\right)^1`$.

To do this, let $`F()`$ denote the distribution function of an individual walker conditioned to stay near $`q_2`$ and define $`q_B=q(B_t)`$ as the solution to the equation

$$\frac{M}{A_t}=1F(q_B)1\mathrm{\Phi }\left((q_Bq_2)\sqrt{V^{\prime \prime }(q_2)}\right),$$

where $`\mathrm{\Phi }()`$ is the distribution function of the standard Gaussian variable. Using the well-known asymptotics for the tail distribution of $`\mathrm{\Phi }()`$,

$$\frac{u}{u^2+1}\frac{e^{u^2/2}}{\sqrt{2\pi }}\frac{1}{\sqrt{2\pi }}_u^{\mathrm{}}e^{x^2/2}𝑑x1\mathrm{\Phi }(u)\frac{1}{u}\frac{e^{u^2/2}}{\sqrt{2\pi }},u\mathrm{},$$

we immediately get

$$q_Bq_2\left(\frac{2\mathrm{log}(A_t/B_t)}{V^{\prime \prime }(q_2)}\right)^{1/2}+\frac{\mathrm{log}(4\pi \mathrm{log}(A_t/B_t))}{\sqrt{8V^{\prime \prime }(q_2)\mathrm{log}(A_t/B_t)}},B_t/M0.$$

Consequently, in the limit of large $`A_t`$, $`M`$ such that $`\rho =1M/A_t`$ is fixed we have

$$q_Bq_2\left(\frac{2\mathrm{log}(\rho ^1)}{V^{\prime \prime }(q_2)}\right)^{1/2}.$$

In this regime of “increasing ranks”, the fluctuations of the $`B_t`$’s order statistic have Gaussian behaviour and their scaling can be derived from \[11, Theorem 2.5.2\]. To do this, consider a small enough $`y`$ such that

$$A_tF(q_B+y)\left(1F(q_B+y)\right)$$

is large for large $`B_t`$ and such that

$$\frac{MA_t\left(1F(q_B+y)\right)}{(B_tM/A_t)^{1/2}}\tau $$

as $`A_t`$, $`M=(1\rho )A_t`$, and $`B_t=\rho A_t`$ are getting large. In view of the definition of $`q_B`$ the LHS expression above equals

$$\sqrt{\frac{B_tA_t}{M}}\left(\frac{F(q_B+y)}{F(q_B)}1\right),$$

where the last ratio can be approximated by

$$\frac{q_2q_B}{q_2q_By}\frac{\mathrm{exp}\left\{(q_B+yq_2)^2V^{\prime \prime }(q_2)/2\right\}}{\mathrm{exp}\left\{(q_Bq_2)^2V^{\prime \prime }(q_2)/2\right\}}1+(q_2q_B)V^{\prime \prime }(q_2)y,$$

assuming that $`(q_2q_B)^2V^{\prime \prime }(q_2)`$ is large enough. As a result, in the limit of large $`M`$ and $`B_t`$, we have

$$\tau \frac{\sqrt{\rho M}(q_2q_B)V^{\prime \prime }(q_2)}{1\rho }y\frac{\sqrt{2\rho MV^{\prime \prime }(q_2)\mathrm{log}(\rho ^1)}}{1\rho }y.$$

It remains to observe that Theorem 2.5.2 from  implies then that the ask price $`p^a(t)`$ satisfies

$$\mathrm{𝐏𝐫}\left(p^a(t)q_B+y\right)\mathrm{\Phi }(\tau ).$$

In other words, the price corresponding to a system of $`A_t`$ such agents has mean $`q_B`$, shifted away from the well $`q_2`$ on a distance of order $`(\mathrm{log}(\rho ^1)/V^{\prime \prime }(q_2))^{1/2}`$ and the variance of the price (ie, the volatility) diverges as $`(MV^{\prime \prime }(q_2))^1(\rho \mathrm{log}(\rho ^1))^1`$ in the limit of small $`\rho =B_t/(M+B_t)`$.

Finally, recalling the hydrodynamic description of $`B_t`$ and the definition of $`T_c`$, we can also describe the time dependence of the mean $`𝐄\rho `$,

$$𝐄\rho =𝐄\rho _t\frac{1\mathrm{exp}\{(w_2+w_1)(T_ct)\}}{1+w_1\mathrm{exp}\{(w_2+w_1)t\}/\left((w_2+w_1)a_0w_1\right)},tT_c.$$

### 3.2. Interacting traders

The simple model introduced in the previous section is rather artificial and simplistic. In reality one would expect that the behaviour of a trader is influenced by the information received from the market as well as external influences. Moreover, the opinion held by a trader with respect to the current price should somehow reflect some his intrinsic psychological characteristics. Finally, the exchange of shares occurring when a transaction takes place should have some visible effect on the time evolution. In the following we suggest some minimal features that should be incorporated in an interacting model to take into account some such effects. We will see that these features correspond to types of interactions that are not commonly considered in the theory of interacting particle systems.

* The derived process of the (ask and bid) price is the most easily accessible piece of information about the trading state of the market for any trader. It is natural that the updating rules should take the current value of this process into account. The simplest and natural modification is to introduce a bias towards the actual price $`(p^a(t),p^b(t))`$ into the distribution of opinion change.
* Traders whose opinion is far from the current price are likely not to pay much attention to what is happening on the market. It is reasonable to assume that they update their opinion less frequently. This feature can be included by reducing the overall transition rates as a function of $`p_i(t)p(t)`$.
* Finally, it is natural to assume that the traders performing a transaction, that is exchange of a share, will update their opinions according to some special rules reflecting the fact that someone buying or selling a share at a given price believes that she has struck a favourable deal, i.e. they attribute a higher value to the share then what they paid, respectively a lower one then what they got.

In the following we describe some concrete framework in which these features are implemented. We describe the construction of the process algorithmically.

Change of opinion: At any time step we first select at random a trader. We will allow this probability to depend on $`p_i`$, and we choose trader $`i`$ with a probability proportional to $`f(p_i(t)p(t))`$, where we define the “current price” via $`p(t)=\left(p^a(t)+p^b(t)\right)/2`$ and the function $`f(x)0`$ has its maximum at zero. The function $`f`$ is responsible for the slow-down phenomenon away from $`p(t)`$. Once a trader has been selected, she changes her opinion from $`p_i`$ to $`p_i^{}`$ with probability proportional to $`q(p_i,p_i^{})`$ which in turn may depend on the entire state of the system. A possible choice for these functions is

$$f(x)=1/(1+|x|)^\alpha ,\alpha >0,$$

and $`q(x,y;𝐩)`$ being, for any fixed $`𝐩`$, a kernel of a random walk. In typical cases, $`q(x,y;𝐩)`$ depends on $`𝐩`$ through $`p^a`$ and $`p^b`$ only. Once $`p_i^{}`$ is chosen, we check whether $`p_i^{}<p^a`$, if $`n_i=0`$, resp. whether $`p_i^{}>p^b`$, if $`n_i=+1`$. If this is the case, we set $`p_i(t+1)=p_i^{}`$, and $`p_j(t+1)=p_j(t)`$ for all $`ji`$, and continue to the next time step. Otherwise, we perform

Transaction: Assume first that $`n_i(t)=0`$ and $`p_i^{}p^a(t)`$. This means that the buyer $`i`$ has decided to buy at the current asked price. Since by definition there is at least one seller who asks only the price $`p^a(t)`$, we select from all these one at random with equal probabilities. Call this trader $`j`$. Then we set

$$p_i(t+1)=p^a(t)+g,p_j(t+1)=p^a(t)g$$

where $`g>0`$ is a fixed or possibly random number. Similarly, if $`n_i(t)=1`$ and $`p_i^{}p^b(t)`$, the seller $`i`$ sells to one of the buyers that offer the price $`p^0(t)`$, and we set

$$p_i(t+1)=p^o(t)g,p_j(t+1)=p^o(t)+g$$

The final state in all cases represents a new stable trading state and the process continues. Note that $`g`$ should be at least as large as to cover the transaction cost.

These additional features make the mathematical analysis of the model much more difficult, but they introduce some interesting effects that are somewhat similar to phenomena observed in real markets. Let us briefly comment on these.

In the ideal gas approximation in the absence of any confining potential all opinions would in the long run spread over all real numbers and the individual opinions would get arbitrarily far from each other. This is avoided by the mechanism of the attraction to the current price.

Slowing down the jump rate of the particles far away from the current price naturally introduces long time memory effects that lead to special features in the distribution of the price process. One of them is possible existence of resistances: if there is a large population of traders in a vicinity of some price $`p`$ far from (say above) the current price $`p(t)`$, then this population tends to persist for a long time, unless the current price approaches this value. If that happens, e.g. due to the presence of an upward drift, one observes a slow-down of the upward movement of the price when it approaches this value. The market has a resistance against increase of the price through this value, reflected in multiple returns to essentially the same extremal values for a long period, see Fig 2.

A further effect of the slowing down is the tendency of the creation of “bubbles” in the presence of strong drifts. In this case one observes a fast motion of the price accompanied by a depletion of the population below this price. Effectively, a few (buying) traders move with the drift, while most are left behind. Such a situation can lead to a crash, if at some moment the drift is removed (due to external effects). Such a scenario was played out in a simulation that is shown in Figure 3.

Note that after the crash there was a strong increase of volatility.

The effect of pushing the opinions from the current price after the transaction  depletes the vicinity of the price and therefore increases the volatility. This effect goes in the opposite direction as the attraction to the price and the interplay of both effects can lead to a non-trivial quasi-equilibrium state. The effect of these mechanisms on the price fluctuation will be studied in forthcoming paper .

To illustrate the influence of different effects discussed above, we present in Fig. 4 and 5 simulation results of the crash-type scenario from Sect. 3.1 where three different parameters of the model were changed independently. All simulations are based on a discretized version of the diffusion model from Sect. 3.1 with the same potential function $`V`$ having two local minima (one metastable and one stable). All simulations start with the same initial condition where all traders are located near the metastable minimum.

Figure 4 shows simulations without the feature of slowing down rates as a function of the distance to the current price: 1) corresponds to the free gas approximation; 2) shows the same scenario with an additional drift towards the current price; 3) and 4) are like 1) and 2) with a trading effect corresponding to $`g`$ being uniformly distributed on the interval $`[3,10]`$.

Figure 5 shows the same sequence of scenarios when the overall rate of updating of trader $`i`$ behaves like $`\left(p_ip(t)\right)^{1.5}`$. Notice the increased volatility compared to the previous picture. The volatility increase before the crash is particularly marked.

## 4. Possible extensions

The basic model we describe above allows for numerous extensions to capture further important features.

### 4.1. Traders of different type.

We have assumed all traders to behave according to the same stochastic rules. It is not difficult to modify this. First, rules can be different between buyers and sellers, generalizing the bias towards the price to some more complicated function. Moreover, we could introduce different species of traders that follow different rules (e.g. optimists vs. pessimists), and study the ensuing effects. Even more challenging, we could try to introduce “intelligent agents” that try to perform arbitrage on the price by estimating future price based on observation of the past price history. This goes beyond the original ideas of the model, but could be interesting when testing some basic principles of financial mathematics in a specific controllable model context.

### 4.2. Coupling to external influences.

To capture the evolution of prices over longer time-scales, it will be indispensable to couple our model to external influences. These should reflect fundamental data on the particular stock considered (such as dividend return, earnings, cash-flow), as well as global macro-economic data (interest rates, growth rates, etc.). Such effects are easily incorporated by making the transition rates $`q(p_i,p_i^{})`$ time dependent. E.g., given the earnings at time $`t`$, one may compute a fictitious “fundamental value” (based e.g. on historic price-earnings ratios), and assume that there should be a certain tendency for market participants to adjust their subjective price towards this value. Changes in earnings (expectations) then induce a change in the transition probabilities. Similarly, other external effects exert their influence most naturally through the transition probabilities of our process. The key question of interest that our model is able to answer is how such external effects are reflected in the evolution of the price of our commodity. Addressing this question via analytical and/or numerical methods may in many respects be the most interesting and promising perspective that our models provide.

## 5. Conclusions

We have presented a class of Markov models that allow a realistic modelling of the price evolution of a commodity under trading. The basic model is a particle system like model for the evolution of a large number of traders whose state space is given by the collection of all opinions of all traders on the current value of the traded commodity. The price process is inferred from this state according to rules analogous to those used in real markets. In the simplest case of independent traders, explicit computations are possible, and we have analysed a crash scenario in a bistable market in this context.

We discussed several basic mechanisms that we think should be taken into account when modelling financial markets. These include attraction to the price, rate dependence from the distance to the price, and repulsion from the price of traders having performed an interaction. These effects lead to interesting properties of the price process which are observed in similar form in reality.

We hope to have motivated that the interacting particle systems have a place in the modelling of financial and economic systems. For this to be fruitful, this requires to chose interactions that take the special features of these systems into account. Moreover, the questions that should be addressed are quite different from what is usually done in the theory of interacting particle systems. In particular the analysis of the price process leads to rather interesting problems regarding order statistics in the particle model. As we will discuss in a forthcoming article , these problems are closely related to the study of interfaces and phase boundaries in particle systems.

## Acknowledgments

A.B. and J.C. are supported in part by the DFG-Research Center “Mathematics for key technologies”. O.H. thanks the Weierstraß Institute, Berlin, for its kind hospitality and the DFG-Research Center “Mathematics for key technologies” for financial support. A.B. thanks the Isaac Newton Institute, Cambridge (UK), for its kind hospitality and financial support within the programme “Interaction and growth in complex stochastic systems”.