format: text

text:
# Photoemission spectroscopy and sum rules in dilute electron-phonon systems

\[

## Abstract

A family of exact sum rules for the one-polaron spectral function in the low-density limit is derived. An algorithm to calculate energy moments of arbitrary order of the spectral function is presented. Explicit expressions are given for the first two moments of a model with general electron-phonon interaction, and for the first four moments of the Holstein polaron. The sum rules are linked to experiments on momentum-resolved photoemission spectroscopy. The bare electronic dispersion and the electron-phonon coupling constant can be extracted from the first and second moments of spectrum. The sum rules could serve as constraints in analytical and numerical studies of electron-phonon models.

\]

Statistical analysis of spectra is an important branch of photoemission spectroscopy. It has been known that the energy moments of the spectral density of atoms can be expressed via the atomic quantum numbers . Such integral relations, or sum rules, provide direct link between experimental data and the internal state of the atoms. Similarly, in crystals analogous sum rules yield valuable information about the Bloch states. This is particularly true in case of two-dimensional systems where the in-plane momentum of the electron does not change in the photoemission process. It allows the spectrum to be recorded for each in-plane momentum individually. In recent years, the angle-resolved photoemission spectroscopy became a powerful experimental technique to study the electronic structure of quasi-two-dimensional compounds such as the high-temperature superconductors and layered colossal magnetoresistive oxides .

This paper discusses statistical analysis of spectra from systems with strong electron-phonon (el-ph) interaction. In such systems, the mobile carriers deform the lattice and form polarons which are bound states of the carriers and the lattice deformations. The polarons then define the low-energy physics of the system, in particular the equilibrium chemical potential. In inverse photoemission (IPES), an electron with momentum k is injected into the system from outside. Initially it feels the bare electronic states in the undeformed lattice and only then begins to deform the lattice around itself. Thus it should be expected that the energy moments of the IPES intensity could provide information about the initial stages of the deformation process and therefore about the strength of the interaction. It turns out that for low density of polarons the moments of IPES intensity can be calculated analytically for a wide class of model el-ph interactions. It leads to a series of exact sum rules which the spectral density, experimental or theoretical, must satisfy. Below I derive the sum rules and discuss their possible applications.

Consider a single band, one optical branch Hamiltonian with an el-ph interaction of arbitrary profile

$`={\displaystyle \underset{𝐤}{}}(\epsilon _𝐤\mu )c_𝐤^{}c_𝐤`$ $``$ $`\sqrt{2}\mathrm{\Omega }{\displaystyle \underset{\mathrm{𝐧𝐦}}{}}g(𝐧𝐦)c_𝐧^{}c_𝐧x_𝐦`$ (1)
$`+`$ $`{\displaystyle \frac{\mathrm{\Omega }}{2}}{\displaystyle \underset{𝐦}{}}(_𝐦^2+x_𝐦^2).`$ (2)
Here $`\mu `$ is the chemical potential, $`x_𝐦=\xi _𝐦\sqrt{\frac{M\mathrm{\Omega }}{\mathrm{}}}`$ is the dimensionless displacement of an oscillator at m-th lattice site, $`M`$ and $`\mathrm{\Omega }`$ are the mass and the frequency of the oscillators, respectively, $`_𝐦/x_𝐦`$, $`c_𝐧=N^{1/2}_𝐤e^{i\mathrm{𝐤𝐧}}c_𝐤`$, and $`N`$ is the number of lattice sites. $`g(𝐧𝐦)`$ is the dimensionless coupling constant proportional to the force with which an electron at site n acts on oscillator m. For simplicity, the fermion spin and multiple phonon polarizations are not included explicitly in Eq.(2). They are not essential for the purposes of the paper and could be added if necessary. The canonical Holstein model  follows from the above Hamiltonian with the choice $`g(𝐧𝐦)=g_0\delta _{\mathrm{𝐧𝐦}}`$.

Within the sudden approximation, direct and inverse photoemission intensities are proportional to the fermion spectral functions $`A_{}(𝐤,\omega )`$ and $`A_+(𝐤,\omega )`$, respectively (hereafter notation $`\mathrm{}=1`$ is used):

$$A_{}(𝐤,\omega )=\frac{2\pi }{Z}\underset{ij}{}e^{\beta _i}|j|c_𝐤|i|^2\delta (\omega +_j_i),$$
(3)
$$A_+(𝐤,\omega )=\frac{2\pi }{Z}\underset{ij}{}e^{\beta _i}|j|c_𝐤^{}|i|^2\delta (\omega +_i_j).$$
(4)
Here $`\omega `$ is the energy relative to the chemical potential, $`|i`$ and $`|j`$ are exact eigenstates of $``$, with energies $`_i`$ and $`_j`$, $`Z=_ne^{\beta _i}`$ is the grand canonical partition function, and $`\beta =(k_BT)^1`$. The validity of the sudden approximation as well as the additional electron-photon matrix element, the Coulomb effects and other complications, are not going to be discussed in this paper, see e.g. . It follows from Eqs. (3) and (4) that the two spectral functions obey the fundamental relation

$$A_+(𝐤,\omega )=e^{\beta \omega }A_{}(𝐤,\omega ).$$
(5)
Thus, $`A_+`$ can be found if $`A_{}`$ is known with sufficient accuracy, and visa versa. In terms of the symmetrized function $`A=A_++A_{}`$ they are expressed as $`A_{}(𝐤,\omega )=(e^{\beta \omega }+1)^1A(𝐤,\omega )`$ and $`A_+(𝐤,\omega )=(e^{\beta \omega }+1)^1A(𝐤,\omega )`$.

Several theoretical results for the polaron spectral function are available. The on-site version of the Holstein model ($`\epsilon _𝐤0`$) is exactly solvable, see e.g. Ref. . For the Holstein model with non-zero kinetic energy the spectral function was calculated by Alexandrov and Ranninger in the Lang-Firsov approximation . Later their result was generalized to the long-range interaction (2) by Alexandrov and Sricheewin . For the more complex Jahn-Teller interaction the polaron spectral function was given by Perebeinos and Allen . The spectral functions mentioned are approximate and all have the form of a sequence of delta-functions.

The formalism of statistical moments is best introduced via the time-dependent correlation functions $`K_\pm (𝐤,t)`$ that are Fourier transforms of $`A_\pm (𝐤,\omega )`$:

$`K_+(𝐤,t)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}e^{i\omega t}A_+(𝐤,\omega )`$ (6)
$`=`$ $`{\displaystyle \frac{1}{Z}}{\displaystyle \underset{j}{}}e^{\beta _j}j|e^{it}c_𝐤e^{it}c_𝐤^{}|j`$ (7)
$``$ $`e^{it}c_𝐤e^{it}c_𝐤^{},`$ (8)
$`K_{}(𝐤,t)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}e^{i\omega t}A_{}(𝐤,\omega )`$ (9)
$`=`$ $`{\displaystyle \frac{1}{Z}}{\displaystyle \underset{j}{}}e^{\beta _j}j|c_𝐤^{}e^{it}c_𝐤e^{it}|j`$ (10)
$`=`$ $`c_𝐤^{}e^{it}c_𝐤e^{it}.`$ (11)
Taking the $`m`$-th derivative with respect to time $`t`$ and setting $`t=0`$ one obtains

$`_m^+(𝐤)`$ $``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}\omega ^mA_+(𝐤,\omega )`$ (12)
$`=`$ $`\underset{m\mathrm{times}}{\underset{}{[[[c_𝐤,],]\mathrm{},]}}c_𝐤^{},`$ (13)
$`_m^{}(𝐤)`$ $``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}\omega ^mA_{}(𝐤,\omega )`$ (14)
$`=`$ $`c_𝐤^{}\underset{m\mathrm{times}}{\underset{}{[[[c_𝐤,],]\mathrm{},]}}.`$ (15)
The moments of the symmetrized function $`A(𝐤,\omega )`$ are defined as $`_m(𝐤)=_m^+(𝐤)+_m^{}(𝐤)`$.

The simplest sum rule is the normalization property of $`A(𝐤,\omega )`$ which is common to any one-particle fermionic spectral function

$$_0(𝐤)=c_𝐤c_𝐤^{}+c_𝐤^{}c_𝐤=1.$$
(16)
The first moment of $`A(𝐤,\omega )`$ is given by

$`[c_𝐤,]`$ $`=`$ $`(\epsilon _𝐤\mu )c_𝐤`$ (17)
$``$ $`{\displaystyle \frac{\sqrt{2}\mathrm{\Omega }}{N}}{\displaystyle \underset{\mathrm{𝐧𝐦𝐤}^{}}{}}e^{i(𝐤^{}𝐤)𝐧}g(𝐧𝐦)x_𝐧c_𝐤^{}.`$ (18)
$`_1(𝐤)`$ $`=`$ $`[c_𝐤,]c_𝐤^{}+c_𝐤^{}[c_𝐤,]`$ (19)
$`=`$ $`(\epsilon _𝐤\mu ){\displaystyle \frac{\sqrt{2}\mathrm{\Omega }}{N}}{\displaystyle \underset{\mathrm{𝐧𝐦}}{}}g(𝐧𝐦)x_𝐦.`$ (20)
The mean displacement of the oscillators $`x_𝐦`$ depends only on the total number of carriers in the systems $`n`$, and can be calculated as follows. Write the displacement operator as $`x_𝐦=x_𝐦+y_𝐦`$ and substitute in $``$. Then minimization of the free energy $`F=\beta ^1\mathrm{ln}e^\beta `$ with respect to $`x_𝐦`$ yields

$`x_𝐦`$ $`=`$ $`y_𝐦+\sqrt{2}{\displaystyle \underset{𝐧}{}}g(𝐧𝐦)c_𝐧^{}c_𝐧`$ (21)
$`=`$ $`\sqrt{2}n{\displaystyle \underset{𝐧}{}}g(𝐧),`$ (22)
because $`y_𝐦=0`$ by definition, and the mean density $`c_𝐧^{}c_𝐧=n`$ does not depend on coordinate $`𝐧`$. Substitution in Eq. (20) results in:

$$_1(𝐤)=(\epsilon _𝐤\mu )2\mathrm{\Omega }n\left[\underset{𝐧}{}g(𝐧)\right]^2.$$
(23)
This is an exact relation for the one-particle spectral function. It is a generalization of the sum rule derived by Mahan for the on-site polaron , and by Perebeinos and Allen for the Jahn-Teller polaron , to non-zero polaron density. The main feature of these results is that the mean value of the spectral function is not affected by interaction with phonons. This conclusion is non-trivial in view of significant changes of the spectral function itself, especially in the strong coupling regime. As follows from Eq. (23), in the many-body case the first moment $`_1(𝐤)`$ acquires a correction that is linear in polaron density. Proceeding with the calculation one obtains for the second moment

$`[[c_𝐤,],]`$ $`=`$ $`(\epsilon _𝐤^2\mu )c_𝐤`$ (24)
$``$ $`{\displaystyle \frac{\sqrt{2}\mathrm{\Omega }}{N}}{\displaystyle \underset{𝐤^{}\mathrm{𝐧𝐦}}{}}e^{i(𝐤^{}𝐤)𝐧}g(𝐧𝐦)`$ (25)
$`\times `$ $`[(\epsilon _𝐤+\epsilon _𝐤^{}2\mu )x_𝐦+\mathrm{\Omega }_𝐦]c_𝐤^{}`$ (26)
$`+`$ $`{\displaystyle \frac{2\mathrm{\Omega }^2}{N}}{\displaystyle \underset{𝐤^{}\mathrm{𝐧𝐦𝐦}^{}}{}}e^{i(𝐤^{}𝐤)𝐧}g(𝐧𝐦)`$ (27)
$`\times `$ $`g(𝐧𝐦^{})x_𝐦x_𝐦^{}c_𝐤^{},`$ (28)
$`_2(𝐤)`$ $`=`$ $`[[c_𝐤,],]c_𝐤^{}+c_𝐤^{}[[c_𝐤,],]=(\epsilon _𝐤\mu )^2`$ (29)
$``$ $`{\displaystyle \frac{\sqrt{2}\mathrm{\Omega }}{N}}{\displaystyle \underset{\mathrm{𝐧𝐦}}{}}g(𝐧𝐦)[2(\epsilon _𝐤\mu )x_𝐦+\mathrm{\Omega }_𝐦`$ (30)
$`+`$ $`{\displaystyle \frac{2\mathrm{\Omega }^2}{N}}{\displaystyle \underset{\mathrm{𝐧𝐦𝐦}^{}}{}}g(𝐧𝐦)g(𝐧𝐦^{})x_𝐦x_𝐦^{}.`$ (31)
Unfortunately, the mean values $`x_𝐦x_𝐦^{}`$ cannot be calculated analytically which limits application of the last expression. In principle, such averages could be computed numerically by various techniques, rendering Eq. (31) a true sum rule. (It is worth noting in this respect that in the case of the Hubbard model, the first and second moments of the spectral function were calculated analytically by Nolting . This finding enabled him to develop the moment approach to the Hubbard model.)

Further progress can be made by going to the low density limit $`n1`$. It follows from definitions (3) and (4) that $`A_{}n`$ while $`A_+1n`$. Therefore $`_m^{}_m^+`$ so that $`_m^{}`$ could be dropped from Eqs. (16), (20), and (31). More importantly, all the phonon operators can now be averaged easily. Indeed, in the low-density limit the lattice is largely undeformed. The function $`A_+`$ describes the process when an electron is created in an empty lattice and removed later. Thus whatever phonon operators appear as a result of multiple commutation with the Hamiltonian, they have to be averaged over the system of undeformed independent harmonic oscillators, which is always possible. For instance, the mean values entering Eq. (31) are $`x_𝐦=_𝐦=0`$, $`x_𝐦x_𝐦^{}=\frac{1}{2}\mathrm{coth}(\frac{1}{2}\beta \mathrm{\Omega })\delta _{\mathrm{𝐦𝐦}^{}}`$. Collecting all the results one obtains the following sum rules for the spectral function $`A_+(𝐤,\omega )`$:

$`_0^+(𝐤)`$ $`=`$ $`1+𝒪(n),`$ (32)
$`_1^+(𝐤)`$ $`=`$ $`(\epsilon _𝐤\mu )+𝒪(n),`$ (33)
$`_2^+(𝐤)`$ $`=`$ $`(\epsilon _𝐤\mu )^2`$ (34)
$`+`$ $`\mathrm{\Omega }^2\mathrm{coth}({\displaystyle \frac{1}{2}}\beta \mathrm{\Omega }){\displaystyle \underset{𝐧}{}}g^2(𝐧)+𝒪(n).`$ (35)
The above exact relations may be used for statistical analysis of inverse photoemission spectra (IPES) from low-density polaron systems. First of all, Eq. (32) allows proper normalization of the spectrum for each electron momentum k. According to Eq. (33), the first moment yields the bare electronic spectrum inside the crystal. However, in an el-ph system the chemical potential $`\mu `$ is lowered by the polaron binding energy $`E_p`$. As a function of $`𝐤`$, $`_1^+(𝐤)`$ never gets smaller than $`E_p`$. In other words, one should expect a gap of size $`E_p`$ in IPES. The second moment $`_2^+(𝐤)`$ also contains information about the strength of el-ph interaction. It is convenient to parameterize the latter in terms of the zero order (in hopping) binding energy of the polaron, $`E_p^{(0)}\mathrm{\Omega }_𝐧g^2(𝐧)`$, see . ($`E_p^{(0)}`$ is employed here as a shorthand notation only. It does not include the effects of the kinetic energy and therefore it is not the true polaron energy, $`E_p^{(0)}E_p`$.) The mean squared deviation of the spectrum, $`\overline{(\mathrm{\Delta }\omega )^2}=_2^+(_1^+)^2=E_p^{(0)}\mathrm{\Omega }\mathrm{coth}(\frac{1}{2}\beta \mathrm{\Omega })`$, is a direct measure of $`E_p^{(0)}`$. It has to scale approximately linearly with the IPES gap. Notice, that $`\overline{(\mathrm{\Delta }\omega )^2}=E_p^{(0)}\mathrm{\Omega }`$ at low temperatures and $`\overline{(\mathrm{\Delta }\omega )^2}2E_p^{(0)}T`$ at high temperatures $`T`$.

One should mention that sum rules on the function $`A_{}(𝐤,\omega )`$ have also been discussed in relation to direct photoemission. In particular, Randeria et al introduced an approximate sum rule for the zeroth moment, $`_0^{}=n(𝐤)`$, where $`n(𝐤)`$ is the mean occupation of the state $`𝐤`$ . This relation may not be generalizable to polaronic systems because in general $`n(𝐤)`$ depends on the strength of the el-ph interaction and therefore is not known a priori.

Calculation of moments can be continued. The low density ensures that any fermionic operator that will arise after the multiple commutation of $`c_𝐤`$ with the Hamiltonian, will factor out. The remaining phonon operators may be complex but they will have to be averaged over the system of free phonons. Such an averaging can always be performed. The conclusion is that any moment of $`A_+(𝐤,\omega )`$ can in principle be calculated analytically. However, the algebra quickly becomes cumbersome. Therefore, I present only two more moments for the Holstein model with nearest neighbor hopping, i.e. for $`g(𝐧𝐦)=g_0\delta _{\mathrm{𝐧𝐦}}`$ and $`\epsilon _𝐤=J_𝐥e^{i\mathrm{𝐤𝐥}}`$, l numbering the nearest neighbors:

$$_3^+(𝐤)=(\epsilon _𝐤\mu )^3+2g_T^2\mathrm{\Omega }^2(\epsilon _𝐤\mu )+g_0^2\mathrm{\Omega }^3+𝒪(n),$$
(36)
$`_4^+(𝐤)`$ $`=`$ $`(\epsilon _𝐤\mu )^4+g_T^2\mathrm{\Omega }^2[3(\epsilon _𝐤\mu )^2+zJ^2]`$ (37)
$`+`$ $`2g_0^2\mathrm{\Omega }^3(\epsilon _𝐤\mu )+[g_T^2+3g_T^4]\mathrm{\Omega }^4+𝒪(n),`$ (38)
where $`g_T^2g_0^2\mathrm{coth}(\frac{1}{2}\beta \mathrm{\Omega })`$, and $`z`$ is the number of nearest neighbors in the lattice. In deriving $`_3^+(𝐤)`$ and $`_4^+(𝐤)`$ the following properties of the harmonic oscillator have been used: $`x_𝐦_𝐦=\frac{1}{2}`$, and $`x_𝐦^4=\frac{3}{4}\mathrm{coth}^2(\frac{1}{2}\beta \mathrm{\Omega })`$ . It is interesting that in the fourth moment the hopping integral $`J`$ appears not only via the bare spectrum $`\epsilon _𝐤`$ but also by itself. Thus $`_4^+(𝐤)`$ is the first moment that explicitly distinguishes the dimensionality of the lattice through the number of nearest neighbors $`z`$. Apparently the higher order moments also distinguish different lattice topologies.

Apart from their relation to photoemission spectroscopy the sum rules may find more “theoretical” applications. First of all, they are valid for the single polaron which continues to draw theoretical interest. Formal transition in the above expressions to the canonical ensemble with one particle is done by setting $`\mu =0`$. Thus the sum rules could be useful checks for polaron spectral functions computed numerically by exact diagonalization , density matrix renormalization group , quantum Monte Carlo , or any other method. Secondly, the sum rules may provide information on the dynamics of polaron formation. For instance, the quantity $`P_+(𝐤,t)=|K_+(𝐤,t)|^2`$ is the probability that the phonon system remains unchanged after time $`t`$ under the disturbance of an electron with momentum $`𝐤`$. On the other hand, $`P_+(𝐤,t)=1[_2^+(𝐤)_1^{+2}(𝐤)]t^2+o(t^2)=1E_p^{(0)}\mathrm{\Omega }\mathrm{coth}(\frac{1}{2}\beta \mathrm{\Omega })t^2+o(t^2)`$. It enables one to interpret $`\tau =[E_p^{(0)}\mathrm{\Omega }\mathrm{coth}(\frac{1}{2}\beta \mathrm{\Omega })]^{1/2}`$ as the average time of phonon emission. This is a non-perturbative result that is valid for any model parameters. Finally, the sum rules may be employed to improve approximate analytic expressions for polaron spectral functions. Consider the two functions ($`\mu =0`$, $`T=0`$):

$`A_+^{\mathrm{low}}(𝐤,\omega )=(2\pi )e^{g_0^2}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{g_0^{2m}}{m!}}\times `$ (39)
$`\delta (\omega +g_0^2\mathrm{\Omega }e^{g_0^2}\epsilon _𝐤m\mathrm{\Omega }),`$ (40)
$`A_+^{\mathrm{high}}(𝐤,\omega )=(2\pi )e^{g_0^2}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{g_0^{2m}}{m!}}\times `$ (41)
$`\delta (\omega +g_0^2\mathrm{\Omega }\epsilon _𝐤m\mathrm{\Omega }).`$ (42)
The first function is the one-polaron limit of Alexandrov and Ranninger’s result for the Holstein model in the Lang-Firsov approximation . In Eq.(40), $`g_0^2\mathrm{\Omega }`$ is the zero-order binding energy $`E_p^{(0)}`$. $`A_+^{\mathrm{low}}`$ correctly describes the low-energy physics of the polaron, most notably the polaron band narrowing. However, it satisfies neither of the sum rules except the zeroth one. In particular, the first moment is equal to $`e^{g_0^2}\epsilon _𝐤`$ instead of the correct $`\epsilon _𝐤`$. The second function $`A_+^{\mathrm{high}}`$ is a “wild guess”. It does satisfy the first and second moments (but not the third and fourth ones) but does not describe the band narrowing. $`A_+^{\mathrm{high}}`$ is wrong at low energies but may be accurate at high energies. This conjecture is supported by the anticipation that at high energies the electron should become “free” from phonons. Thus a good spectral function could be the one that interpolates from $`A_+^{\mathrm{low}}`$ at low energies to $`A_+^{\mathrm{high}}`$ at high energies. With reasonable parameterization, the interpolation parameters could also be deduced from the sum rules.

In conclusion, a series of sum rules on the polaron spectral function have been derived. The sum rules are valid in the low-density limit. The first moment has been generalized to finite density. A connection with the angle-resolved inverse photoemission spectroscopy has been established. The spectrum should display a gap which size correlates with the mean square deviation of the spectral function. Several other applications have been identified. The sum rules are open to further generalizations to dispersive phonons and more complex interactions.

I am grateful to A. M. Bratkovsky, D. Stewart, J. E. Hirsch, and especially to A. S. Alexandrov for very useful discussions on the subject and critical reading of the manuscript.


text:
# On the Properties of Time-Dependent, Force-Free, Degenerate Electrodynamics

## 1 Introduction

In magnetospheres of pulsars and black holes the electromagnetic field is so strong that the inertia and pressure of the plasma can be ignored. As a result the Lorentz force, $`F_{\alpha \beta }J^\alpha `$, vanishes and the transport of energy and momentum is entirely electromagnetic \[Goldreich & Julian 1969, Blandford & Znajek 1977\]. This justifies the name “force-free” to describe the electrodynamics of pulsars and black holes. However, the electrodynamics of the magnetospheres is rather different from that in a vacuum which is obviously also force-free. Indeed, the magnetospheric plasma is plentiful enough to support strong electric currents and screen the electric field \[Goldreich & Julian 1969, Blandford & Znajek 1977\]. An electromagnetic field satisfying this condition is called “degenerate” \[Macdonald & Thorne 1982\]. Thus, the electrodynamics of pulsar and black hole magnetospheres is force-free, degenerate electrodynamics (FFDE).

So far, theoretical studies of pulsar and black hole magnetospheres have assumed a steady-state, indeed the properties of FFDE as a system of time-dependent equation have not been studied systematically. It is only recently that Uchida \[Uchida1997\] has developed a rather elegant theory of this system in which the electromagnetic field is described in terms of two scalar functions, called “Euler potentials”. However, this formulation is not particularly suitable for numerical analysis because its basic equations, when written in components, involve mixed space and time second order derivatives. In this paper, we present another formulation in which the evolution equations have the form of hyperbolic conservation laws. We also describe a simple one-dimensional, upwind numerical scheme based on these analytic results. This scheme can easily be generalized to multidimensions and curved space-time. In fact, we have already constructed a 2D scheme for the Kerr space-time and used it \[Komissarov 2001\] to study the Blandford-Znajek mechanism \[Blandford & Znajek 1977\] for the extraction of rotational energy from black holes.

## 2 Evolution Equations of Degenerate Electrodynamics

Since far away from the relativistic “star” the inertia of particles is bound to become important, the approximation of ideal relativistic magnetohydrodynamics (RMHD) could be considered as more satisfactory in general \[Phinney1982\]. This suggests one way of deriving the time-dependent equations of FFDE, namely by considering the low inertia limit of RMDH. We shell adopt this approach because it reveals the close connection between these two systems. However, the time-dependent equations of FFDE can also be derived entirely within the frame of electrodynamics in a rather obvious way. Besides that, some of the microphysical conditions of MHD, like the small mean free path of plasma particles, do not have to be satisfied in FFDE.

In the covariant form the equations of RMHD are as follows \[Lichnerowicz 1967, Anile 1989\]

Maxwell’s equation:

$$_\alpha ^{}𝑭^{\alpha \beta }=0,$$
(1)
Energy-momentum equation:

$$_\alpha 𝑻^{\alpha \beta }=0,$$
(2)
Continuity equations

$$_\alpha (n_{(i)}𝒖^\alpha )=0,$$
(3)
where $`𝑻^{\alpha \beta }`$ is the total stress-energy-momentum tensor, $`{}_{}{}^{}𝑭_{}^{\alpha \beta }`$ is the dual tensor of electromagnetic field, $`n_{(i)}`$ is the proper number density for species $`i`$ and $`𝒖^\alpha `$ is the fluid 4-velocity. The total stress-energy-momentum tensor

$$𝑻^{\alpha \beta }=𝑻_{(m)}^{\alpha \beta }+𝑻_{(f)}^{\alpha \beta },$$
(4)
is the sum of the stress-energy-momentum tensor of matter, $`𝑻_{(m)}^{\alpha \beta }`$, and the stress-energy-momentum tensor of electromagnetic field

$$𝑻_{(f)}^{\alpha \beta }=𝑭_\mu ^\alpha 𝑭^{\beta \mu }\frac{1}{4}(𝑭_{\mu \nu }𝑭^{\mu \nu })𝒈^{\alpha \beta },$$
(5)
where $`𝒈_{\alpha \beta }`$ is the metric tensor and $`𝑭_{\alpha \beta }`$ is the electromagnetic field tensor. In this paper we assume that Greek indices run from 0 to 3 and Latin indices from 1 to 3.

These differential equations (23) are supplemented by equations of state and the perfect conductivity condition

$$𝑭_{\nu \mu }𝒖^\mu =0.$$
(6)
This condition makes the second group of Maxwell equations redundant. From (6) one obtains

$$𝑭^{}𝑭=0,$$
(7)
$$𝑭𝑭>0,$$
(8)
where $`𝑭𝑭`$ means $`𝑭_{\mu \nu }𝑭^{\mu \nu }`$. Just like classical MHD, RMHD is a hyperbolic system of conservation laws and has waves of four types – entropy, fast, Alfvén, and slow (e.g Anile 1989; Komissarov 1999).

If the contribution of matter to the total tensor of stress-energy momentum is very small then a perturbation technique suggests itself. The zero order equations constitute the system of FFDE:

$$_\mu ^{}𝑭^{\mu \nu }=0,$$
(9)
$$_\mu 𝑻_{(f)}^{\nu \mu }=0,$$
(10)
From (5,9) one has

$$_\mu 𝑻_{(f)}^{\nu \mu }=𝑭_{\nu \mu }𝑱^\mu ,$$
(11)
where

$$𝑱^\mu =_\nu 𝑭^{\mu \nu }$$

and, thus, (10) requires the vanishing of the Lorentz force. Since equations (9,10) do not involve $`𝒖_\alpha `$, the perfect conductivity condition should be replaced with conditions (7,8), which express the degeneracy of electromagnetic field in FFDE. They ensure the existence of observers who detect only magnetic field. Let $`𝒂^\mu `$ to be the velocity of such observer. It easy to see that

$$𝑭_{\mu \nu }𝒂^\nu =0$$

and, thus, $`𝒂^\nu `$ is a zero eigenvector of $`𝑭_{\mu \nu }`$. The second zero eigenvector can then be introduced via

$$𝒃^\mu =^{}𝑭^{\nu \mu }𝒂_\nu .$$

In terms of $`𝒂^\mu `$ and $`𝒃^\mu `$

$$𝑭_{\mu \nu }=𝒆_{\mu \nu \alpha \beta }𝒂^\alpha 𝒃^\beta ,$$
(12)
and

$${}_{}{}^{}𝑭_{}^{\mu \nu }=2𝒃^{[\mu }𝒂^{\nu ]}.$$
(13)
From the last equation it follows that the unit space-like vector $`𝒄^\mu `$ orthogonal to $`𝒂^\mu `$ and $`𝒃^\mu `$ is a zero eigenvector of $`{}_{}{}^{}𝑭_{\mu \nu }^{}`$. The second zero eigenvector, also space-like and orthogonal to $`𝒂^\mu `$ and $`𝒃^\mu `$, is given by

$$𝒅^\mu =𝑭^{\nu \mu }𝒄_\nu .$$

In terms of these vectors

$${}_{}{}^{}𝑭_{\mu \nu }^{}=𝒆_{\mu \nu \alpha \beta }𝒄^\alpha 𝒅^\beta ,$$
(14)
and

$$𝑭^{\mu \nu }=2𝒄^{[\mu }𝒅^{\nu ]}.$$
(15)
From (7) it follows that there are only five independent components of $`𝑭_{\mu \nu }`$. In components, (9) splits into three evolution equations and one differential constraint on the initial data ($`_iB^i=0`$). Thus, if the system (7)–(10) is self-consistent then (10) has only two independent components. Indeed, equations (11,12) show that that in the limit of degenerate electrodynamics, the vector $`_\mu 𝑻^{\nu \mu }`$ always belongs to the two-dimensional vector space generated by $`𝒄^\mu `$ and $`𝒅^\mu `$.

The evolution equations of the system (9,10) have the form of conservation laws. In the following section we shall show that this system is also hyperbolic.

## 3 Hyperbolicity of Degenerate Electrodynamics

In order to demonstrate the hyperbolicity of FFDE, it is sufficient to study its one-dimensional equations in a locally pseudo-Cartesian system of coordinates. In such a frame, the electric and magnetic fields are defined via

$$E_i=F_{i0},B^i=\frac{1}{2}ϵ^{ijk}F_{jk},$$
(16)
where $`ϵ^{ijk}`$ is the 3-dimensional Levi-Civita symbol. Then

$$T^{00}=\frac{1}{2}(B^2+E^2),T^{i0}=ϵ^{ijk}E_jB_k,$$

$$T^{ij}=(E^iE^j+B^iB^j)+\frac{\delta ^{ij}}{2}(E^2+B^2),$$

and conditions (7,8) read

$$\stackrel{\mathbf{}}{𝑬}\mathbf{}\stackrel{\mathbf{}}{𝑩}=0,$$
(17)
$$B^2E^2>0.$$
(18)
Notice that we use units such that the speed of light and $`\pi `$ do not appear in the equations.

Like in MHD, the number of evolution equations and dependent variables for the one-dimensional system is reduced by one. Indeed, if $`/x^2=/x^3=0`$ then $`x^0`$\- ($`t`$-) and the $`x^1`$-components of (9) require

$$B_1=const.$$

Thus, the original system of FFDE is reduced to four independent evolution equations for four dependent variables, $`B_2`$, $`B_3`$, and two components of the electric field. The selection of these two components depends on the direction of magnetic field. For example, if $`B_1=0`$ then one of them must be $`E_1`$. As a result, we have to study three different sets of 1D equations to cover all possible cases. One can avoid such complications by relaxing the degeneracy condition (17) and modifying the equations in a way such that this condition is satisfied automatically if it is satisfied by the initial data. For example, one can select all three space components of the energy-momentum equation and add the term $`B_i(B_iE^i)/t`$ to obtain

$$F_{i\mu }\frac{F^{\mu t}}{t}+F_{i\mu }\frac{F^{\mu 1}}{x^1}+B_i\frac{(B_iE^i)}{t}=0.$$
(19)
Together with the two remaining components of (9) these equations constitute what we shell call the “augmented one-dimensional system of FFDE”. It has five independent equations for five components of $`𝑭_{\nu \mu }`$. Obviously, solutions of this augmented system satisfying the constraint (17) are also solutions of the original system of FFDE. Moreover, contraction of (19) with $`B^i`$ gives

$$(B^iE_i)\frac{E^1}{x^1}+B^2\frac{(B^iE_i)}{t}=0$$

This result shows that if $`B^iE_i=0`$ at $`t=0`$ then all time derivatives of $`B^iE_i`$ at $`t=0`$ also vanish and, therefore, $`B^iE_i=0`$ at $`t>0`$ as well.

All characteristic waves of the FFDE must be present in the augmented system because of the way it has been constructed. It is also expected to have one additional nonphysical wave due to higher number of evolution equations. This extra wave can easily be identified because across this wave $`E_iB^i`$ does not have to stay constant.

In vector form the augmented system reads

$$𝑨\frac{𝑼}{t}+𝑪\frac{𝑼}{x^1}=0,$$
(20)
where

$$𝑼=(B_2,B_3,E_1,E_2,E_3)^t,$$
(21)
is the vector of dependent variables,

$$𝑨=\left(\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ B_1E_2\hfill & B_1E_3\hfill & B_1^2\hfill & B_1B_2B_3\hfill & B_1B_3+B_2\hfill \\ B_2E_2\hfill & B_2E_3\hfill & B_2B_1+B_3\hfill & B_2^2\hfill & B_2B_3B_1\hfill \\ B_3E_2\hfill & B_3E_3\hfill & B_3B_1B_2\hfill & B_2B_3+B_1\hfill & B_3^2\hfill \end{array}\right),$$
(22)
$$𝑪=\left(\begin{array}{ccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ B_2\hfill & B_3\hfill & E_1\hfill & 0\hfill & 0\hfill \\ B_1\hfill & 0\hfill & E_2\hfill & 0\hfill & 0\hfill \\ 0\hfill & B_1\hfill & E_3\hfill & 0\hfill & 0\hfill \end{array}\right).$$
(23)
Given (17,18) the eigenvalue problem

$$(𝑨\mu 𝑪)𝒓=0$$
(24)
has the following solutions

$$\mu _f^\pm =\pm 1,$$
(25)
$$\mu _a^\pm =\frac{B_3E_2B_2E_3\pm \sqrt{B_1^2(B^2E^2)}}{B^2},$$
(26)
$$\mu _n=0.$$
(27)
It is easy to verify that in the limit of FFDE the wave speeds of the fast waves of relativistic MHD \[Anile 1989, Komissarov 1999\] reduce to $`\mu _f^\pm `$ and the wave speeds of Alfvén waves of relativistic MHD reduce to $`\mu _a^\pm `$. Hence, $`\mu _n`$ corresponds to the nonphysical wave of the augmented system.

If $`B_10`$ then the eigenvectors of the fast modes are

$$𝒓_f^\pm =(\eta _f,\zeta _f,0,\mu _f^\pm \zeta _f,\mu _f^\pm \eta _f)^t,$$
(28)
where

$$\eta _f=E_3+\mu _f^\pm B_2,\zeta _f=E_2\mu _f^\pm B_3.$$
(29)
It is shown in Sec.4 that these waves have the same properties as linearly polarised electromagnetic waves in vacuum.

The eigenvectors of the Alfvén modes are

$$𝒓_a^\pm =(\zeta _a,\eta _a,\frac{(\eta _a^2+\zeta _a^2)}{B_1},\mu _a^\pm \eta _a,\mu _a^\pm \zeta _a)^t,$$
(30)
where

$$\eta _a=E_3+\mu _a^\pm B_2,\zeta _a=E_2\mu _a^\pm B_3,$$
(31)
In some applications it might be useful to determine the left eigenvectors as well as the right eigenvectors. Unfortunately, the author was unable to find a concise form for the left eigenvectors and they are therefore not presented here.

Like MHD, the system of FFDE is not strictly hyperbolic, that is it allows multiple eigenvalues under certain conditions. Such degenerate111In this paper three different kinds of degeneracy are encountered: 1) the degeneracy of electromagnetic field itself (see Sec.1), 2) the multiplicity of eigenvalues of the Jacobean matrix (Sec.3), and 3) the so-called linear degeneracy of hyperbolic waves (see Sec.4).  cases require separate treatment because the eigenvectors given by (28) and (30) become either singular or linearly dependent.

1. Provided

$$\stackrel{\mathbf{}}{𝑬}_t=\stackrel{\mathbf{}}{𝒊}_1\mathbf{\times }\stackrel{\mathbf{}}{𝑩},$$
(32)
where $`\stackrel{\mathbf{}}{𝒊}_\mathrm{𝟏}`$ is the unit vector along the $`x_1`$-axis and $`\stackrel{\mathbf{}}{𝑬}_t`$ is the tangential component of electric field, one has

$$\mu _a^+=\mu _f^+=+1.$$

and the corresponding eigenvectors form the two-dimensional vector space generated by

$$(0,1,0,1,0)^t\text{and}(1,0,0,0,1)^t.$$
(33)
Obviously, vacuum electromagnetic waves propagating in the positive direction have the same eigenspace. Thus, in this limit the right Alfvén wave also behaves like a linearly polarised electromagnetic wave in vacuum. Similarly, if

$$\stackrel{\mathbf{}}{𝑬}_t=\stackrel{\mathbf{}}{𝒊}_1\mathbf{\times }\stackrel{\mathbf{}}{𝑩},$$
(34)
one has

$$\mu _a^{}=\mu _f^+=1,$$

and the corresponding eigenvectors form the two-dimensional vector space generated by

$$(0,1,0,1,0)^t\text{and}(1,0,0,0,1)^t.$$
(35)
This eigenspace is the same as that of vacuum electromagnetic waves propagating in the negative direction. We notice that conditions (32) and (34) cannot be satisfied simultaneously unless $`\stackrel{}{𝑬}=0`$ and $`\stackrel{}{𝑩}_t=0`$.
2. If $`B_1=0`$ then

$$\mu _a^+=\mu _a^{}=\mu _a=\frac{E_2}{B^3}=\frac{E_3}{B_2}$$
(36)
with a two-dimensional vector space of eigenvectors. One can use

$$𝒓=(1,0,\frac{B_2(\mu _a^21)}{E_1},0,\mu _a)^t$$
(37)
and

$$𝒓=(0,1,\frac{B_3(\mu _a^21)}{E_1},\mu _a,0)^t$$
(38)
as the basis of the eigenspace if $`E_10`$ or

$$𝒓=(0,0,1,0,0)^t$$
(39)
and

$$𝒓=(B_3,B_2,0,\mu _aB_2,\mu _aB_3)^t$$
(40)
if $`E_1=0`$.

## 4 Properties of Waves of Degenerate Electrodynamics

### 4.1 Fast waves

From (28) one obtains the following system of differential equations describing the fast simple wave:

$$\frac{dB_2}{(E_3+\mu _fB_2)}=\frac{dB_3}{E_2\mu _fB_3}=\frac{dE_1}{0}=\frac{dE_2}{\mu _f(E_2\mu _fB_3)}=\frac{dE_3}{\mu _f(E_3+\mu _fB_2)},$$

where $`\mu _f`$ is either $`\mu _f^+=+1`$ or $`\mu _f^{}=1`$. These are easily integrated to obtain

$$E_1=\text{const}.$$
(41)
$$\eta _f=E_3+\mu _fB_2=\text{const},$$
(42)
$$\zeta _f=E_2\mu _fB_3=\text{const},$$
(43)
$$E_iB^i=\text{const},$$
(44)
If we introduce the vector

$$\stackrel{\mathbf{}}{𝒕}_f=(0,\zeta _f,\eta _f)^t$$

then

$$d\stackrel{\mathbf{}}{𝑬}_t\stackrel{\mathbf{}}{𝒕}_f\text{and}d\stackrel{\mathbf{}}{𝑩}_t\stackrel{\mathbf{}}{𝒕}_f,$$
(45)
where $`\stackrel{\mathbf{}}{𝑩}_t`$ and $`\stackrel{\mathbf{}}{𝑬}_t`$ are the tangential components of the fields. Thus, these are transverse waves and they have the same properties as linearly polarised electromagnetic waves in vacuum.

### 4.2 Alfvén waves

From (30) one obtains the equations of an Alfvén simple wave:

$$\frac{dB_2}{E_2\mu _aB_3}=\frac{dB_3}{E_3+\mu _aB_2}=\frac{B_1dE_1}{\eta _a^2+\zeta _a^2}=\frac{dE_2}{\mu _a(E_3+\mu _aB_2)}=\frac{dE_3}{\mu _a(E_2\mu _aB_3)},$$
(46)
where $`\mu _a`$ is either $`\mu _a^+`$ or $`\mu _a^{}`$. First, we show that

$$B^2E^2=\text{const}$$
(47)
and

$$\mu _a=\text{const}.$$
(48)
Since both these properties are Lorentz invariant it is sufficient to show that they hold in one particular frame. In this case, the most convenient frame is the “fluid” frame, where $`\stackrel{\mathbf{}}{𝑬}=\stackrel{\mathbf{}}{\mathrm{𝟎}}`$ and, hence, $`dE^2=0`$. Moreover, (46) ensures

$$d(B^2)=d(B_3E_2B_2E_3)=0$$

and the results (47,48) are manifest.

Next, equations (46) can easily be integrated to obtain

$$\eta _a=E_3+\mu _aB_2=\text{const},$$
(49)
$$\zeta _a=E_2\mu _aB_3=\text{const},$$
(50)
$$\zeta _aE_2+\eta _aE_3=\text{const},$$
(51)
$$E_iB^i=\text{const},$$
(52)
Since the normal component of electric field varies the Alfvén waves of FFDE are not quite transverse. This is similar to what has been found for Alfvén waves in relativistic MHD \[Komissarov 1997\]. These waves can still be described as linearly polarised. If we introduce the vector

$$\stackrel{\mathbf{}}{𝒕}_a=(0,\zeta _a,\eta _a)^t$$

then

$$d\stackrel{\mathbf{}}{𝑬}_t\stackrel{\mathbf{}}{𝒕}_a\text{and}d\stackrel{\mathbf{}}{𝑩}_t\stackrel{\mathbf{}}{𝒕}_a.$$
(53)
From this equation and (45) one can see that the polarisations of Alfvén and electromagnetic waves propagating in the same direction are not, in general, mutually orthogonal.

If $`B_1=0`$ then from (26) it follows that Alfvén waves propagate with the drift velocity

$$\mu _a=(\stackrel{}{E}\times \stackrel{}{B})_1/B^2,$$

and, therefore, must be somewhat similar to the contact wave of relativistic MHD. Their properties are most transparent in the wave frame, where $`\mu =0`$. From (36) one can see that in such frame

$$\stackrel{\mathbf{}}{𝑬}_t=\stackrel{\mathbf{}}{\mathrm{𝟎}}.$$
(54)
Using this result and (37) one finds that

$$(B_1,B_2,0,0,0)^t\text{and}(E_1,0,B_2,0,0)^t$$
(55)
is a basis of the Alfvén eigenspace. The first vector describes the familiar rotation of $`\stackrel{\mathbf{}}{𝑩}_t`$. The second allows variation of the magnitude of $`\stackrel{\mathbf{}}{𝑩}_t`$, but requires

$$B_t^2E_1^2=\text{const}.$$

Equations (25,48) show that all waves of FFDE are linearly degenerate, that is

$$𝒓\mathbf{}_U\mu =0,$$

and, therefore, this system does not allow the formation of shocks via wave steepening. However, discontinuities may be introduced via discontinuous initial data and/or boundary conditions. For linearly degenerate waves the shock equations coincide with jump equations of simple waves \[Boillat 1982, Anile 1989\].

Finally, figure 1 shows the normal wave speed diagram in the “fluid” frame, where $`\stackrel{}{E}=\stackrel{}{0}`$. Since in this frame

$$\mu _a^\pm =\pm \mathrm{cos}\theta ,$$

where $`\theta `$ is the angle between $`\stackrel{\mathbf{}}{𝑩}`$ and the wave vector, all curves on this diagram are circles.

## 5 1D Upwind Numerical Scheme

Nowadays, the construction of upwind numerical schemes for hyperbolic conservation laws is a well established area of numerical analysis. Here we describe a simple one-dimensional scheme which is very similar to the one constructed by Komissarov \[Komissarov 1999\] for RMHD. The space-time is considered to be flat.

### 5.1 General Structure

The one-dimensional conservation laws of FFDE can be written in the form

$$\frac{𝑸}{t}+\frac{𝑷}{x^1}=0,$$
(56)
where the vector of conserved variables is

$$𝑸=(S_1,S_2,S_3,B_2,B_3)^t,$$

where $`\stackrel{\mathbf{}}{𝑺}`$ is the Poynting flux vector, and the vector of hyperbolic fluxes is

$$𝑷=(T_{11},T_{12},T_{13},E_3,E_2)^t.$$

Notice that we use all three space components of the energy-momentum equation. In addition, we introduce the auxiliary vector of “primitive” variables

$$𝑼=(B_2,B_3,E_1,E_2,E_3)^t.$$

These primitive variables can be converted into the conservative ones and vice versa via

$$\stackrel{\mathbf{}}{𝑺}=\stackrel{\mathbf{}}{𝑬}\times \stackrel{\mathbf{}}{𝑩},$$
(57)
$$\stackrel{\mathbf{}}{𝑬}=\frac{1}{B^2}\stackrel{\mathbf{}}{𝑺}\times \stackrel{\mathbf{}}{𝑩}.$$
(58)
The second equation automatically ensures the degeneracy condition (17). Unfortunately, this does not mean that the second condition, (18), will also be satisfied and extra care would have to be taken in cases where $`B^2E^2`$ is close to zero. At the moment we do not have the means to ensure (17).

Let us define a regular Cartesian grid such that the $`i`$th cell is centered at $`x_i=ih`$ and occupies the region between $`x_{i1/2}=(i1/2)h`$ and $`x_{i+1/2}=(i+1/2)h`$, where $`h`$ is the mesh spacing. Now suppose that we know the solution at $`t=t_n`$ and we want to calculate it at a later time $`t_{n+1}`$. We integrate equations (56) over the $`i`$th cell and from $`t=t_n`$ to $`t=t_{n+1}=t_n+\mathrm{\Delta }t`$ to get

$$𝑸_{i,n+1}𝑸_{i,n}+\frac{\mathrm{\Delta }t}{h}(𝑷_{i+\frac{1}{2},n+\frac{1}{2}}𝑷_{i\frac{1}{2},n+\frac{1}{2}})=0.$$
(59)
Here

$$𝑸_{i,n}=\frac{1}{h}\underset{x_{i\frac{1}{2}}}{\overset{x_{i+\frac{1}{2}}}{}}𝑸(x,t_n)𝑑x$$

is the mean value of $`𝑸`$ in the $`i`$th cell at time $`t_n`$ and

$$𝑷_{i+\frac{1}{2},n+\frac{1}{2}}=\frac{1}{\mathrm{\Delta }t}\underset{t_n}{\overset{t_{n+1}}{}}𝑷(x_{i+\frac{1}{2}},t)𝑑t,$$

is the time average of the flux at the cell interface with $`x=x_i+\frac{1}{2}h`$.

In a first order Godunov-type scheme (Godunov 1959), it is assumed that at the beginning of each time step the solution is uniform within each cell. This implies initial discontinuities at the cell interfaces and the time averaged fluxes can therefore be found by solving the corresponding Riemann problems. In fact, the solution, $`𝑼_{i+1/2}^{}`$, of the problem at the interface with $`x=x_{i+1/2}`$, does not depend on time but only on the initial left, $`𝑼_{i,n}`$, and right, $`𝑼_{i+1,n}`$, states. Then the first order fluxes can be computed from

$$𝑷_{i+\frac{1}{2},n+\frac{1}{2}}=𝑷(𝑼_{i+\frac{1}{2}}^{}).$$

In order to achieve second order accuracy, we use the first order scheme to obtain the solution, $`𝑸_{i,n+\frac{1}{2}}`$ and hence $`𝑼_{i,n+\frac{1}{2}}`$ at the half time, $`t_{n+\frac{1}{2}}=t_n+\mathrm{\Delta }t/2`$. We then use this to compute average gradients of primitive variables in each cell as follows

$$\left(\frac{𝑼}{x}\right)_{i,n+\frac{1}{2}}=\frac{1}{h}\mathrm{av}(\mathrm{\Delta }𝑼_{i,n+\frac{1}{2}},\mathrm{\Delta }𝑼_{i+1,n+\frac{1}{2}}),$$

where

$$\mathrm{\Delta }𝑼_{i,n+\frac{1}{2}}=\left(𝑼_{i,n+\frac{1}{2}}𝑼_{i1,n+\frac{1}{2}}\right),$$

and $`\mathrm{av}(a,b)`$ is a non-linear averaging function whose purpose is to reduce the scheme to first order in space in the neighborhood of discontinuities. Here we will adopt the same averaging function as in \[Falle 1991\]:

$$\mathrm{av}(a,b)=\{\begin{array}{cc}\frac{(a^2b+ab^2)}{(a^2+b^2)}\hfill & \mathrm{if}ab0,\text{ and }a^2+b^20,\hfill \\ & \\ 0\hfill & \mathrm{if}ab<0\text{ or }a^2+b^2=0,\hfill \end{array}$$

These gradients can now be used to set up the left and right states for the second order Riemann problems

$$𝑼_{i+1/2}^l=𝑼_{i,n+\frac{1}{2}}+\frac{1}{2}h\left(\frac{𝑼}{x}\right)_{i,n+\frac{1}{2}},$$

$$𝑼_{i+1/2}^r=𝑼_{i+1,n+\frac{1}{2}}\frac{1}{2}h\left(\frac{𝑼}{x}\right)_{i+1,n+\frac{1}{2}}.$$

The solution, $`𝑼_{i+1/2}^{}=𝑼(𝑼_{i+1/2}^l,𝑼_{i+1/2}^r)`$, to this Riemann problem allows us to compute the second order fluxes

$$𝑷_{i+\frac{1}{2},n+\frac{1}{2}}=𝑷(𝑼_{i+\frac{1}{2}}^{}),$$

which are used to advance the solution through the full time step from $`t=t_n`$ to $`t=t_{n+1}`$ according to (59).

### 5.2 Riemann Solvers

By a Riemann problem we mean the initial value problem for (56) with the initial conditions

$$𝑼(x,0)=\{\begin{array}{cc}𝑼^l& \text{for}x0,\\ 𝑼^r& \text{for}x>0.\end{array}$$
(60)
For a general hyperbolic system this problem has a self-similar solution which involves only shocks and centered rarefactions (e.g. Landau & Lifshitz 1959, Jeffrey & Taniuti 1964). Since all waves of FFDE are linearly degenerate, centered simple waves do not exist and, thus, the solution involves only discontinuities.

#### 5.2.1 Linear Riemann Solver

Since all waves are linearly degenerate we have

$$𝑼_r=𝑼_l+\underset{i=1,4}{}_i𝒓_i,$$
(61)
where $`𝒓_i`$ is the eigenvector of the upstream state of ith wave. For fast waves these states are already known, they are the left and the right states of the Riemann problem. Upstream states of Alfvén waves are different and depend on the amplitudes of fast waves. In our linear Riemann solver we ignore this difference. Then (61) becomes a linear system of five equations for four unknowns, $`_i`$. To handle it, we ignore one of the equations. The corresponding four components of the solution at $`x^1=0`$ are then found via

$$𝑼^{}=𝑼_l+\underset{\mu _i<0}{}_i𝒓_i=𝑼_r\underset{\mu _i>0}{}_i𝒓_i.$$
(62)
The remaining component is determined using the degeneracy condition (17).

#### 5.2.2 Exact Riemann Solver

In most Riemann problems originated from numerical simulations the difference between left and right states is rather small and linear Riemann solvers are sufficiently accurate. However, for handling strong discontinuities one might require an exact solution of the Riemann problem. Thanks to the simplicity of FFDE, an exact Riemann solver can be constructed rather easily.

First we consider the case where none of the states of the Riemann problem has multiple eigenvalues. From the results of Sec.4 it follows that only fast waves change the value of $`p=B^2E^2`$ and only Alfvén waves change the value of $`E_1`$. This implies that the Riemann problem reduces to determining the values of $`p`$ and $`E_1`$ in the region bounded by Alfvén waves, $`\overline{p}`$ and $`\overline{E}_1`$. Given the left state, the values of $`E_1`$ and $`p`$ in the right state, and $`\overline{p}`$ and $`\overline{E}_1`$ the corresponding right state $`𝑽_r(\overline{p},\overline{E}_1)`$ can be found by successive application of the jump equations from left to right:

$$\frac{[B_2]}{\eta _f}=\frac{[B_3]}{\zeta _f}=\frac{[E_1]}{0}=\frac{[E_2]}{\mu _f\zeta _f}=\frac{[E_3]}{\mu _f\eta _f}=\frac{[p]}{2\mu _f(\eta _f^2+\zeta _f^2)},$$
(63)
for fast waves, and

$$\frac{[B_2]}{\zeta _a}=\frac{[B_3]}{\eta _a}=\frac{B_1[E_1]}{\eta _a^2+\zeta _a^2}=\frac{[E_2]}{\mu _a\eta _a}=\frac{[E_3]}{\mu _a\zeta _a}=\frac{[p]}{0}$$
(64)
for Alfvén waves. This gives us the system of nonlinear equations

$$𝑽_r(\overline{p},\overline{E}_1)=𝑼_r.$$
(65)
which can then be solved iteratively. Once $`\overline{p}`$ and $`\overline{E}_1`$ are found the resolved state is determined via (63,64).

In degenerate cases the problem is even easier.

1. If $`\mu _f^+=\mu _a^+`$ in the right state of the Riemann problem, then the right fast and Alfvén waves merge into a single degenerate wave across which both $`p`$ and $`E_1`$ are constant Thus we immediately obtain $`\overline{p}=p_r`$, $`\overline{E}_1=E_{1_r}`$ and then determine the resolved state in the same way as in the non-degenerate case. Similarly, if $`\mu _f^{}=\mu _a^{}`$ in the left state then $`\overline{p}=p_l`$, $`\overline{E}_1=E_{1_l}`$.
2. It is easy to verify that for any of the degenerate waves, the corresponding degeneracy condition holds on both sides of the wave. Thus, if $`\mu _f^+`$ equals to $`\mu _a^+`$ in the right state and $`\mu _f^{}`$ equals to $`\mu _a^{}`$ in the left state then in the resolved state conditions (32) and (34) are satisfied simultaneously. This immediately gives $`\stackrel{}{𝑬}=0`$ and $`\stackrel{}{𝑩}_t=0`$ in the resolved state.
3. If $`B_1=0`$ then the Alfvén waves merge into one contact-like wave across which both $`p`$ and $`\mu _a`$ are constant. This leads to the following two equations:

$$p_l^{}(_{f_l})=p_r^{}(_{f_r}),$$
(66)
$$\mu _{a_l}^{}(_{f_l})=\mu _{a_r}^{}(_{f_r}),$$
(67)
where * indicates states adjacent to the contact wave. Simple manipulations reduce this system to a quadratic equation for $`_{e_l}`$ (or $`_{e_r}`$). It appears that only one of its solutions satisfies the condition (18).

### 5.3 Test calculations

We have tested the scheme against exact analytic solutions for plane electromagnetic and Alfvén waves including all degenerate cases. The exact solutions can easily be constructed using the results of Sec.3 and 4. The profiles of these waves should remain unchanged as they propagate. Figures 23 illustrate the results for some of the test calculations. In these figures the initial and the final exact solutions are shown by continuous lines and markers show the final numerical solution.

Fast wave: Figure 2a shows the results for a fast wave propagating in the positive direction. The initial solution is

$$B_1=1.0,B_3=E_2=0.0,$$

$$B_2=\{\begin{array}{cc}1.0\hfill & \text{if}x^1<0.1\hfill \\ 1.0+\frac{3}{2}(x^1+0.1)\hfill & \text{if}0.1<x^1<0.1\hfill \\ 1.3\hfill & \text{if}x^1>0.1\hfill \end{array}$$

Alfvén wave: Figure 2b shows the results for a non-degenerate Alfvén wave propagating in the negative direction with the wave speed $`\mu =0.5`$. The initial solution is

$$B_1^{}=B_2^{}=1.0,E_2^{}=E_3^{}=0.0,$$

$$B_3^{}=\{\begin{array}{cc}1.0\hfill & \text{if}x^1<0.1\hfill \\ 1.0+\frac{3}{2}(x^1+0.1)\hfill & \text{if}0.1<x^1<0.1\hfill \\ 1.3\hfill & \text{if}x^1>0.1\hfill \end{array},$$

where $`\stackrel{}{B}^{}`$ and $`\stackrel{}{E}^{}`$ are measured in the wave frame which is moving relative to the grid with speed $`\mu =0.5`$.

Degenerate Alfvén wave: Figure (3) shows the results for a rotational degenerate Alfvén wave propagating to the right with $`\mu =0.5`$. The initial solution is

$$\stackrel{}{E}^{}=0,B_1^{}=0,$$

$$B_2^{}=2\mathrm{cos}\varphi ,B_3^{}=2\mathrm{sin}\varphi ,$$

where

$$\varphi =\{\begin{array}{cc}0.0\hfill & \text{if}x^1<0.1\hfill \\ \frac{5\pi }{2}(x^1+0.1)\hfill & \text{if}0.1<x^1<0.1\hfill \\ \frac{\pi }{2}\hfill & \text{if}x^1>0.1\hfill \end{array}.$$

Three waves problem: In this last test the initial discontinuity at $`x=0`$ splits into two fast discontinuities and a stationary Alfvén discontinuity. The initial solution is

$$\begin{array}{ccc}\stackrel{}{𝑩}=(1.0,1.5,3.5)\hfill & \stackrel{}{𝑬}=(1.0,0.5,0.5)\hfill & \text{if}x^1<0,\hfill \\ \stackrel{}{𝑩}=(1.0,2.0,2.(3))\hfill & \stackrel{}{𝑬}=(1.5,1.(3),0.5)\hfill & \text{if}x^1>0.\hfill \end{array}$$

As one can see, numerical diffusion smoothes out the sharp wave fronts of the exact solutions but otherwise the agreement is very good indeed.

Not all Riemann problems with both the left and the right states satisfying the degeneracy conditions (17,18) have solutions within the frame of degenerate electrodynamics. For example, the following problem does not seem have one

$$\begin{array}{ccc}\stackrel{}{𝑩}=(1.0,1.0,1.0)\hfill & \stackrel{}{𝑬}=(0.0,0.5,0.5)\hfill & \text{if}x^1<0,\hfill \\ \stackrel{}{𝑩}=(1.0,1.0,1.0)\hfill & \stackrel{}{𝑬}=(0.0,0.5,0.5)\hfill & \text{if}x^1>0,\hfill \end{array}$$

as iterations of our exact Riemann solver do not converge. To make sure that this is not due to a fault in the solver we have considered a related problem where this initial discontinuity is substituted with a linear transition layer. Numerical simulations (see figure 5) show that within this layer $`B^2E^2`$ decreases in time toward zero leading to the breakdown of the condition

$$B^2E^2>0.$$

As $`B^2E^2`$ tends to zero the drift velocity of plasma

$$\stackrel{}{𝒗}_d=\stackrel{}{𝑬}\times \stackrel{}{𝑩}/B^2$$

tends to the speed of light suggesting that one may no longer neglect the inertia of plasma particles.

## 6 Conclusions

We have shown that the evolution equations of FFDE can be written in the form of conservation laws. One group of equations arise from the energy-momentum conservation. There are only two linear independent dynamic equations (components of one tensor equation) in this group. This is in contrast to MHD, where all four components of the energy-momentum equation are linearly independent, and vacuum electrodynamics, where in general $`𝑭_{\mu \nu }`$ does not have zero eigenvectors and the energy-momentum equation is reduced to

$$_\nu F^{\mu \nu }=0,$$

which has three independent dynamic components. The one-dimensional system of FFDE with plane symmetry, in flat space-time, includes four linearly independent dynamic equations for four independent components of the electromagnetic field tensor (other equations simply require $`B_1=`$const.) This system is hyperbolic, though not strictly hyperbolic, with a pair of fast waves and a pair of Alfvén waves. All characteristic modes are linearly degenerate, which means that this system does not allow the formation of shocks via steepening of continuous waves.

This formulation is particularly useful when it comes to developing numerical schemes for time-dependent FFDE as there has been a great deal of work on numerical methods for hyperbolic conservation laws. Here we have presented a simple one-dimensional scheme based mainly on our linear Riemann solver. Its generalization to multi-dimensions is relatively straightforward. We have already constructed a 2D numerical scheme adapted to Kerr space-time which will be described elsewhere. Such schemes can be useful tools in studying electromagnetically driven jets and winds from black holes and neutron stars and, perhaps, other related phenomena of relativistic astrophysics.

## 7 Acknowledgments

This research is partly supported by PPARC. The author thanks S.A.E.G. Falle for useful discussions and comments on the first version of this manuscript.


text:
# Calculation of 𝚫(𝐤,𝝎) for a 2D t-J Cluster

## Abstract

Using numerical techniques, the diagonal and off-diagonal superconducting one-electron Green’s functions are calculated for a two-dimensional (2D) t-J model on a periodic 32-site cluster at low doping. From these Green’s functions, the momentum and frequency dependence of the pairing gap $`\mathrm{\Delta }(𝐤,\omega )`$ are extracted. It has $`d_{x^2y^2}`$ symmetry and exhibits $`\omega `$-dependent structure which depend upon J/t. We find that the pairing gap persists down to small J/t values. The frequency- and momentum-dependent renormalized energy and renormalization factor are also calculated.

For the last decade, the search for superconductivity in models of strongly-correlated fermions, has been triggered by accumulating experimental evidence in favor of an unconventional (i.e. not based on the usual phonon-mediated interaction) mechanism in high-$`T_c`$ two-dimensional (2D) TK00  or ladder DR96  superconducting cuprates. Soon after the discovery of the layered cuprates, it was proposed that the Hubbard model And87  and its strong coupling limit, the t-J model ZR88 , captured the generic features of these materials. Early RPA MSV86 ; SLH86  and Gutzwiller variational calculations GJR87  supported the notion that the doped state would have $`d_{x^2y^2}`$ pairing. Nevertheless, despite the conceptual simplicity of these Hamiltonians and the results of these approximate calculations, the nature of the basic mechanism responsible for pairing as well as the actual physical properties of these models remained controvertial. In recent years, numerical calculations have provided insight into the second of these questions regarding the actual properties of these models. Specifically, Monte Carlo calculations showed that the ground states of the undoped models had long-range antiferromagnetic order HT89 ; RRY89 . Exact diagonalization studies showed that holes doped into a t-J cluster can form $`d_{x^2y^2}`$ pairs Poi93 ; Leu00 . In addition, numerical calculations support the view that the doped system can have low-lying stripe domain wall states SW00  which can be stabilized by lattice anisotropies KSW01 ; BCS01 . While the static striped phase competes with superconductivity, it has been argued that the addition of a next-near-neighbor hopping WS99  or the use of periodic boundary conditions Sor01  tips the balance in favor of the $`d_{x^2y^2}`$ pairing phase. The point is that a variety of numerical calculations provide evidence that indeed these models do exhibit the basic properties seen in the cuprates. However, in spite of this progress, the question regarding the nature of the basic mechanism responsible for the pairing remains open. Within the traditional BCS framework, one would look for a reflection of the pairing interaction by examining the momentum and frequency dependence of the gap $`\mathrm{\Delta }(𝐤,\omega )`$. Because of the relatively short coherence length and the relatively high energy scales J and t, one can hope to learn about both the momentum and frequency dependence of $`\mathrm{\Delta }(𝐤,\omega )`$ from an exact diagonalization study of a cluster. Here we present a numerical study aimed at doing this for the t-J model and obtain the first results for the $`\omega `$-dependence of the gap.

The 2D t–J Hamiltonian reads,

$`H`$ $`=`$ $`{\displaystyle \underset{i,j}{}}J_{ij}(𝐒_i𝐒_j{\displaystyle \frac{1}{4}}n_in_j)`$ (1)
$`+`$ $`{\displaystyle \underset{i,j,\sigma }{}}t_{ij}(\stackrel{~}{c}_{i,\sigma }^{}\stackrel{~}{c}_{j,\sigma }+h.c.)`$
where the exchange integrals $`J_{ij}`$ and hopping terms are (for simplicity) restricted to nearest neighbor (NN) sites (called hereafter $`J`$ and $`t`$) and $`\stackrel{~}{c}_{i,\sigma }^{}`$ are projected fermion creation operators defined as $`c_{i,\sigma }^{}(1n_{i,\sigma })`$. Hereafter, unless specified otherwise, $`t`$ sets the energy scale. The conclusions drawn in this Letter are supported by Exact Diagonalisation (ED) studies performed on the square cluster of N=32 sites depicted in Fig. 1(a) slightly doped with up to 2 holes note1 . This finite-size system is particularly appealing since it exhibits the full local symmetries of the underlying square lattice as well as all the most symmetric $`𝐤`$-points in reciprocal space as seen in Fig. 1(b). Note that, although the hole doping is quite small, the hole occupation shown in Fig. 1(b) for a few $`𝐤`$-points is rather consistent with a ”large” Fermi surface.

A superconducting ground state (GS) is characterized by Gorkov’s off-diagonal one-electron time-ordered Green function

$$F_𝐤(t)=T\stackrel{~}{c}_{𝐤,\sigma }(t/2)\stackrel{~}{c}_{𝐤,\sigma }(t/2).$$
(2)
Close to half-filling, in a finite system, $`F_𝐤`$ can be computed from,

$$\stackrel{~}{F}_𝐤(z)=N2|\stackrel{~}{c}_{𝐤,\sigma }\frac{1}{zH+\overline{E}_{N1}}\stackrel{~}{c}_{𝐤,\sigma }|N,$$
(3)
defined for all complex $`z`$ (with $`\mathrm{Im}z0`$). Here the number of particles in the initial $`|N>`$ (half-filling) and final $`|N2>`$ (two-hole doped) GS differ by two, reflecting charge fluctuations in a SC state. For convenience, the energy reference $`\overline{E}_{N1}`$ is defined as the average between the GS energies of $`|N>`$ and $`|N2>`$. Note that these states are both spin singlets. In addition, they exhibit s-wave and d$`_{x^2y^2}`$ orbital symmetries (both of even parity) respectively. Owing to these special features, it is straightforward to show that the analytical continuation of $`F_𝐤(\omega )`$ to the real frequency axis is even in frequency and can be expressed as,

$$F(𝐤,\omega )=\stackrel{~}{F}_𝐤(\omega +iϵ)+\stackrel{~}{F}_𝐤(\omega +iϵ),$$
(4)
where $`ϵ`$ is a small imaginary part. The superconducting frequency-dependent gap function $`\mathrm{\Delta }(𝐤,\omega )`$ is directly proportional to $`F(𝐤,\omega )`$ and a simple analysis shows that both functions have d$`_{x^2y^2}`$ orbital symmetry. In particular, they identically vanish on the diagonals $`k_x=\pm k_y`$.

Pairing between holes should also be reflected in the structure of the (time-ordered) diagonal Green function $`G(𝐤,\omega )`$. Here, it is convenient to define the finite size $`G(𝐤,\omega )`$ as the sum of the following electron- (i.e. occupied states for $`\omega <0`$) and hole-like (i.e. empty states for $`\omega >0`$) parts,

$`G(𝐤,\omega )=N|\stackrel{~}{c}_{𝐤,\sigma }^{}{\displaystyle \frac{1}{\omega iϵ+H\overline{E}_{N1}}}\stackrel{~}{c}_{𝐤,\sigma }|N`$ (5)
$`+`$ $`<N2|\stackrel{~}{c}_{𝐤,\sigma }{\displaystyle \frac{1}{\omega +iϵH+\overline{E}_{N1}}}\stackrel{~}{c}_{𝐤,\sigma }^{}|N2>,`$
so that both Green functions $`F`$ and $`G`$ have the same set of energy poles. Here, the well-known continued-fraction method used to compute diagonal correlation functions such as $`G(𝐤,\omega )`$ is extended to deal with off-diagonal ones such as $`F(𝐤,\omega )`$ OSEM94 ; continued\_frac . Data are shown in Fig. 2(a-d). The spectral density $`\mathrm{Im}\mathrm{G}(𝐤,\omega )`$ exhibits sharp quasi-particle-like peaks both above and below the Fermi level ($`\omega =0`$). Contrary to $`G(𝐤,\omega )`$, $`F(𝐤,\omega )`$ has a significant amplitude only at low energy, typically for $`|\omega |<4J`$ as seen in Fig. 2(b), reflecting the energy scale of the pairing interaction. Note that, due to the discreteness of the low energy spectrum of the cluster, a finite value of $`ϵ`$ is necessary to wipe out the irrelevant fast oscillations of the $`F`$ (or $`G`$) Green functions (see e.g. Fig. 2(c-d)). However, the static limit $`\omega 0`$ (together with $`ϵ0`$) is perfectly controled and has a physical meaning.

The frequency-dependent gap function $`\mathrm{\Delta }(𝐤,\omega )`$ can be extracted from the knowledge of the Green functions by assuming generic forms of a SC GS schrieffer  at low energies,

$`G(𝐤,\omega )`$ $`=`$ $`z(𝐤,\omega ){\displaystyle \frac{\omega +ϵ(𝐤,\omega )}{\omega ^2(ϵ(𝐤,\omega )^2+\mathrm{\Delta }(𝐤,\omega )^2)+iϵ}},`$
$`F(𝐤,\omega )`$ $`=`$ $`z(𝐤,\omega ){\displaystyle \frac{\mathrm{\Delta }(𝐤,\omega )}{\omega ^2(ϵ(𝐤,\omega )^2+\mathrm{\Delta }(𝐤,\omega )^2)+iϵ}},`$ (6)
where $`z(𝐤,\omega )`$ and $`ϵ(𝐤,\omega )`$ are the inverse renormalization parameter and the renormalized energy (containing the energy shift) respectively. In the BCS limit, the renormalized energy $`ϵ(𝐤,\omega )ϵ_𝐤`$ and the gap function $`\mathrm{\Delta }(𝐤,\omega )\mathrm{\Delta }_𝐤`$ do not depend on frequency (and $`z(𝐤,\omega )=1`$) so that $`F(𝐤,\omega )`$ exhibits only two poles at $`\pm (ϵ_𝐤^2+\mathrm{\Delta }_𝐤^2)^{1/2}`$. However, it is necessary to assume an explicitely frequency dependence of the gap function in order to reproduce the secondary peaks seen in Fig. 2(d). Using the parity of $`ϵ(𝐤,\omega )`$ and $`\mathrm{\Delta }(𝐤,\omega )`$ in the vicinity of the real frequency axis one gets,

$$\mathrm{\Delta }(𝐤,\omega )=\frac{2\omega F(𝐤,\omega )}{G(𝐤,\omega )G(𝐤,\omega )}$$
(7)
Thus, from a numerical calculation of $`G(𝐤,\omega )`$ and $`F(𝐤,\omega )`$, one can obtain $`\mathrm{\Delta }(𝐤,\omega )`$.

Results for $`\mathrm{\Delta }(𝐤,\omega )`$ with $`𝐤=(\pi ,0)`$ at the $`M`$ point, are shown in Fig. 3(a). Despite the limited resolution in frequency of the ED data, $`\mathrm{\Delta }(𝐤,\omega )`$ exhibits dynamic structure on a scale of energies several times J. We believe that finite-size effects have increased the onset frequency $`\omega _\mathrm{\Delta }`$ and that the corresponding time $`2\pi /\omega _\mathrm{\Delta }`$ should be viewed as a lower bound of the characteristic time scale of the pairing interaction. For comparison, the dynamic spin structure factor $`Im\chi (𝐪,\omega )`$ for the 2-hole doped 32-site cluster with $`𝐪=(\pi ,\pi )`$ is shown in Fig. 3(b).

The quasi-particle gap at the gap edge, where $`ϵ(𝐤,0)=0`$, is determined from the self-consistent solution of

$$Re\mathrm{\Delta }(𝐤,\omega =\mathrm{\Delta }_𝐤)=\mathrm{\Delta }_𝐤$$
(8)
As shown in Fig. 3(a), $`\mathrm{\Delta }(𝐤,\omega )`$ is real and essential constant over an energy region larger than the gap. Furthermore, as we will see, $`ϵ(𝐤,0)`$ for $`𝐤=(\pi ,0)`$ is small so that $`\mathrm{\Delta }_𝐤`$ is to a good approximation given by $`\mathrm{\Delta }(𝐤,0)`$. For $`𝐤=(\pi ,0)`$, the zero frequency limit of the gap function $`\mathrm{\Delta }(𝐤,\omega 0)`$ (which is purely real) computed by numerically taking the $`\omega `$ and $`ϵ0`$ limits of eq. (7) note2  is plotted versus J in Fig. 4(a). In Fig. 4(b) the static $`\omega =0`$ value of the renormalization energy $`ϵ(𝐤,0)`$, the renormalization factor $`z(𝐤,\omega =0)`$, and the equal-time pair amplitude $`\stackrel{~}{c}_{𝐤,\sigma }^{}\stackrel{~}{c}_{𝐤\sigma }^{}`$ for $`𝐤=(\pi ,0)`$ versus J are also shown note2\_bis . The values of $`\mathrm{\Delta }(k,0)`$ and $`z(k,0)`$ are in good agreement with the results Ohta, et. al OSEM94  found by fitting their numerical t-J $`4\times 4`$ and $`\sqrt{18}\times \sqrt{18}`$ cluster calculations to a BCS quasi-particle form in which the frequency dependence of $`\mathrm{\Delta }`$ and $`z`$ were neglected. They found that the magnitude of the gap at the antimode varies as 0.25J to 0.3J at low doping and as seen in Fig. 4(a) we find $`\mathrm{\Delta }(k,0)0.3J`$. It is interesting to note that the data for $`\mathrm{\Delta }(k,0)`$ in Fig. 4(a) do not show any lower bound of J in contrast to the pair binding energy Poi93 ; Leu00 . The spin structure factor is plotted in Fig. 3(b) showing a low energy peak (magnon) and a higher energy background that could be connected to the structures in the gap function at an energy $`J`$. For comparison, the magnon energy and the mean-value note3  of the spin spectral weight are plotted in Fig. 4(a) together with the zero-frequency gap and the onset frequency $`\omega _\mathrm{\Delta }`$ in the gap function. While the similarities in the J-dependence of these quantities suggest that the dynamics underlying the pairing mechanism is related to the spin fluctuations, further work on other clusters, such as 2-leg ladders, are needed to sort out the influence of finite-size effects.

We conclude this investigation with a further discussion of the parameters which enter eq. (6). Note that since these parameters depend upon both $`𝐤`$ and $`\omega `$, the static $`\omega =0`$ results which we will discuss are near the energy shell only for $`𝐤`$ values near the fermi surface and for $`\mathrm{\Delta }_𝐤`$ small compared to the characteristic frequency variations set by J. In fact, even trying to select $`𝐤`$-points near the “Fermi surface” defined by $`ϵ(𝐤,0)=0`$ requires some care since the Fermi surface (FS) cannot be exactly defined on a finite cluster, especially for a small number of holes. However, the smooth variation of the hole occupancy in the BZ \[see Fig. 1(b)\] suggests that this system should still pick up features of a slightly doped AF with a large FS. With this in mind, the static limit of the renormalized energy $`ϵ(𝐤,0)`$ and the renormalization parameter $`z(𝐤,0)`$ were obtained from the calculated $`ϵ,\omega 0`$ limits of $`G(𝐤,0)`$, $`F(𝐤,0)`$, and $`\mathrm{\Delta }(𝐤,0)`$ by using eq. (6) for $`\omega =0`$. These are tabulated for values of $`𝐤`$ both near the nominal FS and away from it in Table 1. The renormalized energy $`ϵ(𝐤,0)`$ has roughly the energy scale of the quasi-particle-like energy poles at the available momenta of the BZ. The small value of $`ϵ(𝐤,0)`$ given by our analysis for $`𝐤=(\pi ,0)`$ show that the excitation that we are probing is relatively close to the FS (defined here as having $`ϵ(𝐤,0)=0`$. Note also that the weight $`z(𝐤,0)`$ of these excitations is small $`(J/2)`$ compared to the large incoherent background \[see Fig. 2(a)\]. Its J-dependence plotted in Fig. 4(b) is consistent with the $`J^\alpha `$ power-law behavior found previously in earlier calculations of a single hole propagating in an AF background PSZ93 .

###### Acknowledgements.

D.J. Scalapino would like to acknowledge support from the US Department of Energy under Grant No DE-FG03-85ER45197. D. Poilblanc thanks M. Sigrist (ETH-Zürich) for discussions and aknowledges hospitality of the Physics Department (UC Santa Barbara) where part of this work was carried out. Numerical computations were done on the vector NEC-SX5 supercomputer at IDRIS (Paris, France).