text:
# Does the Boltzmann principle need a dynamical correction?

## I Introduction

In the notation of Planck, the connection between the thermodynamic entropy $`S`$ of an isolated system and its number of accessible states $`W`$ is often introduced as what Einstein called the Boltzmann principle Cercignani (1998)

$$S=k\mathrm{ln}W.$$
(1)
For classical systems, where the number of states is replaced by continuum integrals over the phase space, two common definitions of entropy for isolated systems are111The prefactors $`1/N!`$ and $`1/h^{3N}`$ are omitted since they are irrelevant in the present context. (see e.g. Huang (1987), or Pearson et al. (1985) and references therein)

$`S_\mathrm{\Phi }`$ $`=k\mathrm{ln}\mathrm{\Phi }(E),`$ (2)
$`S_\mathrm{\Omega }`$ $`=k\mathrm{ln}\mathrm{\Omega }(E),`$ (3)
where

$$\mathrm{\Phi }(E)=_{H(𝐳)<E}𝑑𝐳=𝑑𝐳\theta [EH(𝐳)]$$
(4)
is the volume under the constant-energy surface, and

$$\mathrm{\Omega }(E)=_{H(𝐳)=E}𝑑𝐳=𝑑𝐳\delta [EH(𝐳)]$$
(5)
is the “area” of that surface. Here $`H(𝐳)`$ is the Hamiltonian, $`𝐳`$ denotes a point in the phase space, $`\theta (x)`$ is the Heaviside step function and $`\delta (x)`$ is the Dirac delta function. These two expressions for $`S`$ will be referred to as the volume and area entropies, respectively. It is generally acknowledged that a specific choice of $`S`$ above is immaterial in the thermodynamic limit Huang (1987); Pearson et al. (1985), for therein these expressions differ by terms increasingly small, but in the opposite limit of only a few degrees of freedom they lead to appreciably different observables, most notably the equilibrium temperature (given by the thermodynamic relation $`1/T=S/E`$) Berdichevsky (1997). The temperatures that follow from Eqs. (2) and (3) via this thermodynamic relation will be similarly referred to as the volume ($`T_\mathrm{\Phi }`$) and area ($`T_\mathrm{\Omega }`$) temperatures.

Almost a decade ago, Bianucci et al. (BMWG henceforth) set out to provide a “derivation,” and generalization, of the Boltzmann principle222It is important to stress that what these authors refer to as the Boltzmann principle is the particular case of Eq. (3). \[Eq. (3)\] with particular emphasis on systems far away from the thermodynamic limit Bianucci et al. (1995). Such generalization, which will be concisely reviewed in the next section, is based upon the derivation of an approximate equation of motion for the reduced probability distribution of the degrees of freedom of interest, which is then compared with a postulated form of the Fokker-Planck equation \[Eq. (12) below\]. This comparison allows BMWG to put forward a microscopic expression for the temperature \[Eq. (14) in the present work\], and additional assumptions allow them to relate their temperature expression $`T_B`$ with the phase space area $`\mathrm{\Omega }`$ \[Eq. (18)\]. Since the resulting relationship between $`T_B`$ and $`\mathrm{\Omega }`$ differs from the one that follows from Eq. (3) via the thermodynamic relation $`1/T=S/E`$ by additional terms that only vanish in the thermodynamic limit, Bianucci et al. are led to conclude that they have found a generalization of the Boltzmann principle for finite systems. This fact is confirmed by BMWG through numerical simulations of a particle coupled to a chain of oscillators, in which the authors measure the average of twice the kinetic energy of the particle and compare it with the predicted value from their proposed temperature expression.

It is clear from the above discussion that an implicit assumption in the work of BMWG is the coincidence of their analytic temperature definition $`T_B`$ \[Eq. (18)\] with $`v^2`$ (twice the average kinetic energy, where we temporarily set the mass and the Boltzmann constant to unity), i.e. it is assumed that their temperature definition coincides with, or at least approximates the familiar equipartition one. Since it was actually the latter quantity that was used to assess the quality of the former and not the other way around, one can say that the endeavor of BMWG was to find a dynamical method for evaluating $`v^2`$ in terms of the phase space structure of the system, and that the failure of Eq. (3) in predicting the equilibrium temperature of their numerical experiments is related to the breakdown of the equation $`v^2=(S_\mathrm{\Omega }/E)^1`$ for finite systems. As shown in the present investigation, the BMWG temperature $`T_B`$ is in fact an approximation to the equipartition one, but if one assumes mixing as BMWG did in the final stages of their investigation, the asymptotic value of $`v^2`$ becomes exactly related to $`S_\mathrm{\Phi }`$ \[Eq. (24)\] for any system size, establishing incidentally the equivalence between the equipartition and the volume temperatures. This shows not only that the Boltzmann principle in Eq. (3) is indeed inappropriate to predict the equipartition temperature of small systems Berdichevsky (1997), but also that the dynamical corrections to $`S_\mathrm{\Omega }`$ suggested in the work of Ref. Bianucci et al. (1995) are unnecessary, since an alternate and equally legitimate form of the Boltzmann principle \[Eq. (2)\] is sufficient for that purpose. This observation translates immediately into a negative answer to the question posed in the title of the present work, so long as one answers it in the context of the work of BMWG.

The remaining of this paper is organized as follows. After reviewing the formalism of BMWG (Sec. II), it will be shown that the temperature expression found in Ref. Bianucci et al. (1995) is an approximation to the usual one defined through the equipartition theorem, and hence that the expression derived therein as a generalization of the Boltzmann principle is but a particular approximation to the volume entropy (Sec. III). These findings will be illustrated by means of a Hamiltonian model similar to that of BMWG (Sec. IV), where in addition the distribution of momenta is computed and seen to depart from the canonical one. This provides, in turn, an alternate means of explaining the inability of the area temperature to predict the average of twice the kinetic energy of a particle, as seen in the numerical experiments of Bianucci et al.

## II The dynamical approach of BMWG

In order to understand the origin of the dynamical generalization of the Boltzmann principle mentioned above, it is worth reviewing the scope and results of BMWG. The authors of Ref. Bianucci et al. (1995) were concerned with problems whose degrees of freedom can be separated into “system of interest” and “irrelevant,” the latter playing the role of what is commonly known as heat bath or thermostat in the literature. However, in order to avoid confusion with the usual Fokker-Planck equation approach to their problem, in which the irrelevant degrees of freedom are assumed to be in a specific macroscopic state (usually a canonical one, see e.g. Ford et al. (1965); Zwanzig (1973)), the word booster was adopted instead. Therefore, the general goal of BMWG was to derive the thermostatistical properties of a system coupled to a (finite or infinite) bath whose thermal properties arise naturally from dynamics, bridging the gap between mechanics and thermodynamics much like the ideas of Boltzmann.

Notwithstanding the avoidance of macroscopic assumptions, two fundamental properties concerning the dynamics of the booster had to be assumed in order to make the problem amenable to their analytical approach, namely (i) it is chaotic enough so that its correlation functions decay in a finite time (finite correlation time) and (ii) the average response of the system to an external perturbation is linear (linear response). In addition, the system of interest was taken to be a single particle with coordinates $`x,v`$ linearly coupled to the booster via a “doorway” variable $`\xi `$, so that $`\mathrm{\Delta }\xi `$ can be interpreted as an external force driving the particle (or, vice versa, $`\mathrm{\Delta }x`$ for the booster) whose coupling strength $`\mathrm{\Delta }`$ has to be sufficiently small to justify their perturbative approach (weak coupling condition). Condition (i) implies in particular that the unperturbed autocorrelation function of the doorway variable (with $`_0`$ indicating an equilibrium average in the absence of coupling), viz.

$$\phi (t)\frac{\xi (t)\xi (0)_0}{\xi ^2_0},$$
(6)
vanishes for finite $`t`$ (with $`\xi _0=0`$), which in turn implies that the correlation time

$$\tau _0^{\mathrm{}}𝑑t\phi (t)$$
(7)
is finite. Condition (ii) states that the average behavior of the doorway variable after the introduction of an external time-dependent perturbation $`K(t)`$ to its equation of motion can be written as the following convolution (recall the assumption $`\xi _0=0`$ above)

$$\xi (t)_K=_0^t𝑑uS(u)K(tu)+𝒪(K^2),$$
(8)
where $`S(u)`$ is a function that dictates how the average of $`\xi `$ responds to the small perturbation $`K(t)`$ (response function). Together with (i), condition (ii) furnishes two other finite quantities, the asymptotic susceptibility

$$\chi \underset{t\mathrm{}}{lim}\chi (t)\underset{t\mathrm{}}{lim}_0^t𝑑uS(u),$$
(9)
and the response time

$$\vartheta _0^{\mathrm{}}𝑑tc(t),$$
(10)
where the convenient variable $`c(t)1\frac{\chi (t)}{\chi }`$ was introduced so that $`c(0)=1`$ and $`c(\mathrm{})=0`$.

Relying on the above assumptions (and also on a few others concerning the separation of time-scales) and equipped with a specific projection technique, Bianucci and his collaborators embarked on the rather remarkable endeavor of “integrating out” the irrelevant degrees of freedom of the booster without the precise knowledge of its Hamiltonian structure, ultimately in the search of a differential equation governing the time evolution of the reduced distribution function of the particle, which can be formally written as Bianucci et al. (1995)

$$\sigma (x,v;t)𝑑\xi 𝑑𝝅\rho (x,v,𝝅,\xi ;t),$$
(11)
where $`\rho (x,v,𝝅,\xi ;t)`$ is the distribution function of all the dynamical variables and $`𝝅`$ denotes the remaining degrees of freedom of the booster. The resulting differential equation, Eq. (47) in Ref. Bianucci et al. (1995), is then identified with a general 2D Fokker-Planck equation of the type

$`{\displaystyle \frac{}{t}}\sigma (x,v;t)`$ $`=`$ $`\{_a^{\text{eff}}+{\displaystyle \frac{}{v}}A(x,v)({\displaystyle \frac{}{v}}+{\displaystyle \frac{mv}{kT_B}})`$ (12)
$`+{\displaystyle \frac{}{v}}B(x,v)({\displaystyle \frac{}{x}}+{\displaystyle \frac{U^{}(x)}{kT_B}})\}\sigma (x,v;t),`$
where

$$_a^{\text{eff}}=\frac{U^{}(x)}{m}\frac{}{v}v\frac{}{x}$$
(13)
is the effective time evolution operator (Liouvillian) of the unperturbed particle. From a comparison between the coefficients of their equation and the structure of the above equation, they suggest (despite the implicit dependence on the particle coordinates $`x,v`$) that the following expression should be identified with the temperature of the booster (Eq. (52) in Ref. Bianucci et al. (1995)):

$$kT_B(x,v)\frac{\xi ^2_0}{\chi }\frac{_0^{\mathrm{}}𝑑u\phi (u)\left[\frac{}{x}x_a(tu)\right]}{_0^{\mathrm{}}𝑑uc(u)\left[\frac{}{x}x_a(tu)\right]},$$
(14)
where $`x_a(tu)`$ is the unperturbed evolution of the position of the particle “backwards in time” (being essentially the responsible for the implicit dependence of $`T_B`$ on $`x,v`$). The subscript in $`T_B`$ stands for BMWG. To render the above equation independent of $`x,v`$, BMWG consider the following three cases: (a) the particle is a harmonic oscillator with $`V(x)=m\omega ^2x^2/2`$, (b) $`\phi (u)=c(u)`$, and (c) a “natural time scale” (NTS) condition exists separating the typical timescale of the system of interest from the relaxation TS of the booster (the former being much greater than the latter). Under condition (a), the temperature reduces to

$$kT_B=\frac{\xi ^2_0}{\chi }\frac{\text{Re}\left[\widehat{\phi }(\omega )\right]}{\text{Re}\left[\widehat{c}(\omega )\right]},$$
(15)
with hats indicating a Fourier transform, whereas in the case of NTS (c) one obtains

$$kT_B=\frac{\xi ^2_0\tau }{\chi \vartheta },$$
(16)
which coincides with the one obtained within a Langevin approach (cf. Eq. (22) in Ref. Bianucci et al. (1995)), as expected from the dichotomy of the time scales. Condition (b), however, is peculiar in the sense that it recovers the usual canonical results, which (with the exception of making the connection with Kubo’s linear response theory clear) was not particularly relevant for the remaining part of their work.

To draw a parallel with Boltzmann’s principle, BMWG invoked mixing for the booster, in which case the time-dependent susceptibility can be found via “geometrical” arguments to assume the form Bianucci et al. (1994)

$$\chi (t)=\frac{1}{\mathrm{\Omega }_b(E)}\frac{}{E}\left\{\mathrm{\Omega }_b(E)\xi ^2_0[1\phi (t)]\right\},$$
(17)
where $`\mathrm{\Omega }_b(E)`$ is the phase space area of the unperturbed booster with energy $`E`$. With this result in hand, one can write, for the case (c) above for example,

$$\frac{1}{kT_B}=\frac{}{E}\mathrm{ln}\mathrm{\Omega }_b(E)+\frac{}{E}\mathrm{ln}(\xi ^2_0\tau ),$$
(18)
whereas for case (a) one simply replaces $`\tau `$ with $`\text{Re}\left[\widehat{\phi }(\omega )\right]`$ in this expression. This result should be contrasted with the one obtained from what BMWG calls the Boltzmann principle \[i.e. Eq. (3)\], namely

$$\frac{1}{kT_\mathrm{\Omega }}=\frac{}{E}\mathrm{ln}\mathrm{\Omega }(E).$$
(19)
The additional term in Eq. (18), called by BMWG a “dynamical correction” to Eq. (19), is then argued to become negligible for sufficiently large systems, thus recovering the result above and causing no conflict with the standard macroscopic thermostatistics of Boltzmann, i.e. the Boltzmann principle has been “generalized” to cover systems with only a few degrees of freedom. In the following section it is shown that this generalization is actually an approximate form of the volume temperature $`T_\mathrm{\Phi }`$.

## III Connection with the volume temperature and entropy

Besides heuristic arguments of units, it is not clear how one can identify the “constant” $`T_B`$ \[Eq. (14)\] obtained by BMWG through a comparison with Eq. (12) as an actual temperature. In the usual derivation of the Fokker-Planck equation (or, equivalently, of the Langevin equation) for a Brownian particle Reif (1965); Reichl (1998), the temperature comes into play by invoking the equipartition theorem, viz. $`kTmv^2`$ (this is indeed the case for the Langevin treatment of BMWG, cf. Eq. (22) in Ref. Bianucci et al. (1995)), but in the case of the Fokker-Planck equation derived by BMWG or the generic one in Eq. (12), both with $`x,v`$-dependent coefficients, it is difficult to see whether their temperature definition is compatible with equipartition. Nevertheless, after the simplifications introduced by any of the cases (a), (b), or (c) of the previous section (the only ones considered by BMWG in their generalization of the Boltzmann principle), one can show that the definition of $`kT_B`$ coincides with an approximate form of $`mv^2`$. Indeed, focusing on condition (c) for notational simplicity (the remaining cases (a) and (b) follow by a trivial modification of the procedure below), the time-independent differential equation satisfied by $`\stackrel{~}{\sigma }(x,v)`$ obtained via the projection method reduces to (cf. Eq. (47) of Ref. Bianucci et al. (1995))

$`0=`$ $`v{\displaystyle \frac{\stackrel{~}{\sigma }}{x}}+{\displaystyle \frac{V^{}(x)}{m}}{\displaystyle \frac{\stackrel{~}{\sigma }}{v}}{\displaystyle \frac{\mathrm{\Delta }^2x\chi }{m}}{\displaystyle \frac{\stackrel{~}{\sigma }}{v}}+{\displaystyle \frac{\mathrm{\Delta }^2\xi ^2_0\tau }{m^2}}{\displaystyle \frac{^2\stackrel{~}{\sigma }}{v^2}}`$
$`+{\displaystyle \frac{\mathrm{\Delta }^2\xi ^2_0\eta ^2}{m^2}}{\displaystyle \frac{^2\stackrel{~}{\sigma }}{vx}}+{\displaystyle \frac{\mathrm{\Delta }^2\chi \vartheta }{m}}{\displaystyle \frac{}{v}}(v\stackrel{~}{\sigma })`$
$`+{\displaystyle \frac{\mathrm{\Delta }^2\chi \beta ^2}{m}}{\displaystyle \frac{\stackrel{~}{\sigma }}{v}},`$ (20)
where the constants $`\eta ^2`$ and $`\beta ^2`$ are defined in Eqs. (6) and (13) of Ref. Bianucci et al. (1995). The tilde indicates that the distribution here is not identical to Eq. (11) due to the approximations carried out in the actual process of integration, or “projection” (see below). Multiplying the above equation through by $`v^2`$ and integrating over $`x,v`$, most terms disappear under the assumption that $`\stackrel{~}{\sigma }(x,v)0`$ sufficiently fast for $`x,v\mathrm{}`$ and $`x_{\stackrel{~}{\sigma }}=v_{\stackrel{~}{\sigma }}=0`$ (for simplicity I assume the average position of the particle is centered at the origin, though it is not difficult to adapt this derivation for the general case), so one is left with

$$0=\frac{2\mathrm{\Delta }^2\xi ^2_0\tau }{m^2}\frac{2\mathrm{\Delta }^2\chi \vartheta }{m}v^2_{\stackrel{~}{\sigma }},$$
(21)
where $`_{\stackrel{~}{\sigma }}_{\mathrm{}}^{\mathrm{}}𝑑x_{\mathrm{}}^{\mathrm{}}𝑑v\stackrel{~}{\sigma }(x,v)`$. Solving this equation for $`mv^2_{\stackrel{~}{\sigma }}`$ we obtain the right hand side of Eq. (16), i.e. $`mv^2_{\stackrel{~}{\sigma }}=kT_B`$. This is the desired connection between $`T_B`$ and the equipartition temperature, q.e.d.

It is important to emphasize, however, that the reduced density $`\stackrel{~}{\sigma }`$ (and hence the aforementioned connection) is only approximate. This can be seen in the discussion of paragraph 1, page 3005 of Ref. Bianucci et al. (1995) – apart from conditions (i) and (ii) of the previous section, it is always assumed that there is a large separation of relaxation and intrinsic time scales of the booster. This condition, which amounts to an assumption of weak coupling, was used explicitly between Eqs. (39) and (40) of Ref. Bianucci et al. (1995) (note also that Eq. (39) of Ref. Bianucci et al. (1995) is itself a perturbative result). We now note that the exact average $`mv^2_\sigma `$, obtained with $`\sigma `$ \[Eq. (11)\] before the approximations carried out in the projection method of BMWG which lead to $`\stackrel{~}{\sigma }`$, furnishes the following expression

$$mv^2_\sigma =m𝑑x𝑑v𝑑\xi 𝑑𝝅v^2\frac{\delta [EH(x,v,\xi ,𝝅)]}{\mathrm{\Omega }(E)},$$
(22)
where $`H(x,v,\xi ,𝝅)`$ is the Hamiltonian of the system + booster, and

$$\rho (x,v,\xi ,𝝅)=\frac{\delta [EH(x,v,\xi ,𝝅)]}{\mathrm{\Omega }(E)}$$
(23)
is the asymptotic solution of the Liouville equation for mixing systems (microcanonical distribution, cf. Bianucci et al. (1994)), in accordance with the assumption of mixing in Sec. IV-D of Ref. Bianucci et al. (1995). By means of simple integral manipulations used in the generalized equipartition theorem (see e.g. Huang (1987); Berdichevsky (1997)), we can rewrite Eq. (22) in terms of the volume entropy $`S_\mathrm{\Phi }=k\mathrm{ln}\mathrm{\Phi }(E)`$ as

$$mv^2_\sigma =\left[\frac{}{E}\mathrm{ln}\mathrm{\Phi }(E)\right]^1=k\left(\frac{S_\mathrm{\Phi }}{E}\right)^1kT_\mathrm{\Phi }.$$
(24)
Observe that, contrary to the approximate conclusion made after Eq. (21), the above equation establishes an exact equivalence between the equipartition and the volume temperatures, as already anticipated in Sec. I, showing that the knowledge of the phase space volume alone (or $`S_\mathrm{\Phi }`$) is sufficient to predict the equilibrium value of $`mv^2_\sigma `$. Though not new Huang (1987); Berdichevsky (1997), this result renders the dynamical corrections introduced in Ref. Bianucci et al. (1995) unnecessary. This fact is illustrated in the next section by means of explicit analytical and numerical calculations.

## IV Quartic oscillators: An analytically tractable model

In this last section a specific Hamiltonian model will be analyzed both analytically and numerically. Since BMWG relied on a similar model to illustrate their findings, the present analysis serves the similar purpose of illustrating on an equal footing the results established in the last section. It also shows that, for such small systems where Eqs. (2) and (3) yield appreciably different results, a noncanonical distribution of momenta is a necessary and fundamental property, in contrast to what is argued in Ref. Bianucci et al. (1995).

The model adopted herein is described by the Hamiltonian structure

$$H(𝐳)=H_1(x,p)+H_2(\mathit{\varphi },𝝅)+gH_{1,2}(x,\varphi _0),$$
(25)
where $`x,p`$ and $`\mathit{\varphi }=(\varphi _0,\mathrm{},\varphi _{N1})`$, $`𝝅=(\pi _0,\mathrm{},\pi _{N1})`$ are the particle and booster canonically conjugate coordinates, respectively, $`g`$ is a controllable coupling constant and

$`H_1(x,p)`$ $`={\displaystyle \frac{p^2}{2}}+{\displaystyle \frac{1}{2}}\omega ^2x^2,`$ (26)
$`H_2(\mathit{\varphi },𝝅)`$ $`={\displaystyle \underset{i=0}{\overset{N1}{}}}{\displaystyle \frac{\pi _i^2}{2}}+{\displaystyle \underset{i=1}{\overset{N1}{}}}\left[{\displaystyle \frac{\varphi _i^4}{4}}+{\displaystyle \frac{(\varphi _i\varphi _{i1})^4}{4}}\right],`$ (27)
$`H_{1,2}(x,\varphi _0)`$ $`=x\varphi _0,`$ (28)
are the particle, booster and coupling Hamiltonians, respectively. The intrinsic time scale of the particle is set by $`\omega ^1`$, whereas its coupling with the booster is set linearly through the “doorway” variable $`\varphi _0`$. It should be noted that this model can be rendered unphysical when $`\omega 0`$ (i.e. in the free particle limit) since the energy would no longer be bounded from below (due to the remaining linear coupling term $`H_{1,2}`$) and thence no stable thermodynamical state would exist Ford et al. (1988) (that shall not concern us here since $`\omega `$ will always be non-zero).

The essential difference between the above model and the one adopted by BMWG is the absence of harmonic terms in the booster Hamiltonian, which permits a simple analytic treatment based on scaling arguments, as we shall see in the next sections. While not particularly relevant for the present analysis, it will also be shown that this model retains most of the quantitative features of the model of Ref. Bianucci et al. (1995) (see, in particular, $`q`$ and $``$ below). More important than these quantitative similarities is its mixing behavior which, just as in the model of BMWG, is suggested by the numerical experiments through the observed finite relaxation time towards a unique equilibrium state.

For the numerical simulations to be reported soon, the equations of motion were evolved with an optimized fourth-order symplectic integrator McLachlan and Atela (1992), resulting in numerical energy fluctuations $`\mathrm{\Delta }E/E`$ smaller than one part in $`10^6`$. Initial conditions were picked randomly from the density $`\rho (x,p,\mathit{\varphi },𝝅;t=0)=\delta (x)\delta (p)\delta (H_2E)`$ Bianucci et al. (1995), and no dependence of the observables in equilibrium on these configurations was observed (being already a good indication of ergodicity by Birkhoff’s theorem Berdichevsky (1997)). The typical integration time used was $`t=10^8`$, after allowing a suitable relaxation time of order $`t_{rel}=10^7`$. Since we are not interested in nonequilibrium features, no ensemble average was done (the coincidence of time and phase space averages relying on the Birkhoff-Khinchin theorem for ergodic systems Berdichevsky (1997)). For simplicity, the harmonic frequency of the particle, $`\omega `$, was set to unity (in contrast to the approach of BMWG, the specific value of $`\omega `$ is immaterial for the validity of the method adopted here).

### IV.1 Momentum distribution

As a first probe of ergodicity, one can test whether the observed histogram of the momentum of the particle $`f(p)`$ approximates the distribution law obtained from the structure function $`\mathrm{\Omega }_p(E)`$ of the remaining degrees of freedom $`x`$, $`\mathit{\varphi }`$ and $`𝝅`$ (see below), i.e. for large enough integration times one expects

$$f(p)\frac{\mathrm{\Omega }_p(Ep^2/2)}{\mathrm{\Omega }(E)},$$
(29)
where $`\mathrm{\Omega }(E)`$ is the structure function of the whole system. Most quantities of interest here, including $`\mathrm{\Omega }_p`$, can be obtained from the ($`2N+1`$)-dimensional phase space volume

$$\mathrm{\Phi }_p(E)=𝑑x𝑑\mathit{\varphi }𝑑𝝅\theta \left(E[x^2/2+H_2(\mathit{\varphi },𝝅)]\right),$$
(30)
where $`\theta `$ is the unit step function and the coupling Hamiltonian $`H_{1,2}`$ was neglected (weak coupling approximation). In fact, assuming the above volume to be bounded and sufficiently smooth for all $`E>0`$ (such assumption will be adopted throughout this section), the phase space area $`\mathrm{\Omega }_p`$ can be obtained from $`\mathrm{\Phi }_p`$ via the general relation Khinchin (1949)

$$\mathrm{\Omega }_i(E)=\frac{d\mathrm{\Phi }_i(E)}{dE},$$
(31)
where $`i`$ labels the desired (sub)volume. Though evaluating Eq. (30) can be a hard task for general Hamiltonians, for the choice presented in Eqs. (26)-(28) a simple scaling argument can be used to show that (cf. Ref. Adib et al. (2003)):

$$\mathrm{\Phi }_p(E)=\stackrel{~}{a}E^{(3N+2)/4},$$
(32)
where henceforth the symbols $`\stackrel{~}{a},\stackrel{~}{b},`$ etc. will refer to any function that does not depend on $`E`$. Equation (31) then gives

$$\mathrm{\Omega }_p(E)=\stackrel{~}{b}E^{(3N2)/4},$$
(33)
and therefore, by inserting this relation in Eq. (29), one gets the desired distribution law:

$$f(p)=\frac{1}{Z}\left(1\frac{1}{E}\frac{p^2}{2}\right)^{(3N2)/4},$$
(34)
where $`Z`$ is a normalization constant. This expression is valid for any $`N>0`$, provided the dynamics of the system is sufficiently ergodic (recall that $`f(p)`$ was defined as a histogram). Note also that in the limit $`N\mathrm{}`$ with $`E/N`$ constant one recovers the usual Boltzmann factor $`\mathrm{exp}(\beta p^2/2)`$, with $`\beta =3N/4E`$.

Figure 1 shows the momentum distribution function of the harmonic oscillator obtained from the numerical simulations for the booster model in Eq. (27), from which one obtains good evidence for ergodicity given its clear agreement with Eq. (34). The “departure from Gaussianity” introduced by BMWG was also computed, namely $`q=1p^4/[3p^2^2]`$, to show not only that the distribution is close to canonical (see below however), but also that its numerical value is very close to the one obtained from the model of Ref. Bianucci et al. (1995). Indeed, the observed value of $`q`$ in the present model was $`q=0.1255`$ for $`N=8`$, which should be compared to $`q0.1`$ obtained by BMWG for the same number of particles.

In the following section we will be observing the temperature of the system via the average momentum square of the particle. Here I show that the mismatch between $`T_\mathrm{\Phi }`$ and $`T_\mathrm{\Omega }`$ to be found in that section can alternatively be attributed to the noncanonicity of $`f(p)`$. Indeed, using the integral

$$_{1/\sqrt{a}}^{1/\sqrt{a}}p^n\left(1ap^2\right)^b𝑑p=\frac{\left(1+(1)^n\right)}{2a^{\frac{1}{2}(n+1)}}\frac{\mathrm{\Gamma }(b+1)\mathrm{\Gamma }\left(\frac{1+n}{2}\right)}{\mathrm{\Gamma }\left(\frac{3}{2}+b+\frac{n}{2}\right)},(a,b>0)$$
(35)
where $`\mathrm{\Gamma }(u)`$ is the gamma function and $`n=0,1,2,\mathrm{}`$, we obtain, after using the property $`\mathrm{\Gamma }(u+1)=u\mathrm{\Gamma }(u)`$ and reading the values of $`a,b`$ from Eq. (34),

$$p^2=_{\sqrt{2E}}^{\sqrt{2E}}p^2f(p)𝑑p=\frac{4E}{3N+4}$$
(36)
(compare with Eq. (39) below). Note that this result differs from the one obtained using the Boltzmann factor, namely $`p^2=4E/3N`$ (compare Eq. (40)). This difference, though small for large $`N`$, is ultimately the motivation for BMWG to conclude that their corrected temperature is indeed an improvement over the “standard” one ($`T_\mathrm{\Omega }`$, obtained from Boltzmann’s principle), which shows additionally that, contrary to the discussion that accompanies Fig. 6 of Ref. Bianucci et al. (1995), the small departure from canonicity of the momentum distribution is a pertinent property of these small systems.

### IV.2 Temperature

Using the same weak coupling condition adopted in the last section, one can neglect the interaction Hamiltonian and compute both the phase space volume and area (this time for the whole system) through a procedure similar to the one leading to Eqs. (32) and (33). One then finds, respectively,

$`\mathrm{\Phi }(E)`$ $`=\stackrel{~}{c}E^{(3N+4)/4},`$ (37)
$`\mathrm{\Omega }(E)`$ $`=\stackrel{~}{d}E^{3N/4},`$ (38)
from which the temperatures $`T_\mathrm{\Phi }`$ and $`T_\mathrm{\Omega }`$ follow trivially, viz.

$`T_\mathrm{\Phi }`$ $`=\left({\displaystyle \frac{\mathrm{ln}\mathrm{\Phi }(E)}{E}}\right)^1={\displaystyle \frac{4E}{3N+4}}={\displaystyle \frac{4ϵ}{3}}{\displaystyle \frac{1}{1+4/3N}},`$ (39)
$`T_\mathrm{\Omega }`$ $`=\left({\displaystyle \frac{\mathrm{ln}\mathrm{\Omega }(E)}{E}}\right)^1={\displaystyle \frac{4E}{3N}}={\displaystyle \frac{4ϵ}{3}},`$ (40)
where $`ϵE/N`$ is the total energy per degree of freedom and $`k`$ was set to unity (this harmless simplification will be adopted throughout the remaining part of the paper). In what follows, the numerically observed temperature $`T_{obs}`$ is defined as the time-average $`\overline{v^2}`$, which should coincide with the ensemble average $`v^2`$ \[and hence with $`T_\mathrm{\Phi }`$, cf. Eq. (24)\] for ergodic systems.

In Fig. 2 the observed temperatures from the numerical simulations with $`ϵ=1.0`$, $`g=0.1`$ and various values of $`N`$ are presented, as well as the analytic predictions Eqs. (39) and (40). The agreement with $`T_\mathrm{\Phi }`$ is evident, and we can see clearly from this figure that both temperatures should coincide in the thermodynamical limit Berdichevsky (1997). In order to verify that these results can in fact explain the discrepancy found by BMWG between the “standard” form of Boltzmann’s principle \[which gives rise to Eq. (19)\] and their numerical simulations, the relative error $`=(T_{obs}T_\mathrm{\Omega })/T_\mathrm{\Omega }`$ was computed both for the data available here and in that reference (the temperatures from the latter are $`T_\mathrm{\Omega }14.95`$ and $`T_{obs}11.78`$), giving $`_A15.2\%`$ and $`_B21.2\%`$, respectively. Note that not only the magnitude of these relative deviations are quite close, but also their sign is the same, i.e. the general trend of $`T_\mathrm{\Omega }`$ is to overestimate the final equilibrium temperature, just as BMWG found (see also Bannur et al. (1997)). The additional error $`_A_B6\%`$ can be further explained from the fact that the phase space area in the approximate temperature expression of BMWG involves the degrees of freedom of the booster only, neglecting the contribution of $`x`$ and $`p`$ which tends to increase $`T_\mathrm{\Omega }`$ and thus to decrease $``$. This is the source of the “$`10\%`$” deviation from the observed temperature mentioned before the conclusions of Ref. Bianucci et al. (1995).

## V Conclusions

It is important to emphasize that the purpose of the present contribution was not to determine whether the volume or the area entropy is the correct starting point of thermostatistics for small systems. Although the former bears an important thermodynamic property that is independent of the system size (known as “adiabatic invariance,” cf. Adib (2002) and references therein), these very small and isolated systems such as the ones studied by BMWG are usually far removed from experience, and we do not know whether our empirical and macroscopic thermodynamics acquires a different form in this case. Experiments with small clusters of atoms is a field that can potentially settle this question Freeman and Doll (1996), and a theoretical comparison between the two entropies above has already been carried out in this context Jellinek and Goldberg (2000).

Quite apart from elucidating the ad hoc definition of the BMWG temperature \[cf. discussion following Eq. (21)\], the present contribution has established that one does not have to resort to dynamical quantities such as time correlation functions in order to predict the equilibrium temperature of a small system. Such idea was put forward in the work of Ref. Bianucci et al. (1995) in the form of a correction to the temperature that follows from Eq. (3) for the particular case of a system composed of a particle coupled to a set of “irrelevant” degrees of freedom, leading its authors to conclude that a correction for the Boltzmann principle is necessary when one is dealing with small systems. The key identity that escaped the attention of BMWG is Eq. (24), which establishes an exact relation between the equipartition temperature and the phase space volume, a quantity that can be computed much in the same way as the usual phase space area. This shows that, as long as we are working with the assumptions of Ref. Bianucci et al. (1995), the particular form of the Boltzmann principle in Eq. (2) is sufficient for predicting the equipartition temperature, without any dynamical correction.

###### Acknowledgements.

The author is indebted to Prof. Oliver Penrose and the anonymous referee for their constructive and detailed comments, which largely reshaped the original manuscript. Correspondences or private discussions with Professors Victor Berdichevsky, David Montgomery, Murilo Almeida and Paolo Grigolini are also acknowledged. Financial support and computational resources (the latter through the Field Theory/Cosmology group) were partially provided by the Dartmouth College.


format: text

text:
# Cataclysmic variable evolution: AM Her binaries and the period gap

## 1 Introduction

The disrupted magnetic braking model has enjoyed considerable success in explaining the $`23`$ hour period gap in the period distribution of the cataclysmic variables (CVs). The essential idea is that the CVs with orbital periods $`P_{\mathrm{orb}}10^\mathrm{h}`$ evolve towards shorter periods due to a combination of angular momentum losses from magnetic braking due to a magnetised wind originating from the secondary star, and from gravitational radiation (Verbunt & Zwaan ; Rappaport, Verbunt & Joss ; see Howell, Nelson & Rappaport  for a recent review). Magnetic braking dominates at longer periods. As the donor secondary decreases in mass, its thermal time-scale increases, until thermal relaxation can no longer keep pace with the mass transfer rate being driven by magnetic braking, and it becomes increasingly bloated with respect to thermal equilibrium. When, as the orbital period reaches ($`3^\mathrm{h}`$), the secondary’s convective envelope penetrates to the stellar center, its large-scale magnetic field loses its radiative anchor and, being buoyant, is expelled from the star. This sudden cessation of magnetic braking allows the distended secondary to contract more rapidly than its Roche lobe, and mass transfer ceases until angular momentum losses to gravitational radiation have reduced the binary separation far enough for the secondary to again fill its Roche lobe, at an orbital period of $`2^\mathrm{h}`$. It is important to recognize that, were the secondary to remain in thermal equilibrium throughout its evolution, it would evolve continuously through the period gap with no abatement of mass transfer.

A subset of the CVs are known to be magnetic and classified as either AM Hers or intermediate polars (IPs) . In the AM Hers the magnetic white dwarf is in synchronous rotation with the orbit so that its spin period $`P_{\mathrm{spin}}`$ is equal to the orbital period $`P_{\mathrm{orb}}`$. In the IPs, the magnetic white dwarf is asynchronous (typically $`P_{\mathrm{spin}}0.1P_{\mathrm{orb}}`$). Some 65 confirmed AM Hers are known, of which 37 have direct magnetic field measurements (see Wickramasinghe & Ferrario  for a recent review). The estimated effective dipolar fields range from $`7200`$ MG. The IP class contains fewer confirmed members (26). Indirect estimates of their field strengths suggest dipolar fields in the range $`530`$ MG. Any theory which purports to explain the period distribution of CVs must also account for the properties of the magnetic CVs which, as a group, comprise some $`20\%`$ of all CVs of known orbital period .

The possibility that the magnetic field of the white dwarf may affect the orbital evolution of the binary was first discussed by Wickramasinghe & Wu  in relation to the AM Her type systems. This research was prompted by the early discovery of several AM Herculis systems in the period gap before the results of the enlarged sample of AM Hers from the ROSAT survey became available. Here it was argued that at strong fields, the magnetised wind from the secondary may become trapped by the combined magnetosphere of the primary and the secondary and result in a reduction in the magnetic braking rate. If the braking rate was effectively curtailed, systems will evolve towards shorter periods almost entirely by gravitational radiation. The proposal was that there will be no disruption in magnetic braking, and that no period gap will form.

Li, Wu & Wickramasinghe \[1994a, 1994b, 1995\] extended the Mestel & Spruit  wind theory for single stars and developed a theory for evaluating the braking rate in synchronised AM Herculis type binary systems. These models assume that the secondary star has a dynamo generated field with a dipole moment oriented perpendicular to the orbital plane, and that the magnetic white dwarf is locked in synchronous rotation by magnetostatic torques. The dipole moment of the white dwarf is then oriented anti-parallel to that of the secondary star. The calculations showed that for a fixed orbital period, the braking rates dropped by several orders of magnitude as the polar field strength increased from $`10`$ to $`100`$ MG, and the magnetosphere of the white dwarf engulfed the companion star and trapped the wind. Subsequently, Li & Wickramasinghe  allowed for more general orientations of the dipole moment of the secondary star and showed that the largest reductions in the braking rate occurred for the case where the dipole moment was nearly (within $`20\mathrm{°}`$) aligned with the spin axis.

The ROSAT survey  resulted in a large increase in the number of known AM Hers , and allowed a better statistical analysis of the period distribution of the magnetic and non-magnetic CVs. At long orbital periods, the numbers of AM Hers falls much more rapidly than those of non-AM Hers, a difference attributable to a competition between accretion torques and magnetic torques which favors magnetic torques at low donor masses, short orbital periods, and small orbital separations. As CVs evolve to shorter orbital period, then, a growing fraction achieve synchronism and evolve from non-AM Her to AM Her systems. Meaningful comparisons of AM Her to non-AM Her period distributions are therefore limited to short orbital periods, where statistics are presumably dominated by the influx of systems from longer orbital periods, as opposed to birth within the period interval of interest. Wheatley’s  application of the two-sided Kolmogorov–Smirnov test to systems with orbital periods less than 3$`\stackrel{h}{.}`$2 yielded a $`15\%`$ probability that AM Her and non-magnetic distributions were drawn from the same parent population, a difference he regarded as not significant. (Hellier & Naylor , on the other hand, found a similar figure for systems with orbital period less than $`5^\mathrm{h}`$, and regarded it as significant.) We shall see below that the differences between AM Her and non-AM Her distributions should be most evident within the period gap itself, and that a two-sided Kolmogorov–Smirnov test applied only to systems *within the gap* yields a $`8.5\%`$ probability that they belong to the same parent distribution. This comparison is still hampered by small numbers, but we regard the difference as compelling, if not yet conclusive.

In this paper we reports results of calculations of the orbital evolution of magnetic CVs which allow for a reduction in magnetic braking following the theory of Li et al. \[1994a\]. We use our results to analyse the data on the period distribution of the magnetic CVs, and more generally to re-assess the role played by disrupted magnetic braking in the formation of the period gap.

## 2 Observational constraints

Magnetic field measurements from Zeeman and cyclotron spectroscopy provide estimates either of the effective dipolar photospheric field, or of magnetic fields at accretion spots on the white dwarf surface. We use as our data base for AM Hers the list of field measurements given in  based on these techniques updated to include the recent field measurement of V884 Her . The braking calculations, however, require the magnetic moment which is difficult to estimate from observations. The magnetic moment while being proportional to the magnetic field is also proportional to the cube of the stellar radius. Unfortunately, the latter, which is derived from the mass, is poorly determined. Indeed, even in the few eclipsing systems where the orbital inclination is well constrained, there is no consensus on the component masses, and these are even less well determined in other systems. In our calculations, the magnetic moments have been chosen to be in the range $`10^{32}10^{35}`$ G cm3 which is consistent with the observed effective dipolar field range if allowance is made for a spread in white dwarf mass.

Our list of CVs is drawn from the database of Downes et al.  current on 2002 January 25. We have included the intermediate polars in our ‘non-magnetic’ sample because to a first approximation, we expect them to evolve under the normal braking law. However, it is likely that many IPs become synchronised as they evolve toward shorter periods, becoming AM Her systems. This possibility needs to be considered when comparing observed distributions.

We show in Figure 1 the distributions of AM Hers and other types of CVs as a function of orbital frequency (hr-1), binned into frequency intervals 0.025 hr-1. In these units the canonical ‘period gap’ corresponds to the frequency range $`0.3250.475`$ hr-1. We choose to display the distribution in terms of frequency, rather than more conventionally in terms of period, in order to spread the short-period peak and to make the period gap more apparent. In this figure, and throughout the rest of this paper, we will continue to display results in terms of orbital frequency for ease of comparison; however, we will cast our discussion in the more familiar terms of orbital period, and show at the top of each figure the corresponding period scale.

A glance at Figure 1 shows the following characteristics:

There is an abrupt drop in the number of non-AM Her systems at periods above 2.1 hours, which can be considered the lower end of the gap. At longer periods, there is a gradual rise in their number distribution, with the distribution rising rapidly at $`3`$ hours, which can be considered the upper end of the period gap. This plot clearly shows that the gap is itself not totally devoid of stars. Rather there is a tendency for this canonical period gap to fill in from the higher period end.

The period gap is decidedly less well-defined among the AM Her stars. The average number of systems per unit orbital frequency within the 2.1–3 hour period range is smaller than among systems at shorter period (higher frequency), but the increase in numbers at longer period (lower frequency) is much less pronounced than among the non-AM Her systems. Moreover, among the magnetic systems, the canonical period gap appears to fill from the lower period (higher frequency) end, opposite to the behaviour seen for the other CVs. The Kolmogorov–Smirnov two-sample statistic, applied only to those systems of confirmed orbital period falling within the gap in Figure 1, yields only a $`8.5\%`$ probability that AM Her systems and non-AM Her systems are drawn from the same population.

We show in Figure 2 the distribution of field strength as a function of orbital frequency for the 37 systems for which field estimates are available. The fields correspond to effective dipolar fields for cases where photospheric Zeeman measurements are available, or where cyclotron fields are available for both poles. For other systems, the fields correspond to the value measured at a single pole. This plot shows that there is no significant difference between the mean field of systems within the canonical period gap and outside. Furthermore, there is no obvious correlation between field strength and orbital period for systems within the $`2^\mathrm{h}3^\mathrm{h}`$ period region.

## 3 The magnetic braking of synchronised binaries

### 3.1 The parametrisation

The braking rate of synchronised magnetic binaries depends on the interaction of the stellar wind with the combined magnetospheres of the components stars. It is therefore a function of a number of parameters which fall into two basic categories: those which specify the magnetic interaction, and those which specify the properties of the wind.

There are four parameters belonging to the first category which enter directly into the calculations of the braking rates. These are the masses $`M_1`$ and $`M_2`$ of the primary and the secondary stars respectively, the radius $`R_2`$ of the secondary, and the magnetic moment $`\mu _1`$ of the primary star. The magnetic moment $`\mu _2`$ of the secondary is not a free parameter but is calculated on the assumption that prior to the star becoming fully convective, the mean surface field evolves with angular velocity according to a power law

$$B_2=\left(\frac{\omega }{\omega _{\mathrm{}}}\right)^pB_{\mathrm{}}$$

where $`\omega `$ is the angular frequency of rotation. Mestel & Spruit  considered the cases $`p=1`$ as being most consistent with the observations of late type stars, and we accordingly adopt $`p=1`$ in our models. However, it should be noted that this law does not allow for complexities such as dynamo saturation which may play a role in binary evolution where rapid rotation rates are encountered as the orbit shrinks.

The magnetic moment of the secondary star

$$\mu _2=0.5B_2R_2^3$$

is a strong function of its radius $`R_2`$. Among the CV donors we wish to model, that radius is itself strongly correlated with the stellar mass. We therefore parametrise $`R_2`$ in terms of a pseudo-radius $`_2`$ which removes the lowest-order power-law relationship between $`R_2`$ and $`M_2`$ . Thus,

$$R_2=f_2,$$

where

$$_2=0.094\mathrm{R}_{\mathrm{}}(M_2/0.065\mathrm{M}_{\mathrm{}})^{13/15}.$$

After the star becomes convective we assume that the intrinsic mean field is effectively zero. The secondary star will have an induced dipole moment of

$$\mu _2\mu _1(R_2/a)^3$$

where $`a`$ is the orbital separation. The induced dipole moment is typically a factor $`10100`$ times lower than the intrinsic dipole moment of the primary star, and there will be no open flux to produce braking even if a wind persists. In this regime we can therefore assume that braking is dominated by angular momentum loss by gravitational radiation.

It can be seen from the above, that apart for the wind parameters to be discussed in the next section, the braking rate can be parametrised terms of $`M_1,f,\mu _1,M_2`$.

### 3.2 The calculations of the magnetic braking rate

In the braking of single stars by a magnetised stellar wind, both thermal and centrifugal effects can play an important role in different mass and rotation regimes. The Mestel & Spruit  theory allows a convenient parametrisation of these effects. In the context of the CVs, however, due to their rapid rotation we can assume that centrifugal effects dominate over the entire observed range \[1994a, 1994b\].

The formulation of the magnetic braking problem for centrifugally driven winds has been described in detail by Li et al. \[1994a\]. We therefore present only the basic equations and our method of solution, and refer the reader to this paper for further details.

The basic assumption underlying our calculations is that the dipole moment of the secondary star is oriented perpendicular to the orbital plane – the likely orientation for a dynamo generated field. The magnetic moment of the primary star must then be anti-aligned with that of the secondary in the minimum energy configuration (stable equilibrium). The magnetic moment of the white dwarf rotates through $`\pi `$ radians as the polarity of the secondary star switches during its magnetic cycle. We assume that this adjustment occurs on a time scale that is shorter than the magnetic cycle time scale, and that a given system spends the bulk of its evolution with the dipoles nearly anti-aligned. In this context we note that in the few well studied AM Hers there is a preference for the dipole axes to be aligned with the spin axis . Changes have been observed in the rotational azimuth of the magnetic axis in three of AM Her systems on a time scale of $`25`$ years (see Wickramasinghe & Wu 1991, Wu & Wickramasinghe 1993) but no changes in magnetic polarity have so far been reported.

In what follows we make the explicit assumption that the braking rates that we deduce for anti-aligned systems are representative of the mean braking rate under which the AM Hers evolve.

We consider a coordinate system in which the secondary star is at the origin $`O`$, $`OZ`$ in the direction of the orbital angular momentum, and the primary star is in the negative $`y`$ direction at (0,$`a`$,0). The field components in the plane $`OYZ`$ containing the two dipole moments are

$$B_y=3\mu _2z\left[\frac{y}{r^5}+\frac{k_\mu (a+y)}{r_1^5}\right]$$

and

$$B_z=\mu _2\left[\left(\frac{1}{r^3}+\frac{3k_\mu }{r_1^3}\right)3z^2\left(\frac{1}{r^5}+\frac{3k_\mu }{r_1^5}\right)\right]$$

where

$$k_\mu =\frac{\mu _1}{\mu _2}$$

where $`r_1`$ and $`r`$ are the radial distances of a field point from the primary and the secondary stars respectively.

The dead zones are assumed to be isothermal with sound speed $`c_s^2(d)`$ and to corotate with the secondary star. The condition of magnetostatic equilibrium then yields the density within this region. For a point with polar coordinates ($`r,\theta `$) along a field line starting from the surface of the secondary ($`r=R`$, $`\theta =\theta _0`$), the density is \[1994a\]

$`\rho `$ $`=`$ $`\rho _0\mathrm{exp}[l_d(1{\displaystyle \frac{R}{r}})l_{d1}(1{\displaystyle \frac{r_1(R)}{r}})`$
$`+{\displaystyle \frac{1}{2}}\chi l_d(\mathrm{sin}^2\theta \mathrm{sin}^2\theta _0)+\beta ({\displaystyle \frac{r}{R}}\mathrm{sin}\theta \mathrm{sin}\theta _0)]`$
where $`\rho _0`$ at the density at the base of the dead zone. Here

$$\chi l_d=\frac{R^2\omega ^2}{c_s^2(d)},\beta =\frac{M_1a}{M_1+M_2}\frac{R^2\omega ^2}{c_s^2(d)}$$

and

$$ł_d=\frac{GM_2}{c_s^2(d)R},ł_{d1}=\frac{GM_1}{c_s^2(d)r_1(R)}$$

respectively measure the importance of centrifugal and gravitational effects of the two stars relative to thermal effects on the wind. We trace the magnetic field lines by solving the coupled differential equations

$$\frac{dy}{ds}=\frac{B_y}{B}$$

and

$$\frac{dz}{ds}=\frac{B_z}{B}$$

where $`B`$ is the total field. The critical field lines which separate the wind zones from the dead zones are determined by the pressure balance condition

$$\frac{B^2}{8\pi }=\rho c_s^2(d).$$

For assumed dead zone parameters, which we have evaluated as in Li et al. \[1994a\], these calculations yield the region $`\theta _1\theta \theta _2`$ of open field lines in the plane $`OYZ`$ which is then used to estimate the fractional area $`\mathrm{\Phi }`$ on the surface of open magnetic flux.

$`\mathrm{\Phi }`$ depends primarily on the magnetic moment of the primary star and of the mass of the secondary star, and only weakly on $`M_1`$. The dependence of $`\mathrm{\Phi }`$ on these two parameters is shown in Figure 3. We note in particular the rapid drop in open flux for large $`\mu _1`$ at fixed $`M_2`$ which plays a crucial role in determining the characteristics of the orbital evolution to be discussed in the next section.

As shown by Li et al. \[1994a\], the braking rate takes a particularly simple form when the wind is centrifugally driven. In this limit

$`\dot{J}`$ $`=`$ $`4.63\times 10^{35}\mathrm{dyn}\mathrm{cm}\left({\displaystyle \frac{\dot{J}_0}{4.63\times 10^{35}}}\right)`$
$`\times \left({\displaystyle \frac{5^\mathrm{h}}{P_{\mathrm{orb}}}}\right)^{5/3}\left({\displaystyle \frac{\mathrm{\Phi }}{0.258}}\right)^{5/3}\left({\displaystyle \frac{R_2}{0.5}}\right)^{10/3}`$
The same theory also provides an expression for $`\dot{J}_0`$ which involves wind and dead zone parameters. However, the dependence of these quantities on spectral type is unknown, and as discussed by Li & Wickramasinghe , a simple scaling from solar values is inappropriate. Following previous investigators, we have accordingly chosen to normalise the braking rate so as to reproduce the observed period gap in the limit of zero magnetic field. This normalisation corresponds to a braking rate of $`4.63\times 10^{35}`$ dyn cm for a Roche lobe filling secondary star in a non-magnetic ($`\mu _1=0`$) CV at an orbital period of $`5^\mathrm{h}`$ .

The mass of the white dwarf and its magnetic moment are essentially free parameters in our calculations as are the mass of the secondary and the factor $`f`$ which determines the radius of the donor star. We have constructed a grid of magnetic braking rates as a function of $`M_1,f,\mu _1,M_2`$, and used these to calculate orbital evolution for various initial component masses.

## 4 evolutionary calculations

We have explored the evolutionary consequences of attenuated magnetic braking using the bipolytrope code developed by Rappaport et al. , as modified by Nelson, Rappaport & Joss . This code idealizes the convective envelope-radiative core structure of near-main sequence stars of arbitrary degeneracy as $`n=3/2`$ polytropic envelopes fitted to $`n=3`$ polytropic cores, and reproduces the most important features of mass-losing lower main sequence stars with considerable fidelity. Except as noted below, we assume solar composition for the donor star throughout, with photospheric opacities adopted from Alexander  after Rappaport, Joss & Webbink .

For purposes of computing mass transfer rates, we fit isothermal atmospheres to the bipolytropes, and write the mass transfer rate explicitly in terms of the difference between the fitting point radius for this atmosphere and the Roche lobe radius of the donor star. However, in the calculations presented here, that difference never exceeds $`2\times 10^3`$ of the stellar radius, so we could have achieved the same results for all practical purposes by limiting the stellar radius to the Roche lobe radius as a surface boundary condition.

We assume that the white dwarf mass remains constant throughout the evolution of our model binaries, with all mass lost by the hydrogen-rich donor lost from the system in nova explosions. The ejected mass is assumed to carry off a specific angular momentum equal to the orbital angular momentum per unit mass of the white dwarf. Observationally, nova outbursts are found in non-AM Her systems and AM Her systems alike (*viz.* V1500 Cyg). The heavy-element enhancements typical of the ejecta of well-studied novae  are unmistakable evidence that the underlying white dwarfs do not generally grow in mass during the course of evolution, but may in fact suffer net erosion through nova outbursts. Provided that this erosion is not so severe as to drive the binary toward dynamical instability (roughly, that $`dM_1/dM_2<3/2`$), that the binary in question has not arrived at the onset its cataclysmic variable phase just above the period gap, and that the angular momentum content of nova ejecta has not been seriously underestimated, then the width and position of the period gap are insensitive to our assumption of constant white dwarf mass. Figure 4 illustrates this point, by comparing a family of evolutionary tracks for nonmagnetic CVs, differing in the extent of systemic mass loss, but converging to identical total masses at the upper edge of the period gap.

Figure 5 underscores the point implicit in Figure 4: for a given magnetic braking rate, the position and width of the period gap depends principally on the component masses of the binary as it ceases magnetic braking, and not on its evolutionary history in reaching that point. Since the disrupted magnetic braking model for the period gap associates the cessation of magnetic braking with the donor star becoming fully convective, the principal factors governing the position and width of the period gap (again, for a given magnetic braking rate) are (i) the mass of the white dwarf (or, equivalently, the mass ratio), and (ii) the chemical composition of the donor star.

The mass of the white dwarf influences the period gap through its role in orbital dynamics. Smaller white dwarf masses demand larger mass transfer rates to moderate the rate of Roche lobe contraction induced by angular momentum losses to magnetic braking. These higher mass transfer rates in turn drive the donor star further from thermal equilibrium, expanding the period gap, as seen in Figure 6. However, the effect is only strong for binary mass ratios near the limit of dynamical stability for fully convective donors ($`(M_1/M_2)_{\mathrm{crit}}=3/2`$), and cannot by itself account for CVs found in the gap.

Figure 7 illustrates the effect of varying metallicity of the donor star on the position and width of the period gap. Metal-deficient donors are systematically bluer and more luminous than their solar-composition counterparts of the same mass. Their higher luminosities endow them with shorter thermal time scales, and they therefore adhere more closely to the thermal equilibrium main sequence before turning fully convective than do solar-composition stars. As seen in Figure 6, then, metal-deficient stars encroach into the period gap from the long-period (low-frequency) end. We tentatively suggest, therefore, that the small number of non-AM Her-type CVs spilling into the period gap at long period may be relatively metal-poor. Metal-deficient main sequence stars are also physically smaller at low masses than solar composition stars, and so the short-period (high-frequency) boundary of the period gap migrates toward shorter periods (higher frequencies).

The preceding considerations demonstrate that, within the context of the disrupted magnetic braking model, the period gap is a robust entity. One can account for small numbers of non-AM Her CVs in the gap by appealing to newly-formed systems in which the donor stars are already fully convective main sequence stars with masses in the 0.2–0.3 M range, or by appealing to infiltration by metal-deficient donors at the long-period end of the gap. However, in light of the robustness of the period gap, the difference in orbital period distributions between AM Her systems and non-AM Her systems seen in Figure 1 strongly suggests that the AM Her systems deviate from the disrupted braking law which seems to apply to all other CVs.

As described in Section 3, in synchronous rotators (AM Her binaries), the dipole field of the white dwarf can close field lines from the donor star, suppressing the stellar wind on otherwise open field lines, and thereby suppress angular momentum losses. We have incorporated the attenuated magnetic braking model described there into the evolutionary model described above, under the assumption that the white dwarf magnetic dipole moment, like its mass, remains constant in time. Our model sequences assume that the white dwarf is locked in synchronism with the binary orbit from the very onset of mass transfer, i.e., that the system is of AM Her type throughout its evolution.

We show in Figure 8 the effect of increasing magnetic moment, $`\mu _1`$, of the white dwarf primary on the evolution of a typical AM Her system with initial component masses $`M_1=1.0`$ M for the white dwarf and $`M_2=0.6`$ M for the donor secondary. As the braking torque is suppressed with increasing $`\mu _1`$, mass transfer rates above the period gap subside, and the donor star is no longer driven so far from thermal equilibrium before it becomes fully convective at the top (lower frequency edge) of the period gap. As the donor star is no longer so severely inflated in comparison with its thermal-equilibrium main sequence counterpart, so its contraction following cessation of the magnetic braking driving mass transfer is also less severe. Mass transfer resumes at longer orbital period (lower orbital frequency) than in the disrupted magnetic braking model for non-AM Her-type CVs because the donor is no longer so severely undermassive for its orbital period, and the period gap fills from its shorter period (higher frequency) end. As seen in Figure 5, this filling of the period gap is progressively greater as $`\mu _1`$ increases, and magnetic braking prior to the secondary’s loss of its dipole field is suppressed. For sufficiently high values of $`\mu _1`$ (see $`\mu _1=10^{35}`$ G cm3 in Figure 5), magnetic braking becomes inconsequential, and the system evolves more or less continuously through the period gap (the momentary dip in mass transfer rate at $`P3^\mathrm{h}`$ owing to an inflection point in the thermal equilibrium mass-radius relation at the point stars become fully convective).

Figure 9 illustrates the dependence of the upper and lower edges of the period gap on white dwarf mass $`M_1`$ and magnetic moment $`\mu _1`$. We define the upper (lower-frequency) edge of the period gap as the point at which the donor star becomes fully convective. In those cases where magnetic braking remains strong enough to drive the donor out of thermal equilibrium, this point indeed marks a local minimum in the orbital period, which increases momentarily (if only by a very small increment) as the mass transfer rate subsides. We take the lower (higher-frequency) edge of the period gap as the period maximum (frequency minimum) during the thermal time scale mass transfer episode visible as a slim vertical loop in Figure 5. For sufficiently large white dwarf magnetic moments, where magnetic braking is nearly completely suppressed, this loop disappears, and we then define the lower edge of the period gap as the local maximum in $`\ddot{M}_2`$.

It is evident from Figure 9 that, for white dwarf moments $`\mu _110^{33}`$ G cm3, the gap edges asymptotically approach those of the non-AM Her systems. At the other extreme, for white dwarf moments $`\mu _14\times 10^{34}`$ G cm3, the gap is reduced to a momentary hiccup as the donor becomes fully convective. For intermediate white dwarf moments, the period gap fills from its short-period (high-frequency) end, as indeed observations suggest (see Figure 1). A clear prediction of our model is that, as more AM Her systems are discovered within the period gap, a clear trend of increasing magnetic moment with increasing orbital period (decreasing orbital frequency) should emerge.

To facilitate comparison with observed AM Her stars in the period gap, we have used the Nauenberg  mass-radius relation for cold white dwarfs to convert the dipole moments of Figure 9 to polar fields. The region immediately surrounding the period gap is shown in Figure 10. Unfortunately, empirical masses of white dwarfs in AM Her binaries are generally very poorly-constrained. However, to the extent that our model results can be taken at face value, the positions of observed systems in Figure 10 can be used to set upper limits to their white dwarf masses, since for given dipole moment, the polar field strength decreases rapidly with increasing mass. In this figure, AM Her binaries must fall *outside* the gap corresponding to their white dwarf masses, thus implying that QS Tel, V381 Vel, V516 Pup, and UW Pic must contain white dwarfs of mass $`M_11.1`$ M. The reader is cautioned, however, that this constraint is quite sensitive to the precise location of the short-period (high-frequency) edge of the canonical period gap.

## 5 Discussion and conclusions

The AM Her stars, unlike cataclysmic variables generally, have mostly been discovered through all-sky surveys, where they appear as strong X-ray emitters. We can therefore approximate the selection effects governing their discovery by treating them as an X-ray-brightness limited sample, with X-ray luminosities proportional to the mass transfer rate $`|\dot{M}_2|`$. To the extent that X-ray extinction is negligible, the distance to which systems are sampled is then proportional to $`|\dot{M}_2|^{1/2}`$. If the threshold distance is small compared with the scale height of AM Her systems with respect to the galactic plane, then the apparent distribution should be roughly isotropic, and the number of systems sampled varies as $`|\dot{M}_2|^{3/2}`$; should the distance be large compared with the scale height (but still small compared with the radial scale length of the galactic disk) the number of systems varies as $`|\dot{M}_2|`$ itself. The number of systems in any orbital frequency interval is then simply proportional to the lifetime of the system in that interval, i.e., to $`\dot{\nu }_{\mathrm{orb}}^1`$. As a point of reference, the mean value of $`\mathrm{sin}|b|`$ of the AM Her systems in Figure 1 is 0.51, indistinguishable from that of a spherical source distribution ($`\mathrm{sin}|b|=0.50`$).

Figure 11 illustrates the relative discovery probabilities arising from this simple model of selection effects, in the two limiting cases regarding their spatial distribution, for our fiducial evolutionary sequence. In reality, we expect CVs to form with a range of component masses, and this will affect the overall form of the true distribution. In particular, systems with a given initial secondary mass, $`M_{2\mathrm{i}}`$, contribute to the distribution only at orbital frequencies $`\nu _{\mathrm{orb}}`$ above the frequency at which the CV first reaches interaction (cf. Figure 5). Likewise, magnetic CVs will contribute to the AM Her distribution only at orbital frequencies above those at which they first reach synchronism. Nevertheless, as seen above, the properties of individual systems at the upper and lower edges of the period gap, and the positions of those edges, are insensitive to their initial conditions.

One feature immediately evident from Figure 11 is that the number density per unit orbital frequency of CVs at the upper (low-frequency) edge of the period gap is decidedly greater than at its lower (high-frequency) edge. The ratio of these two number densities, as defined at the very edges of the period gap, amounts to a factor of 5.4 for non-AM Her systems with a spherical source distribution (1.9 in a planar distribution). Among the AM Her systems of increasing magnetic moment $`\mu _1`$, the amount by which this ratio exceeds unity is roughly proportional to the width of the period gap appropriate to that value of $`\mu _1`$.

Reference to Figure 1 shows that the non-AM Her systems have a much higher number density of systems at the upper edge of the period gap than at its lower edge, in at least qualitative agreement with Figure 11. However, among AM Her systems themselves, this ratio of number densities is plainly less than unity. The period gap is for them only partially filled, and the marked increase in number density beyond the lower (high-frequency) edge of the gap should be matched by an equally marked increase at the upper (low-frequency) edge. The clear implication is that more AM Her-type systems emerge from the period gap than enter it. Thus, a significant fraction of AM Her systems below the period gap must have entered the gap in a non-AM Her (asynchronous) state, i.e., as intermediate polars. As asynchronous systems, they suffered full magnetic braking above the gap, and were driven out of thermal equilibrium to the same extent as other non-AM Her systems. They therefore remain detached until resuming mass transfer at the lower edge of the gap as defined by other non-AM Her systems, in transit becoming synchronously locked, but failing to appear either as CVs or AM Her systems within the gap. Systems which achieve synchronism while crossing the period gap thus contribute to the AM Her population only at periods below the gap.

Comparison of Figure 11 with Figure 1 also shows that the theoretical period minimum of AM Her stars lies at significantly shorter period (higher frequency) than observed. This is a problem pervading CV evolution generally, since both AM Her systems and non-AM Her systems are driven through this period minimum by gravitational radiation. Kolb & Baraffe  discuss a number of possible explanations for this discrepancy. Among these possibilities, they note that a four-fold increase in orbital angular momentum loss rate (above that given by gravitational radiation alone) would raise the CV minimum period to that observed (if we discount the anomalous cases of V485 Cen at $`P_{\mathrm{orb}}`$ = 0$`\stackrel{h}{.}`$984 and 1RXS J2329+0628 at $`P_{\mathrm{orb}}`$ = 1$`\stackrel{h}{.}`$070). However, we find that at these ultashort orbital periods, the white dwarf magnetic field should suppress a magnetic stellar wind entirely among AM Her systems. The fact that both AM Her binaries and non-AM Her binaries are subject to the same minimum period (Figure 1) then suggests that stellar winds are not responsible for this discrepancy.

We summarise our principal conclusions as follows:

(a) The width of the period gap, and the position of its lower boundary, is a function of the magnetic moment of the primary star. For sufficiently high magnetic moments ($`\mu _14\times 10^{34}`$ G cm3), the gap disappears, and mass transfer is driven primarily by gravitational radiation.

(b) Observations of AM Her-type systems show that the period gap is less pronounced than that of CVs belonging to all other classes. In particular, the AM Hers appear to fill the gap from the short-period ($`2^\mathrm{h}`$) end, behaviour which supports the reduced braking theory outlined in section 3. In contrast, non-AM Her systems encroach on the period gap from its long-period edge, behaviour which we suggest may reflect lower metallicities of these particular systems.

(c) The differences between the distributions of AM Her binaries and non-AM Her binaries within the period gap in turn provide strong prima facie evidence that magnetic stellar winds are indeed the dominant angular momentum loss mechanism among non-AM Her binaries.

(d) The abundance of AM Hers immediately below the period gap, in comparison with their numbers immediately above the period, requires that a significant proportion of the AM Hers found below $`2^\mathrm{h}`$ entered the period gap as IPs, and become synchronised only as mass transfer ceased within the gap. The smaller accretion torques, and increased magnetic interaction (due to the reduction in orbital separation) allowed these systems to remain synchronised as they emerge from the gap.

## 6 Acknowledgements

RFW thanks the Astrophysical Theory Centre, Department of Mathematics, Australian National University, where this work was begun, for its hospitality, and Lilia Ferrario for her generous assistance. DTW in turn acknowledges the hospitality of the Department of Astronomy, University of Illinois, where this work was completed. Both authors thank an anonymous referee for constructive remarks. Studies of interacting binaries are supported at the University of Illinois by National Science Foundation grant AST 9618462.


text:
# The crossover between lasing and polariton condensation in optical microcavities

## Abstract

We study a model of a photon mode dipole-coupled to a medium of two-level oscillators in a microcavity in the presence of dephasing processes introduced by coupling to external baths. Decoherence processes can be classified as pair-breaking or non-pair-breaking in analogy with magnetic or non-magnetic impurities in superconductors. In the absence of dephasing, the ground state of the model is a polariton condensate with a gap in the excitation spectrum. Increase of the pair-breaking parameter $`\gamma `$ reduces the gap, which becomes zero at a critical value $`\gamma _{C1}`$; for large $`\gamma `$, the conventional laser regime is obtained in a way that demonstrates its close analogy to a gapless superconductor. In contrast, weak non-pair-breaking processes have no qualitative effect on the condensate or the existence of a gap, although they lead to inhomogeneous broadening of the excitations.

There are two apparently distinct phenomena where quantum coherence has been observed on the macroscopic scale: the laser and Bose Einstein condensation (BEC) of massive particles.

The laser is a weak-coupling phenomenon: a coherent state of photons created by stimulated emission from an inverted electronic population due to the strong pumping. The polarisation of the medium is heavily damped and the atomic coherence is very much reduced. A coherent photon field, oscillating at bare cavity mode frequency, is the only order parameter in the system .

Polaritons are strongly coupled modes of light and electronic excitations. A theoretically constructed polariton condensate is a mixture of coherent state of light and coherent state of massive particles in the media. It is characterised by two order parameters: the coherent polarisation and the coherent photon field and exhibit a gap in the excitation spectrum .

The simplest models of both the laser and the polariton condensate start from similar equations of motion, but make different physical approximations. The laser is assumed to be far from equilibrium, and decoherence (dephasing) — arising from collisions, interactions with phonons, defects or impurities, and manifested by homogeneous broadening ($`1/T_2`$) of the electronic excitations — is imposed by hand. In contrast, the polariton condensate has been assumed to be in thermal equilibrium, and dephasing phenomena neglected.

In this letter, we show that a consistent inclusion of decoherence phenomena drives a crossover from the polariton condensate to the lasing regime. Using a model of two-level oscillators coupled to a degenerate optical cavity mode with coupling strength $`g`$, and with further coupling of the internal degrees of freedom to external baths (defined in the Hamiltonian below), we find the generic phase diagram of Fig. 1, where a single pair-breaking parameter $`\gamma `$ encapsulates the strength of electronic dephasing. When the decoherence is weak, we find a strong-coupling regime with a polariton condensate and a gap in the excitation spectrum. The gap is robust against (weak) dephasing processes, a result that can be obtained only by the correct self-consistent treatment of decoherence processes. For larger decoherence, the gap is suppressed to zero at a critical value $`\gamma _{C1}`$ in a similar way to the effect of pair-breaking processes in a BCS superconductor, and the coherent polarisation is reduced. For low excitation densities coherence is completely suppressed at a higher value $`\gamma _{C2}`$. In the limit of high excitation densities (photon dominated) and large decoherence, we recover the traditional model of the laser.

The principal aim of this paper is to establish the correct theoretical framework to discuss the differences and similarities between a condensate and a laser, so we adopt a deliberately simplified model that cannot accommodate all the details of important physical systems. Nevertheless, in closing we attempt estimates of the relevant parameters for some relevant physical systems.

To obtained the phase diagram we have considered the following Hamiltonian: $`H=H_0+H_{SB}+H_B`$, where $`H_0=\omega _c\psi ^{}\psi +_{j=1}^Nϵ_j(b_j^{}b_ja_j^{}a_j)+\frac{g_j}{\sqrt{N}}(b_j^{}a_j\psi +\psi ^{}a_j^{}b_j)`$. $`H_0`$ describes an ensemble of N two-level oscillators with an energy $`ϵ_j`$ dipole coupled to one cavity mode. $`b`$ and $`a`$ are fermionic annihilation operators for an electron in an upper and lower states respectively and $`\psi `$ is a photon bosonic annihilation operator. The sum here is over the possible sites where an exciton can be present (different molecules or localised states in the disorder potential). $`H_{SB}`$ contains all interactions with an environment:

$$\begin{array}{c}H_{SB}=\underset{k}{}g_k^\kappa (\psi ^{}d_k+d_k^{}\psi )+\underset{j,k}{}[g_{j,k}^\gamma _{}(b_j^{}a_jc_{j,k}^\alpha +c_{j,k}^\alpha a_j^{}b_j)\hfill \\ \hfill +g_{j,k}^\gamma _{}(b_j^{}a_jc_{j,k}^\beta +c_{j,k}^\beta a_j^{}b_j)+g_{j,k}^{\gamma _1}(b_j^{}b_j+a_j^{}a_j)(c_{j,k}^\zeta +c_{j,k}^\zeta )\\ \hfill +g_{j,k}^{\gamma _2}(b_j^{}b_ja_j^{}a_j)(c_{j,k}^\theta +c_{j,k}^\theta )]\end{array}$$
(1)
and the $`H_B`$ is a Hamiltonian for the baths. The first term in (1) gives rise to the decay of the photon field from the cavity. The second term describes an incoherent pumping of two-level oscillators. The third term contains all the processes which destroy the electronic excitations such as the decay to photon modes different to the cavity mode. These processes, apart from a dephasing effect, will cause a flow of energy through the system. However, in the steady state the total number of excitations, $`n_{ex}=\psi ^{}\psi +\frac{1}{2}_j(b_j^{}b_ja_j^{}a_j)`$, which is the sum of photons and excited two-level oscillators, is conserved.

Finally, the fourth and the fifth terms describe all the dephasing processes which do not change the total number of excitations in the cavity, for example collisions and interactions with phonons and impurities. They can be divided into a part, which acts in the same way (fifth term) and with an opposite sign (fourth term) on the upper and lower levels. These four terms contain all the essential groups of decoherence processes in the media.

Processes, described by the second, the third and the fourth terms in (1), are analogous to the pair-breaking, magnetic impurities in superconductors and correspond to potentials which vary rapidly in space (on the lengthscale of excitons) or in time. On the contrary the fifth term contains interactions which vary very slowly in space in comparison to the size of excitons and in time and are analogous to normal, non-magnetic, impurities. They do not have any pair-breaking effects and lead only to an inhomogeneous broadening of energies.

For very big pumping and photon decay rates the system would be out of equilibrium. However, if the system is pumped slowly, with a rate slower than the thermalisation rate equilibrium assumption can be justified. In this work we present only the influence of the decoherence on the system leaving the implications of a non-equilibrium system for the future work. This does not put any restrictions on the strength of the decoherence nor the density of excitations. The pair-breaking processes described by the fourth term in (1) give exactly the same decoherence effects as the pumping and damping of two-level oscillators preserving equilibrium for any strength of the interaction. The ratio between the pumping and all the damping processes which can be arbitrary big even for small absolute values of both will be expressed in terms of the excitation density $`\rho _{ex}`$.

Usually the baths are averaged out before the exciton - photon interaction is studied which leads to the well-known quantum Maxwell-Bloch (Langevin) equations for the photon field and polarisation. When there is a coherent polarisation this procedure is not valid ( Langevin equations can be used to study the laser regime but not the polariton or intermediate regime) . Instead, we must use the self-consistent Green’s function techniques similar to the Abrikosov and Gor’kov theory of gapless superconductivity . Decoherence processes are treated by perturbation theory about a BCS - like Green function for the interacting exciton-photon system. The major difference between our calculations and the Abrikosov and Gor’kov theory  is that we have two order parameters, the coherent photon field and polarisation connected through the chemical potential, set by the excitation density. We consider two-level systems in contrast to free propagating electrons model used by Abrikosov and Gor’kov. Within the Markov approximation we can express the influence of the bath by a single parameter $`\gamma _l=g^{\gamma _l}(0)^2N_l(0)`$ where $`l=1,2`$ correspond to the pair-breaking and non-pair breaking decoherence and $`N_l`$ are densities of states for the respective baths. The details of the method are published elsewhere .

We first determine the ground state properties of the system with uniform energies, $`ϵ_j=ϵ`$, for different densities in the presence of the pair-breaking decoherence. Figure 2 shows the ground state behaviour of the coherent fields and inversion. For small values of $`\gamma /g`$, up to some critical value $`\gamma _{C1}`$, $`\psi `$ is practically unchanged while for $`\gamma /g>\gamma _{C1}`$ is damped quite rapidly with the increasing dephasing (see Fig. 2 \- upper panel). This critical value of the decoherence strength, $`\gamma _{C1}`$ is proportional to $`\rho _{ex}`$, suggesting that for higher excitation densities the system is more resistant to the dephasing. At low excitation densities, where $`\rho _{ex}<0`$, there is a second critical value of the decoherence strength, $`\gamma _{C2}`$, where both coherent fields are sharply damped to zero.

The ratio of coherent polarisation to coherent field $`P/\psi `$ (Fig. 2 \- lower panel) for an isolated system, where $`\gamma =0`$, depends on the excitation density. The condensate becomes more photon like as $`\rho _{ex}`$ is increased due to the phase space filling effect. For finite $`\gamma `$, at a given excitation density, this ratio decreases with increasing $`\gamma `$ meaning that $`P`$ is more heavily damped than $`\psi `$.

To understand the behaviour of coherent fields we

study the excitation spectrum of the system (Fig. 3). As $`\gamma `$ increases the two quasi-particle peaks (Fig. 3 a) broaden, which causes the decrease in the magnitude of the energy gap (Fig. 3 b and c). Finally, precisely at $`\gamma _{C1}`$ (shown in Fig. 2), these two broadened peaks join together and the gap closes (Fig. 3 d). $`\gamma _{C1}`$ defines the phase boundary (see Fig 1) between a gapped and gapless condensate at low densities and between a condensate and laser at high densities.

The laser operates in the regime of a very strong decoherence, comparable with the light-matter interaction itself. In the presence of such a large decoherence laser action can be observed only for a sufficiently large excitation density . The coherent polarisation in a laser system is much more heavily damped than the photon field and a gap in the density of states is not observed. Thus the laser is a regime of our system for very large $`\rho _{ex}`$ and $`\gamma `$. It has to be, however, pointed out that there is no formal distinction between a laser and a phase-coherent condensate of polaritons. In the laser, coherence in the media (manifested by the coherent polarisation) although small, is not completely suppressed so the laser can be seen as a gapless condensate with a more photon-like character. One of the possible distinctions between a polariton condensate and laser could be an existence of an energy gap in the excitation spectrum.

The dotted lines in onset of Fig. 2 and in Fig. 3 show the influence of the non-pair breaking processes (the fifth term in (1)). We found that these processes have some quantitative influence on the coherent fields and the gap in the density of states but do not cause any phase transitions. Although the two quasiparticle peaks get very broad, the gap is only slightly affected. Even for much larger values of $`\gamma `$ than presented in Fig. 3 the gap is present until the coherent fields get completely suppressed. This type of processes give similar effect as the inhomogeneous broadening of energy levels.

Finally, we study an influence of the pair-breaking processes on the system of realistic, inhomogeneously broadened two-level oscillators. The coherent fields, $`\gamma _{C1}`$, $`\gamma _{C2}`$ and the energy gap are slightly smaller than in the uniform case but all the regimes and transitions are the same. The broadening of the density of states and the suppression of the energy gap can be observed as $`\gamma `$ is increased (Fig. 4).

General conclusions about the conditions and systems where polariton condensation is likely to be observed can be drawn from our theory. Systems with large dipole coupling $`g`$, and high densities but low decoherence are good candidates. Inhomogeneous broadening of energies have much weaker influence than a homogeneous dephasing. The large value of $`g`$ is particularly important as it defines the energy scale in this problem. Only the relative strength of the dephasing and inhomogeneous broadening with respect to $`g`$ influence the condensate. Thus systems with larger $`g`$ can allow larger values of decoherence and broadening. The second important factor is the excitation density. Since the gap is proportional to the amplitude of the coherent field the condensate is more robust at high densities. Unfortunately the increase in the density usually implies an increase in decoherence as the dephasing time ($`T_2`$) for most of materials increases linearly with density. An additional increase in dephasing comes from the increase in the pump intensity required to achieve higher densities. This could be partially overcome by improving the quality of the mirrors in microcavity systems.

To be more specific we consider in detail different materials. Among semiconductors CdTe, ZnSe and GaN are characterised by large values of $`g`$ (measured by the polariton splitting in the normal state). For CdTe, the $`g`$ was measured to be 29 meV  and the dephasing ($`1/T_2`$), at the exciton density of 0.05 per Bohr radius, to be 2.5 meV (0.08 g) . In the absence of inhomogeneous broadening our theory predicts a gap $`\mathrm{\Delta }=0.7g=20.3`$ meV. In the case of Gaussian broadening with $`\sigma =0.2g=5.8`$ meV reported in some structures the gap is still significant $`\mathrm{\Delta }=0.08g=2.32`$ meV and it gets suppressed for broadening larger then 8.7 meV. Other promising candidates would be excitons in organics, where the reporting polariton splitting was as large as 80 meV , or dilute atomic gases in which the dephasing is measured to be around 1000 times smaller then the dipole coupling  and high excitations can be achieved without an increase in the density of the gas. GaAs — due to the small value of $`g`$ — is not the best material for the observation of the condensate. Polariton splitting in GaAs based structures is usually only around 9 meV  and the recently reported 20 meV splitting in 36 well structure  is probably approaching the absolute limit for these materials.

In this work we have developed a method to include decoherence effects in the systems with an energy gap in the excitation spectrum, where the widely used Langevin equations are not valid. We have studied an equilibrium, incoherently pumped system. However, in the same way it would be possible to examine the influence of an environment on coherently driven condensates. We have shown that in the presence of small decoherence the polariton condensate is protected by an energy gap in the excitation spectrum. The gap is proportional to the coherent field amplitude and thus the excitation density, so the condensate is more robust at high densities. As the decoherence is increased the gap gets smaller and finally is completely suppressed. If the polariton condensate is ever to be observed an increase in density without an increase in decoherence is necessary. Our method allows us to study different regimes in the cavity and thus to examine the crossover from a polariton condensate to a laser. The laser emerges from the polariton condensate at high densities when the gap in the density of states closes for large decoherence and thus is analogous to a gapless superconductor. Our work generalises the existing laser theories to include coherence effects in the media. When this coherence is included the generation of a coherent photon field without population inversion becomes possible.