text:
# Electronic polarization in pentacene crystals and thin films

## I Introduction

Prospective molecular electronic devices are based on organic films deposited on metal or semiconductor surfaces . Due to limited mobility, organic devices are typically restricted to thin (10—100 nm) films. Variations in electronic polarization energies of charge carriers near surfaces and interfaces affect the transport states  and hence device function. We have recently addressed electronic polarization at surfaces and in thin films  using a self-consistent approach in the limit of vanishing intermolecular overlap , when each molecule experiences the non-uniform electrostatic potential of all other molecules. For the prototypical hole conductor , perylenetetracarboxylic acid dianhydride (PTCDA), the calculated transport gap varies by 500 meV between monolayers and thick films, in agreement with photoelectron and tunneling electron spectroscopy . PTCDA showed that electronic polarization is accessible to self-consistent computation in crystalline thin films.

In contrast to PTCDA, whose molecules form one-dimensional stacks and lie almost flat on metallic substrates, many organic molecular crystals of interest have a herring-bone packing, with molecules oriented across molecular layers as sketched in Fig. 1. Both thiophenes  and acenes  have herring-bone packing and are suitable for thin film transistors. The molecules in crystalline thin films are almost upright and charge transport is preferentially parallel to the surface.

In this paper, we analyze electronic polarization in pentacene, considering it as a representative of a wider class of organic materials with electronic applications. We calculate electronic polarization in the bulk to obtain the optical dielectric constant, the transport gap for generating an electron-hole pair at infinite separation, and the electrostatic binding between molecular ions at fixed separation in the lattice. We then consider the experimental situation in Fig. 1 of thin films on a metallic substrate to compute electronic polarization at the surface and at the organic-metal interface, assuming idealized films with structure identical to the bulk. Crystalline thin films exhibit multiple phases that, except in monolayers, are close to the bulk.

Fig. 1 Schematic molecular packing in a pentacene layer on metal. Image charges below the metal surface are also sketched.

The long axis of pentacene or $`\alpha `$-sexithiophene is almost normal to the surface. Image charges across the interface in Fig. 1 then act in the direction of greatest polarizability, which results in contrasting electronic polarization in pentacene and PTCDA films. Electric fields normal to the surface induce large dipoles in pentacenes, but these induced dipoles are parallel and hence repulsive. There is competition between charge redistribution due to image charges and induced dipoles. Our self-consistent calculations of charge redistribution indicate that even a single molecular layer provides a screening environment that reduces repulsion between induced dipoles. The improved iterative procedures of Section II are necessitated by the competition between image charges and induced dipoles.

The paper is organized as follows. In Section II we consider the system of linear equations for electronic polarization and introduce alternative methods based on a minimum principle when the efficient iteration scheme of  fails. Section III contains results for bulk crystals, including the optical dielectric constants, polarization energies of cations and anions, and the binding energies of ion pairs. In Section IV we consider pentacene films of 1—10 molecular layers and report polarization energies at the surface and at the metal-organic interface. The self-consistent pentacene potential at the surface illustrates the competition between image charges and induced dipoles. Section V contains discussion and conclusions.

## II Solution of linear equations

The zero-overlap approximation reduces electronic polarization to charge redistribution on molecules in the electrostatic potential of the solid. The crystal structure fixes all distances. We neglect lattice relaxation whose polarization is estimated to be an order of magnitude smaller in organic molecular crystals . The crystal potential at atoms is readily introduced as site energies in semiempirical theory, which then provides a practical approach to computing self-consistent atomic charges and induced dipoles . We have eight linear equations per atom and need large ($`10^3`$) clusters of molecules for long-range Coulomb interactions. The large number of equations precludes the use of algorithms based on transformations of dense matrices. For example, a cluster of 2000 PTCDA molecules, as in , leads to 608,000 linear equations for 38 partial atomic charges and $`3\times 38`$ components of induced atomic dipoles for each atom of each molecule, and the same number of potentials and potential gradients.

We review briefly some approaches for solving linear equations that appear in polarization problems. We first recast these equations in a formal but transparent form. Let $`q`$ be the vector describing charge redistribution in a cluster of interest. Its components may include partial atomic charges for each atom of each molecule, as well as induced atomic dipoles. Higher atomic multipoles can be included in $`q`$ as well. Let $`\rho ^{(0)}`$ be the “unrelaxed” charge distribution, that is the charge distribution in the individual molecules in the gas phase, the actual charge distribution thus being $`\rho ^{(0)}+q`$.

For zero overlap, charge conservation for individual molecules leads to simple constraints on the components of $`q`$, meaning that $`q`$ belongs to a certain subspace $`𝒬`$ of interest. We note that $`\rho ^{(0)}`$ does not necessarily belong to $`𝒬`$ since the source charges on ions do not sum to zero.

The polarization problem is defined in terms of the “state vector” $`q`$ as minimization of the total energy $`E(q)`$, which can be formally written as

$$E=\frac{1}{2}(qGq)+\frac{1}{2}(\rho ^{(0)}+q|V|\rho ^{(0)}+q).$$
(1)
We use bra- and ket- notation to denote matrix vector products. The first term is the energy of non-interacting molecules. A symmetric positive-definite stiffness matrix $`G`$ ensures that $`q=0`$ in the absence of interactions. It results from solving Schrödinger’s equation for a single molecule in non-uniform external field. The stiffness matrix is block-diagonal, with separate blocks corresponding to individual molecules or ions.

The second term describes intermolecular interactions. The interaction matrix $`V`$ contains charge-charge, charge-dipole, and dipole-dipole interactions, and can also be made symmetric. The matrix $`V`$ is indefinite, reflecting that Coulomb systems are always unstable. In contrast to the stiffness matrix, the diagonal blocks of $`V`$ are zero, since no molecule interacts with itself.

Expression (1) shows a potential for instability, since it may have no minimum when $`V`$ is sufficiently large. In the organic molecular crystals that we consider, $`V`$ is small enough and only results in a shift of the minimum of $`E(q)`$ from $`q=0`$.

Differentiating (1) with respect to $`q`$ we find a set of linear equations on $`q`$:

$$Gq+Vq=V\rho ^{(0)},$$
(2)
subject to the constraint that $`q`$ belongs to $`𝒬`$. While this form of polarization equations is simple conceptually, another form is often more practical. By introducing a set of generalized potentials $`p=Gq`$ we recast (1) in a dual form:

$`p`$ $`=`$ $`V(\rho ^{(0)}+q)`$ (4)
$`q`$ $`=`$ $`\mathrm{\Pi }p,`$ (5)
where $`\mathrm{\Pi }=G^1`$. These are formally written Eqs. (8) and (15) of . The physical meaning of Eq. (II) is simple. Instead of a one-step minimization of $`E(q)`$ as in (2), we optimize the charge distribution of each molecule individually in the self-consistent external field $`p`$ of all other molecules.

The components of the vector $`p`$ contain potentials and potential gradients (and possibly higher derivatives), that couple to atomic partial charges, atomic induced dipoles, and higher atomic multipoles. The symmetric positive-definite polarization matrix $`\mathrm{\Pi }`$ has the same block-diagonal structure as $`G`$ and describes the polarization response of individual molecules to the external field, given by the components of $`p`$. The inverse of $`G`$ is well defined in the subspace of interest $`𝒬`$.

### A Iterative solution of polarization equations

The dual form of Eqs. (II) suggests a natural iterative approach. We start with the bare charge distribution $`\rho ^{(0)}`$ and find “bare” potentials $`p^{(0)}=V\rho ^{(0)}`$. We then repeat updating charge distribution using Eq. (5) and recalculating $`p`$ using Eq. (4) until convergence.

When convergent, this procedure is usually fast. However, it is often unstable, leading to oscillatory behavior. A simple numerical trick has been proven very powerful in suppressing such oscillations. We introduce a damping factor $`0<f<1`$ and iterate, starting with $`q^{(0)}=0`$:

$`p^{(m)}`$ $`=`$ $`V(\rho ^{(0)}+q^{(m)})`$ (7)
$`q^{(m+1)}`$ $`=`$ $`fq^{(m)}(1f)\mathrm{\Pi }p^{(m)}`$ (8)
We increase $`f`$ adaptively when oscillations in $`E(q^{(m)})`$ are detected and decrease it when convergence is stable. All results in Refs.  were obtained using this iteration procedure, which appears to be both remarkably stable and fast. For example, it took just 25 iterations to achieve $`10^7`$ accuracy in the components of $`q`$ for the cluster of 2000 molecules of PTCDA, mentioned above.

Iterations (II A) usually perform quite well. However, they are not guaranteed to converge and sometimes fail. For example, we found intermittent convergence problems for a pentacene cation in a layer on a metal substrate. Image charges in the metal induce strong charge redistribution along the pentacene axis with largest polarizability. On the other hand, the surrounding molecules in the layer have the opposite effect of drawing charge into the inside of the layer, closer to the center of the ion. As follows from the self-consistent solution, the second effect outperforms the first. This clearly leads to difficulties because the initial iterations tend to drive the system away from the convergence point, indicating also that a perturbative solution may be counterintuitive in this case.

Computational difficulties of this kind are rare. We now describe an alternative approach to solving polarization equations, which is formally stable and is guaranteed to converge, while it may not be as efficient as Eqs. (II A).

### B Subspace methods for polarization problems

Let us return to the original Eq. (2). We can treat it as a general linear system

$$Aq=b,$$
(9)
with a symmetric matrix $`A=G+V`$ on the left-hand side. Solution of such a linear system is equivalent to minimizing the energy functional

$$E(q)=\frac{1}{2}(qAq)(bq),$$
(10)
whose minimum exists when $`A`$ is positive-definite.

We present a variationally-stable subspace approach to solving (9) based on the minimum principle (10). By re-scaling the solution we can assume the right-hand side to be normalized, $`(bb)=1`$. We shall be searching for the solution $`q`$ in the Krylov subspace $`𝒦`$ of the matrix $`A`$, generated by the vector $`b`$: $`𝒦=\text{span}(b,Ab,A^2b,\mathrm{},A^{m1}b)`$. At every step $`m`$ we add a new vector to the subspace $`𝒦`$ and construct the best vector $`q`$ which minimizes (10) within $`𝒦`$. This guarantees that each step $`m`$ yields a better solution, and thus the procedure converges.

We use the Hermitian Lanczos recursion  to build an orthonormal basis in $`𝒦`$: $`q^{(1)}=b`$,

$$q^{(m+1)}=\beta _{m+1}^1(Aq^{(m)}\alpha _mq^{(m)}\beta _mq^{(m1)})$$
(11)
where $`\alpha _m`$ and $`\beta _m`$ are chosen at each step $`m`$ to orthogonalize $`q^{(m+1)}`$ with respect to both $`q^{(m)}`$ and $`q^{(m1)}`$, ensuring that $`(q^{(m)}|q^{(n)})=\delta _{mn}`$ for every $`m`$ and $`n`$. Expressing $`q=_mc_mq^{(m)}`$ and solving for the coefficients $`c_m`$ that minimize $`E(q)`$ we find that they obey a linear system $`_n\stackrel{~}{A}_{mn}c_n=\stackrel{~}{b}_m`$, where $`\stackrel{~}{A}_{mn}=(q^{(m)}Aq^{(n)})`$ and $`\stackrel{~}{b}_m=(q^{(m)}b)`$.

The algebraic properties of Lanczos recursion (11) yield readily all the components of $`\stackrel{~}{A}`$ and $`\stackrel{~}{b}`$. The matrix $`\stackrel{~}{A}`$ is symmetric tridiagonal with the diagonal elements $`\stackrel{~}{A}_{mm}=\alpha _m`$ and sub-diagonal elements $`\stackrel{~}{A}_{m,m1}=\beta _m`$. Due to the special choice of $`q^{(1)}=b`$, all $`\stackrel{~}{b}_m=0`$, except $`\stackrel{~}{b}_1=1`$ and $`\stackrel{~}{b}_2=\alpha _1`$.

The method described above is variationally stable, since it is based on a minimum principle. Its drawback, however, is that it is generally slower than iterations (II A). The reason is that Eqs. (II A) follow from the dual form of Eqs. (II), which is based in turn on a physical insight of solving each molecule separately in the external field of other molecules. The subspace approach is not based on such an insight. In practice, it takes 200—300 iterations, rather then 25-35, to achieve convergence with this approach. Thus we use the subspace approach only when iterations fail.

### C Subspace method for non-positive definite cases

The subspace method of the above section being applied to Eq. (2) deals with the stiffness matrix $`G`$. It is preferable numerically to work with $`\mathrm{\Pi }`$, rather than with its inverse. The reason is that whenever there is a small polarizability, e.g. in the direction normal to the plane of a $`\pi `$-conjugated molecule, it translates into a small eigenvalue of $`\mathrm{\Pi }`$ and, correspondingly, into a large eigenvalue of $`G`$. Numerically, it is much easier to handle near-zero than near-infinity.

In a more drastic situation, $`G`$ may have negative eigenvalues, yet the polarization self-consistent equations (II) continue to make sense. For example, we have introduced atomic polarizabilities $`\stackrel{~}{\alpha }`$  as the difference $`\alpha \alpha ^\mathrm{C}`$ between the best available molecular $`\alpha `$ and $`\alpha ^\mathrm{C}`$ based on semiempirical atomic charges. Usually $`\alpha ^\mathrm{C}`$ underestimates $`\alpha `$ by $`10`$%. But the charge-induced polarizability occasionally exceeds $`\alpha `$ in some directions (see e.g. Table I), and the correction $`\stackrel{~}{\alpha }`$ has one or more negative eigenvalues. The self-consistent polarization equations (II) then give a saddle point of the energy functional, which formally has no minimum.

In order to use the subspace method to deal with $`\mathrm{\Pi }`$ rather than $`G`$, we rewrite Eqs. (II) in a matrix form in such a way that the matrix on the left-hand side is symmetric:

$$\left(\begin{array}{cc}V& 1\\ 1& \mathrm{\Pi }\end{array}\right)\left(\begin{array}{c}q\\ p\end{array}\right)=\left(\begin{array}{c}p^{(0)}\\ 0\end{array}\right).$$
(12)
This is again a linear system of type (9), but of twice the size. The matrix is, however, not positive-definite, and so the above method based on the minimization of (10) cannot be applied.

We present a modified subspace approach that is applicable to linear systems of type (9) with symmetric matrix $`A`$ that is not necessarily positive definite. Instead of searching for the minimum of (10) we minimize the residual norm.

$$Z(q)=\frac{1}{2}(Aqb|Aqb),$$
(13)
which is also equal to half the norm of the gradient, $`E(q)/q=Aqb`$. Again, expressing the solution as $`q=_mc_mq^{(m)}`$, but solving now for $`c_m`$ that minimize $`Z(q)`$ we find another linear system, $`_n\stackrel{~}{Z}_{mn}c_n=\stackrel{~}{w}_m`$, where $`\stackrel{~}{Z}_{mn}=(q^{(m)}A^2q^{(n)})`$ and $`\stackrel{~}{w}_m=(q^{(m)}Ab)`$.

After some algebra, we find that the matrix $`\stackrel{~}{Z}`$ is symmetric 5-diagonal with the following nonzero elements: $`\stackrel{~}{Z}_{mm}=\alpha _m^2+\beta _m^2+\beta _{m1}^2`$, $`\stackrel{~}{Z}_{m,m1}=(\alpha _m+\alpha _{m1})\beta _{m1}`$, and $`\stackrel{~}{Z}_{m,m2}=\beta _{m1}\beta _{m2}`$, while the right-hand side has all components zero except $`\stackrel{~}{w}_1=\alpha _1`$ and $`\stackrel{~}{w}_2=\beta _1`$. The above expressions for the elements of $`\stackrel{~}{Z}_{mn}`$ can be found by consequently multiplying Eq. (11) by itself with $`m`$ replaced with $`m1`$, $`m2`$, etc. We notice by inspection that $`\stackrel{~}{Z}=\stackrel{~}{A}^2`$, except for the last row and column.

Throughout this section we have assumed that all vectors are projected onto the subspace $`𝒬`$ to preserve charge conservation relations. We did not discuss preconditioning techniques, such as rescaling the variables, which may improve convergence.

## III Polarization in Pentacene Crystals

### A Gas-phase properties

We performed calculations of the electronic structure of pentacene in gas phase using B3LYP hybrid density functional theory with the extensive 6-311++G(d,p) basis set available in Gaussian 98 . D2h symmetry was assumed and the geometry was separately optimized for the neutral molecule, cation, and anion. Unrestricted Hartree-Fock (UHF) was used for the ions. Limited spin contamination was observed as estimated by the maximum deviation of the total spin $`\overline{S^2}=0.769`$ before annihilation. Comparison of the energies of the ions in the neutral and ionic optimized geometries yields the relaxation energies $`\lambda _+=46`$ meV for the cation and $`\lambda _{}=68`$ meV for the anion, under the assumption of equal zero-point contributions for the neutral molecule and ions. $`\lambda _+`$ compares well with the value of 59 meV derived from experiment .

The vertical ionization potential is found to be $`I=6.229`$ eV, which is significantly below the recently reported experimental value of 6.589 eV . Other values of 6.64 eV  and 6.74 eV  were reported previously. The calculated value of $`I`$ must be corrected by the difference in zero-point energies before the direct comparison can be made, but zero-point differences are probably less than 360 meV.

The experimental value for the electronic affinity is known with an error bar. Our calculated B3LYP value $`A=1.475`$ eV is 105 meV above the recommended average 1.37 eV  of experimental data. Thus, B3LYP values lead to the gas-phase pentacene charge gap $`IA=4.754`$ eV (without zero-point correction) vs. 5.22 eV derived from experiment.

### B Polarization energies

Following the procedure in  we first calculate atom-atom polarizabilities $`\mathrm{\Pi }_{ij}`$ for individual pentacene molecules using the semiempirical INDO/S Hamiltonian . Experimental geometries from the crystal structure data  were used. Positions of hydrogen atoms, not known experimentally, were AM1-optimized using Gaussian 98 with the heavy atoms fixed.

Due to the inequivalence in crystal field for the two molecules in the unit cell, their geometries are slightly different. There are two atom-atom polarizability tensors $`\mathrm{\Pi }_{ij}`$ with $`i`$, $`j`$ ranging over 36 atoms of a pentacene,

$$\mathrm{\Pi }_{ij}=\frac{\rho _i}{\varphi _j}=\frac{^2E}{\varphi _i\varphi _j}$$
(14)
$`E`$ is the INDO/S ground-state energy, $`\varphi _i=\varphi (𝐫_i)`$ is the crystal potential at atom $`i`$ and $`\rho _j`$ are Löwdin atomic charges. Charge redistribution is then described explicitly as

$`\rho _i`$ $`=`$ $`\rho _i^{(0)}{\displaystyle \underset{j}{}}\mathrm{\Pi }_{ij}\varphi _j`$ (16)
$`𝝁_i`$ $`=`$ $`\stackrel{~}{\alpha }_i𝑭_i,`$ (17)
\[cf. Eqs (II)\], where $`\rho _i^{(0)}`$ are gas-phase charges, $`𝝁_i`$ are induced atomic dipoles and $`𝑭_i=\varphi (𝐫_i)`$ are electric fields. $`\stackrel{~}{\alpha }=\alpha \alpha ^\mathrm{C}`$ is a correction  that accounts for the difference in ab-initio molecular polarizability and the “charges-only” part $`\alpha ^\mathrm{C}`$ associated with $`\mathrm{\Pi }_{ij}`$:

$$\alpha _{\alpha \beta }^\mathrm{C}=\underset{ij}{}\mathrm{\Pi }_{ij}r_i^\alpha r_j^\beta .$$
(18)
We computed $`\alpha `$ for the two geometries using B3LYP density functional with 6-311++G(d,p) basis set of Gaussian 98.

Principal values of $`\alpha `$ and “charges-only” $`\alpha ^\mathrm{C}`$ are compared in Table I. Charge redistribution accounts quite well for the in-plane components $`\alpha _{\mathrm{LL}}`$ and $`\alpha _{\mathrm{MM}}`$ of the neutral molecule. The polarizability $`\alpha _{\mathrm{NN}}`$, normal to the molecular plane, is completely “atomic,” as expected, and cannot be described by charge redistribution . Values that are not exactly zero are due to non-planar molecules in the crystal. The INDO/S results for ions are far too large; we return to this point in the Discussion.

Table I. Molecular polarizability of pentacene along principal axes ($`L=\text{long}`$, $`M=\text{medium}`$, $`N=\text{normal}`$)

Molecule $`\alpha _{\mathrm{NN}}`$<sup>3</sup>) $`\alpha _{\mathrm{MM}}`$<sup>3</sup>) $`\alpha _{\mathrm{LL}}`$<sup>3</sup>)     B3LYP/6-311++G\**    Neutral(1) 17.77 37.66 91.43  Neutral(2) 18.02 38.02 99.48  Anion(1) 19.87 40.16 135.47  Anion(2) 20.20 40.81 131.88  Cation(1) 16.31 35.91 130.87  Cation(2) 16.54 36.54 124.14   $`\alpha ^\mathrm{C}`$, Eq. (12)    Neutral(1) 0.06 39.26 87.41  Neutral(2) 0.05 39.61 95.83  Anion(1) 0.06 35.26 194.35  Anion(2) 0.05 39.98 231.24  Cation(1) 0.06 36.03 368.28  Cation(2) 0.05 40.55 282.52

The difference between the B3LYP and $`\alpha ^\mathrm{C}`$ in Table I is distributed over the atoms with the weight proportional to the atomic valence charge . Since charge redistribution overestimates the in-plane molecular polarizability in pentacene, the atomic correction tensor $`\stackrel{~}{\alpha }=\alpha \alpha ^\mathrm{C}`$ is not positive definite. It is still small, however, for neutral molecules, and the self-consistent solution is well defined, even though it does not correspond strictly to the minimum of the total energy (cf. discussion in II.C).

### C Properties of the neutral lattice

Both pentacene  and anthracene  illustrate herring-bone packing with two molecules per unit cell, but their Bravais lattices are different. Pentacene is triclinic . The two molecules are inequivalent and subject to different crystal-field environments. The polarization energy of a cation or anion consequently depends on which molecular site is charged. Since crystalline electric fields are perturbations, we expect comparable polarization energies that differ at most by $`100`$ meV. Such differences have already been computed  in pentacene crystals using the submolecular method  in which the gas-phase molecular polarizability $`\alpha `$ is placed as $`\alpha /5`$ at the ring centers.

Solving self-consistent polarization equations for the neutral lattice of pentacene, we find that polarization contribution to sublimation energy is negligible (1.23 meV per molecule), as in anthracene, since the gas-phase charges are small due to an approximate electron-hole symmetry of the valence shell. We find the dielectric tensor of pentacene to be highly anisotropic, with principal values $`\kappa _1=5.336`$, $`\kappa _2=3.211`$, and $`\kappa _3=2.413`$. For reference purposes, the directional cosines of the principal axes are $`𝐧_1=(0.296,0.314,0.902)`$, $`𝐧_2=(0.021,0.946,0.322)`$, and $`𝐧_3=(0.955,0.076,0.286)`$, respectively, in the Cartesian coordinate system with the pentacene lattice vectors (in Å) $`𝐚=(7.900,0,0)`$, $`𝐛=(0.444,6.044,0)`$, and $`𝐜=(6.153,2.858,14.502)`$. The direction of $`𝐧_1`$ coincides with the direction of long axes of pentacene molecules to the accuracy to which the latter can be defined ($`5^o`$).

The calculated principal axes and values of $`\kappa `$ for anthracene crystals agree with experiment . Measurements of refractive indices or optical dielectric constants are challenging in anisotropic molecular crystals . We are not aware of such data for pentacene crystals. Dielectric data will be necessary for organic devices with improved performance.

### D Charge carriers in the bulk

We compute the polarization energy of charge carriers following , by placing an ion in the infinite neutral lattice and considering an imaginary sphere of $`M`$ molecules surrounding it. We allow only the molecules whose centers fall into the sphere to relax their charge distributions from the self-consistent values for the neutral lattice and monitor convergence as $`M`$ is increased. We use the same molecular geometry for the ion. The small ($`100`$ meV) relaxation energies of large aromatic molecules make the same geometry both an excellent and convenient approximation. Atom-atom polarizabilities $`\mathrm{\Pi }_{ij}`$ for the ion are computed using unrestricted Hartree-Fock (UHF) INDO/S.

Fig. 2 Convergence of $`P_\pm `$ with increasing cluster size $`M`$.

Figure 2 shows linear convergence of polarization energy with $`M^{1/3}`$, which is proportional to the inverse radius of the sphere. The polarization energy extrapolates as $`M^{1/3}`$, because the missing part at large finite $`M`$ can be thought of as polarization of an infinite homogeneous dielectric with a spherical cavity of radius $`RM^{1/3}`$, $`(1\kappa _{\mathrm{eff}}^1)e^2/R`$ . Extrapolation to infinite radius gives the polarization energy in the bulk, $`P=P_++P_{}=2.014`$ eV and $`2.007`$ eV for the ion occupying positions 1 and 2 in the unit cell, respectively. The “charges only” values of $`P=P_++P_{}`$ are $`10`$% less, as also found in other organic crystals. Since $`\kappa _{\mathrm{eff}}`$ for an anisotropic dielectric medium is known , the slope in Fig. 2 at large $`M`$ is directly related to the $`\kappa _\alpha `$ found above, and the slopes agree within a few percent. We note that $`\kappa _{\mathrm{eff}}^1`$ implies identical slopes at large $`M`$ for the anion and cation in either position, as found separately in Fig. 2.

Using B3LYP values for $`I`$ and $`A`$ in the gas phase, we calculate the transport gap of the pentacene crystal $`E_t=IAP=2.740`$ eV in the limit of zero overlap. This correlates quite well with the reported band gap of 2.85 eV, obtained from careful fit of electroabsorption spectra . $`E_t=2.78`$ eV was reported in  also based on electroabsorption spectra of charge-transfer (CT) states. The experimental value for the charge gap, $`IA=5.22`$ eV, used in this context results in $`E_t=3.21`$ eV, which is too high.

### E Ion pairs and CT states

We next report polarization energies $`P_{\mathrm{pair}}`$ of CT states with an ion at the origin and the counter ion at a nearby site. We place the ion pair within an imaginary sphere in a neutral lattice and calculate the effective interaction $`V_{\mathrm{eff}}`$, which is the difference of the total polarization energy and the energy of a pair of well-separated charges, $`V_{\mathrm{eff}}=(P_{\mathrm{pair}}P_+P_{})`$. $`V_{\mathrm{eff}}`$ is the binding energy of CT states in the limit of no overlap and is closely related to the binding energies of Frenkel excitons, when the final state is a molecular excitation.

Fig. 3 shows convergence with the size of the sphere. Since the total ion charge within the sphere is zero, the leading term $`M^{1/3}`$ vanishes. The polarization energy of a dipole in a dielectric medium converges much faster as $`M^1`$. Extrapolating to infinite $`M`$, we obtain the lowest CT states at $`0.719`$ and $`0.679`$ eV for the nearest-neighbors (center-to-center distance 4.799 Å) and $`0.685`$ and $`0.675`$ eV for next-nearest neighbors (5.151 Å), as listed in Table II. These results compare reasonably with submolecular calculations  of $`0.777`$ and $`0.698`$ eV, respectively, with $`\alpha /5`$ at ring centers. Submolecular results do not distinguish between anions and cations and also depend on the precise partitioning of $`\alpha `$, which is left open.

Fig. 3 Convergence of the effective electron-hole interaction energy with cluster size $`M`$.

Table II. Effective interaction energies $`V_{\mathrm{eff}}(𝐫)`$, in eV, of cation-anion pairs. The submolecular data is from . Center-to-center distances $`r`$ are also listed for reference.

Anion Cation $`r`$, Å  $`V_{\mathrm{eff}}(𝐫)`$ $`V_{\mathrm{eff}}^{\mathrm{submol}}(𝐫)`$    $`(0,0,0)`$ $`(\frac{1}{2},\frac{1}{2},0)`$ $`4.7990`$ $`0.719`$ $`0.756`$  $`(\frac{1}{2},\frac{1}{2},0)`$ $`(0,0,0)`$ $`4.7990`$ $`0.679`$ $`0.756`$  $`(0,0,0)`$ $`(\frac{1}{2},\frac{1}{2},0)`$ $`5.1514`$ $`0.685`$ $`0.698`$  $`(\frac{1}{2},\frac{1}{2},0)`$ $`(0,0,0)`$ $`5.1514`$ $`0.675`$ $`0.698`$

It is interesting that the mutual orientation of the cation and ion together with variations in their electrostatic response are enough to reverse the order of states, such that the pair with larger center-to-center distance has greater binding energy (lines 2 and 3 in Table II). This is an effect very similar in nature to the one suggested by Mazur and Petelenz .

## IV Polarization in Pentacene Films

### A Inert metal-organic interfaces

We now address polarization energy in thin pentacene films deposited on metal surface (Fig. 1). We use an approach similar to the one sketched recently for thin films of PTCDA , which we present in greater detail. We assume an ideal metal-organic interface without chemical interactions. Organic molecules are physisorbed and in van der Waals (vdW) contact with the metal. We model the metal as an equipotential surface parallel to the organic layer. The potential at the metal is assumed to be $`\varphi =0`$, as the actual value drops out from the combined polarization energy for the electron and hole, $`P=P_++P_{}`$, even though we perform separate calculations for $`P_+`$ and $`P_{}`$. In fact, any additive constant in the potential cancels out for any neutral entity with an arbitrary but equal numbers of positive and negative charges.

For the organic film we use the bulk crystal geometry, assuming a parallel slab of the crystal cut along molecular layers and placed parallel to the metal surface. The only parameter in the model that allows limited adjustment is the actual distance of the equipotential plane from the slab. We choose it such that the equipotential surface is at the vdW contact distance of 2.80 Å between the closest hydrogen and gold atoms. The molecular arrangement depicted in Fig. 1 is in fact drawn to scale.

As in  the equipotential metallic surface is treated by means of the image charges. For each atomic partial charge and induced dipole we assign and place an opposite charge and dipole in mirror positions when computing potentials and potential gradients. The charges and dipoles within a molecule do not contribute to the potentials and potential gradients that act on that molecule. However, the images of these charges and dipoles do contribute. For example, a single ion exhibits no polarization energy in gas phase, while it acquires polarization energy when placed near the metal surface.

As in bulk calculations, we solve the polarization problem for slabs in two steps. We first consider a neutral film of $`N`$ molecular layers with no ions and with translational symmetry in two dimensions. We solve the self-consistent equations and find ground-state charges $`\overline{\rho }_i^{kj}`$ and dipoles $`\overline{𝝁}_i^{kj}`$, which are different for each molecular layer $`k=1,\mathrm{},N`$ and for each type $`j=1,2`$ of molecule in the unit cell. While the slab is infinite in two dimensions, the number of variables is finite due to translational symmetry.

We then replace one molecule in the surface molecular layer by an ion and solve for the difference $`\delta \rho _i^a=\rho _i^a\overline{\rho }_i^{kj}`$ and $`\delta 𝝁_i^a=𝝁_i^a\overline{𝝁}_i^{kj}`$ for every molecule $`a`$ in the slab. Since the slab is infinite, we define clusters whose shapes are pillboxes with variable radius $`R`$ parallel to the metal and fixed thickness to include $`N`$ molecular layers and their images. We then relax the charge distribution of all $`M`$ molecules whose centers are in the pillbox. We assume $`\delta \rho _i^a=\delta 𝝁_i^a=0`$ for the molecules $`a`$ outside the pillbox, which means that their charge distributions are not relaxed. They contain, however, partial charges and dipoles $`\overline{\rho }_i^{kj}`$ and $`\overline{𝝁}_i^{kj}`$, which contribute to the total polarization energy of the film.

### B Surface and interface polarization

The polarization energy for the ion appears, as in the bulk, as a finite difference of two infinite energies: for the infinite slab with and without the ion. The expression for the energy that takes care of this cancellation is the same as Eq. (27) in , but with the potentials and potential gradients containing contributions from image charges and dipoles. We perform calculations for finite pillbox clusters of increasing $`M`$, and extrapolate to $`M\mathrm{}`$.

Fig. 4 shows the convergence of polarization energies for cations (upper panel) and anions (lower panel), as well as the extrapolated values for pentacene films one to ten molecular layers thick. Two lines for each film thickness correspond to two types of molecules in the unit cell. Convergence with increasing pillbox radius $`R`$ is shown as $`1/M_l`$, where $`M_l=M/N`$ is the average number of pentacenes per layer. The computational effort increases with $`N`$. The largest $`N=10`$ clusters in Fig. 3 contain $`M=2806`$ molecules, each containing 36 atoms, and their images.

The combined extrapolated values $`P^\mathrm{S}=P_+^\mathrm{S}+P_{}^\mathrm{S}`$ are also listed in Table III and plotted vs. $`1/N`$ in Fig. 5. The monolayer ($`N=1`$) and bulk values are almost the same. Image charges and dipoles at vdW separation are sufficiently more polarizable to compensate for a vacuum on the other side. $`P^\mathrm{S}(N)`$ decreases with increasing $`N`$ as less polarizable organic layers are introduced between the surface and the metal. In the limit of a thick film, or a free pentacene $`ab`$ surface, we estimate $`P^\mathrm{S}P0.23`$ eV.

Table III. Variation in polarization energy at the surface ($`P^\mathrm{S}`$) and at the metal/organic interface ($`P^\mathrm{M}`$) in an N-layer thick pentacene film on a metallic substrate. The values are reported for anion at position (1) and cation at position (2) in the unit cell, which correspond to the lowest energy state.

Layers, $`N`$  $`(P^\mathrm{S}P)`$, meV  $`(P^\mathrm{M}P)`$, meV    1 $`6`$ $`6`$  2 $`40`$ $`113`$  3 $`71`$ $`125`$  5 $`108`$ $`129`$  10 $`166`$ $`130`$  thick film $`227`$ $`130`$

In Table III the interface values $`P^\mathrm{M}=P_+^\mathrm{M}+P_{}^\mathrm{M}`$ refer to the polarization energy of an ion next to the metal in a film of $`N`$ molecular layers. The biggest increase of 100 meV at $`N=2`$ comes from the first pentacene overlayer. Additional layers produce small changes and the interface value in a thick films is $`P^\mathrm{M}P=0.13`$ eV. The hole and electron components of $`P^\mathrm{M}`$ are relevant for matching transport levels to Fermi energies for facile injection of carriers. Such optimization for organic devices is largely empirical at present, often without any consideration of polarization. The variation of $`P^\mathrm{S}(N)`$ or $`P^\mathrm{M}(N)`$ with film thickness is less in pentacene than in PTCDA . Pentacene layers are several times thicker, with the long axis normal to the metal, and a single layer already provides an effective polarization environment for the charge.

Fig. 4 Convergence of polarization energies of a cation (upper panel) and anion (lower panel) in the outermost surface layer of pentacene film on gold vs. the average number $`M_l=M/N`$ of molecules in each molecular layer ($`M`$ molecules in $`N`$ layers in the pillbox). Open and solid symbols correspond to the molecules of type (1) and (2) respectively.

Fig. 5 Variation in polarization energy at the surface ($`P^\mathrm{S}(N)P`$) and at the metal/organic interface ($`P^\mathrm{M}(N)P`$) for an $`N`$-layer thick pentacene film on a metal substrate.

### C Molecular potential

We have already noted the competition between image charges and induced dipoles in films of molecules with large polarizability normal to the surface. The self-consistent potential $`\varphi ^a(𝐫)`$ at molecule $`a`$ is due to the polarized densities $`\rho ^b(𝐫)`$ of all other molecules and their images. A single pentacene cation in the surface layer induces an image anion in the metal. Since the interaction is attractive, the cation’s $`\varphi ^{cat}(𝐫)`$ increase monotonically with the distance from the surface. Similarly, the $`\varphi ^{an}(𝐫)`$ of a single anion on the surface decrease monotonically with distance from the surface.

The situation in films is quite different due to the surrounding neutral molecules. Since the large pentacene polarizability is normal to the surface, the induced moments due to the image anion are parallel. The repulsion of parallel dipoles and image dipoles is relieved by redistributing charge towards the middle of the molecular layer. Now the $`\varphi ^{cat}(𝐫)`$ has a minimum at atoms near the center, while the $`\varphi ^{an}(𝐫)`$ has a maximum near the center. The actual $`\varphi ^{ion}(𝐫)`$ in the limit of zero overlap is the solutions of the self-consistent equations.

Figure 6 shows the self-consistent potential $`\varphi ^{cat}(𝐫)`$ for $`𝐫`$ along the axis of a cation in a pentacene monolayer on metal. The gradients yield electric fields $`\pm 10^7`$ V/cm that reverse over a single molecule, directly confirming that fields are strong and non-uniform. The variation of $`\varphi ^{cat}(𝐫)`$ with $`𝐫`$ shows the importance of repulsion between induced dipoles. Image charges account for more negative $`\varphi ^{cat}(𝐫)`$ close to metal. Similar variations of self-consistent potentials are found at ions (1) and at neighboring neutral molecules. The strong shielding of a single layer of pentacene follows from the weak dependence on $`N`$. Similar curves (not shown) for $`N=2`$, 3, and 5 are essentially parallel to the $`N=1`$ curve in Fig. 6, but displaced to less negative (positive) potential for cations (anions).

Fig. 6 Self-consistent potential for a cation created by the polarization field of neutral molecules in the lattice and the image charge distribution in the metal. Data shown for a cation in a monolayer on metal substrate, placed in the position (2) in the unit cell. The sketch of the molecule shows actual position with respect to the metal.

## V Discussion

The bulk polarization energy can be estimated simply in terms of molecular size, $`P_\pm =(e^2/2R)(11/\kappa )`$, for a sphere of radius $`R`$ in an isotropic dielectric medium with $`\kappa 3`$ for organic molecular crystals. Such rough values rationalize much data but preclude accurate positioning of transport or CT states. The effective pentacene radius is 14% greater than that of anthracene, while the polarization energy decreases by only 9%. On the other hand, the radii of pentacene and PTCDA molecules are the same within 3%, but the polarization energy differs by 10%. The greatest contrast to PTCDA films , is that the polarization energy $`P^\mathrm{S}`$ is within 10% of the bulk value in pentacene, compared to $`P^\mathrm{S}/P=0.77`$ in PTCDA.

Direct relations between gas-phase and crystalline energies have long been sought using molecular exciton theory . Accurate evaluation of electronic polarization permits sharper comparisons. There is considerable scatter in reported values for both the gas and solid, as noted in Sections III(A) and III(D) for the ionization potential and electron affinity of pentacene. The transport gap is $`IAP`$ in the limit of zero intermolecular overlap. The experimental values for the pentacene transport gap range from 2.85 eV  to 2.4 eV suggested in , to 2.2 eV used in . This scatter was not “important” when polarization energies had to be estimated in the bulk and were not available at all at surfaces. With the ability to calculate polarization energies, it becomes possible to analyze and connect various sources of experimental data on such quantities as optical and transport gaps, CT binding energies, and work functions.

All molecular calculations, except for the gas-phase properties, were performed using experimental geometries from the crystal structure data . Slight variations in the geometry of the two molecules in the unit cell were consequently taken into account. While the variations in electronic energies are small (tens of meV), which is consistent with the adopted approximation, we believe that crystalline geometries are preferable to gas-phase geometries when vertical electronic transitions are considered in the condensed phase.

We note that INDO/S greatly overestimates the in-plane polarizabilities of ions in Table I. The huge difference of $`\alpha ^\mathrm{C}`$ between cations 1 and 2 is clearly unphysical in view of small distortions in the crystal, and it rationalizes why convergence was particularly difficult for cation 1. The UHF INDO/S approach fails for ions. Radical ions of large conjugated molecules are expected to be more polarizable than the neutral molecules, in line with the B3LYP entries in Table I. While the idea is that $`\stackrel{~}{\alpha }`$ constitutes a small correction, and it is not so for ions, we expect this to have a negligible impact on the results. The reason is that the largest part of the polarization energy comes not from the ion itself but from the polarization of many neutral molecules surrounding it. Moreover, the molecular $`\alpha `$ is kept no matter how it is partitioned between charge redistribution ($`\alpha ^\mathrm{C}`$) and induced dipoles. In the future we may choose to use $`\mathrm{\Pi }_{ij}`$ of the neutral molecule for the ions, and to absorb the polarizability difference in $`\stackrel{~}{\alpha }`$. This requires separate B3LYP calculation of $`\alpha `$ for the ions, as done here, but not in .

Small atomic charges in acenes are related to approximate electron-hole symmetry. The crystal potential $`\varphi _i=\varphi (𝐫_i)`$ at atom $`i`$ due to INDO/S charges is consequently small, as shown by $`1`$ meV per molecule contributions to the sublimation energy. Since the potential due to the B3LYP or any other gas-phase charge distribution can readily be computed at the positions of other atoms in the crystal, we can obtain first-order corrections to polarization energies in general  by combining the best available gas-phase potential with self-consistent charges and dipoles based in INDO/S calculations. The quadrupole moments of acenes lead to opposite shifts of $`P_+`$ and $`P_{}`$, as discussed previously for submolecules . Since such corrections cancel in $`P=P_++P_{}`$ in anthracene , we expect them to be small in pentacene. They are not negligible in PTCDA or conjugated molecules that contain heteroatoms and hence partial charges.

Our values for the the transport gap are likely an overestimate, due to the bandwidth effects for electrons and holes, and the correction can be introduced as $`E_t=IAPW`$. The bandwidth $`W`$ in molecular crystals at room temperature is thought to be of the order of few tens of meV or less, which justifies our zero-overlap approximation. At lower temperatures, however, the bandwidths may be large in certain directions. It has been suggested that the bandwidth in pentacene (at zero temperature) may be anomalously large and reach 500 meV . If we estimate $`W0.5`$ eV, we get $`E_t=2.71`$ eV with the experimental charge gap, while the calculated charge gap leads to the value 2.35 eV, which is too small.

When bandwidths are comparable to polarization energies, the zero-overlap approximation cannot be justified, and we face a challenging problem of combining the “localized” treatment of polarization effects with the “delocalized” language of band theory. The self-consistent polarization approach can provide a proper “zero-order” approximation for the attempts to treat bandwidth effects as perturbation.

We note that our self-consistent approach has the ability to distinguish between the electronic polarization of cations and anions, since the quantum-mechanical structure of molecular response is used to obtain the redistribution of charge. We can also estimate variations of polarization caused by slight changes in molecular geometry between crystal and gas-phase, and especially between inequivalent molecules in the unit cell. These variations appear to be of the order of 20—40 meV, which is on the threshold of the accuracy, but such variations are probably correct on the order of magnitude, and may be important in the analysis of experimental data.

In summary, we have applied the recently developed self-consistent approach based on evaluating charge redistribution in organic molecules to the problem of electronic polarization in pentacene molecular crystals. The power of the approach comes from the combination of semiempirical treatment of charge distribution by means of atom-atom polarizability tensors with accurate ab initio gas-phase calculations of molecular polarizabilities. As a result, we are now able to calculate polarization energies with accuracy on the order of 0.1 eV or better.

ACKNOWLEDGMENTS

It is a pleasure to thank A. Kahn and R.A. Pascal, Jr. for stimulating discussions. We gratefully acknowledge support for work at Rutgers through the Office of Naval Research, grant No. NO00014-01-1-1061, and at Princeton through the MRSEC program under DMR-9400632.


format: text

text:
# Tracking Black Holes in Numerical Relativity

## I Introduction

The idea of a totally collapsed gravitational source, from which nothing — not even light — can escape, is an old one, dating back at least to the work of Laplace (and others) in the eighteenth century . Since those original considerations, research has uncovered a wealth of understanding regarding the physics of black holes. In general, the work has followed two major routes, with an ever - narrowing gap between the approaches. Along one direction, black holes are studied as mathematical solutions in a given theory of gravity including Newtonian theory, post-Newtonian metric theories of gravity, Einsteinian gravity, as well as semi-classical and (more recently) various attempts at full quantum theories of gravity. The second direction considers black holes in the astrophysical and astronomical contexts. Both fields of research have seen increasing activity over the years and have made startling discoveries concerning the physics of black holes; including the gravitational collapse theory of Oppenheimer and Snyder , proofs of uniqueness and stability of black holes in Einstein’s general relativity , thermodynamic properties of black holes including the three laws of black hole mechanics  and mechanisms for black hole radiation and evaporation , discovery of critical phenomena in black hole formation , experimental identification of both astrophysical black hole sources themselves  and their event horizons , experimental signatures for super - massive black holes in galactic centers , and experimental bounds on the distribution and spectrum of black hole sources and collisions .

Numerical relativity, which addresses computational solution of Einstein’s equation , is an active participant in both the mathematical and the astrophysical programmes of research. With the advent of the Laser Interferometric Gravitational Observatory and other similar efforts, considerable attention has focused on the generic binary black hole coalescence (BBHC) problem, expected as the strongest sources for gravitational waves . As discussed in , this is a very difficult and interesting problem that can only successfully be addressed in the computational domain. Consequently, the binary black hole coalescence problem has become an active subject of numerical relativity. One particular problem in the computational domain is the computational definition, detection and tracking of the black hole event horizon itself , , , , , .

The present work completely solves the problem of numerically tracking black hole event horizons. The solution is complete in the sense that a single method is presented such that any one implementation of the method can generically detect arbitrary numbers of black hole event horizons undergoing arbitrarily strong gravitational interactions. For example, using a single computational code of our method we analyse both single black holes and black holes undergoing merger; and no special modifications of our code are required to handle these distinct dynamics.

The present article, describing our generic method for tracking black hole event horizons, is divided as follows. In section II the eikonal equation, the foundation of our method, is described in sufficient detail to be employed in an event horizon tracker. In particular, we focus on the signature behavior of a black hole event horizon in solutions of the eikonal equation. As usual in numerical work, there are a variety of possible implementations of the approach. In section III, several closely related systems of equations are presented and one particular system is singled out for consideration. The system chosen makes use of an explicit second order diffusion, or viscosity, term and we show in section III the relationship between solutions of the diffusive equations of motion and the continuum eikonal equation of interest. Section IV details extraction of the two dimensional sections of the event horizon for each time level of a numerical evolution. Such extraction is crucial for carrying out area, mass, and spin calculations. Section IV also shows the accuracy of our implementation by considering a parameter space survey of single Kerr black holes. Sections V - X describe in detail the first three dimensional application of our (or any) method to the numerical analysis of the asymmetric binary black hole coalescence problem. In particular we search for evidence of any nontrivial topology of the horizon immediately following merger. In distinction to the prediction of  we find a topologically spherical horizon to within the accuracy of our three dimensional mesh. Further, area analysis of the candidate numerical event horizon is carried out in conjunction with analysis of the black hole apparent horizons, which reveals both a time of merger much earlier then estimated using apparent horizons and mass energy estimates much larger those found using apparent horizon tracking methods. Our conclusions are presented in section XI.

## II The Eikonal

To begin, since the Lagrangian $`L=g_{ab}\dot{x}^a\dot{x}^b`$ of null geodesic motion has only kinetic terms it is equal to the associated Hamiltonian***We use early-latin indices to denote spacetime components $`a,b,c=0,1,2,3`$; mid-latin indices denote spatial components $`i,j,k=1,2,3`$.. Legendre transformation

$$H(\tau ,x^a,p_b)=\frac{dx^c}{d\tau }p_cL(\tau ,x^d,\frac{dx^e}{d\tau })$$
(1)
where

$$p_a\frac{L}{\left(\frac{dx^a}{d\tau }\right)},$$
(2)
sets $`L=H`$. The corresponding Hamiltonian equations are

$$\frac{dx^a}{d\tau }=\frac{H}{p_a}=2g^{ab}p_b$$
(3)
$$\frac{dp_b}{d\tau }=_bH=p_cp_d_bg^{cd}.$$
(4)
It is generic that the spacetime metric is independent of the affine parameter $`\tau `$. There is thus a first integral associated to geodesic motion. Making use of this property permits elimination of the affine parameter in favor of coordinate time $`t`$. To see this, it is convenient to adopt the ADM variables

$$g^{tt}=\frac{1}{\alpha ^2},g^{ti}=\frac{\beta ^i}{\alpha ^2},g^{ij}=\gamma ^{ij}\frac{\beta ^i\beta ^j}{\alpha ^2}.$$
(5)
Here $`\gamma ^{ij}`$ is the inverse spatial metric. Substituting the ADM variables directly into the Hamiltonian gives

$$H=p_ag^{ab}p_b=\frac{1}{\alpha ^2}p_t^2+\frac{2}{\alpha ^2}\beta ^ip_i+p_i\left(\gamma ^{ij}\frac{1}{\alpha ^2}\beta ^i\beta ^j\right)p_j.$$
(6)
Since the Hamiltonian is not explicitly $`\tau `$ dependent there is a constant of the motion

$$\frac{1}{\alpha ^2}p_t^2+\frac{2}{\alpha ^2}\beta ^ip_i+p_i\left(\gamma ^{ij}\frac{1}{\alpha ^2}\beta ^i\beta _j\right)p_j=\omega ^2.$$
(7)
For $`\omega ^2<0`$, $`\omega ^2=0`$, or $`\omega ^2>0`$ the motion is said to be timelike, null, or spacelike. Without loss of generality, assuming null geodesic motion, solution of (7) by ordinary algebra yields:

$$p_t=\beta ^ip_i\pm \alpha \sqrt{p_i\gamma ^{ij}p_j}.$$
(8)
With this result, the Hamiltonian can be explicitly factored:

$$H=H_+H_{}=0$$
(9)
where

$$H_\pm =p_t\beta ^ip_i\pm \alpha \sqrt{p_i\gamma ^{ij}p_j}=0.$$
(10)
In the case of either root Hamilton’s canonical equations become

$$\frac{dt}{d\tau }=H_{}\frac{H_\pm }{p_t}=H_{}$$
(11)
$$\frac{dx^i}{d\tau }=H_{}\frac{H_\pm }{p_i}=H_{}\left(\beta ^i\pm \alpha \frac{\gamma ^{ij}p_j}{\sqrt{p_k\gamma ^{kl}p_l}}\right)$$
(12)
$$\frac{dp_t}{d\tau }=H_{}\frac{H_\pm }{t}$$
(13)
$$\frac{dp_i}{d\tau }=H_{}\frac{H_\pm }{x^i}.$$
(14)
According to equation (11), $`\tau `$ can be eliminated in favor of $`t`$. The system of equations becomes simply the equations (12) and (14) with $`t`$ written in place of $`\tau `$ and the factor $`H_{}`$ cancelled everywhere on the right hand side. Explicitly,

$$\frac{dx^i}{dt}=\beta ^i\pm \alpha \frac{\gamma ^{ij}p_j}{\sqrt{p_k\gamma ^{kl}p_l}}$$
(15)
and

$$\frac{dp_i}{dt}=_i\left(\beta ^jp_j\pm \alpha \sqrt{p_j\gamma ^{jk}p_j}\right).$$
(16)
The eikonal, corresponding to the Hamiltonian of equation (6),

$$_aSg^{ab}_bS=0$$
(17)
can be factored similarly. To do so it is a simple matter of making the replacements $`p_t_tS`$, $`p_i_iS`$ in (8) to find the following symmetric hyberbolic partial differential equation

$$_tS=\beta ^i_iS\pm \alpha \sqrt{_iS\gamma ^{ij}_jS}\overline{H}.$$
(18)
Note that a bar is introduced here to distinguish the Hamiltonian used here from the Hamiltonian used in (6). This equation is also used in the method of , although in that work the equation is further reduced to consider the case of a single null surface.

The right hand side of this result is homogeneous of degree one in $`_iS`$. The characteristic curves along which the level sets $`\mathrm{\Gamma }`$ of $`S`$ are propagated, are then

$$\dot{x}^i=\beta ^i\pm \alpha \frac{^iS}{\sqrt{_iS\gamma ^{ij}_jS}}\frac{\overline{H}}{\left(_iS\right)}$$
(19)
$$_i\dot{S}=_i\overline{H}(t,x^j,_jS),$$
(20)
which are the null geodesics of equations (15) and (16). Immediately, the integral curves of the gradients of $`S`$ and $`\mathrm{\Gamma }`$ are also the null geodesics:

$$\frac{dx^i\left(\lambda \right)}{d\lambda }=g^{ia}p_a=g^{ia}_aS=^iS(\lambda ,x^j\left(\lambda \right)).$$
(21)
Hereafter, the bar on $`\overline{H}`$ is dropped with the understanding that the Hamiltonian considered is that of equation (18). This result establishes that the eikonal is technically a Riemann invariant of the null geodesics, a fact that proves useful in establishing the signature of a black hole event horizon in solutions of the eikonal equation. More specifically, since the eikonal is the canonical generator of null geodesics, it can be employed in analysis of black hole event horizons, which are by definition generated by null geodesics having no future end points. To proceed in this manner the equation (18) is employed in an initial value problem and then surveyed for signature features of black hole event horizons.

A black hole event horizon is generated by a congruence of outgoing — but future asymptotically nonexpanding — null geodesics. The scope of the surveys of the eikonal equation that are required to identify a black hole event horizon is then restricted to the space of all outgoing null surfaces. These surveys are greatly reduced by the fact that solutions of the eikonal are categorized by topologically equivalent solutions. To see this, note that

$$_i\psi \left(S\right)=\frac{\psi }{S}_iS=\lambda \left(S\right)_iS,$$
(22)
and

$$_t\psi \left(S\right)=\frac{\psi }{S}_tS=\lambda \left(S\right)_tS.$$
(23)
By homogeneity of the right hand side of (18), if $`S`$ is a solution then $`\psi \left(S\right)`$ is also solution. Thus smoothly related initial data $`S_0S_0^{}=\psi \left(S_0\right)`$ have smoothly related solutions. This feature alleviates the need for surveying over smoothly related initial data.

A further reduction of the scope of solution surveys is provided by the equivalence of ingoing and outgoing solutions under time reversal. Propagation of data for $`S`$ describing an ingoing or outgoing null surface is accomplished by specification of: 1. A definition of the direction of time, 2. A choice of $`\alpha `$ and $`\beta ^i`$, and 3. A choice of the root. With these choices specified, data is then uniquely partitioned into an ingoing type and an outgoing type with the distinction being the gradient of $`S`$.

With the dynamics of any outgoing null surface (including the event horizon) specified by (18), and the scope of solution surveys categorized into topological classes, the task remains of identifying the event horizon within this restricted space of outgoing null surfaces. To do so, the results and approach of  are adopted here and modified to include the eikonal equation as discussed above. The result of this approach is a signature feature of black hole event horizons in the eikonal equation.

An event horizon of a black hole is by definition a critical outgoing surface when tracked into the future. Let $`𝒫`$ be a point interior to the horizon and $`𝒬`$ be a point exterior to the horizon such that $`𝒫`$ and $`𝒬`$ lie on characteristic curves of the eikonal $`\gamma _𝒫`$ and $`\gamma _𝒬`$. At arbitrarily early times let $`\gamma _𝒫`$ and $`\gamma _𝒬`$ pass arbitrarily closely to a point $``$ that lies on the horizon $`\mathrm{\Gamma }`$. Since $`S`$ is a Riemann invariant, at arbitrarily early times the jump of the eikonal at $``$ becomes $`\left[\left[S\left(\right)\right]\right]S_{}\left(\right)S_+\left(\right)S\left(𝒫\right)S\left(𝒬\right)0`$. In the computational domain, where the resolution is finite, this discontinuity will appear generically in finite time. As such, an approximation of the event horizon will appear numerically as the formation of a jump discontinuity in the eikonal for outgoing data that is propagated into the past. This is the numerical signature of black hole event horizons in the eikonal’s solutions.

## III Numerical Methods

Analysis of the continuum properties of the eikonal equation as a Hamilton Jacobi equation identifies three closely related approaches to tracking black hole event horizons:

* System I: The null geodesic equations

$$S(x^j,t)=S_0\left(x^k\right)𝑑tH(t,x^i,p_i)$$
(24)
where the integral is evaluated along the solutions of

$`\dot{x}^i`$ $`=`$ $`\beta ^i\pm \alpha {\displaystyle \frac{p^i}{\sqrt{p_i\gamma ^{ij}p_j}}},`$ (25)
$`\dot{p}_i`$ $`=`$ $`_iH.`$ (26)
* System II: The eikonal

$$_tS=H(t,x^j,_jS)$$
(28)
* System III: The flux conservative form

$$S(x^i,t)=𝑑x^ip_i(x^j,t)$$
(29)
$$_tp_i(x^i,t)+_iH(t,x^j,p_j)=0.$$
(30)

In each of the above three systems of equations the Hamiltonian is given by (18).

In each of the above cases, the structure of the dynamical equations provide certain advantages. For example, in each case the symplectic structure can be employed to identify numerical loss of accuracy in a similar manner to some modern numerical schemes used in Hamiltonian dynamics. Further, as a flux conservative system, high resolution methods from computational fluid dynamics can be applied directly to the third system. In following sections of this article, the second system is considered since this system of equations yields an expedient implementation that is sufficiently accurate for our purposes.

Singular behavior on the eikonal is not specific to only the event horizon and instead, as described in detail by Arnold and Newman , , , the eikonal is known to generically break down on caustic and other sets. Special numerical methods are then required to handle the generic singular behavior of the eikonal; we make use of an explicit viscosity term. In the continuum, addition of our form of numerical viscosity at the level of the finite difference approximation corresponds to replacing the evolution of

$$_tS=H(t,x^i,_jS)$$
(31)
with evolution of the equation

$$_t\psi =ϵ^2^2\psi H(t,x^i,_j\psi )$$
(32)
where $`ϵ`$ is a small quantity which we call the viscosity and $`^2`$ denotes any second order, linear derivative operator. There is a well defined sense in which the solutions $`S`$ relate to the solutions $`\psi `$; it is simply given (when the solutions $`S`$ exist) by the WKB transformation

$$\psi (x^i,t)=\underset{n}{}a_n(x^i,t)ϵ^n\mathrm{exp}\left(S/ϵ\right)A\mathrm{exp}(S/ϵ).$$
(33)
To see the explicit relationship between the solutions $`\psi `$ and the solutions $`S`$ note that

$$_t\psi =\frac{\psi }{ϵ}\left(_tS+ϵ_t\mathrm{log}A\right),$$
(34)
$$_i\psi =\frac{\psi }{ϵ}\left(_iS+ϵ_i\mathrm{log}A\right),$$
(35)
$$^2\psi =\frac{\psi }{ϵ^2}\left(S+ϵ\mathrm{log}A\right)^2+\frac{\psi }{ϵ}\left(^2S+ϵ^2\mathrm{log}A\right).$$
(36)
Assuming that the Hamiltonian is homogeneous of degree one in momentum, and making use of perturbation theory:

$$H(x^i,_j\psi )=\frac{\psi }{ϵ}\left(H(x^i,_jS)+ϵH_1(t,x^i,_j\mathrm{log}A)\right).$$
(37)
Substituting these results into (32) and cancelling an overall factor of $`\psi /ϵ`$ gives at lowest, and first order:

$$_tS=H(t,x^i,_jS).$$
(38)
$`ϵ_t\mathrm{log}A`$ $`=`$ $`ϵH_1(t,x^i,_j\mathrm{log}A)`$ (39)
$`+`$ $`ϵ\left(S+ϵ\mathrm{log}A\right)^2+\mathrm{}`$ (40)
Here $`H_1`$ is the first order linear Hamiltonian obtained from perturbation theory. At zeroth order, $`S`$ then satisfies the eikonal equation, while the first order correction is the linear result of first order perturbation theory and expresses the evolution of $`A`$. Higher order results can be found similarly. Again, these results hold only where solutions of $`S`$ exist; that is, away from discontinuities or other solution singularities in $`S`$ and its derivatives.

The above analysis is a slight modification of usual WKB expansion and expresses the solutions $`\psi `$ in terms of the solutions $`S`$. Such solutions can be inverted using the conventional method of series inversion. To invert the WKB expression, and express the solutions of the eikonal in terms of the solution of the parabolic equation, assume first that the solution of the parabolic equation (32) is given as $`\psi `$ (numerically or otherwise). Writing

$$S=\mathrm{log}\left(\frac{\psi }{\overline{A}}\right),$$
(42)
the procedure is then analogous to that of the WKB expansion: $`\overline{A}`$ is constructed as an asymptotic power series with coefficients depending on $`\psi `$ alone. The result of such an analysis is simply:

$$S(x,t)=\mathrm{log}\psi (x,t)+\overline{ϵ}𝑑t^{}\frac{^2\psi (\overline{x},t)}{\psi }$$
(43)
where the integral is evaluated along

$$\frac{d\overline{x}^i}{dt}=\frac{H(t,\overline{x}^j,_{\overline{x}^k}\mathrm{log}\psi )}{p_i}$$
(44)
and $`\psi `$ is provided independently by (numerical) solution of

$$_t\psi =\overline{ϵ}^2\psi H(t,x^i,_j\psi ).$$
(45)
One advantage of an explicit second order viscosity term is a simple procedure for reducing the error of viscosity solutions by one order of the viscosity $`ϵ`$. The limit $`ϵ0`$ corresponds to a continuous family of zeroth order solutions $`\varphi _ϵ`$; where

$$\varphi _ϵ=S+\eta _ϵ$$
(46)
and $`\eta _ϵ=𝒪\left(ϵ\right)`$. Note that each solution $`\varphi _ϵ`$ is obtained from an analysis similar to that following equation (42), although obtained by neglecting the logarithm.

Given two viscosity solutions $`\varphi _ϵ`$ and $`\varphi _{2ϵ}`$ it is then possible to construct a third improved viscosity solution $`\varphi _I`$ that is accurate to $`𝒪\left(ϵ^2\right)`$. To see this, consider the combination

$$\varphi _I\varphi _ϵ+\left(\varphi _ϵ\varphi _{2ϵ}\right).$$
(47)
$`\varphi _I`$ is accurate to $`𝒪\left(ϵ^2\right)`$ since

$$\varphi _ϵ\varphi _{2ϵ}=\eta _ϵ\eta _{2ϵ}\eta _ϵ2\eta _ϵ=\eta _ϵ.$$
(48)
To make use of these improved viscosity solutions let $`h`$ denote the resolution of the numerical mesh. Any second order finite difference approximation will have then

$$\varphi =\widehat{\varphi }+𝒪\left(h^2\right),$$
(49)
where $`\varphi `$ is the continuum solution and $`\widehat{\varphi }`$ denotes its finite difference approximation.\[We will use a hat throughout the paper to denote the discrete approximation to a continuum object.\] Similarly,

$$S=\widehat{S}+𝒪\left(h^2\right).$$
(50)
Using

$$S=\varphi +𝒪\left(ϵ\right)=\varphi _I+𝒪\left(ϵ^2\right)$$
(51)
gives

$$\widehat{S}=\widehat{\varphi }+𝒪\left(h^2\right)+𝒪\left(ϵ\right)=\widehat{\varphi _I}+𝒪\left(h^2\right)+𝒪\left(ϵ^2\right).$$
(52)
## IV Level Set Extraction

Of crucial interest in the binary black hole coalescence problem are the areas of sections of the black hole event horizon. To find these sections, at any given time level any of the level sets, say $`S=0`$, can be extracted to obtain the surface $`\mathrm{\Gamma }`$, a two dimensional section of the corresponding null surface. This problem of extraction is an inverse problem, since it requires that points $`(x,y,z)`$ are found such that $`S(x,y,z)=0`$. To accomplish this inversion, we find that an ordinary bisection method is sufficient for use with an ordinary second order interpolation scheme. In this method it is assumed that the surface $`\mathrm{\Gamma }`$ can be expressed in spherical coordinates $`(\theta ,\varphi ,u(\theta ,\varphi ))`$, where $`r=u(\theta ,\varphi )`$ is the surface function for a given center $`c^i`$ contained within the surface $`\mathrm{\Gamma }`$. Given a choice for the center, the radial function for $`S\left(\mathrm{\Gamma }\right)=S_o=0`$ can then be approximated via the interpolation and bisection.

To establish the accuracy of our implementation, including the routines that accomplish extraction of the level sets and the accuracy of the viscosity term, we consider stationary, spinning black holes. This case is completely described by the Kerr - Newmann family of axisymmetric solutions of Einstein’s equation. It is convenient to make use of the Kerr-Schild form for the metric  since this form is used in our binary black hole evolution code as well as in our solution of the initial value problem for setting initial data for the evolution. Specifically, the metric in the Kerr-Schild form is

$$g^{ab}=\eta ^{ab}2Hl^al^b.$$
(53)
Here $`\eta _{ab}=\mathrm{diag}(1,1,1,1)`$ is Minkowski’s metric, $`H`$ is a space time scalar, and $`l^a`$ is an ingoing null vector with respect to both the Minkowski and full metric. The Kerr solution is the two parameter family of solutions such that

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

$$l^t=1,$$
(55)
$$l^x=\frac{rx+ay}{r^2+a^2},$$
(56)
$$l^y=\frac{ryax}{r^2+a^2},$$
(57)
$$l^z=\frac{z}{r},$$
(58)
$$r^2=\frac{1}{2}\left(\rho ^2a^2\right)+\sqrt{\frac{1}{4}\left(\rho ^2a^2\right)+a^2z^2}$$
(59)
where

$$\rho ^2=x^2+y^2+z^2.$$
(60)
Finally, the event horizon for the Kerr black hole is located on the ellipsoid $`r=r_+=r(x,y,z)`$ where

$$r_+=M+\sqrt{Ma}.$$
(61)
In figures (1) -(4) we show the evolution backward in time of the eikonal equation (followed by extraction of the $`S=0`$ surface). Errors can arise both in the evolution and in the extraction of the surface. Figure (1) shows the percent errors in the extracted areas for a nonspinning black holes with mass $`M=2`$, $`a=0`$ in a survey over the viscosity parameter. Figure (2) shows a similar study but for the L2 norm of the truncation error in the function $`r_+`$. Note that the function $`r_+`$ is defined for every point of the horizon in the continuum $`r_+=r_+(x,y,z)`$. In the computational domain $`r_+`$ then takes a discrete form $`\widehat{r}_+=\widehat{r}_+(i,j,k)`$, where the integers $`i,j,k`$ span the numerical mesh. The truncation error $`e_{r_+}`$ is then

$$e_{r_+}(i,j,k)=r_+(x,y,z)\widehat{r}_+(i,j,k)$$
(62)
where is it understood that both $`r_+`$ and $`\widehat{r}_+`$ are evaluated at the same point. Bith Figures(1) and (2)show that viscosity parameter $`h^2/8`$ adequately captures the horizon location. While not perfect, it will suffice for the short term horizon tracking reported here.

However these results also suggest that in vanishing viscosity the percentage error in the calculated area is reduced toward a bias. This bias is partly associated to the finite resolution of the computational mesh and partly to accuracy of the extraction routine. These figures were generated using a three dimensional computational domain of $`N^3`$ points with $`N=121`$. The outer boundaries are located at $`[15M,+15M]`$ in the $`x,y,z`$ directions. The resolution of the finite difference mesh for these results is then $`h=M/4`$. Also, a Courant - Friedrichs - Lewy factor of $`\lambda =1/4`$ with an iterated Crank Nicholson scheme  was used as the finite difference approximation of the evolution of the eikonal equation (18). Neumann boundary conditions $`_iS=0`$ on the outer boundary are found to be generically sufficient conditions for stability of the method. The philosophy here is that the primary interest is deep within the bulk of the computational domain where the event horizon of the black hole is located. The outer boundary is then treated only to the degree that the evolution of the interior region remains stable. Further, the interior of the black hole is excised from the computational domain in a sphere of radius $`r=r_0+2dx`$, where $`r_0`$ denotes the radius of the Kerr - Newman ring curvature singularity. Finally, the discrete surface $`\widehat{\mathrm{\Gamma }}`$ was constructed using $`m^2`$ points with $`m=100`$. At such a resolution of the extracted surface any errors in the area must be attributed to all of: the viscosity parameter, the extraction routine, and the resolution of the underlying three dimensional grid.

Figures (3) and (4) show the percent errors in the extracted areas as well as the L2 Norm of the truncation error in the function $`r_+`$ in a survey over the angular momentum parameter $`a`$. Here the viscosity is $`ϵ=h^2`$. These figures both show the evolution from a single null sphere that is completely exterior to the horizon and propagated into the past. Since the event horizon of a spinning black hole is elliptical in its geometry, the spherical data we have chosen for $`t_{}=0`$ corresponds to a percent error in the area and $`r_+`$ that varies with the spin parameter $`a/M`$ at $`t_{}=0`$, explaining why the curves do not intersect at $`t_{}=0`$. (Where appropriate we append a $``$ to $`t`$, thus: $`t_{}`$, indicating evolution into the past; we also sometimes use $`t_{}`$ to emphasize that we mean the forward evolving, usual, time $`t`$. Thus $`t_{}=0`$ corresponds to the late time at which we begin to integrate into the past.)

However, the errors should converge to a constant, which is evident in each of the curves with $`a<0.75`$ in figures (3) and (4). For $`a=0.75`$ the curve does not converge, and we expect that different choice of initial data will exhibit convergence. According to these results, the viscosity level set method does indeed accurately and robustly detect the distorted outermost event horizons of spinning black holes at least when $`a<0.75`$. Note that we study both the accuracy of $`r_+`$ and the accuracy of the calculated areas since the calculations seperately and together establish the accuracy of our area calculation and of our detection of the event horizon.

## V Asymmetric Binary Black Hole Coalescence

Analysis of the event horizon for the binary black hole coalescence problem in the case of head on collision has been considered in detail in . The problem of the event horizon for asymmetric, that is off axis, collision has only been considered analytically  and prior to this work no results for numerically generated sources have been analyzed. Numerical evolution and analysis of an asymmetric binary black hole system was studied in , but at that time the question of the event horizon was not considered.

In this section the results of the previous sections are applied to the first completely numerical analysis of the event horizon for the case of asymmetric collision.

To begin, consider two black holes of mass $`M=1`$ with aligned spins in the positive $`z`$ direction of $`a/M=1/2`$. The computational domain is a grid of $`N^3`$ points with $`N=121`$. The outer boundary is located at $`\pm 15M`$ and the holes are initially positioned at $`(x,y,z)=(+6,+2,0)`$, and $`(x,y,z)=(6,2,0)`$. This computational domain is identical to the mesh used in the previous section to analyze the percent error in area calculations of surfaces extracted from the level set method. The percent error in the calculation of the area of sections of the horizon should have a magnitude of about $`45\%`$. Further, as an order of magnitude estimate, in a flat spatial geometry; e.g., in a Newtonian spacetime, the initial separation $`s`$ of the black hole centers would be $`s=\sqrt{12^2+4^2}12.64`$. This would seem to be an ample initial separation to guarantee that the initial data corresponds to two distinct black holes. However, in the case that the holes are nonspinning, each will have a spherical event horizon of radius $`r=2M`$. Assuming only marginal distortion of the nonspinning event horizons due to spin effects, which is an approximation that is justified by the properties of spin $`a/M=1/2`$ black holes, the nearest separation between the two sections of the black hole event horizon is then approximately $`s_{min}=12.64M4M8.64M`$. Again, this approximation assumes a flat underlying geometry; and so can be considered only as an order of magnitude estimate. These initial data then appear to correspond to a separation of approximately two nonspinning black hole diameters between the surfaces of section of the black holes.

The holes are boosted along the $`x`$ direction with speeds of $`\pm c/2`$. This boost lengthens the nearest separation $`s_{min}`$ of the holes due to Lorentz contraction of the horizons.\[In this coordinate system the horizons undergo contraction in the direction of motion. For a single hole, the area of the horizon does not change under this boost.\] The nearest coordinate separation $`s_{min}`$ between the holes is then expected to lie in the range $`8.64<s_{min}<12.64`$.

The numerical evolution of this collision process is carried out for approximately $`10M`$ of run time with a Courant factor of $`\lambda =dt/dx=1/4`$. The code is the Texas black hole evolution code, a derivative of the Agave code . Apparent horizon finders  locate two distinct apparent horizons of area $`A50M`$ for the initial data and continue to do so until $`t=8M`$, when only a single apparent horizon of area $`A200M`$ can be located. This single apparent horizon persists until approximately $`10M`$, beyond which instability effects, stemming from the outer boundary and the excision boundary, swamp the solution.

## VI Data For the Eikonal

While the run length of this asymmetric collision data is relatively small in units of the black hole masses, the problem of detecting the associated black hole event horizon, or horizons, is not a small computational problem. For example, to analyze this data requires analysis of the lapse $`\alpha `$, the shift vector $`\beta ^i`$, and the three metric $`\gamma ^{ij}`$. By symmetry of the metric, $`\gamma ^{ij}=\gamma ^{ji}`$, there are only 6 independent components. Tracking the event horizon of the associated data then requires analysis of 10 grid functions, where each grid function consists of $`𝒪\left(Nt\times N^3\right)\times 8B=𝒪\left(160\times 121^3\right)\times 8B2.26GB`$ of data. That is, tracking the associated event horizon requires analysis of approximately $`20GB`$ of data.

A more serious difficulty associated to this data set is the relaxation time, $`t4M`$, that is typically required for outgoing data to converge onto the event horizon when followed into the past. Assuming that the collision time is (as suggested by the apparent horizon solvers) near $`t=8M`$, perturbation theory implies that the resulting horizon should undergo quasi normal ringing for another $`t20M`$ from that time. That is, at the time level $`t=10M`$, where the event horizon tracker will begin tracking into the past, it is expected that the event horizon remains highly distorted and far from its stationary regime. The problem is to determine good data for the eikonal at the time level $`t=10M`$, which can then be propagated into the past. In contrast, in analysis of the event horizon for the case of head on collision, researchers made use of approximately $`100M`$ of data , and assumed that the final state at $`t=100M`$ was a stationary or quasi stationary black hole. In such a circumstance using spherically symmetric data at $`t=100M`$ for an event horizon tracking method should prove sufficient. The difficulty with a data set only of length $`t=10M`$ is that using spherically symmetric data for the eikonal at the time level $`t=10M`$ does not accurately approximate sections of the event horizon at that time.

To accommodate data of this $`10M`$ length and to generate sufficient data for the eikonal at the time level $`t=10M`$, a method for finding a candidate section of the event horizon at $`t=10M`$ is proposed here and applied to the problem. This method makes use of the analytic properties of apparent and event horizons.

To motivate the method, note that if the event horizon were stationary at $`t=10M`$ then the expansion of the surface would be vanishing there:

$$\theta \left(t=10M\right)=\frac{1}{A}\frac{dA}{dt}=0.$$
(63)
In such a circumstance the apparent horizon could be used as initial data for the eikonal, which could then be propagated into the past. However, according to the second law of black hole mechanics, in the case that the horizon is nonstationary, which is the situation expected for this asymmetric problem, at $`t=10M`$ the horizon will satisfy

$$\theta \left(t=10M\right)0,$$
(64)
in the forward time direction. In the backwards time direction these dynamics correspond to $`\theta 0`$ at $`t=10M`$. As such, the Taylor series of any compact null surface, including the critical surface of the event horizon, is at least of the form

$$A\left(t+dt\right)=A\left(t\right)+dtA\theta +dt^2\frac{d^2A}{2dt^2}+\mathrm{}$$
(65)
and probably cannot be truncated to below

$$A\left(t+dt\right)A\left(t\right)+dtA\theta +\mathrm{}$$
(66)
The apparent horizon is then a poor estimate for the the event horizon at $`t=10M`$. The objective is to use this behavior (66) in combination with the property that outgoing null data followed into the past converges to the event horizon. The hope is to establish a better approximation for the structure of the section of the event horizon at $`t=10M`$, which can then be used in an event horizon tracking method.

To proceed, consider three compact null surfaces of outgoing data: One completely interior to the horizon $`\mathrm{\Gamma }^i`$, one completely exterior to the horizon $`\mathrm{\Gamma }^e`$, and one that is a surface of section of the event horizon $`\mathrm{\Gamma }^h`$. Let $`t_{}\mathrm{}`$ denote propagation into the past. By the property that outgoing data numerically converges to the horizon when propagated into the past, $`lim_t_{}\mathrm{}\mathrm{\Gamma }^i=lim_t_{}\mathrm{}\mathrm{\Gamma }^e=lim_t_{}\mathrm{}\mathrm{\Gamma }^h`$. Thus in most cases, if the surface $`\mathrm{\Gamma }^i\left(t_{}+dt_{}\right)`$ is pulled back to the time level $`t_{}`$ and compared to the surface $`\mathrm{\Gamma }^i\left(t_{}\right)`$ it will be completely exterior to $`\mathrm{\Gamma }^i\left(t_{}\right)`$. Similarly, if the surface $`\mathrm{\Gamma }^e\left(t_{}+dt_{}\right)`$ is pulled back to the time level $`t_{}`$ it will be completely interior to the surface $`\mathrm{\Gamma }^e\left(t_{}\right)`$. That is, by iteratively pulling the surfaces back to a single time level, spherically symmetric data will approach the numerical event horizon. Exceptions to this general behavior stem from the presence of caustics in the spacetime, where null surfaces intersect, and in neighborhoods of the event horizon. In particular, as discussed above, the horizon of this asymmetric collision data is expected to satisfy $`\theta 0`$. If the surface $`\mathrm{\Gamma }^h\left(t_{}+dt_{}\right)`$ is then pulled back to $`t_{}`$ it will be interior to the surface $`\mathrm{\Gamma }^h\left(t_{}\right)`$. Let $`\mathrm{\Gamma }_ϵ^i\left(t_{}\right)`$ be a perturbation of $`\mathrm{\Gamma }^h\left(t_{}\right)`$ that is arbitrarily close to the horizon $`\mathrm{\Gamma }^h\left(t_{}\right)`$ but completely interior to the surface $`\mathrm{\Gamma }^h\left(t_{}\right)`$. Due to round off and truncation error of the finite difference approximation and the fact that $`lim_t_{}\mathrm{}\mathrm{\Gamma }_ϵ^i=\mathrm{\Gamma }^h`$, if the numerical finite difference approximation $`\widehat{\mathrm{\Gamma }}_ϵ^i\left(t_{}+dt_{}\right)`$, which approximates the continuum section $`\mathrm{\Gamma }_ϵ^i\left(t_{}+dt_{}\right)`$, is pulled back to the time level $`t_{}`$ it will generically intersect $`\widehat{\mathrm{\Gamma }}_ϵ^i\left(t_{}\right)`$ and contain neighborhoods that are both interior and exterior to the finite difference approximation $`\widehat{\mathrm{\Gamma }}_ϵ^i\left(t_{}\right)`$. Similar behavior will hold for numerical data $`\widehat{\mathrm{\Gamma }}_ϵ^e\left(t_{}\right)`$ defined to be a perturbation of $`\widehat{\mathrm{\Gamma }}^h\left(t_{}\right)`$ that is arbitrarily close to the horizon $`\widehat{\mathrm{\Gamma }}^h\left(t_{}\right)`$ , and completely exterior to $`\widehat{\mathrm{\Gamma }}^h\left(t_{}\right)`$.

According to this behavior, data for the eikonal at the time level $`t=10M`$ can be constructed by iteration over that time slice. This approach is similar to treating the event horizon as if it were an apparent horizon, although modified to account for this nonstationary regime. Beginning with spherically symmetric initial data that is well exterior to the horizon $`\widehat{\mathrm{\Gamma }}_1^e\left(t_{}=10M\right)`$ the data is updated, pulled back to the original time level and reset as follows: $`\widehat{\mathrm{\Gamma }}_2^e\left(t_{}\right)=\widehat{\mathrm{\Gamma }}_1^e\left(t_{}+dt_{}\right)`$. The step $`\widehat{\mathrm{\Gamma }}_{n+1}^e\left(t_{}\right)=\widehat{\mathrm{\Gamma }}_n^e\left(t_{}+dt_{}\right)`$ is then repeated for several hundred iterations, which corresponds to $`𝒪\left(10\right)`$ $`e`$ folding times. Note that in the case that the spacetime is stationary the surface $`\widehat{\mathrm{\Gamma }}`$ will converge to the apparent horizon according to this method. Since the apparent horizon of a stationary spacetime coincides with the event horizon, this method will generically find the event horizon in the case of stationary spacetimes. However, in the nonstationary regime, which is the case for the asymmetric collision problem, the result of this procedure is at best an improved initial guess for an event horizon tracker. Further information, such as area analysis or study of the apparent horizon, is required to argue that the resulting surface is a candidate section of the event horizon.

We show in figure (5) the result of such a process, applying $`200`$ iterations over the time slice $`t=10M`$. Figures (5) through (11) show the eikonal function in the $`z=0`$ plane. The location of the determined guess for the surface $`\mathrm{\Gamma }`$ (which we will take to begin our evolution of the horizon into the past) is encoded into the eikonal by the color map. With this $`\mathrm{\Gamma }`$, data for the eikonal can be written in the form:

$$S(0,x^i)=1+\mathrm{tanh}\left(\frac{r_cr}{c}\right)$$
(67)
In equation (67) the first argument of the eikonal is $`t_{}=0`$; $`t_{}`$ will increase into the past. Also, $`r_c`$ denotes the data $`\mathrm{\Gamma }`$ and $`c`$ controls its steepness. Typically we take a transition width $`c`$ on the order of a computational zone. Note that this surface is not considered to be a true section of the event horizon, but instead is a good initial guess or candidate section, which after a few $`M`$ will evolve into better approximation of a true section of the event horizon. By way of comparison figure (5) also shows at $`t=10M`$ the final apparent horizon as a white wire mesh. Note that both the apparent horizon and these data for the eikonal are highly distorted from the stationary case. Figure (6) shows the resulting eikonal function $`S(x,y,z)`$ after $`2M`$ of evolution into the past. Note that the surface $`\mathrm{\Gamma }`$ is not qualitatively changed during the evolution.

## VII Surface Extraction and Apparent Horizons

Figure (7) shows several frames of the evolution of the eikonal using a viscosity solution of $`ϵ=h^2`$ (not of the improved viscosity form). This figure display the value of the eikonal function on the $`z=0`$ plane. Figure (7) shows the eikonal data as an elevation map and also via the color map of (5). Note in those figures that null surfaces interior to the event horizon undergo a change in topology and this topological transition is continuously monitored by the viscosity solutions of the eikonal. In figure (7) $`t=0.562M`$ is shown in the upper left-hand corner, $`t=1.5M`$ is shown in the upper right-hand corner, $`t=2.5M`$ is shown in the lower left-hand corner and $`t=5.0M`$ appears in the lower right-hand corner.

The figures (8) - (11) continue the sequence of figures (5) - (6), and show, for several values of $`t_{}`$, the value of the eikonal in the $`z=0`$ plane; the location of the apparent horizon in the $`3`$ dimensions (the white wire frame); and in the black wire frame, locations of sections of a candidate event horizon $`\widehat{\mathrm{\Gamma }}_c\left(t_{}\right)`$ that is generated by evolution of the eikonal equation from data constructed using the method described in the previous section. In this context, the surfaces $`\widehat{\mathrm{\Gamma }}_c\left(t_{}\right)`$ are extracted from the eikonal data using the technique described in section IV. Note that $`\widehat{\mathrm{\Gamma }}_c`$ completely contains the apparent horizons throughout their evolution. This is a fundamental condition that any numerically constructed black hole event horizon must satisfy. To determine how these results depend on the initial data $`\widehat{\mathrm{\Gamma }}_c\left(t_{}=0\right)`$, choosing initial data $`\widehat{\mathrm{\Gamma }}_\delta \left(t_{}=0\right)`$ of the form

$$u_\delta (\theta ,\varphi )=u_c(\theta ,\varphi )\delta ,$$
(68)
permits survey about the data $`u_c(\theta ,\varphi )`$, where $`u_c(\theta ,\varphi )`$ corresponds to the data $`\widehat{\mathrm{\Gamma }}_c\left(t_{}=0\right)`$. Studies with $`\delta =M/2,M,2M`$ establish that the level sets $`\widehat{\mathrm{\Gamma }}_\delta `$ penetrate both apparent horizons for any $`\delta M/2`$. These results suggest that the true event horizon is contained in a domain parameterized by $`0<\delta <M/2`$.

Figures (12) and (13) show two views of the extracted level set $`S=0`$ at $`t=2M`$. Note that this surface is highly distorted and shows the event horizon just after merger. Also shown in these figures are the apparent horizons for the two black holes. In these figures each color of the color map denotes a level set $`\widehat{\mathrm{\Gamma }}`$ of the eikonal. As such, each color represents a null surface. From these figures it is apparent that at this time there appears one innermost null surface that completely contains both apparent horizons. Thus, the results of the viscosity solutions suggest a merger time much closer to $`2M`$ then the $`8M`$ found with the apparent horizon trackers. \[Note that in figures (12) and (13) the event horizon is shown as a white wire mesh, while the apparent horizons are shown as black wire meshes. This is opposite to the color scheme used in figures 5, 6, 8 \- 11, which was an independent study of the evolution as opposed to the study of the throat geometry considered here.\]

## VIII Change of Topology

As shown in figure (7) the viscosity solutions of the eikonal equation do continuously monitor a change in topology. In the case of asymmetric binary black hole coalescence it is conjectured that the level set $`\mathrm{\Gamma }`$, which corresponds to a section of the event horizon, must take a higher genus topology at merger. To investigate this possibility, figure (14) shows the level set $`\widehat{\mathrm{\Gamma }}`$ viewed along the axis joining the centers of the apparent horizons. In that figure it is apparent that the throat function of the topological transition is elliptical in geometry. Studies indicate that this elliptical throat function persists for all null surfaces (i.e. those slightly inside or slightly outside our candidate event horizon) undergoing the topological transition. Further, for all of our computed transitions of the null surfaces, no higher genus topology is exhibited; instead, the elliptical geometry of the throat persists to the transition. These results suggest that if there is a non trivial topology in the sections of the event horizon as a consequence of the asymmetry of the merger, then that topology change is bounded to occur when the minor axis of the ellipse is within one of our computational zones, or $`h=M/4`$.

## IX Area Analysis

## X Black Hole Areas

Figure (15) shows an area versus time plot for this asymmetric collision. The curve with the lowest area is the result of a viscosity solution with $`ϵ=2h^2`$. The curve with the second lowest area is the result of a viscosity solution with $`ϵ=h^2`$, while the topmost curve is an improved viscosity solution composed of the two higher viscosity solutions. According to these results it is immediately apparent that the viscosity solutions of higher viscosity show a merger time that is later in $`t_{}`$ (and therefore prior in $`t_{}`$) then the merger time found with solutions constructed in the limit of vanishing viscosity. These results then indicate that the error (or bias) in the merger time of the viscosity solutions is directly related to the magnitude of the viscosity. More precisely, for the continuum merger time of $`t_{}^{}`$ and an approximate merger time of $`t_{}^{}{}_{ϵ}{}^{}`$, constructed using a viscosity solution of viscosity parameter $`ϵ`$, the function $`f\left(ϵ\right)=t_{}^{}{}_{ϵ}{}^{}t_{}^{}`$ is increasing in $`ϵ`$. (Here all surface areas are calculated with $`m^2`$ points where $`m=100`$.)

Figure (16) shows several area versus time curves for initial data of the form

$$u_\delta (\theta ,\varphi )=u_c(\theta ,\varphi )\delta $$
(69)
where the data $`u_c(\theta ,\varphi )`$ is that obtained using the method of section VII. From top to bottom, the curves show areas for $`\delta =0,M/4,M/2`$ and with a viscosity parameter of $`ϵ=h^2`$. Recall that studies of the apparent horizons found that the true event horizon is contained in the domain $`0<\delta <M/2`$. This survey over $`\delta `$ is conducted in search for the convergence signature associated to event horizons. Due to the time scale of this data, the time scale of the dynamics, and the relaxation time scale of the event horizon tracking method, the signature is not clearly apparent. However, this study of the area curves does show convergence of the areas, which is expected for null data approaching the horizon when propagated into the past. The curve with $`\delta =M/2`$ shows behavior suggestive of an event horizon since the $`\delta M/2`$ curves all approach that of $`\delta =M/2`$. The $`\delta =M/2`$ curve is therefore considered the best candidate for the numerical event horizon of this study. Note that the sections of this event horizon completely contain the correct apparent horizons for all $`t_{}<9M`$. Further, these $`\delta =M/2`$ data show a bifurcation time at $`t_{}8.3M`$, which corresponds to a merger time of about $`t_{}1.7M`$. This bifurcation time is detected by an algorithm that searches for any points $`\theta ,\varphi `$ of the surface such that $`u(\theta ,\varphi )<h`$. In the circumstance that a point is found such that $`u(\theta ,\varphi )<h`$, the bifurcation is expected to occur in a few more $`dt=M/16`$ in the $`t_{}`$ direction. A few time levels prior to this bifurcation time in $`t_{}`$ (i.e. just after the bifurcation in $`t`$), the single merged horizon component has an area (the area of the level set) computed to be $`A=148.9`$. Based on analysis of area computations for exact solutions (figure(3), we anticipate an error in the horizon area of several percent. We conservatively assign an $`8\%`$ error to the areas. Just prior (in $`t`$) to merger the area of each black hole is then $`A=74.5\pm 6.0M^2`$. These individual areas correspond to a Schwarzschild mass of about $`M=1.48\pm 0.12`$. This result is a substantially larger mass for each hole then determined by the apparent horizon finders at the time level $`t_{}=0`$. Interestingly, studies have found that apparent horizons separated by about $`10M`$ show an increase in their mass due to the effects of binding energy. Individual masses of about $`M1.36`$ are then only a slight departure from studies that account for the binding energy of the holes. Further, at the time of merger the holes have undergone $`1.7M`$ of evolution, during which the holes could accrete any surrounding gravitational radiation present in the initial data. The presence of such radiation would lead to larger masses then those found using apparent horizon finders at $`t_{}=0`$. However, it is important to note that due to the viscosity in the solution the result $`M1.36`$ can only properly be considered as a lower bound on the calculated masses. The most significant contribution to any error in this result must stem from the relatively small time scale of this asymmetric collision data and coupling of that time scale to the $`e`$folding time scale of this event horizon detection method.

## XI Conclusions

In this work we have demonstrated a relatively simple yet robust and (most importantly) generic solution to the problem of numerically tracking black hole event horizons. An implementation of our method made use of an explicit second order diffusion term to regulate the solution singularities associated to caustics. As demonstrated by analysis of analytic sources, this term does introduce numerical error although we demonstrate our control over these effects and the resulting accuracy. But, the use of a second order diffusion term is not required by our method per se; and a variety of other approaches can be employed. Examples of other methods for controlling breakdown of the numerical solution include those classes of high resolution shock capturing numerical schemes that are used extensively in computational fluid dynamics for hyperbolic problems similar to the eikonal equation.

The application of our new method for event horizon tracking method considered the asymmetric binary black hole coalescence problem, including a detailed analysis of areas of the surfaces of sections, the collision time, associated apparent horizons, and the topology of the horizon. Due to the relatively short time scale of the collision data, our method was unable to demonstrate the signature of the black hole event horizon. We believe that this problem is due to the data itself and not due to our method. We anticipate much more accurate and convincing results as more accurate computational simulations of black hole interactions become available.