text:
# Statistical Physics of Traffic Flow

## 1 Introduction

Despite the long history of application of methods from physics to traffic flow problems (going back to the fifties) it has only recently blossomed into a successfull field of “exotic statistical physics” . Until a few years ago most approaches were based on ”classical” methods from physics, especially from mechanics and hydrodynamics. In general one can distinguish microscopic and macroscopic approaches.

In microscopic models individual vehicles are distinguished. A typical example are the so-called car-following theories . For each car one writes an equation of motion which is the analogue of Newton’s equation. The basic philospophy of the car-following approach can be summarized by the equation $`[Response]_n=\kappa _n[Stimulus]_n`$ for the $`n`$th vehicle. Each driver $`n`$ responds to the surrounding traffic conditions which constitute the stimulus for his reaction. The constant of proportionality $`\kappa _n`$ is also called sensitivity. Usually also the reaction-time of the drivers is taken into account. A typical example is the follow-the-leader model  where the stimulus is given by the velocity difference to the next car ahead. Assuming that drivers tend to move at the same speed as the leading car the equations of motion are given by $`\ddot{x}_n(t)=\kappa _n\left[\dot{x}_{n+1}(t)\dot{x}_n(t)\right]`$. In order to obtain realistic behaviour the sensitivity $`\kappa _n`$ has to become a function of the velocity and the distance between the cars, $`\kappa _n=\kappa _n(v_n,x_{n+1}x_n)`$. In recent years the so-called optimal-velocity model  has successfully been used. Here the acceleration is determined by the difference of the actual velocity $`v_n(t)`$ and an optimal velocity $`V^{opt}(\mathrm{\Delta }x_n)`$ which depends on the distance $`\mathrm{\Delta }x_n`$ to the next car. The equations of motion are then of the form $`\ddot{x}_n(t)=\kappa _n[V^{opt}(\mathrm{\Delta }x_n)v_n(t)]`$.

In macroscopic models one does not distinguish individual cars. Instead a ”coarse-grained” fluid-dynamical description in terms of densities $`c(x,t)`$ and flows $`J(x,t)`$ is used. Traffic is then viewed as a compressible fluid formed by the vehicles. Density and flow are related through a continuity equation which for closed systems takes the form $`c/t+J/x=0`$. Since this is only one equation for two unknown functions one needs additional information. In one of the first traffic models Lighthill and Whitham  assumed that $`J(x,t)`$ is determined by $`c(x,t)`$, i.e. $`J(x,t)=J(c(x,t))`$. Inserting this assumption into the continuity equation yields the so-called Lighthill-Whitham equation $`c/t+v_gc/x=0`$ with $`v_g=dJ/dc`$. However, for a more realistic description of an traffic additional equation, the analogue of the Navier-Stokes equation for fluids, describing the time-dependence of the velocities $`v_n(x,t)`$ has to be considered instead of the simple Lighthill-Whitham assumption.

The first approach borrowing ideas from statistical physics is the kinetic theory of vehicular traffic . Here traffic is treated as a gas of interacting particles. The interactions are described by a dynamical equation in phase space which is the analogue of the Boltzmann equation in the kinetic theory of gases.

This is only a a very brief overview over the different approaches for the description of traffic flow. For a more complete exposition of the models and their merits we refer to the literature, e.g. . Nice overviews over the physics of related nonequilibrium models can be found in . In the following we concentrate on a class of models developed in the last few years, namely (probabilistic) cellular automata. These models describe stochastic processes and are therefore of general interest for statistical physics . An extensive review of the application of cellular automata to traffic flow and related systems can be found in . Here we focus on aspects that have not been discussed in much detail in , e.g. the effects of on- and off-ramps.

## 2 Some empirical facts

The most simple empirical facts that should be reproduced by any traffic model are the spontaneous formation of jams and the characteristic form of the flow-density relation, the so-called fundamental diagram. The space-time plot of Fig. 1 shows the formation and propagation of a traffic jam. In the beginning, the vehicles are well separated from each other. Then, without any obvious reason like an accident or road construction, a dense region appears due to fluctuations. This finally leads to the formation of a jam which remains stable for a certain period of time but disappears again without any obvious reason. During its lifetime the front of the jam moves backwards against the driving direction of the cars with a characteristic velocity of approximately 15 km/h.

A more detailed analysis of traffic jams in absence of hindrances has been given by Kerner and Rehborn . They found that the upstream velocity and, therefore, the outflow from a jam is approximately constant. The outflow from a jam and the velocity of the jam fronts are now regarded as two important empirical parameters of highway traffic which can be used for calibrating theoretical models.

Fig. 2 shows a typical fundamental diagram obtained from empirical data. One can clearly distinguish two regions, the free-flow regime at small densities and the congested regime at high densities. In the free-flow regime the cars do not hinder each others motion and therefore the flow grows linearly with density. In the congested regime however, the motion of the vehicles is dominated by the hindrance due to other cars. Here the flow decreases with increasing density.

These considerations only give the basic structure of the fundamental diagram. In Sec. 4 a more detailed analysis of empirical data will be presented which reveals a much richer structure, e.g. the existance of hysteresis and metastable states. Furthermore in the congested regime one can distinguish two different phases, namely stop-and-go traffic and synchronized flow.

## 3 Cellular automata models

The Nagel-Schreckenberg (NaSch) model  is a probabilistic cellular automaton (CA) able to reproduce the basic features of traffic flow as discussed in Sec. 2. In a CA not only space and time are discrete, but also the state variable. In the NaSch model, each cell can be empty or occupied by exactly one car (Fig. 3). The state of the car is then characterized by its velocity $`v`$ which can take one of the $`v_{max}+1`$ allowed integer values $`v=0,1,\mathrm{},v_{max}`$. Suppose, $`x_n`$ and $`v_n`$ denote the position and velocity, respectively, of the $`n`$-th vehicle. Then, $`d_n=x_{n+1}x_n`$, is the gap in between the $`n`$-th car and the car in front of it at time $`t`$. At each time step $`tt+1`$, the arrangement of the $`N`$ cars on a finite lattice of length $`L`$ (i.e. for a global density $`\rho =N/L`$) is updated in parallel according to the following ”rules”:

Step 1: Acceleration.
If $`v_n<v_{max}`$, the speed of the $`n`$-th car is increased by one, i.e.

$$v_n\mathrm{min}(v_n+1,v_{max})$$
$`(\mathrm{U1}).`$
Step 2: Deceleration (due to other cars).
If $`d_nv_n`$, the speed of the $`n`$-th car is reduced to $`d_n1`$, i.e.,

$$v_n\mathrm{min}(v_n,d_n1)$$
$`(\mathrm{U2}).`$
Step 3: Randomization.
If $`v_n>0`$, the speed of the $`n`$-th car is decreased randomly by unity with probability $`p`$, i.e.,

$$v_n\mathrm{max}(v_n1,0)\mathrm{with}\mathrm{probability}p$$
$`(\mathrm{U3}).`$
Step 4: Vehicle movement.
Each car is moved forward according to its new velocity determined in Steps 1–3, i.e.

$$x_nx_n+v_n$$
$`(\mathrm{U4}).`$
The NaSch model is a minimal model in the sense that all the four steps are necessary to reproduce the basic features of real traffic; however, additional rules are needed to capture more complex situations. Step 1 reflects the general tendency of the drivers to drive as fast as possible without crossing the maximum speed limit. Step 2 is intended to avoid collision between the cars. The randomization in step 3 takes into account the different behavioural patterns of the individual drivers, especially, nondeterministic acceleration as well as overreaction while slowing down; this is crucially important for the spontaneous formation of traffic jams. Even changing the precise order of the steps of the update rules stated above would change the properties of the model. E.g. after changing the order of steps 2 and 3 there will be no overreactions at braking and thus no spontaneous formation of jams (see below). The NaSch model may be regarded as stochastic CA . In the special case $`v_{max}=1`$ the deterministic limit of the NaSch model is equivalent to the CA rule $`184`$ in Wolfram’s notation .

The use of a parallel updating scheme (instead of a random-sequential one) is crucial  since it takes into account the reaction-time and can lead to a chain of overreactions. Suppose, a car slows down due the randomization step. If the density of cars is large enough this might force the following car also to brake in the deceleration step. In addition, if $`p`$ is larger than zero, it might brake even further in Step 3. Eventually this can lead to the stopping of a car, thus creating a jam. This simple mechanism of spontaneous jam formation is rather realistic and cannot be modeled by the random-sequential update.

Space-time diagrams showing the time evolutions of the NaSch model demonstrate that no jam is present at sufficiently low densities, but spontaneous fluctuations give rise to traffic jams at higher densities (Fig. 4(a)). From Fig. 4(b) it should be obvious that the intrinsic stochasticity of the dynamics, arising from non-zero $`p`$, is essential for triggering the jams. For a realistic description of highway traffic , the typical length of each cell should be about $`7.5`$ m which is the space occupied by a car in a dense jam. When $`v_{max}=5`$ each time step corresponds to approximately $`1`$ sec of real time which is of the order of the shortest relevant timescale in real traffic, namely the reaction time of the drivers.

### 3.1 Open boundary conditions

For $`v_{max}=1`$ the NaSch model reduces to the totally asymmetric simple exclusion process (ASEP) with parallel update. The ASEP is the simplest prototype model of interacting systems driven far from equilibrium . In the ASEP (Fig. 5) a particle can move forward one cell with probability111Note that conventionally the hopping rate in the ASEP is denoted as $`p`$. Since in the NaSch model $`p`$ is the braking probability the hopping rate in the ASEP (for $`v_{max}=1`$) becomes $`q=1p`$. $`p`$ if the lattice site immediately in front of it is empty. In addition one considers open boundary conditions. If the first cell is empty a particle will be inserted there with probability $`\alpha `$. If the last cell is occupied the particle will be removed with probability $`\beta `$.

By varying the boundary rates $`\alpha `$ and $`\beta `$ one obtains a surprisingly rich phase diagram (see Fig. 6) which is qualitatively the same for all types of dynamics . Three phases can be distinguished by the functional dependence of the current through the system on the system parameters. In the low-density phase A ($`\alpha <\beta ,\alpha _c(p)`$) the current is independent of $`\beta `$. Here the current is limited by the rate $`\alpha `$ which then dominates the behaviour of the system. In the high-density phase B ($`\beta <\alpha ,\beta _c(p)`$) the behaviour is dominated by the output rate $`\beta `$ and the current is independent of $`\alpha `$. In the maximum current phase C ($`\alpha >\alpha _c(p)`$ and $`\beta >\beta _c(p)`$) the limiting factor for the current is the bulk rate $`q=1p`$. Here the current becomes independent of both $`\alpha `$ and $`\beta `$.

The transitions between the phases can be characterized by the behaviour of two correlation lengths $`\xi _\alpha `$ and $`\xi _\beta `$ which only depend on $`p`$ and $`\alpha `$ or $`\beta `$. Apart from $`\xi _\alpha `$ and $`\xi _\beta `$ also a third length $`\xi ^1=|\xi _\alpha ^1\xi _\beta ^1|`$ plays an important role.

The transition from A (B) to C is continuous with diverging correlation length $`\xi _\alpha `$ ($`\xi _\beta `$). The transition from the high- to the low-density phase is of first order. Here both $`\xi _\alpha `$ and $`\xi _\beta `$ are finite, but $`\xi `$ diverges. On the transition line one finds a linear density profile created by the diffusion of a domain wall between a low-density region at the left end of the chain and a high-density region at the right end.

In  a nice physical picture has been developed which explains the structure of the phase diagram not only qualitatively, but also (at least partially) quantitatively. The phase diagram of the open system is completely determined by the fundamental diagram of the periodic system through an extremal-current principle  and therefore independent of the microscopic dynamics of the model.

Fig. 7 shows the full phase diagram of the NaSch model with open boundary conditions and $`v_{max}=4`$. As predicted by the extremal-current principle the same phases as in the case $`v_{max}=1`$ described above appear. Instead of using input- and output-rates at the boundaries here reservoirs have been introduced which induce effective densities $`\rho ^{}`$ and $`\rho ^+`$ near the boundaries. In order to allow for an easier comparison with empirical data obtained near ramps (see Sec. 5.2) the phases are denoted as ”free-flow” and ”congested” instead of ”low-density” and ”high-density”.

## 4 Generalizations and extensions of the NaSch model

### 4.1 Metastability, hysteresis and slow-to-start rules

Measurements on real traffic have revealed that traffic flow can exhibit metastability and hysteresis effects . Fig. 8(a) shows a schematic fundamental diagram with a metastable branch. In the density interval $`\rho _1<\rho <\rho _2`$ the flow is not uniquely determined but depends on the history of the system. Here a metastable high-flow branch exists. Closely related is the observation of hysteresis effects. Fig 8(b) shows a clear hysteresis loop obtained from time-traced measurements of the flow .

The NaSch model in its simplest form as introduced in Sec. 3 does not exhibit metastable states and hysteresis. However, a simple generalization exists which is able to reproduce these effects. It is the so-called Velocity-Dependent-Randomization (VDR) model . Here, in contrast to the original NaSch model, the randomization parameter depends on the velocity of the car, $`p=p(v)`$. The rules of Sec. 3 are supplemented by a new rule,

Step 0: Determination of the randomization parameter. The randomization parameter used in step 3 for the $`n`$-th car is given by $`p=p(v_n(t))`$.

This new step has to be carried out before the acceleration step 1. The randomization parameter used in step 3 depends on the velocity $`v_n(t)`$ of the $`n`$th car after the previous timestep. Most relevant for an understanding of the occurance of metastability and hysteresis are simple slow-to-start rules where one chooses

$$p(v)=\{\begin{array}{cc}p_0\hfill & \text{for }v=0\text{,}\hfill \\ p\hfill & \text{for }v>0\text{,}\hfill \end{array}$$
(1)
with $`p_0>p`$. This means that cars which have been standing in the previous timestep have a higher probability $`p_0`$ of braking in the randomization step than moving cars.

Fig. 9 shows a typical fundamental diagram obtained from computer simulations of the VDR model. It can be clearly seen that it has the same form as Fig. 8(a).

Over a certain density interval $`J(\rho )`$ can take one of the two values depending on the initial state and, therefore, exhibit metastability. The microscopic structure of the states in the free-flow branch is very similar to that of the NaSch model. The metastable homogeneous states have a lifetime after which their decay leads to congested branch (see Fig. 10). In contrast, the microscopic structure of the congested branch differs drastically from the NaSch model. Here one observes phase separation into a macroscopic jam and a macroscopic free-flow regime (see Fig. 10).

It is instructive to compare the fundamental diagram of the VDR model with those of the corresponding NaSch models. The fundamental diagram of the VDR model can be understood from heuristic arguments utilizing the observed structures of the steady-states. For small densities $`\rho 1`$ there are no slow cars in the VDR model since interactions between cars are extremely rare. Every car moves with the free-flow velocity $`v_f=v_{max}p`$ and, therefore, the flux is given by

$$J_{hom}(\rho )=\rho (v_{max}p)$$
(2)
which is identical to the NaSch model with randomization $`p`$. On the other hand, for densities close to $`\rho =1`$, the cars are likely to have velocities $`v=0`$ or $`v=1`$ only and, therefore, the random braking is dominated by $`p_0`$. A simple waiting time argument  shows that the flow in the phase-separated regime is given by

$$J_{sep}(\rho )=(1p_0)(1\rho ).$$
(3)
This corresponds to the NaSch model with randomization $`p_0`$ (see Fig. 9).

The behaviour found in the VDR model is generic for models with slow-to-start rules. Due to the slow-to-start rule the outflow from a large jam is smaller than the maximal possible flow, i.e. the maximum of the fundamental diagram. Therefore the density far downstream is rather small and the vehicles propagate almost freely in the low-density region. The spontaneous formation of jams is suppressed if $`p`$ is not too large. This is the basic mechanism which leads to the formation of the phase-separated state and hysteresis.

### 4.2 Anticipation and synchronized traffic

In the NaSch model as well as in the VDR model drivers only react to the distance to the next car ahead, i.e. the headway. In real traffic often “anticipation” plays an important role. Drivers not only react to a change of their headway, but also take into accout the headway of the preceding car. If this headway is large the preceding car will not brake sharply and therefore much shorter distances between cars are acceptable.

Fig. 11 gives an example for this kind of behaviour. It shows empirical data for time-headways, i.e. the time intervals between the passing of consecutive cars recorded by a detector at a fixed position on the highway. The time-headway distribution in the free-flow regime shows a two peak structure. The peak at $`\mathrm{\Delta }t=1.8`$ sec corresponds to the recommended safe distance222This corresponds to the following rule of thumb taught in (German) driving schools: safe distance (in meters) = half of velocity (in km/h).. What is surprising is the existance of a second even larger peak at $`\mathrm{\Delta }t=0.8`$ sec. On a microscopic level these short time-headways correspond to platoons of some vehicles traveling very fast – their drivers are taking the risk of driving ”bumper-to-bumper” with a rather high speed. These platoons are the reason for the occurrence of high-flow states in free traffic. The corresponding states exhibit metastability, i.e. a perturbation of finite magnitude and duration is able to destroy such a high-flow state .

The NaSch model is not able to reproduce this behaviour. In  a extension of the NaSch model has been introduced which includes anticipation as a crucial ingrediant. The driving strategy of the drivers in this model comprises four aspects:
$`(i)`$ At large distances the cars move (apart from fluctuations) with their desired velocity $`v_{max}`$. $`(ii)`$ At intermediate distances drivers react to velocity changes of the next vehicle downstream, i.e. to “brake lights”. $`(iii)`$ At small distances the drivers adjust their velocity such that safe driving is possible. $`(iv)`$ The acceleration is delayed for standing vehicles and directly after braking events. This is achieved by introducing an interaction horizon and brake lights. The randomization parameter obeys a slow-to-start rule and in addition depends on the state of the brake light. Furthermore a finer discretization of the lattice allows for a more realistic modelling of the acceleration behaviour.

This model is also able to reproduce the typical behaviour found in so-called synchronized traffic. As first pointed out by Kerner and Rehborn  the congested phase can be divided into stop-and-go and synchronized traffic. In the latter, the mean velocity is considerably reduced, but all cars are moving and the flow is rather large.

A statistical observable which allows to identify synchronized traffic is the crosscorrelation

$$cc_{J\rho }=\frac{1}{\sqrt{\sigma (J)\sigma (\rho )}}\left[J(t)\rho (t+\tau )J(t)\rho (t+\tau )\right]$$
(4)
($`\sigma (\xi )`$ denotes the variance of $`\xi `$) between the density $`\rho `$ and the flow $`J`$ almost vanishes . This means that these variables are almost independent of each other and the data points cover randomly a plane in the fundamental diagram . In contrast, in the free-flow phase the flow is completely determined by the density and $`cc_{J\rho }(\tau )1`$.

Another interesting property of synchronized traffic is the existance of strong velocity correlations between the lanes in multilane highways. Lane changes are strongly surpressed and the traffic in the different lanes becomes “synchronized”. For a discussion of further features of synchronized flows and transitions between the different phases we refer to  (and refs. therein).

## 5 Effects of on- and off-ramps

There is large empirical evidence that the activity of on- and off-ramps can trigger transitions from free to congested flow . This view was supported by analysis of macroscopic models, where the introduction of on ramps leads to synchronized flow like behaviour .

In a very recent work , which I will summarize in the following subsection, simulations of microscopic models with different types of ramps have been discussed. In Sec. 5.2 experimental evidance  for the nonequilibrium phase transitions predicted by the ASEP (see Sec. 3.1) will be presented.

### 5.1 Ramps in periodic systems

First we discuss implementations of on- and off-ramps in single-lane systems with periodic boundary conditions where the number of cars is constant. The on- and off-ramps are implemented as connected parts of the lattice where the vehicles may enter or leave the system (see Fig. 12). The activity of the ramps is characterised by the number of entering (or leaving) cars per unit of time $`j_{in}`$ ($`j_{out}`$). A car is only added to the system if the removal of another car at the off-ramp is possible at the same time-step so that $`j_{in}=j_{out}`$. In  input and output are performed with a constant frequency. A stochastic in- and output of cars does not lead to a qualitatively different behaviour. The length chosen for the on- and off-ramps and the distance between them is motivated by the dimensions found on german highways: $`L_{ramp}=25`$ as length of the ramps in units of the lattice constant (usually identified with 7.5 meters). The position, i.e. the first cell of the on-ramp, is located at $`x_{on}=80`$ and that of the off-ramp at $`x_{off}=L80`$ where $`L`$ is the system size. Using periodic boundary conditions the distance of the on-ramp to the off-ramp is given by $`d_{ramp}=x_{on}x_{off}+L`$.

Two different procedures for adding cars to the lattice have been investigated in order to model inconsiderate (type A) and cautious behaviour (type B) of the drivers at the ramp. The only and essential difference in the implementation of the different types of on-ramps lies in the strategy of the cars changing from the acceleration lane to the driving lane. It has been found that all the procedures give qualitatively the same results.

In Fig. 13 the fundamental diagrams of the NaSch model for both types of on-ramps are shown for the input-rate $`j_{in}=j_{out}=\frac{1}{5}`$. For comparison the fundamental diagram without ramps is also shown as solid line.

It is clearly seen that a density regime $`\rho _{low}<\rho <\rho _{high}`$ exists where the flow $`J(\rho )`$ is independent of the density. This so-called plateau value of the flow is lower than the corresponding flow of the model without ramps. An increase of the input-rate $`j_{in}`$ leads to a decrease of the plateau value. If one compares the different types of input strategies it is evident that the plateau value of type B is lower than that of type A. This can be explained by the probability that the cell in front of the inserted car is already occupied. For type B (cautious) this probability is smaller than for type A (inconsiderate). Nevertheless there is no difference between the two types in a qualitative sense.

One can distinguish three different phases depending on the global density. In the high and low density phases the average flow $`J(\rho )`$ of the perturbated system takes the same value as in the system without ramps. For intermediate densities $`\rho _{low}<\rho <\rho _{high}`$ the flow is constant and limited by the capacity of the ramps.

This behaviour can be understood from the form of the density profiles (see Fig. 14). In the high and low density phase one only observes local deviations from a constant density profile. For intermediate densities,however, the system is separated into macroscopic high ($`\rho _{high}`$) and low ($`\rho _{low}`$) density regions. In the region of the ramps the density $`\rho _{ramp}`$ is additionally higher than in the high density region. Thus the ramps act like a blockage in the system that decreases the flow locally. Varying the global density within the phase-separated state the bulk densities in the high and low density region remain constant, only the length of the regions changes. Also the local density at the ramps ($`\rho _{ramp}`$) remains constant. Finally this leads to a constant flow in the segregated phase.

This phase separation in systems with ramps and the form of the fundamental diagram is very similar to the behaviour found in systems with a stationary defect . Typically, in a certain range on the lattice the slowdown parameter $`p`$ is increased to a value $`p_d>p`$ compared to that of the residual lattice sites . For intermediate densities $`\rho _{low}<\rho <\rho _{high}`$ the flow is constant in analogy to the model with ramps. In this regime $`J(\rho )`$ is limited by the capacity of the defect sites. There is no qualitative difference between the effect of on- and off-ramps and that of stationary defects. In both cases one observes plateau formation in the fundamental diagram as well as phase separation in the system. The only difference lies in the nature of the blockage dividing the system into two macroscopic regions. In the case of the ramps it is the local increase of the density that decreases the flow locally. In the model with defects the increased slowdown parameter leads to a local decrease of the flow.

In order to use the equivalence of ramps and defects in realistic simulations of traffic flow one has to determine the relation between the input-rates $`j_{in}`$ and the slowdown parameter $`p_d`$ at the defect sites. Fig. 15 shows the values of $`p_d`$ and $`j_{in}`$ leading to the same fundamental diagrams. It can be seen that a non-trivial relation between the two parameters exists.

This analogy is very important for efficient modelling of complex networks  with a multitude of ramps. Therefore it seems to be possible to achieve a considerable reduction of the computational complexity of network simulation using this kind of implementation of ramps.

### 5.2 Ramps in open systems

In  it has been shown that empirical data obtained in measurements on a German highway near Cologne can be interpreted using the results presented in Sec. 3.1. These experimental data exhibit boundary effects caused by the presence of an on-ramp. Fig. 16 shows a sketch of the relevant part of the highway. No experimental data are available for the density $`\rho ^{}`$ and flow $`j^{}=j(\rho ^{})`$ far upstream from the on-ramp and $`\rho ^+`$ and $`j^+=j(\rho ^+)`$ just before the on-ramp as well as the activity of the ramp. The only data come from a detector located upstream from the on-ramp which measures a traffic density $`\widehat{\rho }`$ and the corresponding flow $`\widehat{ȷ}`$.

Cars entering the motorway cause the mainstream of vehicles to slow down locally. Therefore, the vehicle density just before the on-ramp increases to $`\rho ^+>\rho ^{}`$. Then a shock, formed at the on-ramp, will propagate with mean velocity

$$v_{shock}=\frac{j(\rho ^+)j(\rho ^{})}{\rho ^+\rho ^{}},$$
(5)
obtained from mass conservation. $`v_{shock}`$ is the velocity of a ’domain-wall’ 333In nonequilibrium systems a domain wall is an object connecting two possible stationary states. between two stationary regions of densities $`\rho ^{},\rho ^+`$ The formation of a stable shock is usually a boundary-driven process caused by a ‘bottleneck’ on a road. Bottlenecks on a highway arise from a reduction in the number of lanes and from on-ramps where additional cars enter a road (see Sec. 5.1).

Depending on the sign of $`v_{shock}`$, two scenarios are now possible:
1) $`v_{shock}>0`$ (i.e. $`j^+>j^{}`$): In this case the shock propagates (on average) downstream towards the on-ramp. Only by fluctuations a brief upstream motion is possible. Therefore the detector will measure a traffic density $`\widehat{\rho }=\rho ^{}`$ and flow $`\widehat{ȷ}=j^{}`$.
2) $`v_{shock}<0`$ (i.e. $`j^+<j^{}`$): The shock wave starts propagating with the mean velocity $`v_{shock}`$ upstream, thus expanding the congested traffic region with density $`\rho ^+`$. The detector will now measure $`\widehat{\rho }=\rho ^+`$ and flow $`\widehat{ȷ}=j^+`$.

In order to discuss the transition between these two scenarios suppose that one starts with a situation where $`j^+>j^{}`$ is realized. If now the far-upstream-density $`\rho ^{}`$ increases it will reach a critical point $`\rho _{crit}`$ above which $`j^{}>j^+`$, i.e., the free flow upstream $`j^{}`$ prevails over the flow $`j^+`$ which the ‘bottleneck’, i.e. the on-ramp, is able to support. At this point shock wave velocity $`v_{shock}`$ will change sign (see (5)) and the shock starts traveling upstream. As a result, the stationary bulk density $`\widehat{\rho }`$ measured by the detector upstream from the on-ramp will change discontinuously from the critical value $`\rho _{crit}`$ to $`\rho ^+`$. This marks a nonequilibrium phase transition of first order with respect to order parameter $`\widehat{\rho }`$. The discontinuous change of $`\widehat{\rho }`$ leads also an abrupt reduction of the local velocity. Notice that the flow $`\widehat{ȷ}=j^+`$ through the on-ramp (then also measured by the detector) will stay independent of the free flow upstream from the congested region $`j^{}`$ as long as the condition $`j^{}>j^+`$ holds.

Empirically this phenomenon can be seen in the traffic data taken from measurements on the motorway A1 close to Cologne . Fig. 17 shows a typical time series of the one-minute velocity averages. One can clearly see the sharp drop of the velocity at about 8 a.m.

Also the measurements of the flow versus local density, i.e. the fundamental diagram (Fig. 18), support this interpretation. Two distinct branches can be distinguished. The increasing part corresponds to an almost linear rise of the flow with density in the free-flow regime. This part of the flow diagram is not affected by the presence of the on-ramp at all and one measures $`\widehat{ȷ}=j^{}`$ which is the actual upstream flow. The second branch are measurements taken during congested traffic hours. The transition from free flow to congested traffic is characterized by a discontinuous reduction of the local velocity. However, as predicted above the flow does not change significantly in the congested regime. In contrast, in local measurements large density fluctuations can be observed. Therefore in this regime the density does not take the constant value $`\rho ^+`$ as suggested by the argument given above, but varies from 20 veh/km/lane to 80 veh/km/lane (see Fig. 18). In contrast, spontaneously emerging and decaying jams would lead to the observation of a non-constant flow.

This shows that the empirical data obtained near an on-ramp can indeed be interpreted as a nonequilibrium phase transition. Furthermore it shows exactly the same behaviour as the transition from the low-density to the high-density phase in the ASEP (see Sec. 3.1).

The case of an off-ramp (or expansion of road space etc.) leads to a local decrease $`\rho ^+<\rho ^{}`$ of the density. Here the collective velocity $`v_c=\frac{j(\rho )}{\rho }`$ plays an important role. The theory than predicts a transition from the free-flow regime to the maximal flow (MF) phase . This transition is of second order because $`\widehat{\rho }`$ changes continuously across the phase transition point and has not yet been observed empirically. The existence of a maximal flow phase was not emphasized in the context of traffic flow up to now. At the same time, it is the most desirable phase, carrying the maximal possible throughput of vehicles $`j_{max}`$.

For practical purposes above observations may directly be used in order to operate a highway in the optimal regime. E.g. the flow near a lane reduction could be increased significantly if the traffic state at the entry would allow to attain the maximal possible flow of the bottleneck. This could be achieved by controlling the density far upstream, e.g. by closing temporally an on-ramp, such that the cars still enter the bottleneck with high velocity.

## 6 Summary and conclusion

In this short review only a few examples for the successfull application of CA models for the description of traffic flow have been presented. Despite its conceptual simplicity the NaSch model and its variants are able to reproduce many phenomena observed empirically. Although the same is true for many of the other approaches mentioned here only briefly in the introduction CA have the big advantage of being ideally suited for large-scale computer simulations. It is now possible to simulate large networks of highway  and city traffic  faster than real time. This offers new perspectives for the planning and design of transportation networks. Finally it is worthwhile pointing out that not only traffic science profits from the application of methods from physics. The example of the first empirical observation of a boundary-induced phase transition (see Sec. 5.2) shows that also physics can profit from applications to ”exotic fields”.

Acknowledgements: I like to thank all my collaborators, especially L. Santen and D. Chowdhury for their contribution to the extensive review . I also like to thank the organizers of the 36. Karpacz Winter School for creating a stimulating atmosphere.


format: text

text:
# Dissipative dynamics of vortex arrays in anisotropic traps

## I Introduction

Recent experimental evidence  for vortices in weakly interacting Bose-condensed gases stimulated a number of publications on structure of a rotating state of superfluids (see  for a review). Latest state of the art experimental technology makes possible creation and study of sufficiently large vortex arrays , as well as non-destructive measurements of time evolution in a system with a single vortex . Such possibilities allow, in principle, experimentally clarify the nature of interaction of vortices with thermal excitations. This issue plays an important role for understanding of vortex dynamic in numerous related systems, ranging from superconductors to superfluid cores of neutron stars (see  for a review), and still remains a controversial issue in the literature .

In our recent papers we investigated the dissipative dynamics of a single vortex  and a large vortex array  in a static (non-rotating) trap at finite temperatures. By this we assumed a very fast relaxation of the thermal component towards its static configuration due to, say, the anisotropy of the trap. Since actual experiments are performed in traps, which can be made nearly perfect, the rotation of the thermal component can persist for a very long time. Moreover, in the case of a large vortex array, the decay of vortices reaching the border of the condensate may be accompanied by a substantial heating of the thermal cloud, which, in turn, affects the strength of the friction forces acting on the remaining vortices. Thus, in order to understand an experiment in realistic conditions and use its results for a fine test of vortex dynamics theories, one has to have a solid description of dissipative dynamics including various temperature effects as well as the role of the trap anisotropy.

In this Letter we develop a simple, but yet a feature rich, model of dissipative dynamics of large vortex arrays in realistic anisotropic traps at finite temperatures. We analyze the relaxation kinetics of a system consisting of the rotating thermal component interacting with a condensate the vortices. We show that the trap anisotropy, or any other mechanism of angular momentum dissipation, is absolutely crucial for the onset of the dissipative dynamics. On the basis of our analysis we distinguish the two regimes of the evolution. The first one corresponds to the previously studied static case , i.e. to the situation when the trap is sufficiently anisotropic and the rotation of the thermal cloud quickly ceases. In this case the vortex array decays in a non-exponential fashion (a power law for the number of vortices at sufficiently long times). Remarkably, in this limit the lifetime of vortices does not depend on the parameters of anisotropy. In the opposite case, i.e. when the trap is nearly perfect, the friction of the vortices on the thermal cloud is stronger than the the friction of the thermal cloud on the trap asymmetry. This makes the vortices and the cloud stick to each other: the condensate rotates with the same angular velocity as the thermal particles and the angular velocity of the both components decreases exponentially with a time constant explicitly depending on the parameters of the trap asymmetry.

We discuss the heating of the thermal component due to the vortex array decay. The developed model, together with certain simplifications, allows us to obtain a simple analytical description of dissipative dynamics and compare the results with the experimental data . We believe that the presented model can be used for indirect study of vortex arrays dynamics based on the thermal component thermometry (proposed earlier in ).

## II The description of the model

We consider an evolving vortex array in an anisotropic trap characterized by the radial and the longitudinal frequencies $`\omega _\rho `$ and $`\omega _z`$. The equilibrium vortex configuration in a rotating trap (in the experiment the  this corresponds to the situation when the stirring beam is on) is considered in . Further on we assume that the initial angular velocity of the trap rotation $`\mathrm{\Omega }_0\mathrm{\Omega }_c`$, where $`\mathrm{\Omega }_c`$ is the critical velocity corresponding to the appearance of the first vortex. In this case the equilibrium number of vortices $`N_v1`$, and their spatial distribution is uniform.

As we will see, the relaxation of vortex arrays occurs on a time scale, exceeding the relaxation time in the thermal cloud. Hence, we may consider the evolution of the system in quasiequilibrium approximation. This means, that at any given moment of time, the thermal cloud is characterized by its equilibrium density profile, corresponding to a certain temperature $`T`$, and its angular velocity of the uniform macroscopic rotation $`\mathrm{\Omega }`$ (see Fig. 1). To derive the equations of motion for the whole system, consisting of the condensate with the vortices and the thermal cloud, we first write down the conservation laws. Assuming that the magnetic trap does not dissipate energy, we write down the conservation of the total energy

$$\frac{d}{dt}\left(\frac{I(T)\mathrm{\Omega }^2}{2}+E(T)+E_v\right)=0,$$
(1)
where $`E_v`$ is the energy of the condensate with the vortices (see below),

$$E(T)=3T\frac{\zeta (4)}{\zeta (3)}N(T),$$
(2)
is the internal energy of the non-condensed particles,

$$I(T)=mR_T^2N(T)\frac{\zeta (4)}{\zeta (3)},$$
(3)
is the momentum of inertia and $`R_T=(2T/m\omega _\rho ^2)^{1/2}`$ is the radial size of the thermal cloud. Here

$$N(T)=\frac{\zeta (3)T^3}{\stackrel{~}{\omega }^3}=\frac{T^3}{T_c^3}N$$
(4)
is the number of thermal particles, $`\stackrel{~}{\omega }=(\omega _z\omega _\rho ^2)^{1/3}`$, $`T_c=\stackrel{~}{\omega }(N/\zeta (3))^{1/3}`$ is the BEC transition temperature, and $`N`$ is the total number of particles in the trapped gas. Eqs. (2),(3) imply that $`T\mu `$, where $`\mu =n_0g`$ is the chemical potential, $`g=4\pi \mathrm{}^2a/m`$, $`a`$ is the two-body scattering length, and $`m`$ is the mass of the gas particles. At lower temperatures, $`T\mu `$, the size of the thermal cloud matches the size of the condensate and, hence, the temperature measurements become very difficult. Accordingly, in our quantitative discussions we confine ourself to the case $`T\mu `$.

In the Thomas-Fermi limit ($`\mu \mathrm{}\omega _{z,\rho }`$) the condensate is characterized by the longitudinal $`R_z=(2\mu /m\omega _\rho ^2)^{1/2}`$ and the transverse $`R_c=(2\mu /m\omega _z^2)^{1/2}`$ sizes, so that $`R_{c,z}l_c`$, where $`l_c=(\mathrm{}^2/m\mu )^{1/2}`$ is the correlation (healing) length in the condensate. The velocity of the superfluid flow, generated by a rotating vortex grid, imitates the flow of a normal liquid with the angular velocity $`\mathrm{\Omega }_v`$: $`𝐯_s=[\mathrm{\Omega }_v,𝐫]`$. The corresponding density of the vortex lines is given by

$$n_v=\frac{\mathrm{\Omega }_v}{\pi \kappa },$$
(5)
where $`\kappa =\mathrm{}/m`$ is the vortex circulation (vortices with a larger circulation are unstable ). The angular velocity $`\mathrm{\Omega }_v`$ must not be very large. Indeed, the maximum velocity of the superflow, which coincides with the velocity of the superfluid on the border of the condensate $`\mathrm{\Omega }_vR_c`$ must not reach the velocity of Bogolyubov sound $`c_s=\sqrt{\mu /m}`$, i.e.

$$\mathrm{\Omega }_v\frac{c_s}{R_c}\omega _\rho ,$$
(6)
since otherwise the superflow becomes unstable with respect to the spontaneous creation of excitations already at $`T=0`$ (see Landau arguments in ; the nature of the rotational friction in such a situation is an interesting question by itself and is not considered here, see  for a detailed discussion). Since the angular velocity $`\mathrm{\Omega }_v`$ is directly related with the density of the vortex array (see Eq. (5)), the condition (6) limits the maximum number of vortices in a stable vortex configuration confined in a harmonic trap of a given frequency:

$$N_v\pi R_c^2\frac{\mathrm{\Omega }_v}{\pi \kappa }\frac{\mu }{\mathrm{}\omega _\rho }1,$$
(7)
which is a large number in the Thomas-Fermi limit.

Since the density of the condensate energy is logarithmically singular in the vicinity of the vortex line, the energy of the vortex array, $`E_v`$, consists of two parts:

$$E_v=d^2r\rho _s\frac{v_s^2}{2}=E_0+I_v\frac{\mathrm{\Omega }_v^2}{2},$$
(8)
where $`E_0`$ is the part of the vortex energy independent on the velocities of the vortices, and amounts to $`N_v`$ times energy of a single vortex at rest, $`E_{v0}=\pi mn_0R_z\mathrm{log}(R_c/l_c)`$, see . The quantity $`I_v=mn_0R_z\pi R_c^4/2`$ can be called as the momentum of inertia for the condensate with the vortices (moment of inertia of the vortex array). The ratio of the first to the second term in the r.h.s. of Eq.(8) is $`\mathrm{log}(R_c/l_c)/N_v1`$. This shows that in a sufficiently large vortex array, $`N_v1`$, its rest energy is negligible and the whole internal energy of the condensate coincides with the energy of its macroscopic rotation.

The second equation comes from the angular momentum balance

$$\frac{d}{dt}(I(T)\mathrm{\Omega }+M_v)=\alpha \mathrm{\Omega },$$
(9)
where $`M_v`$ is the angular momentum of the vortex array and the phenomenological coefficient $`\alpha `$ describes the interaction of the thermal particles scattered by the trap anisotropy. In fact, the form of r.h.s. in Eq. (9) corresponds to the dissipation of the angular momentum by a friction force proportional to the velocity of the gas cloud. This relation between the friction force and the gas velocity is a very general one and holds for a friction force acting on an obstacle in a moving media both in the collisionless and hydrodynamic limits , at least for sufficiently small velocities. The choice of the friction coefficient $`\alpha `$ depends on the details of the anisotropy model and requires a microscopic calculation (see below).

Since in the reference frame rotating with the superfluid the trap asymmetry moves with a subsonic velocity (see the discussion leading to Eq.(6)), the condensate particles are not scattered by the trap and, therefore, no term proportional to the condensate angular velocity can appear in the r.h.s. of Eq.(9).

Unlike in the calculation of the condensate energy (8), the expression for the angular momentum of a vortex has no singularity and, hence, the angular momentum of the vortex array coincides with he angular momentum of the superfluid

$$M_v=d^2r\rho _srv_s=I_v\mathrm{\Omega }_v.$$
(10)
Eqs. (1),(9) are still insufficient, since they comprise a system of two equations for the three unknown functions $`\mathrm{\Omega }(t)`$, $`\mathrm{\Omega }_v(t)`$ and $`T(t)`$. The formulation of the model is completed by the addition of the equation of motion for the vortices. It follows from the Magnus law with the friction force originating from the interaction of the vortex lines with the thermal component. In fact it is the continuity equation for the vortex gas moving towards the border of the condensate (see Appendix A for the derivation, which closely follows the arguments given in ):

$$\frac{d\mathrm{\Omega }_v}{dt}+\gamma \mathrm{\Omega }_v(\mathrm{\Omega }_v\mathrm{\Omega })=0,$$
(11)
where $`\gamma =2\kappa \rho _sD/(\kappa ^2\rho ^2+D^2)`$, $`\rho _s=mn_0`$ and $`\rho `$ are the superfluid density and the total density of the gas, and $`D`$ is the mutual friction coefficient, characterizing the interaction of the vortex line with the excitations.

Eqs.(1),(9) and (11) constitute the complete set of equations describing the dissipative dynamics of vortex arrays interacting with thermal particles. The initial condition at $`t=0`$, when the rotation of the trap stops (i.e. when the stirring beam is off), are

$$\mathrm{\Omega }(0)=\mathrm{\Omega }_v(0)=\mathrm{\Omega }_0,T(0)=T_0>\mu .$$
(12)
The formulation of the model should be completed by ascribing a certain value to the friction coefficient $`\alpha `$. This can be done, using the results of the discussion of the rotational properties of trapped gas, presented in . Consider an ideal gas, splined up in an infinitely long trap, characterized by the two close transverse frequencies $`\omega _x`$ and $`\omega _y`$, so that $`\omega _{x,y}=\omega (1\pm ϵ)^{1/2}`$, where $`ϵ`$ is the parameter of anisotropy. Then, after the rotation of the trap stops, the gas continues rotation with the angular velocity, decreasing exponentially with the characteristic time constant given by the expression :

$$\lambda =\mathrm{Re}\frac{1}{4\tau }\left(1\sqrt{1\frac{ϵ^2}{ϵ_c^2}}\right),$$
(13)
where $`ϵ_c=1/4\omega \tau `$ and $`\tau `$ is the characteristic time between the collisions in the thermal cloud. The relaxation time $`\tau `$ can be estimated as $`\tau ^1n\sigma \sqrt{T/m}`$, where $`n`$ is the gas density and $`\sigma `$ is the collisional cross section. To generalize this result for the time $`\tau `$, obtained for a gas sample at $`T>T_c`$, for the case of a Bose condensed gas, we change the total density of the gas to the condensate density (see the argumentation in ).

In our model, when the condensate is absent, or the interaction between the vortices and the thermal component is negligible (this corresponds to $`\gamma D=0`$), the vortices and the thermal excitations evolve independently. From Eq.(1) for the temperature change in the thermal cloud we have

$$\frac{\mathrm{\Delta }T}{T}\frac{I\mathrm{\Omega }^2}{E(T)}\frac{\mathrm{\Omega }^2}{\omega _\rho ^2}1.$$
(14)
The latter inequality follows from Eq.(6) and guarantees that the temperature of the system does not change. In this approximation Eq. (9) reads

$$\frac{d\mathrm{\Omega }}{dt}=\frac{\alpha }{I}\mathrm{\Omega },$$
(15)
and gives the following simple decay law for the angular velocity of the thermal cloud

$$\mathrm{\Omega }(t)=\mathrm{\Omega }(0)\mathrm{exp}(t/\tau _A),$$
(16)
where $`\tau _A=I/\alpha `$ is the relaxation time related to the friction of the moving gas and the trap anisotropy. Introducing the requirement, that the presented model recovers the result of Guery-Odelin , we set $`\tau _A^1=\lambda `$ and find, that $`\alpha =I\lambda `$.

The appearance of imaginary part in the relaxation rate (13) for $`\omega \tau ϵ^1`$ implies an oscillatory behavior in the onset of equilibrium of a free gas cloud (see  for the discussion). As we shall see, the evolution of the combined system of vortices and thermal excitations explicitly depends on the relaxation rate (13) only in the so called “weak anisotropy” limit, when both the vortex array and the non-condensed gas rotate with the same angular velocity. Since, according to Eq. (7), the number of vortices in the vortex array is proportional to the angular velocity $`\mathrm{\Omega }_v`$, such oscillations would imply either periodic oscillations of the vortex density (if the oscillations of $`\mathrm{\Omega }_v`$ are sufficiently small), or a reversible destruction and nucleation of vortices close to the border of the condensate. Well below $`T_c`$, the nucleation of vortices  is slow compared with kinetic times (13), and, thus, the number of vortices, as well as the angular velocity of the condensate rotation, can only decrease in the course of dissipative dynamics. Although the detailed description requires a more thorough discussion, we neglect the oscillations and keep $`\alpha `$ purely real in the remaining discussion. Collective oscillations of vortices and the thermal cloud are reviewed in .

## III Analysis of the dissipative dynamics.

As discussed in the previous section, Eqs. (1),(9) and (11) form a close set of equations describing the evolution of a large ($`N_v1`$) vortex array emerged in a rotating thermal cloud. Although the equations look simple, their analytical solution requires certain simplifications.

First of all we qualitatively analyze the heating of the thermal cloud in the course of dissipation of rotational energy of the gas. According to Eq. (14) the rotational energy of the thermal cloud is always less than its internal energy and therefore makes no contribution to the heating of the system. Thus we can neglect it in the energy conservation law and express the temperature of the thermal cloud from Eq. (1) as a function of angular velocity of the condensate

$$T=T_0\left[1+A\left(1\frac{\mathrm{\Omega }_v^2}{\mathrm{\Omega }_0^2}\right)\right]^{1/4}.$$
(17)
where $`A`$ is the ratio of the rotational energy of the condensate to the internal energy of the thermal cloud. It can be found with the help of Eqs. (2) and (10)

$$A=\frac{I_v\mathrm{\Omega }_0^2}{2E}\left(\frac{\mu }{T_0}\right)^4\frac{\mathrm{\Omega }_0^2}{\omega _\rho ^2}\frac{1}{(n_0a^3)^{1/2}}$$
(18)
and can well be large. The dimensionless parameter $`A`$ indicates how important are the heating effects in the course of the dissipative dynamics. On the other hand, it can be shown that the temperature increase $`T_fT_0T_c`$. This allows us to neglect the temperature dependence of $`\rho _s`$ and $`I_v`$.

Remarkably, in Eq. (9) we can neglect also the time dependence of the momentum of inertia of a thermal cloud. Indeed, the term containing the time derivaty of the temperature, $`\dot{I}\mathrm{\Omega }`$, turns to be always small compared to $`I_v\dot{\mathrm{\Omega }}_v`$: $`\dot{I}\mathrm{\Omega }/I\dot{\mathrm{\Omega }}_v(\mathrm{\Omega }_v\mathrm{\Omega }/\mathrm{\Omega }_0^2)(\omega _\rho /\mathrm{\Omega }_0)^21`$. Using this fact in Eqs. (9) and (11) we obtain

$$\frac{d}{dt}(\mathrm{\Omega }\mathrm{\Omega }_v)=\frac{\alpha }{I}\left[\eta \frac{\mathrm{\Omega }}{\mathrm{\Omega }_0}(\mathrm{\Omega }\mathrm{\Omega }_v)+\mathrm{\Omega }\right]$$
(19)
where the dimensionless parameter

$`\eta ={\displaystyle \frac{(I+I_v)\gamma \mathrm{\Omega }_0}{\alpha }}={\displaystyle \frac{\tau _R}{\tau _v}},`$

characterizes the relative strength of the mutual interaction between the vortices and the thermal particles and the interaction of the thermal cloud and the trap asymmetry. Here $`\tau _v^1=\gamma \mathrm{\Omega }_0`$ and $`\tau _R^1=\alpha /(I+I_v)`$ are the characteristic relaxation times, determined respectively by the interaction of thermal component with the vortices and with the trap anisotropy.

The parameter $`\eta `$ distinguishes two limiting regimes of relaxation. First we turn to the static trap limit, when $`\tau _R\tau _v`$ ($`\eta 1`$). In this case at the initial stage the friction of the thermal cloud on the trap anisotropy is much stronger than the friction on the condensate. Therefore, at this stage the relaxation of the thermal cloud should proceed in the same way as in the absence of vortices, i.e. the angular velocity $`\mathrm{\Omega }`$ exponentially drops on the time scale $`\tau _A\tau _R\tau _v`$, while the angular velocity of a condensate remains almost unchanged, $`\dot{\mathrm{\Omega }}_v0`$. Neglecting the first term in the brackets of r.h.s of Eq. (19) and using Eq. (9) we reduce Eq.(19) to Eq. (15). This qualitative picture remains valid until the thermal clod nearly stops and its angular velocity reaches the very small value $`\eta \mathrm{\Omega }_v^2/\mathrm{\Omega }_0\mathrm{\Omega }_0`$. At this point the trap becomes effectively static and, setting $`\mathrm{\Omega }0`$ in Eq. (11), we find the equation of motion for the condensate angular velocity

$$\frac{d\mathrm{\Omega }_v}{dt}=\gamma \mathrm{\Omega }_v^2.$$
(20)
If the time dependence of the coefficient $`\gamma `$ can be neglected, that either for $`A1`$, or for very small times when $`\mathrm{\Omega }_v\mathrm{\Omega }_0`$, the solution of this equation is given by the expression

$$\mathrm{\Omega }_v(t)=\frac{\mathrm{\Omega }_0}{1+t/\tau _v},$$
(21)
which recovers the result previously obtained in .

To analyze the solutions of Eq. (20) in the strong heating limit, $`A1`$, we note that $`\gamma T`$ and rewrite Eq. (20) with the help of Eq. (17) and the dimensionless units for the condensate angular velocity $`\stackrel{~}{\mathrm{\Omega }}_v=\mathrm{\Omega }_v/\mathrm{\Omega }_0`$ and the time $`\stackrel{~}{t}=t/\tau _R(T_0)`$:

$$\frac{d\stackrel{~}{\mathrm{\Omega }}_v}{d\stackrel{~}{t}}=(1+A(1\stackrel{~}{\mathrm{\Omega }}_v^2))^{1/4}\stackrel{~}{\mathrm{\Omega }}_v^2.$$
(22)
Integrating it we find the exact implicit solution $`\mathrm{\Omega }_v(t)`$:

$`{\displaystyle \frac{(1+AA\stackrel{~}{\mathrm{\Omega }}_v^2)^{3/4}}{(1+A)\stackrel{~}{\mathrm{\Omega }}_v}}+{\displaystyle \frac{A\stackrel{~}{\mathrm{\Omega }}_v}{2(1+A)^{5/4}}}F_{\frac{1}{2};\frac{1}{4};\frac{3}{2}}\left({\displaystyle \frac{A\stackrel{~}{\mathrm{\Omega }}_v^2}{1+A}}\right)=\stackrel{~}{t}+C,`$

where $`F_{a,b,c}(x)`$ is the hypergeometric function and the integration constant $`C`$ is to be determined from the initial condition a $`\stackrel{~}{t}=0`$. In the considered limit of large $`A`$ we find, that

$`C={\displaystyle \frac{\sqrt{\pi }\mathrm{\Gamma }(3/4)}{4\mathrm{\Gamma }(5/4)A^{1/4}}}.`$

Then, for small $`\stackrel{~}{\mathrm{\Omega }}_v`$, which corresponds to large times, when the rotation of the condensate is almost stopped, the hypergeometric function can be expanded in series of its small argument and in the lowest order in $`1/A`$ we obtain the asymptotic expression for the condensate angular velocity

$`\mathrm{\Omega }_v={\displaystyle \frac{\mathrm{\Omega }_0}{0.4+t/\tau _v(T_f)}},`$

which is valid as soon as $`t\tau _v(T_f)`$. This means, that the quantity $`\tau _v(T_f)`$ plays the role of the vortex array life time and is mainly determined by the friction coefficients at the final temperature. This fact is fairly intuitive, since the strength of the interaction of the vortices with the thermal cloud quickly increases with the rise of the temperature and thus most of the relaxation of the system occurs when its temperature is close to its final value. Remarkably, the dynamics of the vortex array in a static trap limit, that is when anisotropy is strong, turns out to be independent from the parameters of the anisotropy.

In the opposite, rotating trap limit, when the anisotropy is weak and $`\tau _R\tau _v`$ ($`\eta 1`$), the friction between the vortices and the thermal cloud is stronger than the friction of the thermal excitation on the trap asymmetry. As the result, the both components stick to each other due to their mutual interaction and dissipate the angular momentum together, as a single rigid body. Neglecting the last term in the brackets of r.h.s. of Eq. (19) for the initial conditions given by (12) we find the solution $`\mathrm{\Omega }=\mathrm{\Omega }_v`$. Using this relation between the angular velocities in Eq.(9) we arrive at

$$\frac{d\mathrm{\Omega }_v}{dt}=\frac{\alpha }{I+I_v}\mathrm{\Omega }_v.$$
(23)
Again, when the heating of the gas sample is small, the coefficient $`\alpha `$ is time independent, and the solution of this equation is simple:

$$\mathrm{\Omega }_v(t)=\mathrm{\Omega }_0\mathrm{exp}(t/\tau _R),$$
(24)
i.e. in this, “rotating trap” approximation, the power law (21) changes into the simple exponential decrease of the angular velocities, with the same time constant for both the vortex array and the thermal cloud.

To discuss the solutions of Eq.(23) in the strong heating limit $`A1`$ we have to make certain assumptions about the trap anisotropy. For example, according to Eq.(13), when the trap asymmetry is small, $`4ϵ\omega \tau 1`$, we have $`\alpha T^2`$. Then Eq.(23) has the solution

$$\mathrm{\Omega }_v=\mathrm{\Omega }_0\frac{2\mathrm{exp}(t/\tau _R(T_f))}{1+\mathrm{exp}(2t/\tau _R(T_f))}.$$
(25)
This means that again, the relaxation time at the finite temperature $`\tau _R(T_f)`$ plays the role of the vortex array life time. At larger times, $`t\tau _v(T_f)`$, the angular velocity is small and the solution has an approximate asymptotic form

$`\mathrm{\Omega }_v=2\mathrm{\Omega }_0\mathrm{exp}(t/\tau _R(T_f)).`$

Remarkably, in logarithmical approximation the solution (25) holds also for a situation, when the trap asymmetry is strong, $`4ϵ\omega \tau 1`$, and $`\alpha T^6`$ (see Eq.(13)).

Using the presented above analysis we now can draw a qualitative picture of dissipative dynamics. To do this, we introduce the dimensionless ratio of the momenta of inertia of the thermal cloud and the vortex array: $`\zeta =I/I_v(T_0/\mu )^4(n_0a^3)^{1/2}`$. Since in a weakly interacting gas the gaseous parameter is small: $`(n_0a^3)^{1/2}1`$, at $`T\mu `$ the condensate is heavier than the thermal cloud, and both the rotational energy and the angular momentum are dominated by the condensate contributions. Thus, in the rotating trap limit ($`\eta 1`$) the relaxation time $`\tau _R`$ is mainly determined by $`I_v`$ and is much larger than the spinning down time of a free thermal cloud $`\tau _A`$ . In the static trap limit the role of the parameter $`\zeta `$ is more subtle. At sufficiently large times, $`t\tau _A`$, we can neglect the term $`I\dot{\mathrm{\Omega }}`$ in Eq. (9), which thus transforms into $`I_v\dot{\mathrm{\Omega }}=\alpha \mathrm{\Omega }`$. Comparing this equation with Eq. (20) we find that $`\mathrm{\Omega }=\eta \mathrm{\Omega }_v^2/\mathrm{\Omega }_0`$, i.e. the thermal cloud angular velocity decays as a power law and closely follows the relaxation of the vortex array. In the same time the heating parameter $`A`$ is related to the parameter $`\zeta `$: strong heating of the sample is only possible when $`\zeta 1`$. As the temperature grows, the number of the condensed particles decreases and at some point the condensate may become so small, that $`\zeta 1`$. One can see, that the crossover between the two regimes occurs at the temperature below $`T_c`$, so the condensate does not disappear. Later on, the heating of the system stops and the system relaxes according to Eqs. (21) and (24) for $`\eta 1`$ and $`\eta 1`$ respectively.

The presented analysis shows that the anisotropy, or any other mechanism of angular momentum dissipation, is absolutely crucial for the onset of dissipative dynamics. It is seen already from the fact, that when $`\tau _R\mathrm{}`$, Eq.(24) gives $`\mathrm{\Omega }_v(t)=\mathrm{\Omega }(t)=\mathrm{\Omega }_0`$, i.e. the rotation of the vortex grid persists forever.

Although the presented analysis is sufficient for a qualitative understanding of the dissipative dynamics in the whole parameters range, the quantitative description requires a numerical solution of the initial Eqs. (1),(9) and (11). For the sake of completeness we provide the results of such calculation based on the experimental parameters of . We use the following values for the model parameters: the condensate density $`n_0=2\times 10^{14}cm^3`$, the trap frequencies $`\omega _\rho =2\pi \times 169s^1`$ and $`\omega _z=2\pi \times 11.7s^1`$, and the temperature $`T_0\mu =80nK`$. Then, the condensate size is $`R_c=4\mu m`$ and, at the reported frequency of the stirring $`\mathrm{\Omega }_0=2\pi \times 135s^1`$, the vortex array consists of $`N_v10`$ vortices.

We provide the results of the numerical calculations for the trap anisotropy parameter $`ϵ=0.06`$ (see Figs. 2 and 3). Our model system turns out to be well in the rotating trap limit ($`\eta 10`$) and, in the full agreement with our qualitative analysis, $`\mathrm{\Omega }\mathrm{\Omega }_v`$ and either of the angular velocities relaxes exponentially (see the solid and the dashed lines on Fig. 2). The solid line on Fig. 3 illustrates the heating of the system in the course of the dissipative dynamics.

The comparison of the exact numerical solution with our approximate analytical results is a bit difficult, since this value of the anisotropy parameter is close to $`ϵ_c`$ (see Eq.(13)). The dashed and the dotted dashed lines on Figs. 2 and 3 represents the analytical results based on Eq. (25) with $`\tau _v`$ calculated in the strong $`ϵϵ_c`$ and the weak $`ϵϵ_c`$ anisotropy limits in Eq.(13) respectively. We note, that all of the presented results are in a good agreement with the experimentally found vortex array life time $`1s`$.

## IV Conclusion

In this Letter we formulate a realistic model, describing the evolution of a large vortex array in a rotating thermal cloud interacting with the trap anisotropy. The model allows both the simple analytical description of the most important limiting cases as well as a very straightforward numerical analysis.

First we show that the angular momentum dissipation mechanism, such as the trap anisotropy is absolutely crucial for the onset of dissipative dynamics of vortices in trapped Bose-condensed gases. We demonstrate, that if the total angular momentum of the system is conserved, the rotation of the vortex array persists forever and, consequently, the number of the vortices in the condensate remains constant.

Since the flow of the condensate is a superflow, the condensate itself can not interact with the trap anisotropy. On the contrary, the thermal excitations are scattered by the trap and can dissipate the angular momentum. Thus, it is just the mutual friction between the vortices and the gas particles that allows the condensate to transfer its angular momentum to the excitations and, subsequently, to the trap. Depending on the relative strength of interaction of the thermal gas with vortices and with the trap anisotropy we distinguish the two regimes of the dissipative dynamics.

First we identify the limit, when the trap anisotropy is very strong and the rotation of the thermal cloud quickly stops. Then the vortex array interacts with the thermal cloud at rest and subsequently looses the vortices due to the mutual friction between the vortices and the excitations (static trap limit). This case has been previously studied in  for a single vortex and in  for a large vortex array.

In the opposite limit, when the trap is nearly perfect, the thermal cloud stops slowly. In this approximation the mutual friction between the vortices and the thermal cloud makes the two components stick to each other, forming a single rotating object with the momentum of inertia equal to the sum of the momenta of inertia of the condensate and the thermal cloud. This composite body is usually much heavier than the thermal cloud and hence its relaxation proceeds on a time scale much larger than that in the Guery-Odelin scenario , where the calculation was done above $`T_c`$ and only the thermal cloud was taken into account.

Another feature of a large vortex array evolution is the heating of the thermal component caused by the energy released by the vortices escaping the condensate. We show, that though the heating itself can be very substantial, it does not change the qualitative picture of the dissipative dynamics. Since the friction coefficients $`\alpha `$ and $`\gamma `$ both grow very quickly with the temperature increase, most of the relaxation processes occurs when the temperature is already high, i.e. when the temperature is close to $`T_f`$. Hence, for estimation purposes, one can use a much simpler calculation with $`T=\mathrm{Const}`$ with all the model parameters taken at the final temperature.

## Acknowledgement.

We acknowledge fruitful discussions with G.V. Shlyapnikov, J. Dalibard and I.E. Shvarchuck. The work was supported by Austrian Science Foundation, by INTAS, and Russian Foundation for Basic Research (grant 99-02-18024).

## A Dynamics of vortex arrays.

To derive the equation of motion for the vortices inside a rotating thermal cloud, we closely follow the derivation presented in . The configuration of the superfluid is described by the superfluid velocity $`𝐯_s`$ and by the two continuous functions: the average vortex density $`n_v(t,𝐫)`$ and velocity $`𝐯_v(t,𝐫)`$. Since the number of vortex lines is a locally conserved quantity, the introduced functions satisfy the continuity equation:

$$\frac{n_v}{t}+\frac{(𝐯_vn_v)}{𝐫}=0.$$
(A1)
The velocity of the superfluid satisfies Maxwell-like equations

$$\mathrm{rot}𝐯_s=2\pi \kappa n_v\widehat{z},\mathrm{div}𝐯_s=0,$$
(A2)
where $`\widehat{z}`$ is a unit vector along the axis of the cylinder . These equations are solved with the boundary condition that the normal component of $`𝐯_s`$ vanishes on border of the condensate. The vortex velocity is related to the superfluid velocity through the Magnus law :

$$\rho _s[𝐯_v𝐯_s,\kappa ]=𝐅,$$
(A3)
where $`\rho _s`$ is the superfluid (mass) density and $`𝐅`$ is the friction force acting on the vortex line (per unit length).

At finite temperatures, the force $`𝐅`$ originates from the scattering of thermal excitations from the moving vortices and consists of two terms:

$$𝐅=D(𝐯_v𝐯_n)+D^{}[\kappa \widehat{z},(𝐯_v𝐯_n)],$$
(A4)
where $`𝐯_n`$ is the velocity of the thermal gas at the position of the vortex line. Here $`D`$ and $`D^{}`$ are the longitudinal and transverse friction coefficients respectively. Both of them are temperature dependent and $`D^{}=\rho _n`$, where $`\rho _n`$ is the density of the normal component . The friction coefficient $`D`$ is also temperature dependent and ranges from $`D\kappa m^{5/2}T^5/\mu ^{7/2}\mathrm{}^3`$ for very small temperatures $`T\mu `$  to $`D0.1\kappa m^{5/2}T\mu ^{1/2}/\mathrm{}^3`$  for $`\mu TT_c`$.

As discussed in , the vortex density is spatially uniform: $`n_v(t,𝐫)=n_v`$ and the superfluid velocity inside the vortex array can be found from Eq.(A2) by using the Stokes theorem. It has no radial component (i.e. it satisfies the necessary boundary condition) and is given by

$$𝐯_s(𝐫)=\pi [\kappa ,𝐫]n_v.$$
(A5)
The solution of Eqs.(A3) and (A4), with $`𝐯_s`$ from Eq.(A5) and $`𝐯_n=[\mathrm{\Omega },𝐫]`$, can be represented in the form: $`𝐯_v=v^{(r)}\widehat{𝐫}+v^{(\varphi )}[\widehat{\kappa },\widehat{𝐫}]`$, where

$$v_v^{(r)}=(v_sv_n)\frac{\kappa \rho _sD}{\kappa ^2\rho ^2+D^2}$$
(A6)
and

$$v_v^{(\varphi )}=\frac{D^2v_n+\rho \kappa ^2(v_s\rho _s+D^{}v_n)}{\kappa ^2\rho ^2+D^2}$$
(A7)
are the radial and tangential components of the velocity field, and $`\rho =\rho _s+\rho _n`$ is the total density. Eqs. (A6) and (A7) show, that if the condensate rotates together with the normal component, $`v_s=v_n`$, then $`v_v^{(\varphi )}=v_n`$ and the radial vortex velocity is absent, i.e. the vortex and the thermal particles move together and are in the equilibrium in the reference frame defined by the rotating cloud.

Using Eqs.(A5) and (A6), we rewrite Eq.(A1) in the form

$`{\displaystyle \frac{dn_v(t)}{dt}}+\gamma n_v(t)(\pi \kappa n_v(t)\mathrm{\Omega })=0,`$

where $`\gamma =2\kappa \rho _sD/(\kappa ^2\rho ^2+D^2)2D/\kappa \rho `$. Finally, expressing the vortex density through the angular velocity of the condensate rotation (see Eq.(5)), we obtain the desired equation of motion for the vortices (11).


text:
# Geometric Phases and Multiple Degeneracies in Harmonic Resonators

\[

## Abstract

In a recent experiment Lauber et al. have deformed cyclically a microwave resonator and have measured the adiabatic normal-mode wavefunctions for each shape along the path of deformation. The nontrivial observed cyclic phases around a 3-fold degeneracy were accounted for by Manolopoulos and Child within an approximate theory. However, open-path geometrical phases disagree with experiment. By solving exactly the problem, we find unsuspected extra degeneracies around the multiple one that account for the measured phase changes throughout the path. It turns out that proliferation of additional degeneracies around a multiple one is a common feature of quantum mechanics.

\]

Geometric phases may show up whenever the system under study depends on several parameters and is transported adiabatically around a closed path in parameters space. Such phases have been predicted and observed in many different systems . In particular, the cyclic phases observed by Lauber et al. for a deformed microwave resonator have been recently accounted for within an approximate theory . However, the closed-path phases are not the only interesting observable quantities: for an open-path evolution, the phase of the projection of a parallel-transported  eigenstate on the initial eigenstate is equally well defined  whenever these two states are not orthogonal. The experiment of Ref.  measures also such phase relations for a number of intermediate points along the path. These phases can be compared with the theoretical outcome of the calculation of Ref. . We observe that the two sets of values disagree at half loop. The geometrical nature of these phases implies that the approximation of Ref.  is missing the topological structure of the simple system at hand.

To find what is missing, we push the perturbative expansion of Ref.  to higher order and demonstrate the presence of unsuspected extra degeneracies which account for the phase discrepancies. We show that similar “satellite” degeneracies close to a multiple degeneracy appear as a systematic, predictable feature of the Laplace operator in distorted domains: they can therefore arise in a broad range of undulatory phenomena, including acoustics, optics and quantum mechanics.

In Ref.  Lauber et al. follow the adiabatic evolution of three eigenstates of a microwave resonator driven along a loop in the two-dimensional space of deformations around the rectangular geometry. We indicate with $`𝐫=(r\mathrm{cos}\theta ,r\mathrm{sin}\theta )`$ the external parameter representing the displacement of the right upper corner of the cavity from the rectangular position $`(a,b)=(\sqrt{3},1)`$ at which the eigenstates $`|\psi _1=\{2,4\}`$, $`|\psi _2=\{5,3\}`$ and $`|\psi _3=\{7,1\}`$ become degenerate \[with $`\{n_x,n_y\}`$ we indicate the wavefunctions $`2\mathrm{sin}(n_xx/a)\mathrm{sin}(n_yy/b)/\sqrt{ab}`$ in the rectangle\]. We define the angle $`\theta `$ so that it goes from $`0`$ to $`2\pi `$ following the path of distortions used in the experiment. In Fig. 1 we report a sequence of adiabatic eigenfunctions $`|\psi _j(\theta )`$ from the original pictures of Ref. .

The cyclic phase factors \[$`\gamma _1(2\pi )=1`$, $`\gamma _2(2\pi )=+1`$ and $`\gamma _3(2\pi )=1`$\] are easily identified in Fig. 1 from the recurrence of the pattern and the sign changes at $`\theta =0`$ and $`2\pi `$. Similarly, the sign change of the central state $`|\psi _2(\theta )`$ at $`\theta =\pi `$ tells us that its Berry-Pancharatnam open-path phase factor is $`\gamma _2(\pi )=1`$.

Let us compare now these results with those of the linear approximation. The eigenstates of the 2-dimensional Laplacian operator represent the normal modes of this cavity. The secular problem in the deformed domain can be conveniently reduced to the secular problem of the following differential operator

$$H(𝐫)=\underset{i,j}{}\frac{}{u_i}M_{ij}(𝐮,𝐫)\frac{}{u_j}+D(𝐮,𝐫)$$
(1)
with vanishing boundary conditions in an undistorted rectangular domain for the transformed variables $`𝐮`$. This form, given explicitly in Ref. , is particularly suitable to set up a perturbation theory in $`r`$

$$H(𝐫)=rH^{(1)}(\theta )+r^2H^{(2)}(\theta )+\mathrm{},$$
(2)
since all the parametric dependence is now in the functions $`M_{ij}(𝐮,𝐫)`$ and $`D(𝐮,𝐫)`$. At first order it suffices to consider the subspace span by the three eigenstates of Lauber et al.’s experiment, and compute the $`3\times 3`$ matrices $`_{ij}^{(1)}(\theta )=\psi _i|H^{(1)}(\theta )|\psi _j`$. The form of this matrix, given in Ref. , is

$$^{(1)}(\theta )=\mathrm{cos}\theta F+\mathrm{sin}\theta F^{},$$
(3)
where $`F`$ and $`F^{}`$ are real symmetric numerical matrices. In a sufficiently small neighborhood of the 3-fold degeneracy, $`^{(1)}(\theta )`$ accounts for the leading contribution to the energy shifts.

Manolopoulos and Child  calculate the (real) transformation matrix $`U_{ij}^{(1)}(\theta )=\psi _i(0)|\psi _j(\theta )`$, which diagonalizes $`^{(1)}(\theta )`$ parallel-transporting the eigenstates. In particular, $`U_{ii}^{(1)}(2\pi )`$ is the cyclic geometric phase of state $`i`$. Each calculated phase agrees with that observed experimentally. More generally, for any $`\theta `$, $`U_{ii}^{(1)}(\theta )/|U_{ii}^{(1)}(\theta )|`$ is the observable (open-path) Berry-Pancharatnam phase factor $`\gamma _i`$. These phase factors are therefore represented by the signs of the curves of Fig. 2b of Ref.. In Fig. 2 we report $`U_{22}^{(1)}(\theta )`$ for the intermediate state $`|\psi _2`$. Near $`\theta =\pi `$, the disaccord with experiment is particularly evident: linear theory has $`U_{22}^{(1)}(\theta )+1`$, while experimentally $`\gamma _2(\pi )=1`$.

Around $`\theta =\pi `$, $`|U_{11}|`$ and $`|U_{33}|`$ are very small ($`1`$): the interesting geometric phase information is contained in the off-diagonal part of the matrix $`U(\pi )`$ . In particular, in this case $`\gamma _{13}(\theta )U_{13}(\theta )U_{31}(\theta )/|U_{13}(\theta )U_{31}(\theta )|`$ is the only observable phase associated to state 1 and 3, since $`|U_{13}(\theta )U_{31}(\theta )|1`$. From Fig. 1 one can read $`U_{13}(\theta )+1`$, and $`U_{31}(\theta )+1`$. Thus the observed $`\gamma _{13}=+1`$. We have verified that the linear theory of Ref.  gives instead $`\gamma _{13}=1`$.

A similar reasoning can be applied to the path $`\mathrm{\Gamma }_2`$ leading back from $`\theta =\pi `$ to $`\theta =2\pi `$. The observed and calculated phase changes for this second path are the same as those for path $`\mathrm{\Gamma }_1`$ going from $`\theta =0`$ to $`\theta =\pi `$: they are summarized in Table I.

The failure of the linear approximation in predicting the correct geometrical phase factors suggests that it misses the topological structure of degeneracies in parameters space. In this Letter, we investigate this issue by considering both the next order in the expansion (2) and an exact numerical solution.

The origin of the problem can be understood from inspection of Fig. 3. Indeed, close to $`\theta =\pi /2`$ and $`3\pi /2`$, the eigenenergies of $`^{(1)}`$ get pairwise very near (the minimum splittings being $`1\%`$ of those in the $`\theta =0`$ direction). In this region, the linear approximation breaks down for $`r`$ much smaller than in the (say) $`\theta =0`$ direction. Second-order terms can close the gaps inducing additional degeneracies in these directions, remaining dominant with respect to third- and higher-order ones. It is the smallness of the first-order splitting coefficient in this direction $`\theta =\pi /2`$ and $`3\pi /2`$ that makes this possible.

We compute therefore the second-order Hamiltonian correction:

$$_{ij}^{(2)}=\psi _i|H^{(2)}|\psi _j+\underset{k1,2,3}{}\frac{\psi _i|H^{(1)}|\psi _k\psi _k|H^{(1)}|\psi _j}{E_iE_k}$$
(4)
To calculate explicitly $`^{(2)}`$, we cutoff the sum of Eq. (4) for $`|E_iE_k|<E_{\mathrm{max}}`$. Convergence is achieved including $`10^2`$ states. We have thus obtained the real symmetric numerical matrices $`G`$, $`G^{}`$, and $`G^{\prime \prime }`$ in $`^{(2)}(\theta )=G+G^{}\mathrm{cos}2\theta +G^{\prime \prime }\mathrm{sin}2\theta `$. This second-order contribution adds to the linear Hamiltonian (3) to constitute the secular problem for the matrix

$$r^{(1)}(\theta )+r^2^{(2)}(\theta ).$$
(5)
We thus obtain the eigenvalues and eigenvectors of (5) as analytical functions of $`r`$ and $`\theta `$. Like in the linear case, the transformation matrix $`U_{ij}^{(2)}(r,\theta )`$ contains all the phase information of the eigenvectors. Now, at second order, all the geometrical phases, both for the closed loop and for the open-path, agree with experiment  (cf. Tab. I and Fig. 2).

To elucidate the differences between first and second order, we show in Fig. 4 the region of validity of the linear approximation. The endpoints W and Z of the open path $`\mathrm{\Gamma }_1`$ considered above lie both well within this region. In addition, there exists a path $`\mathrm{\Gamma }_3`$ connecting W to Z, lying all inside the linear region on the same side of the $`3`$-fold degeneracy. For path $`\mathrm{\Gamma }_3`$ the first-order theory is sufficient to predict the correct phases. Consider now the loop $`\mathrm{\Gamma }_1\mathrm{\Gamma }_3`$ ($`\mathrm{\Gamma }_3`$ is followed backwards from Z to W). If no degeneracies were present inside, the cyclic phase factors $`\gamma _2^{\mathrm{\Gamma }_1\mathrm{\Gamma }_3}`$ would be $`+1`$ for all states. Then the open-path phases $`\gamma _2^{\mathrm{\Gamma }_1}`$ and $`\gamma _2^{\mathrm{\Gamma }_3}`$ would have to coincide. The same applies to the off-diagonal $`\gamma _{13}`$ . Different phases in first and second order calculations indicate additional conical intersections of the second-order Hamiltonian (5) inside the loop.

Indeed, we find that, at the points A ($`r_𝐀0.0126`$, $`\theta _𝐀1.64`$) and B ($`r_𝐁0.0120`$, $`\theta _𝐁4.64`$) of Fig. 4, the two lowest eigenvalues are degenerate. Knowledge of the phases in the linear region around the threefold degeneracy at the origin  together with the presence of the twofold conical intersections at A and B accounts for all (open-path and cyclic) phases. In particular, on the complete loop $`\mathrm{\Gamma }_1+\mathrm{\Gamma }_2`$, the two sign changes that $`|\psi _1`$ and $`|\psi _2`$ acquire encircling A and B cancel out: as a fortuitous coincidence, the observed cyclic Berry phases end up being the same as these obtained in the linear theory.

As well known for a generic real Hamiltonian it is necessary to vary two parameters in order to make a degeneracy occur accidentally, i.e. not on account of symmetry. In other terms, the codimension of the degeneracies is two. Thus, in our two-parameter system, 2-fold degeneracies occur at isolated points. Therefore, a small readjustment of the two parameters can always compensate a sufficiently small perturbation to the Hamiltonian: the conical intersections will only displace slightly upon inclusion of third- and higher-order corrections of expansion (2). The fact that the expansion parameter $`r`$ at the degeneracies is so small indicates that, most likely, higher-order corrections are indeed smaller that the $`r^2`$ term at the degenerate points and even beyond them. In this case, higher-order terms do not change the qualitative configuration of degenerate points and do not affect therefore the geometrical phases calculated at second order. However, without an upper bound on the higher-order terms, nothing guarantees that these terms do not displace significantly or even remove the degenerate points found at second order.

For this reason, we need to verify with a nonperturbative method the actual presence and position of the degeneracies. We thus resort to the numerical technique introduced by Berry and Wilkinson (BW) to solve the Laplace problem in two-dimensional domains. The advantage of this method is to reduce the two-dimensional differential equation to a one-dimensional integral equation for the normal derivative of the wavefunction at the boundary, which can be solved on a discrete mesh. We assess the presence of each degeneracy by encircling it with a loop in the parameters space and calculating the cyclic Berry phases of the involved states with the BW numerical method. Specifically, for the problem at hand we have calculated with the BW technique the Berry phases for the (discretized) loop $`\mathrm{\Gamma }_1\mathrm{\Gamma }_3`$ in Fig. 4, obtaining $`\gamma _1=1`$, $`\gamma _2=1`$, and $`\gamma _3=+1`$. Similar procedure also confirms the presence of the degenerate point B.

At this stage, one may wonder if satellite degeneracies (i.e. degeneracies within the range of validity of perturbation theory, involving irrelevant components on states outside the multiplet) are present only in this specific distorted “quantum billiard”, or if they are a widespread feature. Consider, for example, the 3-fold degenerate multiplet $`\{n_x,n_y\}=\{1,3\}`$, $`\{4,2\}`$, and $`\{5,1\}`$ in the same geometry $`a/b=\sqrt{3}`$: in this case we find no satellite degeneracies in the neighborhood of the 3-fold one, where second-order perturbation theory holds. This is not surprising, since over the whole range $`0<\theta <2\pi `$, the first-order splittings vary much more slowly than those of Fig. 3. Thus, when the quadratic term becomes of the same order of the linear term, also higher-order terms are of comparable size.

If we could establish a general criterion to predict the occurrence of satellite degeneracies for the general case of a $`n`$-fold degenerate multiplet in a billiard of arbitrary side ratio $`a/b`$, then we could easily realize practical classical/quantum systems with clustered conical degeneracies, and thus intricate patterns of parallel-transported wavefunctions. Satellite degeneracies are to be expected whenever (i) the smallest first-order gap $`\delta _{\pi /2}^{(1)}`$ near $`\theta =\frac{\pi }{2}`$ is much smaller than those $`\delta _0^{(1)}`$ at $`\theta =0`$, and (ii) the second-order matrix elements of $`|^{(2)}(\pi /2)|\delta _{\pi /2}^{(1)}`$, so that for $`r\delta _{\pi /2}^{(1)}/|_{ij}^{(2)}(\pi /2)|1`$ first and second order terms are of the same magnitude and at the same time higher-orders are negligible

We first look for the realization of condition (i). $`^{(1)}(\theta )`$ \[Eq. (3)\] determines the first-order splittings $`\delta _\theta ^{(1)}`$. The diagonal matrix elements of $`F`$ determine the linear splittings $`\delta _0^{(1)}`$ since the off-diagonal elements either vanish or are much smaller. The diagonal matrix elements of $`F^{}`$ are instead all equal. Consequently, at $`\theta =\pi /2`$, the mixing and the splittings $`\delta _{\pi /2}^{(1)}`$ of the three states are governed by the off-diagonal elements of $`F^{}`$. A parity selection rule applies to these matrix elements: $`F_{\{n_x,n_y\}\{n_x^{},n_y^{}\}}^{}`$ is zero when

$$(1)^{n_x+n_x^{}}=(1)^{n_y+n_y^{}}=1.$$
(6)
In addition, the the nonzero elements (setting $`a=b`$ for the sake of an estimate) satisfy

$$|F_{\{n_x,n_y\}\{n_x^{},n_y^{}\}}^{}|\frac{2^{9/2}\frac{n_x}{n_x^{}}\frac{n_y^{}}{n_y}}{\left[1\left(\frac{n_x}{n_x^{}}\right)^2\right]\left[1\left(\frac{n_y^{}}{n_y}\right)^2\right]}.$$
(7)
For states widely separated in $`n_xn_y`$ plane ($`n_x/n_x^{}1`$ and $`n_y^{}/n_y1`$) this coupling becomes very small. Small linear gaps occur at $`\pi /2`$ whenever a state $`\{n_x,n_y\}`$ in the degenerate multiplet is weakly coupled to the others, i.e. when all the off-diagonal matrix elements $`F_{\{n_x,n_y\}\{n_x^{}n_y^{}\}}^{}`$ are much smaller than $`\delta _0^{(1)}`$. Indeed in Lauber’s multiplet the linear coupling between states $`\{5,3\}`$ and $`\{7,1\}`$ vanishes because of parity (6), and the coupling determining $`\delta _{\pi /2}^{(1)}`$ is $`F_{\{2,4\}\{7,1\}}^{}0.9`$, compared to $`\delta _0^{(1)}50`$. As we saw, for the multiplet $`\{1,3\}`$, $`\{4,2\}`$, and $`\{5,1\}`$ the ratio $`\delta _{\pi /2}^{(1)}/\delta _0^{(1)}`$ is much closer to one, and indeed, no satellite degeneracies appear there.

We come now to condition (ii). The matrix elements of $`^{(2)}(\pi /2)`$ follow no selection rule. We have verified that these matrix elements are large: for instance $`|_{\{n,n+1\}\{n+1,n\}}^{(2)}(\pi /2)|n^3`$. These matrix elements reduce rapidly as a function of the distance $`(n_xn_x^{})^2+(n_yn_y^{})^2`$ between the states.

In conclusion, whenever in a degenerate multiplet one state is near some states (so that second-order coupling is large) for which selection rule (6) makes first order coupling vanish, and at the same time it is far from all remaining states (so that $`\delta ^{(1)}`$ is small), one expects satellite degeneracies. This rule permits, for instance, to decide that the 3-fold multiplet $`\{1,5\}`$, $`\{4,4\}`$, $`\{6,2\}`$ for $`a/b=\sqrt{5/3}`$ has a structure of satellite degeneracies very similar to that of Lauber’s experiment, as we then verified by BW calculation.

An interesting case emerges from this analysis: when all off-diagonal matrix elements of $`F^{}`$ vanish because of selection rule (6), the linear perturbation does not remove the degeneracy at $`\theta =\pi /2`$. In this direction, the second-order term $`^{(2)}`$ is therefore the leading responsible for the splitting of the degenerate levels, which therefore “kiss” gently, instead of intersecting conically as usual. Many two-fold degenerate multiplets satisfy this condition (for example $`\{2,4\}`$ and $`\{4,2\}`$ for $`a/b=1`$, or $`\{2,3\}`$ and $`\{6,1\}`$ for $`a/b=2`$).

In this work, we have accounted for all experimental findings regarding the eigenstates of Lauber et al.’s deformed resonator. We have found that in deformed rectangular resonators additional degeneracies (with lower multiplicity) may appear very close to multiple degeneracies. These extra degeneracies are responsible for the observed pattern of geometric phases. These results rely only on the properties of the Laplace operator: for this reason they have a wide scope of applicability and are not restricted to simple microwave resonators.

We thank Dirk Dubbers for suggesting us to investigate this system and for useful discussions.