text:
# Dispersion of the dielectric function of a charge-transfer insulator

## I Introduction

The importance of many body effects for the description of dielectric response of insulating solids is generally accepted. In the neighborhood of the fundamental absorption threshold bound exciton lines and continuum excitons drastically change the spectrum. Recently it has been realized that the excitons in charge-transfer insulators (CTI) possess unusual features connected with the strongly correlated character of the antiferromagnetic (AFM) ground state in these compounds. The most pronounced peculiarity consists in the existence of exciton branches with substantially larger dispersion compared with the one-particle excitations. Wang96 ; Fink01  This behavior has a formal analogy with that of the Frenkel exciton Frenkel  that acquires a finite effective mass although both, electron and hole, have infinite masses, Egri  but in CTI it has a completely different origin.

The experimental technique suitable for the observation of the exciton dispersion is the electron energy loss spectroscopy (EELS).Pines  What is actually measured in transmission EELS experiments is the partial cross sectionFink89 ; Schnat79  that may be decomposed into an amplitude factor and a dynamic structure factor

$$\frac{d^2\sigma }{d\mathrm{\Omega }dE}=\frac{4}{(a_0)^2q^4}S(𝐪,\omega ).$$

The dynamic structure factor characterizes the linear response of the whole electronic system on longitudinal electric fields with the momentum $`𝐪`$ and frequency $`\omega `$ ( the ionic contribution may be neglected for the considered frequency range of the order of several eV). Pronounced peaks in $`S(𝐪,\omega )`$ are called excitons. remark  They may correspond to discrete lines in the excitation spectrum of the solid or to resonances in the continuum part of the spectrum.

For the theoretical description of excitonic features in conventional semiconductors and insulators the following scheme is used (The key references are Refs. Egri, ; Knox, ; HankeShm, ; AbInit98, ; Prokof01, ). First, the quasiparticle excitation spectrum is found empirically Egri ; Knox ; Prokof01  or from first principles. HankeShm ; AbInit98  It is essential that the spectrum of the $`N`$-electron system consists mainly of a continuum of electron-hole pairs whose electron or hole quasiparticle excitations are close to eigenstates of the $`(N+1)`$\- or $`(N1)`$-electron system, respectively, with definite quasimomenta and energy. ShRice  The quasiparticle spectrum is usually obtained from the self-consistent-field (SCF) approach. Next, the electron-hole interaction is taken into account. Then the problem for two quasi-particles interacting via the medium is solved. It is crucial that the ground state may be viewed as an occupied valence band that is separated from the first excited state by an energy gap.

For the CTI the above scheme should be revised beginning from the first step. The CTI has an odd number of electrons per formula unit. That is why the SCF calculations usually give a metallic ground state for the CTI and a gapless excitation spectrum. More elaborate SCF methods like the unrestricted Hartree-Fock method or its modern version LDA+U predict an AFM long-range order in the ground state and a gap in the excitation spectrum. Nevertheless, the nature of excitations and the ground state fluctuations are not caught by these approaches. Let us recall that any CTI remains insulating and shows excitonic features in optics and EELS spectra above the Néel temperature, i.e. in the paramagnetic state (with short-range AFM correlations). For strongly correlated systems the electron-electron interaction should be taken into account beyond the mean-field level. This is possible within the framework of Hubbard-like models in the restricted subspace of orbitals close to the Fermi level.

The description of exciton physics in CTI is possible within the framework of one-band models, Penc98 ; Tsut99 ; Essler01 ; Eder02  but more detailed and realistic information may only be obtained from the multiband Hubbard (Emery) model that explicitly includes ligand ion degrees of freedom. For quasi-one-dimensional cuprates the Emery model was considered in Refs. Richt00, ; Hubsch01, ; Atzkern01, ; Moskv03, . The phenomenon of spin-charge separation that is characteristic to one dimension introduces a specific physics into the exciton formation. The mobility of a single electron or hole in 1D is not suppressed by spin correlations and the exciton dispersion is comparable with the one-particle dispersion.

In quasi-two-dimensional cuprates the situation is different. Due to AFM correlations the bandwidth of the one-particle motion is of the order of the AFM exchange integral $`J`$ which is considerably smaller than the bare hopping $`t`$, and the exciton dispersion is of the order of $`t`$. The authors of Ref. ZhNg,  proposed a qualitative physical explanation of the large exciton dispersion in layered cuprates: the propagation of an electron-hole pair does not disturb the AFM background, in contrast to the motion of a single electron or hole. Unfortunately, the calculations of Ref. ZhNg,  do not support this idea since they give no dispersion in the absence of oxygen on-site Coulomb repulsion $`U_p`$ and direct O-O hopping $`t_{pp}`$, which was also pointed out in Ref. Richt00, . One should expect a qualitative description of the EELS spectra Wang96 ; Fink01  already for the ’minimal’ version of the Emery model ($`1/U_d=U_p=t_{pp}=0`$) which can be refined by taking into account additional parameters like e.g. $`t_{pp}`$ or $`U_p`$. On the contrary, the authors of Ref. ZhNg,  had to assume quite unrealistic parameter values to fit the experimental spectra.

In this paper we outline an approach to calculate the longitudinal and the transverse dielectric response (Section II) for CTI within the framework of the multiband Hubbard model (Section III). Introducing an analogy to Wannier’s excitonic representation Wannier  we obtain a set of coupled equations of motion for the two-particle Green’s function that may be solved by means of the projection technique Fulde  (Section IV). It is substantial that the method allows a systematic improvement of approximations. In order to retain only the essential properties of CTI we first consider the ’minimal’ version of the Emery model describing the CuO2 plane of high-Tc superconductors and their parent compounds. The model reflects the main features of CTI: the existence of two kinds of states - strongly correlated ’copper’ states with the prohibition of double occupancy and uncorrelated ’oxygen’ states; the undoped plane has one hole per unit cell, and is a quantum two-dimensional AFM insulator possessing no long-range order at any finite temperature; the charge excitation corresponds to the transfer of a hole from copper to the adjacent oxygen site. In Section V, we discuss the solution of our equations within the small exciton basis and its dependence on various parameters of the extended Emery model. We obtain an appreciable exciton dispersion and a rough, qualitative agreement with the experimental spectra already within the ’minimal’ version of the Emery model.

## II Density-density correlation function

The dynamic structure factor is related to the density-density correlation function

$$S(𝐪,\omega )\frac{1}{2\pi N}_{\mathrm{}}^{\mathrm{}}𝑑t\mathrm{e}^{ı\omega t}\widehat{n}_𝐪(t)\widehat{n}_𝐪(0)=\frac{1}{\pi }\frac{1}{\mathrm{exp}(\beta \omega )1}\mathrm{Im}N_𝐪(\omega ),$$
(1)
where

$$\widehat{n}_𝐪=\frac{1}{\sqrt{N}}\underset{𝐫,s}{}\mathrm{exp}(ı\mathrm{𝐪𝐫})\left(a_{𝐫,s}^{}a_{𝐫,s}a_{𝐫,s}^{}a_{𝐫,s}\right)$$
(2)
is the electronic density operator in the localized basis, the summation runs over all lattice sites $`𝐫`$ and orbital sorts $`s`$ ; $`\mathrm{}`$ means the thermodynamic average. For $`\beta \omega 1`$ we have

$$S(𝐪,\omega )\frac{1}{\pi }\mathrm{Im}N(𝐪,\omega ),$$

where

$$N(𝐪,\omega )\widehat{n}_𝐪|\widehat{n}_𝐪=i_0^{\mathrm{}}𝑑t\mathrm{e}^{ı\omega t}[\widehat{n}_𝐪(t),\widehat{n}_𝐪(0)]$$
(3)
is the retarded Green’s function that defines the inverse dielectric function

$$\epsilon ^1(𝐪,\omega )=1+\frac{4\pi e^2}{v_cq^2}N(𝐪,\omega ),$$
(4)
with $`v_c`$ being the volume of the unit cell, and $`e`$ is the electronic charge. We neglect here the local field effects (coupling with the Fourier components with $`𝐪^{}=𝐪+𝐆`$, $`𝐆`$ being a reciprocal lattice vector). The function $`N(𝐪,\omega )`$ describes the response to the unscreened external potential. It requests the account of the macroscopic electric field, i.e. the long-range part of the Coulomb interaction. The latter is responsible for the splitting into longitudinal and transverse excitons for small wave number which is analogous to the splitting into longitudinal and transverse optical phonons. Knox  The response to the total, screened potential is given byPinesNoz

$$N_s(𝐪,\omega )=\epsilon (𝐪,\omega )N(𝐪,\omega ),$$
(5)
then

$$\epsilon (𝐪,\omega )=1\frac{4\pi e^2}{v_cq^2}N_s(𝐪,\omega ).$$
(6)
In the diagrammatic language the linear response to the total field may be expressed by the polarization operator where only irreducible graphs (which do not contain the contribution of the macroscopic electric field) should be taken into account. AmKo ; AGD  The random-phase approximation results in

$$N(𝐪,\omega )=\frac{N_s(𝐪,\omega )}{1\frac{4\pi e^2}{v_cq^2}N_s(𝐪,\omega )}$$
(7)
which follows from (5) and (4) and which is exact for $`q0`$, as it was shown in Ref. AmKo, .

The construction of Hubbard-like models for strongly correlated systems has an input from LDA band structure calculations where the screening of the long-range part of the Coulomb interaction is already taken into account. The Hubbard terms arise from the short-range residual interaction. Thus, the density response function $`N_H(𝐪,\omega )`$ calculated within the Hubbard model is an approximation to $`N_s(𝐪,\omega )`$. Penc98  In other words, it describes the motion of transverse (or ’mechanical’ by Agranovich’s Agran  terminology) excitons.

Using the spectral representation we may write

$$N_s(𝐪,z)=_0^{\mathrm{}}\left[\frac{1}{\pi }\mathrm{Im}N_s(𝐪,\omega ^{})\right]\frac{2\omega ^{}d\omega ^{}}{z^2\omega ^2}=_0^{\omega _0}+_{\omega _0}^{\mathrm{}}=N_H(𝐪,z)+N_{\mathrm{}}(𝐪,z).$$
(8)
Here we bear in mind that the Hubbard model contributes to transitions in the low frequency region $`\omega <\omega _0`$ with $`\omega _0`$ of the order of the bandwidth, and the electrons of the rest of the solid are excited only at higher energies. In zero approximation we may assume that in the frequency region $`\omega >\omega _0`$ the electronic polarization of the rest of the solid follows immediately the external field

$$N_{\mathrm{}}(𝐪,z)N_{\mathrm{}}(𝐪,0).$$

In other words, the Hubbard model is embedded into the medium with dielectric permeability

$$\epsilon _{\mathrm{}}(𝐪)=1\frac{4\pi e^2}{v_cq^2}N_{\mathrm{}}(𝐪,0).$$

In fact, $`\epsilon _{\mathrm{}}`$ may has its own dispersion and may be quite anisotropic for layered or quasi-one-dimensional compound. In principle, it should be taken from e.g. LDA calculations (we have assumed that the rest of solid is uncorrelated) or from experiment. It is obvious that the peak positions of the loss function

$$L(𝐪,\omega )\mathrm{Im}\left[\epsilon ^1(𝐪,\omega )\right]$$
(9)
and their intensity strongly depend on the value of $`\epsilon _{\mathrm{}}(𝐪)`$. Usually one neglects the $`𝐪`$-dependence and the anisotropy of $`\epsilon _{\mathrm{}}`$ , but it is a crude approximation, as well as another one which assumes $`\epsilon (𝐪,0)=\mathrm{const}`$ (see the discussion in Sec. V). For a quantitative description of EELS experiments the detailed knowledge of $`\epsilon _{\mathrm{}}(𝐪)`$ is necessary. Then the total dielectric function is

$$\epsilon (𝐪,\omega )=\epsilon _{\mathrm{}}\frac{4\pi e^2}{v_cq^2}N_H(𝐪,\omega )\epsilon _{\mathrm{}}\epsilon _H$$
(10)
and its inverse

$$\epsilon ^1(𝐪,\omega )=\epsilon _{\mathrm{}}^1\epsilon _H^1.$$

For the dielectric function of the Hubbard model $`\epsilon _H`$ the usual sum rule holds

$$_0^{\mathrm{}}\omega \mathrm{Im}\epsilon _H(𝐪,\omega )𝑑\omega =_0^{\mathrm{}}\omega \mathrm{Im}\epsilon _H^1(𝐪,\omega )𝑑\omega =\frac{\pi }{2}\frac{4\pi e^2}{\epsilon _{\mathrm{}}v_cq^2}[[\widehat{n}_𝐪,\widehat{H}_H],\widehat{n}_𝐪],$$

where $`\widehat{H}_H`$ is the Hubbard model Hamiltonian and the operator $`\widehat{n}_𝐪`$ acts in the subspace of orbitals which enter into $`\widehat{H}_H`$ (i.e. the summation over $`s`$ in Eq. (2) is restricted to these orbitals). Then for the total dielectric function approximated by Eq. (10) we have

$$_0^{\mathrm{}}\omega \mathrm{Im}\epsilon ^1(𝐪,\omega )𝑑\omega =\frac{\pi }{2}\frac{4\pi e^2}{\epsilon _{\mathrm{}}^2v_cq^2}[[\widehat{n}_𝐪,\widehat{H}_H],\widehat{n}_𝐪]=\frac{1}{\epsilon _{\mathrm{}}^2}_0^{\mathrm{}}\omega \mathrm{Im}\epsilon (𝐪,\omega )𝑑\omega .$$
(11)
The factor $`1/\epsilon _{\mathrm{}}^2`$ arises due to the neglect of the frequency dependence of $`\epsilon _{\mathrm{}}`$.

## III Model Hamiltonian and density operator

As we have mentioned in the Introduction, we consider the ’minimal’ Emery model that exhibits the essential properties of layered cuprates ($`1/U_d=U_p=t_{pp}=0`$). Then the total Hamiltonian in hole notation reads

$$\widehat{H}_H=\widehat{H}_0+\widehat{V},$$
(12)
where

$$\widehat{H}_0=\mathrm{\Delta }\underset{𝐫,\gamma }{}\overline{p}_{𝐫,\gamma }^{}\overline{p}_{𝐫,\gamma },\widehat{V}=t\underset{𝐑,\alpha ,\gamma }{}\left(\overline{p}_{𝐑+𝐚_\alpha ,\gamma }^{}\overline{Z}_𝐑^{0\gamma }+\overline{Z}_𝐑^{\gamma 0}\overline{p}_{𝐑+𝐚_\alpha ,\gamma }\right),$$

and where the Fermi operator $`\overline{p}_{𝐫,\gamma }`$ annihilates a hole at site $`𝐫`$ of the oxygen sublattice with spin projection index $`\gamma `$, the Hubbard projection operator $`\overline{Z}_𝐑^{0\gamma }=\overline{d}_{𝐑\gamma }(1n_{𝐑\overline{\gamma }})`$ annihilates a hole with spin index $`\gamma `$ on a singly occupied copper site, where $`\overline{d}_{𝐑\gamma }`$ is the corresponding Fermi operator. The double occupancy of copper sites is thus excluded from (12). $`\widehat{H}_0`$ includes the on-site energies ($`\mathrm{\Delta }=ϵ_pϵ_d,`$ $`ϵ_d`$ is taken as zero of energy), $`\widehat{V}`$ is the $`p`$-$`d`$ hybridization, $`\alpha =x,x,y,y`$ characterizes the direction of a nearest-neighbor vector $`𝐚`$; the phase factors in $`\widehat{V}`$ are absorbed into the definition of the operators $`\overline{p}_{𝐫,\sigma },\overline{Z}_𝐑^{0\gamma }`$ , they do not change the exciton dispersion.

Taking the limit $`U_d/\mathrm{\Delta }\mathrm{}`$ considerably simplifies the consideration and is a good approximation for weakly doped compounds in a wide range of values $`U_d/\mathrm{\Delta }>2`$ (see e.g. Ref. Stef97, ). The conditions $`U_p=t_{pp}=0`$ are introduced for simplicity and may be easily relaxed (see Appendix), then the Hubbard on-site term for $`p`$-orbitals and the direct O-O hopping of the form $`\widehat{t}_{pp}=t_{pp}_{i,j,\gamma }\overline{p}_{𝐫_i,\gamma }^{}\overline{p}_{𝐫_j,\gamma }`$ is added to $`\widehat{H}_0`$ .

It is well known that the Hamiltonian (12) has an insulating ground state and is equivalent to the nearest-neighbor AFM Heisenberg model in the low-energy region. It means that charge fluctuations in $`\widehat{V}`$ are strongly suppressed and that holes are localized. This fact becomes more apparent if we make a canonical transformation of operators of the form

$$\widehat{A}_{eff}=\mathrm{exp}(\widehat{S})\widehat{A}\mathrm{exp}(\widehat{S})=\widehat{A}+[\widehat{A},\widehat{S}]+\mathrm{},$$

where

$$\widehat{S}=\frac{t}{\mathrm{\Delta }}\underset{𝐑,\alpha =\pm x,\pm y,\gamma }{}\left(p_{𝐑+𝐚_\alpha ,\gamma }^{}Z_𝐑^{0\gamma }Z_𝐑^{\gamma 0}p_{𝐑+𝐚_\alpha ,\gamma }\right).$$

Then $`\widehat{H}_H`$ becomes

$`\widehat{H}_{eff}`$ $``$ $`\widehat{H}_04\tau {\displaystyle \underset{𝐑,\gamma }{}}Z_𝐑^{\gamma \gamma }+\tau {\displaystyle \underset{𝐑,\alpha 1,\alpha 2,\gamma }{}}p_{𝐑+𝐚_{\alpha 1},\gamma _1}^{}p_{𝐑+𝐚_{\alpha 2},\gamma _2}\left(Z_𝐑^{00}\delta _{\gamma _1\gamma _2}+Z_𝐑^{\gamma _2,\gamma _1}\right)`$ (13)
$`\tau {\displaystyle \underset{𝐑,\alpha ,\gamma _1}{}}Z_{𝐑+𝐠_\alpha }^{\gamma _10}Z_𝐑^{0\gamma _1}+\widehat{J}_s,`$
(see also Ref. Kuz98,  for the notation). Here $`p`$ and $`Z`$ mean transformed operators, $`\widehat{J}_s`$ is the AFM copper-copper superexchange interaction, $`𝐠`$ points to neighboring copper sites. Strictly speaking, the Hamiltonian (13) is obtained under the condition $`t/\mathrm{\Delta }1`$, and its parameters are $`\tau =t^2/\mathrm{\Delta }`$ and the AFM exchange $`Jt^4/\mathrm{\Delta }^3`$. Nevertheless, it may be applied in a wider range $`t/\mathrm{\Delta }<1`$ with renormalized values of $`\tau `$ and $`J`$.

The advantage of using the effective Hamiltonian (13) instead of the bare one (12) consists in excluding irrelevant zero-point charge fluctuations. Then the coupling of carriers with spin fluctuations that governs the low-energy physics of CTI becomes apparent. It is essential that the effective Hamiltonian (13) does not contain transitions between $`p`$\- and $`Z`$-states, in other words, it never creates a particle-hole state out of the dielectric state. In this sense it resembles the starting Hamiltonian for the transverse exciton motion in conventional insulators. This allows to introduce an analog of Wannier’s excitonic representation for the description of the electron-hole pair dynamics.

The bare density operator

$$\widehat{n}_𝐪=\frac{1}{\sqrt{N}}\left[\underset{𝐑,\sigma }{}\mathrm{exp}(ı\mathrm{𝐪𝐑})\left(\overline{Z}_𝐑^{\sigma \sigma }\overline{Z}_𝐑^{\sigma \sigma }\right)+\underset{𝐬}{}\mathrm{exp}(ı\mathrm{𝐪𝐬})\left(\overline{p}_{𝐬,\sigma }^{}\overline{p}_{𝐬,\sigma }\overline{p}_{𝐬,\sigma }^{}\overline{p}_{𝐬,\sigma }\right)\right],$$
(14)
($`𝐬`$ runs over O sublattice) is transformed to

$$\widehat{n}_𝐪=\frac{1}{\sqrt{N}}\{\underset{𝐑,\sigma }{}\mathrm{exp}(ı\mathrm{𝐪𝐑})[Z_𝐑^{\sigma \sigma }+\frac{t}{\mathrm{\Delta }}\underset{\alpha =\pm x,\pm y}{}(p_{𝐑+𝐚_\alpha ,\sigma }^{}Z_𝐑^{0\sigma }+Z_𝐑^{\sigma 0}p_{𝐑+𝐚_\alpha ,\sigma })]+$$
(15)
$$\underset{𝐑,\alpha =+x,+y}{}\mathrm{exp}[ı𝐪(𝐑+𝐚)][p_{𝐑+𝐚_\alpha ,\sigma }^{}p_{𝐑+𝐚_\alpha ,\sigma }\frac{t}{\mathrm{\Delta }}[(Z_𝐑^{\sigma 0}+Z_{𝐑+\mathrm{𝟐}𝐚_\alpha }^{\sigma 0})p_{𝐑+𝐚_\alpha ,\sigma }+h.c.]]\}.$$

Note that in the second line of (15) (in the sum over the oxygen sublattice) $`𝐚_\alpha `$ lies only in the same cell as $`𝐑`$. Collecting the terms surrounding a Cu site, we have

$$\widehat{n}_𝐪=\stackrel{~}{n}_𝐪+\frac{1}{\sqrt{N}}\frac{t}{\mathrm{\Delta }}\underset{𝐑,\sigma }{}\mathrm{exp}(ı\mathrm{𝐪𝐑})\underset{\alpha =\pm x,\pm y}{}\left(p_{𝐑+𝐚_\alpha ,\sigma }^{}Z_𝐑^{0\sigma }+Z_𝐑^{\sigma 0}p_{𝐑+𝐚_\alpha ,\sigma }\right)\left(1\mathrm{exp}(ı\mathrm{𝐪𝐚}_\alpha )\right)$$

$$=\stackrel{~}{n}_𝐪+\frac{1}{\sqrt{N}}\frac{t}{\mathrm{\Delta }}\underset{𝐑}{}\mathrm{exp}(ı\mathrm{𝐪𝐑})\underset{\alpha =\pm x,\pm y}{}\left(\psi _{𝐑,𝐚_\alpha }^{}+\psi _{𝐑,𝐚_\alpha }\right)\left(1\mathrm{exp}(ı\mathrm{𝐪𝐚}_\alpha )\right),$$
(16)
where the operator

$$\psi _{𝐑,\alpha }\underset{\gamma }{}Z_𝐑^{\gamma ,0}p_{𝐑+𝐚_\alpha ,\gamma }$$
(17)
annihilates an electron-hole pair with minimal distance, and

$$\stackrel{~}{n}_𝐪\frac{1}{\sqrt{N}}\left[\underset{𝐑,\sigma }{}\mathrm{exp}(ı\mathrm{𝐪𝐑})\left(Z_𝐑^{\sigma \sigma }\overline{Z}_𝐑^{\sigma \sigma }\right)+\underset{𝐬}{}\mathrm{exp}(ı\mathrm{𝐪𝐬})\left(p_{𝐬,\sigma }^{}p_{𝐬,\sigma }\overline{p}_{𝐬,\sigma }^{}\overline{p}_{𝐬,\sigma }\right)\right].$$
(18)
As we have mentioned above, the effective Hamiltonian (13) conserves the number of particles in every band. Therefore $`\stackrel{~}{n}_𝐪=0`$ gives no contribution to $`N_H`$. Having the operator (16) we may proceed with the calculation of the density-density response function (3).

## IV Electron-hole pair dynamics

The problem of the dielectric function (10) calculation is thus reduced to the calculation of the two-particle Green’s function

$$N_H(𝐪,\omega )=\mathrm{\Phi }_𝐪+\mathrm{\Phi }_𝐪^{}|\mathrm{\Phi }_𝐪^{}+\mathrm{\Phi }_𝐪_\omega =\mathrm{\Phi }_𝐪|\mathrm{\Phi }_𝐪^{}_\omega +\mathrm{\Phi }_𝐪|\mathrm{\Phi }_𝐪^{}_\omega $$
(19)
where

$$\mathrm{\Phi }_𝐪\frac{t}{\mathrm{\Delta }}\underset{\alpha }{}\left(1\mathrm{exp}(ı\mathrm{𝐪𝐚}_\alpha )\right)\psi _{𝐪,\alpha },\psi _{𝐪,\alpha }=\frac{1}{\sqrt{N}}\underset{𝐑}{}\mathrm{exp}(ı\mathrm{𝐪𝐑})\psi _{𝐑\alpha },$$
(20)
and where we used again the conservation of particle numbers in the electron and the hole subsystem that excludes anomalous Green’s functions like $`\mathrm{\Phi }_𝐪|\mathrm{\Phi }_𝐪_\omega `$.

The equation of motion

$$\omega \mathrm{\Phi }_𝐪|\mathrm{\Phi }_𝐪^{}_\omega =[\mathrm{\Phi }_𝐪,\mathrm{\Phi }_𝐪^{}]+[\mathrm{\Phi }_𝐪,\widehat{H}_{eff}]|\mathrm{\Phi }_𝐪^{}_\omega $$
(21)
generates more complex operators

$$\xi _{𝐑,\alpha ,\beta }=\underset{\gamma _1,\gamma _2}{}Z_{𝐑+𝐠_\alpha }^{\gamma _1\gamma _2}Z_𝐑^{\gamma _2,0}p_{𝐑+𝐠_\alpha +𝐚_\beta ,\gamma _1},$$
(22)
which annihilate states with an increasing separation between electron and hole, accompanied by spin fluctuations. The set of equations of motion will generate states corresponding to electrons and holes that are more and more separated and dressed by spin fluctuations. These states form a set similar to the excitonic representation for conventional insulators. Wannier ; Knox ; Egri  The complication that arises in CTI consists in the strong interaction of both electron and hole with AFM fluctuations. The effect of this interaction leads to a strong renormalization of the one-particle bandwidth, but it is partially canceled when electron and hole follow each other.

The set of coupled equations (21) may be approximately solved by means of the projection technique. We choose an operator basis $`B_{𝐪,i}`$ and the definition of the scalar product $`[B_{𝐪,i},B_{𝐪,j}^{}]`$. Within the operator subspace spanned by this basis we are looking for the approximate solution of the eigenvalue problem

$$[\mathrm{\Psi }_𝐪,\widehat{H}]=E_𝐪\mathrm{\Psi }_𝐪,\mathrm{\Psi }_𝐪=\underset{i}{}c_i(𝐪)B_{𝐪,i}.$$
(23)
This leads, as usual for a non orthonormal basis, to the generalized eigenvalue problem

$$\underset{i}{}c_i(𝐪)L_{i,j}(𝐪)=E\underset{i}{}c_i(𝐪)S_{ij}(𝐪)$$
(24)
where overlap and Liouvillean matrices

$$S_{ij}[B_i,B_j^{}]=B_i,B_j^{},L_{ij}[[B_i,\widehat{H}],B_j^{}]=[B_i,\widehat{H}],B_j^{}$$
(25)
depend only on spin-spin correlation functions for the system without electron-hole pairs, which is equivalent to the Heisenberg antiferromagnet. The correlation functions may thus be calculated from the Heisenberg model.

The system (24) may be solved numerically. Then we find the eigenvectors

$$\mathrm{\Psi }_𝐪^\lambda =\underset{i}{}c_i^\lambda (𝐪)B_{𝐪,i},$$
(26)
where $`\lambda `$ is the number of the branch in the spectrum. In order to calculate $`N_H(𝐪,\omega )`$ within our basis, we should expand

$$\mathrm{\Phi }_𝐪=\underset{\lambda }{}g^\lambda (𝐪)\mathrm{\Psi }_𝐪^\lambda .$$
(27)
Then one obtains

$$\mathrm{\Phi }_𝐪|\mathrm{\Phi }_𝐪^{}_\omega =\underset{\lambda }{}\left|g^\lambda (𝐪)\right|^2\mathrm{\Psi }_𝐪^\lambda |\left(\mathrm{\Psi }_𝐪^\lambda \right)^{}_\omega =\underset{\lambda }{}\left|g^\lambda (𝐪)\right|^2\left(\omega E_𝐪^\lambda \right)^1,$$
(28)
with

$$\left|g^\lambda (𝐪)\right|^2=\left(\frac{t}{\mathrm{\Delta }}\right)^2\left|\underset{i,\alpha }{}c_i^\lambda (𝐪)S_{i\alpha }\left[1\mathrm{exp}(ı\mathrm{𝐪𝐚}_\alpha )\right]\right|^2,$$
(29)
and finally

$$\epsilon (𝐪,\omega )=\epsilon _{\mathrm{}}\frac{4\pi ^2e^2}{v_cq^2}\underset{\lambda }{}\left|g^\lambda (𝐪)\right|^2\left[\left(\omega E_𝐪^\lambda \right)^1\left(\omega +E_𝐪^\lambda \right)^1\right].$$
(30)
Let us note that the projection technique allows to improve the chosen approximation step by step by enlarging the basis set.

## V Results and discussion

We have restricted ourself to the minimal basis that describes the electron-hole pair with minimal distance. The basis contains operators (17) and (22) with $`\beta =\alpha `$. Then the problem (23) has the dimension $`8\times 8`$. The overlap and Liouvillean matrices are given in the Appendix. Spin-spin correlation functions were taken from the spherically symmetric treatment of an S=1/2 Heisenberg AFM model on the square lattice. BarBer94  For a low temperature $`T=0.1J`$ and a vanishing frustration parameter $`p=0.01`$ they have the following values $`\widehat{𝐒}_𝐑\widehat{𝐒}_{𝐑+𝐠}=0.33`$, $`\widehat{𝐒}_𝐑\widehat{𝐒}_{𝐑+𝐠_x+𝐠_y}=0.20,\widehat{𝐒}_𝐑\widehat{𝐒}_{𝐑+\mathrm{𝟐}𝐠}=0.17`$. For the on-site energy difference and $`p`$-$`d`$ hopping we took the values $`\mathrm{\Delta }=3.6`$eV and $`t=1.3`$ eV, which are characteristic to all cuprates.

Fig. 1 shows the dispersion of the imaginary part of the dielectric function $`\epsilon _2(𝐪,\omega )\mathrm{Im}\epsilon (𝐪,\omega )`$. The oxygen on-site repulsion and the O-O hopping were neglected. We see strongly dispersive branches both in the  and in the  direction.

As we have mentioned above, the comparison with the EELS experiment may have only qualitative character without a detailed knowledge of the background dielectric constant $`\epsilon _{\mathrm{}}(𝐪)`$. Fig. 2 shows the graphs for the loss function (multiplied by $`\left[\epsilon _{\mathrm{}}(\mathrm{𝟎})\right]^2`$ in order to have approximately the same normalization as $`\epsilon _2(𝐪,\omega )`$ according to (11)) under the assumption that ZhNg

$$\epsilon (𝐪,\omega =0)\epsilon (𝐪=\mathrm{𝟎},\omega =0)=4.83.$$
(31)
Then the value of $`\epsilon _{\mathrm{}}`$ was obtained from Eq. (30). In this case the dispersion in the loss function reproduces essentially the dispersion in $`\epsilon (𝐪,\omega )`$. The peaks are slightly shifted to higher energies. In fact, the assumption (31) implies that the dispersion of $`N_{\mathrm{}}(𝐪,\omega )`$ should follow the dispersion of $`N_H(𝐪,\omega )`$ in such a way that

$$\epsilon (𝐪,\omega =0)=1\frac{4\pi e^2}{v_cq^2}\left[N_H(𝐪,0)+N_{\mathrm{}}(𝐪,0)\right]=\mathrm{const}$$

as it follows from (8) and (6).

In general, the interplay of $`N_H(𝐪,\omega )`$ and $`N_{\mathrm{}}(𝐪,\omega )`$ should be more complex. In order to demonstrate the strong dependence on the value of $`\epsilon _{\mathrm{}}`$, we plot in Fig. 3 the same graph assuming a constant value $`\epsilon _{\mathrm{}}=2`$ for all $`𝐪`$. We see a qualitative difference with Fig. 2 and we may conclude that the dependence of the loss function on $`\epsilon _{\mathrm{}}`$ is nontrivial. For quasi-one dimensional compounds a large values of $`\epsilon _{\mathrm{}}8`$ was taken in Refs. Hubsch01, ; Atzkern01, . This means that the rest of the solid strongly screens the long-range part of the Coulomb interaction between electrons that enter the Hubbard model. For this situation the poles of $`\epsilon (𝐪,\omega )`$ are very close to the poles of $`\epsilon ^1(𝐪,\omega )`$ and the shape of the loss function is close to the shape of $`\epsilon _2(𝐪,\omega )`$. Note also that always $`\epsilon (𝐪,0)>\epsilon _{\mathrm{}}`$ as follows from Eq. (10), (30) and with $`\epsilon _{\mathrm{}}=8`$ one will receive unrealistically large $`\epsilon (𝐪,0)`$.

Let us now show some examples for the dependence of the dielectric function on various parameters of the model. For the reasons outlined above, from now on we consider only figures for $`\epsilon _2(𝐪,\omega )`$. As it was mentioned in Section III, the model (12) may be generalized in order to include a finite Hubbard repulsion $`U_p`$ on the oxygen site and a direct oxygen-oxygen hopping $`t_{pp}`$.

The dependence on $`t_{pp}`$ is not very strong for $`t_{pp}<\tau `$. Let us recall that the addition of the O-O hopping to the Hamiltonian (13) within the parameter range $`t_{pp}=(0.3÷0.4)\tau `$ is essential to describe correctly the angle-resolved photoemission spectroscopy (ARPES) of layered cuprates. Kuz98 ; Bar01  Fig. 4 shows $`\epsilon _2(𝐪,\omega )`$ for $`t_{pp}=0.4\tau `$ , and $`U_p=0`$. We see that the main difference to Fig. 1 consists in the redistribution of spectral weight between the different branches of the spectrum.

The spectrum demonstrates a much stronger dependence on $`U_p`$. Fig. 5 displays $`\epsilon _2(𝐪,\omega )`$ for $`U_p=4`$ eV. We have an almost vanishing dispersion of the lower branches. Let us pay attention to the fact that $`U_p`$ does not affect the single hole motion. Due to that reason its value is experimentally not well established. Our results show that the dielectric function dispersion is very sensitive to this parameter.

Now let us discuss a very important peculiarity of our figures, namely the absolute position of the intensive peaks in the spectrum of $`\epsilon _2(𝐪,\omega )`$. At the $`\mathrm{\Gamma }`$ point ($`𝐪=(0,0)`$) we have one peak with the energy $`E_𝚪=\mathrm{\Delta }+4\tau 5.5`$ eV. Let us estimate the edge of the electron-hole continuum. It corresponds to the energy that is needed for the excitation of electron and hole which are independent on each other. The operator that annihilates such a state is

$$\mathrm{\Psi }_𝐪^{eh}=e_{𝐤+𝐪}h_𝐤,[\mathrm{\Psi }_𝐪^{eh},\widehat{H}_{eff}]=(ϵ_{𝐤+𝐪}^e+ϵ_𝐤^h)\mathrm{\Psi }_𝐪^{eh},$$

where $`ϵ_𝐤^e=(4\tau ϵ_𝐤^{tJ})`$ corresponds to the energy of a single ’electron’ quasiparticle, that will be a complex spin-polaron corresponding to the coherent motion in the so called $`tJ`$ model, which describes the $`Z`$ subsystem in (13) in the absence of holes (with $`t^{tJ}=\tau `$). $`ϵ_𝐤^h=\mathrm{\Delta }ϵ_𝐤^s`$ is the energy of the coherent motion of the Zhang-Rice singlet, dressed by spin fluctuations. The minimal energy that one needs to excite such a pair is

$$E_{\mathrm{min}}^{eh}\mathrm{\Delta }ϵ_{\mathrm{min}}^s\left(4\tau ϵ_{\mathrm{min}}^{tJ}\right)\mathrm{\Delta }2.2\tau 2.6\mathrm{eV}.$$

Here we have taken into account that the minimum of the spectrum in the $`tJ`$ model is $`ϵ_{\mathrm{min}}^{tJ}2\tau `$ and that the Zhang-Rice singlet energy is $`ϵ_{\mathrm{min}}^s4.2\tau `$. Kuz98  This result indicates that for layered cuprates the excitonic feature is immersed into the electron-hole continuum and represents a resonance rather than a discrete level. Of course, for the the final conclusion more detailed calculations should be performed within an enlarged basis set.

Summarizing the calculated spectra we may state a rough, qualitative, agreement with the experimental curves Wang96 ; Fink01  already within the ’minimal’ version of the Emery model and using a minimal basis set (Figs. 1-3). One may note a remarkable influence of the background dielectric function $`\epsilon _{\mathrm{}}`$ (Figs. 2 and 3). We found that additional parameters like e.g. $`t_{pp}`$ or $`U_p`$ act in different ways. That might improve a future parameter fit of the experimental curves. But the present accuracy is not sufficient for a reliable fit, which has to be reserved for the future.

In this work we have not taken into account the intersite Coulomb repulsion which leads to electron hole attraction. This term has a twofold influence on the position of the excitonic feature. From one hand, it leads to an effective increase of the fundamental gap, from the other hand it contributes to the electron-hole binding. These tendencies are opposite to each other and should be thoroughly explored in a separate investigation.

## Acknowledgments

This work was supported by DFG (436 UKR 113/49/41) and by RFFR (project N 01-12-16719). R.O.K. and A.F.B. thank for hospitality the IFW Dresden, where the main part of this work has been carried out. We wish to thank S.-L. Drechsler and J. Málek for numerous and useful discussions.

*

## Appendix

Here we give the formulae for the matrices (25) for the Hamiltonian ( 12) that includes also a finite oxygen on-site repulsion $`U_p`$ and a direct O-O hopping $`t_{pp}`$. The effective Hamiltonian may be generally derived by the following procedure. Let be

$$\widehat{H}=\widehat{H}_0+\widehat{V},\widehat{H}_0=\underset{m}{}|mE_mm|,\widehat{V}=\underset{m}{}|nt_{nm}m|,$$
(32)
then the canonical transformation with the operator

$$\widehat{S}=\underset{m}{}|n\frac{t_{nm}}{E_mE_n}m|$$
(33)
gives $`[\widehat{H}_0,\widehat{S}]=\widehat{V}`$, and finally up to the second order

$$\widehat{H}_{eff}=\mathrm{exp}(\widehat{S})\widehat{H}\mathrm{exp}(\widehat{S})=H_0+\frac{1}{2}[\widehat{V},\widehat{S}]=$$

$$=\frac{1}{2}\underset{m}{}|n\left(\frac{t_{nj}t_{jm}}{E_mE_j}+\frac{t_{nj}t_{jm}}{E_nE_j}\right)m|.$$
(34)
Within the operator basis containing (17) and (22) we have the equations of motion for the Hamiltonian containing finite $`U_p,t_{pp}`$

$$[\psi _{𝐑,\alpha },\widehat{H}]=E_0\psi _{𝐑,\alpha }+\tau \underset{\beta }{}\psi _{𝐑,\beta }+\tau _u\xi _{𝐑,\alpha ,\alpha }+[\tau _u\xi _{𝐑+𝐠_\alpha ,\alpha ,\alpha }(\tau _u\tau )\psi _{𝐑+𝐠_\alpha ,\alpha }]t_{pp}\underset{\beta \alpha ,\alpha }{}\psi _{𝐑,\beta }$$

$$[\xi _{𝐑,\alpha ,\alpha },\widehat{H}]=E_0\xi _{𝐑,\alpha ,\alpha }+\tau \xi _{𝐑\alpha ,\alpha }+\tau \underset{\beta \alpha }{}Z_{𝐑+𝐠_\alpha }^{\gamma _1,\gamma _2}Z_𝐑^{\gamma _2,0}p_{𝐑+𝐚_\beta ,\gamma _1}t_{pp}\underset{\beta \alpha ,\alpha }{}Z_{𝐑+𝐠_\alpha }^{\gamma _1,\gamma _2}Z_𝐑^{\gamma _2,0}p_{𝐑+𝐚_\beta ,\gamma _1}+$$

$$+\tau _u\psi _{𝐑,\alpha }+\tau _u\psi _{𝐑+𝐠_\alpha ,\alpha }(\tau _u\tau )\xi _{𝐑+𝐠_\alpha ,\alpha ,\alpha }.$$

Here $`\tau _ut^2/(\mathrm{\Delta }+U_p)`$ and $`E_0=\mathrm{\Delta }+4\tau +2(\tau \tau _u)`$ . It is convenient to introduce the notation

$$\omega _𝐫\underset{\gamma ,\gamma _1,\gamma _2,\mathrm{}}{}Z_0^{\gamma ,\gamma _1}Z_𝐠^{\gamma _1,\gamma _2}\mathrm{}Z_𝐫^{\gamma _r,\gamma }.$$

For $`𝐫`$ up to the third neighbors $`\omega _𝐫`$ is expressed via two point correlation functions

$$\omega _g=\frac{1}{2}+\widehat{𝐒}_𝐑\widehat{𝐒}_{𝐑+𝐠},\omega _{𝐠_\alpha +𝐠_\beta }=\frac{1}{4}+2\widehat{𝐒}_𝐑\widehat{𝐒}_{𝐑+𝐠}+\widehat{𝐒}_𝐑\widehat{𝐒}_{𝐑+𝐠_\alpha +𝐠_\beta }.$$
(35)
The Liouvillean and overlap matrices in $`𝐤`$-space are then

$$\{\xi _{𝐤\alpha ,\alpha },\psi _{𝐤\beta }^{}\}=\delta _{\alpha \beta }\omega _g,$$
(36)
$$\{[\psi _{𝐤\alpha },\widehat{H}],\psi _{𝐤\beta }^{}\}=\delta _{\alpha \beta }\left(\tau +\tau _u\omega _g\right)+\left(\tau t_{pp}\right)\left(1\delta _{\alpha \beta }\right)\left(1\delta _{\alpha ,\beta }\right)$$
(37)
$$+\delta _{\alpha ,\beta }\left[\tau +\mathrm{exp}(ı\mathrm{𝐤𝐠}_\alpha )\left(\tau _u\omega _g\tau _u+\tau \right)\right],$$
(38)
$$\{[\xi _{𝐤\alpha ,\alpha },\widehat{H}],\psi _{𝐤\beta }^{}\}=\delta _{\alpha \beta }\left(\tau _u+\tau \omega _g\right)+\left(\tau t_{pp}\right)\omega _g\left(1\delta _{\alpha \beta }\right)\left(1\delta _{\alpha ,\beta }\right)$$
(39)
$$+\delta _{\alpha ,\beta }\left\{\tau \omega _g+\mathrm{exp}(ı\mathrm{𝐤𝐠}_\alpha )[\tau _u(\tau _u\tau )\omega _g]\right\},$$

$$\{[\xi _{𝐤\alpha ,\alpha },\widehat{H}],\xi _{𝐤\beta ,\beta }^{}\}=[\delta _{\alpha \beta }(\tau +\tau _u\omega _g)+\left(1\delta _{\alpha \beta }\right)\left(1\delta _{\alpha ,\beta }\right)\left(\tau t_{pp}\right)\omega _{𝐠_\beta +𝐠_\alpha }]$$
(40)
$$+\delta _{\alpha ,\beta }\left\{\tau \omega _{2g}+\mathrm{exp}(ı\mathrm{𝐤𝐠}_\alpha )[\tau _u\omega _g(\tau _u\tau )]\right\}.$$


format: text

text:
# INTRODUCTION TO 𝐵 PHYSICS

## 1 Introduction

The completion of a pioneering program on the physics of $`B`$ mesons at the Cornell Electron Synchrotron (CESR) and the highly successful commissioning of asymmetric $`e^+e^{}`$ colliders at PEP-II and KEK-B has led to a wealth of data on the decays of $`B`$ mesons which can shed light on the weak and strong interactions and on the violation of CP symmetry. Hadron colliders are also beginning to utilize specialized triggers to study $`B`$ mesons under conditions of higher background but with the benefit of larger production cross sections. One session of the present workshop is devoted to these results. The present Introduction seeks to put these exciting results in a broader context, showing both what has been learned so far and what the future holds.

After reviewing the charge-changing weak transitions of quarks, described by the Cabibbo-Kobayashi-Maskawa (CKM) matrix, in Sec. 2, we discuss several decays of neutral $`B`$ mesons to CP eigenstates in Sec. 3, and a sampling of $`B\pi K`$ decays in Sec. 4. Sec. 5 concludes.

## 2 The Cabibbo-Kobayashi-Maskawa Matrix

The weak transitions between quarks with $`Q=2/3`$ ($`u,c,t`$) and $`Q=1/3`$ ($`d,s,b`$) are encoded in a $`3\times 3`$ unitary matrix $`V`$ known as the Cabibbo-Kobayashi-Maskawa (CKM) matrix.$`^{\mathrm{?},\mathrm{?}}`$ The $`tb`$, $`cs`$, and $`ud`$ elements of $`V`$ are of approximately unit magnitude, while $`|V_{us}||V_{cd}|0.22`$, $`|V_{cb}||V_{ts}|0.04`$, $`|V_{td}|0.008`$, $`|V_{ub}|0.004`$. Phases in the last two matrix elements can account for the observed CP violation in the kaon system, as first noted by Kobayashi and Maskawa.$`^\mathrm{?}`$ One now seeks tests of this mechanism through $`B`$ decays. A convenient parametrization of the CKM matrix, mentioned elsewhere in these Proceedings, is due to Wolfenstein.$`^\mathrm{?}`$

### 2.1 Unitarity

The unitarity of the CKM matrix implies that the scalar product of any column with the complex conjugate of any other column is zero, for example, $`V_{ub}^{}V_{ud}+V_{cb}^{}V_{cd}+V_{tb}^{}V_{td}=0`$. If one divides by $`V_{cb}^{}V_{cd}`$, this relation becomes equivalent to a triangle in the complex $`\overline{\rho }+i\overline{\eta }`$ plane, with vertices at (0,0) (angle $`\varphi _3=\gamma `$), (1,0) (angle $`\varphi _1=\beta `$), and $`(\overline{\rho },\overline{\eta })`$ (angle $`\varphi _2=\alpha `$). The triangle has unit base and its other two sides are $`\overline{\rho }+i\overline{\eta }=(V_{ub}^{}V_{ud}/V_{cb}^{}V_{cd})`$ (opposite $`\varphi _1=\beta `$) and $`1\overline{\rho }i\overline{\eta }=(V_{tb}^{}V_{td}/V_{cb}^{}V_{cd})`$ (opposite $`\varphi _3=\gamma `$). The result is shown in Fig. 1.

### 2.2 Parameters

Direct constraints on the CKM parameters are obtained from strange particle decays ($`V_{us}0.22`$), $`bc`$ decays ($`V_{cb}0.041`$), and $`bu`$ decays ($`|V_{ub}/V_{cb}|0.080.10`$). Indirect constraints arise from various processes involving flavor-changing box diagrams. The $`\overline{s}d\overline{d}s`$ transition gives rise to CP-violating $`K^0`$$`\overline{K}^0`$ mixing, whose magnitude (expressed through the parameter $`ϵ`$) imposes a constraint on Im($`V_{td}^2`$). The $`\overline{b}d\overline{d}b`$ transition generates $`B^0`$$`\overline{B}^0`$ mixing, from which a constraint on $`|V_{td}||1\overline{\rho }i\overline{\eta }|`$ may be obtained. The $`\overline{b}s\overline{s}b`$ transition generates $`B_s`$$`\overline{B}_s`$ mixing; by comparing the lower limit on this mixing with the observed $`B^0`$$`\overline{B}^0`$ mixing and using estimates of flavor-SU(3)-symmetry breaking in matrix elements, one finds $`|V_{ts}/V_{td}|>4.4`$. The combined impact of these constraints can be conservatively summarized as requiring $`0.08\overline{\rho }0.34`$, $`0.25\overline{\eta }0.43`$,$`^\mathrm{?}`$ though more restrictive analyses (see, e.g., $`^\mathrm{?}`$) appear.

## 3 $`B^0`$ Decays to CP Eigenstates

Consider the decays $`B^0f`$ (amplitude $`A`$) and $`\overline{B}^0f`$ (amplitude $`\overline{A}`$), where $`f`$ is a CP eigenstate with eigenvalue $`\xi _f=\pm 1`$. \[See reviews $`^{\mathrm{?},\mathrm{?}}`$ for more details.\] As a result of $`B^0`$$`\overline{B}^0`$ mixing, a state which is $`B^0`$ at proper time $`t=0`$ will evolve into one, denoted $`B^0(t)`$, which is a mixture of $`B^0`$ and $`\overline{B}^0`$. Consequently, there will be two pathways to the final state $`f`$: one from $`B^0`$ through the amplitude $`A`$ and the other from $`\overline{B}^0`$ through the amplitude $`\overline{A}`$, which acquires an additional phase $`2\varphi _1`$ through $`B^0`$$`\overline{B}^0`$ mixing. The interference of these two amplitudes can be different in the decay $`B^0(t)f`$ from that in $`\overline{B}^0(t)f`$, leading to a time-integrated rate asymmetry

$$𝒜_{CP}\frac{\mathrm{\Gamma }(\overline{B}^0f)\mathrm{\Gamma }(B^0f)}{\mathrm{\Gamma }(\overline{B}^0f)+\mathrm{\Gamma }(B^0f)}$$
(1)
as well as to time-dependent rates

$$\left\{\begin{array}{c}\mathrm{\Gamma }[B^0(t)f]\\ \mathrm{\Gamma }[\overline{B}^0(t)f]\end{array}\right\}e^{\mathrm{\Gamma }t}[1𝒜_f\mathrm{cos}\mathrm{\Delta }mt𝒮_f\mathrm{sin}\mathrm{\Delta }mt],$$
(2)
where

$$𝒜_f\frac{|\lambda |^21}{|\lambda |^2+1},𝒮_f\frac{2\mathrm{I}\mathrm{m}\lambda }{|\lambda |^2+1},\lambda e^{2i\varphi _1}\frac{\overline{A}}{A}.$$
(3)
Note that one must have $`𝒮_f^2+𝒜_f^21`$. I now discuss specific cases.

### 3.1 $`B^0J/\psi K_S`$ and $`\varphi _1=\beta `$

For this decay one has $`\overline{A}/A\xi _{J/\psi K_S}=1`$. One finds that the time-integrated asymmetry $`𝒜_{CP}`$ is proportional to $`\mathrm{sin}(2\varphi _1)`$. Using this and related decays involving the same $`\overline{b}\overline{s}c\overline{c}`$ subprocess, BaBar $`^\mathrm{?}`$ finds $`\mathrm{sin}(2\varphi _1)=0.741\pm 0.067\pm 0.033`$ while Belle $`^\mathrm{?}`$ finds $`\mathrm{sin}(2\varphi _1)=0.719\pm 0.074\pm 0.035`$. The two values are quite consistent with one another; the world average $`^\mathrm{?}`$ is $`\mathrm{sin}(2\varphi _1)=0.734\pm 0.054`$, consistent with other determinations.$`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$

### 3.2 $`B^0\pi ^+\pi ^{}`$ and $`\varphi _2=\alpha `$

Here the situation is more complicated because there are two competing amplitudes contributing to the decay: a “tree” $`T`$ and a “penguin” $`P`$. The decay amplitudes are then

$$A=(|T|e^{i\varphi _3}+|P|e^{i\delta }),\overline{A}=(|T|e^{i\varphi _3}+|P|e^{i\delta }).$$
(4)
The parameter $`\delta `$ is the relative $`P/T`$ strong phase. The asymmetry $`𝒜_{CP}`$ would be proportional to $`\mathrm{sin}(2\varphi _2)`$ if the penguin amplitude could be neglected. However, it cannot, so methods have been developed to deal with its contribution.

An isospin analysis $`^\mathrm{?}`$ makes use of $`B`$ decays to $`\pi ^+\pi ^{}`$, $`\pi ^\pm \pi ^0`$, and $`\pi ^0\pi ^0`$ to separate the contributions of decays involving $`I=0`$ and $`I=2`$ final states. Information can then be obtained on both strong and weak phases. A potential problem with this method is that the branching ratio of $`B^0`$ to $`\pi ^0\pi ^0`$ may be very small, of order $`10^6`$. I shall discuss instead methods $`^{\mathrm{?},\mathrm{?}}`$ in which flavor symmetry is used to estimate the magnitude of the penguin amplitude.$`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$

The tree amplitude for $`B^0(=\overline{b}d)\pi ^+\pi ^{}`$ involves the subprocess $`\overline{b}\pi ^+\overline{u}`$, with the spectator $`d`$ quark combining with $`\overline{u}`$ to form a $`\pi ^{}`$. Its magnitude is $`|T|`$; its weak phase is Arg($`V_{ub}^{})=\varphi _3`$; by convention its strong phase is 0. The penguin amplitude involves the flavor structure $`\overline{b}\overline{d}`$, with the final $`\overline{d}d`$ pair fragmenting into $`\pi ^+\pi ^{}`$. Its magnitude is $`|P|`$. The dominant $`t`$ contribution in the loop diagram for $`\overline{b}\overline{d}`$ can be integrated out and the unitarity relation $`V_{td}V_{tb}^{}=V_{cd}V_{cb}^{}V_{ud}V_{ub}^{}`$ used. The $`V_{ud}V_{ub}^{}`$ contribution can be absorbed into a redefinition of the tree amplitude, after which the weak phase of the penguin amplitude is 0 (mod $`\pi `$). By definition, its strong phase is $`\delta `$.

Knowledge of the time-dependent asymmetries $`𝒮_{\pi \pi }`$ and $`𝒜_{\pi \pi }`$ suffices to specify both $`\varphi _3`$ (or $`\varphi _2=\pi \varphi _1\varphi _3`$) and $`\delta `$, if one has an independent estimate of $`|P/T|`$. This may be obtained by using flavor SU(3) to get $`|P|`$ from $`B^+K^0\pi ^+`$ $`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$ and factorization to get $`|T|`$ from $`B\pi l\nu `$.$`^\mathrm{?}`$ Since the $`T`$ amplitude obtained from factorization is not precisely the same as that which contains a small $`V_{ud}V_{ub}^{}`$ contribution from the penguin, an alternative method $`^{\mathrm{?},\mathrm{?}}`$ makes direct use of the measured ratio of the $`B^+K^0\pi ^+`$ and $`B^0\pi ^+\pi ^{}`$ branching ratios to constrain $`|P/T|`$. I shall discuss the first method since it is simpler.

In addition to $`𝒮_{\pi \pi }`$ and $`𝒜_{\pi \pi }`$, a useful quantity is the ratio of the $`B^0\pi ^+\pi ^{}`$ branching ratio $`\overline{}(\pi ^+\pi ^{})`$ (averaged over $`B^0`$ and $`\overline{B}^0`$) to that due to the tree amplitude alone:

$$R_{\pi \pi }\frac{\overline{}(\pi ^+\pi ^{})}{\overline{}(\pi ^+\pi ^{})|_{\mathrm{tree}}}=12\left|\frac{P}{T}\right|\mathrm{cos}\delta \mathrm{cos}(\varphi _1+\varphi _2)+\left|\frac{P}{T}\right|^2.$$
(5)
One also has

$$R_{\pi \pi }𝒮_{\pi \pi }=\mathrm{sin}2\varphi _2+2\left|\frac{P}{T}\right|\mathrm{cos}\delta \mathrm{sin}(\varphi _1\varphi _2)\left|\frac{P}{T}\right|^2\mathrm{sin}(2\varphi _1),$$
(6)
$$R_{\pi \pi }𝒜_{\pi \pi }=2|P/T|\mathrm{sin}\delta \mathrm{sin}(\varphi _1+\varphi _2).$$
(7)
The value of $`\varphi _1=\beta `$ is specified to within a few degrees; we shall take it to have its central value $`\varphi _1=23.6^{}`$. The value of $`|P/T|`$ (updating $`^{\mathrm{?},\mathrm{?}}`$) is $`0.28\pm 0.06`$. Taking the central value, we plot trajectories in the ($`𝒮_{\pi \pi },𝒜_{\pi \pi }`$) plane as $`\delta `$ is allowed to vary from $`\pi `$ to $`\pi `$. The result is shown in Fig. 2.

The experimental situation regarding the time-dependent asymmetries is not yet settled. As shown in Table 1, BaBar $`^\mathrm{?}`$ and Belle $`^\mathrm{?}`$ obtain very different values, especially for $`𝒮_{\pi \pi }`$. Even if this conflict were to be resolved, however, one sees the possibility of a discrete ambiguity, since curves for different values of $`\varphi _2`$ intersect one another.

The discrete ambiguity may be resolved with the help of $`R_{\pi \pi }`$. The most recent average of branching ratios, including ones from Belle presented at this Conference,$`^\mathrm{?}`$ yields $`\overline{}(B^0\pi ^+\pi ^{})=(4.55\pm 0.44)\times 10^6`$ (see below), while $`^\mathrm{?}`$ $`\overline{}(B^0\pi ^+\pi ^{})|_{\mathrm{tree}}=(7.3\pm 3.2)\times 10^6`$, so that $`R_{\pi \pi }=0.62\pm 0.28`$. One can then plot $`R_{\pi \pi }`$ as a function of $`𝒮_{\pi \pi }`$ for any value of $`\delta `$; some examples are shown in Fig. 3.

If the errors on $`R_{\pi \pi }`$ can be reduced to $`\pm 0.1`$, a distinction between $`\delta =0`$ and $`\delta =\pm \pi `$ will be possible. One possibility for improving this situation would be to measure $`(d\sigma /dq^2)(B\pi l\nu )`$ at $`q^2=m_\pi ^2`$ and to obtain $`|T|`$ using factorization.$`^\mathrm{?}`$ A value of $`|T|`$ obtained $`^\mathrm{?}`$ from $`B^+\pi ^+\pi ^0`$ implies $`\overline{}(B^0\pi ^+\pi ^{})|_{\mathrm{tree}}=(9.0\pm 1.8)\times 10^6`$ and thus would favor $`R_{\pi \pi }`$ significantly below 1. However, this work underestimates the uncertainty due to the color-suppressed contribution to $`B^+\pi ^+\pi ^0`$.$`^\mathrm{?}`$ A value of $`R_{\pi \pi }`$ below 1 would favor large $`\delta `$ (unexpected in factorization approaches) and, referring to Fig. 2, larger values of $`\varphi _2`$. Such a conclusion, in my opinion, is premature.

### 3.3 $`B^0\varphi K_S`$: New Physics?

In the decay $`B^0\varphi K_S`$, governed by the $`\overline{b}\overline{s}`$ penguin amplitude, the standard model predicts the same CP asymmetries as in those processes (like the “golden” $`J/\psi K_S`$ mode) governed by the $`\overline{b}\overline{s}c\overline{c}`$ tree amplitude. In both cases the weak phase is expected to be 0 (mod $`\pi `$), so the indirect CP asymmetry should be governed entirely by $`B^0`$$`\overline{B}^0`$ mixing and thus should be proportional to $`\mathrm{sin}2\varphi _1`$. There should be no direct CP asymmetries (i.e., one expects $`𝒜0`$) in either case. This is certainly true for $`BJ/\psi K`$; $`𝒜`$ is consistent with zero in the neutral mode, while the direct CP asymmetry is consistent with zero in the charged mode.$`^\mathrm{?}`$ However, a different result for $`B^0\varphi K_S`$ could point to new physics in the $`\overline{b}\overline{s}`$ penguin amplitude.$`^\mathrm{?}`$

The experimental situation for the asymmetries in $`B^0\varphi K_S`$ is shown in Table 2. We have included an updated result presented by the BaBar Collaboration at this conference.$`^\mathrm{?}`$ One cannot conclude much about $`𝒜_{\varphi K_S}`$ as a result of the substantial discrepancy between BaBar and Belle. However, the value of $`𝒮_{\varphi K_S}`$, which should be equal to $`\mathrm{sin}2\varphi _1=0.734\pm 0.054`$ in the standard model, is about $`2.7\sigma `$ away from it. If one assumes that the amplitudes for $`B^0\varphi K^0`$ and $`B^+\varphi K^+`$ are equal, as is true in many approaches, one expects the time-integrated CP asymmetry $`A_{CP}`$ in the charged mode to be equal to $`𝒜_{\varphi K_S}`$. The BaBar Collaboration $`^\mathrm{?}`$ has recently reported $`A_{CP}=0.039\pm 0.086\pm 0.011`$.

Many proposals for new physics can account for such a discrepancy.$`^\mathrm{?}`$ I describe a method similar to that $`^{\mathrm{?},\mathrm{?}}`$ used in analyzing $`B^0\pi \pi `$ for extracting a new physics amplitude, developed in collaboration with Cheng-Wei Chiang.$`^\mathrm{?}`$ One uses the measured values of $`𝒮_{\varphi K_S}`$ and $`𝒜_{\varphi K_S}`$ and the ratio

$$R_{\varphi K_S}\frac{\overline{}(B^0\varphi K_S)}{\overline{}(B^0\varphi K_S)|_{\mathrm{std}}}=1+2r\mathrm{cos}\varphi \mathrm{cos}\delta +r^2,$$
(8)
where $`r`$ is the ratio of the magnitude of the new amplitude to the one in the standard model, and $`\varphi `$ and $`\delta `$ are their relative weak and strong phases. For any values of $`R_{\varphi K_S}`$, $`\varphi `$, and $`\delta `$, Eq. (8) can be solved for the amplitude ratio $`r`$ and one then calculates

$`R_{\varphi K_S}𝒮_{\varphi K_S}`$ $`=`$ $`\mathrm{sin}2\varphi _1+2r\mathrm{cos}\delta \mathrm{sin}(2\varphi _1\varphi )+r^2\mathrm{sin}2(\varphi _1\varphi )`$ (9)
$`R_{\varphi K_S}𝒜_{\varphi K_S}`$ $`=`$ $`2r\mathrm{sin}\varphi \mathrm{sin}\delta .`$ (10)
The $`\varphi K_S`$ branching ratio in the standard model is calculated using the penguin amplitude from $`B^+K^0\pi ^+`$ and an estimate of electroweak penguin corrections. It was found $`^\mathrm{?}`$ that $`R_{\varphi K_S}=1.0\pm 0.2`$.

For values of $`\varphi `$ between $`\pi `$ and $`\pi `$, curves of $`𝒮_{\varphi K_S}`$ vs. $`𝒜_{\varphi K_S}`$ are plotted in Fig. 4 as $`\delta `$ varies from $`0`$ to $`\pi `$. Both $`𝒮`$ and $`𝒜`$ are unchanged under $`\varphi \varphi +\pi `$, $`\delta \delta \pi `$, while $`𝒮𝒮`$, $`𝒜𝒜`$ under $`\varphi \varphi +\pi `$, $`\delta \pi \delta `$.

Various regions of $`(\varphi ,\delta )`$ can reproduce the observed values of $`𝒮_{\varphi K_S}`$ and $`𝒜_{\varphi K_S}`$. Some of these are shown in Fig. 5, while others correspond to shifts in $`\varphi `$ and $`\delta `$ by $`\pm \pi `$. As errors on the observables shrink, so will the allowed regions. However, as Fig. 4 makes clear, there will always be a solution for some $`\varphi `$ and $`\delta `$ as long as $`R`$ remains compatible with 1. (The allowed regions of $`\varphi `$ and $`\delta `$ are restricted if $`R1`$.$`^\mathrm{?}`$) Typical values of $`r`$ are of order 1; one generally needs to invoke new-physics amplitudes comparable to those in the standard model.

The above scenario envisions new physics entirely in $`B^0\varphi K^0`$ and not in $`B^+K^0\pi ^+`$. An alternative is that new physics contributes to the $`\overline{b}\overline{s}`$ penguin amplitude and thus appears in both decays. Here it is convenient to define a ratio

$$R^{}\frac{\overline{\mathrm{\Gamma }}(B^0\varphi K^0)}{\overline{\mathrm{\Gamma }}(B^+K^0\pi ^+)},$$
(11)
where $`\overline{\mathrm{\Gamma }}`$ denotes a partial width averaged over a process and its CP conjugate. Present data indicate $`R^{}=0.78\pm 0.17`$. The $`B^0\varphi K^0`$ amplitude contains a contribution from both the gluonic and electroweak penguin terms, while $`B^+K^0\pi ^+`$ contains only the former. Any departure from the expected ratio of the electroweak to gluonic penguin amplitudes would signify new physics. Again, the central value of $`𝒮`$ would suggest this to be the case.$`^\mathrm{?}`$

### 3.4 $`B^0\eta ^{}K_S`$: No New Physics Needed

At present neither the rate nor the CP asymmetry in $`B\eta ^{}K`$ present a significant challenge to the standard model. The rate can be reproduced with the help of a modest contribution from a “flavor-singlet penguin” amplitude, the need for which was pointed out $`^{\mathrm{?},\mathrm{?}}`$ prior to the observation of this decay. One only needs to boost the standard penguin amplitude’s contribution by about 50% via the flavor-singlet term in order to explain the observed rate.$`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$. (Ref. $`^\mathrm{?}`$ instead finds an enchanced standard-penguin contribution to $`\eta ^{}`$ production.) The CP asymmetry is not a problem; the ordinary and singlet penguin amplitudes are expected to have the same weak phase Arg$`(V_{ts}^{}V_{tb})\pi `$ and hence one expects $`𝒮_{\eta ^{}K_S}\mathrm{sin}2\varphi _1`$, $`𝒜_{\eta ^{}K_S}0`$. The experimental situation is shown in Table 3. The value of $`𝒮_{\eta ^{}K_S}`$ is consistent with the standard model expectation at the $`1\sigma `$ level, while $`𝒜_{\eta ^{}K_S}`$ is consistent with zero.

The singlet penguin amplitude may contribute elsewhere in $`B`$ decays. It is a possible source of a low-effective-mass $`\overline{p}p`$ enhancement $`^\mathrm{?}`$ in $`B^+\overline{p}pK^+`$.$`^\mathrm{?}`$

## 4 $`B\pi K`$ Decays and $`\varphi _3=\gamma `$

The decays $`BK\pi `$ (with the exception of $`B^0K^0\pi ^0`$) are self-tagging. For example, the $`K^+\pi ^{}`$ final state is expected to originate purely from a $`B^0`$ and not from a $`\overline{B}^0`$. Since such self-tagging decays do not involve a CP eigenstate, one must contend with both weak and strong phases. Nonetheless several methods permit one to separate these from one another. We give two examples below.

### 4.1 $`B^0K^+\pi ^{}`$ vs. $`B^+K^0\pi ^+`$

The decay $`B^+K^0\pi ^+`$ is a pure penguin ($`P`$) process, while the amplitude for $`B^0K^+\pi ^{}`$ is proportional to $`P+T`$, where $`T`$ is a (strangeness-changing) tree amplitude. The ratio $`T/P`$ has magnitude $`r`$, weak phase $`\varphi _3\pm \pi =\gamma \pm \pi `$, and strong phase $`\delta `$. The ratio $`R_0`$ of these two rates (averaged over a process and its CP conjugate) is

$$R_0\frac{\overline{\mathrm{\Gamma }}(B^0K^+\pi ^{})}{\overline{\mathrm{\Gamma }}(B^+K^0\pi ^+)}=12r\mathrm{cos}\gamma \mathrm{cos}\delta +r^2\mathrm{sin}^2\gamma ,$$
(12)
where the inequality holds for any $`r`$ and $`\delta `$. If $`R_0`$ were significantly less than $`1`$ this inequality could be used to impose a useful constraint on $`\gamma `$.$`^\mathrm{?}`$ On the basis of the latest branching ratios from BaBar,$`^\mathrm{?}`$ Belle,$`^\mathrm{?}`$ and CLEO,$`^\mathrm{?}`$ summarized in Table 4, using the $`B^+/B^0`$ lifetime ratio $`\tau _+/\tau _0=1.073\pm 0.014`$,$`^\mathrm{?}`$, one finds $`R_0=0.99\pm 0.09`$, which is consistent with 1 and does not permit application of the bound. However, using additional information on $`r`$ and the CP asymmetry in $`B^0K^+\pi ^{}`$, one can obtain a constraint on $`\gamma `$.$`^{\mathrm{?},\mathrm{?}}`$

Define a “pseudo-asymmetry” normalized by the rate for $`B^0K^0\pi ^+`$, a process which is expected not to display a CP asymmetry since only the penguin amplitude contributes to it:

$$A_0\frac{\mathrm{\Gamma }(\overline{B}^0K^{}\pi ^+)\mathrm{\Gamma }(B^0K^+\pi ^{})}{2\overline{\mathrm{\Gamma }}(B^+K^0\pi ^+)}=2r\mathrm{sin}\gamma \mathrm{sin}\delta .$$
(13)
One may eliminate $`\delta `$ between this equation and Eq. (12) and plot $`R_0`$ as a function of $`\gamma `$ for the allowed range of $`|A_0|`$. One needs an estimate of $`r`$, whose present value (based on the branching ratios in Table 4 and arguments given in Refs. $`^{\mathrm{?},\mathrm{?}}`$) is $`r=0.17\pm 0.04`$. Here one must take the value of $`T`$ from $`B\pi l\nu `$,$`^\mathrm{?}`$ using flavor SU(3) to relate the strangeness-preserving and strangeness-changing terms. The latest BaBar and Belle data imply $`A_0=0.088\pm 0.040`$, leading us to take $`|A_0|0.13`$ at the $`1\sigma `$ level. Curves for $`A_0=0`$ and $`|A_0|=0.13`$ are shown in Fig. 6. The lower limit $`r=0.134`$ is used to generate these curves since the limit on $`\gamma `$ will be the most conservative.

At the $`1\sigma `$ level, using the constraints that $`R_0`$ must lie between 0.90 and 1.08 and $`|A_0|`$ must lie between zero and 0.13, one can establish the bound $`\gamma \stackrel{>}{}60^{}`$. No bound can be obtained at the 95% confidence level, however. Despite the impressive improvement in experimental precision (a factor of 2 decrease in errors since the analysis of Ref. $`^\mathrm{?}`$), further data are needed in order for a useful constraint to be obtained.

### 4.2 $`B^+K^+\pi ^0`$ vs. $`B^+K^0\pi ^+`$

The comparison of rates for $`B^+K^+\pi ^0`$ and $`B^+K^0\pi ^+`$ also can give information on $`\varphi _3=\gamma `$. The amplitude for $`B^+K^+\pi ^0`$ is proportional to $`P+T+C`$, where $`C`$ is a color-suppressed amplitude. Originally it was suggested that this amplitude be compared with $`P`$ from $`B^+K^0\pi ^+`$ and $`T+C`$ taken from $`B^+\pi ^+\pi ^0`$ using flavor SU(3),$`^\mathrm{?}`$ using a triangle construction to determine $`\gamma `$. However, electroweak penguin amplitudes contribute significantly in the $`T+C`$ term.$`^\mathrm{?}`$ It was noted subsequently $`^\mathrm{?}`$ that since the $`T+C`$ amplitude corresponds to isospin $`I(K\pi )=3/2`$ for the final state, the strong-interaction phase of its EWP contribution is the same as that of the rest of the $`T+C`$ amplitude, permitting the calculation of the EWP correction.

New data on branching ratios and CP asymmetries permit an update of previous analyses.$`^{\mathrm{?},\mathrm{?}}`$ One makes use of the ratios

$$R_c\frac{2\overline{\mathrm{\Gamma }}(B^+K^+\pi ^0)}{\overline{\mathrm{\Gamma }}(B^+K^0\pi ^+)}=1.30\pm 0.15,r_c=|(T+C)/P|=0.20\pm 0.02,$$
(14)
and

$$A_c\frac{\mathrm{\Gamma }(B^{}K^{}\pi ^0)\mathrm{\Gamma }(B^+K^+\pi ^0)}{\overline{\mathrm{\Gamma }}(B^+K^0\pi ^+)}=0.05\pm 0.09.$$
(15)
One must also use an estimate $`^\mathrm{?}`$ of the electroweak penguin parameter $`\delta _{\mathrm{EW}}=0.65\pm 0.15`$. One obtains the most conservative (i.e., weakest) bound on $`\gamma `$ for the maximum values of $`r_c`$ and $`\delta _{\mathrm{EW}}`$.$`^\mathrm{?}`$ The resulting plot is shown in Fig. 7. One obtains a bound at the $`1\sigma `$ level very similar to that in the previous case: $`\gamma \stackrel{>}{}58^{}`$.

## 5 Summary

The process $`B^0J/\psi K_S`$ has provided spectacular confirmation of the Kobayashi-Maskawa theory of CP violation, measuring $`\varphi _1=\beta `$ to a few degrees. Now one is entering the territory of more difficult measurements.

The decay $`B^0\pi ^+\pi ^{}`$ has great potential for giving useful information on $`\varphi _2=\alpha `$. One needs either a measurement of $`(B^0\pi ^0\pi ^0)`$,$`^\mathrm{?}`$ probably at the $`10^6`$ level (present limits $`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$ are several times that), or a better estimate of the tree amplitude from $`B\pi l\nu `$.$`^\mathrm{?}`$ As for the BaBar and Belle experimental CP asymmetries,$`^{\mathrm{?},\mathrm{?}}`$ they will eventually converge to one another, as did the initial measurements of $`\mathrm{sin}2\varphi _1`$ using $`B^0J/\psi K_S`$.

The $`B\varphi K_S`$ decay can display new physics via special $`\overline{b}\overline{s}s\overline{s}`$ operators or effects on the $`\overline{b}\overline{s}`$ penguin. Some features of any new amplitude can be extracted from the data in a model-independent way if one uses both rate and asymmetry information.$`^\mathrm{?}`$

The rate for $`B\eta ^{}K_S`$ is not a problem for the standard model if one allows for a modest flavor-singlet penguin contribution in addition to the standard penguin amplitude. The CP asymmetries for this process are in accord with the expectations of the standard model at the $`1\sigma `$ level or better. Effects of the singlet penguin amplitude may also be visible elsewhere, for example in $`B^+p\overline{p}K^+`$.

Various ratios of $`BK\pi `$ rates, when combined with information on CP asymmetries, show promise for constraining phases in the CKM matrix. These tests have shown a steady improvement in accuracy since the asymmetric $`B`$ factories have been operating, and one expects further progress as the instantaneous and accumulated luminosities increase. In the longer term, hadron colliders may provide important contributions.

## Acknowledgments

I thank C.-W. Chiang, M. Gronau, and Z. Luo for collaborations on the subjects mentioned here; Tom Browder, Steve Olsen, Sandip Pakvasa, Xerxes Tata, and San Fu Tuan for discussions, and the Physics Department of the University of Hawaii for hospitality. This work was supported in part by the United States Department of Energy through Grant No. DE FG02 90ER40560.

## References


text:
# Traffic flow of interacting self-driven particles: rails and trails, vehicles and vesicles

## I Introduction

In the recent years non-equilibrium statistical mechanics has found unusual applications in research on traffic flow of various different types of objects. In this paper we consider mainly three different examples of such traffic, namely, (a) vehicular traffic css ; helbing ; naga , (b) ant-traffic on ant-trails cgns ; ncs ; burd ; couzin , and (c) intra-cellular traffic of molecular motors carrying vesicular cargo moving along cytoskeletal filaments frey . Most of these models are essentially generalizations or extensions of the Asymmetric Simple Exclusion Process (ASEP) spohn , which is, to our knowledge, the simplest model of systems consisting of interacting driven particles; the general aim of these investigations is to understand the interplay of self-organized structures and transport in systems driven far from equilibrium zia ; schutz ; stinchcombe .

The aim of this article is (a) to summarize the main results of recent works on all the three systems mentioned above, elucidating the nature of various types of quenched randomness, (b) to present the challenging open problems, and (c) to indicate the possible trends of future developments in this frontier area of interdisciplinary research.

The common modeling strategy is to represent the motile elements (i.e., vehicles, ants, molecular motors) by “self-propelled” particles which convert chemical energy (derived from fuel or food) into the mechanical energy required for the forward movement. In such generic models, the mutual influences of the motile elements on the movements of each other are captured by appropriate inter-particle interactions. In the spirit of the lattice gas models, the track for the traffic movement (i.e., the highway lane or ant trail or cytoskeletal filaments) are represented as discrete lattice of “cells” each of which can accomodate at most one partcle at a time. The dynamical laws governing the forward movement of the self-propelled particles in such “particle-hopping” models are usually formulated as “update rules” in terms of cellular automata (CA) wolfram .

## II Vehicular traffic

The simplest model of interacting self-driven particles is the so-called totally asymmetric simple exclusion process (TASEP) spohn ; zia ; schutz . Imposition of open boundary conditions leads to richer physics as compared to those for the corresponding model with periodic boundary conditions. The states of the system are updated either in parallel or in a random-sequential manner following rules which will be explained later in this section.

### II.1 TASEP with periodic boundary conditions

TASEP with periodic boundary conditions is sketched schematically in fig.1. In the original formulation of TASEP, the states of the system were updated in a random sequential manner where a particle is picked up randomly and moved forward by one lattice spacing, with the hopping probability $`q`$, provided the site in front of the particle is empty.

The Nagel-Schreckenberg (NS) model ns  is a minimal CA model of vehicular traffic on idealized single-lane highways; the maximum possible (discrete) speed of the vehicles is $`V_{max}`$. However, in the special case $`V_{max}=1`$ this model reduces to the TASEP with parallel updating.

The exact fundamental diagrams of the TASEP with random-sequential updating and parallel updating are given by

$$F_r=qc(1c)$$
(1)
and

$$F_p=\frac{1}{2}\left[1\sqrt{14qc(1c)}\right],$$
(2)
respectively; the two expressions (1) and (2), for identical hopping probability $`q=0.75`$, are compared in fig.2. Note that both the expressions (1) and (2) exhibit particle-hole symmetry, i.e., these are symmetric under the interchange of $`c`$ and $`1c`$. In section III we shall show how this symmetry is broken in a model of ant-traffic on ant-trails.

### II.2 TASEP with open boundary conditions

The open boundary condition is, however, closer to the real vehicular traffic on a stretch of highway. If open boundary conditions are imposed on the TASEP, additional rules must be specified to regulate the entry and exit of the particles at the two boundaries of the finite system. Usually, these are specified as follows: if the first site at the open point of entry is empty it is filled with probability $`\alpha `$ whereas particles occupying the last site at the point of exit hop out of the system with probability $`\beta `$ (see fig.3).

The open boundaries break the translational invariance of the system and give rise to stationary states with non-trivial density profiles. Such model systems have been investigated thoroughly over the last decade from the point of view of fundamental principles of non-equilibrium statistical mechanics. In contrast to equilibrium systems with short-range interactions, such driven non-equilibrium systems can exhibit transitions from one dynamical phase to another, even in one-dimension with only short-range interactions, with the slight change of boundary conditions krug ; css .

The typical phase diagrams of the TASEP with open boundary conditions are sketched in fig.4; the qualitative features of the phase diagram of TASEP is practically independent of the nature of the dynamics. In the low-density phase A the flux is independent of $`\beta `$ and limited only by $`\alpha `$. On the other hand, in the high-density phase B the flux is independent of $`\alpha `$ and determined by $`\beta `$. However, in the maximum flux phase C the current is independent of both $`\alpha `$ and $`\beta `$. Moreover, both the high- and low-density phases can be subdivided into two phases each, namely, AI, AII and BI, BII, respectively; these subphases are distinguished by the asymptotic behaviour of the density profiles at the boundaries (see the insets in fig.4). In section IV we shall see how this phase diagram is modified in models of intra-cellular transport of vesicular cargo by molecular motors.

## III Ant traffic on trails

Ants communicate with each other by dropping a chemical (generically called pheromone) on the substrate as they crawl forward wilson ; camazine . Although we cannot smell it, the trail pheromone sticks to the substrate long enough for the other following sniffing ants to pick up its smell and follow the trail. In our recent papers cgns ; ncs  we have developed a particle-hopping model, formulated in terms of stochastic CA, which may be interpreted as a model of unidirectional ant-traffic on a trail. Rather than addressing the question of the emergence of the ant-trail, we have focussed attention on the traffic of ants on a trail which has already been formed.

Each site of our one-dimensional ant-trail model represents a cell that can accomodate at most one ant at a time (see Fig. 5). The lattice sites are labelled by the index $`i`$ ($`i=1,2,\mathrm{},L`$); $`L`$ being the length of the lattice. We associate two binary variables $`S_i`$ and $`\sigma _i`$ with each site $`i`$ where $`S_i`$ takes the value $`0`$ or $`1`$ depending on whether the cell is empty or occupied by an ant. Similarly, $`\sigma _i=1`$ if the cell $`i`$ contains pheromone; otherwise, $`\sigma _i=0`$. Thus, in contrast to TASEP, we have two subsets of dynamical variables in this model, namely, $`\{S(t)\}`$ and $`\{\sigma (t)\}`$. The instantaneous state (i.e., the configuration) of the system at any time is specified completely by the set $`(\{S\},\{\sigma \})`$.

Since a unidirectional motion is assumed, ants do not move backward. The forward-hopping probability of an ant is higher if it smells pheromone ahead of it. The state of the system is updated at each time step in two stages. In our ant-trail model with parallel dynamics, at each stage the dynamical rules are applied in parallel to all ants and pheromones, respectively.

Stage I: Motion of ants
An ant in cell $`i`$ that has an empty cell in front of it, i.e., $`S_i(t)=1`$ and $`S_{i+1}(t)=0`$, hops forward with

$$\mathrm{probability}=\{\begin{array}{ccc}Q\hfill & \mathrm{if}\sigma _{i+1}(t)=1,\hfill & \\ q\hfill & \mathrm{if}\sigma _{i+1}(t)=0,\hfill & \end{array}$$
(3)
where, to be consistent with real ant-trails, we assume $`q<Q`$.

Stage II: Evaporation of pheromones
At each cell $`i`$ occupied by an ant after stage I a pheromone will be created, i.e.,

$$\sigma _i(t+1)=1\mathrm{if}S_i(t+1)=1.$$
(4)
On the other hand, any ‘free’ pheromone at a site $`i`$ not occupied by an ant will evaporate with the probability $`f`$ per unit time, i.e., if $`S_i(t+1)=0`$, $`\sigma _i(t)=1`$, then

$$\sigma _i(t+1)=\{\begin{array}{ccc}0\hfill & \mathrm{with}\mathrm{probability}f,\hfill & \\ 1\hfill & \mathrm{with}\mathrm{probability}1f.\hfill & \end{array}$$
(5)
We have also considered another version of our ant-trail model where the states of the system are updated in a random-sequential manner rather than in parallel ncs . Note that in both the cases, because of the periodic boundary conditions, the dynamics conserves the number $`N`$ of ants, but not the number of pheromones.

This model is related to several other models. For example, in the limits $`f0`$ and $`\mathrm{}`$ this ant-trail model reduces to TASEP with the hopping probabilities $`q`$ and $`Q`$, respectively. Moreover, the ant-trail model may be regarded as the opposite limit of the bus-route model loan ; cd  (see ref.ncs  for the detailed comparison). Furthermore, the ant-trail model also has some similarities with the particle-hopping models of human trails of pedestrians helbped ; burs .

The typical fundamental diagrams of the ant-trail model cgns  with parallel dynamics are shown in fig.6; the corresponding results in the case of random-sequential updating are qualitatively similar ncs . The unusual shapes of the curves observed over a range of $`f`$ are consequences of the non-monotonic variation of the average speed of the ants with their density on the trail (see fig.7).

Both the ordinary mean-field theory (MFT), which accounts for the exact fundamental diagram of the TASEP with random-sequential updating, and $`2`$-cluster MFT ssni ; css , which succussfully predicted the exact fundamental diagram of the TASEP with parallel updating, fail to capture even the qualitative features of the fundamental diagrams of the ant-trail model shown in fig.6 ncs . However, a heuristic MFT, described in ref.ncs , captures at least the qualitative features of the observed flow properties of our ant-trail model.

In order to develope a quantitative theory for the flow properties of the ant-trail model, we have analyzed the spatial organization of the ants by computer simulations. Analyzing these observations we concluded that in the anomalous regime, loose clusters of ants dominate; the term “loose” means that there are small gaps in between successive ants in the cluster although the cluster appears to be an usual compact cluster if seen from a distance. As shown in fig.8, the fundamental diagram we calculated within the “loose”-cluster approximation (LCA) is in good quantitative agreement with the corresponding data we obtained from computer simulations of the model.

## IV Intra-cellular traffic of molecular motors

Molecular motors are protein molecules that convert the chemical energy, released by the hydrolysis of ATP, into mechanical energy required for its forward movement during intra-cellular transport of vesicular cargo howard . The minimal models developed for explaining the mechanism of directed motion of isolated motor proteins are based on Brownian ratchets reimann . In such models, each motor is represented by a particle. The essential features of the detailed mechano-chemistry of the molecular motor is captured in the Brownian ratchet models by a stochastic sequence of successive attachments and detachments of the motor with the cytoskeletal filamentary track (e,g, microtubule in the case of kinesin and dynein motors). In the simplest versions of these models julicher , in the attached state, the particle representing a motor is subjected to a potential that is spatially periodic, but each period of which is asymmetric. In the detached state the particle executes an unbiased diffusive motion. In spite of its simplicity, such a minimal model can account for the directed, albeit noisy, movement of individual isolated motors.

To our knowledge, the question of the effects of interactions of the motors on the intra-cellular traffic was addressed theoretically for the first time only a few years ago menon . In that work, the filamentary track was discretized in the spirit of the particle-hopping models described above and the motors were represented by field-driven particles. Both forward and backward movement of the particles were possible and the hopping probability of every particle was computed from the local potential. Thus, this model was a generalization of ASEP rather than TASEP where the hopping probabilities were obtained from the local potential which itself was time-dependent. The fundamental diagram of that model, computed imposing periodic boundary conditions, is very similar to those shown in fig.2. This observation indicates that further simplification in the model proposed in ref.menon  is possible to develope a minimal model for interacting molecular motors.

Recently, Parmeggiani et al.frey  have, indeed, developed such a minimal model for interacting molecular motors involved in intra-cellular transport by extending the TASEP with open boundary conditions. In this model, the molecular motors (e.g., kinesin or dynein) are represented by particles whereas the sites for the binding of the motors with the cytoskeletal tracks (e.g., microtubules) are represented by a one-dimensional discrete lattice. Just as in TASEP, the motors are allowed to hop forward, with probability $`q`$, provided the site in front is empty. However, unlike TASEP, the particles can also get “attached” to an empty lattice site, with probability $`A`$, and “detached” from an occupied site, with probability $`D`$ (see fig.9) from any site except the end points. The state of the system was updated in a random-sequential manner.

To my knowledge, this is the first application of TASEP to intra-cellular transport phenomena although it is not the first application of TASEP in the domain of biological systems; for example, a TASEP-like model was considered earlier for protein synthesis protein

Carrying out Monte-Carlo simulations Parmeggiani et al.frey  demonstrated a novel phase where low and high density regimes, separated from each other by domain walls, coexist (see fig.10). Using a MFT, they interpreted this spatial organization as traffic jam of molecular motors.

## V Defects and disorder in particle-hopping models

At least three different types of defects and quenched randomness have been considered so far in the context of the models of interacting particles driven far from equilibrium. (a) First, the randomness may be associated with the track on which the particles move; typical examples are the bottlenecks on the roads (in the context of vehicular traffic) or defects on the microtubules (in intra-cellulat transport), etc. For example, as shown in fig.11(a), normal hopping probability at unblocked sites is $`Q`$ whereas that at the bottleneck is $`q`$ ($`q<Q`$). This type of quenched defect and disorder of the track leads to interesting phase-segregation phenomena (see css  for a review).

(b) The second type of randomness is associated with the hopping particles, rather than with the track. For example, the normal hopping probabilities of the particles may vary randomly from one particle to another (see fig.11(b)); the hopping probabilities are, however, “quenched” random variables, i.e., independent of time. In this case, the system is known to be exhibit coarsening of queues of the particles and the phenomenon has some formal similarities with Bose-Einstein condensation (reviewed in css ).

Note that in case of the randomness of type (a), the hopping probability depends only on the spatial location on the track, independent of the identity of the hopping particle. On the other hand, in the case of randomness of type (b), the hopping probability depends on the hopping particle, irrespective of its spatial location on the track. In contrast to these two types of randomness, the randomness in the hopping probabilities of the particles in some models arises from the coupling of their dynamics with that of another non-conserved dynamical variable. For example, in the ant-trail model, the hopping probability of an ant depends on the presence or absence of pheromone in front of it (see fig.11(c)). Therefore, in such models with periodic boundary conditions, a given particle may hop from the same site, at different times, with different hopping probabilities.

Defects of either (a) the cytoskeletal filaments or (b) the motor proteins or (c) the mechano-chemical coupling can cause malfunctiong of the intra-cellular transport leading to various types of diseases trafdis . In order to get deep insight into the physical origin of such diseases, the recent model developed by Parmeggiani et al. frey  has been extended chowetal . This modeling strategy has opened up a new horizon for further unconventional applications of non-equilibrium statistical mechanics far beyond the traditional borderlines of physics.

## VI Conclusion

In this paper we began with a brief introduction to TASEP zia ; schutz  which is, perhaps, the simplest model of systems of interacting driven particles. TASEP, when updated in parallel, may be regarded as a special case, corresponding to the maximum allowed speed $`V_{max}=1`$, of the NS model ns , the minimal model of vehicular traffic on single-lane highways. We have summarized some of the main known results on the fundamental diagram for the TASEP with periodic boundary conditions and the phase diagram of the TASEP with open boundary conditions. Then, we have shown how these results for TASEP get qualitatively modified by the generalizations or extensions required to model ant-traffic on ant-trails cgns ; ncs  and molecular motor traffic on cytoskeletal filaments frey .

In the context of the ant-traffic on ant-trails, we have established how a combination of analytical and numerical calculations cgns ; ncs  can account for the unusual shape of the fundamental diagram observed in the computer simulations of the ant-trail model. We have also presented the phase diagram obtained by Parmeggiani et al.frey  from studies of their recent model for molecular motor traffic. This phase diagram suggests the possibility of coexistence of high-density regions (traffic jam) and low-density regions (freely flowing traffic), separated from each other by domain walls, in a novel phase. Finally, we have mentioned some ongoing investigations on the effects of defects and disorder on molecular motor traffic chowetal . This trend of research indicates the possibility of further unconventional, but very useful, applications of statistical physics in biological systems.

Acknowledgements: I dedicate this paper to Dietrich Stauffer on the occassion of his 60th birthday. During the two decades of our collaboration he not only taught me the art of computer simulations but also inspired me to try unconventional applications of the conventional tools of statistical physics. I thank E. Frey, V. Guttal, F. Jülicher, A. Kunwar, K. Nishinari and A. Schadschneider for enjoyable recent collaborations and/or discussions. I also thank E. Frey for his kind permission to reproduce figure 10 from his paper frey . Part of this work was supported by the Alexander von Humboldt Foundation (under the re-invitation program for former fellows) and by the German Research Foundation (DFG) through a joint Indo-German research project.


text:
# Nonlocal pseudopotentials and magnetic fields

## Abstract

We show how to describe the coupling of electrons to non-uniform magnetic fields in the framework of the widely used norm-conserving pseudopotential approximation for electronic structure calculations. Our derivation applies to magnetic fields that are smooth on the scale of the core region. The method is validated by application to the calculation of the magnetic susceptibility of molecules. Our results are compared with high quality all electron quantum chemical results, and another recently proposed formalism.

The coupling of the electrons in matter with probing electromagnetic fields or charged particles provides the basis for nearly all known analytical experimental techniques. However, the most computationally efficient schemes for the first principles prediction (and hence interpretation) of the experimental observables require the use of approximate, and crucially, nonlocal Hamiltonians. Most notably, the use of nonlocal pseudopotentials, along with with density functional theory, is often referred to as the Standard Model of modern electronic structure theory. This is not without justification. The computational efficiency, and accuracy of modern first principles pseudopotentials has allowed a vast range of problems to be solved for realistic materials. However, whatever the successes of this method, it is still unclear how nonlocal Hamiltonians should be coupled to arbitrary magnetic fields. In the specific case of uniform magnetic fields described with the symmetric-gauge vector-potential, we have derived the correct Hamiltonian by developing what we called the gauge including projector augmented wave (GIPAW) method Pickard and Mauri (2001). In this Letter we derive the pseudopotential Hamiltonian that describes the coupling between the electrons and a smooth non-uniform magnetic field represented by an arbitrary vector-potential gauge.

We are not alone in this quest. Recently Ismail-Beigi, Chang and Louie Ismail-Beigi et al. (2001), to whom we refer the reader for a more complete summary of earlier work in this field, proposed a scheme (hereafter referred to as the ICL method) which sought to finally provide a rigorous derivation of a closed form for the coupling of nonlocal systems to arbitrary electromagnetic fields. While the widespread use of nonlocal pseudopotentials clearly provided the motivation and applications for their method, they attempted to tackle the problem of a general nonlocal Hamiltonian. In contrast, this Letter focuses exclusively on the class of nonlocal Hamiltonians that arise due to the use of first-principles nonlocal pseudopotentials for electronic structure calculations.

In the pseudopotential approach, and in the absence of a magnetic field, the all-electron (AE) Hamiltonian $`H^{\mathrm{AE}}=p^2/2+V(𝐫)`$, is replaced by its pseudo (PS) equivalent, $`H^{\mathrm{PS}}=p^2/2+V^\mathrm{l}(𝐫)+_𝐑V_𝐑^{\mathrm{nl}}`$. $`V^\mathrm{l}(𝐫)`$ is a local-potential, and $`V_𝐑^{\mathrm{nl}}`$ is a nonlocal operator which acts only within the core region of the atomic sites, at $`𝐑`$:

$$V_𝐑^{\mathrm{nl}}=d^3r^{}d^3r^{\prime \prime }|𝐫^{}𝐫^{\prime \prime }|V_𝐑^{\mathrm{nl}}(𝐫^{},𝐫^{\prime \prime }).$$
(1)
By construction, in the valence-energy-range and to within some controllable error, (i) the eigenvalues of $`H^{\mathrm{PS}}`$ coincide with those of $`H^{\mathrm{AE}}`$, and (ii) the eigenstates of $`H^{\mathrm{PS}}`$ coincide with those of $`H^{\mathrm{AE}}`$ outside the core regions.

Turning on a magnetic field, $`𝐁(𝐫)=\times 𝐀(𝐫)`$, the AE Hamiltonian becomes:

$$H_𝐀^{\mathrm{AE}}=\frac{1}{2}\left[𝐩+\frac{1}{c}𝐀(𝐫)\right]^2+V(𝐫).$$
(2)
The question that we aim to answer in this Letter is: what is the PS Hamiltonian that, in presence of a magnetic field, satisfies the requirements (i) and (ii) stated above? The PS Hamiltonian can be written as:

$$H_𝐀^{\mathrm{PS}}=\frac{1}{2}\left[𝐩+\frac{1}{c}𝐀(𝐫)\right]^2+V^\mathrm{l}(𝐫)+\underset{𝐑}{}V_𝐑^{\mathrm{nl}}+\mathrm{\Delta }H_𝐀,$$
(3)
where $`\mathrm{\Delta }H_𝐀`$ remains to be determined. $`\mathrm{\Delta }H_𝐀`$ cannot be zero for all $`𝐀(𝐫)`$, as $`H_𝐀^{\mathrm{PS}}`$ must be gauge invariant to satisfy (i) and (ii). Its eigenvalues must not depend on the arbitrary choice of the gauge of $`𝐀(𝐫)`$. However, as we will demonstrate, demanding gauge invariance alone is not sufficient to uniquely determine $`\mathrm{\Delta }H_𝐀`$.

To obtain $`H_𝐀^{\mathrm{PS}}`$ we use Blöchl’s projector augmented wave (PAW) theory Blöchl (1994). In the PAW approach the pseudisation procedure is defined as a linear transformation between Hilbert spaces – those of the AE valence wavefunctions, and the computationally convenient PS wavefunctions. This transformation, $`|\mathrm{\Psi }=𝒯|\stackrel{~}{\mathrm{\Psi }}`$, can be applied to obtain PS operators, $`\stackrel{~}{𝒪}`$, which correspond to their AE counterparts, $`𝒪`$:

$`\stackrel{~}{𝒪}`$ $`=`$ $`𝒯^{}𝒪𝒯=𝒪+\stackrel{~}{𝒞}`$ (4)
$`\stackrel{~}{𝒞}`$ $`=`$ $`{\displaystyle \underset{i,j}{}}|\stackrel{~}{p}_i\left[\varphi _i|𝒪|\varphi _j\stackrel{~}{\varphi }_i|𝒪|\stackrel{~}{\varphi }_j\right]\stackrel{~}{p}_j|,`$ (5)
where we have adopted Blöchl’s tilde to denote a PS quantity. By construction, the expectation values of $`\stackrel{~}{𝒪}`$ between PS wavefunctions are equal those of $`𝒪`$ between the corresponding AE wavefunctions. The $`|\varphi _i`$ and $`|\stackrel{~}{\varphi }_i`$ are atomic all electron and pseudo partial waves respectively. The projector functions $`|\stackrel{~}{p}_i`$ act only within some augmentation region – or core radius in the language of pseudopotential theory. The AE and PS partial wave coincide outside this augmentation region, and $`\stackrel{~}{p}_i|\stackrel{~}{\varphi }_i=\delta _{ij}`$. If the norms of the AE and PS partial waves are equal and there is just one partial wave in each angular momentum channel, which we shall assume in the following, then taking Eq. (4) with $`𝒪=H^{\mathrm{AE}}`$ and $`𝐁(𝐫)=\mathrm{𝟎}`$, we obtain the $`H^{\mathrm{PS}}`$ of norm-conserving pseudopotentials in the Kleinman-Bylander form Kleinman and Bylander (1982).

In principle, in order to compute expectation values of physical observables from PS wavefunctions one should use the PS operator $`\stackrel{~}{𝒪}`$. In practice, since the AE and PS partial waves are identical outside the augmentation regions and have the same norm inside, the term $`\stackrel{~}{𝒞}`$ in Eq. (4) can be neglected for a large class of operators. For example, the AE perturbation Hamiltonian describing the coupling with a uniform electric field $`𝐄`$ is $`V_𝐄^{\mathrm{AE}}(𝐫)=𝐫𝐄`$. In norm-conserving PS calculations of the electric-field response-properties, the $`\stackrel{~}{𝒞}`$ term is always neglected, with $`V_𝐄^{\mathrm{AE}}(𝐫)`$ being used as the approximate PS perturbation potential Baroni and Resta (1986); de Gironcoli et al. (1989); King-Smith and Vanderbilt (1993). The $`\stackrel{~}{𝒞}`$ term can also be neglected for operators for which the weight is concentrated away from the atoms or are almost constant in the augmentation regions. This is the case for operators of the form $`[𝐜(𝐫𝐑)]^n`$ in systems with a single PS atom centered at $`𝐑`$. Here, and in the following, $`𝐜`$ is a constant vector and $`n0`$. In addition, the $`\stackrel{~}{𝒞}`$ term is exactly zero for the $`𝐑`$-centered angular-momentum operator $`𝐋_𝐑=(𝐫𝐑)\times 𝐩`$, as the partial waves are eigenstates of $`|𝐋_𝐑|^2`$ and of $`(𝐋_𝐑)_z`$. Following this reasoning the $`\stackrel{~}{𝒞}`$ term is also negligible for operators of the form $`𝐋_𝐑[𝐜(𝐫𝐑)]^n`$. In contrast, $`\stackrel{~}{𝒞}`$ can not be neglected for the kinetic energy operator, $`p^2/2`$, or for operators that are general functions of the $`𝐩`$ and $`𝐫`$ operators. Indeed, by construction, the PS wavefunctions can be expanded on a much smaller set of Fourier components in the momentum space (where $`𝐩`$ is diagonal) than the corresponding AE wavefunctions.

In general, because of the presence of the $`𝐩`$ operator, the $`\stackrel{~}{𝒞}`$ term corresponding to the AE magnetic perturbation Hamiltonian $`H_𝐀^{\mathrm{AE}}H^{\mathrm{AE}}`$,

$$\mathrm{\Delta }H_𝐀^{\mathrm{AE}}=\frac{𝐩𝐀(𝐫)+𝐀(𝐫)𝐩}{2c}+\frac{A(𝐫)^2}{2c^2},$$
(6)
can not be neglected, even if $`𝐁(𝐫)`$ and $`𝐀(𝐫)`$ are smooth on the scale of the core augmentation regions. However, as we shall show, if $`𝐁(𝐫)`$ is smooth, and we consider a system with a single PS atom centered at $`𝐑`$, there is a special vector potential gauge, $`𝐀^{}(𝐫)`$, in which the $`\stackrel{~}{𝒞}`$ term can be neglected. This $`𝐀^{}(𝐫)`$ potential can be defined in terms of the Fourier components $`𝐛_𝐆`$ of the magnetic field,

$$𝐁(𝐫)=\underset{𝐆}{\overset{|𝐆|<G_{\mathrm{max}}}{}}𝐛_𝐆e^{i𝐆𝐫},$$
(7)
where $`𝐛_𝐆𝐆=0`$, and, if the field is smooth on the scale of the core radius $`r_\mathrm{c}`$, then $`r_\mathrm{c}G_{\mathrm{max}}1`$. In particular,

$$𝐀^{}(𝐫)=\underset{𝐆}{\overset{|𝐆|<G_{\mathrm{max}}}{}}\frac{1}{2}𝐛_𝐆\times (𝐫𝐑)f[i𝐆(𝐫𝐑)],$$
(8)
with $`f(x)`$ defined by,

$$f(x)=2\frac{1+xe^xe^x}{x^2}=1+\frac{2}{3}x+O(x^2).$$
(9)
Note that if the field is smooth, in the core region $`x1`$ and $`f(x)`$ can be expanded in powers of $`x`$. Since, for a uniform magnetic field $`𝐀^{}(𝐫)`$ reduces to the symmetric gauge, $`𝐀^{}(𝐫)`$ can be seen as a generalization of the symmetric gauge to non-uniform magnetic fields.

The magnetic coupling Hamiltonian in the $`𝐀^{}(𝐫)`$ gauge is:

$$\mathrm{\Delta }H_𝐀^{}^{\mathrm{AE}}=\underset{𝐆}{\overset{|𝐆|<G_{\mathrm{max}}}{}}\frac{𝐛_𝐆𝐋_𝐑f[i𝐆(𝐫𝐑)]+cc}{4c}+\frac{A^{}(𝐫)^2}{2c^2},$$
(10)
where $`cc`$ stands for the complex conjugate. All the terms in the asymptotic expansion of $`𝐋_𝐑f[x]`$ are of the form $`𝐋_𝐑[𝐜(𝐫𝐑)]^n`$ with $`n0`$. As a result, if $`𝐁(𝐫)`$ is smooth on the scale of the augmentation region, then $`\stackrel{~}{𝒞}`$ term arising from the $`\mathrm{\Delta }H_𝐀^{}^{\mathrm{AE}}`$ operator can be neglected. Thus the total PS Hamiltonian in the special $`𝐀^{}(𝐫)`$ gauge is:

$$H_𝐀^{}^{\mathrm{PS}}=\frac{1}{2}\left[𝐩+\frac{1}{c}𝐀^{}(𝐫)\right]^2+V^\mathrm{l}(𝐫)+V_𝐑^{\mathrm{nl}}.$$
(11)
From this result we can show that the Hamiltonian for an arbitrary gauge $`𝐀(𝐫)`$ is:

$`H_𝐀^{\mathrm{PS}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[𝐩+{\displaystyle \frac{1}{c}}𝐀(𝐫)]^2+V^\mathrm{l}(𝐫)+{\displaystyle }d^3r^{}d^3r^{\prime \prime }|𝐫^{}𝐫^{\prime \prime }|\times `$ (12)
$`V_𝐑^{\mathrm{nl}}(𝐫^{},𝐫^{\prime \prime })e^{\frac{i}{c}_{𝐫^{}𝐑𝐫^{\prime \prime }}𝑑𝐫𝐀(𝐫)},`$
where $`𝐫𝐫^{}`$ indicates a straight line path from point $`𝐫`$ to point $`𝐫^{}`$. To prove Eq. (12) we have just to notice that $`H_𝐀^{\mathrm{PS}}`$ is gauge invariant and reduces to Eq. (11) for $`𝐀(𝐫)=𝐀^{}(𝐫)`$.

Finally, the linearity of the Hamiltonian with respect to the electron-ion potential can be exploited to obtain $`H_𝐀^{\mathrm{PS}}`$ when there is more the one PS atom:

$`H_𝐀^{\mathrm{PS}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[𝐩+{\displaystyle \frac{1}{c}}𝐀(𝐫)]^2+V^\mathrm{l}(𝐫)+{\displaystyle }d^3r^{}d^3r^{\prime \prime }|𝐫^{}𝐫^{\prime \prime }|\times `$ (13)
$`{\displaystyle \underset{𝐑}{}}V_𝐑^{\mathrm{nl}}(𝐫^{},𝐫^{\prime \prime })e^{\frac{i}{c}_{𝐫^{}𝐑𝐫^{\prime \prime }}𝑑𝐫𝐀(𝐫)}.`$
We refer to $`H_𝐀^{\mathrm{PS}}`$ as the GIPAW Hamiltonian for an arbitrary magnetic field, since $`H_𝐀^{\mathrm{PS}}`$ reduces to the Hamiltonian that we have derived in our earlier work Pickard and Mauri (2001), if the magnetic field is uniform and if the symmetric gauge is used. Our new Hamiltonian holds if the magnetic field varies smoothly over the core region, and if the potentials are norm conserving. If the field varies more rapidly, or the norm conservation is relaxed, the $`\stackrel{~}{𝒞}`$ terms must be included in the Hamiltonian.

Our result differs from that of ICL. The Hamiltonian derived in Ref. Ismail-Beigi et al. (2001) can be obtained from Eq. (13) if one replaces our dog-leg path $`𝐫^{}𝐑𝐫^{\prime \prime }`$ with a straight line path $`𝐫^{}𝐫^{\prime \prime }`$. The ICL Hamiltonian is also gauge invariant. In the presence of a magnetic field the value of the integral depends on the path, since $`\times 𝐀(𝐫)\mathrm{𝟎}`$. The two gauge invariant Hamiltonians are therefore different.

To clarify the situation, we examine the consequences of these differences for the calculation of the magnetic susceptibility. To compare the results obtained within the two methods with all-electron calculations, we restrict ourself to molecular systems for which it is possible to compute the susceptibility with quantum chemical approaches. The macroscopic magnetic susceptibility tensor $`\stackrel{}{\chi }`$ is defined from the second derivative of the system energy with respect to the external uniform magnetic field $`𝐁`$:

$$𝐁\stackrel{}{\chi }𝐁=2E^{(2)}$$
(14)
where $`E^{(2)}`$ is the second order variation of the energy with respect to the magnetic field:

$`E^{(2)}`$ $`=`$ $`2{\displaystyle \underset{o}{}}[\stackrel{~}{\mathrm{\Psi }}_o^{(0)}|\stackrel{~}{H}^{(1)}𝒢(ϵ_o)\stackrel{~}{H}^{(1)}|\stackrel{~}{\mathrm{\Psi }}_o^{(0)}`$ (15)
$`+\stackrel{~}{\mathrm{\Psi }}_o^{(0)}|\stackrel{~}{H}^{(2)}|\stackrel{~}{\mathrm{\Psi }}_o^{(0)}],`$
$`|\stackrel{~}{\mathrm{\Psi }}_i^{(0)}`$ and $`ϵ_i`$ are the unperturbed eigenstates and eigenvalues, $`𝒢(ϵ_o)=_e|\stackrel{~}{\mathrm{\Psi }}_e^{(0)}\stackrel{~}{\mathrm{\Psi }}_e^{(0)}|/(ϵ_oϵ_e)`$ and the $`o`$ and $`e`$ sums run over occupied and empty orbitals.

In a uniform magnetic field, with the gauge $`𝐀(𝐫)=𝐁\times 𝐫/2`$, our GIPAW Hamiltonian, Eq. (13), gives rise to the following perturbation Hamiltonians:

$$\stackrel{~}{H}_{\mathrm{GIPAW}}^{(1)}=\frac{1}{2c}\left(𝐋+\underset{𝐑}{}𝐑\times 𝐯_𝐑^{\mathrm{nl}}\right)𝐁,$$
(16)
and

$$\stackrel{~}{H}_{\mathrm{GIPAW}}^{(2)}=\frac{1}{8c^2}\left[(𝐁\times 𝐫)^2+\underset{𝐑}{}𝐁\times 𝐑\stackrel{}{𝐃}_𝐑^{\mathrm{nl}}𝐁\times 𝐑\right],$$
(17)
where $`𝐯_𝐑^{\mathrm{nl}}=[𝐫,V_𝐑^{\mathrm{nl}}]/i`$, $`D_{𝐑,\alpha ,\beta }^{\mathrm{nl}}=[r_\alpha ,[r_\beta ,V_𝐑^{\mathrm{nl}}]]`$, and $`\alpha `$ and $`\beta `$ are Cartesian indexes.

The corresponding perturbation Hamiltonians obtained following the ICL approach are

$$\stackrel{~}{H}_{\mathrm{ICL}}^{(1)}=\frac{1}{2c}\left(𝐫\times 𝐯\right)𝐁,$$
(18)
$$\stackrel{~}{H}_{\mathrm{ICL}}^{(2)}=\frac{1}{8c^2}\left[(𝐁\times 𝐫)^2+\underset{𝐑}{}𝐁\times 𝐫\stackrel{}{𝐃}_𝐑^{\mathrm{nl}}𝐁\times 𝐫\right],$$
(19)
where $`𝐯=𝐩+_𝐑𝐯_𝐑^{\mathrm{nl}}`$.

We compute $`\stackrel{}{\chi }`$ in molecules with both the GIPAW and the ICL approaches. We describe the electronic structure with density functional theory in the local density approximation. We use a large-cubic-periodic supercell of 6000 au3, in order to avoid the interaction of the molecules with their periodic replica, and Troullier-Martins pseudopotentials Troullier and Martins (1991) in the Kleinman-Bylander form Kleinman and Bylander (1982). We expand the wavefunctions in plane-waves with a cutoff of 100 Ry. The position operator is not defined within periodic boundary conditions. We treat it approximately by constructing a periodic saw-tooth like function centered on the molecules. For large cells this operator well approximates the position operator where the electron density is not negligible Pickard and Mauri (2001). The contribution of core electrons to magnetic properties can not be neglected. This contribution is however rigid, i.e. independent of the chemical environment Gregor et al. (1999); Pickard and Mauri (2001). We compute the core contribution with an atomic code.

Both the GIPAW and the ICL approaches are, by construction, gauge invariant. To verify that our numerical implementation and the use of a finite basis set preserve this property, we compute $`\chi `$ for a CH4 molecule as a function of the distance between the gauge origin and the molecular center. The results summarized in Table 1 show that the calculated $`\chi `$ is indeed gauge invariant, while the individual terms due to $`\stackrel{~}{H}^{(1)}`$ and $`\stackrel{~}{H}^{(2)}`$ clearly are not.

In order to obtain an independent assessment of the accuracy of the two methods, we compare to a truly all electron method, the individual gauges for atoms in molecules (IGAIM) method Keith and Bader (1992) as implemented in Gaussian98 Frisch et al. (2001). The magnetic susceptibility converges only slowly with Gaussian or atomic basis sets. The aug-cc-p(C)VxZ basis set series Peterson and T. H. Dunning (2002) has been previously shown exhibit reliable convergence for magnetic response properties Gregor et al. (1999); Helgaker et al. (1999), and we confirm this by converging the magnetic susceptibility for CH4 by using up to the aug-cc-pCV5Z basis for C (and aug-cc-pV5Z for H). For the remaining calculations we use the corresponding quadruple zeta basis sets, at which level the CH4 result is converged to better than 0.1 cm$`{}_{}{}^{3}/`$mole. Indeed, calculations using the largest basis sets rapidly become intractable for even moderately sized molecules. The results are summarized in Table 2, and show that the GIPAW method results in values for the magnetic susceptibility that are consistently closer to the all electron results than those calculated using the ICL method by roughly an order of magnitude. E.g., the ICL results for P2 deviate by almost one third of the total, while the maximum fractional deviation for our method is less than two percent.

We have derived, and demonstrated the practical utility of, a theory for the coupling of nonlocal pseudopotentials to arbitrary electromagnetic fields. While the ICL method may be the best possible for an arbitrary nonlocal Hamiltonian we have shown that by focusing on the nonlocal Hamiltonian most frequently encountered in electronic structure calculations we are able to make considerable improvements. This is of great importance if the results are to be compared quantitatively with experiment.

The work of CJP was supported by the Universités Paris 6 and 7 in France, and the EPSRC in the United Kingdom. We thank the High Performance Computing Facility of the University of Cambridge for access to Hodgkin, their SGI Origin 2000, and the CNRS funded IDRIS supercomputing center. We would also like to thank Thomas Gregor for the molecular geometries, and Kirk Peterson for his kind help with the Gaussian basis functions.