text:
# Dynamical parameter estimation using realistic photodetection

## I Introduction

Quantum parameter estimation originated with the work by Helevo Hol82  and Helstrom Hel76  and is usually formulated as follows. A known quantum state undergoes some evolution. This evolution is usually unitary but need not be ChuNie97 , and is parameterized by one or more unknown parameters. The goal is to estimate these parameters by making a measurement on the final state. In a recent paper, two of us Jay , following Mabuchi MabuchiEst , considered the situation in which an open quantum system parameterized by an unknown dynamical parameter is estimated through continuous-in-time measurements of the environment. Ref. Jay  also focussed on the conditioned state and how the lack of information about the parameter affects its purity. These two papers are part of a small but growing body of work including DohertyStateDeterm  and DohertyOptQPE  (and references therein) that address the problem of reducing classical uncertainties in the evolution of quantum systems by performing continuous-in-time measurements.

In Ref. Jay , the system is taken to be a damped, classically driven two-level atom (TLA) with a driving Hamiltonian ($`\widehat{\sigma }_x`$ in a rotating frame) with an uncertain driving strength $`\mathrm{\Omega }`$. The physical scenario the authors give for this is as follows. An atom is placed in a classical standing wave, thus experiencing a Rabi frequency of

$$\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{max}}\mathrm{sin}(kx),$$
(1)
where $`k`$ is the wavevector for the classical field and $`x`$ is the position of centre of mass of the atom. Gambetta and Wiseman then assume that the TLA is equally likely to be placed at any position in the standing wave. Hence, a range of $`\mathrm{\Omega }_{\mathrm{max}}\mathrm{\Omega }\mathrm{\Omega }_{\mathrm{max}}`$ exists for the possible driving strength. It is also assumed that the TLA has zero translational motion. This can be justified by assuming that the TLA is fixed by a radio frequency Paul trap or similar confining mechanism MabuchiEst . Although it is shown that continual efficient measurement upon the output electromagnetic field of a TLA leads to information gain about the driving strength, the amount of information gained depends strongly on the measurement scheme used. The purity of the conditional state also depends strongly on the measurement scheme, but the performance of schemes by these different criteria are not well correlated at all Jay .

In this paper, we consider only one sort of detection: direct detection, but we consider imperfect detection. This allows us to address the following two questions. First, what (qualitatively) are the aspects of the ideal measurement record that contain the information about the system state and the unknown parameter? Second, in practice, how badly would the amount of information gained be affected by detectors with realistic imperfections. To treat realistic detection we use the theory developed in Refs. WarWisMab02 ; WarWis03a ; WarWis03b . This generalized the theory of quantum trajectories Car93b ; BelSta92 ; Bar93 ; DalCasMol92 ; GarParZol92 ; Wis96a  by determining the state of the system conditioned upon the output of a detector that has dead-time ($`\tau _{\mathrm{dd}}`$), dark counts at rate $`\gamma _{\mathrm{dk}}`$, inefficiency $`\eta `$, and finite bandwidth $`\gamma _\mathrm{r}`$. A diagram of this detector is displayed in Fig. 1.

The reciprocal of the bandwidth $`\gamma _\mathrm{r}^1`$ sets the scale for the uncertainty in the delay between the absorption of photon (which is synonymous with its detection in an ideal detector) and the avalanche of charge which registers this at a macroscopic level in a realistic detector. We find that a finite bandwidth has the greatest effect on our parameter estimation problem. This tells us that it is the timing between detections, not the steady state detection rate, that contains the majority of the information gained. Unless the bandwidth of the detector is substantially greater than $`\mathrm{\Omega }_{\mathrm{max}}`$, the rate of gain of information about $`\mathrm{\Omega }`$ is much reduced, and consequently the purity of the system is also. Like the work of Ref. Jay  this has important implications for the robustness of quantum feedback control (see for example Refs. Dohe2000 ; WisManWan02 ).

This paper is structured as follows. In Sec. II we begin with some numerical results to motivate the theoretical developments in the remainder of the paper. In Sec. III we briefly review the theory of parameter estimation based on linear quantum trajectories in Ref. Jay , and show how it can be applied to realistic detection. In Sec. IV we use this theory to investigate the evolution of the conditional state for an initially unknown parameter, and in Sec. V we investigate the rate of information gain about the parameter. We conclude with a discussion in Sec. VI.

## II Conditional Probability Distribution for Driving Strength

Before providing a theoretical basis, some numerical results are provided to give the reader an insight into the effects of realistic detection on parameter estimation. With an example in mind, the theory contained in the following section should also prove to be more transparent.

We begin by looking at how an observer’s probability distribution for the driving strength evolves given a single measurement record (trajectory). Details of the numerical techniques employed to do this will be given at the end of the next section. The distribution is denoted by $`P(\mathrm{\Omega }|𝐈_{[0,t)})`$, where $`𝐈_{[0,t)}`$ represents the measurement record. Obviously, $`𝐈_{[0,t)}`$ refers to different events for ideal and realistic detection: it is a record of the times of photodetections and of avalanches, respectively. The initial distribution can be found from Eq. (1) with $`x`$ taken to have a flat distribution. The result is

$$P_0(\mathrm{\Omega })=\frac{1}{\pi \sqrt{\mathrm{\Omega }_{\mathrm{max}}^2\mathrm{\Omega }^2}}.$$
(2)
The unconditioned evolution of the atom (which is independent of the scheme, or quality, of detection) is given by the master equation (ME)

$$\dot{\rho }(t)=\frac{i\mathrm{\Omega }}{2}[\widehat{\sigma }_x,\rho (t)]+\mathrm{\Gamma }𝒟[\widehat{\sigma }]\rho (t)=_\mathrm{\Omega }\rho (t).$$
(3)
Here $`\mathrm{\Omega }`$ is assumed known, $`\mathrm{\Gamma }`$ is the spontaneous emission rate (which we will set equal to unity), $`\widehat{\sigma }`$ is the atomic lowering operator, and $`\widehat{\sigma }_x=\widehat{\sigma }+\widehat{\sigma }^{}`$. The superoperator $`𝒟`$ gives the damping of the system into the environment and is defined as WisMil93c

$$𝒟[\widehat{a}]\rho =\widehat{a}\rho \widehat{a}^{}\frac{1}{2}\left\{\widehat{a}^{}\widehat{a}\rho +\rho \widehat{a}^{}\widehat{a}\right\}.$$
(4)
A comparison of ideal and realistic detection is given in Fig. 2. In this figure $`\mathrm{\Omega }_{\mathrm{true}}=4\mathrm{\Gamma }`$, but the observer does not know this, but rather has knowledge described by the prior distribution (2) with $`\mathrm{\Omega }_{\mathrm{max}}=10\mathrm{\Gamma }`$. As time goes on, the detection record reveals more about $`\mathrm{\Omega }_{\mathrm{true}}`$, and this probability is updated in a Bayesian manner. Note that for both ideal and realistic detectors, the distribution always remains symmetric. This is because the sign of $`\mathrm{\Omega }`$ cannot be determined from direct detection. Photoemissions depend only upon the mean of $`\widehat{\sigma }_z`$, and that is independent of whether the TLA state rotates clockwise or counterclockwise around the $`\widehat{\sigma }_x`$ axis.

In Fig. 2, we have chosen a realistic detector whose only imperfection is a finite bandwidth equal to 7$`\mathrm{\Gamma }`$ (see figure caption). It can be seen that the degradation of the quality of measurement is quite large even though our ‘realistic’ detector is well beyond current technology. If we choose the realistic parameters for an APD (avalanche photodioide) used in Refs. WarWisMab02 ; WarWis03b  we obtain results that are so far from those of ideal detection that a comparison lacks interest.

The fact that the chosen detector parameters give results that are still poor compared to ideal measurement is a little surprising. Very few photodetections are missed under realistic detection with a finite $`\gamma _\mathrm{r}`$ being the only imperfection. It illustrates that the rate of information gain is very sensitive to an uncertainty in the state of the TLA after a jump.

The plot for realistic detection, Fig. 2 (B), also shows that it is difficult to distinguish between $`\mathrm{\Omega }_{\mathrm{true}}`$ and an integer multiple of $`\mathrm{\Omega }_{\mathrm{true}}`$, say $`m\mathrm{\Omega }_{\mathrm{true}}`$ (for this figure this occurs only for $`|\mathrm{\Omega }_{\mathrm{true}}|`$ and $`2|\mathrm{\Omega }_{\mathrm{true}}|`$). By $`\mathrm{\Omega }_{\mathrm{true}}`$ we mean the actual driving strength that exists, but which is initially unknown to the observer. This difficulty exists because the evolution for these two proposed values of the driving strength will place the TLA in the excited state at the same time every integer oscillation (at the smaller frequency). Hence, for certain times there will be a strong correlation between the probabilities of decay. Obviously, a photoemission can still occur when a driving strength of $`2|\mathrm{\Omega }_{\mathrm{true}}|`$ would place the TLA very near ground state, thus ruling out that possible value. Whether these decays that are discriminatory in this sense are detected or not depends on the choice (and duration) of the quantum trajectory in question. In Fig. 2 (B) a small but, persistent (at least up to $`t=50`$) probability for $`|\mathrm{\Omega }|=2|\mathrm{\Omega }_{\mathrm{true}}|`$ exists.

It is worth mentioning that the plots in Fig. 2 are based on the ‘same’ photo-absorptions in the ideal and realistic cases. This is possible using the simulation technique of Refs. WarWis03b , in which the realistic quantum trajectory can be assured to be consistent with the ideal quantum trajectory.

## III Parameter and State Estimation Theory

We now briefly present the theory of parameter and state estimation in the context that is relevant to this paper. The specific problem that we wish to address involves calculating the new probability distribution for $`\mathrm{\Omega }`$ and also obtaining the best estimate of the TLA state, based on the measurement record. Bayesian statistics are required, with the new distribution, based on $`𝐈_{[0,t)}`$, given by BayesBook

$$P(\mathrm{\Omega }|𝐈_{[0,t)})=\frac{P(𝐈_{[0,t)}|\mathrm{\Omega })P_0(\mathrm{\Omega })}{P(𝐈_{[0,t)}|\mathrm{\Omega })P_0(\mathrm{\Omega })𝑑\mathrm{\Omega }}.$$
(5)
The best estimate of the TLA state is found from Jay

$$\rho _𝐈(t)=\frac{\stackrel{~}{\rho }_{𝐈,\mathrm{\Omega }}(t)P_0(\mathrm{\Omega })𝑑\mathrm{\Omega }}{P(𝐈_{[0,t)}|\mathrm{\Omega })P_0(\mathrm{\Omega })𝑑\mathrm{\Omega }},$$
(6)
where $`\stackrel{~}{\rho }_{𝐈,\mathrm{\Omega }}(t)`$ is the unnormalised state of the TLA given that there was a measurement record $`𝐈_{[0,t)}`$ and that the driving strength was $`\mathrm{\Omega }`$.

In this paper we leave the state unnormalised even after jumps, so that

$$\rho _{𝐈,\mathrm{\Omega }}(t)=\frac{\stackrel{~}{\rho }_{𝐈,\mathrm{\Omega }}(t)}{P(𝐈_{[0,t)}|\mathrm{\Omega })}$$
(7)
meaning that $`P(𝐈_{[0,t)}|\mathrm{\Omega })`$ is actually equal to the norm of the TLA state, when it is evolved with $`\mathrm{\Omega }`$ under $`𝐈_{[0,t)}`$. Ideally $`𝐈_{[0,t)}`$ would be an experimental record, but for convenience we simulate it using $`\mathrm{\Omega }_{\mathrm{true}}`$ and then ‘forgetting’ it. The integrals in Eqs. (5)–(6) require the simulation of the evolution for the entire range of $`\mathrm{\Omega }`$. This means that the process is very time consuming as many trajectories (for the different $`\mathrm{\Omega }`$’s) have to be run in order to obtain a single trajectory for $`P(\mathrm{\Omega }|𝐈_{[0,t)})`$ or $`\rho _𝐈(t)`$.

In order to treat realistic detection, it is necessary to describe not just the state of the TLA, but also the state of a detector (treated classically) which is correlated with the atomic state WarWisMab02 ; WarWis03a . We do this by introducing a set of three operators $`\rho _{𝐈,s}(t)`$ for this supersystem, where here $`s=0,1,2`$ represents the three classical detector states. These correspond to the ready state ($`s=0`$), the avalanching state ($`s=1`$) and the resetting state ($`s=2`$). Refer to Fig. 1 for details. The state of the TLA is defined as $`\rho _𝐈(t)=\rho _{𝐈,0}(t)+\rho _{𝐈,1}(t)+\rho _{𝐈,2}(t)`$. Equation (6) is trivially generalized to give the best estimate of the supersystem state $`\rho _{𝐈,s}(t)`$:

$$\rho _{𝐈,s}(t)=\frac{\stackrel{~}{\rho }_{𝐈,\mathrm{\Omega },s}(t)P_0(\mathrm{\Omega })𝑑\mathrm{\Omega }}{P(𝐈_{[0,t)}|\mathrm{\Omega })P_0(\mathrm{\Omega })𝑑\mathrm{\Omega }}.$$
(8)
The three equations, Eqs. (5)–(6) and Eq. (8), are unusable in their present form as $`P(𝐈_{[0,t)}|\mathrm{\Omega })`$ and the norm of $`\stackrel{~}{\rho }_{𝐈,\mathrm{\Omega }}(t)`$ (which is equal to $`P(𝐈_{[0,t)}|\mathrm{\Omega })`$) will become extremely small as $`t`$ increases. This will be true even for $`\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{true}}`$, as can be seen by considering the probability of a string of measurement results occurring. The probability of a jump in a given infinitesimal interval is never greater than $`\mathrm{\Gamma }dt`$, so that the norm of the state after $`m`$ jumps would be less than $`(\mathrm{\Gamma }dt)^m`$ (since the no-jump evolution also decreases the norm of the state). In simulations, a time step of $`dt=10^4\mathrm{\Gamma }^1`$ is typical. Thus, numerical error introduced by computer round-off necessitates a different approach which we now explain.

### III.1 Linear Quantum Trajectories

In linear quantum trajectories Wis96a , the state is multiplied by a predetermined factor after every measurement event. The normalisation constants are chosen so that the norm of the states corresponding to the most likely values of $`\mathrm{\Omega }`$ stay relatively close to unity. Each possible measurement result, $`r`$, has an ‘ostensible’ Wis96a  probability, $`\mathrm{\Lambda }(r)`$ associated with it, that we use as the normalising factor. Therefore, we have

$$\overline{\rho }_{𝐈,\mathrm{\Omega },s}(t)=\frac{\stackrel{~}{\rho }_{𝐈,\mathrm{\Omega },s}(t)}{\mathrm{\Lambda }(𝐈_{[0,t)})},$$
(9)
where $`\overline{\rho }_{𝐈,\mathrm{\Omega },s}`$ is the supersystem state when linear quantum trajectories are used and $`\mathrm{\Lambda }(𝐈_{[0,t)})`$ is the ostensible probability for getting $`𝐈_{[0,t)}`$,

$$\mathrm{\Lambda }(𝐈_{[0,t)})=\mathrm{\Lambda }(r_k)\mathrm{\Lambda }(r_{k1})\mathrm{}\mathrm{\Lambda }(r_1).$$
(10)
Here, we have assumed that the time interval $`[0,t)`$ has been divided into a very large number, $`k=t/dt`$, of discrete measurement times. The actual probability of getting $`𝐈_{[0,t)}`$, given $`\mathrm{\Omega }`$ is

$$P(𝐈_{[0,t)}|\mathrm{\Omega })=\mathrm{\Lambda }(𝐈_{[0,t)})\underset{s}{}\mathrm{Tr}[\overline{\rho }_{𝐈,\mathrm{\Omega },s}(t)],$$
(11)
where the trace is over the TLA and the summation is over the detector states.

This allows Eq. (5) and Eq. (8) to be written as Jay

$$P(\mathrm{\Omega }|𝐈_{[0,t)})=\frac{_s\mathrm{Tr}[\overline{\rho }_{𝐈,\mathrm{\Omega },s}(t)]P_0(\mathrm{\Omega })}{_s\mathrm{Tr}[\overline{\rho }_{𝐈,\mathrm{\Omega },s}(t)]P_0(\mathrm{\Omega })d\mathrm{\Omega }}$$
(12)
and

$$\rho _{𝐈,s}(t)=\frac{\overline{\rho }_{𝐈,\mathrm{\Omega },s}(t)P_0(\mathrm{\Omega })𝑑\mathrm{\Omega }}{_s\mathrm{Tr}[\overline{\rho }_{𝐈,\mathrm{\Omega },s}(t)]P_0(\mathrm{\Omega })d\mathrm{\Omega }}.$$
(13)
Since $`\mathrm{Tr}[\overline{\rho }_{𝐈,\mathrm{\Omega },s}(t)]`$ is of order unity for the most likely $`\mathrm{\Omega }`$’s the problems associated with using $`\mathrm{Tr}[\stackrel{~}{\rho }_{𝐈,\mathrm{\Omega },s}(t)]`$ have been avoided. The best estimate of the system state allows the conditional evolution of a single stochastic record to be investigated, while $`P(\mathrm{\Omega }|𝐈_{[0,t)})`$ will allow calculation of the information gain (see following).

In this paper, we choose the ostensible probabilities to be

$$\mathrm{\Lambda }(1)=ϵdt,$$
(14)
where $`r=1`$ indicates the detection of a photon (ideal detector) or an avalanche (realistic detector). No detection ($`r=0`$) has an ostensible probability $`\mathrm{\Lambda }(0)=1ϵdt`$. With these ostensible probabilities, the linear quantum trajectory equations for the supersystem describing realistic detection WarWisMab02 ; WarWis03a  become

$`d\overline{\rho }_0(t)`$ $`=`$ $`dt\left(_\mathrm{\Omega }\gamma _{\mathrm{dk}}\eta \mathrm{\Gamma }𝒥[\widehat{\sigma }]+ϵ\right)\overline{\rho }_0(t)`$ (15)
$`d𝒩(t)\overline{\rho }_0(t)+d𝒩(t)(t\tau _{\mathrm{dd}})\overline{\rho }_2(t)`$
$`d\overline{\rho }_1(t)`$ $`=`$ $`dt\left(_\mathrm{\Omega }\gamma _\mathrm{r}+ϵ\right)\overline{\rho }_1(t)+dt\left(\eta \mathrm{\Gamma }𝒥[\widehat{\sigma }]+\gamma _{\mathrm{dk}}\right)`$ (16)
$`\times \overline{\rho }_0(t)d𝒩(t)\overline{\rho }_1`$
$`d\overline{\rho }_2(t)`$ $`=`$ $`dt\left(_\mathrm{\Omega }+ϵ\right)\overline{\rho }_2(t)+d𝒩(t)\overline{\rho }_1(t)/ϵ`$ (17)
$`d𝒩(t)(t\tau _{\mathrm{dd}})\overline{\rho }_2(t).`$
Here $`d𝒩(t)=0,1`$, the increment in the number of avalanches, is equivalent to the measurement result $`r`$ in the time interval $`[t,t+dt)`$. We maintain the symbol $`𝐈_{[0,t)}`$ for the record of avalanches in the interval $`[0,t)`$, but for clarity we have omitted the subscripts $`𝐈`$ and $`\mathrm{\Omega }`$ for the $`\rho (t)`$’s. The superoperator $`𝒥`$ is defined by WisMil93c

$$𝒥[\widehat{a}]\rho \widehat{a}\rho \widehat{a}^{}.$$
(18)
For simulation purposes we take

$$ϵ=\eta \frac{\mathrm{\Gamma }\mathrm{\Omega }_{\mathrm{true}}^2}{2\mathrm{\Omega }_{\mathrm{true}}^2+\mathrm{\Gamma }^2},$$
(19)
the expected detection rate in steady state (ignoring any dead time) when $`\mathrm{\Omega }`$ is known. This is allowed as the results of the simulation are independent of $`ϵ`$, thus we are free to choose any value. Note that renormalisation of $`\rho (t)`$ is occurring even when the detector is resetting and the probability of an avalanche is strictly zero. This illustrates the flexibility of the ostensible probabilities, where the only concern is maintaining a norm of order unity.

## IV Conditional State Evolution

Now that we have explained how the numerical results of Sec. II were obtained, we return to these simulations to investigate the conditional state evolution in a typical realistic quantum trajectory. In Fig. 3 we see how $`z(t)`$ (the mean of $`\widehat{\sigma }_z`$) behaves under a variety of circumstances. Specifically, we consider known $`|\mathrm{\Omega }_{\mathrm{true}}|`$ and ideal detection (thin dashed line), unknown $`|\mathrm{\Omega }_{\mathrm{true}}|`$ (that is, an initial distribution given by Eq. (2)) and ideal detection (thin solid line), known $`|\mathrm{\Omega }_{\mathrm{true}}|`$ and realistic detection (thick dashed line) and, finally, unknown $`|\mathrm{\Omega }_{\mathrm{true}}|`$ and realistic detection (thick solid line). The reason why a comparison is made between the trajectories for a TLA randomly placed in the standing wave of the field and trajectories of known $`|\mathrm{\Omega }_{\mathrm{true}}|`$, but unknown sign, is that at $`t=\mathrm{}`$ the sign of $`\mathrm{\Omega }_{\mathrm{true}}`$ will still not be known, as discussed in Sec. II. This allows us to reveal how far the initially unknown $`|\mathrm{\Omega }_{\mathrm{true}}|`$ trajectories are from their optimal form, given the constraints of direct detection.

The conditioned state corresponding to a known $`|\mathrm{\Omega }_{\mathrm{true}}|`$ is formed from

$$\rho _\mathrm{I}(t)=(\rho _{\mathrm{I},\mathrm{\Omega }_{\mathrm{true}}}(t)+\rho _{\mathrm{I},\mathrm{\Omega }_{\mathrm{true}}}(t))/2.$$
(20)
That is, it is an equal mixture of the two conditioned states obtained with a driving strength of $`+\mathrm{\Omega }_{\mathrm{true}}`$ and $`\mathrm{\Omega }_{\mathrm{true}}`$. Given that we take the TLA as starting in the ground state, these two conditioned states will have the same $`z(t)`$, but opposite sign for $`y(t)`$ resulting in $`y(t)=0`$ for the mixture as initially $`P_0(\mathrm{\Omega })`$ is an even mix of positive and negative $`\mathrm{\Omega }`$. If the initial state of the TLA (which is presumed to be known) gives a finite value for $`y(t=0)`$ then the time of the first jump will give a very small amount of information about $`\mathrm{\Omega }_{\mathrm{true}}`$. This is because an early jump would tend to support the conclusion that the TLA had been rotated on the Bloch sphere in the direction that takes it towards the excited state first, while a later jump suggests that the TLA was rotated in the opposite direction. However, the initial state is taken to be the ground state in our simulations. Since for direction detection with no detuning $`x(t)=0`$ Car93b ; WisMil93c , we have

$$\rho _\mathrm{I}(t)=(I+z(t)\widehat{\sigma }_z)/2,$$
(21)
where $`z(t)`$ could be obtained from either the $`+\mathrm{\Omega }_{\mathrm{true}}`$ or $`\mathrm{\Omega }_{\mathrm{true}}`$ trajectory. It should now be clear why only $`z(t)`$ is plotted. In fact, the purity can also be inferred from the $`z(t)`$ plot as $`p_\mathrm{I}(t)=(1+z^2(t))/2`$.

As for Fig. 2, the realistic trajectory is forced to be consistent with the ideal trajectory. This is evidenced by the jump times being correlated. Because of the detail of these plots we show two portions of the evolution only, illustrating how improved knowledge over time of the dynamical parameter, $`\mathrm{\Omega }`$, allows closer adhesion to the maximum knowledge trajectory (thin dashed line for ideal detection and thick dashed line for realistic detection).

In the first $`10\mathrm{\Gamma }^1`$ units of time the frequency of oscillation of $`z(t)`$ for the realistic detection state (solid black line) is poorly defined. One can also see that the frequency is faster than the true frequency, due to $`P_0(\mathrm{\Omega })`$ being peaked at $`\pm \mathrm{\Omega }_{\mathrm{max}}`$ Jay . The small amplitude of oscillation is also due to the lack of knowledge about $`\mathrm{\Omega }`$. The range of possible frequencies ensures that when one value of $`\mathrm{\Omega }`$ would place the TLA in the excited state, another would place it near the ground state, thus averaging the amplitude of oscillation of $`z(t)`$ to below unity.

For ideal detection (thin solid line), the frequency of oscillation is also shifting about, but the amplitude of the oscillations have already noticeably increased in size by $`t=10\mathrm{\Gamma }^1`$. Another feature is that jumps take $`z(t)`$ to $`1`$ (a pure state) for ideal detection even if $`|\mathrm{\Omega }_{\mathrm{true}}|`$ is not known. This is because there is no time delay in conveying the decay of the TLA to the observer. This is not true for realistic detection where the jump is slightly delayed and the TLA will have rotated out of the ground state by some unknown amount.

After $`40\mathrm{\Gamma }^1`$ time units, the ideal detection observer has gained enough information about $`|\mathrm{\Omega }_{\mathrm{true}}|`$ that $`z(t)`$ matches well with the ‘ideal’ $`z(t)`$ (Fig. 3 (B)). The realistic observer’s $`z(t)`$ is much closer to the ideal than at the start of the trajectory, but some uncertainty still obviously exists.

## V Information Gain

We now move away from conditional dynamics and look at how the ensemble averaged gain of information concerning $`\mathrm{\Omega }`$ is affected by realistic detection. As in Jay , the measure used to quantify the quality of information gained about $`\mathrm{\Omega }`$ is $`\mathrm{\Delta }I_𝐈`$. It is given, in bits, by

$`\mathrm{\Delta }I_𝐈`$ $`=`$ $`{\displaystyle 𝑑\mathrm{\Omega }P(\mathrm{\Omega }|𝐈_{[0,t)})\mathrm{log}_2P(\mathrm{\Omega }|𝐈_{[0,t)})}`$ (22)
$`{\displaystyle 𝑑\mathrm{\Omega }P_0(\mathrm{\Omega })\mathrm{log}_2P_0(\mathrm{\Omega })}.`$
As the distribution for $`\mathrm{\Omega }`$ becomes more sharply peaked around $`|\mathrm{\Omega }_{\mathrm{true}}|`$, $`\mathrm{\Delta }I_𝐈`$ will increase. Due to the stochastic nature of $`𝐈_{[0,t)}`$, an ensemble of trajectories are run in order to obtain the average, $`\mathrm{\Delta }I`$, of $`\mathrm{\Delta }I_𝐈`$ that is characteristic of the measurement scheme. Each member of the ensemble is formed by picking an $`\mathrm{\Omega }_{\mathrm{true}}`$ randomly (but according to the distribution Eq. (2)) and simulating the evolution based on a stochastic trajectory corresponding to that $`\mathrm{\Omega }_{\mathrm{true}}`$. A single trajectory for each member of an ensemble of $`\mathrm{\Omega }_{\mathrm{true}}`$’s is sufficient as the averaging of the stochasticity is done over the multitude of $`\mathrm{\Omega }_{\mathrm{true}}`$’s.

A major difference between parameter estimation for ideal direct detection and realistic direct detection is found when $`\mathrm{\Omega }_{\mathrm{true}}\mathrm{\Gamma }`$. We investigate this regime by varying the value of $`\mathrm{\Omega }_{\mathrm{max}}`$ used in the simulations. This works as the most likely values for $`\mathrm{\Omega }_{\mathrm{true}}`$ come from close to $`\mathrm{\Omega }_{\mathrm{max}}`$. By increasing $`\mathrm{\Omega }_{\mathrm{max}}`$ we are effectively increasing the values of $`\mathrm{\Omega }_{\mathrm{true}}`$ used. This method corresponds physically to the TLA being placed randomly in standing waves of different amplitude.

In Fig. 4 (A) the information gain, $`\mathrm{\Delta }I`$, averaged over an ensemble size of $`1000`$ (that is, $`1000`$ different $`\mathrm{\Omega }_{\mathrm{true}}`$’s), is shown for ideal detection. The three plots are for $`\mathrm{\Omega }_{\mathrm{max}}=5\mathrm{\Gamma },10\mathrm{\Gamma },20\mathrm{\Gamma }`$. In Fig. 4 (B) the results are shown for the same values of $`\mathrm{\Omega }_{\mathrm{max}}`$ when the monitoring is done with realistic detection. Error bars are included as the ensemble size in this case is only $`100`$, due to the greater computational intensity. As $`\mathrm{\Omega }_{\mathrm{max}}`$ doubles, the information gain curves for ideal detection are separated by approximately $`1`$ bit at large $`t`$. This does not indicate that the detection scheme is more effective at large $`\mathrm{\Omega }_{\mathrm{max}}`$, rather, it is an artifact of the initial information content of the distribution $`P_0(\mathrm{\Omega })`$. This can easily be shown analytically. Consider the integral

$$_{\mathrm{\Omega }_{\mathrm{max}}}^{\mathrm{\Omega }_{\mathrm{max}}}𝑑\mathrm{\Omega }\frac{1}{\pi \sqrt{\mathrm{\Omega }_{\mathrm{max}}^2\mathrm{\Omega }^2}}\mathrm{log}_2\frac{1}{\pi \sqrt{\mathrm{\Omega }_{\mathrm{max}}^2\mathrm{\Omega }^2}}.$$
(23)
If the change $`\mathrm{\Omega }_{\mathrm{max}}^{}m\mathrm{\Omega }_{\mathrm{max}}`$ is made, a change of integration variable to $`\mathrm{\Omega }^{}=m\mathrm{\Omega }`$ allows the information of the new distribution to be expressed in terms of the old one. The result is

$$\mathrm{\Delta }I_{𝐈,(m\mathrm{\Omega }_{\mathrm{max}})}=\mathrm{\Delta }I_{𝐈,(\mathrm{\Omega }_{\mathrm{max}})}\mathrm{log}_2m.$$
(24)
Thus, every time $`\mathrm{\Omega }_{\mathrm{max}}`$ is doubled, the initial information decreases by $`1`$ bit. Since the initial information is being subtracted from the information at time $`t`$, we can conclude that the width of the $`P(\mathrm{\Omega })`$ distributions, for different $`\mathrm{\Omega }_{\mathrm{max}}`$, in the long time limit are roughly equal for ideal detection.

This is qualitatively different from the case of realistic detection where the information gain for $`\mathrm{\Omega }_{\mathrm{max}}=5\mathrm{\Gamma }`$ and $`10\mathrm{\Gamma }`$ are roughly the same and $`\mathrm{\Omega }_{\mathrm{max}}=20\mathrm{\Gamma }`$ sees a significant decrease. Taking into account the fact that as $`\mathrm{\Omega }_{\mathrm{max}}`$ increases less information is being subtracted away from the information at time $`t`$, it is clear that realistic detection is becoming considerably worse at distinguishing between candidates for $`\mathrm{\Omega }_{\mathrm{true}}`$. This is more clearly illustrated in Fig. 4 (C) which shows the difference between ideal and realistic detection information gain for the three values of $`\mathrm{\Omega }_{\mathrm{max}}`$. We now try to understand these trends by discussing what features of a measurement record allow the information to be accrued.

The most simple (and naive) way to determine the actual driving strength is to analyse the mean steady-state photon counting rate $`𝐈_{\mathrm{ss}}`$. For ideal direct detection this is

$$𝐈_{\mathrm{ss}}=\frac{\mathrm{\Gamma }\mathrm{\Omega }^2}{2\mathrm{\Omega }^2+\mathrm{\Gamma }^2}.$$
(25)
By using Eq. (25) the mean flux can be matched with a particular $`\mathrm{\Omega }`$, given that $`\mathrm{\Gamma }`$ is known. However, the flux asymptotes to $`\mathrm{\Gamma }/2`$ at large $`\mathrm{\Omega }`$, essentially making one value of $`\mathrm{\Omega }`$ indistinguishable from another. Because ideal detection still gives information at large $`\mathrm{\Omega }_{\mathrm{max}}`$ we deduce that the specific times at which detections are occurring also play a role in determining $`\mathrm{\Omega }_{\mathrm{true}}`$.

To consider this effect in more detail we look at the waiting time distribution for photon detections. In perfect direct detection wtime

$$w(\tau )=\frac{\mathrm{\Omega }^2}{\mathrm{\Omega }^2(\mathrm{\Gamma }/2)^2}e^{\mathrm{\Gamma }\tau /2}\mathrm{sin}^2\{\frac{1}{2}[\mathrm{\Omega }^2(\mathrm{\Gamma }/2)^2]^{1/2}\tau \},$$
(26)
where $`w(\tau )`$ is the probability density of a photon being emitted by the TLA at a time $`\tau `$ after the previous emission. The waiting time, $`w(\tau )`$, is zero at $`\tau =0`$ as the TLA in the ground state after a detection. A plot of this function for $`\mathrm{\Omega }=5\mathrm{\Gamma }`$, with $`\mathrm{\Gamma }=1`$ is given in Fig. 5 (A). The point that is meant to be taken from this graph is that there are oscillations that make certain emission times more likely (when the TLA is in the excited state) and other times prohibited (when it is in the ground state), given a certain driving strength. These oscillations, which take $`w(\tau )`$ to zero, are present for all $`\mathrm{\Omega }>\mathrm{\Gamma }/2`$, which includes the regime of interest here. This allows the narrowing of the distribution for $`\mathrm{\Omega }`$, based on a measurement record, despite the fact that very little information can be obtained from the mean flux rate.

The waiting time distribution for realistic detection is much more complex. In this case, $`\tau `$ would be the waiting time between avalanches of the photodiode, rather than emissions of the TLA. One of the difficulties lies in the fact that there is not a definite state of the atom after an avalanche, due to the finite response time. Thus, initial conditions for the TLA must be chosen as being characteristic of the post-avalanche state. We chose them to be the solution to the ME (3) averaged over the time for the avalanche to mature. That is,

$$\rho _{\mathrm{aval}.}=\gamma _\mathrm{r}_0^{\mathrm{}}e^{\gamma _\mathrm{r}t}\rho (t)𝑑t,$$
(27)
where $`\rho (t)`$ is the solution to the ME given that at $`t=0`$ the TLA is in the ground state. Of course, this will only be true if the dark counts are ignored, which we will do. The integral in Eq. (27) can be carried out analytically with $`x=0`$ and the $`y`$ and $`z`$ components of the state being

$`y_{\mathrm{aval}.}`$ $`=`$ $`{\displaystyle \frac{2\mathrm{\Omega }(1+\gamma _\mathrm{r})}{(1+3\gamma _\mathrm{r}+2\gamma _\mathrm{r}^2+2\mathrm{\Omega }^2)}},`$ (28)
$`z_{\mathrm{aval}.}`$ $`=`$ $`{\displaystyle \frac{(2\gamma _\mathrm{r}+1)(1+\gamma _\mathrm{r})}{(1+3\gamma _\mathrm{r}+2\gamma _\mathrm{r}^2+2\mathrm{\Omega }^2)}},`$ (29)
where we have set $`\mathrm{\Gamma }=1`$.

If we take the simplest case where $`\tau _{dd}=0`$ and $`\gamma _{dk}=0`$, then we can use these as the initial conditions for the ready state of the detector. (The case where $`\tau _{dd}0`$ is only a little more complicated, because the solutions in Eqs. (28)–(29) can be evolved forward a time $`\tau _{\mathrm{dd}}`$ with the ME to give the state of the TLA when the detector becomes ready to detect photons once more.) A series of $`6`$ linear coupled differential equations must be solved in order to determine $`\gamma _\mathrm{r}\mathrm{Tr}[\stackrel{~}{\rho }_1(\tau )]=w(\tau )_{\mathrm{aval}.}`$, which is equal to the probability of an avalanche at time $`\tau `$ (the waiting time). Here, $`\stackrel{~}{\rho }_1(\tau )`$ is the solution to the unnormalised evolution equations given in Eqs. (15) and (16) with $`ϵ=0`$ and $`d𝒩0`$ (no avalanches) from $`\tau _{\mathrm{dd}}`$ to $`\tau `$.

The $`6`$ real variables that are necessary can be expressed as a vector $`\stackrel{}{r}=(\stackrel{~}{y}_0,\stackrel{~}{z}_0,\stackrel{~}{P}_0,\stackrel{~}{y}_1,\stackrel{~}{z}_1,and\stackrel{~}{P}_1)^\mathrm{T}`$, where the the subscripts refer to the detector state and the $`P`$’s are the probability of occupation. The equations of motion for this vector is $`d\stackrel{}{r}/dt=𝐀\stackrel{}{r}`$ where $`𝐀`$ is the following matrix

$$\left(\begin{array}{cccccc}\frac{1}{2}& 0& \mathrm{\Omega }& 0& 0& 0\\ 0& \frac{1}{2}\gamma _\mathrm{r}& 0& \mathrm{\Omega }& 0& 0\\ \mathrm{\Omega }& 0& 1+\frac{\eta }{2}& 0& 1+\frac{\eta }{2}& 0\\ 0& \mathrm{\Omega }& \frac{\eta }{2}& 1\gamma _\mathrm{r}& \frac{\eta }{2}& 1\\ 0& 0& \frac{\eta }{2}& 0& \frac{\eta }{2}& 0\\ 0& 0& \frac{\eta }{2}& 0& \frac{\eta }{2}& \gamma _\mathrm{r}\end{array}\right),$$
(30)
whereas before we are assuming $`\gamma _{\mathrm{dk}}=0`$ and $`\mathrm{\Gamma }=1`$. The variable of interest is then $`w(\tau )_{\mathrm{aval}}=\gamma _\mathrm{r}\stackrel{~}{P}_1(\tau )`$

Plots of $`w(\tau )_{\mathrm{aval}.}`$ with $`\mathrm{\Omega }=5\mathrm{\Gamma }`$, $`10\mathrm{\Gamma }`$, and $`20\mathrm{\Gamma }`$ are given in Fig. 5 (B), (C) and (D). Most notable is that although there are oscillations, they do not take $`w_{\mathrm{aval}.}(\tau )`$ to zero. The random response time of the photodetector is washing out the peaks and troughs of the ideal detection waiting time distribution, making it much more difficult to distinguish between $`\mathrm{\Omega }`$’s on the basis of the specific times of the avalanches. This effect becomes more pronounced as $`\mathrm{\Omega }`$ increases. The relevant characteristic time scale is the response time of the detector. For $`\mathrm{\Omega }\gamma _\mathrm{r}`$, the TLA may have rotated out of the ground state significantly by the time the avalanche occurs. If the state of the TLA is smeared out and information contained in the higher order moments of the waiting distribution is lost, then there is no way to determine $`|\mathrm{\Omega }_{\mathrm{true}}|`$ when it is large. This is the reason for the drop-off of the information gain for realistic detection that is seen in Fig. 4.

Although we have done our simulations with a dead time $`\tau _{dd}=0`$, the finite response bandwidth leads to an effective dead time of $`\gamma _r^1`$, which leads to an effective drop in the efficiency of the detector. To rule inefficiency as the cause of the loss of structure in $`w(\tau )`$ as $`\mathrm{\Omega }`$ increases we consider the waiting time distribution for a detector of efficiency $`\eta `$. Using the formula given in Ref. wtime  we have calculated $`w(\tau )`$ with $`\mathrm{\Omega }=5\mathrm{\Gamma }`$ , $`10\mathrm{\Gamma }`$, and $`20\mathrm{\Gamma }`$. Plots are given in Fig. 6 (B), (C) and (D) (with $`\eta =70\%`$) while Fig. 6 (A) is a reproduction of Fig. 5 (A) ($`\mathrm{\Omega }=5\mathrm{\Gamma }`$ and $`\eta =100\%`$). Here it is observed that the smearing out of $`w(\tau )`$ (for detectors with only an inefficiency) is due only to this inefficiency. That is, it is independent of $`\mathrm{\Omega }`$. Thus the smearing out the waiting time distribution as a result of increasing $`\mathrm{\Omega }`$ for a realistic detector is brought about by the finite bandwidth, $`\gamma _r`$.

## VI Discussion

We have investigated the effect of imperfections in realistic detectors upon the problem of quantum state and parameter estimation by continuous monitoring of an open quantum system. Specifically, we have re-examined the system of a two-level atom (TLA) with an unknown driving strength introduced in Ref. Jay . Considering only direct photodetection, we find that the most significant effect comes from a finite bandwidth of the detector. This corresponds to a randomness in the response time of the photodiode (the time from photon absorption to the time when the avalanche reaches a pre-set threshold).

We find that unless the bandwidth is significantly greater than the Rabi frequency, the observer’s ability to obtain information about the unknown Rabi frequency, and about the state of the atom, is severely compromised. This is because the waiting time distribution between avalanches is smeared out by the bandwidth, losing the Rabi oscillations that characterize it for an ideal detector (see Fig. 5). This result shows that the Bayesian update method in Ref. Jay  implicitly made use of the structure of the waiting time distribution, rather than simply the mean detection rate. This implies that even for an ideal detector this parameter estimation problem is extremely nonlinear, and requires the full power of Bayesian probability theory. Our results thus have implications for quantum control Dohe2000  in the presence of unknown parameters even for ideal detectors.

A final comment on the question of the physicality of the scenario that we have constructed needs to be made. In order to measure the emissions of an atom with an efficiency anywhere close to unity it would be necessary to place the atom in a micro cavity, with the cavity mode being heavily damped so that the atom can be regarded as emitting predominantly into the cavity output beam RicCar88 ; TurThoKim95 . However, in this paper we have considered a situation in which the position of the TLA is random, but fixed, across one wavelength of a standing wave field. This has the ramification that the effective decay rate is also position-dependent and hence also uncertain, scaling like the driving strength squared. Thus, in reality, there would be uncertainty in the two dynamical parameters $`\mathrm{\Omega }`$ and $`\mathrm{\Gamma }`$. In principle, with ideal detection, both these parameters could be found (apart from the sign of $`\mathrm{\Omega }`$) given a long enough measurement record, as the waiting time distribution for ideal detection will be different for different positions of the TLA in the standing wave. It would be interesting to see how realistic direct detection affects the rate gain of information about $`\mathrm{\Omega }`$ and $`\mathrm{\Gamma }`$.

Another question which would be interesting to answer is, given that homodyne detection of the $`y`$ quadrature in Jay  was shown to provide the greatest gain in information, does this result also apply when realistic detection affects are taken into account. Lastly, an interesting idea would be to investigate this scenario with a frequency-filtered measurement scheme, such as in Ref. WisToo99 . This would require extending the measuring apparatus to include many cavities (one for each frequency) with each cavity having its own realistic detector. Either of these schemes would require far greater numerical resources than the (not insignificant) resources used here.

###### Acknowledgements.

This work was supported by the Australian Research Council.


format: text

text:
# Internet’s Critical Path Horizon

## 1 Introduction

The efficient performance of any communication network is jeopardized by congestion problems, which often show up in unpredictable ways. This seems the case, for example, of so called Internet storms Huberman . Such problems were early identified in different types of engineered networks. Norbert Wiener for example mentions that “a switching service involving many stages and designed for a certain level of failure shows no obvious signs of failure until the traffic comes up to the edge of the critical point, when it goes completely into pieces, and we have a catastrophic traffic jam” Wiener . These observations allow to formulate a number of key questions concerning communication nets, such as: How do critical traffic levels are reached? What is the nature of these thresholds? How appropriate routing algorithms modify this behavior? Are there optimal routing strategies?

Modeling Internet dynamics has been an active area over the last decade. The approaches include detailed simulations Ogielski , simple statistical models and verbal models WillingerPNAS . In this context, in SoleValverde2001  we investigated the existence of a jamming phase transition between the free phase and the congested phase in a model of network traffic over regular meshes. This phase transition depends on traffic density. It was conjectured that large-scale network traffic self-organizes near the critical point of the transition, which is linked to high network efficiency and unpredictability ValverdeSole2002 . At the critical regime, the distribution of congestion duration lengths scales as a power-law Fukuda ValverdeSole2002 .

Such a self-organized scenario is reinforced by a recent study suggesting that Internet fluctuations are a consequence of the internal dynamics of the system Fluctuations . Against previous claims, there is mounting empirical support that network traffic heterogeneity is a consequence of collective dynamics and not because of the high variability injected by external sources.

In this paper we will show that Internet is able to route efficiently its inner flow of packets because a special combination of local routing rules and a particular network architecture. As we will see, Internet routing reaches an optimal, low-cost traffic flow as a result of a trade-off between random and deterministic routing schemes.

## 2 Modelling Internet’s traffic

Following previous approaches SoleValverde2001 ValverdeSole2002  let us consider network defined on a graph $`\mathrm{\Omega }`$ with M nodes (figure 1). The number of links connecting a given node with another (irrespective of their characteristics) will be indicated as $`k_i`$. The subset of closest nodes of $`s_i\mathrm{\Omega }`$ is $`C_i`$. Only a fraction of fixed $`\rho <(0,1]`$ nodes is selected as source and destination traffic endpoints (hosts). Host locations are chosen randomly. The other $`(1\rho )M`$ nodes can only store and forward incoming messages (routers).

Both host and routers can route only one packet at a time (the routing policy is described below). Each node $`s_i`$ is provided with a finite queue of packets waiting for communication resources (free links). Queues are not allowed to store more than $`H`$ packets simultaneously. New packets arriving at an already saturated queue will be simply removed from the system. This provides a very simple dissipative mechanism useful for recovering from heavy traffic congestion. If $`n(s_i,t)`$ is the number of packets at $`s_i`$, the total number of packets in the system will be

$$N(t)=\underset{s_i\mathrm{\Omega }}{}n(s_i,t)$$
(1)
which is the key quantity which has to be analysed here.

Any microscopic host behavior compatible with the fluid scenario must ensure a controlled injection of packets in order to not flooding the network (and thus entering the congested state). Note that in a fluid traffic regime it is unlikely that packets will be alive for an arbitrarily long time. A fluidity requirement will necessarily impose a hard constraint on the maximum time spent by a packet travelling across the network. The simplest way of ensuring fluid packet flow is to stop the sources emitting new packets when detecting local congestion, that is, when there is no empty space in the source neighborhood $`C_i`$ to fill with new packets. If no self-regulation is present, it has been shown that there is a sudden and sharp jamming transition from the free to the congested state depending on the mean (system) packet generation rate $`\lambda `$. At the critical point between the two phases the system displays optimum (global) performance. If self-regulation is allowed, it can be shown ValverdeSole2002  that simple traffic source rules are able to self-organize the traffic network around a definite mean rate $`\lambda `$, which scales as the quotient between the average system load $`<N>`$ and the averaged packet latency $`T`$:

$$\lambda _C\frac{<N>}{<T>}$$
(2)
where the latency is defined as the time comprised from the creation of a packet until its delivering at the destination host. The mean latency $`<T>`$ is averaged over all successfully released packets. Following ValverdeSole2002  we will indicate by $`\xi `$ the number of local congested neighbors:

$$\xi =\underset{kC(i)}{}\theta \left[n(i,t)\right]$$
(3)
where $`\theta [x]=1`$ for $`x>0`$ and zero otherwise. The packet injection rate for host i is updated as follows:

$$\lambda _i(t+1)=\{\begin{array}{cccccccccccccccccccc}\mathrm{min}\{1,\lambda _i(t)+\mu \}& \xi =0& & & & & & & & & & & & & & & & & & \\ \lambda _i(t)& 0<\xi <k_i& & & & & & & & & & & & & & & & & & \\ 0& \xi =k_i& & & & & & & & & & & & & & & & & & \end{array}$$
(4)
where $`\mu `$ is the so-called driving parameter. The second rule allows a particular host to stabilize around a given rate. Traffic rate increases conservatively and drops down to zero when all neighboring nodes are congested. The reader must be aware that the above rules are not intended to be a detailed model of real traffic sources. Traffic sources can not be described with an universal distribution probability. Moreover, it seems that a rich variety of distributions (exponential, bi-modal or log-normal) apply. As we will see later in the paper, it is unlikely that the variability of traffic sources will be the cause of the scaling detected in Internet traffic. Moreover, this model does not introduce any explicit correlations between sources (i.e.: like in active conversation between two hosts) so the dynamics can not be the by-product of pre-defined source correlations. Because the system is poised to criticality (the so-called fluid packet flow regime) the global dynamics emerge from the collective behavior of its components. Ultimate source details are irrelevant in this context. Moreover, although TCP is more realistic than our simplified local rules (also much more complex) our conjecture is that the large-scale dynamics will not be greatly changed by the local behavior of hosts within this fluid flow regime. This idea has been partly answered in a recent work Hohn  in which it has been shown that IP level traffic (packet dynamics) does not depend to any significant extent on the TCP arrival process (host behavior), supporting our view that traffic dynamics is not a consequence of detailed source (host) behavior.

For a homogeneous network with average degree $`<k>`$ (and finite $`<k^2>`$) it can be shown, following ValverdeSole2002  that the time evolution of packet density $`\mathrm{\Gamma }(t)=N(t)/M`$ in the limit of $`H\mathrm{}`$ is defined by the mean field equation:

$$\frac{d\mathrm{\Gamma }}{dt}=\rho \lambda \frac{<k>}{D}\mathrm{\Gamma }\left(1\mathrm{\Gamma }\right)$$
(5)
where $`D`$ indicates the network diameter (also the average transient time, assuming packets jump from node to node with the same average time scale). It is not difficult to show that the equilibrium points of the previous equation are given by $`\mathrm{\Gamma }_\pm =[1\pm (14\rho \lambda D/<k>)^{1/2}]/2`$. For $`\lambda >\lambda _c<k>/4D\rho `$, the fixed points vanish and no finite density exists. In this “congested phase” the density of packets grows without bounds. For $`\lambda <\lambda _c`$, a finite stable density

$$\mathrm{\Gamma }_{}=\frac{1}{2}\left[1\left(1\frac{4\rho \lambda D}{<k>}\right)^{1/2}\right]$$
(6)
is observable (the other fixed point $`\mathrm{\Gamma }_+`$ is unstable). For the particular case $`\rho =1`$ analysed by Fukś and Lawnizak Fuks , we recover their critical point $`\lambda _c=2/L`$ for dynamics taking place on a square lattice, i. e. $`D=L/2`$. If a feedback exists between $`\lambda `$ and $`\gamma `$, some finite equilibrium density $`\mathrm{\Gamma }^{}`$ will be achieved in the previous model and thus no divergence will be allowed to occur. In this case, the mean field model indicates that a scaling relation will be observed, i. e.

$$\lambda \rho ^1$$
(7)
between the (self-organized) packet release rate and the density of hosts (which is here the only relevant external parameter). Such a scaling relation was shown to occur in the previous model.

The previous mean field calculation was done by assuming that a fixed $`\lambda `$ is beign used. The previous rules actually introduce self-regulation of injection rates by traffic. In other words, if traffic level $`N`$ defines an order parameter, it will interact with a control parameter ($`\lambda `$), reducing it when $`N`$ is large and incresing it when low. The feedback between these two key quantities results in a state dominated by fluid, but fluctuating traffic with many characteristics in common with observed Internet’s dynamical patterns.

In order to explictly consider the feedback between order and control parameters, we can consider a new mean field approximation based on the previous local rules. It can be shown, asuming finite $`H`$, that the new set of equations is now:

$$\frac{d\mathrm{\Gamma }}{dt}=\rho \lambda \left(1\frac{\gamma }{H}\right)\frac{<k>}{D}\mathrm{\Gamma }$$
(8)
$$\frac{d\lambda }{dt}=\mu (1\lambda )\frac{\mathrm{\Gamma }}{<k>}$$
(9)
for low density levels (i. e. $`\mathrm{\Gamma }H`$, consistently with a fluid traffic) the single fixed point (obtained from $`d\mathrm{\Gamma }/dt=d\lambda /dt=0`$) is

$$(\mathrm{\Gamma }^{},\lambda ^{})=(\frac{1}{\frac{<k>}{\rho D}+\frac{1}{\mu <k>}},1\frac{\mathrm{\Gamma }^{}}{<k>M})$$
(10)
the Jacobi matrix $`𝐋`$ for the previous set of equations is given by

$$𝐋=\left(\begin{array}{cc}\frac{<k>}{D}& \rho \\ \frac{1}{<k>}& \mu \end{array}\right)$$
(11)
The associated eigenvalues are

$$\mathrm{\Lambda }_\pm =\frac{1}{2}\left[\left(\frac{<k>}{D}+\mu \right)\pm \sqrt{\left(\frac{<k>}{D}+\mu \right)^2\frac{4\rho }{<k>}}\right]$$

both of them real and negative: the point attractor is globally stable. Numerical simulations of the model on a Poissonian graph (where the previous approximation would held) agree with these predicted values. Topological features are included only as averaged quantities, here mean degree $`<k>`$ and diameter $`D`$. However, our interest is to explore the traffic dynamics on a realistic network architecture, in order to provide the closest modeling approach under the previous rules. It has been shown that the Internet does not display the homogeneous architecture assumed by the Poissonian graphs FFF99 . In the next section the network’s topology used in our analysis is presented. Together with an explicit definition of the routing algorithm, they will complete our model’s description.

## 3 Network topology

In previous analyses SoleValverde2001  ValverdeSole2002  the topology of the graph was chosen to be a lattice. Although this might seem a limited topological arrangement, it has been successfully tested on real hardware and the presence of a phase transition fully confirmed Bolding . Moreover, mesh connected topologies may be the most efficient solution at the limit of very large parallel computer architectures Vitanyi Bilardi .

Regular lattices fail to describe Internet topology. Remarkably, Internet displays a scale-free architecture FFF99 . The origins and causes of this topology have been the subject of most discussion and controversy. It has been show that most existing Internet generators fall in a very different region of the phase space where real Internet is located Yook . This phase space is defined by $`(D_f,\sigma ,\alpha )`$, where every parameter defines a major force shaping a different aspect of the large-scale Internet topology. Note that this model does not define all detailed correlations observed in Internet and/or the precise functional form of Internet’s path length and degree distribution (i.e.: the exact exponent of the scale-free distribution). This is a minimal set of universal parameters that any realistic Internet model must satisfy and it is a very good approximation to the large-scale topology.

The spatial distribution of Internet nodes is not random. It has been noticed a strong correlation between fractal distribution of cities and Internet nodes. The measured fractal dimension is $`D_f=1.5`$. When generating the Internet-like network, the position of nodes is obtained by sampling a Rayleigh-Lévy dust of the same fractal dimension. The likelihood of placing a link between two nodes $`s_i`$ and $`s_j`$ depends both on the (euclidean) link length $`w_{ij}`$ and linear preferential attachment:

$$\mathrm{\Pi }(k_j,d_{ij})\frac{k_j^\alpha }{w_{ij}^\sigma }$$

Note that longer links cost more and thus will be selected with less probability. Traditional topology generators based on Waxman model Waxman  wrongly assume exponential decay instead of linear cost decay. On the other side, there is empirical evidence that highly connected nodes will be linked with higher probability. Increasing $`\alpha `$ will favor linking to nodes with higher degree, while a higher $`\sigma `$ will penalize longer links. From real measurements Yook , it has been identified the position of Internet at $`D_f=1.5`$, $`\sigma =1`$ and $`\alpha =1`$. Here, all numerical simulations have been performed using this topological arrangement (figure 1).

## 4 Path horizon and routing tables

An important ingredient when modeling network dynamics is to describe the paths followed by packets towards their destinations, that is, the routing policy. By properly defining this routing algorithm, we will complete our model’s definition.

Real routing protocols do not drive packets at random. Instead, they try to route packets along the most efficient routes (i.e: minimize distance or latency). At the same time, it is unrealistic to assume that all packets follow optimal paths because the large amount of global information replicated at every single node. Clearly, the paths traced by packets in real networks are properly characterized as a trade-off between random diffusion and optimal routing. This is reflected in the real Internet by its two-level routing organization. Nodes are grouped in so-called Autonomous Systems (AS) and different routing rules apply at each level. Intra-AS routing is based on shortest path routing but inter-AS routing does not cleary follow any minimization criteria.

A simple way to explore the cost/efficiency trade-off is by introducing a parameter defining the visibility scope of the node (depth of routing parameter $`m`$ or node domain diameter), that is, a sphere $`\mathrm{\Gamma }^{(m)}(s_i)`$ of radius $`m`$ centered at every node $`s_i\mathrm{\Omega }`$ (figure 1). We allow every node to known every other node at a distance of $`m`$ hops or less but no more. No information will be stored about nodes outside the node domain. This idea is indicated in figure 2, where a given target node $`s_i`$ is shown at the center of its $`m`$-sphere. When a foreign node $`s_j\mathrm{\Omega }\mathrm{\Gamma }^{(m)}(s_i)`$ send a packet towards $`s_i`$, while moving in the outside of the sphere it performs as a random walk. Once the packet hits the boundary of the sphere, $`\mathrm{\Gamma }^{(m)}(s_i)`$, it is routed along the shortest path NOTE1 .

## 5 Efficiency and network’s exploration

The depth of routing parameter m induces a hierarchy of path subsets over the entire set of available network paths (see figure 3). The random $`m=0`$ and deterministic $`m=M`$ routing represent the most general and restrictive subsets, respectively. Increasing m will progressively reduce the randomness at routing decisions. Here, we find a nested collection of subsets corresponding to the intermediate situations $`0<m<M`$. The relative size of each subset can be approximated by measuring the fraction of nodes F visited by packets:

$$F=\frac{1}{M}\underset{i=1}{\overset{M}{}}\theta \left[A_i\right]$$
(12)
where $`A_i`$ is the probability of visiting the node $`s_i`$. This term depends both on the network dynamics and the geometrical situation of the node within the network. Considerable effort has been devoted to characterize the likelihood of visiting a node in a purely static fashion. In this context, two relevant centrality (or ’load’ Noh Goh2003 ) measures are Random Walk Centrality Newman2003  and Betweenness Centrality Freeman77  which certainly represent the two extremes of our hierarchy of path subsets. Anyway, the correlation of these measurements with dynamic centrality $`A_i`$ is weak, to say the least.

In figure 4 the numerically measured $`F(m)`$ for a finite range of $`m`$ values is shown. This fraction equals 1 for random routing $`0m<[d]`$, where $`d`$ is the averaged shortest path length and $`[x]`$ denotes the integer part of x. Random routing does not discard any network node. A sharp transition takes place from this point $`[d]<mM`$ and a large fraction of network nodes (about 45 percent) are never visited by deterministic routing. This severly restricts the diversity of routes and yields a load-insensitive system Lorenz . The system requires a certain degree of noise in order to avoid sending packets through already collapsed nodes while it is enabled to choose less optimal but free routes. Note that depth of routing parameter $`m`$ defines an order parameter because quantifies the degree of order existing in the system. The existence of an order parameter confirms the critical nature of Internet traffic.

Numerical simulations have shown that flow is maximized at the order-disorder transition point $`m=[d]`$. This can be observed in network throughput reaching a maximum at this point (see figure 5). Throughput is defined by the quotient of the number of successfully released packets and the sum of all packets generated during the simulation. In the next section, we will show that several real statistics can be reproduced by the system dynamics poised to criticality.

## 6 Average stretch

The average stretch $`s`$ measures the efficiency of routing by comparing the number of hops $`h`$ traversed by a packet to the shortest path distance $`d`$ between source and destination:

$$s=\frac{h}{d}$$

Compact routing schemes minimize the average stretch while maintaining the size of routing tables small Cowen . Reducing the average strech will progressively raise up the memory requirements at each router. Assuming Thorup-Zwick (TZ) compact routing scheme TZ2001 , the average stretch can be expressed as a function of the distance distribution and the graph size $`M`$ only: $`s=f(d,\sigma _d)`$ Krioukov . This TZ scheme ensures a nearly optimal lower memory upper bound for $`s=3`$ in generic networks. For scale-free networks, TZ achieves lower bounds. In particular, for the Internet interdomain graph (Autonomous Systems) with degree distribution exponent $`\gamma 2.1`$ and $`M=10^4`$, the TZ average stretch $`s1.14`$Krioukov . Also, the average number of entries in the routing tables is approximately 52. It turns out that $`s(d,\sigma _d)`$ surface has unique minimums. Strikingly, the points corresponding to Internet distance distribution are very close to them Krioukov . This suggests that Internet topology is shaped by some hidden optimization criteria. Anyway, TZ scheme is not a realistic Internet interdomain routing scheme because assumes that global topology view is available. We have numerically measured the average stretch and the average fraction of neighbours at distance $`m`$ for our routing scheme (see figure 6). Note that the critical path horizon $`m=[d]`$ is very close to the minimum average stretch.

## 7 Network fluctuations and performance

In order to test the goodness of the model presented here, it is interesting to compare some real Internet statistics with their respective model measurements. In particular, those measurements will be collected for the model at the critical path horizon, that is, when $`m=[d]`$.

When understanding the competition between a network’s internal collective dynamics (i.e: Internet traffic) and external environmental changes (i.e: traffic sources or host behaviour), it is useful to study the relationship between the mean flux and the size of fluctuations around the average Fluctuations . Previous explanations of ocurrence of self-similarity in traffic networks are based on the superposition of many and high-variable (infinite variance) sources Willinger . This point of view discards the effects of the system’s collective dynamics and considers that Internet is an externally driven system. Real data shows that Internet dynamics can not be simply reduced to the behaviour of traffic sources.

Let us define the incoming (outgoing) flux $`f_i`$ as the amount of packets received (forwarded) at router $`s_i`$ during a given and fixed period of time. For every router, compare the average flux $`f_i`$ with the dispersion $`\sigma _i`$ around the mean. It has been noticed that for several real systems the following scaling relation holds:

$$\sigma f^\alpha $$

where $`\alpha `$ is an exponent which can take the values of 1/2 and 1 Fluctuations . This suggests that real systems can be classified in two main classes depending on the value of this exponent. The relevant exponent in this context is $`\alpha =1/2`$, which has been observed in daily traffic measurements at 374 geographically distinct Internet routers Fluctuations . Systems exhibiting the 1/2 exponent are representative of endogeneous dynamics, that is, determined by the system’s internal collective fluctuations. Moreover, measurements on our model reproduce the same exponent (see figure 7).

A measure of Internet end-to-end performance is the normalized latency time $`\tau `$ Vespignani  defined as the quotient between latency time $`L`$ (measured as Round-Trip-Time) and geographical distance $`w`$:

$$\tau =\frac{L}{w}$$

There are several factors governing $`\tau `$. First, propagation rate is finite. There is a minimum delay because packets can not move faster than speed of light. Second, the number of nodes traversed departs from the minimum found along the geodesic path from node to destination. And third, packets spent some time at every intermmediate node because queueing delays. Define $`\tau _{min}=L_{min}/w`$ as the normalized latency without taking into acount queue delays and $`\tau _{av}=L_{av}/w`$ as the normalized latency while considering all factors affecting packet latency.

Their probability distributions have been measured from two years of PingER PingER  data and follow power-law scaling with stable exponents of about -3.0 for $`\tau _{min}`$ and -2.5 for $`\tau _{av}`$ Vespignani . Numerical simulations on the model poised to criticality reproduced the power-law $`P(\tau _{av})`$ probability distribution. The exponent of the distribution is about -2.45, which is very close to real observations (see figure 8). The previous results seems to be quite robust and independent of most model parameters changes. An interesting thing to note is that current model cannot reproduce the existing correlation between $`\tau _{min}`$ and $`\tau _{av}`$ and geographical distance $`w`$ reported by previous studies Vespignani . The reason might be that propagation rate is (unrealistically) assumed to be infinite in the model. It might be that this factor has little influence in shaping the $`\tau _{av}`$ probability distribution. Moreover, it is worth noting that any deviation of the order parameter from the critical regime resulted in an exponential distribution for $`\tau _{av}`$, reinforcing the view of an optimally efficient and self-organized Internet.

## 8 Discussion

In this paper we have shown that a simple model of traffic dynamics incorporating the appropriate Internet’s network topology is able to recover several statistical features of real traffic. More important, we have seen that the routing algorithms can take advantage of the small-world structure of the web by reaching a critical path horizon close to the network’s average path length. In doing so, a highly efficient system is reached at low cost: routing strategies only need to consider a small depth. Once the $`m`$ parameter reaches the network’s diameter, no further information is required to properly reach the target. A full-system deterministic routing strategy is actually unnecesary and would be too costly. Instead, the constraints imposed by network’s architecture allow to exploit the implicit information defined by the small-world architecture.

It might be that Internet evolves in a way that throughput or global performance is maximized. This trend is constrained within the limits of available communication resources. Inside this regime, Internet is shaped in order to provide better response. This is reflected in the optimal and lower average stretch observed in real Internet, which is close to global minimum TZ2001 . This suggest hat some hidden optimization is at work. Clearly, an unresponsive system will be no good. How a distributed collection of designers were able to define this globally efficient infrastructure is a question that deserves attention.

###### Acknowledgements.

This paper is dedicated to the memory of Per Bak. This work was supported by a grant BFM2001-2154 and by the Santa Fe Institute.


text:
# The problem with the Brier score

## 1 Introduction

Users of meteorological forecasts need to be able to judge which forecasts are the best in order to decide which to use. We distinguish two cases. The first case is one in which the user plans to use the forecast for making a certain specific decision the details of which can be specified entirely in advance. The second is one in which the user plans to use the forecast for making one or more decisions which cannot be specified in detail in advance.

In the first case it may be possible to decide which forecast is the best by analysing the effect of using different forecasts on the quality of the final decisions made (for an example of this situation see Richardson (2000)). In the second case, however, the user cannot convert forecasts into decisions ahead of time because they do not know what decisions they are going to have to make. By the time they know what decision they are going to have to make, they do not have time to re-evaluate the available forecasts and potentially switch to a different forecast provider. In this second case forecasts have to be analysed and compared on their own merits, rather than on the merits of the decisions that can be based on them. In such a situation, the forecast user needs standard measures which can distinguish between forecasts at a general level222A good example is the root mean square error, which is a general measure used for comparing forecasts for the expectation. It is this second case that we will consider.

Forecasts can be divided into forecasts of the expectation of future outcomes and probabilistic forecasts that give probabilities of different outcomes. Probabilistic forecasts can then be divided into continuous and discrete probabilistic forecasts. A continuous probabilistic forecast gives a continuous density for the distribution of possible outcomes. We have discussed how to measure the skill of such forecasts in Jewson (2003b) and have applied the measures we propose to the calibration and the comparison of forecasts in a number of studies such as Jewson et al. (2003) and Jewson (2003a).

Discrete probabilistic forecasts give probabilities for a number of discrete events. Any number of events can be considered, but in this article we will restrict ourselves to the case of only two events, which we call a binary probabilistic forecast. We will address the question of how binary probabilistic forecasts can be compared.

One of the standard tools used by meteorologists to answer the question of which of two binary probabilistic forecasts is the better is the Brier score, first used over 50 years ago (Brier, 1950), and still in use today (see, for example, Vitart (2003), page 25). Nevertheless, we are going to argue that the Brier score is flawed. This is not something that can be proven *mathematically*, of course. Our arguments will be based on an appeal to intuition: we will present a simple case in which we believe it is intuitively clear which of two forecasts is the better, and we will show that the Brier score then comes to the opposite conclusion to our intuition i.e. it gives the *wrong answer*. We will then present an alternative score that overcomes this problem, that has a definition that accords more clearly with intuition, and that is also more firmly grounded in standard statistical theory.

## 2 The Brier Score

The Brier score for a binary event is defined as:

$$b=<(fo)^2>$$
(1)
where $`f`$ is a forecast for the probability that an event $`X`$ will happen, and $`o`$ is an observation which takes the value 1 if the event happens and 0 otherwise. Lower values of the Brier score indicate better forecasts. A detailed discussion is given in Toth et al. (2003).

We can expand the Brier score as:

$$b=<f^2>2<fo>+<o^2>$$
(2)
When we are comparing two forecasts on the same observed data set the difference in the Brier score is given by:

$$b_2b_1=<f_2^2>2<f_2o><f_1^2>+2<f_1o>$$
(3)
where the $`<o^2>`$ term has cancelled because it is the same for both forecasts. If this difference is positive ($`b_2>b_1`$) then we conclude that $`b_1`$ is the better forecast.

A particularly simple case is where the forecast probabilities have constant values, giving:

$$b_2b_1=f_2^22f_2<o>f_1^2+2f_1<o>$$
(4)
A further simplification is possible if the event occurs with a constant probability $`p`$, in which case $`<o>=p`$ and

$$b_2b_1=f_2^22f_2pf_1^2+2f_1p$$
(5)
## 3 A simple example

We now consider a very simple example, with constant event probability and constant forecast probabilities. We set $`p=\frac{1}{10}`$, and consider the forecasts $`f_1=0`$ and $`f_2=\frac{1}{4}`$.

In this case the difference between the Brier scores is given by:

$`b_2b_1`$ $`=`$ $`f_2^22f_2pf_1^2+2f_1p`$
$`=`$ $`\left({\displaystyle \frac{1}{4}}\right)^22.{\displaystyle \frac{1}{4}}.{\displaystyle \frac{1}{10}}`$
$`=`$ $`{\displaystyle \frac{1}{16}}{\displaystyle \frac{1}{20}}`$
$`=`$ $`{\displaystyle \frac{1}{80}}`$
The Brier score leads us to conclude that forecast $`f_1`$ is the better forecast. However, this does not agree with our intuition. Forecast $`f_1`$ is a disaster: it predicts a zero probability (a very strong statement!) for something that happens not infrequently. Forecast $`f_1`$ is completely invalidated whenever event $`X`$ actually occurs (on average, 1 in every 10 trials). Forecast $`f_2`$, on the other hand, is not so bad. It gives a lowish probability for something that does indeed occur with a low probability. Its only fault is that the probability is not exactly correct.

The reason that the Brier score makes this mistake is that it does not penalise forecasts that predict a zero probability strongly enough when they are wrong, even though our intuition tells us that they should be heavily penalised. More generally, the Brier score does not penalise forecasts that give very small probabilities when they should be giving larger probabilities to the same extent that we penalise such forecasts with our intuition. This is because the Brier score is based on a straight difference between $`f`$ and $`o`$. Our intuition, on the other hand, considers the difference between probabilities of 0% and 10% to be very different from the difference between probabilities of 40% and 50%. Intuition apparently uses fractional or logarithmic rather than absolute differences in probability.

One can easily construct other examples that illustrate this point. The more extreme the events considered, the more striking is the problem with the Brier score. Consider, for example, $`p=\frac{1}{1000}`$, $`f_1=0`$ and $`f_2=\frac{1}{400}`$. Again the Brier score prefers $`f_1`$, while our intuition considers $`f_1`$ to be a failure, and $`f_2`$ to be a reasonably good attempt at estimating a very small probability.

We conclude that the Brier score cannot be trusted to make the right decision about which of two forecasts is better. It should also not be used to calibrate forecasts or evaluate forecasting systems since it will over-encourage prediction of very small or zero probabilities. We need a different measure.

## 4 The likelihood score

The standard measure used in classical statistics for testing which of two distributions gives the best fit to data is the likelihood $`L`$ defined as the probability (or probability density) of the observations given the model and the parameters of the model (Fisher, 1922). In our case this becomes the probability of the observations given the forecast.

We advocate the likelihood as the best metric for calibrating and comparing *continuous* probabilistic forecasts (see the previous citations) mainly on the basis that it is very intuitively reasonable: the forecast that gives the highest probability for the observations is the better forecast. We also advocate the likelihood as the best metric for calibrating and comparing binary forecasts. In this case the likelihood is given by:

$$L=p(x|f)$$
(7)
where $`x`$ is the full set of observations and $`f`$ is the full set of forecasts. If we assume that the forecast errors are independent in time then this becomes:

$`L`$ $`=`$ $`\mathrm{\Pi }_{i=1}^{i=n}p(x_i|f_i)`$
$`=`$ $`\mathrm{\Pi }_{i=1}^{i=n}o_if_i+(1o_i)(1f_i)`$
We can also use the log-likelihood, which gives a more compressed range of values, and is given by:

$`l`$ $`=`$ $`lnL`$
$`=`$ $`ln[\mathrm{\Pi }_{i=1}^{i=n}o_if_i+(1o_i)(1f_i)]`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{i=n}{}}}ln[o_if_i+(1o_i)(1f_i)]`$
If we put all cases of event $`X`$ occuring into set $`A`$, and all cases of event $`X`$ not occuring into set $`B`$ then:

$$L=\mathrm{\Pi }_Af_i\mathrm{\Pi }_B(1f_i)$$
(10)
and

$$l=\underset{A}{}f_i+\underset{B}{}(1f_i)$$
(11)
If we now consider the special case in which $`f`$ is constant then:

$$L=f_i^a(1f_i)^b$$
(12)
and

$$l=alnf_i+bln(1f_i)$$
(13)
where $`a`$ is the number of occurences of $`X`$, $`b`$ is the number of occurences of not $`X`$, and $`b=na`$.

If any of the predictions $`f`$ are 0 or 1 (i.e. are completely certain) then $`L=0`$ and $`l=\mathrm{}`$. If not, then $`L>0`$ and $`l>\mathrm{}`$. We see that use of the likelihood penalises the use of probability forecasts with values of 0 or 1 very heavily. Such forecasts get the worst possible score, as they should (since one can *never* be completely certain).

In our simple example the difference in likelihoods for the two forecasts is:

$$L_2L_1=f_2^a(1f_2)^b$$
(14)
Since this is positive for all samples we see that the likelihood concludes that forecast 2 is better, in line with our intuition.

## 5 Summary

Meteorologists have used the Brier score to compare binary probabilistic forecasts for over 50 years. However, we find that in simple cases it makes the wrong decision as to which is the better of two forecasts (where we define *wrong* in terms of our intuition). We reach this conclusion independently of any detailed analysis of the preferences of the user of the forecast.

We advocate scores based on the likelihood as a replacement for the Brier score. On the one hand the likelihood is conceptually simpler than the Brier score: it decides which forecast is better simply according to which forecast gives the higher probability for the observed data, which seems immediately reasonable. On the other hand the likelihood accords with our intuition in the simple example that we present, and punishes forecasts that give probabilities of 0 and 1 appropriately.

We conclude that use of the Brier score should be discontinued, and should be replaced by a score based on the likelihood.

## 6 Acknowledgements

The author would like to thank Rie Kondo and Christine Ziehmann for useful discussions, and Christine for reading the manuscript and making some helpful comments.

## 7 Legal statement

The author was employed by RMS at the time that this article was written.

However, neither the research behind this article nor the writing of this article were in the course of his employment, (where ’in the course of his employment’ is within the meaning of the Copyright, Designs and Patents Act 1988, Section 11), nor were they in the course of his normal duties, or in the course of duties falling outside his normal duties but specifically assigned to him (where ’in the course of his normal duties’ and ’in the course of duties falling outside his normal duties’ are within the meanings of the Patents Act 1977, Section 39). Furthermore the article does not contain any proprietary information or trade secrets of RMS. As a result, the author is the owner of all the intellectual property rights (including, but not limited to, copyright, moral rights, design rights and rights to inventions) associated with and arising from this article. The author reserves all these rights. No-one may reproduce, store or transmit, in any form or by any means, any part of this article without the author’s prior written permission. The moral rights of the author have been asserted.


text:
# Why magnesium diboride is not described by anisotropic Ginzburg-Landau theory

## Abstract

It is well established that the superconductivity in the recently discovered superconducting compound MgB2 resides in the quasi-two-dimensional band ($`\sigma `$-band) and three-dimensional band ($`\pi `$-band). We demonstrate that, due to such band structure, the anisotropic Ginzburg-Landau theory practically does not have region of applicability, because gradient expansion in the $`c`$ direction breaks down. In the case of dirty $`\pi `$-band we derive the simplest equations which describe properties of such superconductors near $`T_c`$ and explore some consequences of these equations.

Ginzburg-Landau (GL) theory is the most powerful and widely used phenomenological theory of superconductivity (see, e. g., Refs. GL, ; Tinkham, ). It describes practically all known superconductors in the vicinity of transition temperature. GL theory is fully microscopically justified and all its parameters can be derived from the microscopic BCS theory Gorkov . GL theory provides the basis for such elaborated fields as vortex physics Abrikosov  and the theory of fluctuation phenomena Tinkham .

The recently discovered superconductor MgB2 Akimitsu  gives an example of a superconductor which *is not described by the anisotropic GL theory*. This very unusual feature is a consequence of a specific band structure of this compound. It is reliably established that superconductivity in MgB2 resides in two families of bands: strongly superconducting quasi-two-dimensional $`\sigma `$-bands and weakly superconducting three-dimensional $`\pi `$-bands (see, e.g., Ref. PhysCSpecial ). Both bands are characterized by their intrinsic coherence lengths, and the $`c`$ axis coherence length in the $`\sigma `$-band is much smaller than $`c`$ axis coherence length in the $`\pi `$-band. Typically, the strong band forces the order parameter in the weak band to change in the $`c`$ direction at distances smaller than the intrinsic $`c`$ axis coherence length in this band. This means that almost in the whole temperature range the effective coherence length in the $`c`$ direction, $`\xi _z(T)`$, is smaller than the intrinsic coherence length in the $`\pi `$-band, $`\xi _{\pi ,z}`$. The crossover to the GL region takes place only when $`\xi _z(T)`$ exceeds $`\xi _{\pi ,z}`$, which occurs in the very close vicinity of $`T_c`$. In this narrow region the $`\pi `$-band strongly increases the $`c`$ axis coherence length. Beyond the narrow region, the variations of the order parameter in the $`c`$ directions are not described by the anisotropic GL theory. Important consequencies of GL theory breakdown are the strong temperature dependence of the $`H_{c2}`$ anisotropy Hc2Clean ; Hc2Exp ; GurevichPRB03 ; GK-Hc2  and strong deviations of the $`H_{c2}`$ angular dependence from the simple ”effective mass” law GK-Hc2 ; Rydh03 .

Obviously, the breakdown of the anisotropic GL theory has numerous consequences and it would be desirable (i) to trace the reason of this breakdown and (ii) to derive the simplest model, which replaces the GL model near $`T_c`$. This paper addresses these issues. For illustration, we use the simplest microscopic model, multiband generalization of the Usadel theory, describing a dirty two-band superconductor with weak interband scattering KG ; GurevichPRB03 . However, the main conclusions are very general and do not depend much on the intraband scattering strength. In the model we use the GL expansion for the $`\sigma `$-band and keep the microscopic description for the $`\pi `$-band, i.e., only ”dirtyness” of the $`\pi `$-band is essential for a particular form of equations.

We consider a dirty two-band superconductor with weak interband scattering. Such superconductor is described by Usadel equations for the impurity averaged normal and anomalous Green’s functions, $`G_\alpha `$ and $`F_\alpha `$, $`G_\alpha ^2+\left|F_\alpha \right|^2=1`$, and the pair potentials $`\mathrm{\Delta }_\alpha `$,

$$\omega F_\alpha _j\frac{𝒟_{\alpha ,j}}{2}\left[G_\alpha D_j^2F_\alpha F_\alpha _j^2G_\alpha \right]=\mathrm{\Delta }_\alpha G_\alpha $$
(1)
where $`\alpha =1,2`$ is the band index, $`j=x,y,z`$ is the coordinate index, $`D_j_j(2\pi i/\mathrm{\Phi }_0)A_j`$, $`𝒟_{\alpha ,j}`$ are diffusion constants, and $`\omega =2\pi T(s+1/2)`$ are Matsubara frequencies. Bearing in mind the application to MgB2, in our notations index 1 corresponds to $`\sigma `$-bands and index 2 to $`\pi `$-bands, $`𝒟_{1,j}𝒟_{\sigma ,j}`$ and $`𝒟_{2,j}𝒟_{\pi ,j}`$. All bands are isotropic in the $`xy`$ plane, $`𝒟_{\alpha ,x}=𝒟_{\alpha ,y}`$ and anisotropic in the $`xz`$ plane with the anisotropy ratios $`\gamma _\alpha =\sqrt{𝒟_{\alpha ,x}/𝒟_{\alpha ,z}}`$. Self-consistency conditions can be written as KG

$`W_1\mathrm{\Delta }_1W_{12}\mathrm{\Delta }_2`$ $`=2\pi T{\displaystyle \underset{\omega >0}{}}\left(F_1{\displaystyle \frac{\mathrm{\Delta }_1}{\omega }}\right)+\mathrm{\Delta }_1\mathrm{ln}{\displaystyle \frac{1}{t}}`$ (2a)
$`W_{21}\mathrm{\Delta }_1+W_2\mathrm{\Delta }_2`$ $`=2\pi T{\displaystyle \underset{\omega >0}{}}\left(F_2{\displaystyle \frac{\mathrm{\Delta }_2}{\omega }}\right)+\mathrm{\Delta }_2\mathrm{ln}{\displaystyle \frac{1}{t}}`$ (2b)
where $`tT/T_c`$ and the matrix $`W_{\alpha \beta }`$ is related to the matrix of coupling constants $`\mathrm{\Lambda }_{\alpha \beta }`$ as

$`W_1`$ $`={\displaystyle \frac{\mathrm{\Lambda }_{}+\sqrt{\mathrm{\Lambda }_{}^2+\mathrm{\Lambda }_{12}\mathrm{\Lambda }_{21}}}{\mathrm{Det}}},W_2={\displaystyle \frac{\mathrm{\Lambda }_{}+\sqrt{\mathrm{\Lambda }_{}^2+\mathrm{\Lambda }_{12}\mathrm{\Lambda }_{21}}}{\mathrm{Det}}},`$
$`W_{12}`$ $`=\mathrm{\Lambda }_{12}/\mathrm{Det},W_{21}=\mathrm{\Lambda }_{21}/\mathrm{Det},`$
$`\mathrm{\Lambda }_{}(\mathrm{\Lambda }_{11}\mathrm{\Lambda }_{22})/2,`$ $`\mathrm{Det}\mathrm{\Lambda }_{11}\mathrm{\Lambda }_{22}\mathrm{\Lambda }_{12}\mathrm{\Lambda }_{21}`$, $`W_1W_2=W_{12}W_{21}`$. The supercurrent components are given by

$$j_j=4\pi eT\underset{\alpha }{}\underset{\omega >0}{}N_\alpha 𝒟_{\alpha ,j}\mathrm{Im}\left[F_\alpha ^{}D_jF_\alpha \right],$$
(3)
where $`N_\alpha `$ are the partial densities of states

We start with the derivation of the GL equations from the Usadel equations in the close vicinity of $`T_c`$ following a standard route. In the lowest approximation $`G_\alpha ^{(0)}=1`$ and $`F_\alpha ^{(0)}\mathrm{\Delta }_\alpha /\omega `$. When $`\mathrm{\Delta }_\alpha `$ are small and change slowly in space (the exact criterion will be derived below) one can keep only the leading nonlinear and gradient corrections

$$F_\alpha \frac{\mathrm{\Delta }_\alpha }{\omega }\frac{\mathrm{\Delta }_\alpha ^3}{2\omega ^3}+_j\frac{𝒟_{\alpha ,j}}{2\omega ^2}D_j^2\mathrm{\Delta }_\alpha .$$
(4)
Substituting this expansion into the self-consistency conditions (2), we obtain coupled GL equations for two gap parameterscleanGL

$`W_1\mathrm{\Delta }_1W_{12}\mathrm{\Delta }_2`$ $`={\displaystyle _j}\xi _{1,j}^2D_j^2\mathrm{\Delta }_1b\mathrm{\Delta }_1^3+\tau \mathrm{\Delta }_1`$ (5a)
$`W_{21}\mathrm{\Delta }_1+W_2\mathrm{\Delta }_2`$ $`={\displaystyle _j}\xi _{2,j}^2D_j^2\mathrm{\Delta }_2b\mathrm{\Delta }_2^3+\tau \mathrm{\Delta }_2`$ (5b)
with $`\xi _{\alpha ,j}^2=(\pi /8T)𝒟_{\alpha ,j}`$, $`b=7\zeta (3)/(8\pi ^2T^2)`$, and $`\tau =\mathrm{ln}(1/t)(T_cT)/T_c`$.
Near $`T_c`$ the right hand sides of Eqs. (5a,b) are small. This allows us to reduce Eqs. (5a,b) to a single GL equation by looking for solution for $`\mathrm{\Delta }_2`$ in the form

$$\mathrm{\Delta }_2\frac{W_{21}}{W_2}\mathrm{\Delta }_1+\delta _2\text{,}$$
(6)
From Eqs. (5a,b) we obtain

$`W_{12}\delta _2`$ $`={\displaystyle _j}\xi _{1,j}^2D_j^2\mathrm{\Delta }_1b\mathrm{\Delta }_1^3+\tau \mathrm{\Delta }_1`$ (7a)
$`W_2\delta _2`$ $`={\displaystyle _j}\xi _{2,j}^2D_j^2\mathrm{\Delta }_2b\mathrm{\Delta }_2^3+\tau \mathrm{\Delta }_2.`$ (7b)
The second equation indicates that $`\delta _2`$ is a small correction, $`\delta _2\mathrm{\Delta }_2`$, and one can use $`\mathrm{\Delta }_2\left(W_{21}/W_2\right)\mathrm{\Delta }_1`$ in the right hand side of this equation. Excluding $`\delta _2`$ and introducing the band-averaged order parameter
$$\mathrm{\Delta }^2=\frac{W_2\mathrm{\Delta }_1^2+W_1\mathrm{\Delta }_2^2}{W_2+W_1}\frac{W_{12}W_2^2+W_{21}W_1^2}{W_{12}W_2\left(W_2+W_1\right)}\mathrm{\Delta }_1^2$$

we finally obtain the anisotropic GL equation for $`\mathrm{\Delta }`$

$$_j\xi _j^2D_j^2\mathrm{\Delta }+b\mathrm{\Delta }^3\tau \mathrm{\Delta }=0,$$
(8)
with the average coherence lengths

$$\xi _j^2=\frac{W_2\xi _{1,j}^2+W_1\xi _{2,j}^2}{W_2+W_1}.$$

For the supercurrent, using relation $`W_{21}/W_{12}=\mathrm{\Lambda }_{21}/\mathrm{\Lambda }_{12}=N_1/N_2`$, we derive

$$j_j4eNP\xi _j^2\mathrm{Im}\left[\mathrm{\Delta }^{}D_j\mathrm{\Delta }\right]$$
(9)
with $`NN_1+N_2`$ and

$$P\frac{N_1N_2\left(W_2+W_1\right)^2}{\left(N_1+N_2\right)\left(N_2W_2^2+N_1W_1^2\right)}.$$

From Eqs. (8) and (9) we derive the components of the London penetration depth

$$\lambda _j^2\frac{32\pi ^2eNP\xi _j^2\tau }{c\mathrm{\Phi }_0b}.$$

For the parameters of MgB2, $`W_1W_2`$, $`\xi _{1,z}\xi _{2,z}`$, $`\xi _{1,x}\xi _{2,x}`$, the dominating effect of the $`\pi `$-band is the renormalization of the $`c`$ axis lengths

$`\xi _z^2`$ $`\xi _{1,z}^2+S_{12}\xi _{2,z}^2,`$ (10)
$`\lambda _z^2`$ $`{\displaystyle \frac{32\pi ^2e\tau N_1}{c\mathrm{\Phi }_0b}}\left(\xi _{1,z}^2+S_{12}\xi _{2,z}^2\right).`$ (11)
with $`S_{12}W_1/W_21`$. The influence of the $`\pi `$-band on properties not related with the variations of the order parameter along the $`c`$ axis are weak and can be treated perturbatively.

We obtain now the applicability criterion for the GL expansion. The gradient expansion is justified if $`\xi _{\alpha ,j}^2_j^2\mathrm{\Delta }_\alpha <\mathrm{\Delta }_\alpha `$ for all $`\alpha `$ and $`i`$. Because a typical scale of the spatial variations is the temperature-dependent GL coherence length $`\xi _j(T)`$, this condition simply means

$$\xi _j(T)>\xi _{\alpha ,j}$$
(12)
The most restraining inequality is the one for $`\alpha =2`$ and $`i=z`$ which gives

$$(T_cT)/T_c<\xi _{1,z}^2/\xi _{2,z}^2+S_{12}$$
(13)
Because $`\xi _{1,z}\xi _{2,z}`$ and $`S_{12}1`$, the applicability of the GL approach is limited to extremely narrow temperature range near $`T_c`$, i.e., the situation is very different from conventional superconductors. For parameters of MgB2 this condition implies $`\left(T_cT\right)/T_c0.05`$. On the other hand, near $`T_c`$ the fluctuation effects become important. This means that the mean-field GL theory practically does not have a region of applicability.

We derive now the simplest theory which replaces the GL theory in the conventional GL region $`\left(T_cT\right)/T_c<1`$. As the gradient expansion actually breaks down only for the $`\pi `$-band, in the vicinity of $`T_c`$ we can proceed with the expansion (4) for the $`\sigma `$-band, $`\alpha =1`$. Substituting this expansion into the self-consistency conditions, we obtain Eq. (5a). The $`\pi `$-band only weakly renormalizes the nonlinear term and we can use the linear approximation in this band

$$\omega F_2_j\frac{𝒟_{2,j}}{2}D_j^2F_2=\mathrm{\Delta }_2.$$
(14)
The $`\pi `$-band order parameter can again be represented by Eq. (6) with $`\delta _2`$ being a small correction. Finding this correction from Eq. (2b) and substituting it into Eq. (2a) together with the GL expansion for $`F_1`$, we derive coupled equations for $`\mathrm{\Delta }_1`$ and reduced $`\pi `$-band $`F`$-function $`f_s`$, $`f_s(2\pi TW_2/W_{21})F_2(\omega _s)`$criterion

$`\left(1+S_{12}\right)\tau \mathrm{\Delta }_1+b\mathrm{\Delta }_1^3{\displaystyle _j}\xi _{1,j}^2D_j^2\mathrm{\Delta }_1`$
$`S_{12}{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}\left(f_s{\displaystyle \frac{\mathrm{\Delta }_1}{s+1/2}}\right)`$ $`=0,`$ (15a)
$`(s+1/2)f_s{\displaystyle \frac{2}{\pi ^2}}{\displaystyle _j}\xi _{2,j}^2D_j^2f_s`$ $`=\mathrm{\Delta }_1,`$ (15b)
and the expression for the supercurrent

$$j_j4eN_1\xi _{1,j}^2\mathrm{Im}\left[\mathrm{\Delta }_1^{}D_j\mathrm{\Delta }_1\right]+\frac{8e}{\pi ^2}N_1S_{12}\xi _{2,j}^2\underset{s=0}{\overset{\mathrm{}}{}}\mathrm{Im}\left[f_s^{}D_jf_s\right].$$

These equations replace the GL equations in the case of dirty $`\pi `$-band. Note that the same equations are also valid in the case of clean $`\sigma `$-band but with different definition of the coherence length $`\xi _{1,j}`$, $`\xi _{1,j}^2=7\zeta (3)v_{1,j}^2/(4\pi T)^2`$.

In the case of weak superconductivity in the $`\pi `$-band, $`S_{12}1`$, and for $`\xi _{1,z}`$ $`\xi _{2,z}`$, one can neglect the in-plane gradients in the Eq. (15b) and obtain an even simpler set of equations which describe only *the dominating strong effects*, related to inhomogeneities of the gap parameter along the $`c`$ axis, and neglect small renormalizations of the coefficients by the weak $`\pi `$-band

$`\tau \mathrm{\Delta }_1`$ $`+b\mathrm{\Delta }_1^3{\displaystyle _j}\xi _{1,j}^2D_j^2\mathrm{\Delta }_1`$
$`S_{12}{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}\left(f_s{\displaystyle \frac{\mathrm{\Delta }_1}{s+1/2}}\right)=0,`$ (16a)
$`\left(s+1/2\right)f_s{\displaystyle \frac{2}{\pi ^2}}\xi _{2,z}^2D_z^2f_s=\mathrm{\Delta }_1,`$ (16b)
$`j_j`$ $`4eN_1\xi _{1,j}^2\mathrm{Im}\left[\mathrm{\Delta }_1^{}D_j\mathrm{\Delta }_1\right]`$
$`+\delta _{j,z}{\displaystyle \frac{8e}{\pi ^2}}N_1S_{12}\xi _{2,z}^2{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}\mathrm{Im}\left[f_s^{}D_zf_s\right],`$ (16c)
We explore now some consequences of these equations. To define an effective coherence length, we consider the response of the order parameter to the weak $`z`$-dependent variation of $`T_c`$, $`\tau \tau (z)=\tau +\delta \tau (z)`$. In linear approximation with respect to $`\delta \tau (z)`$ Eqs. (16a) and (16b) can be solved by Fourier transform yielding $`\mathrm{\Delta }_1=\mathrm{\Delta }_1^{(0)}+\delta \mathrm{\Delta }_1(z)`$

$`\delta \mathrm{\Delta }_1(z)`$ $`={\displaystyle G(zz^{})\delta \tau (z^{})𝑑z^{}},`$
$`G(z)`$ $`={\displaystyle \frac{dk}{2\pi }\frac{\mathrm{exp}(ikz)}{2\tau +\xi _{1,z}^2k^2+S_{12}g[(2/\pi ^2)\xi _{2,z}^2k^2]}},`$
where $`g(u)\psi (1/2+u)\psi (1/2)`$ and $`\psi (u)`$ is the digamma-function. In contrast to the GL model, the decay of the perturbation $`\delta \mathrm{\Delta }_1(z)`$ is not exponential. Using the last equation, one can introduce the effective coherence length $`\xi _z`$, which determines the scale of spatial variations of the order parameter in the $`z`$ direction,

$$\xi _{1,z}^2/\xi _z^2+S_{12}g\left[(2/\pi ^2)\xi _{2,z}^2/\xi _z^2\right]=\tau $$
(17)
The dependence $`\xi _z(T)`$ computed from this equation using parameters $`S_{12}=0.034`$ KG  and $`\xi _{2,z}^2=300\xi _{1,z}^2`$ is shown in Fig. 1.

Consider the relation between the supercurrent $`j_z`$ and supermomentum $`p_z=_z\varphi (2\pi /\mathrm{\Phi }_0)A_z`$, which determines the $`c`$ axis London length and depairing current. From Eqs. (16a) and (16b) we obtain

$`j_z(p_z)`$ $`=4eN_1\mathrm{\Delta }_1^2(p_z)p_z`$ (18)
$`\times \left(\xi _{1,z}^2+{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(2/\pi ^2)S_{12}\xi _{2,z}^2}{\left(s+1/2+(2/\pi ^2)(\xi _{2,z}p_z)^2\right)^2}}\right)`$
$`\mathrm{\Delta }_1^2(p_z)`$ $`=\left(\tau \xi _{1,z}^2p_z^2S_{12}g\left[(2/\pi ^2)\xi _{2,z}^2p_z^2\right]\right)/b.`$
In the linear regime $`j_z(4e\tau N_1/b)\left(\xi _{1,z}^2+S_{12}\xi _{2,z}^2\right)p_z`$. This means that in the whole range $`(T_cT)/T_c1`$ the $`z`$ component of the London length is given by the GL formula (11). In conventional superconductors the dependence $`j_z(p_z)`$ is nonmonotonic and its maximum gives the well-known GL result for depairing current, $`j_{dp}=c\mathrm{\Phi }_0/(12\sqrt{3}\pi ^2\lambda ^2\xi )\tau ^{3/2}`$ for $`\tau 0`$.Tinkham  In our case the situation is different. The amplitude of the order parameter is suppressed at $`p_z1/\xi _z(T)`$. However, in the region $`\xi _z(T)\xi _{2,z}`$ the dependence $`j_z\left(p_z\right)`$ becomes nonlinear at much smaller $`p_z`$, $`p_z1/\xi _{2,z}`$. The shape of this dependence is determined by the parameter $`S_r=S_{12}\xi _{2,z}^2/\xi _{1,z}^2`$. The dependencies $`j_z(p_z)`$ for $`S_r=6`$ and different temperatures are plotted in the left panel of Fig. 2. For large values of $`S_r`$ the dependence $`j_z(p_z)`$ has *two maxima* within some temperature range, where first (second) maximum corresponds to the suppression of supercurrent in the $`\pi `$-($`\sigma `$-)band. At low temperatures the depairing current $`j_{dp}`$ is given by the second maximum and is determined mainly by the $`\sigma `$-band. At certain temperature near $`T_c`$ global maximum switches to the first maximum (see Fig. 2). The temperature dependence of $`j_{dp}`$ has a kink at this temperature (see the right panel of Fig. 2). For $`S_r>6`$ the local maximum of $`j_z\left(p_z\right)`$ at $`p_z1/\xi _{2,z}`$ exists even in the limit $`\xi _z(T)\xi _{2,z}`$.

As another example, we compute from Eqs. (16) the in-plane upper critical field near $`T_c`$, $`H_{c2,a}(T)`$. Experiment Hc2Exp  shows strong upward curvature of $`H_{c2,a}(T)`$, leading to the temperature dependent anisotropy factor. Microscopic calculations reproduce this feature, both in clean Hc2Clean  and dirty GK-Hc2  cases, but require rather heavy numerical computations. Our model allows to trace origin of the upward curvature in a simple way. Selecting the gauge $`A_z=Hx`$ and introducing reduced variables $`h=H/H_{c2}^{(1)}`$ with $`H_{c2}^{(1)}\mathrm{\Phi }_0/(2\pi \xi _{1,x}\xi _{1,z})`$, $`x\sqrt{h}x/\xi _{1,x}`$, $`r_z=𝒟_{2,z}/𝒟_{1,z}`$, we write the linear equation for determination of the upper critical field, $`h=H_{c2}/H_{c2}^{(1)}`$, as

$`{\displaystyle \frac{S_{12}}{h}}{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\mathrm{\Delta }_1}{s+1/2}}f_s\right)_x^2\mathrm{\Delta }_1+x^2\mathrm{\Delta }_1`$ $`={\displaystyle \frac{\tau }{h}}\mathrm{\Delta }_1,`$ (19a)
$`(s+1/2)f_s+{\displaystyle \frac{2}{\pi ^2}}r_zhx^2f_s`$ $`=\mathrm{\Delta }_1.`$ (19b)
Excluding $`f_s`$, we obtain the Schroedinger equation for $`\mathrm{\Delta }_1`$ with non-parabolic potential

$$_x^2\mathrm{\Delta }_1+\left(x^2+\frac{S_{12}}{h}g\left[\frac{2r_zhx^2}{\pi ^2}\right]\right)\mathrm{\Delta }_1=\frac{\tau }{h}\mathrm{\Delta }_1$$
(20)
Only in the limit $`h\sqrt{1+S_{12}r_z}/r_z1`$ this equation reduces to the usual oscillator equation. In this limit, using expansion $`g(u)(\pi ^2/2)u`$, we reproduce the GL result, $`h_{c2}=\tau /\sqrt{1+S_{12}r_z}`$. The inequality $`h_{c2}\sqrt{1+S_{12}r_z}/r_z`$ reproduces criterion (13) for the validity of the GL theory. In the opposite limit, $`\sqrt{1+S_{12}r_z}/r_zh1`$, one can use the asymptotics $`g(u)\mathrm{ln}\left(u\right)+\gamma _E+2\mathrm{ln}2`$ for $`u1`$ with $`\gamma _E0.577`$ being the Euler constant, and obtain

$`_x^2\mathrm{\Delta }_1+\left(x^2+{\displaystyle \frac{S_{12}}{h}}\mathrm{ln}(x^2)\right)\mathrm{\Delta }_1=\alpha (h)\mathrm{\Delta }_1`$
$`\tau =\alpha (h)h+S_{12}\left(\mathrm{ln}\left[{\displaystyle \frac{8hr_z}{\pi ^2}}\right]+\gamma _E\right)`$
This gives the following equation for the upper critical field

$$h_{c2}+S_{12}\mathrm{ln}\left[Ch_{c2}r_z\right]=\tau $$

with $`C1`$ ($`C(8/\pi ^2)\mathrm{exp}\left[\mathrm{ln}\left(x^2\right)+\gamma _E\right]=2/\pi ^2`$ for $`hS_{12}`$). In this limit the $`\pi `$-band gives only small logarithmic correction to the upper critical field. As we can see, the upper critical field has a strong upward curvature in a narrow region neat $`T_c`$: the slope $`dh_{c2}/d\tau `$ changes from $`1/\sqrt{1+S_{12}r_z}`$ to 1 near $`\tau =S_{12}+1/r_z`$, in agreement with microscopic calculations and experiment.

In conclusion, we demonstrated that the properties of magnesium diboride are not described by the anisotropic GL theory. We derived a simple model which replaces this theory in the vicinity of $`T_c`$, and explored some consequences of this model.

We would like to thank U. Welp for critical reading of the manuscript. This work was supported by the U.S. DOE, Office of Science, under contract # W-31-109-ENG-38.