text:
# Electronic spin precession in semiconductor quantum dots with spin-orbit coupling

## Abstract

The electronic spin precession in semiconductor dots is strongly affected by the spin-orbit coupling. We present a theory of the electronic spin resonance at low magnetic fields that predicts a strong dependence on the dot occupation, the magnetic field and the spin-orbit coupling strength. Coulomb interaction effects are also taken into account in a numerical approach.

In recent years the spin properties of non-magnetic semiconductors have attracted an increasing attention, not only for the fundamental physics behind the subject but also for the future technological applications of the electronic spin in spin-based devices . The available experimental techniques allow for a precise observation of spin dynamics in a wide range of semiconductor structures. Actually, spin precession can be monitored with femtosecond resolution using time resolved Faraday rotation, as reported in Ref.  for GaAs quantum wells; and in Ref.  for CdSe excitonic quantum dots. Another exciting possibility comes from spatially resolved spin detection, achieved in Ref.  for organic molecules on a graphitic surface by combining the spatial resolution of scanning-tunneling microscopy (STM) with the spin sensitivity of electron-spin resonance (ESR).

In this work we report a theoretical study of the spin precessional properties of electrons confined to a model GaAs quantum dot, including spin-orbit (SO) coupling. This mechanism gives rise to a rich variety of spin precessional frequencies, depending on the orbital state of the electrons, even in the absence of a vertical magnetic field; as compared to the Larmor frequency for systems without SO coupling.

In order to model a GaAs quantum dot with SO coupling for the standard (001) plane , we add to the Hamiltonian the Dresselhaus term originating from the bulk inversion asymmetry:

$$_D=\frac{\lambda _D}{\mathrm{}}\underset{i=1}{\overset{N}{}}\left[P_x\sigma _xP_y\sigma _y\right]_i,$$
(1)
where the $`\sigma `$’s are the Pauli matrices and $`𝐏=i\mathrm{}+\frac{e}{c}𝐀`$ represents the canonical momentum containing the vector potential $`𝐀`$ —within the symmetric gauge for a vertical magnetic field $`B`$ one has $`𝐀=\frac{B}{2}(y,x)`$. The intensity of the SO term depends on the effective dot height $`z_0`$ as  $`\lambda _D\gamma (\pi /z_0)^2`$ where $`\gamma `$ is a material dependent constant that for GaAs takes the value $`\gamma =27.5`$ eVÅ3 . In the present work we shall consider $`\lambda _D`$ parameters in the range $`[\mathrm{\hspace{0.17em}0.44},\mathrm{\hspace{0.17em}1.08}]\times 10^9\mathrm{eVcm}`$, in the same order of magnitude of those found in the literature for GaAs/AlGaAs heterojunctions  ($`0.25\times 10^9\mathrm{eVcm}`$ for electrons, $`0.6\times 10^9\mathrm{eVcm}`$ for holes).

In the following, the dot vertical extent only determines the SO coupling strength, the electronic motion is otherwise considered bidimensional in a lateral confinement potential with circular symmetry, $`V(r)`$. Taking into account the Zeeman $`_Z`$ and spatial $`_0`$ energies and neglecting for the moment electron-electron interaction, the full Hamiltonian reads $`=_0+_D+_Z`$, where

$`_0`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left[\text{}{\displaystyle \frac{𝐏^2}{2m^{}}}+V(r)\text{}\right]_i,`$ (2)
$`_Z`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left[\text{}{\displaystyle \frac{1}{2}}g^{}\mu _BB\sigma _z\text{}\right]_i.`$ (3)
In Eqs. (2) and (3) $`\mu _B`$ denotes the Bohr magneton while $`m^{}`$ and $`g^{}`$ are, respectively, the effective mass and gyromagnetic factor for the conduction band of bulk GaAs, i.e., $`m^{}=0.067m_e`$ and $`g^{}=0.44`$.

For $`_0_D_Z`$ a diagonalization in spin space to order $`\lambda _D^2`$ can be obtained by means of a unitary transformation  to a new Hamiltonian $`\stackrel{~}{}=U^+U`$, with

$`\stackrel{~}{}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}[\text{}{\displaystyle \frac{𝐏^2}{2m^{}}}+V(r)+\lambda _D^2{\displaystyle \frac{m^{}}{\mathrm{}^3}}(xP_yyP_x)\sigma _z`$ (4)
$`+`$ $`{\displaystyle \frac{1}{2}}g^{}\mu _BB\sigma _z\text{}]_iN\lambda _D^2{\displaystyle \frac{m^{}}{\mathrm{}^2}}+O(\lambda _D^3).`$ (5)
In the new, intrinsic, reference frame the eigenstates are orbitals with well defined spin and spatial angular momentum in $`z`$ direction, i.e., $`\phi _{n\mathrm{}\pm }(𝐫)\chi _\pm (\eta )`$; where $`\eta =,`$ and $`\sigma _z\chi _\pm =\pm \chi _\pm `$. Note that the spatial parts also depend on the spin label since the effective radial confinement in Eq. (4) is different for $`\chi _+`$ and $`\chi _{}`$ orbitals. When these eigenstates are transformed back to the laboratory frame, spin and angular momentum become ill defined, but they deviate little from the well defined intrinsic values. Therefore, we shall retain the intrinsic labels $`(n\mathrm{}\pm )`$ to characterize the laboratory-frame eigenstates

$`\chi _{n\mathrm{}+}(𝐫,\eta )`$ $``$ $`\phi _{n\mathrm{}+}(𝐫)\left({\displaystyle \genfrac{}{}{0pt}{}{1}{i\lambda _D\frac{m^{}}{\mathrm{}^2}re^{i\varphi }}}\right)`$ (6)
$`\chi _n\mathrm{}(𝐫,\eta )`$ $``$ $`\phi _n\mathrm{}(𝐫)\left({\displaystyle \genfrac{}{}{0pt}{}{i\lambda _D\frac{m^{}}{\mathrm{}^2}re^{i\varphi }}{1}}\right).`$ (7)
Using this approach, in Ref.  we showed that the static spin of odd-$`N`$ quantum dots alternates between up and down states as a consequence of the SO coupling when the magnetic field and/or the SO coupling strength are varied. Other static properties have been studied by Governale  using a SO coupling of the Rashba type, unitarily equivalent to the Dresselhaus contribution . Here we shall focus on the dynamical spin evolution in a model quantum dot when all the system parameters are kept fixed.

In order to excite the electronic spin precession one needs to perturb the ground state spin configuration. The usual way to achieve this consists in applying a horizontal magnetic field for a certain time interval that rotates the spin and triggers the precessional motion. By performing a spin rotation about an arbitrary horizontal axis of the above given spinors and decomposing the result in the stationary basis (6), we observe some interesting features. In absence of SO coupling ($`\lambda _D=0`$) the only allowed transitions are the spin flips $`(n\mathrm{}+)(n\mathrm{})`$, leading to an in-plane spin precession at the usual Larmor frequency $`\omega _L=|g^{}|\mu _BB/\mathrm{}`$. This is a well-known result valid even when spin-independent interactions are present . When SO coupling is considered, besides the pure spin flips other transitions involving additional changes in $`\mathrm{}`$ and/or $`n`$ are allowed. In addition to monopolar $`\delta \mathrm{}=0`$, dipolar $`\delta \mathrm{}=\pm 1`$ and quadrupolar $`\delta \mathrm{}=2`$ spin flip transitions also contribute with different weights. As we shall show below the $`\delta \mathrm{}=0`$ transitions of the pure Larmor mode are still the dominant ones in the precessional spectrum with SO coupling; the dipolar ones are weaker by more than an order of magnitude while the quadrupolar spin flip excitations turn out to be negligible in all cases studied. It is worth to point out that even transitions between orbitals with different $`n`$ could contribute because of the non orthogonality of $`+`$ and $``$ radial functions, although these will normally involve high energies and low strengths due to the small deviation from pure orthogonality.

To quantify the SO coupling effect we shall assume a parabolic confinement potential whose eigenenergies and eigenstates are analytically known:.

$$V(𝐫)=\frac{1}{2}m^{}\omega _0^2(x^2+y^2).$$
(8)
With the above analysis the lowest $`\delta \mathrm{}=0`$ spin-flip mode, that we shall call the SO precessional mode $`\omega _𝒫`$, has a frequency, in the low $`\lambda _D`$ limit,

$`\omega _𝒫=|\text{}\omega _L`$ $`+`$ $`2\mathrm{}\lambda _D^2{\displaystyle \frac{m^{}}{\mathrm{}^4}}`$ (9)
$``$ $`(2n+|\mathrm{}|+1){\displaystyle \frac{\lambda _D^2\omega _c}{\sqrt{\omega _0^2+\frac{\omega _c^2}{4}}}}{\displaystyle \frac{m^{}}{\mathrm{}^4}}\text{}|.`$ (10)
The SO precessional mode $`\omega _𝒫`$ is the dominant excitation in the spin rotation spectrum and it can be considered in a natural way as the modification of the pure Larmor mode $`\omega _L`$ by the SO coupling. Note that in Eq. (9) we have introduced the usual cyclotron frequency $`\omega _c=eB/(m^{}c)`$ and that $`\mathrm{}`$ is the angular momentum in the intrinsic frame of the precessionally active orbital, i.e., that with spin flips allowed by the Pauli principle; normally the highest or the second highest occupied level in an odd-$`N`$ dot. Even $`N`$ systems will generally not possess a net spin at low enough magnetic fields $`_D_Z`$.

Equation (9) already allows us to point out several interesting predictions: a) At $`B=0`$ the SO precessional frequency does not vanish when $`\mathrm{}0`$, with the offset indicating the SO coupling intensity. b) In general, positive and negative $`\mathrm{}`$ orbitals will display different $`B`$ dispersions for a fixed $`\lambda _D`$. c) When the precessionally active orbital changes due to an internal rearrangement the SO precessional frequency will display a discontinuity.

It is worth to mention that a $`B=0`$ offset similar to the one mentioned above was observed in GaAs quantum wells already in 1983  and that a zero-field level splitting in quantum dots was also discussed in Ref. . Note that since the SO coupling term is time reversal invariant Kramers degeneracy is preserved at $`B=0`$ and, therefore, the precessional offset can be found only when the transition involves non-conjugate states.

The validity of the preceding $`O(\lambda _D^2)`$ analysis can be tested with direct numerical calculations which avoid the approximate diagonalization procedure. To this end, we have implemented in a spatial grid the solution to the time-dependent Schrödinger equation, labelling the discrete set of orbitals by an index $`j`$,

$$i\mathrm{}\frac{}{t}\chi _j(𝐫,\eta )=\underset{\eta ^{}}{}h_{sp}(𝐫,\eta \eta ^{})\chi _j(𝐫,\eta ^{});$$
(11)
where we have defined the single-particle Hamiltonian from

$$\underset{i=1}{\overset{N}{}}h_{sp}(𝐫_i,\eta _i\eta _i^{}).$$
(12)
A ‘real time’ simulation of the precession can be performed by taking the stationary solutions to Eq. (11), rotate in spin space with a given horizontal axis (say the $`x`$ axis) and use the resulting perturbed spinors as starting point for the time evolution. Figure 1 displays one such simulation for a $`N=7`$ quantum dot. The analysis is based on the $`x`$-component of the total spin, in time (lower panel) and energy (upper panel) domains. The different features discussed above in the analytical model are nicely manifested by the numerical signals. From Fig. 1 we can also see quantitatively the strength of the SO precessional mode with respect to the doubly split upper and lower branches of the dipole modes. The minor energy differences between the numerical and analytical peak positions can be attributed to effects beyond $`O(\lambda _D^2)`$ and, also, to a slight departure from the $`_Z_D`$ limit for a finite $`B`$.

A shortcoming of the time simulation technique is found in the determination of very low energy excitations. Low frequency signals may require extremely long simulation times, exceeding the limit of computational feasibility; either by excessive computing time or by accumulated numerical error. In our case, we have estimated this limit at $`T_{\mathrm{m}\mathrm{a}\mathrm{x}}100\mathrm{ps}`$ and the corresponding minimum frequency $`\omega _{\mathrm{m}\mathrm{i}\mathrm{n}}/2\pi 10\mathrm{GHz}`$. Nevertheless, in the noninteracting case one can directly compute the perturbative strength function from the stationary ground state,

$`S_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}}(\omega )`$ $`=`$ $`{\displaystyle \underset{ij}{}}(1f_i)f_j\left|\chi _i|\sigma _x|\chi _j\right|^2`$ (13)
$`\times `$ $`\delta (\epsilon _i\epsilon _j\mathrm{}\omega ),`$ (14)
where $`i`$ and $`j`$ span the whole single particle set while the $`f_i`$’s and $`\epsilon _i`$’s give the orbital occupations and energies, respectively. We have checked that the perturbative and time simulation methods yield the same results when both are feasible, whereas the sub-10 GHz points in Fig. 2 have been computed using Eq. (13).

A systematics of the lowest peak energy, i.e., the SO precessional mode, is gathered in Fig. 2 as a function of the magnetic field for two different $`\lambda _D`$’s. Focussing first on the non-interacting results, we note that for the smaller $`\lambda _D`$ the agreement between analytical and numerical values is excellent, proving the equivalence of both methods; while the slight differences for $`\lambda _D=1.1\times 10^9\mathrm{eVcm}`$ can be understood on the basis of the previous discussion. The already mentioned offset at $`B=0`$ with respect to the Larmor frequency is clearly seen in Fig. 2 for $`N=7`$ and 11, as well as the diferent slopes for different $`\mathrm{}`$ and $`\lambda _D`$ values.

A better understanding of the precessional mode systematics is obtained from Fig. 3, which displays the level scheme as a function of the magnetic field for $`\lambda _D=1.1\times 10^9\mathrm{eVcm}`$. In this figure the active level for the same electron numbers of Fig. 2 are marked with thick dots and dashes. We note the clear correspondence of the discontinuities in Fig. 2 with the crossings in Fig. 3, which correspond to changes in the precessional level.

In the rest of the work we present numerical results for the SO precessional frequencies when the electron-electron Coulomb interaction is added to the model. In this case we rely exclusively on the real-time simulation method since the perturbative treatment equivalent to Eq. (13), known as the random-phase approximation, becomes extremely demanding in the present context of spinors without good angular momentum. Electronic exchange and correlation effects will be approximated within density-functional theory in the local-spin density approximation, as in Refs. . The Hamiltonian $`h_{sp}(𝐫;\eta \eta ^{})`$ of Eq. (12) is thus extended to include the dynamical terms

$`V_H(𝐫;t)`$ $`=`$ $`{\displaystyle \frac{e^2}{\kappa }}{\displaystyle 𝑑𝐫^{}\frac{\rho (𝐫^{};t)}{|𝐫^{}𝐫|}},`$ (15)
$`V_{\mathrm{x}\mathrm{c}}(𝐫,\eta \eta ^{};t)`$ $`=`$ $`{\displaystyle \frac{\delta E_{xc}[\rho _{\eta \eta ^{}}]}{\delta \rho _{\eta \eta ^{}}(𝐫;t)}};`$ (16)
i.e., the Hartree and exchange-correlation contributions, respectively. In Eq. (15) we have used the spin density matrix $`\rho _{\eta \eta ^{}}`$, the total density $`\rho `$ and exchange-correlation energy $`E_{xc}`$, as well as the dielectric constant $`\kappa =12.4`$.

Figure 4 is the analog of Fig. 1 within LSDA. The time signal has a similar large period, but the lower period modulations are manifestly different. Accordingly, the Fourier transform (upper panel) shows a similar low energy precessional mode but the distribution of minor peaks is rather different. The dipole peaks of Fig. 1 are washed out and, instead, new excitations at $`\mathrm{}\omega 1.7`$ and $``$ 4 meV appear. As was discussed in Ref.  in the context of the far-infrared absorption, these excitations are collective spin oscillations known as dipole magnons. Figure 2 also shows the LSDA numerical calculations for the higher SO coupling constant. The characteristics of the precessional mode discussed above are qualitatively retained within LSDA, although with the important difference that the discontinuity points are changed because of the interaction-induced orbital rearrangements.

In general, along with the transverse magnetic field the system will be probed by an electric field. This modifies the relative strength of the precessional mode with respect to the plasmon and magnon peaks of Figs. 1 and 4. We have checked this numerically by using an initial charge translation, simulating the effect of the electric field at $`t=0`$, simultaneously to the spinor rotation. The corresponding spectra display the same peaks of Fig. 1 for the non-interacting case and of Fig. 4 in LSDA, but with different heights. Therefore, only the strength, not the energy of the precessional mode depends on the coupling with the electric field.

In summary, the theory of electronic spin precession in GaAs quantum dots with SO coupling predicts a rich behaviour of the precessional frequencies with the electron number, the magnetic field and the intensity of the coupling. In this way, the spin precessional channel reveals information not only about the SO coupling intensity but also about the intrinsic level structure. It also opens the possibility to control the magnetic dynamical properties through the nanostructure parameters or, simply, by changing the number of electrons in the quantum dot, which can actually be varied one by one in GaAs nanostructures using electric gates. Although some of the precessional characteristics of GaAs dots we have discussed seem large enough to be experimentally accessible, the relevance for real samples of many effects beyond the ideal system considered here deserve more theoretical work. In particular, we mention the possible dephasing mechanisms in dot samples, influence of the temperature on the electronic precession and role of the coupling with nuclear spins.

This work was supported by Grant No. BFM2002-03241 from DGI (Spain), and by COFINLAB from Murst (Italy).


format: text

text:
# 1 INTRODUCTION

## 1 INTRODUCTION

Both from an operations and from a design perspective, the capacity of a road is an important quantity. Clearly, if demand exceeds capacity, queues will form, which represent a cost to the driver and thus to the economic system. In addition, such queues may impact other parts of the system, for example by spilling back into links used by drivers who are on a path that is not overloaded.

For a variety of reasons, however, capacity is not a deterministic fixed quantity. It is possible that one day a queue forms and the next day not, and this may even happen in spite of demand being larger on the second day. In consequence, any definition of capacity needs to take its stochastic nature into account.

For example, one could measure flow in 15-min intervals, say from 6am to 6:15am, from 6:15am to 6:30am, etc. One could then take the daily maximum of these values, and average the result over many typical workdays.

As an alternative, one could measure flow as a function of density irrespective of any other state variables. One could then obtain the average flow for each density interval, and the maximum of these flow-values would represent capacity.

All of these measurements have the property that they result in an expected value, i.e. in a number that, for a given day, can be exceeded or not be reached. In consquence, it is useful to develop models of traffic which reflect the stochastic nature of traffic. Clearly, the stochasticity can come at many different levels: demand can vary; road conditions can vary; driving behavior can vary; etc. These different contributions to stochasticity will have different influences, which need to be debated. In this paper, we want to concentrate on road capacity. We understand that there is active research to eventually include aspects of stochastic transitions into the Highway Capacity Manual (1).

The starting point for our work are simple single-lane car following models. These models are typically either of the type $`v(t+\tau _v)=f(g,\mathrm{\Delta }v,\mathrm{})`$ or of the type $`a(t+\tau _a)=h(g,\mathrm{\Delta }v,\mathrm{})`$, where $`v(t)`$ is the velocity of a car at time $`t`$, $`\mathrm{\Delta }v`$ is the velocity difference to the car ahead, and $`a`$ is the acceleration. $`g`$ is the gap to the car ahead, where $`g=\mathrm{\Delta }x\mathrm{}`$, with $`\mathrm{\Delta }x`$ the front-buffer-to-front-buffer distance, and $`\mathrm{}`$ is the space the car occupies in a jam. These models can for example be found in (2)

$$v(t+\tau _v)=V(g(t))\text{ with }V(g)=v_fv_f\mathrm{exp}(\lambda g/v_f),$$

in (3)

$$a(t)=\alpha \left(V(g(t))v(t)\right),\text{ with }V(g)=v_f(\mathrm{tanh}(g+\mathrm{})\mathrm{tanh}(\mathrm{})),$$

or in (4)

$$a(t+\tau _a)\frac{[v(t+\tau _a)]^l}{[\mathrm{\Delta }x(t)]^m}\mathrm{\Delta }v(t).$$

Additional parameters here are $`v_f`$ (the free speed), $`\lambda `$, $`l`$, and $`m`$.

When these models are implemented on a computer, they need to be discretized in time, and one has to concern oneself with the size of the integration time step, $`\mathrm{\Delta }t`$. A typical discretization is

$$a(t)=\text{given by the model}$$

$$v(t)=v(t\mathrm{\Delta }t)+\mathrm{\Delta }ta(t),$$

and

$$x(t+\mathrm{\Delta }t)=x(t)+\mathrm{\Delta }tv(t).$$

These discretizations are meant to approach the original coupled differential equations for $`\mathrm{\Delta }t0`$, and there is a whole body of literature available for this. Once time delays (via $`\tau >0`$) are introduced into such equations, numerical treatment becomes considerably more difficult, and much less is known about efficient numerical integration.

Given that state of affairs, it makes sense to look for computational models which are not based on the limit $`\mathrm{\Delta }t0`$, but which generate useful results also for relatively large time steps of, say, one second. The model that we will use in this paper has been introduced by Krauss (5); it is a variant of a model used by Gipps (6). The Krauss model has been shown to be free of collisions (i.e. that $`g<0`$ never occurs).

In addition to being crash free at large time steps, the Krauss model is also stochastic. The important parameter for our study is a noise amplitude $`ϵ`$, which we will vary from 0.5 to 2. For $`ϵ<0`$ or $`ϵ2`$ the model leaves the range of where it is plausible for traffic.

Our main results are the following:

* For medium $`ϵ`$, there are three states of traffic, which we will call laminar, mixed, and jammed. The state depends on the density. Laminar, occuring at low density, means that there are no stopped cars; jammed, occuring at high density, means that no car is driving at full speed; and mixed, occuring at intermediate density, means that the system is a mix of laminar and jammed traffic. In the mixed state, traffic is strongly inhomogeneous; the mixed state is often called the coexistence state.

It is important to note that there are three states (laminar, jammed, mixed) but only two phases (laminar, jammed). The phases refer to homogeneous sections of the system; the state refers to the system as a whole.
* For large $`ϵ`$, there is only one phase of traffic and therefore only one state. When going from low to high density, cars move closer and closer together, but traffic remains homogeneous at all times.
* At some $`ϵ`$ in between, there is a transition from the 1-phase to the 2-phase regime.
* In the Krauss model, changes of $`ϵ`$ also change the average acceleration. This is an unfortunate coincidence, and we believe that our general results regarding the number of phases are not related to this effect.
* Deterministic models, formulated either as car following models or as fluid-dynamical models, can display 1-phase or 2-phase behavior. They can however *not* display stochastic transitions between the phases.

The results are important for model building as well as for understanding field measurements. In a 2-phase model, theory predicts that there can be a hysteretic transition from the laminar to the mixed state *without a change in density*. This means that, at a given density, traffic can operate in the laminar flow state for long times, until it will eventually “break down” and switch to the mixed state. In a 1-phase model, this is impossible, and there is only one phase for any given density.

A direct consequence of this is that, if traffic follows a 1-phase model, any initial jam will “smear out” and thus eventually go away, *even with unchanged traffic conditions*. Conversely, in a 2-phase model, jams “pull themselves together” and thus stabilize themselves.

This paper starts with Sec. 2 which describes the general idea of a gas-liquid transition. Sec. 3 describes the general simulation setup including the car following model that is used, discusses space-time plots of the resulting dynamics, and investigates transients vs. the steady state. Sec. 4 then establishes how a mixed state can be nuemrically detected for a given model. Sec. 5 discusses how these results relate to deterministic model; the paper is concluded by a discussion and a summary.

## 2 PHASES IN TRAFFIC

At a first glance, our results may seem largely irrelevant to traffic operations. As already discussed above, we contend that it is not. There are thoughts about extending the Highway Capacity Manual so that it includes the concept of stochastic traffic breakdown (1). This could for example mean that, for certain flow levels, one would include a curve describing the probability that traffic flow has not broken down as a function of time. In order for such a description to make sense, the existence of spontaneous flow break-down needs to be established. Some researchers however have doubts that flow breakdown truly happens (7, 8). This paper will contribute to the theoretical understanding of the issue, by showing how one model can move from displaying breakdown to not displaying breakdown.

In fact, the analogy between a gas-liquid transition and the laminar-jammed transition of traffic was pointed out many times (e.g. (9, 3)). The description of traffic in the well-known 2-fluid-model (10) implicitely assumes the existence of two phases; and all simulation models which use spatial queues (e.g. (11, 12, 13)) will display two phases because of the definition of the dynamics. The two phases in models with queues are however much easier to understand than the phases in more realistic models.

In a gas-liquid transition, one observes the following (see also Fig. 1 left):

* In the gas phase at low densities, particles are spread out throughout the system. Distances between particles vary, but the probability of having two particles close to each other is very small.
* In the liquid phase at high densities, particles are close to each other. There is no crystalline structure as in solids, but the density is similar and in some cases (e.g. in water) even high in the liquid than in the solid phase. Because of the fact that the particles are so close to each other, it is difficult to compress the fluid any further.
* In between, there is the so-called coexistence state, where gas and liquid coexist. In typical experiments in gravity, the liquid will be at the bottom and the gas will be above it. Without gravity, as well as for example within clouds, droplets form within the gas and remain interdispersed. In clouds, small droplets will eventually merge together into bigger droplets, which will fall out of the cloud as rain. Without gravity, the droplets will just merge but never fall out. The final state of the system is having one big droplet of liquid, surrounded by gas.

The kinetics of the droplet formation (e.g. (14)) is ruled by a balance between surface tension and vapor pressure. Since surface tension pulls the droplet together, it has a tendency to push particles out of the droplet. Vapor pressure is the balancing force – it pushes particles into the droplet.

Surface tension depends on the droplet radius – the smaller the droplet, the larger the surface tension. The result is that, when coming from small densities, there is a regime where large droplets would already be stable, but small droplets are not. In equilibrium, this density would be in the coexistence state, but coming from low density, the homogeneous gaseous phase can survive for some time. This way, meta-stability can be explained in the gas-liquid transition, and by extension also in traffic models.

When compressing the system beyond this point, at density $`\rho _1`$, the following can be observed: Pressure will *not* go up any further. Instead, more droplets will form, and the density outside the droplets will remain constant. Pressure will not go up until all the space is used up by droplets, i.e. all particles are moved into the liquid phase. Only at this point, at density $`\rho _2`$, will pressure go up again when compressing the system.

The state between $`\rho _1`$ and $`\rho _2`$ is called coexistence state. In it, gas and liquid coexist, and changes on the global density are absorbed by the fraction of space and particles that the liquid occupies. The interface between droplets and gas remains the same no matter what the global density is.

In addition to the coexistence, when keeping density constant one will observe coagulation as already mentioned earlier: The many original small droplets will coagulate to fewer and fewer large droplets; in the limit there will only be one large droplet.

A direct consequence of the above behavior is hysteresis:

* When coming from low densities, it is possible to go beyond $`\rho _1`$ and still remain in the gaseous phase. Only after some waiting time will, by a fluctuation, some particles get close enough to each other to start the formation of a droplet.
* When coming from densities above $`\rho _1`$, it is possible that the droplet survives for some time even below $`\rho _1`$.

When increasing temperature with the above model, the 2-phase structure will eventually go away. This happens via $`\rho _1`$ and $`\rho _2`$ approaching each other and eventually meeting. That is, depending on the temperature, a fluid system will either display transitions from gas to coexistence and from coexistence to liquid, or there will be *no transition at all* (Fig. 1 right).

We will now move on to describe the supporting evidence for our claims. As is typical in computational statistical physics, our evidence is based on computer simulations. It is however backed up by generic knowledge about phase transitions as they are well understood in physics.

## 3 THE SIMULATIONS

### 3.1 Krauss Model

The velocity update of the Krauss model (15, 16) reads as follows:

$`v_{\mathrm{safe}}`$ $`=`$ $`\stackrel{~}{v}(t)+{\displaystyle \frac{g(t)\stackrel{~}{v}(t)\tau }{\overline{v}/b+\tau }}`$ (1)
$`v_{\mathrm{des}}`$ $`=`$ $`\mathrm{min}\{v(t)+a\mathrm{\Delta }t,v_{\mathrm{safe}},v_{\mathrm{max}}\}`$ (2)
$`v(t+\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{max}\{0,v_{\mathrm{des}}ϵa\eta \}.`$ (3)
$`\stackrel{~}{v}`$ is the speed of the car in front, $`\overline{v}=(v+\stackrel{~}{v})/2`$ is the average velocity of the two cars involved, $`v_{\text{max}}`$ is the maximum allowed velocity, $`a`$ is the maximum acceleration of the vehicles, $`b`$ their maximum deceleration for $`ϵ=0`$, $`ϵ`$ is the noise amplitude, and $`\eta `$ is a random number in $`[0,1]`$. The meaning of the terms is as follows:

* Eq. (1): Calculation of a “safe” velocity. This is the maximum velocity that the follower can drive when she wants to be sure to avoid a crash. The main assumption is that the car ahead will never decelerate faster than $`b`$, and that the car of the follower can also decelerate with up to $`b`$.
* Eq. (2): The desired velocity is the minimum of: (a) current velocity plus acceleration, (b) safe velocity, (c) maximum velocity (e.g. speed limit).
* Eq. (3): Some randomness is added to the desired velocity.

After the velocities of all vehicles are updated, all vehicles are moved.

The deterministic limit $`ϵ=0`$ of the Krauss model has been proven to be free of crashes for numerical time steps $`\mathrm{\Delta }t`$ smaller than or equal to the reaction time, $`\tau `$. We will use $`\mathrm{\Delta }t=\tau =1`$ as has conventionally been used for the Krauss model. Even for $`ϵ>0`$, as was used for our studies, we never observed vehicles getting closer than their minimum distance. We further use $`a=0.2`$, $`b=0.6`$, $`v_{\text{max}}=3`$ for all simulations.

The model is free of units; this is a property that it has inherited from the cell-based cellular automata models. A reasonable calibration is: time steps correspond to seconds and cells correspond to $`7.5`$ meters. The reaction time then is assumed to be $`1`$ sec, and $`v_{\text{max}}=3`$ corresponds to $`22.5`$ m/s or $`81`$ km/h. $`a=0.2`$ corresponds to a maximum acceleration of $`1.5`$ m/s per second or $`5.4`$ km/h per second. $`b=0.6`$ corresponds to a maximum deceleration of $`16.2`$ km/h per second.

All simulations are done in a 1-lane system of length $`L`$ with periodic boundary conditions (i.e. the road is bent into a ring). Let $`N`$ be the number of cars on the road. The (global) density is $`\rho =\frac{L}{N}`$.

### 3.2 Pictures

Before analysing the Krauss-model numerically, it is instructive to look at the space-time plots in Fig. 1. Space-time plots are pictures of the time evolution of the system. In Fig. 1, vehicles drive to the right and time goes down. Each row of pixels is a “snapshot” of the state of the road. In principle, one can reconstruct the trajectory of a particular car by connecting the corresponding pixels. At the displayed resolution this is however close to impossible and it is mostly the larger scale traffic jam structure that one observes. Traffic jams move *against* the direction of driving.

* The laminar state: All cars drive at high speed. The available space is shared evenly among the cars. The traffic is very homogeneous.
* The mixed state: The slow cars are all together in one big jam. On the rest of the road, the cars drive at high speed. The traffic is very inhomogeneous.
* The jammed state: The density is so high that not a single car can drive fast. As in a), the traffic is very homogeneous.
* The single phase at high $`ϵ`$: Many small jams are distributed over the whole system. There is neither a larger area of free flow, nor a major jam. The traffic is homogeneous.

Note that “homogeneous” here means “homogeneous on large scales”. What this means is that there is a spatial measurement length $`\mathrm{}`$ above which all density measurements return the same value.111More precisely, the fluctuations from one measurement to the next are the same as in case (i).  If a system goes from a 2-phase to a 1-phase model, then even in the regime which technically has only one phase, structure formation on small scales is still possible. Fig. 1 bottom right is indeed an example for this. Only with larger distance from the 2-phase model will these structures go away.

### 3.3 Defining a jam

In order to make progress, one needs to define where a jam starts and where it ends. Our definition of homogeneity (see later) will not depend on this, however. A *jam* is a sequence of cars driving with speed less or equal $`v_{\text{max}}/2`$. The cars between two neighbouring jams are in *laminar* flow.

This definition is very simple, but will not always correspond to our natural understanding of the word jam. Thus, whether a car is jammed or not according to this definition is just a starting point and not the final answer.

### 3.4 Initial Condition And Relaxation

For many parameters of the Krauss model, there is a unique equilibrium state, which the system will attain after a finite time $`t_{\text{relax}}`$, no matter how it was started. However, deciding when the equilibrium is reached is not trivial (running the simulation for $`t\mathrm{}`$ clearly is not an option).

Let $`r_t`$ be the state of the road at time $`t`$. To find the equilibrium value of some property, $`E[f(r_{t_{\text{relax}}})]`$, we use the following idea: For small $`t`$, $`E[f(r_t)]`$ will depend on the initial condition. With increasing time, $`E[f(r_t)]`$ converges towards the equilibrium value. Assume the convergence is from above. Now we need another initial condition that approaches the equilibrium value from below. Once these two sequences are close enough together, the equilibrium value is found. Unfortunately, it cannot be guaranteed that the value thus obtained really is the equilibrium value.

We use the following two initial conditions:

* laminar: The cars are positioned equidistant over the road with speed zero.
* jammed: all cars are cramped together in a big jam without any gap. Their speed is zero.

An example of this method is shown in Fig. 2. For $`f(.)`$ the number of jams was used. Since both initial conditions start with $`v=0`$, the criterion finds one large jam. Vehicles then accelerate, but because of interaction will form small jams. For that reason, the laminar start leads to many jams very quickly. From then on, the number of jams goes down, because jams coagulate. In contrast, when starting with a large single jam, than that jam remains the only one in the system for large times. In Fig. 2, we see that for $`ϵ=1.0`$, the system eventually goes to a state where it has, in the average, about 1.8 jams. In contrast, with $`ϵ=1.5`$, the system converges to an average of more than 20 jams. Also, the figure shows that the system goes to those long-run states no matter how it starts.

## 4 ESTABLISHMENT OF A PHASE DIAGRAM VIA A MEASURE OF INHOMOGENEITY

In this section, a criterion is established that distinguishes homogeneous from mixed states. As pointed out before, mixed states are characterized by the coexistence of laminar and jammed traffic. As we have also said, if we wait long enough, the phases will coagulate, leading to exactly one laminar and one jammed section in the system. When approaching the boundaries of the mixed regime, this characterization will become less clear-cut, and it may be possible to have more than one jam. Typically, there will be one major jam and many small ones, but for many measurement criteria this will cause enough problems to no longer be able to differentiate between the mixed and a homogeneous state. Our criterion will only show a gradual decrease in differentiating power.

As already discussed earlier, it should be noted that some states that are called “homogeneous” in this paper may appear inhomogeneous to an observer. An example for this is Fig. 1 bottom right. As said before, these states are “homogeneous on large scales”, which is the important criterion here. Essentially, this means that for system size $`L\mathrm{}`$ and measurement interval $`\mathrm{}\mathrm{}`$ (but $`\mathrm{}L`$), all density measurements will eventually return the same value. This will not be the case for “mixed” states.

The criterion is defined as follows: Partition the road into segments of length $`\mathrm{}`$ (for simplicity let $`\mathrm{}`$ divide $`L`$ without remainder). For each segment the local density $`\rho _{\mathrm{}}`$ can be computed as the number of cars in that segment divided by $`\mathrm{}`$. An interesting value is the variance of the local density:

$$\mathrm{Var}[\rho _{\mathrm{}}]=\frac{1}{L/\mathrm{}}\underset{i=1}{\overset{L/\mathrm{}}{}}(\rho _{\mathrm{}}(i)\mathrm{E}[\rho _{\mathrm{}}])^2,$$
(4)
where $`E[.]`$ is the expected value, which in our case is the same as the systemwide density. Note that since the density lies within $`[0,1]`$, the variance cannot exceed $`1/4`$.

What this value picks up is how much each individual measurement segment of length $`\mathrm{}`$ deviates, in terms of its density, from the average density. Assume a system consisting of jammed and laminar traffic. If there is a jam in one segment, then the segment’s density will be much higher than the average density. Conversely, if there is only laminar traffic in a segment, then the segment’s density will be much lower than the average density. $`\mathrm{Var}[\rho _{\mathrm{}}]`$ takes the average over the square of these deviations.

Fig. 3 shows this value as a function of the global density $`\rho `$ and the noise parameter $`ϵ`$. Look at it for fixed $`ϵ`$, say $`ϵ=1`$. One sees that at densities up to $`\rho 0.2`$, the value of $`\mathrm{Var}[\rho _{\mathrm{}}]`$ is close to zero, indicating a homogeneous state, which is in this case the laminar state. Similarly, for densities above $`\rho 0.8`$, $`\mathrm{Var}[\rho _{\mathrm{}}]`$ is again close to zero, indicating a homogeneous state, which is in this case the jammed state. In between, for $`0.2\rho 0.8`$, the value of $`\mathrm{Var}[\rho _{\mathrm{}}]`$ is significantly larger than zero, indicating a mixed state.

Now slowly increase $`ϵ`$. We see that the laminar regime ends at smaller and smaller densities, while the jammed regime starts at smaller and smaller densities. The latter means that for large $`ϵ`$, the jammed phase has many relatively small holes, which reduce the density, but do not break the jam. At $`ϵ1.7`$, the mixed phase completely goes away; for larger $`ϵ`$, we do not pick up any inhomogeneity at *any* density. Compare this to the theoretical expectation in Fig. 1, where for increasing $`T`$ the two densities eventually merge and thus the different phases go away. Note that close to the transition the system still looks like it possesses different phases (see Fig. 1(iv) and locate the corresponding $`ϵ=1.8`$ and $`\rho =0.2`$ in Fig. 3). These structures do however exist *on small scales only*; when averaging over larger segments, then all segments contain exactly the same density. A segment length of $`\mathrm{}=62.5`$, as used in the figure, is already sufficient in order to not measure any inhomogeneity for the state in Fig. 1(iv).

Remember again that $`ϵ`$ is a model parameter while $`\rho `$ is a traffic state. That is, once one has settled for an $`ϵ`$, the model behavior is fixed, and one has decided if one can encounter a second phase or not. *If* one can encounter a second phase, it will come into existence through changing traffic demand throughout the day – traffic can move from the laminar into the mixed and potentially into the jammed state and back.

As a side remark, let us note that there is also another 1-phase regime for $`ϵ0`$. Albeit potentially interesting, this is outside the scope of this paper.

The maximum of the density variance is at $`\rho 0.5,ϵ1`$. This can be explained as follows: $`ϵ=1`$ produces a sharp separation of the jam (high density) and the laminar structure (low density). With $`\rho =0.5`$ it turns out that these two phases have the same length, thus $`\mathrm{Var}[\rho _{\mathrm{}}]`$ is maximal. Increasing (decreasing) $`\rho `$ will increase (decrease) the length of the jam and therefore in both cases decreases $`\mathrm{Var}[\rho _{\mathrm{}}]`$.

The elliptical shape with its diagonal axis (instead of a vertial axis as in Fig. 1) can be explained as follows: Increasing $`ϵ`$ decreases, via Eq. 3, the acceleration. This leads to lower density in the laminar region. Since $`\mathrm{Var}[\rho _{\mathrm{}}]`$ is maximal when both phases have the same length, and since a lower density in the laminar region makes the jam occupy more space, one needs to reduce the density in order to go back to the state where they occupy equal space.

In summary, one obtains, for the traffic model, a phase diagram as in Fig. 1, which is the schematic phase diagram for a gas-liquid transition in fluids. Again, the important feature of this phase diagram is that there are three states (laminar/gas; mixed/coexistence; jammed/liquid) for low temperatures. For higher temperatures, the coexistence range becomes more and more narrow, while the density of the gas phase and the density of the liquid phase approach each other. Eventually, they become equal, and the coexistence state dies out. The only important difference is that for our traffic model the phase diagram is bent to the left with increasing $`ϵ`$.

There are other criteria which can be used to understand these types of phase transitions. In particular, one can look at the gap distribution between jams, and one would expect a fractal structure at the point where the 2-phase and the 1-phase model meet, i.e. at $`\rho 0.2`$ and $`ϵ1.7`$. This is still under investigation but looks promising. This should be followed by an investigation in how far cell-based traffic models (17) models display a similar transition.

## 5 PHASE TRANSITIONS IN DETERMINISTIC MODELS

Only stochastic models can display spontaneous transitions between homogeneous and coexistence ($`=`$ mixed) states. The nature of the transition can however also become clear in deterministic models. We will discuss these similarities first for a deterministic car following model and then for deterministic fluid-dynamical models.

### 5.1 Car Following Models

For the model of Eq. (1), it has been shown (18) that the homogeneous solution of the model is linearly stable for all densities for $`V^{}<\alpha /2`$, and linearly unstable for certain densities for $`V^{}>\alpha /2`$, where $`V^{}`$ is the first derivative of the function $`V(.)`$. The instability sets in for certain intermediate densities, that is, for low and high densities *all* models are stable in the homogeneous (laminar or jammed) state. For intermediate densities, one can select the curve $`V(g)`$ and the parameter $`\alpha `$ such that the model is either stable or not stable.

If all parameters including the density are such that the homogeneous solution is not stable, then the system rearranges itself into a pattern of stop-and-go traffic, corresponding to the coexistence state. The density of the laminar and the jammed phase in the coexistence state are independent from the average system density, that is, if in that state system density goes up, it is reflected in the jammed phase using up a larger fraction of space.

The type of the instability is similar to the better-known instability of Eq. (1). However, once the instability is triggered in Eq. (1), it will just grow exponentially, and no stable 2-phase solution is found (e.g. (19)).

### 5.2 Fluid-Dynamical Models

Standard Lighthill-Whitham theory, of the type

$$_t\rho +_xQ(\rho )=0$$

with a strictly convex flow-density-curve $`Q(\rho )`$, results in a 1-phase model. When $`Q(\rho )`$ has linear sections, then in that range shock waves are marginally stable, in the sense that disturbances to those shocks are neither amplified nor dissipated away. (In a 1-phase model, disturbances are dissipated away and the final state is always homogeneous; in a 2-phase/3-state model, there is a density regime where disturbances to shock waves are dissipated away in the sense that the shock wave will regain its standard profile again.) When $`Q(\rho )`$ has concave sections, then the situation becomes more complicated (e.g. (20)).

Fluid-dynamical theory, of the type

$$_t\rho +_x(\rho v)=0$$

and

$$_tv+v_xv=\frac{1}{\tau }\left(V(\rho )v\right)+\alpha (\rho )_x\rho +\nu (\rho )_x^2v$$

*can*, depending on the choice of parameters including the $`V(\rho )`$-curve, either be a 1-phase/1-state or a 2-phase/3-state model (21).

As pointed out before, these models are deterministic, so in no situation will these models display stochastic transitions.

## 6 DISCUSSION

There is no general agreement if measurements show 1-phase/1-state or 2-phase/3-state traffic (or possibly even three phases (22)). There is some evidence for hysteresis in Germany (23), manifesting itself in transitions from high to lower flow values at the same density. Hysteresis, which was also found earlier (24), is a strong indication for a 2-phase model. However, even in Germany, most measurements indicate highly variable traffic at intermediate densities, which does not correspond to any clear-cut picture.

In this context, one should note that a 1-phase model which is close to a 2-phase model would also display highly variable traffic at intermediate densities, although it would be homogeneous at large scales as discussed in Sec. 4. This variability is however a property of stochastic models only and for that reason it is not well integrated into current theory development. A precise investigation of these relations is beyond the scope of this paper. It seems however impossible to us to clarify the question if traffic displays several phases or not – and therefore, if breakdown probability should be entered into the Highway Capacity Manual or not – without having understood how different phases are generated by stochastic models. The present paper fills exactly this gap.

## 7 SUMMARY

This paper shows, via numerical evidence, that a specific stochastic car following model can either display 1-phase/1-state or 2-phase/3-state traffic, depending on the choice of parameters. With 2-phase parameters, the two phases are laminar and jammed, which also corresponds to two of the three states. Those states are homogeneous. The third state, at intermediate densities, is a coexistence or mixed state, consisting of sections with jammed and sections with laminar traffic.

The transition to a 1-phase/1-state model happens via the densities of the laminar and of the jammed phase approaching each other until they become the same. Beyond this point, there is only one homogeneous phase of traffic.

Some of these findings can be understood by looking at deterministic models for traffic, either car-following or fluid-dynamical. However, the stochastic elements of the transition cannot be explained by deterministic models. An important stochastic element is that structure formation and strong variability can also happen in a 1-phase model as long as the parameters are close to the 2-phase model – a deterministic model would converge to a homogeneous solution here.

In our view, it is important to understand this possibility of stochastic models to be in different regimes if one considers to enter discussions of traffic breakdown probabilities into the Highway Capacity Manual. If traffic is best described by a 1-phase model, then there is, in our view, no theoretical justification for such probabilities. If, however, traffic is best described by a 2-phase model, then the 2-phase model could even give theoretical predictions for breakdown probabilities. A discussion of breakdown probabilities in 2-phase models can be found in Ref. (25).

## ACKNOWLEDGMENTS

We thank F. Wagner for intensive discussions. Also D. Stauffer and M. Droz have significantly helped in defining good computational criteria for the detection of multiple phases. P. Wagner and C. Daganzo have helped indirectly (and maybe inadvertently), via the discussion of a separate but related paper. A significant portion of the computing times on 192-CPU cluster Xibalba and on 44-CPU cluster Linneus was used for some of the computational results.