text:
# Generic Tracking of Multiple Apparent Horizons with Level Flow

## I Introduction

Our goal is to investigate the strong field regime of general relativity. In particular we wish to focus on the study of coalescing black hole binaries. Over the last three decades since the pioneering work of Cadez, Smarr, Eppley and others, advances in computing technology, numerical algorithms and techniques and our understanding of the underlying physics have advanced to a point where we are able to carry out simulations of binary black hole collisions in 3+1 dimensions. One of the outcomes of such simulations will be an understanding of the underlying physics of the problem; and, therefore, a prediction and understanding of the gravitational radiation content. A detailed knowledge of how the resultant gravitational waveforms relate to the physical parameters of the binaries that produce them will be of importance to gravitational wave observatories (such as LIGO, VIRGO, TAMA, GEO600) now under construction around the world. With the building of these observatories, we stand at the epoch of the first direct observations of astrophysical sources that involve strong field general relativity.

The orbit and merger of two black holes is one candidate source for ground based detection of gravitational waveforms. This is of great interest to the relativity community. The binary black hole problem is a two-body problem in general relativity. It is a stringent dynamical test of the theory. However, studying spacetimes containing multiple black holes involves solving the Einstein equations, a complex system of non-linear, dynamic, elliptic-hyperbolic equations intractable in closed form. The intractability of the problem has led to the development of numerical codes capable of solving the Einstein equations.

Our approach to numerically solving the Einstein equations involves reformulating them as an initial value problem. In this 3+1 formulation , spacelike hypersurfaces parameterized in time foliate the spacetime. The resulting equations are coupled elliptic and hyperbolic differential equations of the 3-metric, $`g_{ij}`$, and the extrinsic curvature, $`K_{ij}`$. The initial value problem is solved by specifying a hypersurface at an instant of time, say $`t=0`$, and evolving to the next hypersurface at time $`t=\delta t`$ with the evolution equations to obtain $`g_{ij}`$ and $`K_{ij}`$ at the next time $`t=\delta t`$.

One vital issue in numerically solving the Einstein equation describing spacetimes containing black holes is handling the physical singularities. As one approaches the singularity, the values of the fields being computed approach infinity; therefore, a region containing the singularity must be avoided to keep the computation from halting. Excision techniques are promising in avoiding the singularity. Excising the singularity involves locating a region interior to the event horizon containing the singularity on each evolving hypersurface. This region is then “masked” from the computation. The derivatives at the boundary between the masked region and the computational domain are handled using causal differencing, a differencing scheme that respects the causal structure of the spacetime.

In deciding where in space to excise we use the apparent horizon as opposed to the event horizon. By its very nature, the event horizon is a global construct depending on the entire spacetime. The apparent horizon, a local, i.e. spacelike 2-surface is more suitable. Following a suggestion by Unruh , the apparent horizon is used to define the excised region to be masked at each time during the evolution. Apparent horizons are defined locally for each time, and exist at, or inside of, the event horizon. In some spacetimes, choices of foliation may lead to the absence of an apparent horizon. When the discussions in this refer to a hypersurface, it is assumed an apparent horizon exists on that hypersurface.

Recent three-dimensional horizon locator codes are capable of finding the location of an apparent horizon in generic single black hole spacetimes  and two  are capable of finding disjoint multiple apparent horizons in the special case of conformally flat binary, time-symmetric black hole spacetimes. Multiple apparent horizon finding algorithms will be necessary in simulations of generic binary black hole spacetimes. The method presented in this paper, called the level flow method, is capable of detecting multiple apparent horizons in generic spacetimes. The level flow algorithm has two advantages: 1) It is independent of a good initial guess and 2) It is capable of following the surface through a change in topology. In level flow, the apparent horizon equation is reformulated as a parabolic equation and a set of surfaces are flowed with speeds dependent on the expansion of the outgoing null vectors normal to each surface.

The purpose behind developing the level flow method of tracking apparent horizons is to have a method capable of detecting multiple apparent horizons on any given hypersurface without a good guess. Specifically, we want a tracker that can detect the transition from a double to a single apparent horizon in single time step without a prior knowledge of the transition.

In the rest of this paper we discuss apparent horizons in general and current 3D work in $`\mathrm{\S }\text{II}`$ and $`\mathrm{\S }\text{III}`$. In $`\mathrm{\S }\text{IV}`$ we describe the level flow method in detail and give a brief description of the numerical method involved in solving the apparent horizon is in $`\mathrm{\S }\text{V}`$. We demonstrate the use of the level flow tracker on model data and a binary black hole grazing collision in $`\mathrm{\S }\text{VI}`$ and $`\mathrm{\S }\text{VII}`$.

## II Apparent Horizons

$`M`$ is the spacetime with metric $`{}_{}{}^{4}g_{ab}^{}`$ foliated by hypersurfaces $`\mathrm{\Sigma }`$ parameterized by $`t`$ with 3-metric $`g_{ab}`$. Let $`S`$ be a surface with $`S^2`$ topology on $`\mathrm{\Sigma }`$. The apparent horizon is the outermost marginally trapped surface in $`\mathrm{\Sigma }`$, a surface with zero expansion. The zero expansion of the surface is defined in terms of outgoing null vectors to $`S`$, $`k^a`$, such that $`k^a`$ have zero divergence

$$_ak^a=0,$$
(1)
where $`_a`$ is the covariant derivative associated with $`{}_{}{}^{4}g_{ab}^{}`$. Referring to fig. (1), $`k^a`$ is defined in terms of the spacelike normals to $`S`$, $`s^a`$, and the future directed timelike normals, $`n^a`$, such that:

$$k^a=s^a+n^a.$$
(2)
The expansion of the outgoing normals, $`_ak^a`$, can be rewritten in terms of quantities defined on the hypersurface:

$$\kappa D^as_aK+s^as^bK_{ab}.$$
(3)
$`D_a`$ denotes covariant differentiation with respect to $`g_{ab}`$, $`K_{ab}`$ is the extrinsic curvature, and $`K`$ is $`g^{ab}K_{ab}`$. In fact, there is a level set of surfaces in $`\mathrm{\Sigma }`$ parameterized by $`\kappa `$. Each surface in the level set is defined by the constant expansion of its null vectors such that

$$\kappa =c_n,$$
(4)
where $`c_n`$ are constants labeled by the positive integer $`n`$. Marginally trapped surfaces are members of this set for

$$\kappa =0.$$
(5)
Eqn.(5) is called the apparent horizon equation since the apparent horizon is the outermost surface with $`\kappa =0`$ in $`\mathrm{\Sigma }`$.

The $`S^2`$ topology of the apparent horizon naturally lends itself to characterization via spherical coordinates. The function,

$$\psi =rh(\theta ,\varphi )$$
(6)
is a level set of 2-spheres in $`\mathrm{\Sigma }`$, where $`h(\theta ,\varphi )`$ is called the apparent horizon shape function. A marginally trapped surface has $`\psi =0`$. The spacelike normals to $`S`$ are defined from eqn. (6) such that

$$s^i=g^{ij}_j\psi /\sqrt{g^{kl}_k\psi _l\psi }$$
(7)
is the spacelike vector field at every point of $`S`$. The apparent horizon equation in spherical coordinates ($`h(\theta ,\varphi ),\theta ,\varphi `$) is a 2-dimensional problem in $`\theta `$ and $`\varphi `$.

## III Current 3-D Methods

The approaches to solving the apparent horizon equation on a three-dimensional hypersurface can be addressed roughly in two categories: methods that solve the apparent horizon equation directly and methods that solve it by first recasting it as a parabolic equation. This paper does not address spherical and axi-symmetric approaches. One of the first three-dimensional apparent horizon solvers was published by Nakamura, Kojima, and Oohara . They directly solve the apparent horizon equation by using spherical harmonic basis functions to expand the apparent horizon shape function, $`h(\theta ,\varphi )`$ into

$$h(\theta ,\varphi )=\underset{l=0}{\overset{l_{max}}{}}\underset{m=l}{\overset{l}{}}a_{lm}Y_{lm}(\theta ,\varphi ).$$
(8)
This method is called the pseudospectral method. A finite number of the coefficients, $`\{a_{lm}\}`$ parameterize the horizon shape function, and the maximum $`l_{max}`$ depends on the computation and deviations from a sphere. The apparent horizon equation can then be solved by writing it as

$$\kappa (a_{lm})=0,$$
(9)
and using functional integration routines to find the coefficients $`a_{lm}`$. Others have used similar methods .

In another approach to direct solutions of the apparent horizon equation is to treat it as a boundary value problem. One notes that a discretization of this equation leads to a system of algebraic equations which can then be solved via Newton’s method. Thornburg  discusses applications of Newton’s method to this problem in general and shows results in axisymmetry. Huq  has implemented a similar algorithm based on Newton’s method that utilizes Cartesian coordinates to difference the equations.

Tod  first suggested the use of curvature flow in solving the apparent horizon equation by recasting it as a parabolic equation. Bernstein  implemented Tod’s suggestion in axisymmetry. Gundlach  introduced a fast flow method which combines the ideas of the flow method with the pseudospectral method. Pasch  and Diener  implemented a similar method, a level-set method, in three-dimensions and found discrete apparent horizons in multi-black-hole spacetimes; however these spacetimes were confined to be conformally flat and time-symmetric.

Each of the approaches briefly described above, solving the apparent horizon equation directly or solving it via a parabolic equation, has its advantages. Direct methods tend to be faster while flow methods do not rely on “good” initial guesses. However, none of these methods are applicable to the generic, multi-black-hole problem. Herein lies the motivation behind the level flow method. The level flow method is the only method designed to locate discrete apparent horizons in generic spacetimes containing one or two discrete horizons.

## IV Level Flow Method

### A Curvature Flow

The level flow method is a hybrid flow$`/`$level$``$set method. The previous section mentioned the flow method, this section gives more detail on the flow method which is the foundation of the level flow method. The flow approach, as suggested by Tod, is to rewrite the apparent horizon equation as the speed function in a parabolic equation. In the case of a time symmetric hypersurface, in which $`K_{ab}=0`$, the apparent horizon equation reduces to the condition for a minimal surface, $`D_as^a=0`$. In this case, the surface $`S`$ is at a local extremum of the area. The starting surface, $`S(\lambda =0)`$, is parameterized by coordinates $`x^a`$ and evolved in terms of a parameter $`\lambda `$. The equation of motion is

$$\frac{x^a}{\lambda }=Hs^a$$
(10)
where $`x^a/\lambda `$ is a vector field, and $`H`$ is the mean curvature, which is the trace of extrinsic curvature associated with embedding $`S`$ in $`\mathrm{\Sigma }`$ given by

$$H=D_as^a.$$
(11)
Eqn.(10) is the gradient flow for the area functional. The area decreases monotonically with increasing $`\lambda `$. Grayson  has shown that a surface deforming under its gradient field (Eqn.(10)) will evolve to a stable minimum surface (surface is local minimum of area) if there is one, otherwise to a point.

In numerical relativity, we are interested in the generic case, with $`K_{ab}0`$, for which the marginally trapped surfaces differ from minimal surfaces, the surfaces are not extrema of the area. However, Tod suggests an equation similar to Eqn.(10) as a curvature flow:

$$\frac{x^a}{\lambda }=F(\kappa )s^a$$
(12)
using $`F(\kappa )=\kappa `$ where $`\kappa =D_as^a+s^as^bK_{ab}K`$ as in eqn.(3). We have found eqn. (12) to be a successful practical implementation of the flow method.

### B Level Flow

Eqn. (12) gives us an initial value problem. Given information about the system at some initial $`\lambda `$, eqn. (12) will describe the system for all its future propagation in $`\lambda `$. Directly solving eqn. (12) will lead to a successful detection of single apparent horizons; however, solving it directly does not ensure correct handling of a topology change which is necessary for detection of multiple apparent horizons. By combining the flow method with a level-set idea however, this topology change can be effected and multiple apparent horizons can be tracked starting from a single guess surface.

First eqn. (12) is recast from an equation governing the motion of the coordinates parameterizing $`S`$, namely $`x^a`$, to an equation governing the motion of the surface $`\psi `$. Noting that $`\psi `$ is parameterized by $`\lambda `$,

$`{\displaystyle \frac{\psi }{\lambda }}={\displaystyle \frac{x^a}{\lambda }}{\displaystyle \frac{\psi }{x^a}}`$ (13)
by the chain rule, and multiplying eqn.(12) by $`\frac{\psi }{x^a}`$ gives

$$\frac{\psi }{\lambda }=F(\kappa )s^a\psi /x^a.$$
(14)
Using

$$s^a=g^{ab}\frac{\psi }{x^b}/\psi $$
(15)
and

$$\psi =\sqrt{\frac{\psi }{x^a}\frac{\psi }{x^b}g^{ab}},$$
(16)
the test surface’s flow is given by:

$$\frac{\psi }{\lambda }=F(\kappa )\psi .$$
(17)
Eqn.(17) is a reformulation of eqn. (12) and will flow the surface, $`\psi `$, to a marginally trapped surface at $`\psi =0`$ when $`\kappa =0`$.

The strength of the flow methods is their ability to locate a surface with $`\kappa =0`$ given any non-pathological initial surface. For example, the apparent horizon in a spherically symmetric spacetime (Schwarzschild) was found by flowing an initial surface shaped as a leaf, see Fig. (2). This ability is especially important when tracking horizons during evolutions of binary black hole spacetimes. In this case, finding apparent horizons for two discrete apparent horizons in each $`\mathrm{\Sigma }(t)`$ involves flowing two initial guesses, one for each horizon . On $`\mathrm{\Sigma }(t=0)`$, the location of the apparent horizons may be known; however, as the black holes accelerate the task of guessing the locations of the two horizons becomes more difficult. Further, the two horizons merge into a single horizon at a single instant of time, rendering a good initial guess difficult. Some way of determining when two horizons merge into a single horizon is necessary. The level flow method takes care of these issues by not requiring a good initial guess $`(\psi (\lambda =0))`$ and by detecting multiple apparent horizons from a single guess $`(\psi (\lambda =0))`$.

The level flow method differs from the flow method in the specification of the speed function, $`F(\kappa )`$. $`F(\kappa )=0`$ determines when the propagation of the surface stops. The flow method is in the form of eqn. (12), in which $`F(\kappa )=\kappa `$. A good choice since $`F(\kappa )=\kappa =0`$ indicates a marginally trapped surface; but this choice will not flow $`\psi `$ though a fission. In general the scheme fails as the surface pinches due to ill-defined normals at the surface. The level flow method alleviates this problem by looking for indications that the surface topology is about to change before the pinching occurs. (Another method which was introduced by Sethian and Osher  for non-relativistic problems is to flow a higher dimensional surface in which $`\psi `$ is embedded. This higher dimensional surface does not fission. This has only been implemented in a time-symmetric spacetime  and requires more computational power due to the extra dimension in the problem.)

The level flow method flows a set of two-dimensional surfaces in $`\mathrm{\Sigma }`$ parameterized by $`\kappa `$. We call the set of surfaces a level set and label the set $`S(c_n)`$. Each surface has a constant value of $`\kappa =c_n`$ everywhere on it. The set of surfaces is defined by varying $`c_n`$ as the flow progresses

$$c_{n+1}=c_n\pm \mathrm{\Delta }c,$$
(18)
where $`(+)`$ indicates outward flow, $`()`$ inward flow, and $`\mathrm{\Delta }c\kappa _2`$. Each surface obeys the equation of motion given in eqn. (17) with $`F(\kappa )`$ defined to flow to multiple surfaces. We choose two options for the speed function:

$`F(\kappa )`$ $`=`$ $`\kappa c_n`$ (19)
$`F(\kappa )`$ $`=`$ $`(\kappa c_n)\mathrm{arctan}^2({\displaystyle \frac{\kappa c_n}{\kappa _o}}).`$ (20)
As $`\kappa c_n0`$, both functions are solving for a particular surface in the level set, $`S(c_n)`$. The second function, eqn.(20), behaves similarly to the first but allows for larger time steps near a fissioning surface because it moves points further from the $`\kappa c_n=0`$ surface faster than the points closer to this surface.

Eqn. (17) is an initial value problem requiring $`\psi `$ to be specified at $`\lambda =0`$. To initialize the starting surface, we need only supply an origin and radius. Taking into account that there may be more than one marginally trapped surface, it is best to start with an initial surface larger than the expected horizon. The values of $`g_{ij}`$ and $`K_{ij}`$ are required everywhere on the surface to evaluate $`\kappa `$. These functions can be known analytically or generated by evolution codes. As the flow velocity approaches zero, $`F(\kappa c)0`$, $`\kappa c_n`$ and a surface $`S(c_n)`$ is found within a tolerance ($`ϵ_\kappa `$). When $`\kappa =0`$, the located surface describes a marginally trapped surface.

Fig.(4) shows the level set found in a spacetime containing two black holes with coordinate locations $`(0.954,0,0.3)M`$ and $`(0.954,0,0.3)M`$. Each 2-surface has a constant value of $`\kappa `$. We monitor the topology of the deforming surface by computing the radial component of the gradient of $`\kappa `$ with respect to the normals of each surface in the level set. The gradient is defined as:

$$\frac{\kappa _{n1}\kappa _n}{\psi _{n1}\psi _n},$$
(21)
where $`\psi `$ is the function given in eqn.(6). A sharp increase in the gradient indicates the existence of multiple surfaces. To ensure that we do not erroneously abandon a single surface, we also monitor the maximum of the $`l_2`$-norm of $`\kappa `$. If $`\kappa `$ is no longer decreasing, we are no longer finding a solution to eqn.(4); otherwise the single surface is retained. The level flow method is essentially a special set of surfaces with properties that let us determine when to break. If we only flowed to $`\kappa =0`$, we would not form the collection of $`\kappa =`$constant surfaces.

Once a topology change is indicated, the radii and origins for each of the new surfaces are found (note that these four parameters for each surface are all that is needed to initiate two new trial surfaces). These origins and radii are determined using the location of the last of the single surfaces. Using this last single surface, we can find the origin of the last surface and the location on the surface with minimum gradient of the value of $`\kappa `$. This occurs at the farthest points from the pinching in the surface. Picture (a) in fig. (5) shows the last single surface with an arrow drawn from the origin to the point on the surface with a minimal change in $`\kappa `$. The arrow indicates a chosen direction. All points lying in this direction are collected and averaged to find a radius and center of mass. All points lying in the opposite direction are also collected and used to calculate the radius and center of mass for the second surface using the dotted arrows in picture (b) of fig.(5). These two sets of radii and centers are the initial starting parameters for each new trial surfaces. The tracker then flows the two new surfaces depicted in (c) of fig.(5) until $`\kappa =0`$ within $`ϵ_\kappa `$.

The level flow code is only started by the user once, the subsequent flowing to multiple apparent horizons is done automatically.

The advantage to the level flow method is its capability to detect apparent horizons in generic, multiple black hole spacetimes from a single reasonable initial guess. The drawbacks are the dependence of $`\mathrm{\Delta }\lambda `$ on the spatial grid size, $`\mathrm{\Delta }\lambda N^2`$ where $`N^2=N_\theta N_\varphi `$ is the number of grid points, and the fact that we flow to a speed of zero (the flow speed approaches zero as $`\kappa `$ approaches zero). When using apparent horizon tracking in our evolution code, we will not require knowledge of the apparent horizon location to high precision; in fact we can find a surface with $`\kappa 0`$ to remove the singularity thus alleviating some of the speed issues. Nonetheless, we plan to improve the speed of this algorithm. Some improvements have already been made to increase the efficiency of the current algorithm. The addition of the $`\mathrm{arctan}^2`$ function, eqn.(20), speeds up the algorithm during the fissioning process. Monitoring the number of steps needed to complete a Crank-Nicholson iteration (see $`\mathrm{\S }`$V)) has also proven useful in increasing efficiency.

## V Numerical Method

The previous section described how the level flow method is used to solve the apparent horizon equation. The resulting parabolic equation is updated using an iterative Crank-Nicholson method updating the variables at every $`\lambda `$-step. Iterative Crank-Nicholson converges to an exact solution of the implicit problem. However, the detailed behavior of this convergence  shows that the Crank-Nicholson solution at a particular iteration has an amplification factor $`|𝒜^{(n)}|`$ that oscillates around unity. The behavior varies in pairs: $`|𝒜^{(n)}|<1`$ for $`n=2,3`$; $`|𝒜^{(n)}|>1`$ for $`n=4,5`$, etc. while $`||𝒜^{(n)}|1|0`$ monotonically as $`n\mathrm{}`$. $`n`$ is counting the number of iterations it takes to get $`\widehat{\kappa }=\kappa `$ within the specified tolerance. For the data presented here, a Crank-Nicholson iteration of $`n=2`$ or $`n=3`$ was maintained for errors less than the spatial grid spacing squared, $`h^2`$.

The spatial derivatives are approximated to second order in truncation error using centered finite differencing molecules. To verify the convergence of the level flow code, we include a plot of the convergence factor versus the number of iterations taken in fig. (6).

The convergence factor is given by

$$C_f\frac{\widehat{\kappa }_{2h}\widehat{\kappa }_{4h}}{\widehat{\kappa }_h\widehat{\kappa }_{2h}},$$
(22)
where $`\widehat{\kappa }`$ is the discretized $`\kappa `$ and $`h`$ is the spatial grid spacing. For a second order scheme, the convergence factor in eqn. (22) is $`C_f=4+O(h^2)`$.

For the closed form solutions detailed in the next section, the data are given by evaluating the closed form analytically on the two-surface. However, the goal is to use the level flow method during an evolution including evolutions involving a region excised from computational consideration. The approach we take to evaluate $`\kappa `$ from a Cartesian grid of data $`(g_{ij},K_{ij})`$ is the same as that used and described in Huq . This approach discretizes the apparent horizon equation using Cartesian coordinates on 3d-stencils centered on points on the surface. These stencils are not aligned with the 3d-lattice from which we obtain $`g_{ij}`$ and $`K_{ij}`$ data. Our apparent horizon surfaces are embedded in such lattices and as a result interpolations must be carried out to obtain the metric data on the surface as it evolves. The algorithms and methodology for evaluating $`\kappa `$ are described in detail in .

## VI Testing the Method with Closed Form Solutions

The level flow method of tracking apparent horizons has been designed to locate apparent horizons in single and multiple black hole spacetimes. To test the level flow tracker, we locate apparent horizons in Schwarzschild, Kerr, and Brill-Lindquist data. In particular, we also demonstrate the level flow method’s ability to detect binary black holes in the Brill-Lindquist data.

### A Schwarzschild Data

The Kerr-Schild metric provides a closed-form description of both the Schwarzschild and the Kerr solutions to the Einstein equation and is given by:

$$g_{ab}=\eta _{ab}+2Hl_al_b,$$
(23)
where $`\eta _{ab}`$ is the Minkowski metric, $`\eta _{ab}=`$diag$`(1,1,1,1)`$. $`H`$ is a scalar function of the coordinates and $`l_a`$ is an ingoing null vector with respect to both the Minkowski and full metrics; that is $`l_a`$ satisfies the relation:

$$\eta ^{ab}l_al_b=g^{ab}l_al_b=0.$$
(24)
For the ingoing Eddington-Finkelstein form of the Schwarzschild solution, the metric given in eqn.(23) has the scalar function, $`H`$, given by:

$$H=\frac{M}{r}$$
(25)
and the components of the null vector:

$`l_t`$ $`=`$ $`1`$ (26)
$`l_x`$ $`=`$ $`{\displaystyle \frac{x}{r}}`$ (27)
$`l_y`$ $`=`$ $`{\displaystyle \frac{y}{r}}`$ (28)
$`l_z`$ $`=`$ $`{\displaystyle \frac{z}{r}}`$ (29)
where we have adopted rectangular coordinates ($`t,x,y,z`$) with $`r=\sqrt{x^2+y^2+z^2}`$, and $`M`$ the mass of the black hole.

We track the apparent horizon in this situation for a single black hole of mass $`M`$. The area and radius of the event horizon for the Schwarzschild solution of the Kerr-Schild metric is known in closed form  given by:

$$A=4\pi r_+^2$$
(30)
where $`r_+`$ is the event horizon radius and equals $`2M`$. In the slicing we have chosen, the apparent horizon coincides with the event horizon. Using the level flow method we found the apparent horizon to converge to the closed form solution giving a $`0.35\%`$ relative error at a course resolution ($`17\times 17`$ grid). The area of the tracked apparent horizon is computed by

$$A_{num}_S\sqrt{h}𝑑x𝑑y,$$
(31)
and converges to the closed form solution, eqn.(30). In eqn.(31) $`h`$ is the determinant of the 2-metric $`h_{ab}`$ on the apparent horizon surface, and $`x`$ and $`y`$ are surface coordinates. The numerical area is determined from eqn.(31) by calculating the determinant at every point in the grid and using a trapezoidal integration scheme .

### B Kerr Data

The Kerr solution is a second solution given by the Kerr-Schild metric, eqn.(23). The Kerr solution is the solution for a spinning black hole, i.e. a black hole with an internal angular momentum per unit mass given by $`a`$. In rectangular coordinates $`(t,x,y,z)`$, the scalar function and null vector are given by:

$$H=\frac{Mr^3}{r^2+a^2z^2}$$
(32)
and

$$l_\mu =(1,\frac{rx+ay}{r^2+a^2},\frac{ryax}{r^2+a^2},\frac{z}{r}),$$
(33)
where $`\mu =(t,x,y,z)`$, $`M`$ is the mass of the black hole, $`a=J/M`$ is the angular momentum per unit mass of the black hole in the z-direction, and r is obtained from:

$$\frac{x^2+y^2}{r^2+a^2}+\frac{z^2}{r^2}=1:$$
(34)
$$r^2=\frac{1}{2}(\rho ^2a^2)+\sqrt{\frac{1}{4}(\rho ^2a^2)^2+a^2z^2},$$
(35)
with $`\rho =\sqrt{x^2+y^2+z^2}`$. The difference here is the addition of angular momentum. We test two cases, $`a=0.5M`$ and $`a=0.9M`$. Fig.(7) presents a Schwarzschild ($`a=0M`$) case, together with the $`a=0.5M`$ and $`a=0.9M`$ cases. The solid line is the $`\theta =\pi /2`$ slice and the dashed line is the $`\varphi =\pi `$ slice. We find the expected result, that the deformation in the $`\varphi =\pi `$ slice increases with $`a`$.

The radius of the event horizon is given by

$$r_+=M+\sqrt{M^2a^2}.$$
(36)
The solution to eqn.(36) for $`a=0.5M`$ is $`r_+=1.87M`$ and the numerical solution we obtain for the horizon radius is $`r_{num}=1.87M`$, with an error of $`0.06\%`$. In the $`a=0.9M`$ case, $`r_+=1.44M`$ and $`r_{num}=1.46M`$, with a $`1.39\%`$ error. The area of the horizon for each case can be calculated using

$$A=4\pi (r_+^2+a^2)$$
(37)
(generalizing eqn.(30)), and compared to a numerical eqn.(31). computation using eqn.(31). In the $`a=0.5M`$ case, eqn.(30) gives $`A=46.89M^2`$, numerically we obtain $`A_{num}=46.88M^2`$, resulting in a $`0.21\%`$ error. In the $`a=0.9M`$ case, $`A=36.09M^2`$ and $`A_{num}=36.39M^2`$ with a $`0.83\%`$ error. The errors will decrease as $`\kappa `$ is driven closer to zero.

### C Brill-Lindquist Data

In this section, we study a binary black hole system using Brill-Lindquist data \[IX\]. These data are useful to us for two reasons: We can verify previous results of the critical separation, and study an example of how the tracker works in finding multiple apparent horizons. The 3-metric is time symmetric, $`K_{ab}=0`$, and is conformally flat:

$$g_{ab}=\varphi ^4\eta _{ab}$$
(38)
where

$$\varphi =1+\underset{i=1}{\overset{N}{}}\frac{M_i}{2r_i}$$
(39)
and $`N`$ is the number of holes (here $`N=2`$), $`M_i`$ is the mass of the ith black hole and the $`r_i`$ are the radial distances from the centers of the black holes

We use isotropic coordinates to express the metric as

$$ds^2=\varphi ^4(dr^2+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2)$$
(40)
with

$$r_i=\sqrt{r^2+d_i^22d_ir\mathrm{cos}\theta },$$
(41)
where $`d_i`$ is the distance between the holes and the center of the coordinate system. When they are far apart, each hole has an apparent horizon of radius $`M/2`$ in these coordinates. The area of each of the holes when they are well separated is $`16\pi M^2`$.

The limiting separation of the holes between single and double horizons was found by Brill and Lindquist  to be $`1.56M`$, Cadez $`1.534M\pm 0.002M`$ , $`[1.5M,1.6M]`$ by Alcubierre et al. , and $`1.535M`$ by Huq . We found a critical separation 1.53(5)M. The apparent horizon at the critical separation of 1.535M is shown in fig.(8) using the level flow code with $`33^2`$ grid points.

The horizon found for a separation of $`d=1.5M`$ which is less than the critical separation, is shown in fig.(9). Fig.(10) is a plot of the $`l_2`$-norm of the maximum of $`\kappa (\theta )`$ for the separation $`d=1.5M`$ at each iteration plotted versus the number of $`\lambda `$-steps. This is one of the checks in the level flow code to ensure that the apparent horizon equation is still being solved. We expect the expansion to continue to decrease if we have started outside the apparent horizon and are flowing inward. As we will see in fig.(14), the expansion increases as fission occurs in a data set with separated holes.

As we increase the separation between the two holes to a separation greater than the critical separation, we can test the apparent horizon tracker in the case of multiple apparent horizons. We demonstrate with a separation of $`d=2.0M`$. The initial surface flows to the point of fissioning where the topology of the surface changes from a one surface into two surfaces. Fig.(11) is a plot of the initial surface that begins the flow. The level set found during this flow is depicted in fig.(12). Each of the surfaces in fig.(12) has a constant expansion, $`\kappa =c_n`$ and was used to indicate a topology change in the test surface. The values for the expansion are $`c_1=0.14`$, $`c_2=0.12`$, $`c_3=0.1`$, and $`c_4=0.08`$. The last single surface just before the topology change is not a surface in the level set; it is plotted in fig.(13). At this point the tracker begins to flow two surfaces.

In contrast to a separation of $`d=1.5M`$ where there is no fission, here as fissioning becomes imminent, the $`\kappa `$ begins to increase. Fig.(14) is a plot of the absolute value of the maximum across the surface of the expansion, $`\kappa `$, versus $`\lambda `$ up to the point of fission. The increase in the expansion is one of the signals of imminent fission. As the algorithm tries to find a surface with $`\kappa =0`$ everywhere, it is driven into two surfaces. Once the new surfaces are found, the maximum of the expansion begins a monotonic decrease as in fig.(10).

The exaggerated peanut shape in fig.(15) and fig.(16) is taken for the same $`\lambda `$-value as fig.(13).

Once the fissioning is detected by the code, it automatically begins flowing two new surfaces of the same resolution as the parent surface. The series of snapshots shown in fig.(15) and fig.(16) is a subset of the set of surfaces found by the apparent horizon tracker as it follows the fission of the trial surface into two surfaces. The tracker starts with a spherical starting surface that deforms along the gradient field.

As the points defining the surface flow, the distance between the points can become too small for the finite difference scheme at that resolution. Redistribution of the points on the surface is taken care of automatically by updating the center and radius.

## VII Apparent Horizons in a Grazing Collision

As stated in the Introduction, one of the main motivators of this work is to have an apparent horizon finder that can locate disjoint horizons during the evolution. This entails 1) finding the horizons without a good initial guess, and 2) detecting the topology change from two disjoint horizons to one horizon. To demonstrate the level flow’s ability to carry out 1) and 2), we report the results of apparent horizon tracking in the particular case of a short evolution of two spinning, Kerr-Schild black holes using the Binary Black Hole Grand Challenge Alliance Cauchy code . A future paper  will discuss the details of the evolution.

The evolution is free, i.e. the momentum and Hamiltonian constraint equations are only used as checks during the evolution. Since we cannot hold infinity on the grid in this formalism, we must specify outer boundary conditions for the dynamic variables, $`g_{ij}`$ and $`K_{ij}`$. For the following work, we specify analytic outer boundary conditions with blending  between the analytic and numeric regions.

To specify initial data for two spinning, boosted black holes we use superposed Kerr-Schild black holes. We chose a Kerr-Schild metric \[IX\] for two reasons: 1) The metric is well defined at the event horizon, and 2) The metric is Lorentz form-invariant in a simple sense, under boosts ($`v>0`$). The superposed data are constructed in an approximate manner by a conformal method based on the superposition of two isolated, boosted Kerr black holes.

The initial data follows from Matzner et al.  and was first implemented by Correll . The data is the superposition of two, isolated, boosted Kerr-Black holes with individual metrics given by eqn.(23). The resulting superposed metric is:

$$\widehat{g}_{ij}={}_{(1)}{}^{}g_{ij}^{}+{}_{(2)}{}^{}g_{ij}^{}\eta _{ij}$$
(42)
with the $`\widehat{}`$ symbol indicating a quantity conformally related to the physical metric, $`g_{ij}=\varphi ^4\widehat{g}_{ij}`$.

$`{}_{(1)}{}^{}g_{ij}^{}`$ $`=`$ $`\eta _{ij}+{}_{(1)}{}^{}H(r_1){}_{(1)}{}^{}l_{i}^{}{}_{(1)}{}^{}l_{j}^{}\text{ and}`$ (43)
$`{}_{(2)}{}^{}g_{ij}^{}`$ $`=`$ $`\eta _{ij}+{}_{(2)}{}^{}H(r_2){}_{(2)}{}^{}l_{i}^{}{}_{(2)}{}^{}l_{j}^{}`$ (44)
are the the isolated Kerr-Schild metric forms with $`l_i`$ and $`H`$ corresponding to the single black holes. The two holes have comparable masses, $`M_1M_2`$, coordinate separation $`r_{12}`$, and velocities $`v_1`$ and $`v_2`$ assigned to them. For the argument of $`H`$ and $`l_j`$, we use

$`r_{1}^{}{}_{}{}^{2}`$ $`=`$ $`(xx_1)^i(xx_1)^j\delta _{ij}\text{ and}`$ (45)
$`r_{2}^{}{}_{}{}^{2}`$ $`=`$ $`(xx_2)^i(xx_2)^j\delta _{ij}`$ (46)
with $`x_{1}^{}{}_{}{}^{i}`$ and $`x_{2}^{}{}_{}{}^{j}`$ the coordinate positions of the holes on the initial slice.

The extrinsic curvatures of the two isolated black holes are added to obtain a trial $`\widehat{K}_{ab}`$ for the binary black hole system given as:

$${}_{(0)}{}^{}\widehat{K}_{ij}^{}={}_{(1)}{}^{}\widehat{K}_{ij}^{}+{}_{(2)}{}^{}\widehat{K}_{ij}^{}.$$
(47)
The subscript $`0`$ indicates that this is an approximation to the true extrinsic curvature of the binary black hole spacetime. $`{}_{(1)}{}^{}\widehat{K}_{ij}^{}`$ and $`{}_{(2)}{}^{}\widehat{K}_{ij}^{}`$ are the individual extrinsic curvatures associated with the isolated Kerr-Schild metric and their indices are raised and lowered by their individual metrics, eqn.(43).

For the data we describe here (holes center initially separated by a coordinate distance exceeding $`10M`$ where $`M`$ is the mass of one of the holes), we expect that an initial value solution will be $`10\%`$ in error on the domain outside of the excision volume. See further discussion in . We set the lapse function to:

$$\alpha =\alpha _1+\alpha _21,$$

and the shift vector to:

$$\beta ^i=\beta _1^i+\beta _2^i.$$

The run presented in this paper has a grid $`81\times 81\times 81`$ in Cartesian coordinates $`(x,y,z)`$ with a domain of $`(\pm 10M,\pm 10M,\pm 10M)`$ resulting in a spatial resolution of $`M/4`$. The data represent two black holes in a grazing collision. The holes are set initially at $`(5M,1M,0M)`$ and $`(5M,1M,0M)`$ in Cartesian coordinates with a boost speed of $`\pm 0.5\widehat{x}`$ toward one another and each has an angular momentum per unit mass of $`a=0.5M`$ in the $`()`$z-direction. Fig.(17) is the initial configuration of this run; note that a naive sum of the spin and the orbit angular momentum yields zero for this configuration.

We post-process the data obtained from the evolution. For the purposes of this paper, we track the apparent horizons at three specific times during the evolution, namely $`t=0M`$, $`2.8M`$, and $`3.4M`$. At $`0M`$, the apparent horizons of the initial data are found, at $`2.8M`$ two disjoint apparent horizons are found; and finally, at $`3.4M`$ a single apparent horizon is found. For the horizons shown here, the level flow method used a sphere of radius $`8M`$ to initialize each run. Fig. 18 is a plot of the horizons with time going up the page. The lowest plot is of $`t=0M`$, with each horizon being a sphere centered at coordinates $`(\pm 5M,\pm 1M,0M)`$. The middle plot shows the horizons at a later time, $`t=2.8M`$. Here the deviation in shape as the horizons accelerate towards each other is seen. The final plot at the top of fig. 18 is the first single apparent horizon that envelops both black holes at $`t=3.4M`$.

The areas for the apparent horizons at $`t=0M`$ are $`A=43.6M^2`$ for each hole. At $`t=2.8M`$, the horizons have deviated from a spherical shape and the areas for each hole are $`A=44.2M^2`$, giving a measure of the accuracy to which we can maintain their areas constant. The area of the merged apparent horizon at $`t=3.4M`$ is $`A_{merged}=184M^2`$. According to the black hole area theorem of Hawking and Ellis , the area of the merged event horizon must equal at least the combined area of the individual event horizons. Although no strong statements can be made about the area of an apparent horizon, we do find $`A_{merged}>A_1+A_2`$ in a consistent manner. We can further surmise that the final maximum area we could expect based on the initial configuration should be that of a Schwarzschild black hole of mass $`2M`$, giving an area of approximately $`201M^2`$. In some sense the area predicted by the Schwarzschild case is an upper bound. We see a $`8.5\%`$ deviation from that “idealized” case. In view of this upper limit, $`8.5\%`$ may be an indication of the greatest amount of gravitational radiation up to the time of merger ($`t=3.4M`$) given our approximate initial data, gauge condition, and boundaries.

## VIII Conclusion

Apparent horizon location and tracking constitute an important part of numerical evolutions of black hole spacetimes using excision techniques. We have demonstrated a method for finding apparent horizons in situations where the location of the apparent horizon may not be known; hence a good initial guess for the finder may not be possible. The method we have discussed works with generic 3-metric and extrinsic curvature data and with an arbitrary initial starting surface. Furthermore, the method is capable of detecting a topology change as the finder flows towards the apparent horizon. This ability is important for situations where there are multiple apparent horizons in the data. This allows us to locate apparent horizons in binary black hole evolutions without knowing where the apparent horizons are; and it allows us to locate the first single apparent horizon that forms at the merger of two black holes. The level flow method is successful at locating the apparent horizons in generic spacetimes as demonstrated by the Schwarzschild and Kerr data. It also found multiple apparent horizons starting from a single starting surface as demonstrated with the Brill-Lindquist data. Most importantly, the level flow method has been successful at identifying apparent horizons in a binary black hole evolution involving two Kerr-Schild black holes. Beginning with a single guess surface, two discrete apparent horizons were found at early times, and the later single merged horizon was found.

One of the drawbacks of this method currently is its slow convergence property due to the parabolic nature of the equation solved. This, however does not pose a problem since the level-flow method can be used in conjunction with other methods which may be more efficient given a good initial guess. The level-flow method has the definite advantage of being capable of finding multiple surfaces in the data. It can be used to get extremely good initial guesses for other methods that converge more quickly close to the solution. We are currently using the level flow method in this manner in numerical evolutions of black hole collisions.

## IX Acknowledgments

We thank Pablo Laguna for suggesting this topic to DMS, Randy Correll for supplying the initial data routine used in the evolution, and Luis Lehner for discussions on the tracker and evolution. This work was supported by NSF ASC/PHY9318152, NSF PHY9800722, NSF PHY9800725 to the University of Texas and NSF PHY9800970 and NSF PHY9800973 to the Penn State University.


format: text

text:
# Modified BCS mechanism of Cooper pair formation in narrow energy bands of special symmetry II. Matthias rule reconsidered

## I Introduction

In an early paper, Matthias pointed out that the superconducting transition temperature $`T_c`$ depends crucially on the number $`n`$ of valence electrons per atom. Matthias proposed that $`n=5`$ and $`n=7`$ valence electrons per atom are particularly favorable for superconductivity. By means of this empirical rule, which is often referred to as Matthias rule, a great number of superconducting alloys composed of transition as well as nontransition elements have been discovered.

Later, Hulm and Blaugher presented experimental data on the variation of $`T_c`$ with composition over the entire range of binary body-centered cubic solid solutions in the transition metal alloy system. They found two distinct maxima of $`T_c`$, one located between group 4 and group 5 close to the valence number $`n4.7`$, and the second located between group 6 and group 7 close to the valence number $`n6.4`$; see Fig. 1. Though the positions of the observed maxima differ from the values of $`n=5`$ and $`n=7`$ originally proposed by Matthias, the results of Hulm and Blaugher confirm the essential physical substance of the Matthias rule, namely that $`T_c`$ strongly depends on the number $`n`$ of valence electrons per atom and that this $`n`$ dependence is similar in every alloy system.

As is well-known, there is overwhelming theoretical and experimental evidence that an electron-phonon mechanism is responsible for superconductivity in the transition elements. In particular, the variation of $`T_c`$ with the number $`n`$ of valence electrons per atom is understood within the McMillan theory of superconductors with strong electron-phonon coupling. Nevertheless, it cannot be excluded that this variation of $`T_c`$ can also be understood in a new way provided that this new way does not contradict the standard theory of superconductivity.

In this paper, the Matthias rule is reconsidered within the nonadiabatic Heisenberg model which is defined by three new postulates given, e.g., in Ref. .

Niobium possesses a narrow, roughly half-filled energy band of special symmetry which in part I of this paper was called “superconducting band” ($`\sigma `$ band). It was shown that if the postulates of the nonadiabatic Heisenberg model are satisfied within this $`\sigma `$ band, then the electrons of this band form Cooper pairs at zero temperature. Just as in the Bardeen-Cooper-Schrieffer (BCS) theory, the coupling of the electrons is mediated by phonons. However, the formation of Cooper pairs is constrained in an absolutely new way by the conservation law of spin angular momentum: within a narrow, partly filled $`\sigma `$ band the crystal-spin angular momentum is conserved at zero temperature if and only if the electrons form Cooper pairs.

Within the nonadiabatic Heisenberg model, the superconducting transition temperature $`T_c^j`$ of the $`j`$th $`\sigma `$ band of a metal may be calculated in the BCS limit from the modified BCS equation

$$T_c^j=1.14\theta e^{1/N_j^\sigma (E_F)V};$$
(1.1)
see Sec. IIIC of part I. $`V`$ and $`\theta `$ denote the effective electron-phonon interaction and the Debye temperature, respectively. The nonadiabatic Heisenberg model modifies the interpretation of the BCS formula in that the parameter $`N_j^\sigma (E_F)`$ now is only a partial density of states of the electrons at the Fermi level, namely the density of states of the electrons of the $`j`$th $`\sigma `$ band.

In part I of this paper, only bcc metals with one atom per unit cell have been considered. The present paper, however, deals with materials with more than one atom in the unit cell. Therefore, single $`\sigma `$ bands must be replaced by “$`\sigma `$-band complexes” containing just as many bands as there are atoms in the unit cell; see appendix B. In equation (1.1) the index $`j`$ now labels the $`j`$th $`\sigma `$-band complex. In this paper, $`j=L`$ and $`j=U`$ denote the $`\sigma `$-band complex belonging to the lower and upper maximum of $`T_c`$, respectively.

In Sec. II free-electron energy bands will be constructed which in this paper shall approximate the $`s`$-$`d`$ energy bands of the binary bcc and hcp solid solutions in the transition metal alloy system. That means, a free-electron band structure is constructed consisting only of linear combinations of plain waves with the symmetry of a $`s`$-$`d`$ band group. Because $`\sigma `$-band complexes can already be identified within these free-electron bands, in Sec. III the parameters $`N_L^\sigma (E_F)`$ and $`N_U^\sigma (E_F)`$ will be calculated within this free-electron model. With these $`N_{L,U}^\sigma (E_F)`$, the positions of the maxima of $`T_c^L`$ and $`T_c^U`$ will be calculated by means of Eq. (1.1) without any fitting parameter, whereas the constant $`V`$ will be adjusted such that the maximum $`T_c`$ amounts to about $`10K`$. This adjustment does not influence the positions of the maxima noticeably.

In appendix B the definition of a $`\sigma `$-band complex will be given in a metal with any given space group and any given number of atoms per unit cell, and in appendix C the related spin dependent Wannier functions (spin dependent Wfs) will be defined.

## II Free-electron energy bands with s-d symmetry

The energy bands of the valence electrons of the transition elements are often referred to as $`s`$-$`d`$ bands since the valence electrons in isolated atoms of these metals have $`s`$ or $`d`$ symmetry. Consider a metal with the space group $`G`$ and $`\mu `$ atoms in the unit cell. In this paper, $`6\mu `$ energy bands of this metal are called “$`s`$-$`d`$ bands” if, first, the Bloch functions of these bands can be unitarily transformed into Wfs centered at the atoms and if, secondly, at each atom there are six Wfs with one of them having $`s`$-like and the other five having $`d`$-like symmetry. $`s`$\- and $`d`$-like means that in Eq. (1.8) of Ref.  the Wfs belong to a representation of $`G`$ which is compatible with the representation $`D_0+D_2`$ of the three-dimensional rotation group $`O(3)`$.

The band structure of the transition metals shall be approximated by a set of free-electron bands which satisfy the following three conditions.
(1) The set shall consist of $`5\mu `$ bands with the remaining $`\mu `$ bands connecting it with higher lying bands.
(2) Within the representation domain of the Brillouin zone the energy $`E_i(\stackrel{}{k})`$ of the $`i`$th band shall be calculated by

$$E_i(\stackrel{}{k})=\frac{\mathrm{}^2}{2m}(\stackrel{}{k}\stackrel{}{K}_M)^2,$$
(2.1)
where $`\stackrel{}{K}_M`$ denotes the wave vector of the reciprocal lattice point $`M`$ related according to Table I to the $`i`$th band.
(3) The bands shall be labeled by the representations given in Table II and III, respectively. Since, however, the set consists of $`5\mu `$ bands only, some representations may be absent.

### A Body-centered cubic structure

In the bcc structure we have $`\mu =1`$. Table II gives all the representations of a set of six $`s`$-$`d`$ bands in the Brillouin zone of the bcc lattice as derived from Table 2.7 of Ref. .

In the band structure of the typical transition element Nb \[the calculation of Mattheis is depicted, e.g., in Fig. 1 of part I of this paper\], the majority of the Bloch functions has the symmetry given in Table II. However, Nb and the other transition elements do not possess an ideal set of $`s`$-$`d`$ bands because one level near the Fermi energy at point $`N`$ is labeled by the $`p`$-like representation $`N_1^{}`$.

The symmetry of the free-electron wave functions (which are symmetrized plane waves) can be determined by standard group theoretical methods.

It turns out that no set of free-electron bands exists which complies with all the three conditions and which satisfies the compatibility relations in the Brillouin zone. Therefore, as suggested by the band structures of the bcc transition elements, I replace one of the two $`N_1`$ levels by the $`p`$-like $`N_1^{}`$ level. With this modification of the third condition we get the one free-electron band structure given in Fig. 2 which is used in this paper as approximation for the energy bands of the bcc solid solution alloys of the transition elements.

The symmetry structure of the free-electron bands has a close resemblance to what is found in the band structures of the bcc transition elements, whereas, of course, the free-electron energy $`E_i(\stackrel{}{k})`$ deviates considerably from the actual energy. By “symmetry structure” I mean the occurence and the connections of the representations between each other. For instance, in the free-electron case as well as in the actual bands, the $`H_{12}`$ state occures and is connected with the $`P_4`$ and the $`N_1`$ and $`N_1^{}`$ states.

### B Hexagonal close-packed structure

Table III gives the representations of a set of $`s`$-$`d`$ bands in the hcp structure as determined by the methods derived in Refs.  and . Now the set consists of twelve bands since there are $`\mu =2`$ atoms in the unit cell.

\[The representations $`R_\stackrel{}{k}`$ given in Table III satisfy Eq. (B4) if we put $`d(\alpha )=\chi _0(\alpha )+\chi _2(\alpha )`$, where $`\chi _0(\alpha )=1`$ and $`\chi _2(\alpha )`$ denote the characters of $`\alpha `$ in the representations $`D_0`$ and $`D_2`$, respectively, of the three-dimensional rotation group $`O(3)`$.\]

The free-electron band structure which can be constructed by starting from the three conditions at the top of this section, is depicted in Fig. 3. Other sets of free-electron bands which also comply with these conditions, only differ unessentially from the band structure depicted in Fig. 3. For instance, the upper $`\mathrm{\Gamma }_1^+`$ and $`\mathrm{\Gamma }_4^{}`$ states may be interchanged. Hence, I use the band structure in Fig. 3 as approximation of the energy bands of the hcp solid solution alloys of the transition elements.

In the case of the hcp structure, it is not necessary to use $`p`$-like states in the free-electron $`s`$-$`d`$ band structure in order for the compatibility relations to be fulfiled. In the real band structures of the hcp transition elements, however, $`p`$-like states are present and influence the symmetry-structure of the bands. In rhenium, e.g., a $`p`$-like state, namely $`\mathrm{\Gamma }_3^+`$, is situated near the Fermi level; see Fig. 2 of Ref. . Therefore, an additional band exists in the band structure which causes a change of the symmetry structure. In this paper, this pertubation by the $`p`$-electrons will be ignored since it is difficult to include $`p`$-like states in an unequivocal way.

## III Calculation of the positions of the maxima of $`T_c`$.

In this section I shall calculate the positions of the maxima of the superconducting transition temperature $`T_c^j`$ from equation (1.1). The density of states $`N_j^\sigma (E_F)`$ at the Fermi level is determined exactly for the free-electron bands given in Figs. 2 and 3, the Debye temperature $`\theta `$ is arbitrarily put equal to $`300K`$ and the effective spin-phonon interaction $`V`$ is chosen in such a way that, in each case, the maximum transition temperature amounts to about $`10K`$. A variation of $`V`$ or $`\theta `$ does not modify the positions of the maxima of $`T_c`$ noticeably.

Changing the Fermi energy $`E_F`$ stepwise over the entire energy range of the considered set of $`s`$-$`d`$ bands, I calculate for each $`E_F`$ first the total number $`n`$ of electrons per atom by integration and then, for each $`i`$, the density of states $`N_i(E_F)`$ of the $`i`$th band at the Fermi level. The calculated $`n`$ is referred to as the average number of valence electrons per atom.

The density of states $`N_j^\sigma (E_F)`$ is given by

$$N_j^\sigma (E_F)=\underset{i}{}d_{ij}N_i(E_F),$$
(3.1)
where

$$d_{ij}=\{\begin{array}{cc}\hfill 1& \text{if the }i\text{th band belongs}\hfill \\ & \text{to the }j\text{th }\sigma \text{-band complex}\hfill \\ \hfill 0& \text{else.}\hfill \end{array}$$
(3.2)
In the following, we put $`j=U`$ and $`j=L`$ for the upper and lower maximum of $`T_c`$, respectively.

### A Body-centered cubic structure

For an integer number $`n`$ of valence electrons, the bcc solid solution alloys studied by Hulm and Blaugher are elements with the space group $`G=O_h^9`$ and one atom per unit cell. $`\sigma `$-band complexes in this structure are simple $`\sigma `$ bands consisting of only one energy band.

If $`n`$ is not integer, the space group of the alloy is a subgroup $`G^{}`$ of $`O_h^9`$. An energy band recognized as $`\sigma `$ band in $`O_h^9`$ stays a $`\sigma `$-band (complex) in the subgroup $`G^{}`$, i.e., it stays a $`\sigma `$-band (complex) when the symmetry is lowered on the fixed atomic bcc lattice. On the other hand, an energy band not forming a $`\sigma `$ band in $`O_h^9`$ may become a $`\sigma `$-band complex when it is folded into the Brillouin zone of the subgroup $`G^{}`$.

#### 1 Lower maximum of $`T_c`$

Table IV lists the representations of all the $`\sigma `$ bands in the bcc structure. These representations are so-called extra or double-valued representations of the group $`G_\stackrel{}{k}`$ of the wave vector $`\stackrel{}{k}`$.

The only $`\sigma `$ band which can be found among the energy bands in Fig. 2 is band 3 denoted by the dotted line. This $`\sigma `$ band can be identified by means of a compatibility table for double-valued and single-valued representations. From the tables in Ref.  we get

$`D_{1/2}\times \mathrm{\Gamma }_{25}^{}`$ $`=`$ $`\mathrm{\Gamma }_7^++\mathrm{\Gamma }_8^+,`$ (3.3)
$`D_{1/2}\times H_{25}^{}`$ $`=`$ $`H_7^++H_8^+,`$ (3.4)
$`D_{1/2}\times N_2`$ $`=`$ $`N_5^+,`$ (3.5)
$`D_{1/2}\times P_4`$ $`=`$ $`P_7+P_8.`$ (3.6)
These equations show that the spin dependent Bloch functions $`\varphi _{\stackrel{}{k}qm}(\stackrel{}{r},t)`$ \[given in Eq. (C3)\] of a band with the single-valued representations

$$\mathrm{\Gamma }_{25}^{},H_{25}^{},P_4,\text{ and }N_2$$
(3.7)
can be transformed in such a way that at each symmetry point two functions form a basis of the double-valued representations

$$\mathrm{\Gamma }_7^+,H_7^+,P_7,\text{ and }N_5^+$$
(3.8)
given in the second line of Table IV. Band 3 in Fig. 2 is a $`\sigma `$ band because the representations (3.7) can be found among the representations related to band 3.

Fig. 4 shows $`T_c^L`$, as calculated with

$$N_L^\sigma (E_F)=N_3(E_F),$$
(3.9)
as function of the number $`n`$ of valence electrons per atom. The maximum at $`n4.9`$ coincides with surprising accuracy with the experimental value $`n4.7`$ of Hulm and Blaugher. Obviously, the $`T_c(n)`$ curve is determined to a large extend by the symmetry structure of the $`s`$-$`d`$ bands rather than by the real values of the energy $`E(\stackrel{}{k})`$. This idea is corroborated by the calculated band structure of Nb \[see Fig. 1 of part I\], showing that a band with the symmetry (3.7) is roughly half-filled. This band has already been considered in Fig. 1 of Ref. .

#### 2 Upper maximum of $`T_c`$

The upper maximum of $`T_c`$ is not present in the bcc structure. Band 4, as denoted in Fig. 2 by the bold line, is characterized by the (double-valued) representations

$$\mathrm{\Gamma }_7^+,H_7^+,P_8,\text{ and }N_5^+,$$
(3.10)
indicating that at point $`P`$ this band is connected with other bands, since the representation $`D_{1/2}\times P_3=P_8`$ is four-dimensional. Therefore, $`P_8`$ cannot be found in Table IV and, hence, band 4 does not form a $`\sigma `$ band. This result, too, is confirmed by Hulm and Blaugher since the curve of the upper maximum does not go back to the integer value $`n=6`$; see Fig. 1. Later measurements have shown that both Mo and W (having six valence electrons per atom) become superconducting at the relatively low temperatures $`0.92K`$ and $`0.0012K`$, respectively. Their superconducting state is clearly related to band 3 because both metals have a band with the symmetry (3.7) in their band structure which is far from half, but not completely filled; see Fig. 1 in Ref. .

Fig. 5 shows a plot of $`T_c^U(n)`$ as calculated with

$$N_U^\sigma (E_F)=N_4(E_F),$$
(3.11)
e.g., this curve is calculated by assuming that, nevertheless, band 4 is a $`\sigma `$ band. The calculated maximum of $`T_c^U(n)`$ at $`n6.5`$ coincides with the experimental value $`n6.4`$ with high accuracy. This fact together with the steep rise of the experimental $`T_c(n)`$ curve for $`n>6`$ in Fig. 1 suggests that band 4 becomes a $`\sigma `$-band complex when the ideal bcc symmetry of the electron system is disturbed, i.e., when band 4 is folded into the Brillouin zone of the space group $`G^{}`$ of the Mo-Re or W-Re alloy.

Unfortunately, in the literature I have not found any details of the symmetry of this alloy in the bcc solid solution range. Possibly, the atoms are completely irregularly distributed over the bcc lattice sites. In this case, band 4 becomes clearly a $`\sigma `$ band since the problem with the $`P_3`$ state disappears. More likely, however, the atoms in these alloys are ordered or have at least a short-range order. Therefore, I have examined band 4 with the symmetry (3.10) in four ordered structures in a body-centered cubic solid solution.

First, consider the CsCl structure with the cubic space group $`O_h^1`$ and two atoms per unit cell. The double-valued representations of the four possible $`\sigma `$-band complexes (of two bands each) in the CsCl structure are listed in Table V. Fig. 6 shows band 4 folded into the Brillouion zone of $`O_h^1`$. This set of two bands does not form a $`\sigma `$-band complex because the representations at point $`R`$, namely $`R_8^\pm `$, cannot be found in Table V. The same result we get for the NaTl structure with the space group $`O_h^7`$. Consequently, we can exclude that the CsCl as well as the NaTl structure is realized in bcc solid solutions of Mo-Re or W-Re since in both structures the measured upper maximum of $`T_c`$ cannot exist.

The MoSi2 structure with the tetragonal space group $`D_{4h}^{17}`$ has been observed rather frequently, but often with distortions of the atomic bcc lattice. Table VI lists the (double-valued) representations of the four possible $`\sigma `$-band complexes in the MoSi2 structure and Fig. 7 shows band 4 folded into the Brillouin zone of $`D_{4h}^{17}`$. The resulting set of three bands forms a $`\sigma `$-band complex as indicated by the second row of Table VI.

Other conceivable structures on an atomic bcc lattice have the hexagonal space group $`D_{3d}^3`$, as, e.g., the CdI2 structure. \[The International Tables give $`\pm (\frac{1}{3},\frac{2}{3},z)`$ for the positions of the I atoms. The atomic lattice has bcc structure for $`z=\frac{2}{3}`$.\] A hexagonal space group is suggested by the neighbouring “pseudo hexagonal” $`\sigma `$ structure (see Fig. 1). It turns out that, as in the MoSi2 structure, band 4 becomes a $`\sigma `$-band complex in the CdI2 structure when it is folded into the Brillouin zone of $`D_{3d}^3`$.

To summarize, in all the considered structures on an atomic bcc lattice, except for the CsCl and NaTl structure, the upper maximum of the transition temperature as depicted in Fig. 5 may exist.

### B Hexagonal close-packed structure

In the hcp structure we have two atoms in the unit cell and, hence, any $`\sigma `$-band complex consists of two energy bands. All the possible $`\sigma `$-band complexes in this structure are listed in Table VII. By means of the compatibility tables in Ref.  it can be verified that the set of band 5 and 6 as well as the set of band 7 and 8 denoted in Fig. 3 by the dotted and bold lines, respectively, form $`\sigma `$-band complexes. Figs. 8 and 9 show $`T_c^L`$ and $`T_c^U`$, respectively, as function of $`n`$ calculated by Eq. (1.1) with

$$N_L^\sigma (E_F)=N_5(E_F)+N_6(E_F)$$
(3.12)
and

$$N_U^\sigma (E_F)=N_7(E_F)+N_8(E_F).$$
(3.13)
Again the two maxima of $`T_c`$ at $`n4.7`$ and $`n6.7`$ coincide almost exactly with the measured values of $`n4.7`$ and $`n6.4`$, respectively. It should be noted that other possible definitions of the single free-electron bands yield slightly shifted positions of the maxima. This paper defines the bands as given in Eq. (2.1).

## IV Discussion

Within the nonadiabatic Heisenberg model, the electrons of a narrow $`\sigma `$-band complex are constrained to form Copper pairs at zero temperature because in any normal conducting state the conservation of crystal-spin angular momentum would be violated; see part I of this paper. On the basis of this modified BCS mechanism of Cooper pair formation, a theoretical interpretation of the precise measurements of Hulm and Blaugher , see Fig. 1, becomes possible within an appropriate free-electron model and within the BCS theory.

Hence, there is great evidence that the BCS theory yields the absolute energy minimum in the Hilbert space only if the formation of Cooper pairs is additionally constrained by the constraining forces existing in a narrow, partly filled $`\sigma `$-band complex. Under this assumption, relevant results (but, of course, not all the results) of the strong coupling theory of superconductivity can already be obtained within the weak coupling (BCS) limit. Within the approximations of this paper, I was able to calculate the numbers $`n`$ of valence electrons per atom related to the maxima of the transition temperature $`T_c`$, but not the values of $`T_c`$.

The calculated positions of the maxima of $`T_c`$ at valence numbers approximately equal to 4.9 and 6.5 in bcc solid solution (see Figs. 4 and 5) are in good agreement with the measured values of 4.7 and 6.4 by Hulm and Blaugher; see Fig. 1. Obviously, the positions of the maxima of $`T_c`$ are determined to a large extend by the symmetry structure of the $`s`$-$`d`$ electrons rather than by the special form of the periodic potential of the ion cores. This theoretical result emphasizes the universal character of the Matthias rule.

In order to test the influence of the atomic structure on the transition temperature, the related maxima have been calculated in hcp solid solutions, too. Here we get nearly the same positions as in bcc solid solutions, namely $`n4.7`$ and $`n6.7`$; see Figs. 8 and 9. This result suggests that the Matthias rule reflects a property of $`s`$-$`d`$ electrons which is not influenced much by the lattice symmetry.

This statement, however, should be modified: A too high symmetry of the lattice, i.e., of the periodic potential, may suppress superconductivity. A remarkable finding is that the upper maximum of $`T_c`$ at $`n6.5`$ exists neither in the transition elements with bcc structure nor in alloys with CsCl or NaTl structure, while it may occur in several ordered structures in a bcc solid solution alloy of two components. These are, e.g., the tetragonal MoSi2 structure and several atomic arrangements with the hexagonal space group $`D_{3h}^3`$. Obviously, this fact is responsible for the abrupt increase of the $`T_c`$ curve in Fig. 1 between $`n=6`$ and $`n=6.3`$.

The upper maximum in Fig. 1 is markedly narrower than the lower. This might indicate that a (short-range) order with the composition MoRe2 with $`n6.33`$ is present in the solid solution. I believe that the ordered structure has the highest $`T_c`$, though, at this stage, I cannot give a theoretical foundation for this speculation.

The actual band structures of the transition metals contain also $`p`$-like Bloch states. In this paper I have taken into account $`p`$-like states in the bcc structure but have ignored them in the hcp structure; see Sec. II B. Therefore, I cannot exclude that the calculated maxima of $`T_c`$ in the hcp structure are considerably shifted in reality. At this stage, I only can state that a pure $`s`$-$`d`$ band structure in hcp solid solutions leads to positions of the maxima which are similar to what is found in bcc materials.

###### Acknowledgements.

I wish to thank Ernst Helmut Brandt for stimulating discussions and his great help in completing this work. I thank Max Wagner and Günther Tatzel for critical comments on the manuscript and Uwe Eßmann for continuing support.

## A General

Consider a metal with the space group $`G`$, the point group $`G_0`$, and $`\mu `$ atoms in the unit cell. The positions of the atoms are written as

$$\stackrel{}{T}=\stackrel{}{R}+\stackrel{}{\varrho }_i,$$
(A1)
where $`\stackrel{}{R}`$ and $`\stackrel{}{\varrho }_i`$ ($`i=1`$ to $`\mu `$) denote the lattice points and the positions of the $`i`$th atom within the unit cell, respectively. The elements

$$a=\{\alpha |\stackrel{}{t}\}$$
(A2)
of $`G`$ consist of a point group operation $`\alpha `$ and a translation

$$\stackrel{}{t}=\stackrel{}{\tau }(\alpha )+\stackrel{}{R},$$
(A3)
which is the sum of the nonprimitive translation $`\stackrel{}{\tau }(\alpha )`$ associated with $`\alpha `$ and a primitive translation $`\stackrel{}{R}`$.

The operators $`P(a)`$ act on a function of position, $`f(\stackrel{}{r})`$, according to

$$P(a)f(\stackrel{}{r})=f(\alpha ^1\stackrel{}{r}\alpha ^1\stackrel{}{t}),$$
(A4)
and on Pauli’s spin functions,

$$u_s(t)=\delta _{st},$$
(A5)
according to

$$P(\alpha )u_s(t)=\underset{s^{}}{}d_{s^{}s}(\alpha )u_s^{}(t).$$
(A6)
The matrices $`[d_{s^{}s}(\alpha )]`$ are the representatives of the two-dimensional double-valued representation $`D_{1/2}`$ of the three-dimensional rotation group $`O(3)`$, and $`s=\pm \frac{1}{2}`$ and $`t=\pm \frac{1}{2}`$ denote the spin quantum number and the spin coordinate, respectively.

The effect of $`K`$ is given by the equations

$$Kf(\stackrel{}{r})=f^{}(\stackrel{}{r})$$
(A7)
and

$$Ku_s(t)=g_su_s(t),$$
(A8)
with

$$g_{\pm 1/2}=i.$$
(A9)
## B Superconducting band complexes

The Bloch functions belonging to a $`\sigma `$-band complex can be unitarily transformed into spin dependent Wfs which are best localized, symmetry adapted to the paramagnetic group

$$M^P=G+KG,$$
(B1)
and situated at the atoms in the sence that there is at each atom the center of symmetry of exactly one Wannier function. $`G`$ denotes the space group and $`K`$ stands for the operator of time inversion. A $`\sigma `$-band complex contains just as many bands as there are atoms in the unit cell.

In this section, I shall give the group-theoretical condition to identify such a $`\sigma `$-band complex in a given calculated band structure.

Consider a set of $`\mu `$ bands in the given metal and assume the Bloch functions in the points of symmetry $`P_\stackrel{}{k}`$ to be labeled by double-valued $`2\mu `$-dimensional representations $`R_\stackrel{}{k}^d`$ of the group $`G_\stackrel{}{k}`$ of $`P_\stackrel{}{k}`$. $`R_\stackrel{}{k}^d`$ may be irreducible or is the direct sum of irreducible double-valued representations. This set of $`\mu `$ bands forms a $`\sigma `$-band complex if two conditions are fulfiled.
(1) In a $`\sigma `$-band complex, only representations $`R_\stackrel{}{k}^d`$ occur which can be written as a Kronecker product of a single-valued representation $`R_\stackrel{}{k}`$ with the representation $`D_{1/2}`$,

$$R_\stackrel{}{k}^d=D_{1/2}\times R_\stackrel{}{k}.$$
(B2)
Again, $`R_\stackrel{}{k}`$ may be irreducible or is the direct sum of irreducible single-valued representations.
(2) In a $`\sigma `$-band complex, the matrix representatives $`D_\stackrel{}{k}(a)`$ of the representation $`R_\stackrel{}{k}`$ in Eq. (B2) satisfy, at each symmetry point $`P_\stackrel{}{k}`$, the equation

$`\text{trace }D_\stackrel{}{k}(a)`$ (B3)
$`=`$ $`d(\alpha )e^{i\alpha \stackrel{}{k}\stackrel{}{t}}{\displaystyle \underset{i=1}{\overset{\mu }{}}}n_i(a)e^{i\stackrel{}{\varrho }_i(\stackrel{}{k}\alpha \stackrel{}{k})}\text{ for }aG_\stackrel{}{k},`$ (B4)
where $`d(\alpha )`$ stands for the representatives of any real one-dimensional single-valued representation, say $`\mathrm{\Gamma }_0`$, of $`G_0`$ and

$$n_i(a)=\{\begin{array}{cc}\hfill 1& \text{if }\alpha \stackrel{}{\varrho }_i+\stackrel{}{t}=\stackrel{}{\varrho }_i+\stackrel{}{R}\hfill \\ \hfill 0& \text{else}\hfill \end{array}$$
(B5)
with $`\stackrel{}{R}`$ being a primitive translation (which, of course, may be equal zero).
Eq. (B4) is derived from Eqs. (1.8) and (4.28) of Ref.  (cf. Sec. 3 of Ref. ) assuming that there is one Wannier function with $`\mathrm{\Gamma }_0`$ symmetry at each atom. Eq. (B4) is satisfied for all $`aG_\stackrel{}{k}`$ if it is satisfied for one representant of each class of $`G_\stackrel{}{k}`$.

## C spin dependent Wannier functions

In this section the spin dependent Wfs belonging to $`\sigma `$-band complexes will be defined. Consider a $`\sigma `$-band complex consisting of $`\mu `$ bands. The Bloch functions $`\phi _{\stackrel{}{k}q}(\stackrel{}{r})`$ of this band complex can be unitarily transformed into spin dependent Wfs

$`w_{im}(\stackrel{}{r}\stackrel{}{R}\stackrel{}{\varrho }_i,t)`$ (C1)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{\stackrel{}{k}}{\overset{BZ}{}}}{\displaystyle \underset{q=1}{\overset{\mu }{}}}e^{i\stackrel{}{k}(\stackrel{}{R}+\stackrel{}{\varrho }_i)}g_{iq}(\stackrel{}{k})\varphi _{\stackrel{}{k}qm}(\stackrel{}{r},t),`$ (C2)
where the spin dependent Bloch functions $`\varphi _{\stackrel{}{k}qm}(\stackrel{}{r},t)`$ have $`\stackrel{}{k}`$-dependent spin directions,

$$\varphi _{\stackrel{}{k}qm}(\stackrel{}{r},t)=\underset{s=\frac{1}{2}}{\overset{+\frac{1}{2}}{}}f_{sm}(q,\stackrel{}{k})u_s(t)\phi _{\stackrel{}{k}q}(\stackrel{}{r}).$$
(C3)
The first sum in Eq. (C2) runs over the $`N`$ wave vectors of the Brillouin zone (BZ), $`q=1`$ to $`\mu `$ denotes the band index, the $`\mu `$-dimensional matrices $`[g_{iq}(\stackrel{}{k})]`$ are, for each $`\stackrel{}{k}`$, unitary, and the two-dimensional matrices $`[f_{sm}(q,\stackrel{}{k})]`$ are, for each q and $`\stackrel{}{k}`$, unitary. $`i=1`$ to $`\mu `$ labels the Wfs at different atoms in the unit cell and the crystal-spin quantum number $`m=\pm \frac{1}{2}`$ distinguishes between the functions at the same atom.

In a $`\sigma `$-band complex, that means, if Eqs. (B2) and (B4) are satisfied, the matrices $`[g_{iq}(\stackrel{}{k})]`$ and $`[f_{sm}(q,\stackrel{}{k})]`$ can be chosen in such a way that
(1) the spin dependent Bloch functions $`\varphi _{\stackrel{}{k}qm}(\stackrel{}{r},t)`$ vary smoothly through the whole $`\stackrel{}{k}`$ space and
(2) the Wfs are symmetry adapted to $`M^P`$ according to

$`P(a)w_{\stackrel{}{T}m}(\stackrel{}{r},t)`$ (C4)
$`=`$ $`d_i(\alpha ){\displaystyle \underset{m^{}=\frac{1}{2}}{\overset{+\frac{1}{2}}{}}}d_{m^{}m}(\alpha )w_{\stackrel{}{T}^{}m^{}}(\stackrel{}{r},t)\text{ for }aG`$ (C5)
and

$$Kw_{\stackrel{}{T}m}(\stackrel{}{r},t)=\pm w_{\stackrel{}{T}m}(\stackrel{}{r},t),$$
(C6)
with the abbreviations

$$w_{\stackrel{}{T}m}(\stackrel{}{r},t)w_{im}(\stackrel{}{r}\stackrel{}{R}\stackrel{}{\varrho }_i,t)$$
(C7)
and

$$\stackrel{}{T}^{}=\alpha (\stackrel{}{R}+\stackrel{}{\varrho }_i)+\stackrel{}{t}.$$
(C8)
The plus in Eq. (C6) is defined to belong to $`m=+\frac{1}{2}`$ and the minus to $`m=\frac{1}{2}`$. As a consequence of Eq. (B2), the two-dimensional matrices $`[d_{m^{}m}(\alpha )]`$ belong to the representation $`D_{1/2}`$ of $`O(3)`$ and the coefficients $`d_i(\alpha )`$ are c-numbers with $`|d_i(\alpha )|=1`$. The smoothness of the spin dependent Bloch functions $`\varphi _{\stackrel{}{k}qm}(\stackrel{}{r},t)`$ guarantees that the spin dependent Wfs are optimally localizable.

\[If in Eq. (B5) $`n_i(\alpha )=1`$, the c-numbers $`d_i(\alpha )`$ are equal to $`d(\alpha )`$ as given in Eq. (B4). In all cases we have $`d_i(\alpha )=\overline{D}_{ji}(\alpha )`$ with $`j`$ labeling the non-vanishing element in the $`i`$th column of the matrix $`[\overline{D}_{ji}(\alpha )]`$ given in Eq. (2.18) of Ref.  which has only one non-vanishing element in each column.\]

If in Eq. (C3) we have

$$f_{sm}(q,\stackrel{}{k})=\delta _{sm},$$
(C9)
the Wfs are usual Wfs. We have spin dependent Wfs if the matrix $`[f_{sm}(q,\stackrel{}{k})]`$ is not independent of $`\stackrel{}{k}`$.

Usual Wfs could be constructed only if the $`\sigma `$-band complex would not be connected with other bands (if we still demand that the Wfs are best localized and symmetry adapted). In the metals, however, any energy band is connected which other bands at the point and lines of symmetry. Therefore, the matrix $`[f_{sm}(q,\stackrel{}{k})]`$ cannot be chosen independent of $`\stackrel{}{k}`$.

The abbreviation in Eq. (C7) is possible since there is one Wannier function at each atom in a $`\sigma `$-band complex. Therefore, the description of the nonadiabatic Heisenberg model given in Sec. II of part I of this paper need not be changed in the case of $`\sigma `$-band complexes consisting of more than one band.


text:
# The rf control and detection system for PACO the parametric converter detector

## 1 Introduction

In this technical note we describe the rf control and detection system for a detector of small harmonic displacements based on two coupled microwave cavities (PACO). This experimental configuration was initially proposed by Bernard, Pegoraro, Picasso and Radicati ,  and later put in practice by Reece, Reiner and Melissinos , , and has been discussed in some detail in previous papers , , , . Here we just remind the basic principles underlying the detector operation.

The detector consists of an electromagnetic resonator with two resonant frequencies $`\omega _s`$ and $`\omega _a`$ both much higher than the characteristic frequency of the harmonic perturbation $`\mathrm{\Omega }`$, and which satisfy the resonance condition $`|\omega _s\omega _a|=\mathrm{\Omega }`$. In the scheme proposed in ,  the two resonant modes are obtained coupling two identical resonators; $`\omega _s`$ is the frequency of the symmetric (even) resonant mode, while $`\omega _a`$ is the frequency of the antisymmetric (odd) one. If some energy is initially stored in the symmetric mode, an harmonic perturbation can induce the transition of some energy to the initially empty level that can be extracted at frequency $`\omega _a`$. The electromagnetic power in the antisymmetric mode is proportional to the square of the amplitude of the perturbation.

In this experiment we use two cylindrical cavities coupled trough an axial iris and we change the cavity length applying an harmonic voltage to a piezo-electric crystal. To increase the sensitivity of the detector a resonator geometry and field configuration with high geometrical factor and high quality factor $`Q`$ are preferred. To avoid electron field emission from the cavity surface rf modes with vanishing electric field at the surface are mandatory. For these reasons we have chosen TE mode (TE011) superconducting cavities. The choice of frequency was imposed by the maximum dimension of the resonator that can be housed in a standard vertical cyostat; in our case the inner diameter is 300 mm, giving us enough room for a 3 GHz, TE011 resonator.

The system can achieve a great sensitivity to small harmonic displacements (up to $`\delta l/l10^{20}`$ @ $`\mathrm{\Omega }=`$1 MHz, and can be an interesting candidate for the detection of small, mechanically coupled, interactions (e.g. high frequency gravitational waves).

## 2 PACO rf system

The three main functions of rf control and measurement system of the PACO experiment are shown in figure 1.

The system performs three separate functions listed in the following.

1. The first task of the system is to lock the rf frequency of the master oscillator to the resonant frequency of the even mode of the cavity and to keep constant the energy stored in the mode.
2. The second task is to increase the detector’s sensitivity by driving the coupled resonators purely in the even mode and receiving only the rf power up- converted to the odd mode by the perturbation of the cavity walls. This goal can be obtained by rejecting the signal at the even mode frequency tacking advantage of the symmetries in the field distribution of the two modes.
3. The third task is the detection of the up-converted signal pushing the detector’s sensitivity to the limit set by the contribution of the noise sources at the operating frequency. The various noise sources have been described and discussed in a previous paper .

## 3 The rf control loop

The output of the master oscillator (HP4422B) is fed into the cavity trough a TWT amplifier giving a saturated output of 20 Watt in the frequency range of 2-4 GHz. The stored energy in the cavity is adjusted at the operating level by controlling the output of the master oscillator via the built in variable attenuator.

The output signal is sampled via a 3 dB power splitter. The $`A`$ output of the splitter is sent to the TWT amplifier, the $`B`$ output is sent, trough the phase shifter (PS), to the local oscillator (LO) input of a rf mixer acting as a phase detector (PD). The output of the rf power amplifier is fed to the resonant cavity trough a double directional coupler, and a $`180^\mathrm{o}`$ hybrid ring acting as a magic tee. The rf power enters the magic tee via the sum arm, $`\mathrm{\Sigma }`$, and is split in two signals of same amplitude and zero relative phase, coming out the tee co-linear arms 1 and 2.

The rf signal, reflected by the input ports of the cavity, enters the magic tee trough the co-linear arms. The two signals are added at the $`\mathrm{\Sigma }`$ arm and sampled by the directional coupler to give information about the energy stored in cavity allowing for the measurement of the coupling factor, quality factor, stored energy. While driving the cavity on the even mode no reflected signal is shown at the $`\mathrm{\Delta }`$ port of the magic tee where the signals coming from the co-linear arms are algebraically added to zero due to the $`180^\mathrm{o}`$ phase shift between the two signals coming by the co-linear arms.

To get the maximum of the performances of the magic tee we need to have equal reflected signals (phase and amplitude) at the cavity input ports. The equal amplitude goal can be achieved by a careful design of the input couplers; a good design guarantees us also very nice phase equality at the ports. To preserve the signal integrity we use matched output lines (in phase and amplitude) inside the cryostat. The cable used for the output lines is the best rf cable (as far as attenuation and phase stability are concerned) money can buy, nevertheless no characterisation at cryogenic temperatures exists. Because the phase shift is very sensitive to temperature inhomogeneities between the two cables and the phase difference between the two co-linear arms of the magic tee gives a quite strong signal at the $`\mathrm{\Delta }`$ port, we need to compensate for differential thermal contractions of the cables inside the cryostat, leading to phase unbalance in the feed lines. To do that we insert a phase shifter in one of the lines to reduce to a minimum the leakage of the unwanted modes on the two ports. As we will show in the following section, mode leakage of the even mode to the $`\mathrm{\Delta }`$ port sets a limit to the system sensitivity increasing the overall noise level of the detector.

Mode leakage of the odd mode to the $`\mathrm{\Sigma }`$ port reduces the system sensitivity by reducing the signal level available for detection. The converted rf power (odd mode) coupled to the input $`\mathrm{\Sigma }`$ port is lost forever. The rf system is symmetric on both the input and output ports of the detector cavity. The output ports of the cavity are coupled for a maximum output signal on the odd mode (detection mode) and the magic tee is used to reject the rf power at the frequency of the even mode. The up-converted signal (odd mode) comes out at the $`\mathrm{\Delta }`$ port of the magic tee. A fraction of the signal at the $`\mathrm{\Sigma }`$ port is fed to the rf input of the phase detector PD via a low noise rf amplifier. The intermediate frequency (IF) output of the phase detector PD is fed back to the rf master oscillator to lock the output signal to the resonant frequency of the resonator. The total phase shift around the loop is set trough the phase shifter PS, to have the maximum of energy stored in the detector. A carefully design of the servo loop amplifier (SLA) guarantees the stability of the system and the rejection of the residual noise of the master oscillator up to one MHz. The same fraction of the $`\mathrm{\Sigma }`$ output of the output magic tee is used to keep constant, to 100 ppm, the energy stored in the cavity feeding back an error signal to drive the electronically controlled output attenuator of the master oscillator.

Great deal of care is needed in tailoring the frequency response of both controls because the two loops can interact producing phase-amplitude oscillations in the rf fields stored in the cavities.

## 4 Sensitivity enhancement using the mode symmetry

The two modes of the detector cavity have (as in the case of two coupled pendulums) opposite symmetries of the fields. The rf detection system can greatly be improved if the mode symmetry information is used in the rf feed and detection system to reduce the residual components of the noise produced by the master oscillator line width. A further improvement is obtained tacking advantage of the very high selectivity of the resonator’s modes having a quality factor $`Q10^{10}`$.

Using two separate sets of ports to drive the cavity and to receive the up-converted signal (at the frequency of the odd mode), the cavity acts as a very sharp filter, with an high rejection of the signal noise (coming from the master oscillator) at the frequency of the up-converted signal. The resulting attenuation is given by the shape of the cavity resonance, a Lorentz shape with a full width half maximum $`FWHM=f/Q`$, where $`f`$ is the resonant frequency. This already low residual noise, can be even more reduced with the two magic tees. Using two equal couplers at the cavity input, driven by the two co-linear arms of a magic tee fed via the $`\mathrm{\Sigma }`$ port, we store energy in the cavity only on the even mode. Receiving the up-converted signal at the $`\mathrm{\Delta }`$ port the of the second magic tee, rejects the even mode components from the cavity by an amount given by the magic tee insulation.

In the case of an ideal magic tee the mode rejection is infinite as the tee insulation. No even mode component is transmitted trough the system and there will be no signal at the output port if the cavity is driven purely in the even mode. In the ideal case this results is obtained also in the more simple scheme used by Melissinos and Reece , , measuring the up converted power coming out of the detector along the input lines.

Our scheme gives better sensitivity and performances in the real case. The first obvious gain is the sum of the $`\mathrm{\Delta }`$ to $`\mathrm{\Sigma }`$ port insulation of the two tees, plus the possibility of adjusting separately the input and output lines to get better mode rejection. In a commercial magic tee the insulation is specified to be $`2535`$ dB over the whole bandwidth of the device (usually an octave in frequency). The reason for this quite low insulation is mainly due to the difficulty of balancing on a large range of frequency the phases of the signals coming from the two co-linear arms of the tee.

Consider two equal signals in phase

$$V=V_0\mathrm{cos}(\omega t+\varphi )$$
(1)
entering the co-linear arms of the magic tee; suppose a phase unbalance of $`\alpha `$ degrees between the two path to the $`\mathrm{\Delta }`$ port : the resulting signals at the $`\mathrm{\Delta }`$ are

$`V_\mathrm{\Delta }=V_0\mathrm{cos}(\omega t+\varphi )V_0\mathrm{cos}(\omega t+\varphi +\alpha )=`$
$`=V_0\mathrm{cos}(\omega t+\varphi )V_0\mathrm{cos}(\omega t+\varphi )\mathrm{cos}\alpha +V_0\mathrm{sin}(\omega t+\varphi )\mathrm{sin}\alpha `$ (2)
or, for a small value of $`\alpha `$

$$V_\mathrm{\Delta }=\alpha V_0\mathrm{sin}(\omega t+\varphi )$$
(3)
Now a phase unbalance as small as five degrees reduces the insulation from $`\mathrm{\Delta }`$ to $`\mathrm{\Sigma }`$ port to only 25 dB. This fairly low insulation in magic tees comes from the difficulty of designing a 90 degrees rf broad band phase shifter. A custom designed magic tee can be optimized to give a better insulation (35-40 dB) on a more limited frequency interval. Since the bandwidth-insulation product is roughly constant, pushing the requirements to high insulation (70-80 dB) reduces too much the useful bandwidth, giving not so much flexibility in the rf system design . Any mismatch between the frequencies of the cavity modes and of the magic tees results in a very severe degradation of the noise rejection for the whole system at the detection frequency, spoiling the ultimate sensitivity of the detector.

Our electronic scheme allows for an independent compensation of the magic tee phase mismatch either at the feed frequency and at the detection frequency in a flexible way: the phase mismatch is compensated using a variable phase shifter at the input of one of the co-linear arms. Getting the optimum phase at the input side will results in a pure excitation of the even (drive) mode of the two cavity system, keeping the power at the frequency of the odd (detection) mode 70 dB (the tee insulation) below the level of the drive mode. Adjusting the phase at the output will couple to the output only the odd mode components rejecting the even mode component by 70 dB. The total even-odd mode rejection of the system is the sum of the attenuation we can obtain from the two $`90^\mathrm{o}`$ hybrids; in the real world a practical limit is set by the cross talk effects between input and output limiting to 120-130 dB the maximum achievable insulation.

The price we pay for the mode selectivity improvement in our set-up is affordable and gives us the possibility of a separate and more controlled tuning of the input and output circuits. The use of two separated ports for the input and output (with separated magic tees) add some complication to the rf system. The coupling coefficient of the input and output port of the two cell cavity need to be critically coupled ($`\beta `$ = 1) to the rf source and the rf detection system. In this way we will have the optimum transfer of power to the even mode (a maximum of stored energy) and to the odd mode (a maximum in the detector output). Because the frequency and field distribution of the two modes are quite close, the input and output ports will be critically coupled to both modes. For that reason 50$`\%`$ of the even mode signal will be coupled to the idle $`\mathrm{\Sigma }`$ port at the output magic tee, and symmetrically 50$`\%`$ of the odd mode signal is coupled to the idle $`\mathrm{\Delta }`$ port at the input magic tee. The new couplings will greatly affect the detector’s performances, affecting the loaded quality factor of the cavity and changing in a substantial way the way to couple the cavity to the rf source.

At the glance it is clear that closing the two idle ports with a matched load will worsen the detector performance due to the following reasons:

1. Half of the converted power at the odd mode frequency will be lost worsening of a factor 2 the detector sensitivity.
2. The loaded quality factor of the cavity is reduced by a factor roughly 2 with the same reduction in the system sensitivity.
3. The coupling to the rf generator is reduced with a reduction of the detector efficiency (due to the input mismatch the incoming power is partially reflected and more power is needed from the rf system to store the same energy in the cavity).

A detailed analysis showed that closing the idle ports with an impedance giving a reflection coefficient amplitude one and phase zero solves the problem reflecting back to the cavity all the rf power coupled to the idle ports. Using that termination of the idle ports the two rf arms at the input and output are de-coupled (at least by the mode selectivity factor equal to the sum of the $`\mathrm{\Delta }`$ and $`\mathrm{\Sigma }`$ ports insulation of the two hybrids). An rf impedance having reflection coefficient equal to one is an open circuit or a short circuit plus a $`180^\mathrm{o}`$ phase shift. For this reason we close the two idle ports using a phase shifter and a short. Again this solution gives the possibility to fine tune the system and compensate for the effect of non-ideal elements in the rf system.

At the $`\mathrm{\Sigma }`$ port of the detection arm we insert a directional coupler to sample a tiny amount of the even mode power coming from the cavity. This transmitted power is fed to the frequency-amplitude servo loop used to lock the frequency of the master oscillator to the cavity frequency and to keep constant the energy content in the cavity. This choice gives to our rf control system a better reliability over the Melissinos-Reece control system using the reflected power from the cavity for the frequency control loop. The effect of this sampling behaves in our system as the sampling antenna used to monitor and control the rf power fed to a superconducting cavity.

As a final remark our system, despite some complexity, guarantees the following improvement over the one used in the previous experiment:

* A better rejection of the phase amplitude noise of the master oscillator obtained using the sharp resonance of the resonator itself.
* A better insulation of the drive and detector ports obtained by using separate drive and detection arms of the rf system.
* The possibility of an independent adjustment of the phase lag in the two arms giving a better magic tee insulation at the operating frequencies.
* A greater reliability for the frequency amplitude loop using the transmitted power, instead than the reflected, coming from the cavity.

## 5 Detection of the converted signal

The signal converted to the odd mode by the interaction between the mechanical perturbation and the rf fields is coupled to the $`\mathrm{\Delta }`$ port of the detection arm of the rf system and amplified by the low noise rf amplifier LNA. Great deal of care must be used in the choice of the LNA, because the noise figure of the amplifier greatly affects the detector’s sensitivity in the detection region above the mechanical resonances of the cavity (typically spanning the frequency interval 1-10 kHz) The LNA we choose is a commercial, room temperature, low noise amplifier with a 48 dB gain and quite large bandwidth of 500 MHz centered on our operating frequency of 3 GHz. The LNA Noise figure is 0.8 dB corresponding to a noise temperature of 360 K.

The converted rf signal amplified by the LNA is fed to the rf input port (RF) of a low noise double balanced mixer M1; the local oscillator port (LO) of M1 is driven by the even mode rf power (at frequency $`\omega _s`$) transmitted by the cavity. The LO input level is adjusted for a minimum noise contribution from the mixer. As shown in the previous section the input spectrum to the RF port of the mixer is composed by two signals: the first at frequency $`\omega _s`$ coming from the rf leakage of the even mode trough the detection system (even if greatly reduced by the aforementioned double rejection of the carefully tuned magic tees); the second is the converted energy on the odd mode at frequency $`\omega _a`$. Both signal are down converted by the M1 mixer giving to the IF port a DC signal proportional to the even mode leakage and the signal at frequency $`\mathrm{\Omega }`$ proportional to the odd mode excitation.

The down-converted IF output is further amplified using a low noise audio preamplifier (Stanford Research SR560). The combination of tunable built in filters of the SR560 amplifier allows us to further reject (if necessary) the DC component coming from the even mode leakage. Last the output of the audio amplifier is fed to the lock in amplifier (Stanford Research SR844) used as synchronous detector driven by the low frequency synthesizer used to drive the detector cavity trough a piezoelectric ceramic at frequency $`\mathrm{\Omega }`$.

The detection electronics for the detection of harmonic, mechanically coupled interactions, at frequency $`\mathrm{\Omega }`$ (as gravitational waves) is slightly different. Since the exact frequency and phase of the driving source is not known, we can’t perform a synchronous detection; we need to perform a auto correlation of the detector output - or to cross-correlate the outputs of two different detectors - to detect the down converted component at $`\mathrm{\Omega }`$. The output of the SR560 preamplifier is fed to a fft signal analyzer (Stanford Research SR780) able to average the input signal on a bandwidth as small as 0.1 mHz, giving us a very comfortable margin to get the ultimate sensitivity foreseen for the feasibility study of our detector.

The outlined scheme of electronic detection gives us the non-marginal benefit of being (at the first order at least) self compensating against perturbations changing in the same way the frequency of the two modes of the detector. This type of perturbation is usually produced by changes of the cavity walls due to changes of pressure of the helium bath used to keep the cavity at the operating temperature or by the thermal contraction of the cavity walls produced by changes of the operating temperature. Furthermore similar effects are produced by the radiation pressure of the electromagnetic energy stored in the even mode; this pressure acts in a symmetrical way on the cavity walls producing equal deformation of both the detector’s cells. The result of the deformation is a frequency shift of both the modes of a same amount of frequency; this kind of deformation does not affect (at the first order at least) the coupling of the cells fixing the amount of splitting between the even and odd modes. As a result the frequencies of the two modes $`\omega _s`$ and $`\omega _a`$ will change under the effect of pressure changes on the cavity walls but the difference $`\mathrm{\Omega }`$ (related to the coupling) will not.

In our scheme the local oscillator (used for the down conversion of the rf signal coming from the detector) is the rf signal transmitted by our cavity at the frequency of the even mode. In this way the common mode frequency drift of the resonators is automatically cancelled (at the first order at least) and the mixer M1 output is a signal exactly at the perturbation frequency $`\mathrm{\Omega }`$. This effect gives us the possibility of using a quite simple refrigeration scheme (without any complex and cumbersome servo loop) to keep constant the helium bath temperature and pressure.

## 6 Experimental results

The electromagnetic properties of the cavity have been measured in a vertical cryostat after careful tuning of the two cells frequencies. The symmetric mode frequency was $`\omega _s/(2\pi )=3.03431`$ GHz and the mode separation was $`\mathrm{\Omega }=1.3817`$ MHz. The unloaded quality factor at 4.2 K was $`Q=510^7`$, and no significant improvement was found lowering the helium bath temperature at 1.8 K. Even after a second chemical polishing, performed at CERN, which removed approximately $`300\mu \mathrm{m}`$ of niobium from the surface, no improvement was observed. We believe that this very low $`Q`$ value is due to hot spots on the surface caused by welding problems occurred during the cavity fabrication. Adjusting the phase and amplitude of the rf signal entering and leaving the cavity, the arms of the two magic tees were balanced to launch the even mode at the cavity input and to pick up the odd mode at the cavity output. With 30 dBm of power at the $`\mathrm{\Sigma }`$ port of the first magic tee, -90 dBm were detected at the $`\mathrm{\Delta }`$ port of the second one, giving an overall attenuation of the symmetric mode of 120 dB. The energy stored in the cavity with 30 dBm input power was approximately 1 mJ. The signal emerging from the $`\mathrm{\Delta }`$ port of the output magic tee was amplified by the LNA and fed into a spectrum analyzer. The signal level at frequency $`\omega _a`$ was 120 dBm in a 1 Hz bandwidth. System sensitivity at this stage is given by

$$h_{min}6.5510^{18}\mathrm{Hz}^{1/2}$$
(4)
This value is quite far from our goal of $`h_{min}=10^{20}`$ (Hz)-1/2. Since, to our knowledge, this is the first example of a parametric detector operated in trasmission, and since this configuration requires very careful adjustments of the input and output ports balancing, we believe that significant improvements are obtainable. Furthermore the new cavity under construction should show the high quality factor needed to reach high sensitivity.