text:
# Theoretical Study of Electronic Structure and Superconductivity in 𝑁𝑏_{1-𝑥}𝐵₂ Alloys

## 1 INTRODUCTION

The discovery of superconductivity in $`MgB_2`$ and the fact that the carbides and the nitrides of $`Nb`$ show superconductivity have renewed interest in finding superconductivity in the borides of $`Nb`$, in particular, the diborides. However, the search for superconductivity in $`NbB_2`$ and $`NbB_2`$-based alloys is not new. Earlier work by Zeigler $`etal`$  and Hulm $`etal`$  showed $`B`$ deficient $`NbB_2`$, namely $`NbB_{1.94}`$ to be superconducting with superconducting transition temperature, $`T_c`$, of around $`1.2K`$, while Schriber $`etal`$  found $`NbB_x`$ (synthesized under pressure) with $`x2`$ to be superconducting with $`T_c`$ of $`9.2K`$.

Recently, Yamamoto $`etal`$  and Kotegawa $`etal`$  have reported superconductivity in $`Nb_{1x}B_2`$ alloys when synthesized under pressure for $`x`$ ranging from $`0`$ to $`0.48`$. They find that $`Nb_{1x}B_2`$ samples become superconducting for $`x=0.04`$, and the $`T_c\text{ }`$ increases to around $`9.2K`$ for $`x=0.24`$, thereafter $`T_c`$ remains essentially unchanged up to $`x=0.48`$.

In an effort to theoretically understand the changes in the electronic structure of $`Nb_{1x}B_2`$ and $`NbB_{1.94}`$ alloys and their superconducting properties as vacancies are introduced either in the $`Nb`$ plane or the $`B`$ plane, we have carried out first-principles electronic structure calculations of these alloys as a function of vacancy concentration. We have used Korringa-Kohn-Rostoker coherent-potential approximation  in the atomic-sphere approximation (KKR-ASA CPA) method for taking into account the effects of disorder, Gaspari-Gyorffy formalism  for calculating the electron-phonon coupling constant $`\lambda `$, and Allen-Dynes equation  for calculating $`T_c`$ in $`Nb_{1x}B_2`$ and $`NbB_{1.94}`$ alloys.

Here, we like to point out that the recent theoretical work of Shein $`et`$ $`al`$  to model disordered $`Nb_{0.75}B_2`$ alloy using a 12-atom supercell may not be adequate. This can be seen as follows. Each metal atom in $`AlB_2`$-type diborides has $`D_{6h}`$ symmetry and thus has 12 $`B`$ neighbors at the vertices of a hexagonal prism. In addition, the central metal atom also sees $`8`$ others through the faces of the $`B_{12}`$ prism at a distance comparable with the metallic diameter. Thus a minimum of $`24`$-atom supercell would have been more appropriate. However, their conclusion that the superconductivity in $`Nb_{1x}B_2`$ alloys is governed more by the changes that take place in the lattice properties when vacancies are introduced in the metal plane of the diborides, is consistent with our results. Before describing our results, we provide some of the computational details of the present approach.

## 2 Computational Details

### 2.1 The densities of states

The charge self-consistent electronic structure of $`Nb_{1x}B_2`$ and $`NbB_{1.93}`$ alloys has been calculated using the KKR-ASA CPA method with $`x`$ ranging from $`0`$ to $`0.48.`$ The lattice constants for $`Nb_{1x}B_2`$ alloys are taken from the experimental work of Yamamoto $`etal`$. To determine the electronic structure of $`NbB_2`$ using the present approach, we have used the lattice constants $`a`$ and $`c`$ from our earlier full-potential linear muffin-tin orbital calculations. The absence of $`Nb`$ and $`B`$ atoms in the $`Nb`$ and $`B`$ planes, respectively, were modelled by introducing empty spheres at those sites. The disorder was accounted for by the coherent-potential approximation , which has been very successful in describing reliably many physical properties of disordered alloys . The exchange-correlation potential was parametrized as suggested by Perdew-Wang within the generalized gradient approximation . The Brillouin zone (BZ) integration was carried out with $`1215`$ k-points in the irreducible part of the BZ. For densities of states (DOS) calculations, we added a small imaginary component of $`1`$ $`mRy`$ to the energy and used $`4900`$ k-points in the irreducible part of the BZ. For valence states the angular momentum cut-off was kept as $`l_{max}=3`$, and the core states were recalculated. For $`Nb`$ and $`B`$ the ratio of the atomic sphere radius to the Wigner-Seitz radius was kept fixed at 1.294 and 0.747, respectively. The empty spheres inherited the radii of the corresponding atomic spheres. The present calculations were carried out using the scalar-relativistic Schroedinger equation. The Green’s function is calculated in a complex plane with a $`20`$ point Gaussian quadrature.

### 2.2 The electron-phonon coupling and the superconducting transition temperature

The superconducting transition temperature $`T_c`$ was calculated with the Allen-Dynes equation,

$$T_c=\frac{\omega _{ln}}{1.2}\mathrm{exp}\left[\frac{1.04(1+\lambda )}{\lambda \mu ^{}(1+0.62\lambda )}\right].$$
(1)
In our calculation, the value of the phonon frequency $`\omega _{ln}`$ for $`NbB_2`$ was taken from Ref. , $`\mu ^{}`$ was set equal to 0.09, and the electron-phonon coupling constant $`\lambda `$ was calculated using the Gaspari-Gyorffy formalism with the charge self-consistent potentials of $`Nb_{1x}B_2`$ and $`NbB_{1.93}`$ alloys obtained with the KKR-ASA CPA method. This entails writing $`\lambda `$ in terms of Hopfield parameter, $`\eta `$, as

$$\lambda =\underset{i}{}\frac{\eta _i}{M_i<\omega _i^2>},$$
(2)
where the sum is over the basis atoms in the primitive cell, $`M`$ is the atomic mass and $`<\omega ^2>^{1/2}`$ is an average phonon frequency. The average phonon frequency $`<\omega ^2>^{1/2}`$ is taken to be the root-mean- square frequency $`\omega _{rms}=<\omega ^2>^{1/2}`$. In the Gaspari-Gyorffy formalism the spherically averaged Hopfield parameter is given by

$$\eta =2\underset{l}{}\frac{(l+1)}{(2l+1)(2l+3)}M_{l.l+1}^2\frac{N_l(E_F)N_{l+1}(E_F)}{N(E_F)}.$$
(3)
The total DOS at the Fermi energy $`E_F`$, $`N(E_F)`$, and the $`l`$resolved DOS, $`N_l(E_F)`$ as well as the electron-phonon matrix elements,

$$M_{l,l+1}=\varphi _l(E_F)\varphi _{l+1}(E_F)\left[(D_l(E_F)l)(D_{l+1}(E_F)+l+2)+(E_FV(S))S^2\right]$$
(4)
are calculated from the charge self-consistent potentials obtained with the KKR-ASA CPA method. In Eq. (4), $`\varphi _l(E_F)`$ is the amplitude of the $`l`$th partial wave at the sphere boundary, $`S,`$ evaluated at $`E_F,`$ $`D_l(E_F)`$ is the corresponding logarithmic derivative, and $`V(S)`$ is the one-electron potential at $`S`$.

## 3 Results and Discussion

In the following we describe the results of our calculations for $`NbB_2`$ and its alloys in terms of (i) the densities of states, (ii) the electron-phonon coupling and (iii) the superconducting transition temperature. We first describe how the densities of states, the electron-phonon coupling and the superconducting transition temperature change for low concentrations of vacancies in either the $`Nb`$ or the $`B`$ plane. Then, we describe the changes in these properties as a function of increasing vacancy concentration.

### 3.1 The densities of states of $`NbB_2`$, $`Nb_{0.99}B_2`$ and $`NbB_{1.93}`$

#### 3.1.1 $`NbB_2`$

In Fig. 1 we show the calculated total DOS and the partial DOS together with the $`l`$-resolved partial DOS of $`NbB_2`$. The overall structure of the DOS is consistent with the previous work. The calculations based on the tight-binding LMTO method yields a value of total DOS at the Fermi energy equal to $`4.65`$ states/Ry-cell , while the corresponding value from the full-potential LMTO calculations is found to be $`4.28`$ states/Ry-cell. The present calculations, which is based on the KKR-ASA method, gives $`N(E_F)`$ to be equal to $`4.41`$ states/Ry-cell. At $`E_F`$ the contributions of $`Nb`$ and $`B`$ to the total DOS is found to be $`3.44`$ and 0.97 states/Ry-cell, respectively.

The characteristic features in the DOS of $`NbB_2`$, as shown in Fig. 1, are $`(i)`$ a low lying broad peak around $`0.900.75`$ Ry below $`E_F`$ due to $`B`$ $`s`$ and $`Nbd`$states, $`(ii)`$ another broad peak at around $`0.450.25`$ Ry comprising of $`Bp`$, $`Nbp`$ and $`Nbd`$ states, $`(iii)`$ a relatively constant DOS from $`0.70.5`$ Ry primarily of $`Bp`$ and $`Nbd`$ states, and $`(iv)`$ a pseudo-gap at $`0.14`$ Ry, which appears in many transition-metal diborides . Thus the dominant presence of $`Nbd`$ states at $`E_F`$ in $`NbB_2`$ makes it quite different from $`MgB_2`$ which has largely $`B`$ $`p`$states at $`E_F`$.

#### 3.1.2 $`Nb_{0.99}B_2`$

To understand how low concentrations of vacancies in the $`Nb`$ plane of $`NbB_2`$ would affect its electronic properties, we have used KKR-ASA CPA method to calculate the electronic structure of dilute $`Nb_{1x}B_2`$ with $`x`$$`=`$$`0.99`$ alloy. In this calculation the lattice constants were kept the same as that of $`NbB_2`$. In Fig. 2, we show the calculated $`l`$-resolved partial DOS at the empty sphere site in $`Nb_{0.99}B_2`$. It is seen that upon $`1`$% creation of vacancy in the $`Nb`$ plane, the total DOS at $`E_F`$ for $`Nb_{0.99}B_2`$ is reduced by $``$$`4`$% relative to that of $`NbB_2`$. We find that the creation of vacancies in the $`Nb`$ plane do not lead to appreciable change in the $`B`$ partial DOS at $`E_F`$. Thus the metal deficiency introduced in $`NbB_2`$ enhances the ratio of $`B`$ states to that of the total DOS at $`E_F`$.

#### 3.1.3 $`NbB_{1.93}`$

Motivated by the work of Otani $`etal`$ on $`NbB_{1.93}`$ alloy, which showed that a depletion of $`B`$ in $`NbB_2`$ led to an increase in the lattice parameter $`c`$ from $`6.16`$ to $`6.25`$ a.u. and a decrease in the lattice parameter $`a`$ from $`5.89`$ to $`5.84`$ a.u., we have calculated the electronic structure of $`NbB_{1.93}`$ alloy. The calculation is carried out by disordering the $`B`$ sub-lattice with empty spheres and keeping the $`Nb`$ sub-lattice ordered. The calculated densities of states at $`E_F`$ are summarized in Table I, where we have also listed the lattice constants used in the present work.

Our analysis of the DOS of $`NbB_{1.93}`$ shows that the creation of vacancies in the $`B`$ plane largely reduces the $`d_{3z^21}`$ component ($``$12.5%) of the $`Nb`$ $`d`$-states and the $`B`$ $`s`$-states ($``$8.9%) at $`E_F`$. The creation of vacancies in the $`B`$ sub lattice of $`NbB_2`$ thus leads to a major reduction in those states which are directed along the $`c`$axis of the unit cell.

### 3.2 The electron-phonon coupling and the superconducting transition temperature of $`NbB_2`$, $`Nb_{0.99}B_2`$ and $`NbB_{1.93}`$ alloys

The superconducting transition temperatures of $`NbB_2`$, $`Nb_{0.99}B_2`$ and $`NbB_{1.93}`$ were calculated with the Allen-Dynes equation. The factor representing the Coulomb interaction $`\mu ^{}`$ was 0.09. For $`\omega _{rms}`$ and $`\omega _{ln}`$, we used the values $`464cm^1`$ and $`494K`$, respectively. Our results for the Hopfield parameter $`\eta `$, electron-phonon coupling constant $`\lambda `$ and the superconducting transition temperature are summarized in Table II. From Table II, we see that $`\eta _B`$, the Hopfield parameter due to $`B`$, is relatively small in comparison to $`\eta _{Nb}`$, which reflects the fact that $`B`$ contribution to the total density of states at $`E_F`$ is small. Since we use the same $`\omega _{rms}`$ for $`NbB_2`$, $`Nb_{0.99}B_2`$ and $`NbB_{1.93}`$ the changes in $`\lambda `$ are similar to that of $`\eta `$. The calculated $`T_c`$ for $`NbB_2`$, $`Nb_{0.99}B_2`$ and $`NbB_{1.93}`$ are found to be $`0.31`$K, $`0.13`$K and $`0.004`$K, respectively. In the case of $`NbB_{1.93}`$ alloy, we find that the deficiency of $`B`$ largely affects $`\eta _{Nb}`$ rather than $`\eta _B`$, which reduces $`T_c`$ significantly.

### 3.3 The densities of states of $`Nb_{1x}B_2`$ alloys

In Fig. 3, we show the DOS of $`Nb_{0.96}B_2`$, $`Nb_{0.84}B_2`$, $`Nb_{0.76}B_2`$ and $`Nb_{0.52}B_2`$ alloys using the KKR-ASA CPA method. The lattice constants for these alloys were taken from Ref. . The changes in the DOS that arise in $`Nb_{1x}B_2`$ alloys as a function of vacancies are $`(i)`$ the inward movement of $`E_F`$ with increasing vacancy concentration, leading to a decrease in the total DOS at $`E_F`$ for $`x<0.16`$ and an increase for $`x>0.16`$, and $`(ii)`$ the broadening of peaks with increasing vacancy concentration due to disorder. We find that the movement of $`E_F`$ follows the trend as suggested by the rigid-band model, and the DOS of $`Nb_{0.84}B_2`$ alloy is similar to $`ZrB_2`$.

### 3.4 The electron-phonon coupling and the superconducting transition temperature of $`Nb_{1x}B_2`$ alloys

From the resulting electronic structure we have calculated the variation in the Hopfield parameter, electron-phonon coupling constant and the superconducting transition temperature of $`Nb_{1x}B_2`$ alloys as a function of vacancy concentration. Our results for $`\eta `$ and $`\lambda `$ are summarized in Table III. We find that both $`\eta `$ and $`\lambda `$ decrease with increasing vacancy concentration up to $`x0.25`$, thereafter they start to increase. In Fig. 4, we show the calculated $`T_c`$ as a function of $`x`$ in $`Nb_{1x}B_2`$ alloys, together with the experimentally reported values. The parameters used in the calculation of $`T_c`$ were $`\mu ^{}=0.09`$, $`\omega _{rms}=400cm^1`$, and $`\omega _{ln}=494K`$. The calculated trend in the variation of $`T_c`$ is in disagreement with the experimental results as shown in Fig. 4. In our opinion, the disagreement is very likely due to the use of a constant $`\omega _{rms}`$ value in the Allen-Dynes equation for the whole range of $`x`$ from $`0`$ to 0.48 for $`Nb_{1x}B_2`$ alloys.

In the present approach, the calculation of $`T_c`$ incorporates two important approximations, (i) use of Gaspari-Gyorffy formalism for estimating the electron-phonon matrix elements, which is known to underestimate the strength of the coupling, and (ii) use of a constant phonon frequency over the whole range of $`x`$. Using Gaspari-Gyorffy formalism with concentration-weighted averages of the phonon frequencies, we have been able to describe reliably the trend in the variation of $`T_c`$ in many $`MgB_2`$-based alloys . Since the electron-phonon coupling in $`NbB_2`$ is dominated by $`Nb`$ and not by $`B`$, as is the case in $`MgB_2`$, the Gaspari-Gyorffy formalism may be more suitable for $`NbB_2`$ than for $`MgB_2`$. This can be seen by comparing the $`T_c`$ as calculated in the present approach with that of the linear response calculation of Ref. . The lack of agreement between the present calculations and the observed $`T_c`$ in $`Nb_{1x}B_2`$ alloys as a function of $`x`$ indicates the importance of incorporating the lattice-dynamical effects in these calculations.

## 4 Conclusions

We have studied the electronic structure of $`NbB_2`$, $`NbB_{1.93}`$ and $`Nb_{1x}B_2`$ alloys using the KKR-ASA CPA method. We find that the vacancies on the metal plane reduce the number of $`Nb`$ $`d`$-states at $`E_F`$. The changes in the density of states in $`Nb_{1x}B_2`$ alloys as a function of $`x`$ reflect a rigid-band picture. The variation in $`T_c`$ in $`Nb_{1x}B_2`$ alloys as a function of $`x`$, calculated using the Gaspari-Gyorffy formalism disagrees with the observed trend. In our opinion, this disagreement arises due to the use of a constant $`\omega _{rms}`$ in the Allen-Dynes equation over the whole concentration range in our calculations.


format: text

text:
# Wigner’s little group and Berry’s phase for massless particles

## Abstract

The “little group” for massless particles (namely, the Lorentz transformations $`\mathrm{\Lambda }`$ that leave a null vector invariant) is isomorphic to the Euclidean group E2: translations and rotations in a plane. We show how to obtain explicitly the rotation angle of E2 as a function of $`\mathrm{\Lambda }`$ and we relate that angle to Berry’s topological phase. Some particles admit both signs of helicity, and it is then possible to define a reduced density matrix for their polarization. However, that density matrix is physically meaningless, because it has no transformation law under the Lorentz group, even under ordinary rotations.

Eugene Wigner considered his paper “On unitary representations of the inhomogeous Lorentz group”  as his most important contribution to physics . The key feature in that article was the introduction of a little group, namely a subgroup under which a standard vector $`s_\mu `$ is invariant. For example, a timelike $`s_\mu `$ is (1,0,0,0) and the little group is the familiar rotation group SO(3). A null standard vector can be taken as (1,0,0,1). Spacelike standard vectors have no physical interest.

To find explicitly the little group that corresponds to the null $`s_\mu `$, let us introduce an auxiliary complex null vector $`m_\mu `$ such that

$$m_\mu s^\mu =m_\mu m^\mu =0,$$
(1)
$$m_\mu ^{}m^\mu =1,$$
(2)
and a real null vector $`n_\mu `$ that satisfies

$$n_\mu m^\mu =0,$$
(3)
$$n_\mu s^\mu =1.$$
(4)
All these properties are manifestly Lorentz invariant. Moreover they still hold under the transformation

$$m_\mu e^{i\xi }m_\mu +\gamma s_\mu ,$$
(5)
$$n_\mu n_\mu +e^{i\xi }\gamma ^{}m_\mu +e^{i\xi }\gamma m_\mu ^{}+|\gamma |^2s_\mu ,$$
(6)
where $`\gamma =\alpha +i\beta `$ is any complex number. Therefore the above transformation of $`m_\mu `$ and $`n_\mu `$ is a subgroup of the Lorentz group. We thereby obtain a representation of the generic little group element

$$g=S(\alpha ,\beta )R_z(\xi ),$$
(7)
where the $`S(\alpha ,\beta )`$ form a subgroup which is isomorphic to the translations of a Euclidean plane and $`R_z(\xi )`$ is a rotation around the origin of that plane, which in this case is also a rotation around the $`z`$-axis. We see that the little group is E2, as shown in various ways by other authors .

In the case of infinitesimal transformations in Minkowski space, we have

$$g=\mathrm{𝟏}+\alpha A+\beta B+\xi M_{12},$$
(8)
where $`M_{12}J_3`$, and $`A`$ and $`B`$ are commuting generators of translations in the 12-plane, for example ,

$$A=M_{01}+M_{31}\text{and}B=M_{02}+M_{32}.$$
(9)
In general, let $`\mathrm{\Lambda }_\nu ^\mu `$ denote the Lorentz transformation matrix and $`k^\nu `$ an arbitrary momentum. The task is to find the little group matrix $`W_\nu ^\mu `$ that corresponds to $`\mathrm{\Lambda }_\nu ^\mu `$ and $`k^\nu `$, namely

$$W(\mathrm{\Lambda },k)=L^1(\mathrm{\Lambda }k)\mathrm{\Lambda }L(k),$$
(10)
where spacetime indices were omitted for brevity, and $`L(k)`$ is the standard Lorentz transformation (next equation) that converts $`s`$ to an arbitrary momentum $`k`$. For example, in the case of massive particles for which the little group is SO(3), if $`\mathrm{\Lambda }`$ is an ordinary rotation then $`W`$ is the same rotation irrespective of $`k`$. However if $`\mathrm{\Lambda }`$ is a boost, then $`W`$ does depend on $`k`$. Likewise the various terms in (7) may depend on $`k`$.

In this Letter, we shall examine the group properties of massless particles. For $`s=(1,0,0,1)`$, the standard Lorentz transformation, as defined in , is

$$L(k)=R(\widehat{𝐤})B_z(|𝐤|),$$
(11)
where $`B_z(|𝐤|)`$ is a boost along the $`z`$-axis with velocity $`(1𝐤^2)/(1+𝐤^2)`$, and $`R(\widehat{𝐤})`$ is the standard rotation that carries the $`z`$-axis into the direction of the unit vector $`\widehat{𝐤}`$. Again following , if the direction of $`\widehat{𝐤}`$ is given by spherical angles $`\theta `$ and $`\varphi `$, the standard rotation $`R(\widehat{𝐤})`$ consists of a rotation $`\theta `$ around the $`y`$-axis, followed by a rotation $`\varphi `$ around the $`z`$-axis:

$$R(\widehat{𝐤})=R_z(\varphi )R_y(\theta ).$$
(12)
In these formulas, all Lorentz transformations are passive, even if we occasionally use an active wording.

We now prove by a classical geometric argument that, if the Lorentz transformation is a pure rotation ($`\mathrm{\Lambda }=`$), then it follows from Eq. (10) that $`S(\alpha ,\beta )`$ simply is the unit matrix, and therefore $`\alpha =\beta =0`$. We also give a simple expression for the rotation angle $`\xi `$. From the definition of the little group element we have

$$W(,𝐤)=B_z^1(|𝐤|)R^1(\widehat{𝐤})R(\widehat{𝐤})B_z(|𝐤|).$$
(13)
Since the action of $`R^1(\widehat{𝐤})R(\widehat{𝐤})`$ leaves the $`z`$-axis invariant, it is equivalent to some rotation $`R_z(\varpi )`$ around that axis,

$$W(,𝐤)=B_z^1(|𝐤|)R_z(\varpi )B_z(|𝐤|)=R_z(\varpi ),$$
(14)
and since $`W(,k)`$ is a special case of (10), it follows that in Eq. (7) we have $`S_\nu ^\mu =\delta _\nu ^\mu `$ and $`\alpha =\beta =0`$. Equation (7) also gives $`\xi =\varpi `$, and the remaining problem is to find the explicit value of this angle.

To clarify the origin of the phase $`\xi `$ in Eq. (7), we note that any rotation in a three-dimensional space can be described by two angles that give the direction of the rotation axis, and a third angle that gives the amount of rotation around that axis. A rotation from $`𝐤`$ to $`𝐪=𝐤`$ can be performed in many ways (denoted below by $`R_{\mathrm{𝐪𝐤}}`$), in addition to the given $``$ that we are seeking to decompose. Since all such rotations satisfy

$$𝐤=R_{\mathrm{𝐪𝐤}}𝐤,^1𝐪=R_{\mathrm{𝐪𝐤}}^1𝐪,$$
(15)
the difference between them is a rotation that preserves $`\widehat{𝐪}`$, if done after $`R_{\mathrm{𝐪𝐤}}`$, or a rotation that preserves $`\widehat{𝐤}`$, if done before $`R_{\mathrm{𝐪𝐤}}`$. In particular, $`𝐪=R_{\mathrm{𝐪𝐤}}^1𝐪`$, so that

$$R_{\mathrm{𝐪𝐤}}^1=R_{\widehat{𝐪}}(\omega ),$$
(16)
where $`R_{\widehat{𝐪}}(\omega )`$ is a rotation around $`\widehat{𝐪}`$. Among the infinity of possible $`R_{\mathrm{𝐪𝐤}}`$ we choose

$$R_{\mathrm{𝐪𝐤}}=R(\widehat{𝐪})R^1(\widehat{𝐤}),$$
(17)
where $`R(\widehat{𝐪})`$ and $`R(\widehat{𝐤})`$ are standard rotations, as in Eq. (12). It follows that

$$=R_{\widehat{𝐤}}(\omega )R(\widehat{𝐤})R^1(\widehat{𝐤}),$$
(18)
where $`R_{\widehat{𝐤}}(\omega )`$ is a rotation around $`\widehat{𝐤}`$, while $`R(\widehat{𝐤})`$ and $`R(\widehat{𝐤})`$ are the standard rotations that carry the $`z`$-axis to $`\widehat{𝐤}`$ and $`\widehat{𝐤}`$, respectively. We can thus consider Eq. (16) as the definition of $`R_{\widehat{𝐤}}(\omega )`$.

Substituting this decomposition into Eq. (14), we obtain

$$W(,k)=R^1(\widehat{𝐤})R_{\widehat{𝐤}}(\omega )R(\widehat{𝐤})=R_z(\varpi ),$$
(19)
and we conclude that $`\xi =\varpi =\omega `$.

To obtain the rotation angle under a general Lorentz transformation, we decompose the latter into two rotations and a standard boost $`B_z`$ along the $`z`$-axis :

$$\mathrm{\Lambda }=_2B_z(u)_1.$$
(20)
As shown below, $`B_z`$ alone does not lead to a phase rotation. Therefore,

$$\xi =\omega _1+\omega _2,$$
(21)
where both $`\omega _1`$ and $`\omega _2`$ are due to the rotations and are given by Eq. (18).

We now prove that $`B_z`$ alone induces no phase rotation . Consider a pure boost along the $`z`$-axis, $`\mathrm{\Lambda }=B_z(u)`$, and a generic null vector $`k=(|𝐤|,𝐤)`$, where

$$𝐤=|𝐤|(\mathrm{sin}\theta \mathrm{sin}\varphi ,\mathrm{sin}\theta \mathrm{cos}\varphi ,\mathrm{cos}\theta ).$$
(22)
We define $`q=B_z(u)k=(|𝐪|,𝐪)`$ where

$$𝐪=|𝐪|(\mathrm{sin}\theta ^{}\mathrm{sin}\varphi ,\mathrm{sin}\theta ^{}\mathrm{cos}\varphi ,\mathrm{cos}\theta ^{}).$$
(23)
Note that the angle $`\varphi `$ is the same for $`𝐤`$ and for $`𝐪`$. Thus

$$R^1(\widehat{𝐪})=R_y^1(\theta ^{})R_z^1(\varphi ).$$
(24)
We now substitute Eqs. (12) and (24) into Eq. (10), to obtain

$$W(B_z(u),k)=B_z^1(|𝐪|)R_y^1(\theta ^{})B_z(u)R_y(\theta )B_z(|𝐤|),$$
(25)
where we used $`R_z^1(\varphi )B_z(u)R_z(\varphi )=B_z(u)`$.

Consider now the effect of the little group element (25) on the spacelike vector $`y=(0,0,1,0)`$. That vector is not affected by a boost in the $`z`$ direction, nor by a rotation around the $`y`$-axis. Therefore

$$W(B_z(u),k)y=y,$$
(26)
so that in this case $`\xi `$ is either $`0`$ or $`2\pi `$. Since for $`u=0`$ we expect $`\xi =0`$, by continuity $`\xi =0`$ for all $`u`$.

Note that although $`B_z(u)`$ alone does not lead to a phase rotation, it can affect the value of $`\omega _2`$, since it indirectly appears in the definition of $`_2`$. Indeed, if we decompose $`_2`$ as in Eq. (18), we obtain

$$_2=R_{_2\widehat{𝐤}_2}(\omega _2)R(_2\widehat{𝐤}_2)R^1(\widehat{𝐤}_2),$$
(27)
where $`𝐤_2`$ is defined by

$$k_2=(|𝐤_2|,𝐤_2)=B_z(u)_1k.$$
(28)
Thus we see that $`B_z(u)`$ appears in the decomposition of $`_2`$ and therefore affects $`\omega _2`$.

Up to this point, the discussion and the formalism were purely classical. In quantum theory, one needs the unitary representations of the little group, from which those of the complete Lorentz group can be derived. Each irreducible representation corresponds to some species of elementary particles. According to Eq. (10), the general transformation law is

$$U(\mathrm{\Lambda })|𝐤,\sigma =\underset{\sigma ^{}}{}D_{\sigma ^{}\sigma }[W(\mathrm{\Lambda },k)]|𝐪,\sigma ^{},$$
(29)
where $`D_{\sigma ^{}\sigma }`$ is a unitary representation of the little group and

$$|𝐤,\sigma |𝐤|\sigma $$
(30)
is an appropriate basis. The helicity $`\sigma =𝐉\widehat{𝐤}`$ of a massless particle is Lorentz invariant, so that if we use it for labelling basis states, then the sum in Eq. (29) consists of a single term, and

$$D_{\sigma ^{}\sigma }=\mathrm{exp}(i\xi \sigma )\delta _{\sigma ^{}\sigma },$$
(31)
where, for a Lorentz transformation which is a pure rotation, $`\xi `$ is a function of $``$ and $`k`$ which is explicitly given by Eq. (18).

It is experimentally known that some particles, like neutrinos (if they are indeed massless), come with only one sign of helicity. Others, like photons, may have it with both signs, but then the phase angle $`\xi `$ for them is different. In general, different values of $`|\sigma |`$ refer to different species of particles, such as photons and gravitons. Within the present formalism, we cannot offer an explanation why, for a given species of particles, half-integral helicities appear with a definite sign. Of course, there can be no helicity-statistics theorem, since we deal with a single particle. These properties must follow from quantum field theory for interacting fields.

An application of the above results is a direct derivation of the Berry phase for massless particles. Soon after the introduction of Berry’s phase  the latter was derived for photons in the adiabatic approximation  and then  for arbitrary changes in momentum. Finally, derivations that are based on the analysis of connections on Lie groups, and the Poincaré group in particular, were given in .

Consider a sequence of rotations that eventually restores the particle momentum to its original value. Its net effect is some active rotation around the momentum’s direction, $`R_{\widehat{𝐤}}(\omega )`$. According to Eqs. (19) and (31), the helicity eigenstates acquire phases $`\omega \sigma `$, where the minus sign arises from the fact that the transformation law (31) is for passive rotations. To relate this phase to the area on the unit sphere that is enclosed by the orbit of $`\widehat{𝐤}`$, consider an auxiliary unit vector $`\widehat{𝐯}`$ in the plane perpendicular to $`\widehat{𝐤}`$. It is tangent to the sphere at the endpoint of $`\widehat{𝐤}`$. Let the orthonormal triad $`\widehat{𝐤}`$, $`\widehat{𝐯}`$, and $`\widehat{𝐰}=\widehat{𝐤}\times \widehat{𝐯}`$ be parallel-transported along that orbit. When the latter closes, the orthonormal triad does not return to itself, but ends up as another triad at the same point, which is rotated by the angle $`\omega `$ in the $`\widehat{𝐯}\widehat{𝐰}`$ plane. Owing to properties of the holonomy group , the rotation angle $`\omega `$ is related to the spherical angle $`\mathrm{\Omega }`$ by

$$\omega =K𝑑S=𝑑\mathrm{\Omega }=\mathrm{\Omega },$$
(32)
where $`K`$ is the Gaussian curvature, which equals 1 for the unit sphere. The above integral is over the area enclosed by the trajectory of $`\widehat{𝐤}`$, and the sign of $`\mathrm{\Omega }`$ depends on the orientation of the trajectory. We thus obtain

$$\omega =\mathrm{\Omega }.$$
(33)
As another application, consider the concept of reduced density matrices  which is fundamental in quantum information theory. Their properties are significantly modified by relativistic effects . For massless particles that admit both signs of helicity, such as photons, a generic one-particle state is

$$|\mathrm{\Psi }=𝑑\mu (k)\underset{\sigma }{}f_\sigma (𝐤)|𝐤,\sigma ,$$
(34)
where $`d\mu (k)=d^3𝐤/(2\pi )^3(2|𝐤|)`$ is a Lorentz-invariant measure. Then the reduced density matrix for helicity, according to the usual rules, would be

$$\rho _{\sigma \tau }=𝑑\mu (k)f_\sigma (𝐤)f_\tau ^{}(𝐤).$$
(35)
However, since $`\xi `$ in Eq. (31) depends on the photon momentum (even for ordinary rotations) the standard density matrix given by Eq. (35) has no transformation rule at all. This makes the standard density matrix a useless concept, even when only a fixed reference frame is considered, since any POVM that describes an experimental setup must have definite transformation properties at least under ordinary rotations. It is only possible to define an “effective” density matrix which depends on the detection method . This behavior contrasts with that of massive particles, for which the little group is SO(3) and reduced density matrices behave properly under rotations, while there is no transformation law only under boosts .

The absence of any Lorentz transformation law for $`\rho `$ is due to the fact that the momenta $`k`$ transform linearly, but the law of transformation of helicity depends explicitly on $`k`$. When we compute $`\rho `$ by summing over momenta, all knowledge of them is lost and it is then impossible to obtain the new $`\rho `$ by transforming the old one. There is an analogous situation in classical statistical mechanics: a Liouville function can be defined in any Lorentz frame , but it has no definite transformation law from one frame to another. Only the complete dynamical system has a transformation law .

In summary, we have shown how apparently disparate notions — Wigner’s little group and Berry’s phase — are closely related. It is curious that the proof made repeated use of ordinary rotations, namely the SO(3) group, which is by itself another little group of the Lorentz group.

Work by AP was supported by the Gerard Swope Fund and the Fund for Promotion of Research.


text:
# On the detection of pre-low-mass X-ray binaries

## 1 Introduction

Ongoing deep X-ray surveys of selected Galactic regions have sparked a renewed interest in low-luminosity Galactic X-ray sources such as wind-accreting neutron stars. The XMM–Newton Galactic Plane Survey (e.g. Warwick 2002), the XMM-Newton Serendipitous Survey (Watson 2001; Motch et al. 2002; Motch, Herent & Guillout 2003; Watson et al. 2003) and the serendipitous Chandra Multi-wavelength Plane Survey (Grindlay et al. 2003) will all probe the faint Galactic X-ray point source population.

The discovery by Chandra of more than 500 previously undetected point sources with luminosities above $`10^{35}`$ erg/s in a $`2.0^{}\times 0.8^{}`$ field near the Galactic centre (Wang, Gotthelf & Lang 2002) dramatically illustrates the potential of such surveys and has already prompted complementing theoretical studies. Pfahl, Rappaport & Podsiadlowski (2002) considered the population of neutron stars that accrete from the wind of an intermediate- or high-mass star (mass $`3M_{}`$). They found that these wind-driven intermediate-mass and high-mass X-ray binaries may indeed account for a dominant fraction of the point source population found by Wang et al. (2002).

Independently, Bleach (2002) investigated the possibility of detecting detached close binaries consisting of a neutron star and a red dwarf, which may be the progenitors of low-mass X-ray binaries (LMXBs). The proposed detection mechanism consists of looking for red dwarfs with X-ray luminosities larger than those expected from intrinsic X-ray emissions from the star’s corona. The excess could then possibly be ascribed to the emission of X-rays by a neutron star companion accreting mass from the wind of the red dwarf.

Our aim in this paper is to elaborate on the idea by Bleach (2002) and examine the detection probability of pre-LMXBs, i.e. detached neutron-star systems with low-mass companions (mass $`2M_{}`$), in the context of a binary population synthesis study.

## 2 pre-LMXBs and LMXBs

In close binary systems of stars, mass can be transferred from one component to the other when one of the stars loses mass in a stellar wind or when one of the stars fills its Roche lobe. If the companion star is a compact object which is able to accrete some of the transferred mass, the liberation of gravitational potential energy gives rise to an accretion luminosity

$$L_{\mathrm{acc}}=\frac{GM\dot{M}_{\mathrm{acc}}}{R},$$
(1)
where $`G`$ is the Newtonian constant of gravitation, $`M`$ and $`R`$ are the mass and the radius of the accreting compact object, and $`\dot{M}_{\mathrm{acc}}`$ is the mass-accretion rate. We assume that the accretion process yields an X-ray luminosity $`L_\mathrm{X}L_{\mathrm{acc}}`$.

In the particular case of a wind-accreting neutron star, the spin evolution of the neutron star and the possible inhibition of accretion by strong centrifugal forces may play an important role in the study of the accretion process (e.g. Pringle & Rees 1972; Illarionov & Sunyaev 1975; Stella, White, Rosner 1986; Urpin, Geppert & Konenkov 1998). The treatment of these inhibiting effects is however beyond the scope of this initial investigation. Instead, we assume that at the onset of mass transfer via the stellar wind, the neutron star’s radiation pressure and rotation rate is low enough for mass accretion to be possible.

The evolution of binaries containing a neutron star is driven by the nuclear evolution of the companion star or by the loss of orbital angular momentum via magnetic braking and/or gravitational radiation. When the combined effect of these processes causes the companion to become larger than its Roche-lobe, an X-ray binary may be formed. The further evolution of the binary depends critically on the orbital period and the mass of the donor star at the onset of Roche-lobe overflow. The formation and evolution of X-ray binaries with low-mass donor stars have been discussed extensively in the past by, e.g., Verbunt (1993); Verbunt & van den Heuvel (1995); Iben & Tutukov (1995); Kalogera & Webbink (1996, 1998); Kalogera (1998); Kalogera, Kolb & King (1998); and more recently by Podsiadlowski, Rappaport & Pfahl (2002). We therefore refer to these studies and references therein for detailed discussions on the various formation channels and here consider the population of X-ray binaries descending from the different channels as a whole.

For the purpose of this investigation, we define a candidate pre-LMXB as a detached binary in which a neutron star accretes mass from the stellar wind of its companion. We assume the latter to be on the main sequence, in the Hertzsprung gap, or on the giant branch, and to have a mass smaller than or equal to $`2M_{}`$. The luminosity associated with the accretion process is required to be larger than $`10^6L_{}`$ or ten times the mass donor’s coronal X-ray luminosity, whichever is larger. We note that our definition of pre-LMXBS differs from that of Bleach (2002) in that we consider stars with masses up to $`2M_{}`$ and binaries that are expected to become semi-detached due to the nuclear evolution of the secondary as well as those that are expected to become semi-detached due to angular momentum losses via magnetic braking and/or gravitational radiation.

Similarly, we define a LMXB as a semi-detached system in which a neutron star accretes mass from a Roche-lobe filling main-sequence, Hertzsprung-gap, or giant-branch star with a mass smaller than or equal to $`2M_{}`$. In order to be classified as an X-ray binary, we furthermore require the mass-accretion rate to be high enough to produce an X-ray luminosity in excess of $`10^1L_{}`$.

For the remainder of the paper, we denote the mass and the radius of the neutron star by $`M_{\mathrm{NS}}`$ and $`R_{\mathrm{NS}}`$, and the mass and the radius of the donor star by $`M_\mathrm{d}`$ and $`R_\mathrm{d}`$, respectively.

## 3 Assumptions and computational technique

We explore the possibility of detecting and identifying LMXBs prior to the onset of Roche-lobe overflow from the neutron star’s companion by means of the BiSEPS binary population synthesis code described by Willems & Kolb (2002). The code is based on the analytic approximation formulae for the evolution of single stars derived by Hurley, Pols & Tout (2000) and follows the main lines of the binary evolution algorithm described by Hurley, Tout & Pols (2002). The binary orbits are assumed to be circular and the stellar rotation rates are kept synchronised with the orbital motion at all times. We furthermore restrict ourselves to Population I stellar compositions.

For the purpose of this investigation, the mass-loss rates from stellar winds given by Hurley et al. (2000), which are limited to evolved stars and massive main-sequence stars, must be supplemented by a prescription for the winds from low-mass main-sequence stars. Although these winds are not fully understood yet, it is generally accepted that they depend sensitively on both the mass and the age of the donor star (e.g. Wood et al. 2002, and references therein). Young solar-type stars and M dwarfs, for instance, are thought to have mass-loss rates which may be as high as $`10^{12}M_{}`$ per year, while older G dwarfs such as the Sun have mass-loss rates of the order of $`10^{14}M_{}`$ per year (Hartmann 1985, Lim & White 1996, van den Oord & Doyle 1997, Wargelin & Drake 2001, Sackmann & Boothroyd 2003). In view of these still existing uncertainties, we here simply adopt a constant mass-loss rate of $`10^{13}M_{}`$ per year for low-mass main-sequence stars (see also Bleach 2002). The resulting mass-loss rates are shown in Fig. 1 for donor stars with masses up to $`2M_{}`$ and evolutionary stages from the zero-age main sequence up to the tip of the giant branch. The lines labelled ZAMS, TMS, and BGB represent the radii of the stars at the zero-age main sequence, the terminal main sequence, and the base of the giant branch, respectively. The maximum evolutionary age considered for the calculation of the mass-loss rates is 10 Gyr, so that only stars with masses larger than $`1M_{}`$ significantly evolve away from the zero-age main sequence.

We estimate the mean mass-accretion rate onto a star orbiting in the stellar wind of its companion by means of a standard Bondi-Hoyle-Lyttleton formalism (Hoyle & Lyttleton 1941, Bondi & Hoyle 1944). The accretion rate for binaries with circular orbits then depends on the orbital separation $`a`$, the binary mass ratio $`q`$, the donor star’s radius $`R_\mathrm{d}`$, the wind mass-loss rate $`\dot{M}_\mathrm{d}`$, and the wind velocity $`v_\mathrm{w}`$. We parametrise the latter as a fraction $`\beta _\mathrm{w}`$ of the escape velocity $`v_\mathrm{e}=\left(2GM_\mathrm{d}/R_\mathrm{d}\right)^{1/2}`$ at the surface of the mass-losing star (see, e.g., Hurley et al. 2002). The mean wind mass-accretion rate can then be cast in the form

$$\dot{M}_{\mathrm{acc}}=\frac{3}{16}\left(\frac{R_\mathrm{d}}{a}\right)^2\frac{q^2}{\beta _\mathrm{w}^4}\left(1+\frac{1+q}{2\beta _\mathrm{w}^2}\frac{R_\mathrm{d}}{a}\right)^{3/2}\left|\dot{M}_\mathrm{d}\right|.$$
(2)
In the particular case of a wind-accreting neutron star the mass ratio is given by $`q=M_{\mathrm{NS}}/M_\mathrm{d}`$.

The detection mechanism for pre-LMXBs proposed by Bleach (2002) furthermore requires a prescription for the intrinsic X-ray luminosity of the neutron star’s companion. Since this luminosity is thought to have its origin in dynamo processes generated by the interaction between convection and rotation, it can be expected to be correlated with the star’s rotational angular velocity (e.g. Belvedere, Chiuderi & Paterno 1982; Noyes et al. 1984). The existence of such a correlation was explored extensively from an observational point of view by, e.g., Pallavicini et al. (1981); Marilli & Catalano (1984); Vilhu (1984); Doyle (1987); Fleming, Gioia & Maccacaro (1989); Hempelmann et al. (1995); Randich et al. 1996; Singh et al. (1999); Gondoin (1999). The numerous investigations however show a clear lack of consensus on whether or not the correlation exists and, if it exists, what the precise relationship is between the rotational angular velocity and the coronal X-ray luminosity. This deficiency may be related to the use of different stellar samples and observational selection criteria (Hempelmann et al. 1995, Gondoin 1999).

A point of lesser dispute seems to be the existence of an upper limit on the intrinsic stellar X-ray luminosity, although its origin remains to be unambiguously identified (Walter 1982, Byrne et al. 1984, Vilhu & Walter 1987, Jardine & Unruh 1999, Singh et al. 1999). The upper limit may be related to the finite size of the stellar surface in the sense that the magnetic activity can only increase with increasing rotational angular velocities until the whole stellar surface is covered with magnetically active regions. Based on this idea, Fleming et al. (1989) showed that the coronal saturation luminosity $`L_{\mathrm{cor}}`$ could be related to the square of the stellar radius $`R_\mathrm{d}`$ as

$$\mathrm{log}L_{\mathrm{cor}}=2.9+2\mathrm{log}R_\mathrm{d},$$
(3)
where both $`L_{\mathrm{cor}}`$ and $`R_\mathrm{d}`$ are expressed in solar units (see their Fig. 3).

Rather than face the uncertainties regarding the relation between stellar rotation and coronal X-ray activity, we here opt to use a prescription for the upper limit on the coronal X-ray luminosity. We use the mass of the convective envelope, $`M_{\mathrm{conv}}`$, to distinguish between coronally active and non-active stars and, somewhat arbitrarily, set the limit for coronal X-ray activity at $`M_{\mathrm{conv}}=0.01M_\mathrm{d}`$. Stars with smaller convective envelopes are assumed to have negligible coronal X-ray luminosities, while stars with larger convective envelopes are assumed to have coronal X-ray luminosities given by Eq. (3). The upper limit may significantly overestimate the coronal X-ray luminosity of giant-type stars with slow rotation rates and correspondingly weak dynamos (Schröder, Hünsch & Schmitt 1998; Gondoin 1999). The slow rotation, however, only applies to giants in wide binaries where tidal forces are negligible. Giants in closer binaries are subjected to tidal torques which tend to synchronise the giant’s rotation with the orbital motion of the companion. A significant dynamo action may therefore be sustained in these stars during much more evolved stages of stellar evolution than for single stars or stars in wide binaries. Dempsey et al. (1993, 1997) have shown that in the particular cases of RS Canum Venaticorum (RS CVn) and BY Draconis (BY Drac) binaries the coronal X-ray luminosity of the giant-type component may still reach values up to $`0.1L_{}`$. To be on the safe side, we therefore still use the upper limit given by Eq. (3) for stars that have evolved beyond the main-sequence. Where relevant, we will discuss the possible implications of this on our results. In any case, we expect the consistent use of the upper limit for the coronal X-ray luminosity to imply that our selection criteria for the identification of pre-LMXBs are systematically too strict.

## 4 Population synthesis

We evolved a large number of binaries with initial component masses ranging from $`0.1`$ to $`60M_{}`$ and with initial orbital periods between 0.1 to 10 000 days. We used logarithmically spaced grids consisting of 60 grid points for the initial stellar masses and 200 grid points for the initial orbital periods. The maximum evolutionary age considered was 10 Gyr. If a binary undergoes an asymmetric supernova explosion leading to the birth of a neutron star, 200 random kick velocities are generated from a Maxwellian velocity distribution for the magnitude of the kick velocities. The evolution of the binaries surviving the explosion is then followed until the imposed age limit of 10 Gyr.

The outcome of the binary evolution calculations depends on the parameters adopted for the various evolutionary processes a binary may be subjected to. We therefore vary the input parameters one at a time and compare the resulting populations of pre-LMXBs with that of a reference model. For the latter, hereafter referred to as model A, we adopt the same set of input parameters as in Willems & Kolb (2002), i.e. the common-envelope ejection efficiency parameter takes the value $`\alpha _{\mathrm{CE}}=1.0`$, and the Maxwellian kick-velocity dispersion takes the value $`\sigma _{\mathrm{kick}}=190`$ km/s. In addition, we here adopt a standard value of 1.0 for the wind-velocity parameter $`\beta _\mathrm{w}`$. The different parameters adopted in the other population synthesis models are summarised in Table 1. In the case of the different supernova-kick models, the average kick velocities $`\overline{v}_{\mathrm{kick}}`$ adopted in the Maxwellian velocity distributions are listed in addition to the velocity dispersions $`\sigma _{\mathrm{kick}}`$. We also note that the variations in the wind-velocity parameter $`\beta _\mathrm{w}`$ considered in models MA1 and MA2 are equivalent to considering variations in the wind mass-loss rate $`\dot{M}_\mathrm{d}`$, but that the parameter $`\beta _\mathrm{w}`$ is a more appropriate population synthesis parameter since it has a greater impact on the determination of the mean mass-accretion rate. This is readily illustrated by means of Eq. (2) which shows that, to a first approximation, an increase in $`\beta _\mathrm{w}`$ by a factor of $`2`$ corresponds to a decrease in $`\dot{M}_\mathrm{d}`$ by a factor of the order of $`10`$. Models COR1, COR2, and COR3 are introduced to assess the sensitivity of our results to the prescription for the coronal X-ray luminosity. In particular, the use of a constant value for $`L_{\mathrm{cor}}`$ in models COR1 and COR2 allows us to account for the possible effects of overestimating the coronal X-ray luminosity of giant-type stars by the use of Eq. (3).

The contribution of a system to the population of pre-LMXBs depends on the probability distribution functions of its initial parameters and on the time the system spends as a member of the population. We assume the initial mass $`M_1`$ of the primary to be distributed according to the initial mass function

$$\xi \left(M_1\right)=\{\begin{array}{cc}0\hfill & M_1/M_{}<0.1,\hfill \\ 0.38415M_1^1\hfill & 0.1M_1/M_{}<0.75,\hfill \\ 0.23556M_1^{2.7}\hfill & 0.75M_1/M_{}<\mathrm{},\hfill \end{array}$$
(4)
the initial mass ratio $`q=M_2/M_1`$ according to

$$n(q)=\{\begin{array}{cc}1\hfill & 0<q1,\hfill \\ 0\hfill & q>1,\hfill \end{array}$$
(5)
and the initial orbital separation $`a`$ according to

$$\chi (a)=\{\begin{array}{cc}0\hfill & a/R_{}<3\text{ or }a/R_{}>10^4,\hfill \\ 0.12328a^1\hfill & 3a/R_{}10^4\hfill \end{array}$$
(6)
(see Willems & Kolb 2002, and references therein).

Each time a binary from our initial grid of parameters evolves into a pre-LMXB, its weighted contribution is added to the probability density function (PDF) describing the distribution of newborn pre-LMXBs in the $`(t,M_\mathrm{d},P_{\mathrm{orb}},L_{\mathrm{acc}})`$-space at the time $`t=t_b`$ of its formation. The finite lifetime of the pre-LMXBs is taken into account by subtracting the same weighted contribution from the distribution function at the time $`t=t_d`$ when the system ceases to be a pre-LMXB. In doing so, we implicitly assume that the stellar and orbital parameters of a pre-LMXB do not change significantly once it is formed. We will discuss the circumstances under which this assumption breaks down and the possible implications thereof in the final section of the paper. For now, we remark that this assumption does not affect the determinations of the birthrates and the total number of systems obtained after integration over all system parameters.

The resulting distribution function for the population of pre-LMXBs is subsequently convolved with a constant star-formation rate and normalised so that the integral over all systems found is equal to one.

## 5 Donor star masses and orbital periods

The normalised distribution function of pre-LMXBs in the $`(\mathrm{log}M_\mathrm{d},\mathrm{log}P_{\mathrm{orb}})`$-plane is displayed in Fig. 2 for the most illustrative of the population synthesis models listed in Table 1. Bins containing systems that actually evolve into LMXBs within the imposed age limit of 10 Gyr are outlined with a black solid line.

There are three reasons why systems satisfying our pre-LMXB criteria do not show up in the corresponding LMXB population. Donor stars with mass $`M_\mathrm{d}1M_{}`$ evolve too slowly to fill their Roche lobe as a result of their nuclear evolution within the imposed age limit of 10 Gyr. They can therefore only become LMXBs if their orbital period is shorter than the bifurcation period of $`1`$ day, which separates diverging from converging systems (Pylyser & Savonije 1988, 1989). Donor stars with mass $`M_\mathrm{d}1.4M_{}`$ in systems with $`P_{\mathrm{orb}}2.5\mathrm{days}`$ on the other hand may fill their Roche lobe when the donor star ascends the first giant branch, but the resulting mass transfer is dynamically unstable so that a common-envelope phase ensues instead of a LMXB phase. Systems with $`P_{\mathrm{orb}}800\mathrm{days}`$, finally, are too wide for the donor star to fill its Roche lobe prior to the AGB stage.

The lack of pre-LMXBs in the upper-left corner of each of the plots presented in Fig. 2 is related to our treatment of the coronal X-ray luminosity of the donor stars. Since stars less massive than $`1M_{}`$ are essentially unevolved and have appreciable convective envelopes, the adopted upper limit for the coronal X-ray luminosity is proportional to the square of the zero-age main sequence radii \[Eq. (3)\]. The decrease of the latter with decreasing mass then implies that binaries with lower-mass donor stars can satisfy the criterion $`L_{\mathrm{acc}}>10L_{\mathrm{cor}}`$ up to longer orbital periods. The exact position of the lower border of the wedge in the upper-left corner of the figures introduced by this criterion depends on the wind-velocity parameter $`\beta _\mathrm{w}`$ and on the assumptions made for the determination of the coronal X-ray luminosity $`L_{\mathrm{cor}}`$. Increasing the wind-velocity parameter (model MA2) yields lower mass-accretion rates and thus lower accretion luminosities. Correspondingly, the bottom of the wedge shifts to lower orbital periods. Adopting a constant upper limit for the coronal X-ray luminosity (models COR1 and COR2, not shown) instead of Eq. (3) renders the border of the wedge introduced by the $`L_{\mathrm{acc}}>10L_{\mathrm{cor}}`$ criterion independent of the mass of the donor star.

Stars more massive than $`1M_{}`$, on the other hand, may evolve away from the ZAMS within the imposed age limit of 10 Gyr. However, the deep convective envelopes and large radii developed on the giant branch yield high upper limits for the coronal X-ray luminosity so that only few wind-accreting neutron stars manage to top this limit with their accretion luminosity. Most of the pre-LMXBs with $`M_\mathrm{d}1M_{}`$ therefore also have main-sequence donor stars. The rather sharp vertical border of the wedge in the upper-left corner of the figures at $`M_\mathrm{d}1M_{}`$ is then related to the assumed lower limit $`M_{\mathrm{conv}}=0.01M_\mathrm{d}`$ for the mass of the convective envelope of coronally active stars. Main-sequence stars more massive than $`1M_{}`$ have smaller or no convective envelopes and are thus assumed to have negligible coronal X-ray luminosities. The pre-LMXB luminosity criterion therefore reduces to $`L_{\mathrm{acc}}>10^6L_{}`$, which is easily satisfied. If the minimum convective envelope mass necessary for X-ray activity is increased to, e.g., $`M_{\mathrm{conv}}=0.1M_\mathrm{d}`$ (model COR3, not shown), the sharp vertical border of the wedge shifts to $`M_\mathrm{d}0.8M_{}`$.

The group of systems with orbital periods longer than $`1000`$ days and donor star masses between 1.0 and 1.5 $`M_{}`$ in the low kick-velocity dispersion models shown in Fig. 2 corresponds to systems with initial orbital separations wide enough to avoid any kind of Roche-lobe overflow prior to the formation of the neutron star. These systems easily survive the supernova explosion of the neutron star’s progenitor when the average kick velocity is small compared to the relative orbital velocity of the component stars, but generally get disrupted for higher average kick velocities.

The narrow horizontal ridge of systems near $`P_{\mathrm{orb}}10`$ days in the low kick-velocity dispersion models ($`\sigma _{\mathrm{kick}}100`$ km/s) is associated with systems that spiralled-in during a common-envelope phase triggered by Roche-lobe overflow from the neutron star’s progenitor prior to its supernova explosion. The range of orbital periods available for this evolutionary channel is limited by the requirements that the pre-common-envelope orbit must be small enough to initiate Roche-lobe overflow prior to the supernova explosion and wide enough to provide enough orbital energy to expel the common envelope before the core of the Roche-lobe overflowing star merges with the spiralling-in companion. After the expulsion of the envelope, the change in the orbital configuration caused by the supernova explosion depends on the amount of mass lost from the system and on the kick velocity imparted to the neutron star at birth. For a detailed account of the effects of asymmetric supernova explosions on the orbital characteristics of binary stars, we refer to Kalogera (1996).

In the case of small supernova kick-velocity dispersions, the effect of the mass loss dominates over the effect of the kick. The small spread of orbital periods available for the common-envelope scenario is then preserved by the one-to-one relation between the post-supernova orbital period and the mass lost from the system (e.g. Verbunt 1993). Furthermore, since binaries in which more than half of the total mass is lost from the system become unbound after a symmetric supernova explosion, the probability of survival is larger for binaries with more massive donor stars.

In the case of higher supernova kick-velocity dispersions the effect of the kick dominates over the effect of the mass loss. This facilitates the survival of binaries with lower-mass donor stars and increases the spread of the post-supernova orbital periods. The horizontal ridge around $`P_{\mathrm{orb}}10`$ days therefore widens and shifts to lower donor star masses.

The common-envelope ejection efficiency parameter $`\alpha _{\mathrm{CE}}`$ only affects systems with orbital periods less than $`1000`$ days. A decrease of $`\alpha _{\mathrm{CE}}`$ (model CE1) yields smaller orbital separations at the end of the common-envelope phase and thus shifts the horizontal ridge near $`P_{\mathrm{orb}}10`$ days to lower orbital periods. This effectively widens the gap between the long- and short-period systems. Increasing the common-envelope efficiency $`\alpha _{\mathrm{CE}}`$ has the opposite effect.

We conclude this section by noting that the position of the high-density areas in Fig. 2 also depends on the adopted initial mass-ratio distribution. Mass ratio distributions favouring lower initial mass ratios (e.g. $`n(q)=1/q`$) shift the high-density areas to lower donor star masses, while mass ratio distributions favouring higher initial mass ratios (e.g. $`n(q)=q`$) shift the high-density areas to higher donor star masses.

## 6 X-ray luminosities

The normalised distribution functions of pre-LMXBs in the $`(\mathrm{log}P_{\mathrm{orb}},\mathrm{log}L_{\mathrm{acc}})`$\- and the $`(\mathrm{log}M_\mathrm{d},\mathrm{log}L_{\mathrm{acc}})`$-plane are displayed in Figs. 3 and 4, respectively. Since the qualitative differences between the models listed in Table 1 are small, we limit the figures to models KM50, A, and MA2.

The bulk of the systems shown in Fig. 3 have logarithmic accretion luminosities that are linearly correlated with the logarithm of the orbital period. The slope of -4/3 arises from Kepler’s third law and the factor $`1/a^2`$ in Eq. (2) for the mean mass-accretion rate $`\dot{M}_{\mathrm{acc}}`$. Since most of the systems following the correlation have main-sequence donor stars, the small spread around the correlation is the result of the limited range of donor star masses and radii and the adopted constant wind mass-loss rate of $`10^{13}M_{}`$ per year \[see Eq. (2)\]. In the case of low supernova kick-velocity dispersions, the two distinct higher-density regions along the $`L_{\mathrm{acc}}P_{\mathrm{orb}}^{4/3}`$ relation correspond to the two groups of systems identified in the low kick-velocity dispersion models shown in Fig. 2. The high accretion luminosities around $`10^{31}`$ erg/s stem from systems on the horizontal ridge around $`P_{\mathrm{orb}}10`$ days, while the lower accretion luminosities around $`10^{28}`$ erg/s stem from the long-period systems with $`P_{\mathrm{orb}}1000`$ days. For higher kick-velocity dispersions, only the first group of systems survives. The accretion luminosities of the pre-LMXBs on the $`L_{\mathrm{acc}}P_{\mathrm{orb}}^{4/3}`$ furthermore decrease with increasing values of the wind-velocity parameter $`\beta _\mathrm{w}`$.

The branch of systems located to the right of the $`L_{\mathrm{acc}}P_{\mathrm{orb}}^{4/3}`$ correlation corresponds to pre-LMXBs with giant branch donor stars. Most of these already appeared as pre-LMXBs when the donor star was still on the main-sequence, but the growth of the convective envelope and the associated increase in the coronal X-ray activity during the post-main-sequence evolution temporarily concealed the X-ray accretion luminosity emitted by the wind-accreting neutron star. The functional form of the mean mass-accretion rate $`\dot{M}_{\mathrm{acc}}`$ \[Eq. (2)\] implies that systems with longer orbital periods require more evolved donor stars, i.e. stars with larger radii and higher wind mass-loss rates, to satisfy the imposed $`L_{\mathrm{acc}}>10L_{\mathrm{cor}}`$ criterion. However, the larger radii in turn imply a larger upper limit for the coronal X-ray luminosity $`L_{\mathrm{cor}}`$ \[Eq. (3)\], so that at the time of reappearance the $`10L_{\mathrm{cor}}`$ threshold itself will be higher for systems with wider orbital separations. Hence, the accretion luminosity at which the systems reappear as pre-LMXBs increases with increasing orbital periods. Their reappearance also depends critically on whether or not the accretion luminosity increases faster than the coronal X-ray luminosity.

The position of the giant-branch stars in the $`(\mathrm{log}P_{\mathrm{orb}},\mathrm{log}L_{\mathrm{acc}})`$-plane clearly depends on the assumptions adopted for the determination of $`L_{\mathrm{acc}}`$ and $`L_{\mathrm{cor}}`$. Since higher wind velocities yield lower mass-accretion rates, an increase in the wind-velocity parameter $`\beta _\mathrm{w}`$ as in model MA2 requires even more evolved donor stars to satisfy the $`L_{\mathrm{acc}}>10L_{\mathrm{cor}}`$ criterion than in model A. The reappearance threshold of these systems consequently moves to higher accretion luminosities. In the case of a constant coronal X-ray luminosity (models COR1 and COR2, not shown), the accretion luminosity at which systems with giant branch donor stars reappear as pre-LMXBs is independent of the orbital period. The temporary concealment of the accretion luminosity probably does not occur for the longer period systems ($`P_{\mathrm{orb}}100`$ days) in which tidal forces are too weak to keep the donor star’s rotation synchronised with the orbital motion of the companion. The decrease of the donor star’s rotational angular velocity during its post-main-sequence evolution is here more likely to be accompanied by a decrease in the coronal X-ray activity so that the wind-accretion luminosity in these binaries has no problem topping the coronal X-ray luminosity at any time after they first appear as a pre-LMXB.

The linear correlation between $`\mathrm{log}L_{\mathrm{acc}}`$ and $`\mathrm{log}P_{\mathrm{orb}}`$ for the bulk of the systems presented in Fig. 3 essentially makes Fig. 4 a mirror image of Fig. 2. The low kick-velocity dispersion models hence again show two distinct groups of systems at high and low accretion luminosities, which are associated with the horizontal ridge and the long-period systems observed in Fig. 2. In our standard model (model A), even systems with donor stars less massive than $`1.0M_{}`$ are seen to give rise to accretion luminosities up to $`10^{32}`$ erg/s. An increase in the wind-velocity parameter $`\beta _\mathrm{w}`$ shifts all systems to lower accretion luminosities.

Finally, the one-dimensional distribution functions of the pre-LMXB accretion luminosities are presented in Fig. 5 for the same population synthesis models as in Fig. 2. We distinguish between systems that actually evolve into LMXBs within the imposed age limit of 10 Gyr and those that do not by means of light and dark grey shadings, respectively. The distribution of pre-LMXB accretion luminosities in our standard population synthesis model shows a primary peak around $`10^{31}`$ erg/s and a small secondary peak around $`10^{28}`$ erg/s. The latter becomes more important with decreasing kick-velocity dispersions and dominates when $`\sigma _{\mathrm{kick}}=25`$ km/s. This shift towards lower accretion luminosities is related to the increasing contribution of the long-period systems with $`P_{\mathrm{orb}}1000`$ days. The width of the peak at $`10^{31}`$ erg/s depends primarily on the common-envelope ejection efficiency parameter $`\alpha _{\mathrm{CE}}`$. The peak becomes narrower with decreasing values of $`\alpha _{\mathrm{CE}}`$ due to the shrinking range of post-supernova orbital periods around $`P_{\mathrm{orb}}10`$ days. Increasing the wind-velocity parameter $`\beta _\mathrm{w}`$ moves the distribution function to lower X-ray luminosities, while a decrease in $`\beta _\mathrm{w}`$ has the opposite effect. A similar tendency was found by Pfahl et al. (2002). It is also interesting to note that the peak luminosity of $`10^{31}`$ erg/s is comparable to the luminosity originating from non-magnetic white dwarfs accreting mass from a Roche-lobe filling companion (cataclysmic variables). Such systems may be distinguishable from pre-LMXBs based on the harder X-ray spectrum expected for wind-accreting neutron stars. In addition, white dwarfs can possibly be detected in the optical or UV wavelengths, provided that they are bright enough.

## 7 Relative and absolute number of pre-LMXBs in the Galaxy

The relative number of pre-LMXBs evolving into LMXBs within the imposed age limit of 10 Gyr are listed in the first column of Table 2. In general, about 20 to 30 per cent of all systems satisfying our pre-LMXB criteria evolve into a LMXB in a time span of less than 10 Gyr. This fraction decreases significantly for kick-velocity dispersions smaller than $`50`$ km/s (model KM50) due to the growing contribution of systems with orbital periods longer than $`1000`$ days. A decrease of the common-envelope efficiency parameter $`\alpha _{\mathrm{CE}}`$ (model CE1) on the other hand increases the relative number of pre-LMXBs evolving into LMXBs since it boosts the contribution of the short-period systems. Similarly, increasing the wind-velocity parameter $`\beta _\mathrm{w}`$ (model MA2) decreases the accretion luminosities of the longer period systems below the $`L_{\mathrm{acc}}>10L_{\mathrm{cor}}`$ threshold so that again the relative number of pre-LMXBs evolving into LMXBs increases. If Eq. (3) is replaced by a constant coronal X-ray luminosity of $`10^2L_{}`$ (model COR1), a large fraction of the low-mass systems with $`M_\mathrm{d}1M_{}`$ can no longer satisfy the $`L_{\mathrm{acc}}>10L_{\mathrm{cor}}`$ criterion. The disappearance of these systems causes an increase of the relative number of pre-LMXBs that will evolve into LMXBs.

Although a detailed population synthesis study of the formation and evolution of LMXBs is beyond the scope if this investigation, it is interesting to consider the relative number of LMXBs that are expected to evolve through a phase satisfying our pre-LMXB criteria. The relative numbers are listed in the second column of Table 2 and typically range 30 to 50 per cent. The LMXBs that did not evolve through a pre-LMXB phase all descend from intermediate mass X-ray binaries (IMXB) with initial donor star masses outside the mass range considered in this investigation (see also King & Ritter 1999; Podsiadlowski & Rappaport 2000; Kolb et al. 2000; Pfahl, Rappaport & Podsiadlowski 2003). A substantial fraction of LMXBs thus evolves through an X-ray bright wind-accretion phase prior to the onset of Roche-lobe overflow from the neutron star’s companion.

For conclusion, we estimate the order of magnitude of the birthrate and the absolute number of pre-LMXBs currently populating the Galaxy. To this end, we assume all stars to be in binaries and we normalise the star formation rate to produce a Galactic type II supernova rate of $`10^2`$ per year (Cappellaro, Evans & Turatto 1999). The resulting birthrates and absolute numbers are listed in Table 3. Based on these calculations, we expect the current population of pre-LMXBs in the Galactic disk to exist of about $`10^4`$$`10^5`$ systems, so that they may account for a non-negligible part of the X-ray output in large scale surveys performed by Chandra and XMM.

## 8 Discussion

We have investigated the population of candidate pre-low-mass X-ray binaries according to the detection mechanism proposed by Bleach (2002). The mechanism is based on the search for low-mass stars with X-ray luminosities that are too large to be attributed to the stellar coronal activity, and on the subsequent identification of neutron star binaries amongst these stars. The X-ray excess in these binaries can then possibly be ascribed to the accretion of mass by the neutron star from the stellar wind of its companion.

In our study, we used the BiSEPS binary population synthesis code described by Willems & Kolb (2002). We limited ourselves to neutron star binaries with main-sequence, Hertzsprung-gap, or giant-branch donor stars of mass $`M_\mathrm{d}2M_{}`$ and simulated the X-ray excess by the requirement that the X-ray luminosity generated by the wind-accreting neutron star be at least 10 times larger than the coronal X-ray luminosity of the mass-losing star. The X-ray luminosity associated with the accretion process was determined under the assumptions that no accretion-inhibiting effects such as strong centrifugal forces are present and that all gravitational potential energy of the infalling matter is converted into X-ray radiation. The maxiumum coronal X-ray luminosity for a star of a given mass and radius was assumed to be proportional to the square of the star’s radius.

It follows that, in the case of low supernova kick-velocity dispersions, the population of pre-LMXBs can be divided into two groups of systems: a long-period group with $`P_{\mathrm{orb}}1000`$ days and a short-period group with $`P_{\mathrm{orb}}10`$ days. The division yields a bimodal distribution for the accretion luminosities of the candidate pre-LMXBs. The long-period systems produce a peak in the distribution around $`10^{28}`$ erg/s, while the short-period systems produce a peak near $`10^{31}`$ erg/s. For kick-velocity disperions larger than $`100`$ km/s only the peak near $`10^{31}`$ erg/s remains. The luminosity at which the distribution peaks furthermore depends sensitively on the velocity of the donor star’s wind.

These results are derived under the assumption that the stellar and orbital parameters, and in particular the orbital period and the accretion luminosity, of pre-LMXBs remain unchanged after their formation. This assumption breaks down for binaries with $`M_\mathrm{d}1M_{}`$ and $`P_{\mathrm{orb}}1`$ day in which the evolution is governed by the loss of angular momentum via gravitational radiation and/or magnetic braking, and for binaries with $`M_\mathrm{d}1M_{}`$ in which the evolution is governed by the nuclear evolution of the donor star.

First, the decrease of the orbital period resulting from the angular momentum loss in binaries with $`M_\mathrm{d}1M_{}`$ and $`P_{\mathrm{orb}}1`$ day yields an increase of the mean mass-accretion rate and thus of the accretion luminosity generated by the neutron star. However, from Eqs. (1) and (2) it follows that an increase of the accretion luminosity by a factor of 10 requires a decrease of the orbital period by a factor of $`5.6`$. For a pre-LMXB with an orbital period of $`1`$ day, this comes down to a decrease of the period to $`4.3`$ hours. Since neutron star binaries with secondaries more massive than $`0.5M_{}`$ becomes semi-detached before they can reach such short periods, we can expect that the neglect of the orbital evolution does not affect our estimates for the accretion luminosities by more than an order of magnitude.

Second, the nuclear evolution of donor stars more massive than $`1M_{}`$ may yield a significant increase of the star’s radius and wind mass-loss rate within the age of the Galaxy. The effect of both these evolutionary aspects is to increase the mean mass-accretion rate of the neutron star. However, for main-sequence stars less massive than $`2M_{}`$, the increase of the stellar radius increases the mean mass-accretion rate by less than an order of magnitude \[see Eqs. (1) and (2)\]. The overall effect on our results can therefore again be expected to be small. The main uncertainty for main-sequence donor stars consequently stems from their poorly understood wind mass-loss rates (e.g. Wargelin & Drake 2001). Pre-LMXBs with donor stars that have evolved beyond the main-sequence, on the other hand, have better understood mass-loss rates, but contribute too little to the total population of newborn pre-LMXBs to have a significant effect on the distribution functions presented. In addition, systems with an orbital period of the order of a few days will evolve into LMXBs shortly after the secondary leaves the main sequence. A significant increase in the X-ray accretion luminosity is therefore only possible for wider systems which do not initiate Roche-lobe overflow until the secondary is well on its way towards the tip of the giant branch or for systems which do not initiate Roche-lobe overflow at all. For these systems, we may underestimate the maximum accretion luminosity by several orders of magnitude, but since the life time of a star on the giant branch is much shorter than its main sequence life time, the contribution of the systems to the pre-LMXB population during this phase is smaller than their contribution when the donor star was still on the main sequence. In addition, the wider orbital separations required for this scenario somewhat counteract the increase in the accretion luminosity associated with the evolution of the donor star.

Finally, we estimated the total number of pre-LMXBs currently populating the Galaxy to be of the order of $`10^4`$$`10^5`$. For comparison, the estimated number of RS CVn type binaries in the Galaxy is of the order of $`10^6`$ (Watson 1990) and the theoretically predicted number of cataclysmic variables in the Galaxy is of the order of $`10^610^7`$ (Kolb & Willems 2003, in preparation). Both these types of systems have X-ray luminosities comparable to the $`10^{31}`$ erg/s peak luminosity found in our models for candidate pre-LMXBs. In addition, the distribution of X-ray luminosities at $`10^{28}`$ erg/s may be expected to be dominated by coronally active stars. Our predicted number of pre-LMXBs satisfying the X-ray excess criterion is therefore small in comparison to the number of these other low-luminosity X-ray sources, but nevertheless high enough to possibly account for an interesting fraction of the low-luminosity discrete X-ray sources observed by Chandra and XMM. A systematic search for candidate pre-LMXBs based on the X-ray excess criterion can in principle be initiated by comparing exisiting X-ray (ROSAT, Einstein, Chandra and XMM) and optical (Hipparcos, Tycho, HST GSC) source catalogs. Simultaneous X-ray and optical observations by XMM provide an attractive alternative.

## Acknowledgements

We thank Jarrod Hurley, Onno Pols, and Chris Tout for sharing their SSE software package; and Philipp Podsiadlowski, Firoza Sutaria and Robin Barnard for useful discussions. The referee, Peter Wheatley, is thanked for his constructive remarks. This research was supported by the British Particle Physics and Astronomy Research Council (PPARC).


text:
# Comment on ”Negative heat capacities and first order phase transitions in nuclei” by L.G. Moretto et al

## Abstract

In a recent paper L.G. Moretto et al Moretto  claim that the negative heat capacities presented in our previously published paper PRL2000  are ”artifacts” coming from the use of periodic boundary conditions in the Lattice-Gas calculations. We stress in this comment that this claim is wrong: in ref. PRL2000  we did not use periodic boundary conditions and anyhow the boundary conditions are irrelevant for the statistical ensemble used in PRL2000 . The second claim of Moretto  is that, because of the Coulomb repulsion, systems ”with $`A>60`$ should present no anomalous negative heat capacities”. We show that this conclusion is contradicted by exact Lattice-gas simulations including Coulomb forces which present negative heat capacities even for $`A>200`$.

Let us start with the discussion about the boundary conditions. It is clearly stated in ref. PRL2000  that”$`L^3`$ \[the lattice size\] is large enough (typically greater than $`20^3`$ lattice sites) so that the boundary conditions do not affect the calculations with a constraining $`\lambda `$.” This means that the results of ref. PRL2000  do not depend upon the conditions used at the boundary, and in fact the calculations were made without periodic boundary conditions. At that time we checked the independence of the boundaries comparing $`N=8000`$ and $`N=27000`$ lattices. To see a sizeable effect of the boundary, we have decreased the size of the lattice. A fast calculation of $`50000`$ events for $`216`$ particles in a $`N=5832`$ lattice at an energy $`E=0.4ϵ`$ confined by a Lagrange multiplier $`\mathrm{log}\lambda =8`$ , gives a temperature $`T=(0.687\pm 0.004)ϵ`$, a heat capacity $`C=16.3\pm 0.2`$, and a kinetic energy fluctuation $`\sigma _k^2/T^2=1.64\pm 0.03`$. If periodic boundary conditions are imposed, the temperature becomes $`T=(0.682\pm 0.004)ϵ`$, the heat capacity $`C=16.9\pm 0.4`$ and the fluctuation $`\sigma _k^2/T^2=1.63\pm 0.07`$. The temperature decrease in the phase transition region is $`\mathrm{\Delta }T=(7.6\pm 0.4)10^3ϵ`$ ($`\mathrm{\Delta }T=(6.7\pm 0.6)10^3ϵ`$) with (without) periodic boundary conditions. This shows that the boundary conditions do not affect the thermodynamics. Indeed, in ref. PRL2000  we have analyzed an ensemble of particles in the vacuum for which only the average spatial extension of the system is defined (isobar ensemble). Then, provided that the lattice in which the system has still to be discretized for technical reasons is big enough, boundary conditions are irrelevant since the particles never explore the outer region.

In a previous paper PRL1999  we have considered a canonical Lattice Gas model in a box of constant volume (isochore ensemble), and there we have used periodic boundary conditions. In these calculations negative compressibilities are reported but the heat capacity is always positive. In fact, in the canonical ensemble the heat capacity is proportional to the energy variance and can never be negative Shro .

Concerning the isochore microcanonical ensemble, the published results Richert  show that the heat capacity at constant volume $`C_V`$ is always positive. Therefore, contrary to what is claimed in ref. Moretto , published results dealing with Lattice Gas calculations with typical nuclear size systems and periodic boundary conditions Richert ; IMFM  do not show any negative heat capacity.

We may incidentally note that recently a debate has started in the statistical physics community after the observation of a negative heat capacity branch for very large systems in the the constant-magnetization microcanonical Ising model Pliemling  (equivalent to the microcanonical isochore Lattice Gas model). In these cases however, the densities involved are so low that once again boundary conditions are irrelevant. Moreover, the sizes discussed are orders of magnitude larger than the nuclear ones.

Finally we would like to stress that we have recently demonstrated that negative heat capacities in finite systems are the origin of first order phase transitions with a finite latent heat: to present an energy discontinuity at the thermodynamic limit, finite systems must present a negative microcanonical heat capacity if the number of particles is large enough zeroes .

Let us now comment about the effect of Coulomb. The conclusion of ref.Moretto  is based on a simple model of a unique cluster in equilibrium with a gas. However, this configuration is not the most probable one both in statistical models smm ; gross  and in nuclear physics experiments at high excitation energy data . In particular, the negative heat capacity region is characterized by a large ($`m3`$) multiplicity of fragments palluto . Therefore any conclusion drawn considering only the specific channel of multiplicity 1 is not valid. If we would follow the simple model of ref. Moretto  and consider partitions with $`m=35`$ fragments of similar mass data , this would lead to a figure analogous to figure 5 of ref. Moretto  but with an abscissa approximately multiplied by $`m`$ (or more, since all the light particles have to be included), because the inter-fragment coulomb interaction can be neglected if the volume is large. One would then find $`A200`$ as the limit of negative heat capacities.

Of course fragments do not have the same size and to get a quantitative result one needs to compute the relative weight of all the possible channels; this is indeed what is done in nuclear statistical models gross ; smm ; raduta  which show that the intermediate state between the compound nucleus and the vaporized system is multifragmentation: the presence of many drops in equilibrium.

In the case of the Lattice Gas model the different partitions and their relative weight can be calculated without any approximation. In the negative heat capacity energy regime there are in average three fragments of size greater than 4 PRL2000 . The sudden opening of this multifragment channelgross  causes a convex intruder in the entropy that is responsible of the negative heat capacity. The effect of the Coulomb interaction in the Lattice Gas model can be appreciated from oFigure 1 which shows a caloric curve at constant $`\lambda =\beta p`$ for a system of $`A=207`$ particles and charge $`Z=82`$ vero . A clear backbending is visible. In the charged Lattice-gas model, the energy interval corresponding to the backbending decreases with increasing charge vero , showing that the idea that Coulomb tends to suppress the negative heat capacity is qualitatively correct. However, when all the available channels are correctly weighted, the quantitative effect is very different then the one reported in ref. Moretto . Only a slight effect of the Coulomb interaction has also been reported in the isochore framework carmona .

Obviously the result of Figure 1 is model dependent. Different macroscopic models show a higher raduta  or lower smm  sensitivity to the Coulomb depending on the detailed implementation of the Coulomb and nuclear interactions. On the other hand molecular dynamics (which can be solved exactly as the Lattice Gas model) leads to caloric curves which are almost independent of Coulomb claudio . Whether this can be interpreted as an evidence of metastable long time tails, the relaxation time of long range interactions being excessively long, is a subject of debate tsallis . In any case these examples and the results of Figure 1 show that a definitive understanding of the effect of a non saturating long range interaction on a first order phase transition cannot certainly be achieved through the oversimplified model of ref. Moretto .


text:
# Numerical Analysis of the Stub Transistor

## I Introduction

Many devices are being proposed to substitute the MOSFET, as microelectronics reaches the nanoscale. Among them are the stubbed waveguide and the stub transistor Sols1988 ; Sols1989 ; DJ1994 ; DM2000 ; DRVRPM2000 ; WVP2002 . The stubbed waveguide consists of a nanowire with a set of periodically spaced stubs. The stub transistor is a nanowire with an electric potential applied to the stub, which acts as a gate voltage, see Fig. 1. Such devices can be fabricated with dimensions below $`100nm`$, which, at low temperatures, is smaller than the length scale over which the electron preserves its phase coherency. This length is known as the electron phase-coherence length, $`L_\varphi `$, and the transport through the device is ballistic. These kind of structures are classified as mesoscopic systems and cannot be described by the usual semi-classical transport theory, the wave nature of the electron needs to be taken explicitly into account.

The conductance, $`G`$, of the stub transistor displays an oscillatory behavior, as the electric potential applied to the stub, $`V_G`$ is varied. Therefore, the channel can be opened and closed. A minimum of $`G`$ is a reflection resonance or antiresonance and a maximum of $`G`$ is a transmission resonance or just resonance DJ1994 . The transmission zeroes occur as function of the length of the stub, $`L_{stub}`$, given by $`kL_{stub}=n\pi `$, and the resonances are given by $`kL_{stub}=(n+1/2)\pi `$. The use of multiple periodically placed stubs transforms the transmission zeroes in blocked transmission bands DJ1994 . Such devices have been fabricated on GaAs/AlGaAs by MBE DRVRPM2000 .

To understand the behavior of the device, many theoretical techniques have been used, such as, the recursive Green function technique Sols1989  to calculate the scattering matrix, the transfer matrix technique DJ1994 ; WVP2002 . Santos Edval2002  has developed a recursive Green function program to simulate ballistic quantum devices.

The reason for the great interest in the stub transistor is the potential for ultrafast signal processing without giving prohibitive energy dissipation. The traditional FET transistor works using a “brute” force switching voltage, because the electric field has to deplete electrons under its gate to change the conductance. The quantum interference devices, QIDs, on the other hand, work by changing the electron interference pattern. This quantum phenomenon requires a lower switching voltage, therefore a lower energy dissipation. The reduced size of the QIDs decreases the time that the electron takes to travel through it, as a consequence, its switching frequency can be increased up to the terahertz frequency region MKS2000 .

Strong analogies exist between ballistic quantum devices and similar devices in optics and microwave engineering. One, however, has to take care when try to implement optical or microwave devices using matter electronic waves. These analogies are not perfect because electrons interact strongly and obey Fermi statistics, while photons are bosons.

In this paper, the quantum stub transistor is examined by applying the recursive Green function technique. The paper is divided in five sections, this introduction is the first, next the quantum stub transistor is described. Third an overview of the simulation method is presented. Next, the results and analysis are presented. Finally, the conclusions.

## II The quantum stub transistor

A schematic view of the quantum stub transistor is shown in Fig. Numerical Analysis of the Stub Transistor, where $`V_G`$ is the gate voltage and $`V_D`$ is the drain voltage. The drain voltage is related to the incident electron energy by $`E=eV_D`$. If the stub length, $`L_{stub}`$, is much larger than the inelastic mean free path $`L_{ine}`$, then this structure represents just two intersecting wires and the transconductance, $`g_m=I_D/V_G`$ is zero, where $`I_D`$ is the drain current. For $`L_{stub}L_{ine}`$, however, the quantum nature of the electronic transport makes the transconductance $`g_m`$ different from zero and the situation changes completely. What happens is similar to the conventional metal waveguide, where one has to move a piston inside a cavity to change the interference patterns of the electromagnetic field. In the case of the quantum stub transistor, the gate voltage acts as a piston by varying the length of the stub. This changes the penetration of the electron inside the stub, and therefore the interference pattern. To represent this change of the stub length with the gate voltage, $`L_{eff}`$ is defined. The stub length $`L_{stub}`$ for which quantum effects become important vary widely depending on temperature and material used. In ultra-pure semiconductor material and at liquid Helium ($`4.2K`$) temperatures and below, coherent effects are observable over distances of $`1\mu m`$ Tpalm1993 .

The confinement of the propagating electrons from source to drain in the transverse direction allows only discrete energy levels or modes. The number of transverse modes, $`M`$, for a quantum stub transistor with width $`W`$ is given by the ratio between the width and the Fermi wavelength, $`\lambda _F`$, as shown in Equation 1.

$$M=Int\left(\frac{W}{\frac{\lambda _F}{2}}\right)$$
(1)
where $`Int(x)`$ represents the integer that is just smaller than $`x`$.

If more than one mode propagates through the device, each mode will have a different penetration length in the stub, as a consequence the device will display an erratic interference pattern. Hence, for a well behaved device the incident energy, $`E`$, must be in the range $`E_1E<E_2`$, where $`E_1`$ and $`E_2`$ are the energy of the first and second mode respectively. The energy of the $`n`$-th mode can be estimated with the infinite potential box approximation.

$$E_n=\frac{\mathrm{}^2}{2m^{}}\frac{n^2\pi ^2}{W^2}$$
(2)
where $`m^{}`$ is the electron effective mass.

The electron may follow many different paths, as it propagates through the device. The maximum length difference occurs between a path going straight from source to drain and another going through the stub and then to the drain. This difference is approximately $`2L_{stub}`$. So the number of minima expected in the interference pattern is $`m+1`$, where $`m`$ is the maximum integer that satisfies the inequality below

$$L_{stub}(m+\frac{1}{4})\lambda _F.$$
(3)
For example, for $`L_{stub}=30nm`$ and $`\lambda _F=20nm`$, one gets $`m+1=2`$.

## III The simulation method

To simplify the analysis, it is assumed that the potential inside the device is equal to zero, and the borders of the device are defined by a high potential barrier. Even with these simplified boundary conditions, the two dimensional nature of the devices makes the problem of finding the total transmission probability, $`T`$, of the electron through the device, analytically intractable. To obtain the conductance of the device using the Laudauer formula ($`G=\frac{2e^2}{h}T`$), it is necessary to use a numerical method to solve Schrödinger equation ($`H\mathrm{\Psi }=E\psi `$). Hence, the geometry has to be described with a model Hamiltonian, this is done, by using the nearest-neighbor tight-binding model. To implement such model, the structure showed in Fig. Numerical Analysis of the Stub Transistor is filled with a sufficient dense square lattice with periodicity $`a`$. The plane wave motion can be emulated by a tight-binding Hamiltonian, when the lattice parameter, $`a`$, is much smaller than all other length scales of the problem. Proceeding as indicated, one may reasonably expect to correctly reproduce the continuous dynamics.

To better illustrate the method, an example of a square lattice is shown in Fig. Numerical Analysis of the Stub Transistor. Each slice is stored as a column of the device matrix, where the matrix entry is related to the potential energy confining the electron. Each column (slice) of the device matrix is used to construct the Hamiltonian or slice matrix for numerical calculation by using the tight-binding Hamiltonian Model. The $`n`$th slice matrix is shown in Equation 4.

$$H_n=\left(\begin{array}{ccccc}V_14t& t& 0& \mathrm{}& 0\\ t& V_24t& t& \mathrm{}& 0\\ 0& t& V_34t& \mathrm{}& 0\\ .& .& .& & .\\ .& .& .& & .\\ .& .& .& & .\\ 0& 0& 0& \mathrm{}& V_M4t\end{array}\right)$$
(4)
where $`V_i`$ in the diagonal is the local potential energy profile at a given slice, and $`t`$ is a coupling term, given by Equation 5.

$$t=E_F\left(\frac{1}{2\pi }\right)^2\left(\frac{\lambda _F}{a}\right)^2$$
(5)
where $`E_F`$ is the Fermi energy level.

The Schrödinger equation is now solved by using the Green function method, i.e., the equation to the solved has the form:

$$(E\underset{¯}{I}\underset{¯}{H})\underset{¯}{G}=\underset{¯}{I}$$
(6)
where $`\underset{¯}{I}`$ is the identity matrix, and $`\underset{¯}{G}`$ is the Green function matrix. For simulations with practical interest, the direct inversion of the matrix $`(E\underset{¯}{I}\underset{¯}{H})`$ requires huge amounts of memory and processing power. This is circumvented by employing an iterative method for the computation of the Green function. The details of the program can be found in Reference Edval2002 .

The data input process is time consuming, as the device has to be described in matrix format. To facilitate data input, a front-end that converts the device geometry from a standard bit-mapped format into the matrix format that the computation algorithm understands has been developed GS2002 . The relation between the width $`W`$ of the device and the number of pixels $`P_W`$ used to represent this dimension is given by

$$W=(P_W+1)a.$$
(7)
For $`a=0.5nm`$ and $`W=10nm`$, $`P_W=19`$.

## IV Results and analysis

The program uses normalized units, $`E_F`$ is the unit of energy and $`2e^2/h`$ is the unit of conductance. In the simulations, it is assumed that the gate voltage, $`V_G`$, drops linearly along the stub111The program does not have a Poisson solver for the calculation of the gate potential distribution, this feature will be added in the near future.. The effective mass is $`m^{}=0.07m_0`$ (GaAs), the wire width is $`W=10nm`$ and the lattice parameter is $`a=0.5nm`$.

The normalized coupling term is calculated with Equation 5. Considering the single mode device, $`t=40.52`$. Using an incident energy equal to the first allowed mode $`E_1=53meV`$, the plot of the conductance as a function of the normalized incident energy is presented in Figure Numerical Analysis of the Stub Transistor. The result confirms that the device has only one mode. The gate voltage is set equal to zero, $`L=73.3nm`$ and $`L_{stub}=11nm`$.

Considering now a device with four open transmission channels, the simulation yields a conductance with interference patterns, at each conductance step DRVRPM2000 . The plot of the normalized conductance as a function of the normalized incident energy is presented in Figure Numerical Analysis of the Stub Transistor. For this simulation, $`W=10nm`$, $`P_W=19pixels`$, $`L_{stub}=50nm`$, $`M=4`$ (four modes), and $`t=2.533`$.

Next, the incident energy (drain voltage) is fixed at a position to get just one open channel, and the gate voltage is varied. The normalized conductance as a function of the gate voltage are presented in Figure Numerical Analysis of the Stub Transistor and Numerical Analysis of the Stub Transistor. The results show that the conductance varies between one and zero, as the gate voltage varies. In Figure Numerical Analysis of the Stub Transistor, there are two conductance dips or antiresonances (zero conductance). In Figure Numerical Analysis of the Stub Transistor, there are three antiresonances. This can be used to test the correctness of Equation 3, two values $`L_{stub}`$ are used. For this calculation, the incident energy is $`53meV`$ and $`\lambda _F=20nm`$. Hence, for $`L_{stub}=29nm`$, one gets $`m=1`$, and the number of minima is $`m+1=2`$. For $`L_{stub}=50nm`$, one gets $`m=2`$, and the number of minima is $`m+1=3`$. The calculated number of minima are in accordance with the simulation.

Another configuration is the gate-all-around stub transistor. This transistor has a double stub, one opposite to the other. The calculated conductance as a function of the gate voltage is presented in Figure Numerical Analysis of the Stub Transistor. The presence of the gate-all-around stub narrows the antiresonance.

These are interesting results, because one can control the channel transmission, by applying a control voltage, which is useful in digital logic, for instance. The main practical difficulty is that the width of the antiresonance is very narrow, and one would have to set the voltage with great precision. However, It is known that in a stubbed waveguide the presence of a series of stubs transforms the antiresonance in blocked bands DJ1994 . The idea is to increase the number of stubs to get a wider blocked transmission band. The simulation is now carried out varying the number of stubs along the device. The double stub transistor is shown in Figure 8. The result, presented in Figure 9, shows that the antiresonance is replaced by double antiresonances. By adding three stubs, one gets a triplet. Therefore, increasing the number of stubs may be a practical solution for the realization of stub transistors, as the blocked transmittance gets less sensitive to a particular gate voltage level.

## Conclusions

Besides the potential for high-speed operation, the appeal of the quantum stub transistor is its very small dimensions and the possibility of operation at very low power levels, which could allow extremely high-level of integration. But, as pointed by Laundauer, like the others (QIDs), this kind of device requires very specific values of the control signal and of the devices parameters for a conductance maximum or minimum. This leaves little room for device fabrication errors. However, the results in this paper show that the use of multiple stubs reduces this problem by widening the antiresonance. This may be a reasonable approach for the practical realization of stub transistors.

## Acknowledgments

The authors would like to acknowledge the support of CNPq under the “Instituto do Milênio” Initiative.

Figure captions

Fig. 1. Schematic view of the quantum stub transistor, where $`V_G`$ is the voltage applied to the stub, $`V_D`$ is the voltage applied to the nanowire, $`W`$ is the width of the nanowire, $`L_{stub}`$ is the stub length, and $`L_{eff}`$ is the stub effective length.

Fig. 2. Quantum device divided in slices for the simulation of quantum transport. The device is connected to two reservoirs.

Fig. 3. The quantum stub transistor with only one open transmission channel. $`W=10nm`$, $`a=0.5nm`$, $`P_W=19pixels`$, $`L_{stub}=11nm`$, $`M=1`$ (one mode), and $`t=40.52`$.

Fig. 4. The quantum stub transistor with four modes. $`W=10nm`$, $`P_W=19pixels`$, $`L_{stub}=50nm`$, $`M=4`$ (four modes), and $`t=2.533`$.

Fig. 5. Normalized conductance as a function of the gate voltage $`V_G`$, for $`L_{stub}=29nm`$.

Fig. 6. Normalized conductance as a function of the gate voltage $`V_G`$, for $`L_{stub}=50nm`$.

Fig. 7. Conductance of the gate-all-around stub transistor as a function of the gate voltage $`V_G`$, $`L_{stub}=50nm`$ at each side.

Fig. 8. Geometry of the double stub stub transistor, the stub separation is $`10nm`$, $`L_{stub}=50nm`$.

Fig. 9. Conductance of the double stub stub transistor as a function of the gate voltage $`V_G`$, the stub separation is $`10nm`$, $`L_{stub}=50nm`$ at each side.