text:
# BACKGROUND DISCRIMINATION CAPABILITIES OF A HEAT AND IONIZATION GERMANIUM CRYOGENIC DETECTOR EDELWEISS COLLABORATION

## 1 Introduction

The search for dark matter is one of the major challenges of contemporary physics. Despite extensive scans for MAssive Compact Halo Objects (MACHOs) initiated by several groups , it appears at present that the local density of dark matter can hardly be composed of baryonic matter essentially. Indeed, recent results from these experiments indicate that a significant fraction of the MACHOs observed in the direction of the Small Magellanic Cloud (SMC) are probably due to deflectors internal to the SMC, further reducing the possible amount of dark matter in the form of MACHOs in the halo itself .

The coincidence of the electroweak interaction scale (SUSY theory) with that required for Weakly Interacting Massive Particles (WIMPs) to contribute significantly to the solution of the dark matter problem is a further motivation for attacking the WIMPs direct detection challenge. In addition these hypothetical particles represent an attractive solution to the problem of galaxy formation. Experimental efforts are underway to detect these particles, either indirectly, by searching for products of their annihilation in the core of the Sun or of the Earth, or directly, by detecting the interactions of the WIMPs themselves in ordinary matter, as first suggested in the mid 1980s  .

The sensitivity of present experiments appears to be limited by the radioactive background rate of the detectors  and by the systematics of the rejection scheme using statistical identification methods based on pulse shape discrimination techniques (PSD) in NaI crystals  . The motivation for the development of cryogenic detectors is the possibility, when the heat measurement is coupled with a measurement of ionization  (or of scintillation ), to discriminate much more reliably the main source of radioactive background, producing electron recoils, from the nuclear recoils expected from WIMP interactions. Indeed, WIMPs scatter off nuclei, which ionize (or produce scintillation) less efficiently than electron recoils for a given energy deposit in the absorber.

In the following, we present the discrimination capabilities of this new generation of detectors in the EDELWEISS (Expérience pour DEtecter Les Wimps En SIte Souterrain) underground experiment and we show that a 70 g high-purity Ge crystal can already reach sensitivity levels typical of the best existing experiments despite its small mass and a relatively high radioactivity level in the detector environment. The present limitations of these detectors are then explored and the first attempts to counteract them are discussed.

## 2 The EDELWEISS experiment

The EDELWEISS cryogenic experiment operates in the Laboratoire Souterrain de Modane (LSM), an underground laboratory off the Fréjus highway tunnel beneath the French-Italian Alps. The 4600 mwe of rock reduce the muon flux to one part in two million of its surface value, that is to about 4.5 muons.m-2.day-1. The fast ambient neutron flux is measured to be $`\mathrm{4.\; 10}^6\mathrm{cm}^2.\mathrm{s}^1`$ .

The cryogenic setup and shielding are presented in . In particular, the dilution refrigerator has a cooling power of 100 $`\mu `$W at 100 mK, a useful volume for bolometer installation of about one liter and can reach temperatures as low as 10 mK. The optimized readout electronics are described in . The results presented here concern a prototype germanium detector of only moderate radioactive cleanliness, without near Roman lead shielding, radon removal or ambient neutron moderator.

### 2.1 The cryogenic detector

The absorber of the detector consists in a 70 g germanium monocrystal of high purity ($`n_Dn_A5.10^9\mathrm{cm}^3`$). The crystal shape is a disc 8 mm thick and 48 mm in diameter with bevelled edges.

Both charge and heat signals resulting from a particle interaction are measured by this detector. Charge collection is achieved by an electric field applied between two electrodes realized by boron implantation on one face of the monocrystal and by phosphorus implantation on the other. The applied bias voltages are of the order of a few volts. The p-i-n structure has been chosen to decrease as much as possible the injection of carriers in the detector volume and to allow large bias voltages. The doses of implanted ions are two orders of magnitude above the insulator-metal transition. The implantation energies have been chosen in order to reduce the thickness of the non-metallic implanted region and to ensure that the metallic layer reaches the detector surface. Recrystallization in the implanted region has been achieved by a fast thermal annealing and its quality checked by ellipsometry and by SIMS.111Second Ionisation Mass Spectroscopy. These electrode characterizations have been done with a previous detector by J.-P. Ponpon of the PHASE Laboratory (Strasbourg, France). In addition, sheet resistances have been measured on this last detector to verify the metallic behavior at very low temperatures . Several improvements of the detector design result from a study of the trap ionization and neutralization mechanisms in the bulk of the detector. The very low density of ionized traps we reached allows a noticeable improvement of the ionization channel time stability and energy resolution \[20-24\]. In order to minimize edge effects and to maximize charge collection, the thickness of the monocrystal has been reduced to 4 mm along its contour to increase the field strength, thereby reducing trapping and recombination.

The heat measurement is based on the variation, according to the Mott-Efrös-Schklovskii law, of the resistivity of a Neutron Transmutation Doped Germanium (NTD). This thermal sensor is glued to one face of the detector. The size of this NTD is 2x1x0.8 mm3. The heat sink consists of three copper wires 50 $`\mu `$m in diameter and about 2 cm in length. The operating temperature of the NTD sensor is about 20 mK for a mixing chamber temperature around 10 mK.

The whole detector lies on sapphire balls to ensure a good thermal decoupling from the detector holder made of low radioactivity copper and brass. A detailed description of the detector manufacturing and performances is given in \[19-24\].

### 2.2 Off-line data analysis

The small size of the heat and ionization signals (typically a temperature increase of $`10^6`$ K, and of the order of a thousand electron-hole pairs) exposes them to various electronic noises, be they fundamental or instrumental (Figure 1). Therefore, analysis of the data has been performed off-line using various methods including the optimal filtering technique in Fourier space  . This involves a least-square fit in frequency space between the actual event and a model event (obtained by averaging away the noise on several real events) while inversely weighting each frequency by the standard deviation of the noise at that frequency.

Neglecting pileup and assuming pulse shapes independent of pulse amplitude, each event $`s`$ can be represented as the sum of a scaled and shifted model $`m`$ and of a random noise $`n`$, $`s(t)=Am(tt_s)+n(t)`$. Further assuming the noise to be gaussian in frequency space, we can construct a $`\chi ^2`$ as the weighted difference between signal and scaled and shifted model :

$$\chi ^2=\underset{\omega }{}\frac{|N(\omega )|^2}{\sigma _\omega ^2}=\underset{\omega }{}\frac{|S(\omega )Ae^{i\omega t_s}M(\omega )|^2}{\sigma _\omega ^2},$$

where $`\sigma _\omega `$ is the standard deviation of the noise at frequency $`\omega /2\pi `$ and capital letters represent the Fourier transform of their lowercase counterparts. Minimization with respect to $`A`$ and $`t_s`$ yields an estimate of these two parameters. Parseval’s identity allows us to simply maximize over time shifts the scalar product $`A=_tf(tt_s)s(t)`$ where the function $`f`$ can be calculated from its Fourier transform once for all events

$$F(\omega )=\frac{e^{i\omega t_s}M(\omega )}{\sigma _\omega ^2}/\underset{\omega }{}\frac{|M(\omega )|^2}{\sigma _\omega ^2}.$$

Calculating the $`\chi ^2`$ to judge the quality of the fit is however more tractable in frequency space, and this remains computationally feasible thanks to Fast Fourier Transforms.

With respect to time fit techniques, this method selectively filters out noisy frequencies and thus is advantageous in the many cases of non-white noise spectra (Figure 1). Like straightforward time fit techniques however, it is ill-suited to data presenting pileups - though given current masses of bolometers this is not yet a problem when looking for WIMPs.

## 3 In situ detector calibrations

### 3.1  Energy normalization and recoil energy determination

The heat signal results from the sum of the heat deposited by the incident particle and of the heat generated by the charge carriers during their drift (the so-called Luke-Neganov effect ). This last term is given by the product of the bias voltage times the collected charge. Thus as bias voltage increases, the fraction of the heat signal directly correlated to the ionization signal increases, hampering separation of the electron recoils from the nuclear recoils. In practice, although bias voltages up to -12 V have been used, most of the data have been accumulated under bias voltages of -2 V and -6 V. There is in fact an additional electric field to that applied; it is caused by the 0.7 eV gap in Ge which creates a -0.7 V reverse field in the p-i-n diode. However, only the applied field ultimately counts for the Luke-Neganov effect : the energy equal to the additional field times the collected charge is indeed released when the charges are collected, but merely compensates the energy used to create the pairs in the first place .

Moreover it is standard procedure to normalize heat (phonon) energies (keV heat) to equivalent electron energies (keV e.e.). Volume electron recoils (with almost complete charge collection) are used to calibrate the charge and phonon channels. For such a volume electron recoil of a given energy (e.g. 122 keV for a 57Co calibration), both charge and phonon amplitudes ($`E_{ch}`$ and $`E_{ph}`$ respectively) are normalized to the deposited energy. In order to reconstruct this energy deposited in the detector, or recoil energy $`E_{rec}`$, from these quantities, the following formula  is applied:

$$E_{rec}=E_{ph}(1+\frac{eV}{ϵ_\gamma })E_{ch}\frac{eV}{ϵ_\gamma },$$

where $`V`$ is the bias voltage for charge collection (in Volts), $`e`$ the elementary charge and $`ϵ_\gamma `$ the average energy per electron-hole pair (in eV) for an electron recoil. The recoil energy resolution becomes thus bias-dependent, as shown in Figure 6. In the absence of Luke-Neganov effect ($`V=0`$) $`E_{rec}=E_{ph}`$; for a photon interaction where all the charge is collected (volume electron recoil) $`E_{rec}=E_{ph}=E_{ch}`$.

It is shown in  how the $`ϵ_\gamma `$ parameter can be measured using, in X- and $`\gamma `$-ray line calibrations, the small fraction of events which exhibits an incomplete charge collection. The measured value is close to the standard value $`ϵ_\gamma `$ 3 eV for the 77 K germanium diode detectors .

### 3.2 Electron recoil thresholds and resolutions

Electron recoil energy resolutions and thresholds have been studied using a 57Co source which provides essentially photoelectric interactions from its main peak at 122 keV. The ionization signal being roughly 1000 times faster than the heat signal, the former has usually been used as trigger for the acquisition with an effective threshold at 4 keV equivalent electron (e.e.) on the charge channel. Threshold efficiency has been studied by comparing a 60Co calibration with a Monte Carlo simulation. It is found to rise rapidly to 100% by 6 keV , see Figure 2.

The heat and ionization channels exhibit energy FWHM resolutions of $``$ 1 and $``$ 1.2 keV e.e. at 122 keV. According to the Erec expression of section 3.1, the corresponding resolutions on the recoil energy are $``$ 1.8 keV at -2 V and $``$ 3.8 keV at -6 V.

### 3.3 Signal and background calibrations

The detector has been calibrated with both neutron and gamma sources, in order to check its power of separation between nuclear recoil (”signal”) and electron recoil (background) events. High energy photons provide Compton interactions throughout the crystal rather than biased towards its surface.

A 252Cf neutron source, of about 2.5 $`\mu `$Ci activity, has been used for this calibration, at two different values of the reverse bias voltage (-2 V and -6 V). The branching ratio for neutron emission in the decay of the 252Cf isotope is about 12% and the average kinetic energy of the emitted neutrons is around 2 MeV. The 252Cf isotope is also a gamma source with a branching ratio of about 60%, and about 80% of these gammas have an energy less than 1 MeV. So, in principle this source alone is sufficient to obtain simultaneously both neutron and gamma calibrations. In order to improve the statistics without too long an exposure of the detector to the neutron source and the associated risk of neutron activation, we use also a 60Co source to mimick the background noise of electron recoils. Results are presented in Figure 3 for both polarization voltages, showing a remarkably good separation between gamma and neutron recoil events down to low energy. While most electron recoil events appear well behaved, and neatly line up along the main diagonal in the heat-ionization plane, approximately 5% suffer from incomplete charge collection.

The recoil energy $`E_{rec}`$ is derived according to section 3.1. Merging the data obtained at the two bias voltages yields the quenching factor $`(E_{ch}/E_{rec})`$ versus $`E_{rec}`$ in Figure 4. The curve $`4\mathrm{k}\mathrm{e}\mathrm{V}/E_{rec}`$ , which represents the threshold (at 4 keV e.e.), is also given. A ”neutron line” (the quenching factor for nuclear recoils) can be parameterized using the measured spreads and mean values of the charge over recoil ratio, calculated in discrete recoil energy intervals above 20 keV (under 20 keV a threshold bias is clearly visible). We obtain $`E_{ch}/E_{rec}=(0.16\pm 0.07)E_{rec}^{0.18\pm 0.1}`$, a result in agreement (within error limits) with those of the Heidelberg-Moscow collaboration , parameterized in  by $`E_{ch}=0.14E_{rec}^{1.19}`$, and in turn in agreement with the Lindhard theory  and most of the previous experiments (see  and references therein). The comparison with the results of the CDMS experiment , displayed in the Figure 4 inset, shows a very good agreement.

A neutron zone of known acceptance can be constructed by adjusting a power law of the recoil energy to the points at +$`n\sigma `$,-$`m\sigma `$ above and below the mean charge over recoil ratio. For a given background, the choice of the $`n`$ and $`m`$ values is a matter of optimizing the signal-to-noise ratio. On the Figure 4, the 95% acceptance neutron zone ($`\pm 2\sigma `$) is displayed as an example. This neutron calibration induces a slight bias in the determination of the quenching factor: the non-negligible number of multiple diffusions in the crystal leads to a merging of the quenching factor values at different energies. The width of the nuclear recoil band is overestimated accordingly for WIMPs for which it should be defined by single interactions only.

We now turn our attention to the separation capabilities of the detector based upon the two types of calibrations just discussed. In Figure 5 we apply the neutron zone method to the data of the 60Co calibration taken at -2 V. In the 15-100 keV recoil energy interval, 3 events (out of 3229) are found within the 95%-acceptance nuclear recoil band: (99.90 $`\pm `$ 0.05) % of the photon background is rejected. At -6V we find a rejection of (99.6$`\pm `$ 0.2) %. These results are very promising; however it must be kept in mind that they have been obtained from calibrations and that operational values may differ as we will see later.

## 4 A realistic case of data taking

Following the benchmarks just discussed, 0.65 kg.day and 1.17 kg.day of data were taken at the bias voltages of -2 V and -6 V respectively. Both spectra show a rate of roughly $`3540\mathrm{event}.\mathrm{kg}^1.\mathrm{day}^1.\mathrm{keV}_{\mathrm{recoil}}^1`$ in the 15-45 keV recoil energy interval, below the 46.5 keV peak due to 210Pb contamination.The recoil energy versus ionization plots of the data (Figure 6) show a large number of off-axis events. This population, interpreted as incomplete charge collection of electron surface events, is discussed in  together with another population which appears at recoil energies of the order of 80-100 keV and is attributed to surface nuclear recoils.

Given the large number of events with incomplete charge collection, we only retain the (0$`\sigma `$, -2$`\sigma `$) neutron zones as discussed in the previous section (Figure 7). This improves the signal-to-noise ratio and halves the acceptance to 0.475. The samples recorded respectively at -2 V and -6 V bias contain 6 events (respectively 7 events) in the 15-45 keV recoil range. We assume they are nuclear recoils, and use them to derive a conservative upper limit on a signal. For the respective exposures and acceptance, they yield upper limits at the 90% confidence level of 1.1 $`\mathrm{event}.\mathrm{kg}^1.\mathrm{day}^1.\mathrm{keV}_{\mathrm{recoil}}^1`$ at -2 V, and 0.7 $`\mathrm{event}.\mathrm{kg}^1.\mathrm{day}^1.\mathrm{keV}_{\mathrm{recoil}}^1`$ at -6 V. It will appear below that a large fraction of the events observed in the nuclear recoil band can be attributed to neutron interactions inside the detector. One can thus only set a lower limit of 98% for the rejection factor of the electron recoil background.

When we compare the mean ionization over recoil energy ratio of the off-axis events (electron surface events) at both bias voltages ($`0.46\pm 0.01`$ at -2 V and $`0.539\pm 0.006`$ at -6 V), we find that the higher bias voltage brings the off-axis events significantly closer to the well-collected events.

We have also run the detector using a heat trigger. Because the ionization signal is roughly 1000 times faster than the heat signal, it is not always possible to correlate the heat pulse with the charge pulse. Nonetheless, this trigger mode has brought to light a distinct second type of thermal event, one with a rise time dominated by electronics and practically never correlated with an ionization event. For these two reasons we attribute these events to interactions in the heat sensor itself. The energy spectrum of these events is compatible with a small amount of tritium activity in the NTD thermometer.

The events surviving in the neutron zone are not sufficiently separated from the population of electron surface events to be attributed with certainty to nuclear recoil events, induced in the Ge detector by ambient neutrons (the observed rate being much higher than what would be expected from WIMPs). The flux and energy spectrum of the fast ambient neutrons originating for the most part from internal radioactivity in the surrounding rock, have been measured in the LSM by our collaboration using conventional scintillation detectors. The flux was found to be $`4.0\pm \mathrm{0.1\; 10}^6\mathrm{cm}^2.\mathrm{s}^1`$ (statistical error only) . A Monte-Carlo simulation was then performed to predict the expected number of neutron-induced nuclear recoils in our 70 g Ge detector. It uses the GEANT Monte-Carlo simulation  in conjunction with a package specially developed  to handle low energy neutrons (from a few eV to a few MeV). The expected and measured numbers of recoils in the nuclear recoil band (47.5 % acceptance) over the recoil energy interval 15 - 100 keV are given in Table 1.

By inspection of Table 1 it appears that at the -6 V polarization a large fraction, and possibly all the events observed in the nuclear recoil region can be attributed to neutron interactions inside the detector. At -2 V however, the separation between nuclear recoils and electron surface events is, as already noted, slightly less favorable, and accordingly the simulation indicates a pollution of the nuclear recoil band by surface events.

## 5 Conclusion and outlook

Our prototype heat and ionization bolometer provides an excellent separation of the calibrated background and signal down to an energy of 4 keV e.e.. However data taken during background runs contain an important population of off-axis events (electron surface events) essentially not seen during the calibrations. Despite this, the rejection capabilities of the detector remain satisfactory, still reducing the number of electron recoil events by nearly two orders of magnitude. Sensitivity limits on WIMPs  obtained with this prototype are only a factor of 10 higher than the best existing limits obtained with longer-established techniques  despite the small mass (70 g) of our detector and the relatively high radioactivity level in the detector environment. Moreover the present sensitivity of the experiment appears to be limited, principally for the runs with a -6V bias voltage, not by the gamma background or by parasitic electron surface events, but by the ambient neutrons in the underground laboratory. Given the systematic uncertainties in the Monte Carlo neutron simulations a background subtraction is not attempted but, clearly, shielding against neutrons is required.

A new generation of detectors is being tested in the LSM, with new implantation schemes for the electrodes and a rigorous selection procedure of the crystal holder material . Several improvements have been undertaken in the detector environment: nitrogen flushing for radon removal, improved passive protection against gamma radioactivity (including close Roman lead shielding) and 30 cm removable paraffin shielding against neutrons. In an initial test run, two 70 g prototypes show background levels before rejection of $`2\mathrm{event}.\mathrm{keV}^1.\mathrm{kg}^1.\mathrm{day}^1`$ below 50 keV, a factor of 10 improvement over previous data, with a comparable progression for electron surface events. The number of events remaining in the neutron zone is now extremely limited and could be due to residual neutron interactions inside the present shielding. Assuming a good performance of the neutron shielding and that the energy resolution and background rejection are the same as those reported here, a few months of data taking with these detectors should allow a test of the entire parameter region claimed by the DAMA experiment  and the entrance into the relevant supersymmetry parameter space.

## 6 Acknowledgments

The support of the technical staff of the Laboratoire Souterrain de Modane and of the participating laboratories is gratefully acknowledged. We thank J.P. Passerieux, A. Le Coguie and B. Cahan for the quality of their contribution to the drawing and realization of electronic components. This work has been partially funded by the EEC-Network program under contract ERBFMRXCT980167.


format: text

text:
# 1 Introduction

## 1 Introduction

In recent years there has been a persistent interest in obtaining accurate estimates of critical parameters of the three-dimensional (3D) Ising-like systems through high performance simulations and perturbative expansions . Our aim here is to enlarge the knowledge about critical behavior of Ising-like 3D systems by calculating not only the critical exponent of the correlation length $`\nu `$, but also by tackling the much harder problem of calculating the first correction to scaling $`w`$.

A common way to extract information on phase transitions from Monte Carlo simulations is by means of finite-size scaling (FSS), for instance by analyzing the partition function zeros in the complex temperature plane. This approach is not restricted to Ising-like systems and was recently even used to study structural transitions in bio-molecules . However, separation of universality classes can be tricky for 3D systems, and requires high-precision data. An important tool for obtaining such data was provided by Ferrenberg and Swendsen  who revived reweighting techniques introduced by Salsburg et al.  forty years ago. In the sequence, multiple-histogram  and multicanonical simulations  were proposed for a reliable numerical determination of the density of states and extensively checked in two-dimensions where exact results are available. Only by combining one of these sophisticated new simulation techniques (exhaustive multicanonical Monte Carlo simulations for $`L32`$) with Itzykson’s finite-size scaling relation for the partition function zeros , we were able to obtain the results presented in this article.

Before proceeding further, we give the outline of the paper. In the next section we describe the numerical evaluation of the partition function. Section 3 is concerned with the crucial task of data analysis. There we calculate the critical exponent $`\nu `$ and the correction to scaling. The exponent $`\nu `$ was first obtained by a four-parameter fitting. Next, we used the Bulirsch-Stoer algorithm  to extrapolate results from finite lattices to the thermodynamical limit. That algorithm has a free-parameter $`w`$ which is exactly the wanted correction to scaling. We further corroborate on results recently obtained by other techniques. Finally, we test the performance of our approach, and specially the BST algorithm, for the 2D Ising model where the exact zeros are known.

## 2 Multicanonical simulation and partition function zeros

In the multicanonical algorithm  conformations with energy $`E`$ are assigned a weight $`w_{mu}(E)1/n(E)`$. Here, $`n(E)`$ is the density of states. A simulation with this weight will lead to a uniform distribution of energy:

$$P_{mu}(E)\rho (E)w_{mu}(E)=\mathrm{const}.$$
(1)
This is because the simulation generates a 1D random walk in the energy, allowing itself to escape from any local minimum. Since a large range of energies are sampled, one can use the reweighting techniques  to calculate thermodynamic quantities over a wide range of temperatures by

$$<𝒜>_T=\frac{{\displaystyle 𝑑x𝒜(x)w^1(E(x))e^{\beta E(x)}}}{{\displaystyle 𝑑xw^1(E(x))e^{\beta E(x)}}},$$
(2)
where $`x`$ stands for configurations.

It follows from equation 1 that the multicanonical algorithm allows us to calculate estimates for the spectral density:

$$\rho (E)=P_{mu}(E)w_{mu}^1(E).$$
(3)
We can therefore construct the partition function of the three-dimensional Ising model

$$Z(\beta )=\underset{E}{}\rho (E)u^E,$$
(4)
where we define $`u=e^{4\beta }`$.

In the present article, the estimates of the partition function rely on averages over NRUN simulations of NSWEEP Monte Carlo updates for 3D Ising models of linear size $`L`$. In order to allow the system to thermalize, additional sweeps in the beginning were performed and discarded. Table 1 lists the respective values for NSWEEP and NRUN.

Once, we have calculated reliable estimates for the partition function we can calculate the zeros. In the polynomial form, equation (4) has a large number of coefficients for large lattice sizes. For this reason we apply the method presented in  to obtain the leading complex zeros $`u_1^0(L)`$. In table 1 we collect the so calculated zeros as obtained from our multicanonical simulations.

## 3 Finite size scaling analysis

The standard FSS approach for the zero $`u_1^0(L)`$ closest to the real positive axis neglects, for sufficiently large $`L`$, corrections to scaling ,

$$u_1^0(L)=u_c+AL^{1/\nu }\left[1+O(L^w)\right].$$
(5)
Since we are interested in analysing $`\nu `$ in function of the correction to scaling exponent $`w`$, we follow ref.  and fit our data to a four parameter scaling relation,

$$|u_1^0(L)u_c|=a_1L^{1/\nu }+a_2L^{a_3}.$$
(6)
However, due to the small number of available lattices, such a multi-parameter fit is not appropriate. To circumvent this difficulty we proceed including statistics already available for smaller lattice sizes. In ref. , table 1, zeros for 3D Ising model are presented for $`L=6,8,10`$ and $`14`$, together with data from ref.  for $`L=3,4,5,6,8`$ and $`10`$. In addition, we decided to take into account that information and our results in table 1, weighted with the corresponding available precision. The final estimates for each lattice is a linear combination of available data with normalized weight factors which are taken as the reciprocals of the corresponding empirical variances for each data .

In figure 1 we show the corresponding fit for all lattice sizes, whose parameters are found by monitoring the goodness of fit Q . Here we choose the recently obtained critical value $`u_c=0.41204684(25)`$ , to apply the least-square method to equation (6), although the fit is not very sensible to the precision in the value of critical temperature $`u_c`$. The error bar of those data is included, but it is hardly seen in that scale. The goodness of fit $`(Q=0.89)`$ reveals a very good agreement with the data. We obtain $`\nu =0.62853(35)`$ and $`a_3=4.861(84)`$.

If we discard the smallest sizes $`L=3,4`$ and 5, we obtain $`\nu =0.6280(15)`$ ($`Q=0.84`$), corresponding to $`y_t=1/\nu =1.5924(38)`$, which is in remarkable agreement with previous results . However, this fit is less stable with relation to the parameters in the second term. This is because of the presence of rather large lattice sizes.

### 3.1 Correction to scaling from RG transformation

In this section we now evaluate the correction to scaling by a method which is similar to the “finite size phenomenological renormalization group (RG)” analysis by Binder  (which in turn was based on Nightingale’s finite size RG transformation  for the correlation length $`\xi _L`$). Our method to evaluate $`w`$ was previously used  to analyse the 2D Ising model, and it is briefly recalled here.

In order to consider the scaling relation for the longitudinal correlation length $`\xi _L(\beta )`$, we assume that the system is of finite length scale $`L`$ in one direction and infinite in all other directions. The standard expression for the correlation exponent $`\nu `$ is now given by

$$1+\frac{1}{\nu _{L,L^{}}}=\mathrm{ln}\left(\frac{\xi _L^{}/\beta }{\xi _L/\beta }\right)_{\beta _c}/\mathrm{ln}(\frac{L^{}}{L}).$$
(7)
This expression is obtained from a linearization around the fixed critical point $`\beta _c`$ for fixed scaling transformation $`LL^{}`$. The scaling equation for the finite size longitudinal correlation length is given by

$$\xi _L=LY_\xi ((\beta \beta _c)L^{1/\nu },hL^{y_H},\stackrel{~}{u}L^{y_3}).$$
(8)
This differentiable equation includes corrections due to the leading bulk irrelevant scaling field $`\stackrel{~}{u}`$ with exponent $`y_3<0`$, and a magnetic field dependence for the sake of completeness.

From equation (7) and equation (8) one obtains, for $`h=0`$,

$$\frac{1}{\nu _{L,L^{}}}=\frac{1}{\nu }+a_0\frac{L^{y_3}L^{y_3}}{\mathrm{ln}(L^{}/L)}+b_0\frac{L^{2y_3}L^{2y_3}}{\mathrm{ln}(L^{}/L)}+\mathrm{}$$
(9)
where $`a_0`$ and $`b_0`$ include derivatives like $`Y_\xi (y,z)/y|_{y=0,z=0}`$. With the introduction of the rescaling factor $`s=L^{}/L`$ in equation (9), we can now evaluate $`y_3`$.

However an important point remains to be answered: how to estimate finite size dependence of $`\nu `$ on lattice sizes $`L\mathrm{and}L^{}`$ ? This can be achieved from large enough pairs of lattices $`L`$ and $`L^{}`$, with $`L^{}>L`$, through the folowing expression for the partition function zeros:

$$\frac{1}{\nu _{L,L^{}}}=\mathrm{ln}\left(\frac{|u_1^0(L^{})u_c|}{|u_1^0(L)u_c|}\right)/\mathrm{ln}(\frac{L}{L^{}}).$$
(10)
This equation defines our finite size estimators $`\nu _{L,L^{}}`$ from data in table 1. A second estimate can be obtained with the replacement $`|u_1^0u_c|`$ by its imaginary part Im$`(u_1^0)`$ in equation (10). For large enough systems we have Re$`(u_1^0)u_c`$, and both approaches should lead to the same result.

In table 2 we present sequences of those two possible estimates $`\nu _{L,sL}`$ as a function of the fixed rescaling factor $`s=2`$, where the second column stands for the results of equation (10) and the third one for the replacement $`|u_1^0(sL)u_c|/|u_1^0(L)u_c|`$ by $`\mathrm{Im}u_1^0(sL)/\mathrm{Im}u_1^0(L)`$.

As $`L`$ increases, the values obtained by matching pairs of lattices are expected to converge to a limiting value. If we match our largest lattices $`L=24`$ and 32 we obtain from equation (10) $`\nu =0.6260(66)`$ whereas $`\nu =0.6292(70)`$ is obtained by using only the imaginary part of the zeros.

Looking at the values for the crossings $`(L,2L)`$ with $`(12,24)`$ and $`(16,32)`$ in the second column of table 2, we realize that our values for the real part of the corresponding zeros are not precise enough to obtain an ordered sequence towards its critical value $`\nu `$ as $`L\mathrm{}`$. On the other hand the estimate based on the imaginary part of zeros seems to preserve this trend.

Since equations (9) and (10) are valid only for large enough lattice sizes, we have to discard the smallest values in table 2. For this reason we prefer not to evaluate equation (9) by a multi-parameter fit. On the other hand, we note that this relation with $`L^{}=sL`$, is an equation of type

$$T(h)=T+a_1h^w+a_2h^{2w}+\mathrm{},$$
(11)
where we identify,

$`y_3`$ $`=`$ $`w`$
$`h`$ $`=`$ $`1/L`$
$`T(h)`$ $`=`$ $`1/\nu _{L,2L}`$
and $`T=1/\nu `$ the asymptotic limit for $`h0`$. Hence, equation (11) is in the proper form to be analysed by the so called BST approximants, on which we elaborate in the next section.

### 3.2 BST extrapolation

Bulirsch and Stoer  developed an algorithm to extrapolate a sequence $`h_N`$, $`(N=0,1,2,\mathrm{})`$ converging to zero as $`N\mathrm{}`$. See also ref.  for a recent discussion on BST algorithm.

The BST algorithm approximates tabulated data $`T(h)`$ by a sequence of rational functions . The limiting value $`T`$ is computed from a table of recurrent relations defined from:

$`T_1^{(N)}`$ $`=`$ $`0`$
$`T_0^{(N)}`$ $`=`$ $`T(h_N)`$ (13)
and

$`T_m^{(N)}=T_{m1}^{(N+1)}+(T_{m1}^{(N+1)}T_{m1}^{(N)})`$ (14)
$`\left[\left({\displaystyle \frac{h_N}{h_{N+m}}}\right)^w\left(1{\displaystyle \frac{T_{m1}^{(N+1)}T_{m1}^{(N)}}{T_{m1}^{(N+1)}T_{m2}^{(N+1)}}}\right)1\right]^1.`$
Here $`w`$ plays the role of a free parameter. If one defines $`\epsilon _m^{(i)}=2(T_m^{(i+1)}T_m^{(i)})`$, it is expected $`|T_m^{(i)}T|\epsilon _m^{(i)}`$ in the limit $`i\mathrm{}`$. The above remark gives a criterion  for choosing $`w`$ as the value to minimize $`\epsilon _m^{(i)}`$, in order to have a fast and reliable convergence.

Our aim is to extrapolate the finite size sequence in table 2 by the BST algorithm. However, before proceeding with our analysis, we would like to explore first how the BST extrapolant approach performs for exact data. For this end, we come to the 2D Ising model. In ref.  the exact values for $`\nu _{L,2L}`$ are presented. We show in figures 2, 3, and 4 the BST extrapolants obtained from different sets of sequences $`T(h_N)`$.

In figure 2 we show the BST estimates of the critical exponent $`\nu ^{\mathrm{BST}}`$ from a sequence for lattices of lengths $`L=h^1`$, with $`h=1/4,\mathrm{\hspace{0.17em}1}/6,\mathrm{\hspace{0.17em}1}/8,`$ $`1/10,\mathrm{\hspace{0.17em}1}/12,\mathrm{\hspace{0.17em}1}/16,\mathrm{\hspace{0.17em}1}/20,\mathrm{\hspace{0.17em}1}/24`$ and $`1/32`$. This figure presents a pole behavior at $`w1.580`$ for the extrapolated results (full lines) to $`h0(L\mathrm{})`$ for this sequence. In the neighbourhood of that pole we note the corresponding large values for the systematic error (dashed lines) according to the scale at the right hand side of that figure. Here we define the error as the difference between the extrapolated value $`T`$ and a value of last but one interaction: $`TT_{m1}^{(2)}`$.

We also observe that the known value $`\nu =1`$ is obtained for $`w1`$ with error $`0`$. Moreover, it is remarkable that for this sequence of lengths $`h`$, $`\nu ^{\mathrm{BST}}`$ is weakly dependent on $`w`$. For instance, we obtain $`\nu `$ with $`0.1\%`$ precision, $`\nu ϵ[0.999,\mathrm{\hspace{0.17em}1.001}]`$ for a large range of $`w`$ ($`wϵ[0.1888,\mathrm{\hspace{0.17em}1.5482}]`$) before the pole.

In figure 3 we restrict the available sequence to higher lattice sizes, $`h=1/12,\mathrm{\hspace{0.17em}1}/16,\mathrm{\hspace{0.17em}1}/20,\mathrm{\hspace{0.17em}1}/24`$ and $`1/32`$. In that case the extrapolation with $`0.1\%`$ precision is also compatible with a still large range of values for $`w`$: $`wϵ`$$`[0.2816,\mathrm{\hspace{0.17em}1.1969}]`$, before the pole at $`w1.2078`$. In figure 4, we restrict further the sequence to $`h=1/16,\mathrm{\hspace{0.17em}1}/20,\mathrm{\hspace{0.17em}1}/24`$ and $`1/32`$, and we obtain $`wϵ[0.411,\mathrm{\hspace{0.17em}1.3220}]`$, and a pole at $`w1.3871`$.

Therefore, as we restrict our sequences to higher lattice sizes, the effects of the correction to scaling term become more pronounced, leading to a smaller range for acceptable values of $`w`$. This effect has a stronger counterpart in the criterium of minimum error: the acceptable range for $`w`$ is actually narrower than stated above, mainly for figures 3 and 4.

We are finally now at a point where we can use the BST algorithm to analyse our 3D Ising model data of table 2. Since we have to discard smaller lattices in order to accomplish with the expected validity conditions to derive equations (9) and (10), we are left with the sequence in the third column in table 2. Figure 5 presents our so obtained BST estimates for the critical exponent $`\nu `$ from a sequence for lengths $`h=1/6,\mathrm{\hspace{0.17em}1}/8,\mathrm{\hspace{0.17em}1}/12,`$ and $`1/16`$. Since our sequences for $`1/\nu _{L,2L}`$ are restricted to lattice sizes from $`L=6`$ up to $`16`$ we take the four parameter fit estimate $`\nu =0.6280`$ and its statistical error $`0.0015`$ as our input condition to find out $`w`$, as exemplified for the 2D Ising model. This statistical error leads to a range in $`w`$ given by $`w=0.745(74)`$, marking the region of minimum systematic error.

Finally, we remark that this approach yields a value for $`w`$ which is in full agreement with recent results from various other and different approaches: from scaling relations of observables related to the magnetization  it was estimated $`w=0.87(9)`$, from perturbative expansion at fixed $`D=3`$ dimension followed $`w=0.799(11)`$, and $`w=0.814(18)`$ was estimated from $`\epsilon `$-expansion .

## 4 Conclusions

In summary, we have described a new calculation of the correction to scaling exponent $`w`$ based on Monte Carlo multicanonical simulations and finite size scaling theory for the partition function zeros. A four-parameter fitting was done in order to find out the correlation length critical exponent $`\nu `$. Next, the result for $`\nu `$ including the statistical error was used to obtain the acceptable values for the renormalization exponent $`y_3=w`$ by means of the BST algorithm. This algorithm helped us to overcome the difficulties in the straight application of the multiparameter fit (9) to few data points and rather large lattice sizes.

Our results for $`\nu `$ and $`w`$ are in good agreement with recently obtained estimates by Ballesteros et al.  as well as perturbative expansion calculations by Guida and Zinn-Justin .

It is tempting to assume that an accurate value for $`w`$ could be pursued by increasing the significant precision of the complex partition function zeros of the 3D Ising model. This would account for a more precise calculation of $`\nu _{L,L^{}}`$, by evaluating crossings between lattice sizes $`L`$ and $`L^{}`$. However as we have seen from figures 2 to 4, where we used exact values for $`\nu _{L,L^{}}`$, high precision in $`\nu `$ does not necessarily lead to a smaller range in $`w`$. Hence, as an overall conclusion, we note that the large range of $`w`$ is not just a matter of a lack of statistical precision but demonstrates that it is necessary to go to much larger lattice sizes. In particular, for 2D Ising model even $`L=64`$ seems to be not large enough.

Acknowledgements

U. Hansmann gratefully acknowledges support by the National Science Foundation (CHE-9981874), and a generous travel grant by FAPESP (Brazil) which allowed him to visit the Ribeirão Preto Campus of the Universidade de São Paulo, Departamento de Física e Matemática. N.A. Alves and J.R. Drugowich de Felicio acknowledge support by the brazilian agencies FAPESP, CAPES and CNPq.

Table Captions:

Table 1: First partition function zeros for the 3D Ising model on a cubic lattice.

Table 2: Sequence of estimates for the critical exponent $`\nu _{L,2L}`$ from pairs of lattices $`(L,2L)`$. The third column is obtained by replacing $`|u_1^0(2L)u_c|/|u_1^0(L)u_c|`$ by $`\mathrm{Im}u_1^0(2L)/\mathrm{Im}u_1^0(L)`$ in equation (10).

Figure Captions:

Figure 1: Four parameter fit for $`|u_1^0(L)u_c|`$ in function of $`L^{1/\nu }`$ in the range $`L=332`$. The least square method gives $`\nu =0.62853(35)`$ with $`Q=0.89`$.

Figure 2: BST estimates of the critical exponent $`\nu `$ (full lines) and systematic error (dashed lines), as a function of the free parameter $`w`$, for 2D Ising model. The extrapolation is obtained from the finite size sequence with lattice sizes from $`L=4`$ up to $`L=32`$. The right side scale refers to the systematic error.

Figure 3: As in figure 2, BST estimates of the critical exponent $`\nu `$ (full lines) and systematic error (dashed lines) for 2D Ising model, with lattice sizes from $`L=12`$ up to $`32`$. The right side scale refers to the systematic error.

Figure 4: As in figure 2, BST estimates of the critical exponent $`\nu `$ (full lines) and systematic error (dashed lines) for 2D Ising model, with lattice sizes from $`L=16`$ up to $`32`$. The right side scale refers to the systematic error.

Figure 5: BST estimates of the critical exponent $`\nu `$ (full lines) and systematic error (dashed lines), as a function of the free parameter $`w`$ for 3D Ising model. The sequence is obtained with lattice sizes $`L=6,8,12`$ and $`16`$. The right side scale refers to the systematic error.


text:
# The Innermost Stable Circular Orbit of Binary Black Holes

## I INTRODUCTION

Binary black holes are among the most promising sources of gravitational radiation for the new generation of gravitational wave detectors such as the Laser Interferometric Gravitational Wave Observatory (LIGO), VIRGO, GEO and TAMA. This has motivated an intense theoretical effort to predict the gravitational waveform emitted during the inspiral and coalescence of two black holes .

Because of the circularizing effects of gravitational radiation damping, we expect the orbits of close binary systems to have small eccentricities. The inspiral of a binary black hole system then proceeds adiabatically along a sequence of quasicircular orbits up to the innermost stable circular orbit (hereafter ISCO), where the evolution is expected to change into a rapid plunge and coalescence. The ISCO therefore leaves a characteristic signature in the gravitational wave signal, and knowledge of its location and frequency is thus very important for the prospect of future observations.

While various approximations may be adequate to model the adiabatic inspiral up to the ISCO, it is generally expected that only numerical simulations in full general relativity can accurately model the dynamical plunge and merger and predict the gravitational signal from that phase. It is therefore desirable to construct initial data for numerical evolution calculations describing binary black hole pairs at the ISCO, which adds another motivation for determining the location of the ISCO.

Various approaches have been adopted to locate the ISCO in compact binaries, including first order post-Newtonian approximations , variational principles , second order post-Newtonian methods combined with a “hybrid” approach , a Padé approximation  and an effective-one-body approach , and numerical solutions to the constraint equations of general relativity . Unfortunately, however, the results differ significantly and yet have to show any sign of convergence (see Table II below). It would clearly be desirable to understand the origin of these differences. In this paper, we revisit binary black hole solutions to the constraint equations, and evaluate how some of the choices which have to be made in this approach affect the location of the ISCO.

Before the constraint equations of general relativity can be solved, a background geometry and topology have to be chosen. In the conformal-imaging approach adopted by Cook , a conformally flat (spatial) background metric is chosen together with a two-sheeted manifold structure (see Sec. II A). It has been suggested that these choices may affect the location of the ISCO, and may explain the difference between these and the more recent post-Newtonian results.

In this paper, we combine the methods of Cook  and Brandt and Brügmann  to introduce a new approach to constructing binary black holes in quasicircular orbit. We follow Cook  and choose a conformally flat background metric, but do not assume an inversion-symmetry as is done in the conformal-imaging approach. This considerably simplifies the solution of the momentum constraint (see Sec. II B), and produces binary black holes in a three-sheeted manifold structure as opposed to the two-sheeted structure in the conformal-imaging approach. Moreover, adopting the “puncture” approach of Brandt and Brügmann , the Hamiltonian constraint can be solved very easily numerically on R3 without having to impose boundary conditions on interior boundaries (see Sec. II C). We locate the ISCO, and find that its physical parameters agree very well with those found with the conformal-imaging approach of Cook . We therefore conclude that the choice of the underlying manifold structure has a very small effect on the location of the ISCO. Our new approach, which is significantly simpler than the conformal-imaging approach, may also provide a framework in which the conformal-flatness assumption may be relaxed, and its effect on the ISCO be evaluated.

The paper is organized as follows. In Sec. II, we introduce the basic equations and explain how binary black holes in quasicircular orbit can be constructed. We discuss our numerical implementation in Sec. III. In Sec. IV we present our results and compare with those from other approaches. We briefly summarize in Sec. V.

## II SETUP OF THE PROBLEM

### A The Initial Value Problem

A framework for constructing initial data describing binary black holes has been provided by Arnowitt, Deser and Misner’s $`3+1`$ decomposition of Einstein’s equations  and York’s conformal decomposition .

The $`3+1`$ decomposition splits Einstein’s equations into evolution and constraint equations for the metric $`\gamma _{ij}`$ of a spatial hypersurface $`\mathrm{\Sigma }`$, and the extrinsic curvature $`K_{ij}`$, which describes the embedding of the hypersurface $`\mathrm{\Sigma }`$ in the full spacetime. The physical metric $`\gamma _{ij}`$ can now be decomposed into a conformal factor $`\psi `$ and a conformal background metric $`\widehat{\gamma }_{ij}`$,

$$\gamma _{ij}=\psi ^4\widehat{\gamma }_{ij}.$$
(1)
It is also convenient to decompose the extrinsic curvature $`K_{ij}`$ into its trace $`K`$ and a trace-free conformal background extrinsic curvature $`\widehat{A}_{ij}`$ according to

$$K_{ij}=\psi ^2\widehat{A}_{ij}+\frac{1}{3}\gamma _{ij}K.$$
(2)
The Hamiltonian constraint then reduces to an equation for the conformal factor $`\psi `$,

$$8\widehat{}^2\psi \psi \widehat{R}\frac{2}{3}\psi ^5K^2+\psi ^7\widehat{A}_{ij}\widehat{A}^{ij}=0,$$
(3)
and the momentum constraint can be written

$$\widehat{D}_j\widehat{A}^{ij}\frac{2}{3}\psi ^6\widehat{\gamma }^{ij}\widehat{D}_jK=0.$$
(4)
Here $`\widehat{D}_i`$ is the covariant derivative compatible with the conformal background metric, $`\widehat{}^2`$ the Laplacian, and $`\widehat{R}`$ is the Ricci scalar.

Binary black hole initial data cannot be constructed uniquely, because the constraint equations of general relativity determine neither the background geometry nor the topology of the spacetime, both of which have to be chosen before the constraint equations can be solved. Loosely speaking, these ambiguities correspond to different amounts of gravitational radiation in the initial data sets. In this paper, we will follow Cook et.al.  and choose a flat background geometry, but we will choose a three-sheeted topology as opposed to the two-sheeted topology of Cook.

Choosing the background geometry amounts to choosing the conformal background metric $`\widehat{\gamma }_{ij}`$. Following Cook et.al. , we choose the conformal background geometry to be flat so that $`\widehat{\gamma }_{ij}=f_{ij}`$, where $`f_{ij}`$ is the flat metric in a so far arbitrary coordinate system. The covariant derivative $`\widehat{D}_i`$ then becomes the flat-space covariant derivative, and the Ricci scalar $`\widehat{R}`$ vanishes. We will later specialize to cartesian coordinates, $`\widehat{\gamma }_{ij}=\delta _{ij}`$, for which $`\widehat{D}_i`$ reduces to a partial derivative. We also take the hypersurface $`\mathrm{\Sigma }`$ to be maximally embedded in the spacetime so that $`K=0`$. With these choices, the constraint equations simplify to

$$\widehat{}^2\psi =\frac{1}{8}\psi ^7\widehat{A}_{ij}\widehat{A}^{ij}$$
(5)
and

$$\widehat{D}_j\widehat{A}^{ij}=0.$$
(6)
Note that maximal slicing $`K=0`$ automatically decouples the momentum constraint from the Hamiltonian constraint.

Choosing the topology of the spacetime is less straight-forward (compare the discussion in ). Since we are interested in isolated black-holes systems, it is natural to assume the hypersurface $`\mathrm{\Sigma }`$ to be asymptotically flat. Constructing black hole data in vacuum, however, necessarily involves non-trivial topologies. This can be illustrated by a $`t=const`$ slice of the Schwarzschild geometry in isotropic coordinates, where every point inside the black hole’s throat can be mapped into a point outside the throat and vice versa. Moreover, such a mapping can be accomplished with an isometry, which maps the metric into itself, implying that the physical fields at a point inside the throat are identical to those at a point outside the throat. In particular, the geometry near the center is identical to the geometry near infinity. We can therefore think of this solution as describing two identical, asymptotically flat “universes” or “sheets”, which are connected by a throat or Einstein-Rosen bridge .

There is no unique generalization of this topology to the case of multiple black holes . For two black holes, the two throats could either connect to the same asymptotically flat sheet, or else to two separate asymptotically flat sheets. The former approach results in a two-sheeted topology, the latter in a three-sheeted topology.

Cook et.al.  implemented a “conformal-imaging” formalism, which adopts a two-sheeted topology together with the additional demand that the two sheets are related by an isometry so that their physical fields are identical (cf. ). It has been argued that this choice is the “most faithful generalization of the Schwarzschild geometry to the case of multiple holes” . Moreover, the isometry conditions on the throats can be used as boundary conditions in numerical implementations, so that singularities inside the throats can be eliminated from the numerical grid. The computational disadvantage of this method is that boundary conditions have to be imposed on fairly complicated surfaces. In finite difference algorithms, this can be accomplished either with bispherical or Čadež coordinates , designed such that a constant coordinate surface coincides with the throat, or else with fairly complicated algorithms in cartesian coordinates. Both approaches, together with a spectral method, have been compared in .

In this paper, we choose instead a three-sheeted topology and do not assume an inversion-symmetry across the throats, which simplifies the problem in two respects. The analytical solution to the momentum constraint becomes very simple, since we no longer need to construct inversion-symmetric solutions (see Sec. II B). Moreover, the singularities inside the black holes can be removed analytically using a “puncture” method recently suggested by Brandt and Brügmann  (see Sec. II C). The problem can then be solved quite easily on R3 in cartesian coordinates, without having to impose interior boundary conditions. The only added complication is that one now has to locate apparent horizons in the numerically constructed hypersurface.

### B Solving the Momentum Constraint for Binary Black Holes

For maximally sliced hypersurfaces, the momentum constraint decouples from the conformal factor, and analytical solutions to (6) can be given. Moreover, for conformally and asymptotically flat data, the total (physical) linear momentum

$$P^i=\frac{1}{8\pi }_{\mathrm{}}\widehat{A}^{ij}d^2S_j$$
(7)
and the total (physical) angular momentum

$$J_i=\frac{ϵ_{ijk}}{8\pi }_{\mathrm{}}x^j\widehat{A}^{kl}d^2S_l$$
(8)
can be determined from $`\widehat{A}^{ij}`$ without having to solve the Hamiltonian constraint (5) (see ).

Analytical solutions to the momentum constraint (6) describing single boosted or spinning black holes have been given by Bowen and York . A solution $`\widehat{A}^{ij}`$, describing a single black hole at the coordinate location $`𝐂`$ with linear momentum $`𝐏`$ is given by

$$\widehat{A}_{\mathrm{𝐂𝐏}}^{ij}=\frac{3}{2r_𝐂^2}\left(P^in_𝐂^j+P^jn_𝐂^i+(f^{ij}+n_𝐂^in_𝐂^j)P_kn_𝐂^k\right).$$
(9)
Here $`r_𝐂=x^iC^i`$ is the coordinate distance to the center of the black hole and $`n_𝐂^i=(x^iC^i)/r`$ is the normal vector pointing away from that center. Additional terms have to be added in the conformal-imagine approach for an isometry condition to hold across the throat. Note that we have only included linear momentum terms in this expression, and that we are therefore restricting our analysis to non-spinning black holes.

Since the momentum constraint (6) is linear, we can construct binary black hole solutions by superposition of single solutions

$$\widehat{A}^{ij}=\widehat{A}_{𝐂_1𝐏_1}^{ij}+\widehat{A}_{𝐂_2𝐏_2}^{ij}.$$
(10)
From (7) and (8) we find the that total momentum of this solution is $`𝐏=𝐏_1+𝐏_2`$, and the angular momentum about the origin of the coordinate system

$$𝐉=𝐂_1\times 𝐏_1+𝐂_2\times 𝐏_2.$$
(11)
Note that constructing inversion-symmetric solutions for multiple black holes in the conformal-imaging approach is fairly complicated. There, the components of the extrinsic curvature are expressed in terms of an infinite series of recursively defined quantities (see  and Appendix A of ). Relaxing the inversion symmetry, so that the extrinsic curvature can be written as a simple superposition of two solutions (9), therefore greatly simplifies the problem.

### C Solving the Hamiltonian Constraint for Binary Black Holes

Solutions to the Hamiltonian constraint (5) can be constructed by generalizing the Schwarzschild solution in isotropic coordinates for a static (i.e. $`\widehat{A}_{ij}=0`$) and spherically symmetric black hole at coordinate location $`𝐂`$,

$$\psi =1+\frac{}{2r_𝐂}$$
(12)
(note that asymptotic flatness demands $`\psi 1`$ as $`r\mathrm{}`$). Solutions describing multiple static black holes can be constructed by adding contributions $`/(2r_𝐂)`$ for each black hole. To establish an inversion-symmetry, additional terms would again have to be added .

In the “puncture” method suggested by Brandt and Brügmann , a general nonstatic solution to the Hamiltonian constraint is written as a sum of the static, analytic contribution plus a term correcting for finite $`\widehat{A}_{ij}`$. Adopting their notation, we write

$$\psi =u+\frac{1}{\alpha },$$
(13)
where $`\alpha `$ is defined by

$$\frac{1}{\alpha }=\frac{_1}{2r_{𝐂_1}}+\frac{_2}{2r_{𝐂_2}}.$$
(14)
The Hamiltonian constraint then becomes an equation for the correction term $`u`$

$$\widehat{}^2u=\beta (1+\alpha u)^7,$$
(15)
where we have abbreviated

$$\beta =\frac{1}{8}\alpha ^7K_{ij}K^{ij}.$$
(16)
For asymptotic flatness, we impose a Robin boundary condition $`(r(u1))/r=0`$ at large distances from the black holes. The existence and uniqueness of solutions $`u`$ on R3 has been established in . The beauty of this approach is that the poles at the center of the black holes have been absorbed into the analytical terms. The corrections $`u`$ are regular everywhere and can be solved for very easily on a simple computational domain, without having to impose boundary conditions on the throats.

Once the conformal factor $`\psi `$ has been determined, the ADM mass of the solution can be found from

$`E`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}}\widehat{}^i\psi d^2S_i`$ (17)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}}\widehat{}^i\left({\displaystyle \frac{1}{\alpha }}\right)d^2S_i{\displaystyle \frac{1}{2\pi }}{\displaystyle \widehat{}^2u𝑑V}`$ (18)
$`=`$ $`_1+_2+{\displaystyle \frac{1}{2\pi }}{\displaystyle \beta (1+\alpha u)^7𝑑V}.`$ (19)
Note that the integral extends over all space.

### D Constructing Equal-mass Binary Black Holes in Quasicircular Orbit

We now specialize to equal mass black holes with $`_1=_2`$. In the center-of-mass frame of the binary system, we have

$$𝐏𝐏_1=𝐏_2.$$
(20)
Binaries in quasicircular orbit should furthermore satisfy $`𝐏𝐂=0`$, where we have defined

$$𝐂𝐂_1𝐂_2.$$
(21)
Without loss of generality, we can then take $`𝐂`$ to be aligned with the $`z`$-axis, $`𝐏`$ to be aligned with the $`x`$-axis, and place the origin of the cartesian coordinate system at the center between the two black holes.

The problem has now been reduced to a three-dimensional parameter space with the free parameters $``$, $`C𝐂`$ and $`P𝐏`$. For every configuration, we compute several physical quantities. We determine the total ADM mass $`E`$ from (17) and the total angular momentum $`JJ_y=PC`$ from (11). Since we have restricted our analysis to non-spinning black holes, the mass of each individual black hole can be identified with the irreducible mass

$$M=M_{\mathrm{irr}}\left(\frac{A}{16\pi }\right)^{1/2},$$
(22)
where $`A`$ the proper area of the black hole’s apparent horizon . We now define the effective potential as the binding energy

$$E_b=E2M.$$
(23)
Lastly, we compute the proper separation $`l`$ between the two horizons along the line connecting the centers of the two apparent horizons, which is a very good approximation to the shortest proper separation between the two horizons.

Quasicircular orbits can then be found quite easily (see ) by computing the effective potential $`E_b`$ as a function of separation $`l`$ along a sequence of constant black hole mass $`M`$ and angular momentum $`J`$ and locating turning points

$$\frac{E_b}{l}|_{M,J}=0.$$
(24)
A minimum corresponds to a stable quasicircular orbit, while a maximum corresponds to an unstable orbit. For a quasicircular orbit, the binary’s orbital angular velocity $`\mathrm{\Omega }`$ as measured at infinity can then be determined from

$$\mathrm{\Omega }=\frac{E_b}{J}|_{M,l}$$
(25)
(see footnote  for a Newtonian illustration).

## III NUMERICAL IMPLEMENTATION

We adopt a finite difference approach to solve eq. (15) in cartesian coordinates. The numerical code is implemented in a parallel, distributed memory environment using DAGH software , and the Laplace operator in (15) is inverted using PETSc software . We linearize eq. (15) and iteratively solve for corrections to approximate solutions until convergence to a desired accuracy has been achieved. Since the components of the extrinsic curvature (10) are either symmetric or antisymmetric across the coordinate planes $`x=0`$, $`y=0`$ and $`z=0`$, it is sufficient to solve the Hamiltonian constraint in only one octant.

In addition to verifying second order convergence of our code, we have performed tests in the linear regime by comparing with the linear analytic solution for black holes boosted towards each other (eq. (16) in ), and in the nonlinear regime by comparing with the “A2B8” dataset , for which values of the ADM mass in a three-sheeted manifold structure have been given in Table I of . Note, however, that those masses have erroneously been calculated for black hole spins with signs opposite to those given in that table, “$`𝐒_{1,2}=𝐒_{1,2}`$” .

Given a solution $`\psi `$ for a set of parameters $``$, $`C`$ and $`P`$, we can locate an apparent horizon and determine the black hole mass $`M`$ from the horizon’s proper area using the algorithm described in . This algorithm expresses the location of a closed surface in terms of symmetric trace-free tensors, and varies the expansion coefficients until an outer-most trapped surface has been found. We found that for the horizons in this problem, which are fairly spherical and only very mildly deformed, an expansion up to order $`\mathrm{}=4`$ is adequate.

Before constructing sequences of constant black hole mass $`M`$, it is convenient to rescale all variables with respect to that desired value of the black hole mass $`M`$. We introduce the sum of the black hole masses

$$mM_1+M_2=2M$$
(26)
and the reduced mass

$$\mu \frac{M_1M_2}{M_1+M_2}=\frac{M}{2},$$
(27)
and define the dimensionless parameters $`\overline{}/m`$ $`\overline{C}C/m`$ and $`\overline{P}P/\mu `$. We also rescale the angular momentum and angular velocity as $`\overline{J}J/\mu m`$ and $`\overline{\mathrm{\Omega }}m\mathrm{\Omega }`$, and identify the dimensionless effective potential with the rescaled binding energy $`\overline{E}_bE_b/\mu `$.

Sequences of constant $`\overline{J}`$ can now be constructed by setting $`\overline{P}=\overline{J}/\overline{C}`$ for a set of different values of $`\overline{C}`$, and by iterating over $`\overline{}`$ until the (dimensionless) numerical value of the black hole mass $`\overline{M}M_{\mathrm{num}}/M`$ has converged to unity within a desired accuracy for each $`\overline{C}`$ .

For typical cases of interest, the difference between the ADM mass $`E`$ and the sum of the black hole masses $`m`$ is quite small. According to (23), the binding energy is therefore much smaller than those masses, and its relative numerical error much larger. In order to reliably locate a minimum in the binding energy, we therefore have to determine the masses to very high accuracy (cf. ). Moreover, numerous models have to be calculated to construct a sufficient number of sequences over a sufficient range of separations, which is only affordable for a fairly moderate maximum grid size. For a uniform, cartesian grid of a given size, a compromise then has to be found between extending the computational grid to large enough distances and sufficiently resolving the individual black holes.

The location of the outer boundary of the computational grid affects the results through the Robin boundary condition on $`u`$, which is correct only asymptotically, and the energy integral (17), which should extend over all space. The latter effect can be improved by extending the energy integral beyond the numerical grid, where $`\alpha `$ and $`\beta `$ can be evaluated analytically, and where $`u`$ can be estimated from its value on the surface of the computational grid and its $`1/r`$ falloff. In the top panel of Fig. 1, we show the binding energy for different locations of the outer boundary $`\overline{Z}_{\mathrm{out}}Z_{\mathrm{out}}/m`$ for a typical configuration of interest ($`\overline{C}=2.5`$, $`\overline{J}=3.0`$) with a fixed grid resolution $`\overline{h}h/m=0.03125`$. For all calculations presented in this paper we use $`\overline{X}_{\mathrm{out}}=\overline{Y}_{\mathrm{out}}=\overline{Z}_{\mathrm{out}}/2`$. For the dashed line, only contributions to the ADM mass from inside the computational grid have been taken into account, and for the solid line we have expanded the volume for the energy integral by a factor of six in each dimension. Obviously, the latter converges much more rapidly and yields a more accurate value for all locations of the outer boundary.

The resolution of the individual black holes affects the accuracy with which their apparent horizons can be located, and hence the accuracy of their masses $`\overline{M}`$. The effect of the resolution on the binding energy is demonstrated in the bottom panel of Fig. 1, where we show $`\overline{E}_b`$ for different grid resolutions for the same configuration, this time with the outer boundary fixed at $`\overline{Z}_{\mathrm{out}}=6`$.

From Fig. 1, we find that the binding energy $`\overline{E}_b`$ can be determined to within at most a few percent error for $`\overline{Z}_{\mathrm{out}}=6`$ and $`\overline{h}=0.03125`$, corresponding to a numerical grid of size $`96\times 96\times 192`$. The iteration to construct one model then takes approximately 3 CPU hours on the NCSA Origin2000, which makes an extensive survey of parameter space affordable. We use these grid specifications for all results presented in the following section.

## IV RESULTS

In Fig. 2, we show contour plots of the conformal factor $`\psi `$ for a configuration close to the ISCO ($`\overline{C}=2.25`$ and $`\overline{J}=2.95`$). The apparent horizons, marked by the thick dashed lines, are dragged along by the black holes and lag slightly behind in their counter-clockwise orbit. This effect has been discussed for single boosted black holes in . Note that we compute the proper separation $`\overline{l}l/m`$ between the horizons along the line connecting the centers of the apparent horizons. This is a coordinate-dependent quantity, but a very good approximation to the (coordinate-independent) shortest proper separation between the horizons.

We now construct sequences of constant angular momentum for various values of $`\overline{J}`$, and plot the effective potential $`\overline{E}_b`$ along these sequences as a function of $`\overline{l}`$ in Fig. 3. A minimum in the effective potential corresponds to a stable quasicircular orbit. The bold line connecting these minima in Fig. 3 represents a sequence of quasicircular orbits. This sequence terminates at the ISCO, where the adiabatic, quasicircular inspiral of the two black holes is expected to change into a rapid plunge and merger . Since this transition leaves a characteristic signature in the gravitational wave signal, the knowledge of the location and frequency of the ISCO are of great importance for future observations with the new generation of gravitational wave detectors.

In Table I, we list our results for the physical parameters of the ISCO (top line) and compare with the results of Cook  (bottom line). We tabulate the proper separation between the horizons $`\overline{l}`$, the binding energy $`\overline{E}_b`$, the angular momentum $`\overline{J}`$, and the angular velocity $`\overline{\mathrm{\Omega }}`$ as well as the the linear momentum of each black hole $`P/a`$ and the coordinate separation $`C/a`$, where $`a`$ is the (coordinate) radius of the black holes. The latter is not well-defined in our calculation, but since the black holes are nearly spherical it is very reasonable to estimate an average radius from the $`\mathrm{}=0`$ monopole term in the multipole expansion for the horizon.

We conclude that all quantities agree fairly well with those of Cook  within our estimated numerical error of a few percent. As the most significant deviation, we find that in our calculation the binary is slightly more tightly bound at the ISCO, and correspondingly has a slightly larger angular velocity. However, even these quantities differ by less than $`5\%`$, which may be caused by numerical effects. We conclude that the choice of the underlying manifold structure has a very small effect on the location of the ISCO.

In Table II we compare the binding energy, the angular momentum and the angular velocity at the ISCO from various different calculations. For a test particle orbiting a Schwarzschild black hole, the ISCO can be located analytically, which yields $`\overline{E}_b=\sqrt{8/9}10.0572`$, $`\overline{J}=2\sqrt{3}3.464`$, and $`\overline{\mathrm{\Omega }}=1/6^{3/2}0.0680`$. Clark and Eardley  adopted a first order post-Newtonian argument to approximately estimate the location of the ISCO. We list their values for nonrotating neutron stars. Blackburn and Detweiler  adopted a variational principle and assumed a periodic solution to Einstein’s equations. At the ISCO, for which the approximations of this approach fail, they find extremely tightly bound binaries. Kidder, Will and Wiseman  adopted a second order post-Newtonian approximation together with a “hybrid” approach and found an extremely weakly bound ISCO. However, several authors have cast doubt on the robustness and consistency of the hybrid approach . Baumgarte et.al.  constructed fully relativistic models of corotating binary neutron star in quasiequilibrium, albeit assuming conformal flatness, and found that the ISCO depends on the compaction of the neutron stars. In Table I we list their results for $`n=1`$ polytropes of compaction $`M/R=0.2`$, where $`M`$ and $`R`$ are the mass-energy and areal radius which the stars would have in isolation. Naively extrapolating their results to $`M/R=0.5`$ yields values which are very similar to our result for binary black holes. Damour, Iyer and Sathyaprakash  combined a second order post-Newtonian approximation with a Padé approximation, which yields a slightly more tightly bound ISCO than that for a test-particle orbiting a Schwarzschild black hole. A similar result is found by Buonanno and Damour, who combine a second order post-Newtonian approximation with an effective-one-body method .

## V SUMMARY AND DISCUSSION

Since an accurate knowledge of the ISCO is very important for possible future gravitational wave observations, it is very unsettling that different approaches to computing the ISCO lead to very different results (compare Table II). One of these approaches, namely constructing binary black hole solutions to the constraint equations of general relativity, involves choosing the background geometry and topology, and it would be very desirable to know how much the results depend on these choices.

In this paper, we introduce a new method to construct solutions to the constraint equations of general relativity describing binary black holes in quasicircular orbit. We combine the approaches of Cook  and Brandt and Brügmann  to construct binary black holes in a three-sheeted manifold structure, as opposed to the two-sheeted topology in the conformal-imaging approach adopted by Cook . We locate the ISCO and find that its physical parameters are very similar to those found by Cook . Our results confirm those earlier results and imply that the underlying manifold structure only has a very small effect on the ISCO. The latter is perhaps not entirely surprising, since it reflects the fact that the strength of the imaged poles in the conformal-imaging approach is smaller than the strength of the poles themselves .

Our new approach is considerably simpler than the conformal-imaging approach of Cook . The analytic solution to the momentum constraint simplifies because no inversion-symmetric solutions have to be constructed, and the numerical solution to the Hamiltonian constraint simplifies because we can adopt the “puncture” method of Brandt and Brügmann . In particular, we can solve the Hamiltonian constraint in cartesian coordinates on $`R^3`$ without having to impose interior boundary conditions. One disadvantage of our approach is that the apparent horizons have to be located numerically, which we do with the algorithm developed in .

In this paper, we follow Cook  and choose a conformally flat background metric. Accordingly, we cannot address the dependence of the ISCO on the choice of the background geometry. However, since our new method is significantly simpler than the conformal-imaging approach, it may provide a useful framework to relax the assumption of conformal flatness and to construct binary black holes in quasicircular orbit for more general background geometries.

###### Acknowledgements.

It is a pleasure to thank G. B. Cook, S. L. Shapiro and M. Shibata for various very useful conversations. The author gratefully acknowledges the hospitality at the Institute for Theoretical Physics, Santa Barbara, where this research was initiated during the mini-program on “Colliding Black Holes”. Calculations were performed on SGI CRAY Origin2000 computer systems at the National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign. This work was supported by NSF Grants PHY 99-02833 at Illinois and PHY 94-07194 at Santa Barbara.


text:
# Fermi level quantum numbers and secondary gap of conducting carbon nanotubes

## Abstract

For the single-wall carbon nanotubes conducting in the simplest tight binding model, the complete set of line group symmetry based quantum numbers for the bands crossing at Fermi level are given. Besides linear ($`k`$), helical ($`\stackrel{~}{k}`$ and angular momenta, emerging from roto-translational symmetries, the parities of $`U`$ axis and (in the zig-zag and armchair cases only) mirror planes appear in the assignation. The helical and angular momentum quantum numbers of the crossing bands never vanishes, what supports proposed chirality of currents. Except for the armchair tubes, the crossing bands have the same quantum numbers and, according to the non-crossing rule, a secondary gap arises, as it is shown by the accurate tight-binding calculation. In the armchair case the different vertical mirror parity of the crossing bands provides substantial conductivity, though $`k_\mathrm{F}`$ is slightly decreased.

thanks: E-mail: yqoq@afrodita.rcub.bg.ac.yu

Since their discovery  conducting properties of carbon nanotubes have been extensively studied, due to expected technological applications . The electronic bands have been found within the tight binding approximation and assigned either by the linear $`k`$ and the angular $`m`$ quasi momenta quantum numbers , or by the quantum numbers $`\stackrel{~}{k}`$ and $`\stackrel{~}{m}`$ of the helical and ”pure” (i.e. the remnant being not coupled to translations) angular momentum . The both sets of quantum numbers reflect translational and rotational tube symmetries , forming the group $`𝑳^{(1)}`$. The local curvature has been partly neglected in these calculations, assuming the same interaction of carbon atom with its three neighbors. This model predicts  that tube $`(n_1,n_2)`$ is conducting whenever $`n_1n_2`$ is multiple of 3, what is well experimentally justified . It is also suggested  that even in this case the local curvature opens a small secondary gap (less then 0.01 eV and decreasing with the square of the tube diameter ) when $`n_1n_2`$ is non-vanishing; this feature seems to be still beyond the experimental precision.

Here, we present the full set of symmetry quantum numbers of the bands crossing at the Fermi level within the mentioned simple model. They are used to show that the symmetry gives the most profound and easy prediction of the opening of the secondary gap: non-crossing rule prevents conductivity in all but the armchair tubes. The conductivity of the later is due to the extra parity contained in the recently determined  full symmetry group, which makes non-crossing rule inapplicable. Finally, the secondary gap is estimated within a distortion sensitive tight-binding model.

The tube $`(n_1,n_2)`$ is translationally periodic along the tube axis, with period $`a=a_0\sqrt{3q/2n}`$ ($`a_0=2.461`$Å), where $`n`$ is the greatest common divisor of $`n_1`$ and $`n_2`$, $`q=2(n_1^2+n_1n_2+n_2^2)/n`$ and $`=3`$ if $`(n_1n_2)/n`$ is multiple of 3 while $`=1`$ otherwise. The tube is also invariant under rotations for $`2\pi /n`$ around the tube axis. These roto-translational symmetries form the group $`𝑳^{(1)}`$. In addition, chiral ($`𝒞`$) tubes have horizontal $`U`$-axis (rotation for $`\pi `$ around the $`x`$-axis; Fig. 1), doubling the former group:

$$𝑳_𝒞=𝑻_q^r𝑫_n=𝑳^{(1)}+U𝑳^{(1)}.$$
(1a)
Here, the parameter defining the screw-axis (helical) group $`𝑻_q^r`$ is
$$r=\frac{q}{n}\mathrm{Fr}\left[\frac{n}{q}(32\frac{n_1n_2}{n_1})+\frac{n}{n_1}\left(\frac{n_1n_2}{n}\right)^{\phi (\frac{n_1}{n})1}\right]$$

(Fr denotes the fractional part and $`\phi `$ the Euler function). Zig-zag ($`𝒵`$, with $`=1`$) and armchair ($`𝒜`$, $`=3`$) tubes, $`(n,0)`$ and $`(n,n)`$ also have vertical mirror planes $`\sigma _x`$ (the other symmetry elements are generated from the mentioned ones), which doubles the group $`𝑳_𝒞`$ (in these cases $`q=2n`$ and $`r=1`$):
$$𝑳_{𝒵𝒜}=𝑻_{2n}^1𝑫_{nh}=𝑳_𝒞+\sigma _x𝑳_𝒞.$$
(1b)
The additional $`U`$ and $`\sigma _x`$ symmetries yield new conserved quantum numbers: $`+`$ and $``$ ($`A`$ and $`B`$) denote even and odd states with respect to the $`U`$-axis ($`\sigma _x`$ plane). The electronic bands assignation by the complete set of conserved quantum numbers makes easier theoretical predictions of their degeneracies and selection rules relevant for the various physical processes. Here we perform this task for the bands necessary to discuss the conductivity.

The tight binding hamiltonian with single atomic $`p`$-orbital $`i`$ per $`i`$-th site encountering only the nearest neighbors interaction is $`H=\frac{1}{2}_s_{i=1}^3V_iss,i`$, with $`i`$ running over the neighbors of the $`s`$-th C atom (Fig. 1). Neglecting the local curvature, i.e. taking $`V_i=V`$ (estimated between $`3`$eV and $`2.5`$eV), the band structure of such simplified hamiltonian $`H^0`$ can be found easily . The assignation by the full line group quantum numbers has been performed recently . In particular, for the tubes with integer $`(n_1n_2)/3`$, there are pairs of bands crossing at Fermi level $`E_\mathrm{F}=0`$. These bands and their intersection point are characterized by the quantum numbers given in table 1.

For $`𝒞`$ tubes, the simpler to handle with set of quantum numbers, $`\stackrel{~}{k}_\mathrm{F}`$ and $`\stackrel{~}{m}_\mathrm{F}`$, is used. Introducing the angles

$`\stackrel{~}{\psi }_1=\stackrel{~}{k}a{\displaystyle \frac{n_2}{q}}+4r\pi {\displaystyle \frac{2n_1+(1+r)n_2}{3q}},`$ (2a)
$`\stackrel{~}{\psi }_2=\stackrel{~}{k}a{\displaystyle \frac{n_1}{q}}+4r\pi {\displaystyle \frac{(1r)n_1+2n_2}{3q}},`$ (2b)
$`\stackrel{~}{\psi }_3=\stackrel{~}{\psi }_2\stackrel{~}{\psi }_1,`$ (2c)
the Fermi bands for $`𝒞`$, $`𝒵`$ and $`𝒜`$ conducting tubes are:

$`E_𝒞(\stackrel{~}{k})=\pm V\sqrt{3+2{\displaystyle \underset{i}{}}\mathrm{cos}\stackrel{~}{\psi }_i},`$ (3a)
$`G_𝒵(k)=\pm V\sqrt{2(1\mathrm{cos}{\displaystyle \frac{ka}{2}})},`$ (3b)
$`E_𝒜^{A/B}(k)=\pm V(12\mathrm{cos}{\displaystyle \frac{ka}{2}}).`$ (3c)
In the irreducible domain $`\stackrel{~}{k}[0,q\pi /na]`$ of the Brillouine zone , the both bands $`E_𝒞(\stackrel{~}{k})`$ correspond to the same two dimensional irreducible representation $`{}_{\stackrel{~}{k}}{}^{}E_{\stackrel{~}{m}_\mathrm{F}}^{}`$ of the group (1a), . On the contrary, two bands $`E_𝒜^A(k)`$ and $`E_𝒜^A(k)`$ (in the irreducible domain $`k[0,\pi /a]`$) correspond to different two dimensional representations $`{}_{k}{}^{}E_{n}^{A}`$ and $`{}_{k}{}^{}E_{n}^{B}`$ of the group (1b), . As for the $`𝒵`$ tubes, the both bands $`G_𝒵(k)`$ correspond to the four fold representation $`{}_{k}{}^{}G_{2n/3}^{}`$ of the group (1b). These two bands cross at $`k_\mathrm{F}=0`$, giving raise to four fold degeneracy of Fermi level being spanned by two pairs of states carrying representations $`{}_{0}{}^{}E_{2n/3}^{+}`$ and $`{}_{0}{}^{}E_{2n/3}^{}`$ of different $`U`$ parities.

The non-crossing rule  asserts that the bands assigned by the same values of complete set of symmetry based quantum numbers may not cross (or touch). In fact, it can be easily derived from the Wigner-Eckart theorem that in such a case more accurate approach (retaining geometrical symmetry of the model) results in the gap opening at the crossing points. Therefore, in a higher order approximation, $`𝒞`$ tubes attain a gap at $`\stackrel{~}{k}_\mathrm{F}`$; the same argument prevents crossing in any other $`\stackrel{~}{k}`$. Further, although at $`k_\mathrm{F}=0`$ the non-crossing rule at first glance does not obstruct the conductivity of $`𝒵`$ tubes, the finite gap between the bands $`G_𝒵^\pm `$ along the zone implies, by continuity argument, that a gap arises at $`k_\mathrm{F}=0`$, too. On the contrary, the relevant $`𝒜`$ tube bands have different parities along Brillouine zone, being thus beyond the scope of the non-crossing rule. In fact, in all these cases the ”accidental” degeneracy at $`k_\mathrm{F}`$ reflects the increased symmetry tacitly imposed by modeling the interatomic interactions by mutually equal constants $`V_i`$: the hamiltonian became additionally invariant under the permutations of the neighboring atoms. As the following analysis shows, this artefact of the model is to disappear in a more precise approach obeying the original symmetry.

The high symmetry of the nanotubes enables to get the exact dispersion relations for the tight-binding calculations  with arbitrary chosen interaction coefficients $`V_i`$:

$`E_𝒞(\stackrel{~}{k})=\pm \sqrt{{\displaystyle \underset{i}{}}(V_i^2+2{\displaystyle \frac{V_1V_2V_3}{V_i}}\mathrm{cos}\stackrel{~}{\psi }_i)},`$ (4a)
$`G_𝒵(k)=\pm \sqrt{V_1^2+V_2^22V_1V_2\mathrm{cos}{\displaystyle \frac{ka}{2}}},`$ (4b)
$`E_𝒜^{A/B}(k)=\pm (V_32V_1\mathrm{cos}{\displaystyle \frac{ka}{2}}).`$ (4c)
To make the model realistic, the coefficients should be proportional to the overlap integrals  of the site $`p^{}`$ orbitals, meaning that the first order correction to $`V_i`$ is proportional to the bond length $`d_i`$ variation. Especially, $`V_0=V_1<V_2=V_3`$ for $`𝒵`$ and $`V_1=V_2<V_3`$ for $`𝒜`$ tubes. The calculation made with such class of interaction confirm the appearance of the gap (fig. 2). It decreases with tube diameter $`D`$ as $`\beta /D^2`$, with the coefficient $`\beta `$ depending on the chiral angle. The armchair tubes preserve conductivity, since their bands $`E_𝒜^A(k)`$ and $`E_𝒜^B(k)`$ cross in the new point $`k_\mathrm{F}=(2/a)\mathrm{arccos}(V_3/2V_1)<2\pi /3a`$.

To conclude, only for the armchair tubes the bands yielding the conductivity within the simplest tight binding model have different set of the conserved quantum numbers. The difference, the opposite vertical mirror parity, can be observed only with help of the full symmetry line group. In the view of the non-crossing rule, exclusively armchair tubes conductivity is preserved, as illustrated by more exact model. Regularity in the width of the opened gap in $`𝒞`$ and $`𝒵`$ tubes, as well as the shift of $`k_\mathrm{F}`$ in $`𝒜`$ tubes, can be used to resolve some ambiguities in the tubes identification. It should be emphasized that the established complete set of bands quantum numbers enables to calculate exhaustive set of selection rules for the various processes involving interband transitions. In this context it is relevant that the point $`k_\mathrm{F}`$ given in the table 1 becomes weak van Hove singularity for $`𝒞`$ and $`𝒵`$ tubes in the presented model. Concerning $`𝒜`$ tubes, the new crossing point carries nontrivial helical $`\stackrel{~}{k}`$ and total angular momentum $`n\mathrm{}`$, although the ”pure” angular momentum vanishes; this property of the conducting band profoundly confirms proposed chirality of the currents .