text:
# ICANOE and OPERA experiments at the LNGS/CNGS

## 1 INTRODUCTION

Only a few years ago, neutrinos were believed to be massless. Today a picture is emerging in which three (possibly more) neutrino mass eigenstates exist and mix in different amounts to form the flavor eigenstates relevant to the weak interaction. Each year brings its load of new information on neutrinos. But a complete neutrino pattern is still missing. Several new experiments will further fill in the neutrino pattern, the final goal being the comprehensive elucidation of the neutrino masses and mixings. This goal will necessarily be remembered as a fundamental milestone in particle physics, astrophysics and cosmology.

Focusing on experiments detecting neutrinos produced in the Earth atmosphere, the data consistently show that muon-type neutrinos are depleted while passing through the Earth compared to those incident on detectors directly from above. The best agreement with the data is obtained by assuming a muon-type disappearance probability given by

$$P_{disapp}=1\mathrm{sin}^22\theta \mathrm{sin}^2(\beta LE^n)$$
(1)
where $`\beta =\mathrm{\Delta }m^2/4`$ and $`n=1`$. The $`\mathrm{sin}^22\theta `$ is the mixing amplitude squared, $`\mathrm{\Delta }m^2`$ is the difference of the squares of two mass eigenstates, $`L`$ is the path length between the source and the detector and $`E`$ the neutrino energy. The energy of the involved neutrinos and the reconstructed path lengths between production points and detectors can be fitted assuming the disappearance expression of Eq. (1) to extract $`\mathrm{\Delta }m^2`$ and the mixing amplitude. The most stringent constraint comes from SuperK which yields the following 90%C.L. region:

$`\mathrm{\Delta }m^2`$ $``$ $`(26)\times 10^3\mathrm{eV}^2`$ (2)
$`\mathrm{sin}^22\theta `$ $``$ $`0.81`$ (3)
While current data are extremely convincing, in particular that of SuperK, there is a certain number of new measurements that can further improve our understanding of the atmospheric neutrino effect:

1. Only muon “disappearance” has so far been observed; convincing signal for flavor oscillation is the detection of an “appearance” effect; the presence of matter effects disfavors transitions to sterile neutrinos; since maximal $`\nu _\mu \nu _e`$ is excluded by current data, this means detecting $`\nu _\mu \nu _\tau `$ appearance.
2. Given the tau threshold, tau appearance is most easily performed with high energy neutrinos produced at accelerators.
3. Evidence for $`\nu _\mu \nu _\tau `$ appearance in the atmospheric data should also be attempted; this requires large exposures to accumulate enough events at sufficiently high energy, good event reconstruction at high energy and good detector granularity to separate tau decays from backgrounds.
4. The predicted $`L/E`$ dependence of the flavor oscillation has so far been observed with rather poor resolution; in particular, other models of disappearance with “exotic” decays cannot be completely excluded; a convincing signal is the observation of at least a full $`L/E`$ “oscillation” which requires sufficient $`L/E`$ “range” and resolution.

In addition, most of the current analyses have been performed within the so-called “two-family” mixing scenario. Since we know with certainty that three flavor eigenstates exist, we should at least consider the mixing with three mass eigenstates $`m_1`$, $`m_2`$ and $`m_3`$. For the case of atmospheric region, this opens the interesting possibility that oscillations involving tau-type and electron-type take place simultaneously, though with different probabilities. Current best constraints come from reactor experiments reported by Gratta, which limit the element of mixture between electron-type neutrino and $`m_3`$ to a mere $`\mathrm{sin}^22\theta _{13}<0.1`$. One would clearly like to improve this result.

## 2 PROPOSED PROGRAM AT LNGS

Two experiments have been proposed within the context of long-baseline and atmospheric neutrino experiments at CERN: ICANOE and OPERA. A third experiment, MONOLITH, has also been proposed and discussed by Geiser.

The joint ICANOE and OPERA program aims at contributing to the questions described in the introduction. It is really based on these points:

1. a continuation of the observation of atmospheric neutrinos with the improved technique of ICANOE/ICARUS;
2. a sensitive $`\nu _\mu \nu _e`$ and $`\nu _\mu \nu _\tau `$ appearance program with the accelerator neutrinos coming from CERN (CNGS) from a distance of 730 km.

In addition, ICANOE is sensitive to proton decays in the lifetime range of $`>10^{33}`$ years. While this sensitivity is not sufficient to explore convincingly proton decay, ICANOE will be the test-bed to confirm that absolutely background-free searches can be achieved thanks to the excellent imaging of the detector. Hence, ICANOE could open the road to more massive detectors (in the range of 30 kt) capable of exploring proton decay in the range of $`>10^{34}`$ years.

The two proposals, while technologically challenging compared to “traditional” massive neutrino detectors, are both based on many years of R&D undertaken in Europe. They can be considered as natural follow-ups of the NOMAD and CHORUS experiments performed at CERN between 1994 and 1998.

The CERN council has so far only approved in December 1999 the construction of the CNGS beam. The two experiments have not yet been approved, but progress is expected within a short time, since the current schedule foresees the commissioning of the CNGS beam in the year 2005. It is currently planned that experiments should have access to the LNGS hall already in the year 2002 in order to be ready for 2005.

Clearly, the part of the program dedicated to atmospheric neutrinos should start as soon as possible.

## 3 ICARUS/ICANOE

The ICARUS technology provides an excellent general purpose large scale neutrino detector: it offers within the same volume a tracker in three dimensions with high spatial resolution and particle identification, and also, because of the high density of the liquid medium, precise homogeneous calorimetry.

This detector technology combines the characteristics of a bubble chamber with the advantages of the electronic read-out. The $`LAr`$ TPC is based on the fact that ionization electrons can drift over large distances (meters) in a volume of purified liquid Argon under a strong electric field. If a proper readout system is realized (i.e. a set of fine pitched wire grids), it is possible to obtain a massive “electronic bubble chamber” with superb 3D imaging.

The characteristics are the following: (a) it has excellent imaging capabilities, i.e. a “bubble-chamber” like device; (b) the target is fully isotropic and homogeneous; (c) it is a tracking device, capable of dE/dx measurement. The high dE/dx resolution allows both good momentum measurement and particle identification for soft particles; (d) electromagnetic and hadronic showers are fully sampled. This allows to have a good energy resolution for both electromagnetic, $`\sigma (E)/E3\%/\sqrt{E/GeV}`$, and hadronic contained showers, $`\sigma (E)/E15\%/\sqrt{E/GeV}`$; (e) It has excellent electron identification and $`e/\pi ^0`$ discrimination thanks to the ability to distinguish single and double m.i.p. by ionization and to the bubble chamber quality space resolution; (f) calorimetry allows full kinematics reconstruction of “contained” events; (g) muon momentum can be determined by multiple scattering $`\mathrm{\Delta }p/p20\%`$ for long tracks; (h) it is continuously sensitive, self-triggerable, cost effective and simple to build in modular form, sufficiently safe to be located underground.

The ICARUS technique was studied and demonstrated by an extensive R&D program which included ten years of studies on small volumes (proof of principle, purification methods, readout schemes, mixtures of argon-methane, diffusion coefficients, electronics) and five years of studies with several detectors at CERN (purification technology, real events, pattern recognition, event simulations, long duration tests, doping, readout technology). The largest of these devices has a mass of 3 tons. A 50 liters prototype has been exposed to the CERN neutrino and recently tracks with electron drift paths of about $`140\mathrm{cm}`$ have been obtained.

With the cooperation of specialized industries, it is now conceivable to build very large scale devices. AirLiquide has been chosen for the construction of the large cryostat and for the Argon purification system. BREME Tecnica has been selected for the internal detector mechanics. CAEN will produce the readout electronics in large scales. After several years of intense R&D and prototyping, the ICARUS Collaboration is now realizing the first 600 ton module, which will be installed at the Gran Sasso Laboratory in the year 2001.

A recent major step in the R&D program has been the construction and operation of a fully industrial module of 15 ton (T15 prototype). This device has provided (a) a full-scale test of the cryostat technology; (b) a test of the “variable geometry” wire chambers; (c) a test of the liquid phase purification system; (d) a test of the trigger via the scintillation light; (e) a large scale test of the final readout electronics. This test can be considered as the first operation of a 15 ton $`LAr`$ mass as an actual “detector”. More details can be found in Ref..

## 4 OPERA

OPERA is the design of a massive detector able to operate on a medium or long baseline location, to explore $`\nu _\mu \nu _\tau `$ oscillations based on the emulsion technique. In OPERA, emulsions are used as high precision trackers, unlike in CHORUS where they compose the active target. The extremely high space resolution of the emulsion should cope with the peculiar signature of the short lived $`\tau `$ lepton, produced in the interactions of the $`\nu _\tau `$. Since the emulsion does not have time resolution, there are electronic detectors after each module in order to correlate the neutrino interactions to the brick where they occur and to guide the scanning.

## 5 PHYSICS WITH THE CERN BEAM

### 5.1 The CNGS beam

The general strategy was to opt for a wide band neutrino beam based on the experience gathered at CERN with the design and the operation of the WANF.

The primary proton have 400 GeV and assuming 200 days per year at peak intensity 4.8 $`10^{13}`$ per cycle, with 55% overall efficiency, one expects 4.5 $`10^{19}`$ pot/year. In dedicated mode, up to 7.6 $`10^{19}`$ pot/year could be achieved.

The expected event rates per kton and year in shared mode are listed in Table 1.

A first optimization of the beam has been carried out with the goal of maximizing the $`\nu _\tau `$ CC interactions at LNGS for appearance experiments. The optimization of the tau event rates introduces conflicting requirements. Indeed, given the fixed baseline of 730 km and the wish to optimize the rate in the region $`\mathrm{\Delta }m^210^310^2\mathrm{eV}^2`$, the probability of oscillation scales like $`1/E^2`$. But the tau kinematical suppression requires high energy. This rate has been optimized by adjusting the focusing system of the beam and the following rates for tau appearance with maximal mixing $`\mathrm{sin}^22\theta `$ = 1 in shared running mode are obtained in Table 2.

### 5.2 $`\nu _\mu \nu _\tau `$ appearance

It should be stressed that for most of the $`\mathrm{\Delta }m^2`$ region allowed by Superkamiokande, the rate of $`\nu _\tau `$ CC events is so high as to give a statistical excess even prior or with mild selection cuts. This is one of the main reason to perform this experiment at long baseline!

In the case of ICANOE, the expected number of $`\nu _\tau `$ CC with $`\tau e\nu \nu `$ before any selection cuts for $`\mathrm{\Delta }m^2=3.5\times 10^3\mathrm{eV}^2`$ is $`110`$ events in four years of “shared” running, while the background from $`\nu _e,\overline{\nu _e}`$ CC amounts to about 470 events. Such an excess can be seen prior to kinematical cuts, for example in the visible energy distribution of events as shown in Figure 1. The actual cuts will therefore be imposed a posteriori in order to optimize the sensitivity for a given $`\mathrm{\Delta }m^2`$.

Because of the high resolution on measuring kinematical quantities, the $`\nu _\tau `$ appearance search is based on the kinematical suppression of the background using similar techniques to those of the NOMAD experiment. The electron channel in the liquid Argon provides a golden sample. The leading electron is identified and precisely measured. To reconstruct the hadronic jet, the detector is used as an homogeneous calorimeter. The main background of the intrinsic $`\nu _e`$ CC component of the beam is suppressed to a few events, while keeping a $`30\%`$ efficiency for the signal.

In OPERA, the $`\tau `$’s produced in $`\nu _\tau `$ CC interactions, are detected by measuring their decay kink. Current estimates of the number of background events expected yield values below 1 in four years. Thus OPERA is believed to be essentially free of backgrounds.

The number of fully identified tau events from $`\nu _\mu \nu _\tau `$ oscillations as a function of $`\mathrm{\Delta }m^2`$ and for maximal mixing, and the remaining background after selection cuts for ICANOE and OPERA are listed in Table 3. These excellent performances justify the label “direct tau appearance searches”.

### 5.3 $`\nu _\mu \nu _e`$ appearance in three family mixing

We note that for three neutrinos, the three squared mass differences $`\mathrm{\Delta }m_{12}^2`$, $`\mathrm{\Delta }m_{13}^2`$ and $`\mathrm{\Delta }m_{23}^2`$ are linearly dependent, so one has only two free squared mass difference scales. The standard assignment is as such: the smallest one (say $`\mathrm{\Delta }m_{12}^2`$) is attributed to the solar neutrino deficit and $`\mathrm{\Delta }m_{23}^2\mathrm{\Delta }m_{13}^2`$ is set at the atmospheric neutrino anomaly scale.

In the “standard” parameterization of the mixing matrix analogous to the CKM matrix in the quark sector, the amount of $`\nu _e\nu _\mu `$ oscillations is determined by the mixing angle known as $`\theta _{13}`$. Indeed, one has:

$`P(\nu _e\nu _e)`$ $`=`$ $`1\mathrm{sin}^2(2\theta _{13})\mathrm{\Delta }_{32}^2`$ (4)
$`P(\nu _e\nu _\mu )`$ $`=`$ $`\mathrm{sin}^2(2\theta _{13})\mathrm{sin}^2(\theta _{23})\mathrm{\Delta }_{32}^2`$ (5)
$`P(\nu _e\nu _\tau )`$ $`=`$ $`\mathrm{sin}^2(2\theta _{13})\mathrm{cos}^2(\theta _{23})\mathrm{\Delta }_{32}^2`$ (6)
$`P(\nu _\mu \nu _\tau )`$ $`=`$ $`\mathrm{cos}^4(\theta _{13})\mathrm{sin}^2(2\theta _{23})\mathrm{\Delta }_{32}^2`$ (7)
where $`\mathrm{\Delta }_{32}^2=\mathrm{sin}^2\left(\mathrm{\Delta }m_{32}^2L/4E\right)`$ and identical expressions for anti-neutrinos. The negative reactor data sets the limits $`\mathrm{sin}^22\theta _{13}<0.1`$.

The excellent electron identification of ICANOE allows to look for electron excess coming from $`\nu _\mu \nu _e`$ oscillations. At high energy, we expect in addition $`\nu _\mu \nu _\tau `$ with $`\tau e\nu \nu `$ to distort the electron sample. Kinematical analysis can be used to separate the two contributions.

In terms of events, one can see in Table 4 the different contributions expected for a normalization of 4 years running in “shared” mode. Different values of the mixing angles $`\theta _{13}`$ are shown in the first column. The third, fourth and fifth columns list the contribution from $`\nu _e`$ CC (intrinsic beam contamination), the amount of $`\nu _\tau `$ CC with electron decay of the tau for $`\mathrm{\Delta }m^2=3.5\times 10^3\mathrm{eV}^2`$ and maximal mixing, and the amount of $`\nu _\mu \nu _e`$ oscillated events. The last column shows the statistical significance of the $`\nu _\mu \nu _e`$ excess. Hence, ICANOE will test the $`\theta _{13}`$ angle down to a few degrees.

## 6 ATMOSPHERIC NEUTRINOS

ICANOE can give fundamental contributions to the study of atmospheric neutrinos in different aspects, thanks to its unique performances in terms of resolution and precision. In order to achieve a better understanding of neutrino phenomenology, it is fundamental to have coherent results from different kind of measurements, showing that all of them can be interpreted within an unique model. Therefore, the possibility of measuring atmospheric neutrinos with the same detector operating in the long–baseline beam is considered a fundamental feature in order to establish a robust confidence on the results. The perspective, as far as atmospheric neutrinos are concerned, is to provide redundant, high precision measurement and minimize as much as possible the systematics uncertainties of experimental origin which affect the results of existing experiments. Improvements over existing methods are expected in

1. neutrino event selection
2. identification of $`\nu _\mu `$, $`\nu _e`$ and $`\nu _\tau `$ flavors
3. identification of neutral currents

To illustrate the physics reach, we consider the so-called $`L/E`$ analysis. In order to verify that atmospheric neutrino disappearance is really due to neutrino oscillations, an effective method consists in observing the modulation given by the characteristic oscillation probability of Eq. 1. This modulation will be characteristic of a given $`\mathrm{\Delta }m^2`$, when the event rate is plotted as a function of the reconstructed $`L/E`$ of the events when compared to theoretical predictions. The ratio of the observed and predicted spectra has the advantage of being quite insensitive to the precise knowledge of the atmospheric neutrino flux, since the oscillation pattern is found by dips in the $`L/E`$ distribution while the neutrino interaction spectrum is known to be a slowly varying function of $`L/E`$. Such a method is in principle capable of measuring $`\mathrm{\Delta }m^2`$ exploiting atmospheric neutrino events. A smearing of the modulation is introduced by the finite $`L/E`$ resolution of the detection method. Precise measurements of energy and direction of both the muon and hadrons are therefore needed in order to reconstruct precisely the neutrino $`L/E`$. Figure 2 shows the survival probabilities as a function of $`L/E`$ for the neutrino decay hypothesis (see Ref. ), and the oscillation hypothesis with $`\mathrm{sin}^22\theta =1`$ and $`\mathrm{\Delta }m^2=1.0\times 10^3`$, for an exposure of 50 kton$`\times `$year. Both hypothesis are distinguishable from each other at around 2000 km/GeV within the statistical errors.

More details can be found in Ref. . In particular, tau appearance in the atmospheric “beam” is very attractive. Statistically significant signals require exposures of more than $`20\mathrm{kt}\times \mathrm{years}`$, but still are achievable in 4-5 years running.

Finally, ICANOE can also measure upward–going muons from neutrinos interacting in the surrounding rock. In fact, the design and size of ICANOE are such to make it also a large area detector.

## 7 NUCLEON DECAY

Nucleon decay is likely to occur at some level; if experimentally observed it can be used to determine fundamental properties of the nature and the structure of the Unifying Gauge Theory at a scale of $`10^{32}\mathrm{cm}`$. From the experimental point of view, the challenge matches the importance of the issue.

Future nuclear decay experiments have to be able to combine a large mass, the capability of distinguishing between several possible decay channels and a good background discrimination, in order to increase their sensitivity linearly with the mass.

As a specific example of the capabilities of ICANOE, Figure 3 shows a simulated proton decay in the mode $`p\nu K^+`$. This is one of the favored SUSY decays, and is quite typical due to the presence of a strange meson in the final state. A liquid Argon detector can profit from its very good particle identification capabilities to tag the kaon and its decay products. We recall here that the kaon is typically below the Cerenkov threshold for water, therefore it can be seen in a large water detector only via its decay products.

Background-free searches can be performed given the excellent imaging of the detector. Single events will provide evidence for signal. More details can be found in Refs. .

## ACKNOWLEDGMENTS

I thank A. Ereditato for material on OPERA. The help of A. Bueno and J. Rico in preparing this talk is greatly acknowledged. I also thank the organizers of the Neutrino-2000 conference.


format: text

text:
# On pinning and snapping of elastic strands

## Introduction

This paper addresses a seemingly paradoxical result uncovered when assessing the thermodynamic consistency of Yamamoto’s transient network model (Yamamoto 1956,1957,1958, Tanaka and Edwards 1992a-d, Palierne 2000): Consider an elastic strand, i.e. a chainlike molecule or part of molecule, with both ends pinned to fixed sticky sites by reversible bounds. It is fairly obvious that the probability per unit time that a pinned end snaps back and becomes free is an increasing function of the strand extension, hence of the separation of the sticky sites. It is less obvious, however, that the probability per unit time that a free strand becomes pinned when encountering a sticky site should increase in the same proportion with the strand extension. This fact stems from the detailed balance condition for the connection/disconnection reaction formulated in Palierne (2000), the violation of this condition being at the root of the thermodynamic inconsistencies in Yamamoto’s model and in its offspring (Palierne 2000). The argument is as follows.

We take one end of a strand fixed at the origin of space, and consider the other end. Its position is $`𝐡`$, and we consider its pinning to a sticky site, a process hereafter called the connection of a free strand, as well as the reverse process called disconnection, whereby the pinned end snaps back and the connected strand becomes free. Connection and disconnection are instantaneous and take place without change of the connector $`𝐡`$. (hereafter, notations close to those of Palierne 2000 will be used, except for introducing the one-strand distribution $`\psi _\alpha (𝐡)=N^1\mathrm{\Gamma }_{\alpha 𝐡}`$, where $`\alpha =F`$ or $`C`$ denotes the connection state). In the ideal systems we consider, the strand chemical potential $`\mu _\alpha (𝐡)`$ is the sum of the distribution entropy $`k_\mathrm{B}T\mathrm{ln}\psi _\alpha (𝐡)`$ and a part $`\zeta _\alpha (𝐡)`$ depending on the number and the energy of the strand conformations compatible with the $`𝐡`$ value of the connector. The chemical potential at equilibrium is $`\alpha `$ and $`𝐡`$independent, $`\mu _\alpha (𝐡)=\mu ^{eq}`$, and the equilibrium distribution reads

$$\psi _\alpha (𝐡)=Ke^{\zeta _\alpha (𝐡)/k_\mathrm{B}T},K^1=\underset{\alpha =F,C}{}e^{\zeta _\alpha (𝐡)/k_\mathrm{B}T}d^3𝐡$$
(1)
where $`K`$ ensures the normalisation $`(\psi _C(𝐡)+\psi _F(𝐡))d^3𝐡=1`$. The $`𝐡`$-dependence of $`\zeta _\alpha (𝐡)`$ only involves the configuration of the strand, irrespective of its connection state, therefore the connection potential, defined as the difference

$$\mathrm{\Delta }\zeta =\zeta _C(𝐡)\zeta _F(𝐡),$$
(2)
can not depend on $`𝐡`$, and $`\psi _C(𝐡)`$ must be proportional to $`\psi _F(𝐡)`$: $`\psi _C(𝐡)=e^{\mathrm{\Delta }\zeta /k_\mathrm{B}T}\psi _F(𝐡)`$ For instance, freely-jointed strands comprising $`n`$ segments of length $`a`$ have $`\zeta _F(𝐡)=\zeta _C(𝐡)\mathrm{\Delta }\zeta =(3k_\mathrm{B}T/2na^2)𝐡^2`$ and both the free and the connected strands are Gaussian: $`\psi _C(𝐡)\psi _F(𝐡)e^{3𝐡^2/2na^2}`$ (Palierne 2000).

Consider now the kinetics of the connection/disconnection reaction. Let $`\beta _C(𝐡)`$ denote the connection rate, i.e. the probability per unit time that a free strand connects, and let accordingly $`\beta _F(𝐡)`$ denote the disconnection rate. At equilibrium, the distributions $`\psi _F(𝐡)`$ and $`\psi _C(𝐡)`$ of free and connected strands are such that the connection/disconnection reaction rate

$$v(𝐡)\beta _C(𝐡)\psi _F(𝐡)\beta _F(𝐡)\psi _C(𝐡)$$
(3)
vanishes, resulting in the detailed balance condition

$$\frac{\beta _C(𝐡)}{\beta _F(𝐡)}=\frac{\psi _C(𝐡)}{\psi _F(𝐡)}=e^{\mathrm{\Delta }\zeta /k_\mathrm{B}T}\text{independent of }𝐡$$
(4)
The inescapable consequence of this relation is, if we admit that the disconnection rate $`\beta _F(𝐡)`$ increases with $`𝐡`$, then the connection rate $`\beta _C(𝐡)`$ must increase in the same proportion. In the following section, both $`\beta _F(𝐡)`$ and $`\beta _C(𝐡)`$ will be derived from a microscopic model, and relation 4 will be checked.

## microscopic model

We are interested in the interaction between a sticky site and one strand’s end we call the sticker. The other end will therefore be kept fixed at the origin of space. Let $`\{c,𝐡\}`$ denote the set of degrees of freedom, $`𝐡`$ being the position of the non-fixed end and $`c`$ denoting the other degrees of freedom. $`c`$ and $`𝐡`$ are assumed to vary independently and to fluctuate according to an equilibrium Gibbs distribution at temperature $`T`$. The energy of configuration $`\{c,𝐡\}`$ is the sum of the strand internal energy $`ϵ(c,𝐡)`$ and the interaction energy

$$ϵ_s(𝐡)=w\theta (a|𝐡𝐬|)$$
(5)
describing a spherical well centered on $`𝐬`$, of depth $`w`$ and radius $`a`$, attracting the sticker and no other part of the strand ($`w<0`$ for an attracting well). The link between this microscopic model and macroscopic thermodynamics is provided by the quantity

$$\zeta (𝐡)=k_\mathrm{B}T\mathrm{ln}\underset{c}{}e^{\frac{ϵ(c,𝐡)+ϵ_s(𝐡)}{k_\mathrm{B}T}}=\zeta _0(𝐡)+ϵ_s(𝐡)$$
(6)
$$\text{ where }\zeta _0(𝐡)=k_\mathrm{B}T\mathrm{ln}\underset{c}{}e^{\frac{ϵ(c,𝐡)}{k_\mathrm{B}T}}$$
(7)
where $`^c`$ designates teh statistical sum over all configurations keeping a given $`𝐡`$. The conector force

$$𝐟(𝐡)=\frac{\zeta _0(𝐡)}{𝐡}$$
(8)
is the statistical average of the microscopic force $`ϵ(c,𝐡)/𝐡`$ over the Gibbs distribution of $`c`$. The radius $`a`$ will be considered small with respect to the length scale of $`ϵ`$, say $`a\sqrt{ϵ/^2ϵ}`$, but the product $`a𝐟(𝐡)`$ need not be smaller than $`k_\mathrm{B}T`$. Expanding $`\zeta _0(𝐡)`$ about $`𝐬`$ according to $`\zeta _0(𝐡)=\zeta _0(𝐬)+𝐟(𝐬)(𝐡𝐬)`$, we write down for further use the following integrals

$$I_0=e^{\frac{\zeta _0(𝐡)}{k_\mathrm{B}T}}d^3𝐡$$
(9)
$$I(𝐬)=_{|𝐡𝐬|<a}e^{\frac{\zeta _0(𝐡)}{k_\mathrm{B}T}}d^3𝐡=4\pi a^3e^{\frac{\zeta _0(𝐬)}{k_\mathrm{B}T}}(u^2\mathrm{cosh}uu^3\mathrm{sinh}u)$$
(10)
$$J(𝐬)=_{|𝐚|=a}e^{\frac{\zeta _0(𝐬+𝐚)}{k_\mathrm{B}T}}d^2𝐚=4\pi a^2e^{\frac{\zeta _0(𝐬)}{k_\mathrm{B}T}}u^1\mathrm{sinh}u$$
(11)
$$K(𝐬)=e^{\frac{\zeta _0(𝐡)+ϵ_s(𝐡)}{k_\mathrm{B}T}}d^3𝐡=I_0+(e^{\frac{w}{k_\mathrm{B}T}}1)I(𝐬)$$
(12)
where $`u=a|𝐟(𝐬)|/k_\mathrm{B}T`$. $`I_0`$ and $`K(𝐬)`$ are integrals over the whole $`𝐡`$ space, $`I(𝐬)`$ over the sphere $`|𝐡𝐬|<a`$, and $`J(𝐬)`$ over the surface of the same sphere. The equilibrium Gibbs distribution of $`𝐡`$ reads

$$\varphi (𝐡)=\frac{e^{\frac{\zeta _0(𝐡)+ϵ_s(𝐡)}{k_\mathrm{B}T}}}{K(𝐬)}$$
(13)
The probability $`p`$ of finding the sticker within the well can be written

$$p=\frac{e^{\frac{w}{k_\mathrm{B}T}}I(𝐬)}{K(𝐬)}$$
(14)
We now introduce the kinetic properties of our model. The sticker is assumed to diffuse in a Brownian fashion both within and outside the well. The surface $`|𝐡𝐬|=a`$ is assumed to be a potential barrier such that $`\lambda _{\mathrm{out}}`$, the probability per unit time and per unit area that the sticker crosses the barrier from the inside to the outside at point $`𝐬+𝐚`$ reads

$$\lambda _{\mathrm{out}}(𝐬+𝐚)=t_{\mathrm{out}}\varphi _{\mathrm{in}}(𝐬+𝐚)=t_{\mathrm{out}}\underset{\eta 0^{}}{lim}\varphi (𝐬+(1+\eta )𝐚),$$
(15)
$`t_{\mathrm{out}}`$ being the outwards transmission factor, and the probability of the reverse crossing can be written accordingly

$$\lambda _{\mathrm{in}}(𝐬+𝐚)=t_{\mathrm{in}}\varphi _{\mathrm{out}}(𝐬+𝐚)=t_{\mathrm{in}}\underset{\eta 0^+}{lim}\varphi (𝐬+(1+\eta )𝐚),$$
(16)
$`t_{\mathrm{out}}`$ and $`t_{\mathrm{in}}`$ are intrinsic properties of the barrier, independent of the position on the sphere. According to eq. 13, the limits $`\varphi _{\mathrm{in}}`$ and $`\varphi _{\mathrm{out}}`$ are such that

$$\frac{\varphi _{\mathrm{in}}(𝐬+𝐚)}{\varphi _{\mathrm{out}}(𝐬+𝐚)}=e^{\frac{w}{k_\mathrm{B}T}}$$
(17)
independent of the position on the sphere. The principle of detailed balance states that

$$\lambda _{\mathrm{in}}=\lambda _{\mathrm{out}}$$
(18)
at equilibrium, so the transmission factors are linked to the well depth by the customary relation

$$\frac{t_{\mathrm{in}}}{t_{\mathrm{out}}}=e^{\frac{w}{k_\mathrm{B}T}}$$
(19)
The escape probability per unit time of a sticker in a well reads then

$$\mathrm{\Omega }(𝐬)=\frac{1}{p}_{|𝐚|=a}\lambda _{\mathrm{out}}(𝐬+𝐚)d^2𝐚=\frac{t_{\mathrm{out}}J(𝐬)}{I(𝐬)}=\frac{t_{\mathrm{out}}}{a}\frac{u}{(u)}$$
(20)
where $`u=a|𝐟(𝐬)|/k_\mathrm{B}T`$, and $`(u)=\mathrm{coth}uu^1`$ is the Langevin function. One has $`u/(u)=3+u^2/5+\mathrm{}`$ for $`u0`$ and $`u/(u)=u+1+o(1)`$ for $`u\mathrm{}`$ ; and $`u/(u)(1+2(1+u^2/4)^{\frac{1}{2}})`$ over the whole range of $`u`$ to within 3% precision.

## distribution of strands and wells

We model the network as a system of non-interacting strands in a statistically homogeneous distribution of attracting wells. The wells number density $`n`$ is such that many wells enter the volume spanned in the sticker’s Brownian motion. The results of the preceding section need but little modification: using the identity $`I(𝐬)d^3𝐬=\frac{4}{3}\pi a^3I_0`$, the normalisation integral can be rewritten

$$\stackrel{~}{K}=I_0+(e^{\frac{w}{k_\mathrm{B}T}}1)nI(𝐬)d^3𝐬=(1x+xe^{\frac{w}{k_\mathrm{B}T}})I_0$$
(21)
where $`x=\frac{4}{3}\pi na^3`$ is the volume fraction occupied by the wells. The distribution for a given strand becomes

$$\varphi (𝐡)=\{\begin{array}{c}\frac{1}{\stackrel{~}{K}}e^{\frac{\zeta _0(𝐡)}{k_\mathrm{B}T}}\mathrm{if}𝐡\text{outside a well (probability }1x)\\ \frac{1}{\stackrel{~}{K}}e^{\frac{w}{k_\mathrm{B}T}}e^{\frac{\zeta _0(𝐡)}{k_\mathrm{B}T}}\mathrm{if}𝐡\text{inside a well (probability }x)\end{array}$$
(22)
The ensemble of strands is thus described by the connected strand distribution

$$\psi _C(𝐡)=\frac{x}{\stackrel{~}{K}}e^{\frac{w}{k_\mathrm{B}T}}e^{\frac{\zeta _0(𝐡)}{k_\mathrm{B}T}}$$
(23)
and the free strand distribution

$$\psi _F(𝐡)=\frac{1x}{\stackrel{~}{K}}e^{\frac{\zeta _0(𝐡)}{k_\mathrm{B}T}}$$
(24)
satisfying the one-strand normalisation of eq. (1), $`(\psi _C(𝐡)+\psi _F(𝐡))d^3𝐡=1`$. The ratio

$$\frac{\psi _C(𝐡)}{\psi _F(𝐡)}=\frac{x}{1x}e^{\frac{w}{k_\mathrm{B}T}}$$
(25)
takes the form (4) upon introducing the connection chemical potential

$$\mathrm{\Delta }\zeta =w+k_\mathrm{B}T\mathrm{ln}\frac{1x}{x}$$
(26)
as the sum of the energy $`w`$ and the distribution entropy $`k_\mathrm{B}T\mathrm{ln}\frac{1x}{x}`$ accounting for the different volume fraction offered to the connected and to the free stickers. The present model thus satisfies the equilibrium properties announced in the introduction.

The disconnection rate $`\beta _F(𝐡)`$ is the escape probability per unit time for a connected strand, i.e. a strand whose sticker lies in a well centered on $`𝐬`$ such that $`|𝐡𝐬|<a`$. At lowest order in $`a`$, we simply have $`\mathrm{\Omega }(𝐬)=\mathrm{\Omega }(𝐡)`$, and

$$\beta _F(𝐡)=\frac{1}{x}\mathrm{\Omega }(𝐡)$$
(27)
where the fraction $`1/x`$ compensates for the normalisation factor $`x`$ appearing in eq. (23), so that $`\beta _F(𝐡)\psi _C(𝐡)=\mathrm{\Omega }(𝐡)\varphi (𝐡\text{ in a well})`$. The conection rate $`\beta _C(𝐡)`$ relates the probability per unit time of falling into the well to the free strand distribution, so one has similarly

$$\beta _C(𝐡)=\frac{e^{\frac{w}{k_\mathrm{B}T}}}{1x}\mathrm{\Omega }(𝐡)$$
(28)
It must be noted that our $`\beta _F(𝐡)`$ has a much weaker dependence in $`𝐡`$ than the exponential form proposed by Tanaka and Edwards (1992b), whereas their $`\beta _C(𝐡)`$ is a constant, resulting in a connected strand distribution (their eq. 2.23) incompatible with the Gibbs distribution.

Our $`\beta _F(𝐡)`$ and $`\beta _C(𝐡)`$ are found to satisfy the detailed balance (4):

$$\frac{\beta _C(𝐡)}{\beta _F(𝐡)}=\frac{x}{1x}e^{\frac{w}{k_\mathrm{B}T}}$$
(29)
In fact, relation(4) expresses the detailed balance at the global level of connection/disconnection reaction, whereas relation (29) is a consequence of the detailed balance at the local level of crossing the barrier at the well border.

## references

Palierne J.F. (2000) Rheothermodynamics of transient networks. submitted to Rheol. Acta= (2000). posted at http://xxx.lanl.gov/abs/cond-mat/0007368
Tanaka, F. and Edwards S.F.
(1992a) Viscoelastic Properties of Physically Cross-Linked Networks. Transient Network Theory. Macromolecules 25: 1516-1523
(1992b) Viscoelastic properties of physically cross-linked networks Part 1. Non-linear viscoelasticity. J. Non-Newtonian Fluid Mech. 43: 247-271
(1992c) Viscoelastic properties of physically cross-linked networks Part 2. Dynamic mechanical moduli. J. Non-Newtonian Fluid Mech. 43: 273-288
(1992d) Viscoelastic properties of physically cross-linked networks Part 3. Time-dependent phenomena. J. Non-Newtonian Fluid Mech. 43: 289-309
Yamamoto M. J.
(1956) The Visco-elastic Properties of Network Structure I. General formalism. J. Phys. Soc. Jpn 11, 413-421,
(1957) The Visco-elastic Properties of Network Structure II. Structural Viscosity. J. Phys. Soc. Jpn 12, 1148-1158,
(1958) The Visco-elastic Properties of Network Structure III. Normal Stress Effects (Weissenberg Effect). J. Phys. Soc. Jpn 13, 1200-1211


text:
# Density functional formalism in the canonical ensemble

## Abstract

Density functional theory, when applied to systems with $`T0`$, is based on the grand canonical extension of the Hohenberg-Kohn-Sham theorem due to Mermin (HKSM theorem). While a straightforward canonical ensemble generalization fails, work in nanopore systems could certainly benefit from such extension. We show that, if the asymptotic behaviour of the canonical distribution functions is taken into account, the HKSM theorem can be extended to the canonical ensemble. We generate $`N`$-modified correlation and distribution functions hierarchies and prove that, if they are employed, either a modified external field or the density profiles can be indistinctly used as independent variables. We also write down the $`N`$-modified free energy functional and prove that its minimum is reached when the equilibrium values of the new hierarchy are used. This completes the extension of the HKSM theorem.

Density functional theory (DFT) is, undoubtedly, one of the more reliable and established tools in condensed matter theory. It has successfully been used in a great variety of classical systems  as, e.g., uniform and non uniform systems in simple  and general  fluids, confined fluids , melting and freezing , interfaces , etc. as well as in the calculation of electronic properties in all sort of systems . This impressive amount of work directly descends from the pioneering work of Hohenberg, Kohn and Sham  and its extension to non-zero temperature for systems described in the grand canonical ensemble (GCE) by Mermin . Roughly speaking, the Hohenberg-Kohn-Sham-Mermin (HKSM) theorem states that, either the external potential or the density profile can indistinctly be used as independent variables and that the thermodynamic grand potential reaches its minimum value when the equilibrium density profile is used. Therefore, it is in the foundation of all sort of variational principles. The failure to implement a straightforward canonical ensemble (CE) extension is already well known, referred to in , and can be traced back to the fixed N constraint (see eq. (8)). As the GCE is an extension of the CE, an obvious question is where the need for such an implementation exists. In relation to that question, we can mention that experiments done in porous glasses (with a mean pore radius $`r2030\AA `$ ), simulations  and DFT studies  show the interest of having a mesoscopic DFT in the CE in order to study the statistical mechanics of finite closed systems like fluids in spherical pores . In this letter we consider a general mixture of $`p`$ components (labeled by greek indices) with $`N_\alpha `$ particles in a box of volume $`V`$ and under the influence of one body, species dependent, external potentials $`V_\alpha ^{(1)}(𝐱)`$. Following Ramshaw’s work  we will show how, by stripping the canonical correlation functions off of its asymptotic behaviour, an N-modified set of distribution and correlation function hierarchies (which we also show can be defined by functional derivatives of the canonical partition function $`Q_N`$ with respect to a modified external potential) is introduced and it paves the way to prove the HKSM theorem in the CE. More specifically, we do that by introducing an $`N`$-modified free energy which is minimized by the $`N`$-modified density profiles. Here we present the essential results with the pair correlation function case as a prototype and a complete derivation will be published in a paper in preparation.

As usual , the distribution functions and full distribution functions, the truncated correlation functions associated to them and direct correlation functions hierarchies ($`n_𝜶^{(s)}(\{𝐱\})`$, $`\widehat{n}_𝜶^{(s)}(\{𝐱\})`$, $`t_𝜶^{(s)}(\{𝐱\})`$, $`\widehat{t}_𝜶^{(s)}(\{𝐱\})`$, ($`c_𝜶^{(s)}(\{𝐱\})`$ respectively) can be defined by functional derivatives as

$`n_𝜶^{(s)}(\{𝐱\})`$ $`=`$ $`{\displaystyle \frac{_{k=1}^se^{\mathrm{\Phi }_{\alpha _k}(𝐱_k)}}{Q_𝐍}}{\displaystyle \frac{\delta ^sQ_𝐍}{_{k=1}^s\delta e^{\mathrm{\Phi }_{\alpha _k}(𝐱_k)}}}`$ (1)
$`\widehat{n}_𝜶^{(s)}(\{𝐱\})`$ $`=`$ $`{\displaystyle \frac{1}{Q_𝐍}}{\displaystyle \frac{\delta ^sQ_𝐍}{_{k=1}^s\delta \mathrm{\Phi }_{\alpha _k}(𝐱_k)}}`$ (2)
$`t_𝜶^{(s)}(\{𝐱\})`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{s}{}}}e^{\mathrm{\Phi }_{\alpha _k}(𝐱_k)}{\displaystyle \frac{\delta ^s\mathrm{ln}Q_𝐍}{_{k=1}^s\delta e^{\mathrm{\Phi }_{\alpha _k}(𝐱_k)}}}`$ (3)
$`\widehat{t}_𝜶^{(s)}(\{𝐱\})`$ $`=`$ $`{\displaystyle \frac{\delta ^s\mathrm{ln}Q_𝐍}{_{k=1}^s\delta \mathrm{\Phi }_{\alpha _k}(𝐱_k)}}={\displaystyle \frac{\delta \widehat{t}_{\alpha _1\mathrm{}\alpha _{s1}}^{(s1)}(\{𝐱_1,\mathrm{},𝐱_s\})}{\delta \mathrm{\Phi }_{\alpha s}(𝐱_s)}}`$ (4)
$`c_𝜶^{(s)}(\{𝐱\})`$ $`=`$ $`\beta {\displaystyle \frac{\delta ^sF^{exc}}{_{k=1}^s\delta n_{\alpha _k}^{(1)}(𝐱_k)}}={\displaystyle \frac{\delta c_{\alpha _1\mathrm{}\alpha _{s1}}^{(s1)}(\{𝐱_1,\mathrm{},𝐱_s\})}{\delta n_{\alpha s}^{(1)}(𝐱_s)}}`$ (5)
where $`𝜶=(\alpha _1,\mathrm{},\alpha _s)`$, $`\{𝐱\}=(𝐱_1,\mathrm{},𝐱_s)`$ are the species and coordinate sets, $`\beta =(kT)^1`$, $`F^{exc}`$ is the excess free energy and $`\mathrm{\Phi }_\alpha (𝐱)=\beta V_\alpha ^{(1)}(𝐱)`$. The GCE definition is formally the same with the grand partition function $`\mathrm{\Xi }`$ replacing $`Q_N`$. We will consider that the full distribution and correlation functions $`\widehat{n}^{(s)},\widehat{t}^{(s)}`$ are the fundamental quantities (together with $`c^{(s)}`$) in the sense that they not only pinpoint the problems in the CE, is through these functiones that we solve them and no new information can be gotten from the other functions. It is obvious that $`n_\alpha ^{(1)}=\widehat{n}_\alpha ^{(1)}=t_\alpha ^{(1)}=\widehat{t}_\alpha ^{(1)}`$ and that

$$\delta n_\alpha ^{(1)}(𝐱)=\underset{\lambda }{}\widehat{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)\delta \mathrm{\Phi }_\lambda (𝐲)𝑑𝐲$$
(6)
The link with the more usual notation is $`\widehat{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)=t_{\alpha \lambda }^{(2)}(𝐱,𝐲)+n_\alpha ^{(1)}(𝐱)\delta _{\alpha \lambda }\delta (𝐱𝐲)`$ where $`t_{\alpha \lambda }^{(2)}(𝐱,𝐲)=n_\alpha ^{(1)}(𝐱)n_\lambda ^{(1)}(𝐲)h_{\alpha \lambda }^{(2)}(𝐱,𝐲)`$ and, in this case, the independent variables are the external fields. If the density profiles $`\{n_\alpha ^{(1)}\}`$ can also be used as independent variables, then an inverse kernel $`\widehat{t}_{\eta \alpha }^{(2)1}(𝐳,𝐱)=\delta \mathrm{\Phi }_\eta (𝐳)/\delta n_\alpha ^{(1)}(𝐱)`$ must exist such that

$$\underset{\alpha }{}\widehat{t}_{\eta \alpha }^{(2)1}(𝐳,𝐱)\widehat{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)𝑑𝐱=\delta _{\eta \lambda }\delta (𝐳𝐲)$$
(7)
Using the definition of the pair distribution function in both ensembles  , it can be obtained that

$$\widehat{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)d𝐲=\{\begin{array}{c}\text{0}\text{CE}\\ n_\alpha ^{(1)}(𝐱)\delta _{\alpha \lambda }+n_\alpha ^{(1)}(𝐱)\frac{\mathrm{ln}\left[\mathrm{\Xi }n_\alpha ^{(1)}/z_\alpha \right]}{\mathrm{ln}z_\lambda }\frac{\mathrm{ln}\mathrm{\Xi }}{\mathrm{ln}z_\alpha }\frac{\mathrm{ln}\mathrm{\Xi }}{\mathrm{ln}z_\lambda }\hfill & \text{GCE}\hfill \end{array}$$
(8)
Eqs. (7) and (8) are obviously incompatible in the CE while there is no incompatibility in the GCE; because of the fluctuations, the nucleus $`\widehat{t}_{\alpha \lambda }^{(2)}`$ is invertible. Therefore, the Ornstein-Zernike (OZ) equation (obtained when $`c_{\eta \alpha }^{(2)}(𝐳,𝐱)`$ is defined by $`\widehat{t}_{\eta \alpha }^{(2)1}(𝐳,𝐱)=c_{\eta \alpha }^{(2)}(𝐳,𝐱)+\delta _{\eta \alpha }\delta (𝐳𝐱)/n_\alpha ^{(1)}(𝐱)`$) is, mathematically speaking, undefined in the CE and, within that framework, no rigorous DFT is possible. In other words, any variational principle formulated in the CE with the canonical functions hierarchies and having the density profiles as independent variables does not have a mathematically sound foundation. Specifically, this is the problem alluded to in .

The way out we present in this letter is through stripping the correlation functions $`\widehat{t}^{(s)}`$ off of its asymptotic behaviour. This stripping defines the new, $`N`$-modified distribution and correlation functions, $`\stackrel{~}{n}^{(s)},\stackrel{~}{t}^{(s)}`$ respectively. We then show that $`\stackrel{~}{t}^{(s)}`$ is invertible and, from $`\stackrel{~}{t}^{(s)1}`$ we define $`N`$-modified direct correlation functions $`\stackrel{~}{c}^{(s)}`$, write down OZ equations that link $`\stackrel{~}{t}^{(s)}`$ and $`\stackrel{~}{c}^{(s)}`$ and, by proving that they can be also obtained by functional diferentiation of the CE partition function with respect to modified external fields $`\stackrel{~}{\mathrm{\Phi }}_\lambda `$, we conclude that $`\stackrel{~}{n}^{(s)},\stackrel{~}{t}^{(s)},\stackrel{~}{c}^{(s)}`$ constitute a new hierarchical set of distribution and correlation functions in the CE without the mathematical restrictions that the more conventional ones have. Lastly, we write down a $`N`$-modified free energy in terms of these $`N`$-modified functions and prove that it reaches its minimum value when the equilibrium $`N`$-modified density profiles $`\{\stackrel{~}{n}_\alpha ^{(1)}\}`$ are used. This completes the proof of the extension of the HKSM theorem to the CE. Here we will sketch the main steps using the pair correlations case as a prototype and write down the general results. The complete derivation will be reported in a paper in preparation.

The irreducible two bodies behaviour is exclusively contained in $`h_{\alpha \lambda }^{(2)}`$ and we will analyze its asymptotic behaviour. The conditional probability of finding an $`\alpha `$ particle in $`𝐱`$ when a $`\lambda `$ particle is fixed in $`𝐲`$ can be written, when $`𝐱`$ and $`𝐲`$ are very far away, as

$$n_{\alpha \lambda ;\mathrm{}}(𝐱|𝐲)=n_\alpha ^{(1)}(𝐱)+\frac{n_\alpha ^{(1)}(𝐱)}{\rho _\alpha }\mathrm{\Delta }\varphi _\lambda (𝐲)$$
(9)
Here $`\mathrm{\Delta }\varphi _\lambda (𝐲)`$ is an unknown proportionality factor (obviously $`\{\lambda ,𝐲\}`$ dependent); it is not the difference in density due to anchoring a $`\lambda `$ particle. Notice also that the derivative is with respect to $`\rho _\alpha `$, not $`\rho _\lambda `$. The reasons are not only the particle exchange symmetry, they are also related to the fact that one of our goals is to obtain a full hierarchy of $`N`$-modified functions (otherwise, in eq. (16) a modified external field could not be defined). Therefore, the asymptotic behaviour of $`h_{\alpha \lambda }^{(2)}`$ can be written in a fully symmetrical way

$$h_{\alpha \lambda ;\mathrm{}}^{(2)}(𝐱,𝐲)=\mathrm{\Gamma }_{\alpha \lambda }\frac{\mathrm{ln}n_\alpha ^{(1)}(𝐱)}{\rho _\alpha }\frac{\mathrm{ln}n_\lambda ^{(1)}(𝐲)}{\rho _\lambda }$$
(10)
$`\mathrm{\Gamma }_{\alpha \lambda }`$ is an as yet unknown constant. Next, we define an $`N`$-modified correlation function $`\stackrel{~}{h}_{\alpha \lambda }^{(2)}`$ as the correlation function stripped off of its asymptotic behaviour. More precisely,

$$h_{\alpha \lambda }^{(2)}(𝐱,𝐲)=\{\begin{array}{cc}\stackrel{~}{h}_{\alpha \lambda }^{(2)}(𝐱,𝐲)+\mathrm{\Gamma }_{\alpha \lambda }\frac{\mathrm{ln}n_\alpha ^{(1)}(𝐱)}{\rho _\alpha }\frac{\mathrm{ln}n_\lambda ^{(1)}(𝐲)}{\rho _\lambda }\hfill & 𝐱𝐲\hfill \\ \stackrel{~}{h}_{\alpha \lambda }^{(2)}(𝐱,𝐲)=1\hfill & 𝐱=𝐲\hfill \end{array}$$
(11)
It must be noticed that, as the stripping is done through a separation of variables the irreducible two bodies component is not affected, only the long range behaviour (due to the fixed $`N`$ constraint) is isolated. As a consequence, we want to emphasize that $`\stackrel{~}{h}_{\alpha \lambda }^{(2)}`$ is the correlation function with a truly irreducible two-bodies behaviour, not $`h_{\alpha \lambda }^{(2)}`$. Also, the excluded volume effects are not altered and the exchange symmetry imposes that $`\mathrm{\Gamma }_{\alpha \lambda }=\mathrm{\Gamma }_{\lambda \alpha }`$.Therefore, we can define an $`N`$\- modified full truncated correlation function

$$\stackrel{~}{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)=\{\begin{array}{cc}\widehat{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)\mathrm{\Gamma }_{\alpha \lambda }\frac{\mathrm{ln}n_\alpha ^{(1)}(𝐱)}{\rho _\alpha }\frac{\mathrm{ln}n_\lambda ^{(1)}(𝐲)}{\rho _\lambda }\hfill & 𝐱𝐲\hfill \\ \widehat{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)\hfill & 𝐱=𝐲\hfill \end{array}$$
(12)
and verify that it satisfies

$$\stackrel{~}{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)𝑑𝐱𝑑𝐲=\mathrm{\Gamma }_{\alpha \lambda }V^2$$
(13)
This shows that $`\stackrel{~}{t}_{\alpha \lambda }^{(2)}`$ is indeed invertible. Also, $`\mathrm{\Gamma }_{\alpha \lambda }`$ can be written as a functional of $`\stackrel{~}{h}_{\alpha \lambda }^{(2)}`$

$$\mathrm{\Gamma }_{\alpha \lambda }=\frac{1}{V^2}𝑑𝐱n_\alpha ^{(1)}(𝐱)\left[\delta _{\alpha \lambda }+𝑑𝐲n_\lambda ^{(1)}(𝐲)\stackrel{~}{h}_{\alpha \lambda }^{(2)}(𝐱,𝐲)\right]$$
(14)
i.e., a sort of averaged compressibility. This recovers the classical Lebowitz-Percus results  in a form suitable for our purposes. Replacing eq. (12) in (6) and using the chemical potential definition, it can be obtained that

$$\delta n_\alpha ^{(1)}(𝐱)=\underset{\lambda }{}\stackrel{~}{t}_{\alpha \lambda }^{(2)}(𝐱,𝐲)\delta \stackrel{~}{\mathrm{\Phi }}_\lambda (𝐲)𝑑𝐲$$
(15)
where $`\stackrel{~}{\mathrm{\Phi }}_\lambda `$ is the modified external field

$$\delta \stackrel{~}{\mathrm{\Phi }}_\lambda (𝐲)=\delta \mathrm{\Phi }_\lambda (𝐲)+\beta \delta \mu _\lambda $$
(16)
If we write that, by definition,

$$\stackrel{~}{t}_{\eta \alpha }^{(2)1}(𝐳,𝐱)=\frac{\delta \stackrel{~}{\mathrm{\Phi }}_\eta (𝐳)}{\delta \stackrel{~}{n}_\alpha ^{(1)}(𝐱)}=\stackrel{~}{c}_{\eta \alpha }^{(2)}(𝐳,𝐱)+\frac{\delta _{\alpha \eta }\delta (𝐳𝐱)}{\stackrel{~}{n}_\alpha ^{(1)}(𝐱)}$$
(17)
the $`N`$-modified functions $`\stackrel{~}{t}_{\eta \alpha }^{(2)1},\stackrel{~}{t}_{\alpha \lambda }^{(2)}`$ obey eq. (7) and, therefore, $`\stackrel{~}{h}_{\alpha \lambda }^{(2)}`$ and $`\stackrel{~}{c}_{\eta \alpha }^{(2)}`$ are linked by an OZ equation with $`\stackrel{~}{c}_{\eta \alpha }^{(2)}`$ playing the role of an $`N`$-modified direct correlation function.

Now we define an $`N`$-modified intrinsic free energy functional by

$$\beta \stackrel{~}{F}\left[\{\stackrel{~}{n}_\alpha ^{(1)}\}\right]=\beta (K+U\mathrm{ln}P_N)$$
(18)
where $`K,U,P_N`$ are the kinetic, potential energy and $`N`$-bodies probability distribution function in the CE respectively with the average being a CE one and minimize the functional

$$\beta \stackrel{~}{\Omega }\left[\{\stackrel{~}{n}_\alpha ^{(1)}\}\right]=\underset{\alpha }{}\stackrel{~}{n}_\alpha ^{(1)}(𝐱)\stackrel{~}{\mathrm{\Phi }}_\alpha (𝐱)𝑑𝐱+\beta \stackrel{~}{F}\left[\{\stackrel{~}{n}_\alpha ^{(1)}\}\right]$$
(19)
It is easy to prove, using the Gibbs-Bogoliubov inequality , that the minimum is reached when the equilibrium profiles $`\{\stackrel{~}{n}_\alpha ^{(1)}\}`$ are used. This completes the proof of the extension of the HKSM theorem to the CE. It can also be seen that the $`N`$-modified direct correlation function hierarchy $`\{\stackrel{~}{c}^(s)\}`$ can be started by

$`\stackrel{~}{c}_\alpha ^{(1)}(𝐱)=\beta {\displaystyle \frac{\delta \stackrel{~}{F}^{exc}\left[\{\stackrel{~}{n}_\alpha ^{(1)}\}\right]}{\delta \stackrel{~}{n}_\alpha ^{(1)}(𝐱)}}`$

And in particular, we have

$`\stackrel{~}{c}_{\alpha \lambda }^{(2)}(𝐱,𝐲)={\displaystyle \frac{\delta \stackrel{~}{c}_\alpha ^{(1)}(𝐱)}{\delta \stackrel{~}{n}_\lambda ^{(1)}(𝐲)}}={\displaystyle \frac{\delta _{\alpha \lambda }\delta (𝐱𝐲)}{\stackrel{~}{n}_\alpha ^{(1)}(𝐱)}}{\displaystyle \frac{\delta \stackrel{~}{\mathrm{\Phi }}_\alpha (𝐱)}{\delta \stackrel{~}{n}_\lambda ^{(1)}(𝐲)}}`$

Therefore, we have proved that all of our $`N`$-modified functions can be generated by functional differentiation with the formal structure of eqs. (2), (4) and (5).

Summarizing, we have proved that the HKSM theorem of the GCE can be extended to the CE and, in the process of doing that, we have defined a new set of $`N`$-modified distribution and correlation functions that can be built using the same rules and formal structure of the conventional canonical functions with the proviso that the set of modified external fields must be used. A paper with a more detailed derivation together with some applications is in preparation.

We acknowledge support from the National Science Foundation through grants CHE-95-13558, Epscor OSR-94-52893, by the DOE-EPSCoR grant DE-FCO2-91ER75674 and CONICET grant PIP 0859/98.


text:
# A polarized-neutron scattering study of the Cooper-pair moment in Sr2RuO4

## Abstract

We report a study of the magnetization density in the mixed state of the unconventional superconductor Sr2RuO4. On entering the superconducting state we find no change in the magnitude or distribution of the induced moment for a magnetic field of 1 Tesla applied within the RuO2 planes. Our results are consistent with a spin-triplet Cooper pairing with spins lying in the basal plane. This is in contrast with similar experiments performed on conventional and high-Tc superconductors.

Sr2RuO4 has attracted attention since it was discovered to be a superconductor. The superconductivity of this compound is interesting because it is isostructural with the high-$`T_c`$ material La2-xSrxCuO4 and because the superconducting state appears unconventional (i.e. not of the $`s`$-wave singlet type). The low-temperature normal-state of Sr2RuO4 is a quasi-2D Fermi liquid with enhanced quasiparticles. Soon after the discovery of superconductivity in Sr2RuO4, it was suggested that the superconducting state might be unconventional. This suggestion is now supported by a number of experiments including: the observation of a very strong dependence of $`T_c`$ on impurities, a temperature-independent 17O Knight shift on entering the superconducting state, a muon-spin rotation study indicating broken time-reversal symmetry, Andreev reflection, and the observation of power law $`T`$-dependences in the superconducting state for the electronic heat capacity and NQR $`T_1^1`$.

While there is general agreement that the superconducting state of Sr2RuO4 is unconventional, the nature of the wavefunction is still controversial . A knowledge of the spin susceptibility in the superconducting state provides constraints on the pairing wavefunction of a superconductor. Such information can be obtained indirectly by nuclear-resonance techniques through the measurement of the polarization of the $`s`$ electrons on a given site. Alternatively, neutron scattering can directly measure the magnetization density induced by an applied magnetic field. This technique was first used by Shull and Wedgewood to study V3Si and more recently it has been applied to heavy-fermion and high-$`T_c`$ superconductors. In this letter we report a study of the induced magnetization of Sr2RuO4 through the superconducting transition. On entering the superconducting state we find no change in the magnitude or distribution of the induced moment. Our results are in contrast to similar observations on conventional and high-$`T_c`$ superconductors and are consistent with triplet spin-pairing in Sr2RuO4.

The single crystal of Sr2RuO4 used in this study (C117) was prepared by a floating-zone method  in an infrared image furnace. A piece of approximate dimensions 1.5 mm$`\times `$ 2mm $`\times `$ 5 mm was cut using a diamond saw. A.c. susceptibility measurements indicated a sharp superconducting transition with $`T_c`$=1.47 K and $`B_{c2}(T=100\text{mK})`$=1.43 T for $`𝐁`$ $``$ $`[1\overline{1}0]`$. To perform our neutron scattering experiments, the sample was glued to a copper stage using Stycast 2850FT. The copper support was connected to a dilution refrigerator via two 1mm2 diameter copper wires. In the present experiment we applied the magnetic field along the $`[1\overline{1}0]`$ direction. The large anisotropy in $`B_{c2}`$ means that the mutual alignment of the magnetic field and $`[1\overline{1}0]`$ is crucial. Accurate alignment was achieved by mounting the sample and copper stage on a micro-goniometer inside a 2.5 T magnet. In order to verify that the crystal was correctly aligned and at low temperature, we performed an in-situ a.c. susceptibility measurement by mounting two co-axially wound coils near the sample but out of the neutron beam.

The sensitivity of the present magnetic-moment measurement can be dramatically increased by the use of polarized neutrons. Our measurements were performed on the IN20 spectrometer at the Institut Laue-Langevin, Grenoble. A beam of neutrons with a polarization greater than 93% and with energy $`E_i`$= 34.8 meV was produced using the (111) Bragg reflection of a Heusler monochromator. We searched for possible depolarization of the beam caused by the presence of the vortex lattice by measuring the polarization of neutrons scattered by the (002) Bragg reflection using a Heusler analyzer. However, the presence of the vortex lattice produced no detectable depolarization for fields B $``$ 1 T under the experimental conditions used.

In the present experiment we measure the magnetization density $`𝐌(𝐫)`$ in the presence of an applied magnetic field $`𝐁`$. Because of the periodicity of the crystal, an applied magnetic field induces a magnetization density with spatial Fourier components $`𝐌(𝐆)`$ corresponding to reciprocal lattice vectors $`𝐆`$, where

$`𝐌(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{\nu _0}}{\displaystyle \underset{𝐆}{}}𝐌(𝐆)\mathrm{exp}(i𝐆𝐫)`$ (1)
and $`\nu _0`$ is the unit-cell volume. The Fourier components of the magnetization density are given by

$`𝐌(𝐆)`$ $`=`$ $`{\displaystyle _{\text{cell}}}𝐌(𝐫)\mathrm{exp}(i𝐆𝐫)𝑑𝐫.`$ (2)
A diffraction experiment allows these spatially varying components of the magnetization to be measured, even in the superconducting state.

Neutrons interact with condensed matter both through the strong nuclear interaction and through the electromagnetic interaction. If the neutron momentum transfer $`𝜿=𝐤_i𝐤_f`$ equals a reciprocal lattice vector $`𝐆`$, i.e. we satisfy the Bragg condition, then scattering occurs both because of the periodicity of the nuclear density and because of the microscopic periodicity of the magnetization density. The two scattered waves interfere. For neutrons with initial and final spin polarizations $`𝝈_i`$ and $`𝝈_f`$, the total cross section is,

$$\frac{d\sigma }{d\mathrm{\Omega }}_{\sigma _i\sigma _f}\left|\sigma _i|\frac{\gamma r_0}{2\mu _B}𝝈\widehat{𝜿}\times \left\{𝐌(𝜿)\times \widehat{𝜿}\right\}+F_N(𝜿)|\sigma _f\right|^2$$
(3)
where $`\gamma r_0`$= 5.36$`\times 10^{15}`$ m, $`F_N(𝜿)`$ is the nuclear structure factor (in units of m f.u.-1) and $`𝐌(𝜿)`$ has units of $`\mu _B`$ f.u.-1. The sign of the magnetic term in Eq. 3 is controlled by the polarization of the neutrons. By reversing the polarization of the incident neutrons we are able to isolate the term due to the interference between magnetic and nuclear contributions, yielding the magnetization density. Experimentally we measure the flipping ratio $`R`$, defined as the ratio of cross sections for initial neutron-spin states which are parallel or antiparallel to the applied magnetic field and with arbitrary final spin state. In the present experiment, the applied field, and hence neutron polarization, is perpendicular to the scattering vector $`𝜿`$. Under these circumstances the flipping ratio $`R`$ evaluated from Eq. 3 is only sensitive to the component of magnetization parallel to the applied field, $`M_{}(𝜿)`$. Because the induced moment is small, the experiment is carried out in the limit $`(\gamma r_0/2\mu _B)|𝐌(𝜿)|/|F_N(𝜿)|0.0011`$. In this limit, the flipping ratio derived from Eq. 3 is,

$$R=\frac{|F_N(𝜿)(\gamma r_0/2\mu _B)M_{}(𝜿)|^2}{|F_N(𝜿)+(\gamma r_0/2\mu _B)M_{}(𝜿)|^2}1\frac{2\gamma r_0}{\mu _B}\frac{M_{}(𝜿)}{F_N(𝜿)}.$$
(4)
As the nuclear structure factors $`F_N(𝜿)`$ are known from the crystal structure, Eq. 4 directly gives the magnetization $`𝐌(𝐆)`$. Fig. 1 shows the susceptibility $`\chi (𝜿,0)=𝐌(𝜿)/B`$ for a number of wavevectors $`𝜿`$ determined from the measured flipping ratios $`R`$ and converted into magnetic moment using Eq. 4. The measured flipping ratios were corrected for imperfect beam polarization and finite flipper efficiency. Other possible corrections including extinction, absorption, incoherent scattering, the neutron spin-orbit interaction  , half-wavelength contamination in the incident beam were estimated to be small. In order to allow us to enter the mixed state, measurements were generally made with an applied field of 1 Tesla. However, Fig. 1 shows that measurements in the normal state at $`𝜿=(002)`$ and $`B`$=2 Tesla yielded consistent results.

Neutron scattering measures the total moment, thus Fig. 1 includes contributions from diamagnetic, orbital and spin components of the susceptibility. The fall off of the susceptibility $`\chi (𝜿,0)`$ with increasing $`|𝜿|`$ is due to the finite extent of the induced moment in space. Fig. 1 shows the Ru form factor (solid line)  scaled to the measured bulk susceptibility  of Sr2RuO4, $`\chi _{ab}`$=0.9$`\times 10^3`$ e.m.u. mole-1. The poor agreement of the $`𝜿`$=(110) component with the Ru form factor may suggest the presence of a significant induced moment on the oxygen atoms.

Before discussing our results in the superconducting state of Sr2RuO4, we will briefly discuss the same measurement of the spin susceptibility in the conventional superconductor V3Si. In a conventional superconductor with spin-singlet pairing, the spin susceptibility is suppressed on entering the superconducting state because electrons with anti-parallel spins pair up. For $`BB_{c2}`$ the temperate dependence is described by the Yosida function . Wedgewood and Shull  performed a polarized-neutron measurement of the susceptibility on the conventional superconductor, V3Si, and observed a reduction of the susceptibility, due to the formation of spin singlets ($`S=0`$) \[see Fig. 2(a)\]. We have reproduced the Wedgewood-Shull result using the same experimental set-up as for our Sr2RuO4 measurements. Our results are consistent with Wedgewood and Shull and are shown as open circles in Fig. 2(a).

Having observed the induced moment in the normal state of Sr2RuO4, we proceeded to investigate the effect of the superconductivity. Because of the low $`T_c`$ and strongly anisotropic $`B_{c2}`$ it was important to verify that the sample was in good thermal contact with the dilution refrigerator and well aligned with the applied magnetic field. Fig. 2(c) shows an in-situ measurement of the a.c. susceptibility, made using the balance output of an inductance bridge, plotted against applied field $`B`$ for T=100 mK. The kink corresponds to $`B_{c2}`$=1.43 Tesla, which is indistinguishable from the published value , demonstrating that the sample was well-aligned. In order to ensure penetration of the magnetic field throughout the sample the sample was always “field-cooled”.

Fig. 2(b) shows the temperature dependence of the induced moment corresponding to the (002) Bragg position for $`B`$=1 Tesla. This component was chosen for detailed study because its amplitude is proportional to the sum of the moments induced on the in-plane ruthenium and oxygen atoms $`(m_{\text{Ru}}+2m_\text{O})`$. For a 1 Tesla field applied parallel to $`[1\overline{1}0]`$, $`T_c=`$1 K (See Fig. 2(d)). On entering the superconducting state, we find that there is no change in this component of the induced moment within the experimental error. In contrast to the V3Si measurement, we investigated Sr2RuO4 at relatively high fields, $`B/B_{c2}=0.68`$, thus the presence of normal vortices leads to a significant density of quasiparticles and finite spin susceptibility in the mixed state. Using the measured linear heat capacity in the superconducting state we estimate the zero-temperature spin susceptibility $`\chi (T0,B=1T)=0.45\chi _{\text{normal}}`$. The dashed line in Fig. 2(b) is a Yosida function modified to include the finite susceptibility in a field: this prediction is still at variance with the data. Thus, the absence of a change in spin susceptibility is not compatible with spin-singlet or even-parity pairing.

The absence of a change in the spin susceptibility can be explained if Cooper pairs form from electrons with parallel spins. Such “equal-spin pairing” (ESP) was first proposed in the context of 3He by Anderson and Morel . Within an ESP scenario the superconducting state is a superposition of the two possible ($`S=1`$) parallel paired states. In an applied magnetic field, one state is favored yielding a net spin moment and the same susceptibility as in the normal state . An EPS-type pairing implies an odd-parity or spin-triplet state, thus the present experiment supports the notion that the superconducting wavefunction in Sr2RuO4 has an odd-parity representation. There are five unitary odd-parity representations ($`\mathrm{\Gamma }_{15}^{}`$) of the order parameter for the crystal point group $`C_{4h}`$ . Of the allowed representations only the degenerate $`\mathrm{\Gamma }_5^{}`$ ($`E_u`$) or $`𝐝(\widehat{𝐤})=\widehat{𝐳}(\widehat{k}_x\pm i\widehat{k}_y)`$ state is expected to show no change in its spin susceptibility for magnetic fields applied in the basal plane.

So far we have discussed the spin part of the Cooper-pair wavefunction. The $`\mathrm{\Gamma }_5^{}`$ state is special in that the orbital part of the wavefunction suggests that the paired electrons have relative angular motion with orbital angular momentum $`L_z=\pm 1`$. $`\mu `$SR measurements  reveal a spontaneous internal magnetic field in the superconducting state which is thought to be associated with the internal orbital moment. The present experiment is sensitive to the total electronic moment ($`𝐒+𝐋`$) in the superconducting state. Under the present experimental conditions, we measure magnetic moments parallel to the applied field. The absence of a change in the orbital moment measured by the present experiment is entirely consistent with the $`\mathrm{\Gamma }_5^{}`$ assignment above. Firstly, because the moment is expected to be parallel to the $`c`$-axis and, secondly, because the bulk orbital moment is expected to be small.

At first sight our results appear to contradict recent heat capacity  and other experiments  which suggest that the superconducting gap is strongly anisotropic, possibly having nodes for certain directions. However, a strongly anisotropic gap function is still allowed within the $`\mathrm{\Gamma }_5^{}`$ representation. For example, the anisotopic $`f`$-states which have recently been proposed  still have the $`\mathrm{\Gamma }_5^{}`$ representation and are consistent with our interpretation.

Our results complement and contrast with recent NMR measurements  of the 17O Knight shift in the mixed state of Sr2RuO4. These do not detect a reduction in the spin susceptibility. The Knight shift measures the polarization of the $`s`$-electrons at a given site: electrons in other orbitals and on other sites are probed because of the overlap of orbitals. In contrast, the present neutron-scattering measurement directly measures the spatially-averaged total moment.

Measurements of the induced moment in the mixed state have also been performed on other unconventional superconductors. In the high-temperature superconductor YBa2Cu3O$`_{\text{6+x}}`$  a suppression of the spin susceptibility is observed which is consistent with an even-parity or singlet pairing. In contrast, the heavy-fermion superconductors UPt3 and UBe13  show no reduction in the spin susceptibility on entering the superconducting state suggesting that they have odd-parity pairing.

In summary, we have used a spin-polarized neutron-scattering technique to measure the magnetization in the mixed state of the unconventional superconductor Sr2RuO4. We find that for a 1 Tesla field applied parallel to the $`[1\overline{1}0]`$ basal-plane direction, there is no detectable change in the component of the moment parallel to the applied field. Our results strongly support the identification of the paired state of Sr2RuO4 as the $`\mathrm{\Gamma }_5^{}`$ or $`𝐝(\widehat{𝐤})=\widehat{𝐳}(\widehat{k}_x\pm i\widehat{k}_y)`$ state.

We are grateful to D. Bintley, N. Clayton, and J.-L. Raggazoni for their help in preparing this experiment and G. Aeppli, J. Annett, E. M. Forgan and A. P. Mackenzie for useful discussions.


text:
# Simple excision of a black hole in 3+1 numerical relativity

\[

## Abstract

We describe a simple implementation of black hole excision in 3+1 numerical relativity. We apply this technique to a Schwarzschild black hole with octant symmetry in Eddington-Finkelstein coordinates and show how one can obtain accurate, long-term stable numerical evolutions.

\]

The simulation of a black hole inspiral collision is one of the most important open problems facing numerical relativity. Traditional techniques using singularity avoiding slicings will not be able to follow such a collision since problems associated with the stretching of the slice typically cause simulations to crash or to become extremely inaccurate in time scales far shorter than the orbital time scale. Black hole excision techniques (also known as “apparent horizon boundary condition” ) appear to be the most promising way of eliminating the problem of the slice stretching, thus in principle allowing numerical simulations to follow the inspiral from well separated holes through the merger and the ring-down phase.

Black hole excision was first attempted successfully by Seidel and Suen in spherical symmetry , and was later studied in more detail by Anninos et.al. . However, the original idea is older, and Thornburg  has attributed it to a suggestion by Unruh from 1984. The idea consists of two parts: First, one places a boundary inside the black hole and excises its interior from the computational domain. Second, one uses a shift vector that keeps the horizon roughly in the same coordinate location during the evolution (“horizon tracking”, see ). Since no information can leave the interior of the black hole, excision should have no effect on the physics outside. Ideally, one would like to know the position of the event horizon which marks the true causal boundary, but the global character of its definition means that in principle one can only locate it once the whole evolution of the spacetime is known. The apparent horizon, on the other hand, can be located on every time slice and is guaranteed to be inside the event horizon. In practice one therefore needs to find the apparent horizon and excise a region contained inside it.

Though black hole excision has been successful in spherical symmetry , it has been difficult to implement with a 3+1 approach in three-dimensions (3D) , where instabilities typically plague the evolutions (but some progress has been made, see ). Black hole excision using a characteristic formulation, on the other hand, has been very successful in 3D, allowing stable evolutions of perturbed black holes for thousands of $`M`$’s . However, such characteristic formulations are likely to have problems with the development of caustics in the case of extremely distorted or colliding black holes, so the search for a stable 3+1 excision implementation is still of great importance.

Here we present a 3+1 approach to black hole excision in 3D that has allowed us to obtain long-term stable, accurate evolutions of a single black hole spacetime. These results are currently limited to simulations in octant symmetry as discussed below.

I. SIMPLE BLACK HOLE EXCISION

Though conceptually simple, black hole excision in 3D is a complicated problem numerically. First, one has to cut a hole in the computational domain that has a spherical topology and is therefore not well adapted to the Cartesian coordinates typically used. Second, one has to apply some condition at the boundary of the excised region that is stable and respects the causality of the physical system. As the excised region is inside a black hole, no boundary condition should be needed since all the information required to update the boundary comes from outside the excised region. However, achieving this “boundary without a boundary condition” (BWBC)  in 3D is difficult, particularly if one uses a formulation of the evolution equations that is not hyperbolic. The way this problem is usually approached is by using “causal differencing”  or “causal reconnection” , where the computational molecules are adapted to follow the causal structure. The mixture of these issues makes it difficult in practice to identify what particular element of an algorithm is responsible for causing a numerical simulation to go unstable.

In our approach we have simplified the algorithm as much as possible, separating out what we believe is essential to the excision problem. Our algorithm is based on the following simplifications:

* Excise a region adapted to Cartesian coordinates, i.e. excise a cube contained inside the horizon.
* Do not attempt to fulfill the BWBC ideal, and use instead a simple but stable boundary condition at the excision boundary.
* Do not use causal differencing. Use instead centered differences in all terms except the advection terms on the shift (terms that look like $`\beta ^i_i`$). For these terms use upwind along the shift direction (we use the standard 1D second-order upwind stencil in each of the Cartesian coordinate directions based on the sign of the corresponding shift component at each point). This is very important, as it is the only place where any information about causality (i.e. the direction of the shift) enters our scheme. Using a centered approximation for these terms results in an unstable scheme.

One can worry that excising a cube will introduce artifacts into the evolution, but as long as the boundary condition used at the sides of the cube is consistent those artifacts will converge away with increased resolution. Similarly, one can argue that applying a boundary condition instead of using causal differencing is inconsistent with the physics, but since this condition is applied well inside the horizon, any error introduced is unlikely to propagate outside the hole.

II. STATIC BLACK HOLE SPACETIME

As the first test of our excision algorithm we have considered a single static black hole written in “3+1 Eddington-Finkelstein”(EF) coordinates. These 3+1 EF coordinates are a simple transformation of the standard ingoing EF coordinates  to a 3+1 form. The resulting metric has no coordinate singularities, penetrates the event horizon, reaches the physical singularity, and is manifestly time independent. This makes it ideal for excision tests where one can excise the physical singularity and try to keep the numerical evolution stable and close to static. The 3+1 EF metric has the form

$`ds^2`$ $`=`$ $`\left(12M/r\right)dt^2+\left(4M/r\right)dtdr`$ (2)
$`+\left(1+2M/r\right)dr^2+r^2d\mathrm{\Omega }^2.`$
with $`M`$ the black hole mass and $`d\mathrm{\Omega }`$ the solid angle element. From this metric one can read the values of the 3-metric, lapse and shift. The extrinsic curvature can then be obtained in a straightforward way.

III. EVOLUTION EQUATIONS

Formulation— We comment briefly on the formulation used for the simulations described below. Our simulations have been performed using a formulation of the 3+1 evolution equations developed by Baumgarte and Shapiro  (BS), based on previous work of Shibata and Nakamura  (SN). The motivation for using this BSSN formulation comes from the fact that it has shown remarkable stability properties when compared to the Arnowitt-Deser-Misner (ADM) formulation  in a wide range of numerical simulations .

The BSSN variables are defined in terms of the spatial metric $`\gamma _{ij}`$ and the extrinsic curvature $`K_{ij}`$ as: $`\varphi =\text{ln}(\text{det}\gamma _{ij})/12`$, $`\stackrel{~}{\gamma }_{ij}=e^{4\varphi }\gamma _{ij}`$, $`\text{tr}K=\gamma ^{ij}K_{ij}`$, $`\stackrel{~}{A}_{ij}=e^{4\varphi }(K_{ij}\gamma _{ij}\mathrm{tr}K/3)`$, and $`\stackrel{~}{\mathrm{\Gamma }}^i=\stackrel{~}{\gamma }^{jk}\mathrm{\Gamma }_{jk}^i`$ (note that $`\text{det}\stackrel{~}{g}=1`$ and $`\text{tr}\stackrel{~}{A}=0`$). See  for the explicit form of the evolution equations, and  for an analysis that indicates why the BSSN formulation should be superior to ADM at least for linearized perturbations of flat space.

In order to obtain the stable evolutions described below, we have found it necessary to add the following ingredients to the BSSN formulation:

1. As discussed in , we actively force the trace of the conformal-traceless extrinsic curvature $`\stackrel{~}{A}_{ij}`$ to remain zero during our simulations by subtracting it after each time step.
2. We use the independently evolved “conformal connection functions” $`\stackrel{~}{\mathrm{\Gamma }}^i`$ only in terms where derivatives of these functions appear. Whenever these functions are undifferentiated, we recompute them from the conformal Christoffel symbols. We have found this to be very important to achieve long-term stability, but at the moment we lack a theoretical understanding as to why this is so.

Slicing conditions— As a first approach to evolving the solution described above, one could think of using the exact value of the lapse. It turns out that it is difficult to keep the evolution stable if the lapse is not allowed to adapt to the (numerically induced) evolution of other dynamical quantities, particularly the trace of the extrinsic curvature. In order to obtain stable evolutions we have found it crucial to use a “live” slicing condition. What is required is a slicing condition that is well adapted to the exact solution in the sense that for this solution it recovers the exact lapse. For this we start from the Bona-Massó family of slicing conditions

$$_t\alpha =\alpha ^2f\left(\alpha \right)\mathrm{tr}K,$$
(3)
with $`f(\alpha )>0`$. As it is, this condition does not reproduce our exact solution for which $`\mathrm{tr}K0`$, but $`_t\alpha `$=0. However, one can easily see that for zero shift Eq. (3) implies $`_t\alpha _t(\text{det}g)`$. For this to hold also with non-zero shift Eq. (3) must be generalized to

$$_t\alpha =\alpha f\left(\alpha \right)\left[\alpha \mathrm{tr}K_i\beta ^i\right].$$
(4)
For any static solution Eq. (4) implies $`_t\alpha `$=0.

Another natural slicing condition to consider is $`_t\mathrm{tr}K`$=0. For initial data with $`\mathrm{tr}K`$=0 this condition leads to maximal slicing, but $`_t\mathrm{tr}K`$=0 is a gauge choice that can be made in general, even if $`\mathrm{tr}K0`$, as is the case for the constant time slices of the black hole in EF coordinates. This “K freezing” condition leads to an elliptic equation for the lapse,

$$\mathrm{\Delta }\alpha \alpha K_{ij}K^{ij}\beta ^i_i\mathrm{tr}K=0.$$
(5)
In the numerical implementation, we solve this equation for the lapse but we hold $`\text{tr}K`$ constant in time by hand. In  in the context of the evolution of strong waves we have found that otherwise a drift away from the initial value due to numerical errors can lead to an instability. Such drifts were one of the reasons that led us to consider trace-split formulations like BSSN, because here $`\text{tr}K`$ is evolved as an independent variable which makes it trivial to enforce $`_t\text{tr}K=0`$.

Shift conditions— In contrast to the experience with the lapse, we have found that using a static (exact) shift does allow us to get long-term stable evolutions. However, this is not useful in general, so we have considered also live shift conditions. Live shifts have been studied before for black hole spacetimes in , where a minimal distortion shift condition  led to limited stability ($`t100M`$) for a single excised black hole.

In our case a good choice was a conformal version of the 3-harmonic shift . 3-harmonic shifts play a natural role in mixed elliptic-hyperbolic systems . The condition we impose in the BSSN system is $`_t\stackrel{~}{\mathrm{\Gamma }}^k`$=0 (“Gamma freezing” condition, note that $`\stackrel{~}{\mathrm{\Gamma }}^k0`$), or

$`\stackrel{~}{\gamma }^{jk}_j_k\beta ^i+{\displaystyle \frac{1}{3}}\stackrel{~}{\gamma }^{ij}_j_k\beta ^k\stackrel{~}{\mathrm{\Gamma }}^j_j\beta ^i+{\displaystyle \frac{2}{3}}\stackrel{~}{\mathrm{\Gamma }}^i_j\beta ^j+\beta ^j_j\stackrel{~}{\mathrm{\Gamma }}^i`$ (6)
$`2\stackrel{~}{A}^{ij}_j\alpha 2\alpha \left({\displaystyle \frac{2}{3}}\stackrel{~}{\gamma }^{ij}_j\mathrm{tr}K6\stackrel{~}{A}^{ij}_j\varphi \stackrel{~}{\mathrm{\Gamma }}_{jk}^i\stackrel{~}{A}^{jk}\right)=0.`$ (7)
As mentioned before, $`_k\stackrel{~}{\mathrm{\Gamma }}^i`$ is computed from the independent variable $`\stackrel{~}{\mathrm{\Gamma }}^i`$, in other terms we use $`\stackrel{~}{\mathrm{\Gamma }}^i=\stackrel{~}{\gamma }^{jk}\stackrel{~}{\mathrm{\Gamma }}_{jk}^i`$ (notice that the momentum constraint was used to replace $`_j\stackrel{~}{A}^{ij}`$ in (LABEL:eq:Gammafreezing)). Equation (LABEL:eq:Gammafreezing) is an elliptic equation for the shift vector. For the solution of (5) and (LABEL:eq:Gammafreezing) we have used the multi-grid solver from BAM . As in the case of K freezing, we explicitly hold the value of $`\stackrel{~}{\mathrm{\Gamma }}^i`$ constant in time in order to prevent this quantity from drifting due to numerical errors. As shown in Section IV, allowing $`\stackrel{~}{\mathrm{\Gamma }}^i`$ to drift results in an unstable evolution.

We have also looked at shift prescriptions given by evolution equations instead of elliptic conditions. One way to do this is to transform an elliptic equation into a parabolic one by making $`_t\beta ^i`$ proportional to the given elliptic operator (“driver” conditions, see ). As an example we considered the following evolution equation for the shift obtained from the Gamma freezing condition (a “Gamma driver” condition)

$$_t\beta ^i=k_t\stackrel{~}{\mathrm{\Gamma }}^i,(k>0).$$
(9)
Boundary conditions— There are two very different boundaries to consider in our simulations: the outer boundary of the numerical grid, and the inner boundary of the excised region. In principle there should be a rigorous treatment of numerical boundaries at finite radii (starting e.g. from , the first analytic treatment of the initial boundary value problem). Here we are looking for simple numerical methods that are sufficient for the evolution of excised black holes.

At the outer boundary we have attempted to keep all fields equal to their exact values, but have found that this introduces late time instabilities. Using a live boundary condition allows us to eliminate these instabilities. The boundary condition we use is a radiative boundary condition applied to the difference between a given variable and its exact value: $`ff_{\mathrm{exact}}=u(rt)/r`$. We apply this condition to all fields (even to the lapse and shift in the case of the algebraic gauge conditions) except the $`\stackrel{~}{\mathrm{\Gamma }}^i`$ which we leave fixed to their exact values at the boundary. Applying this condition to the $`\stackrel{~}{\mathrm{\Gamma }}^i`$ causes a drift away from the exact solution that eventually crashes the simulation (the origin of this drift is not well understood, but it seems to be related to the shift choice and is not present if one uses the Gamma driver shift described above).

As to what boundary condition to use at the sides of the excision cube, we have experimented with many different conditions and have finally settled on one that simply copies the time derivative of every field at the boundary from its value one grid-point out along the normal direction to the cube (at edges and corners we define the normal direction as the diagonal). This condition is perfectly consistent with evolving a static solution, where the time derivatives are supposed to be zero. Even in a dynamical situation, this condition is still consistent with the evolution equations since it is equivalent to just calculating the source term one grid point away. This means that our boundary condition should introduce a first order error, but as mentioned above, we do not expect this error to affect the solution outside the horizon. One could in principle argue that nothing prevents gauge modes and constraint violating modes from propagating outside the horizon, thus spoiling the second order convergence of the exterior scheme. We have looked carefully at the convergence of our simulations, and have found no evidence that this happens in practice.

IV. NUMERICAL RESULTS

We now present some results of our numerical simulations. As discussed above, our simulations have been done with a live lapse condition, and we have considered both a static shift, and several live shift conditions. In our runs we have always taken $`M`$=1, so the horizon is a sphere of radius $`r`$=2, and we excise a cube of side 1 (we have in fact also excised cubes of different size, but the results discussed below are not affected by this). The numerical integration is carried out using the so-called iterative Crank-Nicholson scheme with 3 iterations (counting the initial Euler step as iteration 1). Because of the spherical symmetry of the problem typically only one octant was evolved (with positive $`x`$, $`y`$, and $`z`$). However, as discussed at the end of this section, an unstable mode appears when the same simulations are performed on the corresponding full grids.

Static shift— We first consider the case when the shift remains equal to its exact value. Figure 1 shows a log plot of the root mean square (r.m.s.) of the change in the lapse between consecutive time steps for two simulations using slicing condition (4) with $`f=1/\alpha `$ (“1+log” slicing ), a grid spacing $`\mathrm{\Delta }x=0.4`$, and a time step $`\mathrm{\Delta }t=0.1`$. Figure 1a shows the results of a simulation using $`53^3`$ grid points, with the outer boundaries at $`20M`$. The change in the lapse drops as an exponentially damped oscillation until at $`t3500M`$ it reaches the level of round-off error ($`10^{16}`$) and settles down (other functions show a similar behavior). The evolution was stopped at $`t=4000M`$, but it is clear that it could have continued. Figure 1b shows a simulation with the same resolution, but using $`103^3`$ grid points, with the outer boundaries now at $`40M`$. The simulation goes past $`t2000M`$, and seems to have settled on an exponentially decaying oscillating pattern. (This simulation took 100 hours running on 16 processors of an Origin 2000 SGI machine. If the pattern continues, round-off error level would be reached by $`t12000M`$, requiring another $`500\times 16`$ hours of computer time). The most obvious differences between the run with the boundaries at $`20M`$ and that with the boundaries at $`40M`$ is the fact that the period of the oscillations increases and the rate of decay decreases. The period increases by a factor of 3.4 as we double the distance to the outer boundaries, so the oscillation time scale is not given directly by the light travel time from the boundary (which would approximately double). We do not know exactly what fixes this time scale, but the fact that when we look at individual metric components we see that the oscillations behave like standing waves (and not travelling pulses) would seem to indicate that we are looking at different modes of oscillation of the whole system (interior plus boundaries).

These simulations are not only stable for very long times, they are also exceedingly accurate. We have located the apparent horizon every 50 time-steps (using the 3D finder described in ), measured its area $`A`$ and computed its mass $`M=\sqrt{A/(16\pi )}`$. Figure 2 shows the behavior of the horizon mass as a function of time. In both cases, after an initial transient, the mass settles on a stationary value with an error of less than $`1\%`$.

In Figure 3 we consider the convergence of our simulations by looking at the late time value of the Hamiltonian constraint along the $`x`$ axis for simulations with $`28^3`$, $`53^3`$ and $`103^3`$ grid points and resolutions of $`\mathrm{\Delta }x`$=0.4,0.2,0.1 respectively (boundaries at $`10M`$). The Hamiltonian constraint for the higher resolution runs has been multiplied by factors of 4 and 16. The fact that the three lines coincide indicates second order convergence.

Elliptic shifts— We now consider results with elliptic shifts, such as those that we expect will be needed in a 3D black hole merger simulation. Figure 4 shows two stable and three unstable runs up to $`t=400M`$, and Figure 5 shows those three runs that lasted longer up to $`t=3000M`$. Second order convergence has been checked using two grids with $`19^3`$ and $`35^3`$ points with the outer boundary at $`7M`$. For 1+log slicing a radiative boundary condition is applied to the lapse, while lapse and shift for the elliptic conditions are held fixed at the exact values.

Stable runs are obtained for Gamma freezing shift with either 1+log or K freezing slicings. Referring to Figures 4 and 5, for 1+log slicing $`\mathrm{\Delta }\alpha _{\mathrm{rms}}`$ falls below $`10^{16}`$ at $`t1500M`$ after four oscillations (run 1), while for K freezing there are more than fifteen oscillations, which damp out at around $`10^{10}`$ followed by a straight line decay (run 2).

The 1+log, Gamma freezing run becomes unstable if the boundary values of all fields are static (run 3, crashing at $`t1500M`$, Figure 4b), or if $`_t\stackrel{~}{\mathrm{\Gamma }}^i=0`$ is not set to zero identically and is allowed to drift because of numerical errors (run 4, crashing at $`375M`$). We also tested 1+log slicing with a minimal distortion shift  computed from the ADM variables, but this run fails already at $`27M`$ (run 5).

Algebraic shifts— Finally, we consider a simulation using 1+log slicing and a Gamma driver shift with $`k=0.1`$. Figure 6 shows the r.m.s. of the change in the lapse and the horizon mass for a simulation with $`\mathrm{\Delta }x=0.4`$, $`\mathrm{\Delta }t=0.1`$ and $`53^3`$ grid points. After $`t2500M`$ the solution becomes static up to round-off error.

Discussion— The above results demonstrate that stable 3D black hole runs can be obtained with the simple excision technique that we introduced in this paper, with a variety of different gauge conditions. However, repeating these runs on a full grid as opposed to just one octant, with otherwise identical parameters, uncovers an unstable mode. Figure 7 shows as an example the situation for 1+log slicing and static shift, although the problem appears for all the gauge conditions considered here. Tracing the growth of the unstable mode back in time suggests that it has started as numerical round-off error of around $`10^{14}`$ at $`t=0`$. Increasing the grid resolution appears to have no significant effect on the growth rate of the unstable mode, but the simulation now crashes slightly sooner. However, we do see good second order convergence at early times, before the instability becomes apparent. The situation does not improve if we impose the exact data at the excision boundary (imposing exact data at the excision boundary in octant mode works well and leads to stable simulations). Also, the presence of a horizon does not seem to be the cause of the problem since when we excise a cube that contains the horizon, as opposed to being contained by it, the instability is still present although it becomes somewhat milder (not surprising since we have excised a region with stronger data). While the achievable run times of about $`500M`$ are roughly 10 times larger than for singularity avoiding slicings, we have found that introducing an artificial asymmetry on the full grid by simply off-setting the excision box one grid point in all directions makes the runs fail much sooner. Although the slope of the blow-up is not significantly affected when this artificial asymmetry is introduced, the exponential growth becomes evident from the very beginning. On the other hand, the full grid runs can be stabilized by setting certain terms in the BSSN equations to their analytic values. In particular, freezing the evolution of the $`\stackrel{~}{\mathrm{\Gamma }}^i`$ while keeping the shift static suffices to obtain stability. In conclusion, the instability appears to be more directly linked to the system of evolution equations than to the boundary condition, and we will investigate different variations of the evolution system in the future.

We have also repeated the above simulations using the ADM equations with the same gauge and boundary conditions, and the same numerical techniques, but these runs fail already at $`t30M`$ even in octant mode.

V. CONCLUSIONS

We have described a black hole excision technique in 3+1 numerical relativity that has allowed us to obtain accurate, long-term stable evolutions of a black hole spacetime in 3D. The main limitation is that the transition from octant symmetry to full grids introduces an unstable mode, which is currently under investigation. Our implementation of excision is based on the idea of simplifying all ingredients of the excision algorithm as much as possible. In our case this means: 1) excising a cube naturally adapted to the underlying Cartesian coordinates, 2) imposing a simple but stable boundary condition on the sides of this cube, and 3) using an upwind scheme instead of causal differencing. Crucial for obtaining our long-term stable evolutions has been the use of a live slicing condition and a radiation outer boundary condition. Although keeping a static shift does not appear to have a detrimental effect on the stability of our simulations, we have also experimented with several live shift conditions, both algebraic and elliptic, that can be generalized to more interesting physical situations. We consider these results a necessary first step towards the development of excision techniques capable of evolving the full inspiral collision of two black holes in an accurate and stable way.

###### Acknowledgements.

We would like to thank J. Baker, D. Pollney, E. Seidel, W.-M. Suen and J. Thornburg for useful discussions and comments. The numerical experiments were implemented using BAM and the Cactus Computational Toolkit . All computations were performed at the Max-Planck-Institut für Gravitationsphysik.