format: text

text:
# A numerical canonical transformation approach to quantum many body problems

## I Introduction

The vast majority of current methods for numerically solving quantum many body systems fall into a few broad categories. One type of approach involves producing approximate representations of wavefunctions. Examples of this approach include configuration interaction, coupled cluster methodscc , and the density matrix renormalization group (DMRG)dmrg ; WhiteMartin ; mitrushenkov . A second type of approach is based on perturbation theory. Some versions of perturbation theory are closely related to wavefunction approaches, while others, for example utilizing finite temperature imaginary time Green’s functions, are more closely related to path integrals. A third type of approach is quantum Monte Carlo, which also may involve representations of wavefunctions (e.g. Green’s function Monte Carlo) or of path integrals (e.g. determinantal/auxiliary field methods).

Much less explored than these approaches are methods based on similarity or unitary transformations of the Hamiltonian $`H`$. In these approaches the primary focus is on $`H`$ and transformed versions of it; wavefunctions play a much more minor role. In this paper we develop such an approach in the context of ab initio quantum chemistry calculations in a finite basis. This approach is based on unitary canonical transformations (CTs) of $`H`$ written in second-quantized operator form. Such transformations have been used in analytical work for very long timevanvleck . A well known example in condensed matter physics is the Schrieffer-Wolff transformation of the Anderson model into the Kondo modelschrieffer ; another is the well known mapping of the Hubbard model in the large $`U/t`$ limit into the $`t`$-$`J`$ model. Often these transformations are performed once, and relate one model system to another, simpler system, with fewer degrees of freedom, in an approximate way. In some cases special CTs can produce exact solutions for certain model systems.

Recently, substantial progress has been made by the development of continuous unitary transformations, in which a set of differential equations is solved to perform the CT.wegner ; glazek  This method, which was developed independently by Wegner and by Glazek and Wilson, is known by the names “flow equation method” wegner  and “similarity renormalization” glazek . A key advantage of this approach is that one does not need to know in advance the transformation operator to be used; it is determined implicitly by the solution of the differential equations. Another advantage is that once the differential equations are set up, there is no operator algebra to be performed in the course of the numerical solution of the differential equations.

This approach can be performed in a semianalytic context, where typically one is deriving one model from another. For example, an improved treatment of the Schrieffer-Wolff transformation of the Anderson model has been performed. kehrein  It can also be used to obtain ground state energies and dispersion relations (i.e. excitation energies) for infinite lattice systems. knetter

Here we extend and develop the CT approach in a quantum chemical context. We first introduce a new version of the CT approach which is closely related to the Jacobi method of diagonalizing matrices. In this case, rather than solving a differential equation, one performs a sequence CTs, each involving the smallest possible number of operators. This approach is designed to solve a finite quantum system, namely a molecule or cluster in a standard quantum chemical basis, as an alternative to other standard approaches, such as configuration interaction or coupled cluster methods. We demonstrate that our method can determine ground states in ab initio chemical systems with excellent precision, utilizing tests on a water molecule for which exact results are available. We expect similar performance for low lying excited states. We then present a version of the “flow equation” approach for use in the ab initio context. Utilizing ideas developed in our Jacobi method, we present a new version of the differential equations which make their numerical solution particularly efficient. We demonstrate that this method also works very well for the water molecule.

Perhaps the most important aspect of our approaches is the ability to remove both low energy (i.e. core) and high energy virtual orbitals from the problem, leaving a system with a small number of “active” orbitals. This approach is especially useful for “strongly correlated” systems, i.e. those with open shells, breaking bonds, etc., where a single reference approach fails. In these systems, the strong correlation is generally confined to a relatively small number of orbitals. In wavefunction-based approaches, it can be awkward and expensive to deal with the strongly correlated part of the problem with a powerful method (e.g. full diagonalization of the active space) and the simple high-energy part of the problem with another, simpler method. By working with the Hamiltonian directly, one can separate the two parts of the problem simply and efficiently. First one transforms the whole Hamiltonian into a form in which the high and low energy orbitals are either completely unoccupied or complete occupied (for core orbitals). Then these orbitals can be thrown out of the problem completely, leaving a smaller system of partially occupied orbitals which can then be solved with any of a variety of nonperturbative techniques. Here we demonstrate this powerful hybrid approach on a stretched water molecule, solving the smaller systems with DMRG.

Another advantage in working with the Hamiltonian directly is in obtaining excitations. In transforming the Hamiltonian, we make progress in solving for excited states at the same time as we obtain the ground state. Essentially the same approach is used to obtain both ground and excited states. Although we have not performed tests yet for excited states, CT approaches in other contexts have obtained excellent results. knetter

We will call the set of methods we present here canonical diagonalization (CD). The term diagonalization indicates our intent to solve the system fully, rather than just transforming to a simpler model. The transformations involved are both unitary and canonical (which is defined in the next section); we choose only the term canonical to emphasize that the method works in the space of second quantized operators. We further distinguish the Jacobi CD method (JCD) and the flow equation CD method (FECD). All CD approaches share the feature that the object that one is manipulating is the second quantized Hamiltonian, as a collection of abstract operator terms with specific numerical coefficients.

These approaches seem particularly suited to ab initio quantum chemistry, which are characterized by a very general quantum Hamiltonian, containing almost all possible one and two-electron terms. The CTs generate additional terms involving one, two, three, and more particles. Since general one and two particle terms are already present in the Hamiltonian, no extra inconvenience arises from these terms. Three and more particle terms are more inconvenient, but most such terms can be neglected, to an excellent approximation. CD is size-consistent, and many-particle terms which may be left out involve the simultaneous interaction of three or more (dressed) electrons, so that neglecting them is analogous to neglecting connected clusters involving triple and higher excitations in coupled cluster methods. Note that canonical transformations also appear in the theory of the coupled cluster methodnooijen , although the method remains largely a wavefunction approach. Note also that although one does not need to write any wavefunctions explicitly, CD in its simpler forms implicitly expresses the ground state using the exponential of an operator acting on a reference state. Further links to coupled cluster methods are made in the Discussion Section.

CD fits naturally into a renormalization group (RG) framework. First, one can remove (“integrate out”) higher energy orbitals, one at a time if one wishes, leaving a system where the effects of the removed orbitals are incorporated into an effective Hamiltonian for the remaining orbitals. Thus each step resembles a transformation in a typical RG calculation in statistical mechanics, although unlike in that case one cannot continue indefinitely and there are no fixed points. Second, even if one is not integrating out orbitals, the transformation of the Hamiltonian, like in RG methods, occurs in a sequence of steps, with truncation of higher order terms occuring at each step. Third, the differential “flow” equation form of CD, in which a time-like variable controls the evolution of the Hamiltonian operator towards a more diagonal form, closely matches Wilson’s original conception of the RG approachwilsonrg .

CD is a natural complement to DMRG, and this was a principle motivation in developing it. When applied to quantum chemistry problems, DMRG does very well in describing non-dynamical correlations (strong correlations associated with partially occupied orbitals), but it is inefficient in describing dynamical correlations, since it describes high energy virtual orbitals on the same footing as partially occupied strongly interacting orbitalschan . CD has complementary behavior. It can be used to remove the nearly unoccupied orbitals, leaving a smaller Hamiltonian involving strongly interacting orbitals for DMRG to solve.

## II Jacobi CD

We begin with the Jacobi CD approach. In the ordinary Jacobi method for diagonalizing matrices, one applies a large number of unitary transformations to a Hermitian matrix to bring it into diagonal formnumrecipe . A unitary transformation gives a new matrix which has the same eigenvalues as the old. In the Jacobi method, each unitary transformations consists of a rotation of two rows and columns to zero out a single off-diagonal element $`H_{ij}`$. The part of the unitary transformation matrix $`\mathrm{exp}(A)`$ corresponding to rows and columns $`i`$, $`j`$ is

$$\mathrm{exp}\left(\begin{array}{ccc}0& \theta & \\ \theta & 0& \end{array}\right)=\left(\begin{array}{ccc}\mathrm{cos}\theta & \mathrm{sin}\theta & \\ \mathrm{sin}\theta & \mathrm{cos}\theta & \end{array}\right)$$
(1)
The angle $`\theta `$ which removes the term $`H_{ij}`$ is given by

$$\theta =\frac{1}{2}\mathrm{tan}^1[2H_{ij}/(E_iE_j)]$$
(2)
where $`E_i=H_{ii}`$ and the transformation is applied as

$$\mathrm{exp}(A)H\mathrm{exp}(A).$$
(3)
Any unitary transformation can be written as an exponential of an antiHermitian matrix $`A`$wagner . For real symmetric matrices and operators $`H`$ and for what follows, it suffices to use real antisymmetric matrices and operators $`A`$. In the Jacobi method, one traverses the matrix repeatedly, rotating away off-diagonal elements, starting with the largest off-diagonal elements for efficiency.

In Jacobi CD we construct unitary transformations to successively remove off-diagonal terms of a second quantized Hamiltonian. We consider a quantum chemical system in a Hartree Fock basis, with Hamiltonian

$$H=\underset{ij\sigma }{}T_{ij}c_{i\sigma }^{}c_{j\sigma }+\frac{1}{2}\underset{ijkl\sigma \sigma ^{}}{}V_{ijkl}c_{i\sigma }^{}c_{j\sigma ^{}}^{}c_{k\sigma ^{}}c_{l\sigma }.$$
(4)
Here $`T_{ij}`$ contains the electron kinetic energy and the Coulomb interaction between the electrons and the nuclei, while $`V_{ijkl}`$ describes the electron-electron Coulomb interaction. Indices such as $`i`$ denote spatial Hartree-Fock orbitals, and $`\sigma `$ is a spin index. Later on we shall sometimes use $`i`$ to denote a spin-orbital, in which case the context should make the usage self-evident. We use a computer representation of $`H`$ as a sum of abstract operator terms, each with a coefficient and a product of $`c`$ and $`c^{}`$ operators involving various orbitals. In our program a $`c`$ or $`c^{}`$ operator is described by a single byte, which encodes the orbital index involved, the spin, and whether it is $`c`$ or $`c^{}`$. This implementation is thus limited to systems with at most 64 orbitals. A complete operator term is stored as an array of such bytes plus a floating point coefficient. We developed a set of C++ routines to take products and commutators of such operators, putting the result in normal ordered form using the anticommutation relations. Having these operations be reasonably efficient is crucial to the method. This approach, using formal operator terms to describe $`H`$, rather than specific matrix expressions for terms of various orders, is both simple and general. In the future, in order to implement specific approximations within CD, more efficient code could be produced by deriving and implementing the relevant matrix expressions.

An off-diagonal term which can be rotated away is simply any term which is distinct from its Hermitian conjugate. Self-adjoint terms constitute the diagonal elements. They can always be written as products of density operators $`n_{i\sigma }=c_{i\sigma }^{}c_{i\sigma }`$. Consider, as a specific example, the term

$$V_\alpha =a\overline{V}_\alpha =ac_i^{}c_j^{}c_kc_l,$$
(5)
where $`a`$ is a numerical coefficient. Let

$$A=\theta (\overline{V}_\alpha \overline{V}_\alpha ^{})$$
(6)
Given the proper choice of $`\theta `$, the transformation $`H\mathrm{exp}(A)H\mathrm{exp}(A)`$ will remove $`V_\alpha +V_\alpha ^{}`$ from $`H`$, introducing other terms instead. We will consider the choice of $`\theta `$ momentarily. These additional terms in general have smaller coefficients than $`V_\alpha `$, making $`H`$ more diagonal. If one continued the process indefinitely, one would have a “classical” Hamiltonian where every term was diagonal. In this case any Slater determinant $`c_i^{}c_j^{}\mathrm{}|0`$ is an eigenstate and all eigenvalues can be read off essentially by inspection. After the particle-hole transformation of the occupied orbitals described below, normally the ground state energy is simply the constant term in $`H`$.

This unitary transformation is a canonical transformation, in that when applied to the operators $`c`$ and $`c^{}`$, the anticommutation relations are preserved, e.g.

$$\{e^Ac_ie^A,e^Ac_j^{}e^A\}=e^A\{c_i,c_j^{}\}e^A=\{c_i,c_j^{}\}.$$
(7)
One can view a CT as a complicated many-particle change of basis. In that sense CD is similar to the DMRG method, the biggest difference being that DMRG works in a wavefunction basis. In the Discussion section, we comment more on applying the CTs to the $`c`$ operators.

In order to carry out the transformations, it is convenient to use the well-known formula

$`e^AHe^A`$ $`=`$ $`H+[A,H]+{\displaystyle \frac{1}{2!}}[A,[A,H]]`$ (8)
$`+`$ $`{\displaystyle \frac{1}{3!}}[A,[A,[A,H]]]+\mathrm{}`$
A commutator between an $`N_1`$-particle term and an $`N_2`$-particle term gives up to $`(N_1+N_21)`$-particle terms. Hence, if $`A`$ is a one particle term, it does not generate any higher order terms. Indeed, rotations by one-particle $`A`$’s correspond to single particle changes of basis. CD using one-particle $`A`$’s can be used to perform a Hartee-Fock calculation, if the initial orbitals are not the HF orbitals. To evaluate Eq. (8), we generally treat each term in $`H`$ separately. For most terms $`V_\beta `$ the transformation has no effect, since $`[A,V_\beta ]=0`$. For relevant terms the expansion can be carried out order by order until all terms of a given order are neglible (i.e. their coefficients are very small). At each order the terms should be put in normal ordered form. The number of distinct terms generated from a single term $`V_\beta `$ is a modest finite number, since a particular $`c_i`$ or $`c_i^{}`$ can appear at most once in a term in its normal ordered form.

For several reasons it is convenient to perform a particle hole transformation on the occupied orbitals. Thus we define $`d_{i\sigma }=c_{i\sigma }`$ if orbital $`i`$ is unoccupied in the HF state, and $`d_{i\sigma }=c_{i\sigma }^{}`$ if it is occupied. The $`d`$’s have identical anticommutation relations to the $`c`$’s. After this transformation, anticommutation relations are used to put the terms in normal ordered form, putting $`d`$ to the right of $`d^{}`$. The resulting anticommutators generate new lower order terms, after which the HF energy appears as a constant term in $`H`$. In terms of the $`d`$ operators, the HF state is the vacuum $`|0`$. There are no off-diagonal single particle terms such as $`d_i^{}d_j`$. There are terms which appear to violate particle conservation, such as $`d_i^{}d_j^{}d_k^{}d_l^{}`$, but which in fact do not in terms of real particles. As one performs the canonical transformations, the vacuum state approaches the exact ground state.

We now consider the choice of $`\theta `$. An operator term $`V_\alpha `$ connects an exponentially large number of pairs of states $`l`$,$`r`$ together, $`l|V_\alpha |r0`$. However, one of these pairs of states can be considered the most important, namely, the one closest to the HF vacuum. This pair has the fewest number of $`d^{}`$’s operating on $`|0`$ which generate states not destroyed by $`V_\alpha `$. As a specific example, let $`V_\alpha =0.1d_i^{}d_jd_kd_m`$. Then the most important pair of states is

$`|r`$ $`=`$ $`d_j^{}d_k^{}d_m^{}|0`$
$`|l`$ $`=`$ $`d_i^{}|0.`$ (9)
We will call these states simply the left and right states of $`V_\alpha `$. Other pairs of states, considered to be less important, have additional $`d^{}`$’s, which do not appear in $`V_\alpha `$, applied to both $`|l`$ and $`|r`$. For example, one could take the pair $`d_n^{}|l`$ and $`d_n^{}|r`$. Define for a term $`V_\alpha `$ the left energy $`E_l=l|H|l`$ and similarly the right energy $`E_r=r|H|r`$. We now choose to use the Jacobi formula for $`\theta `$ to attempt to eliminate the off-diagonal term in the Hamiltonian connecting $`l`$ and $`r`$. Thus, if $`a`$ is the coefficient of $`V_\alpha `$, we choose

$$\theta =\frac{1}{2}\mathrm{tan}^1[2a/(E_lE_r)].$$
(10)
This does not eliminate $`V_\alpha `$ exactly, since the operator $`A`$ connects many different states, not just $`l`$ and $`r`$, so it is not a $`2\times 2`$ transformation, as it would be for a matrix. Nevertheless, we find that this choice generally works well, typically reducing the size of the coefficient of $`V_\alpha `$ by a few orders of magnitude. Note that the degenerate case $`E_l=E_r`$ is nonsingular, generating an angle of $`\pm \pi /4`$ (either angle can be chosen). Such a large transformation angle should be avoided if possible, however, since it generates high order terms in the transformation of $`H`$.

A more common choice in analytic work using CTs is to choose to eliminate $`V_\alpha `$ to first order in the expansion Eq. (8), namely choosing $`\theta `$ to set the coefficient of $`\overline{V}_\alpha `$ in $`[A,H^D]`$ to $`a`$, where $`H^D`$ is the diagonal part of the Hamiltonian. This is closely related to our approach: note that $`l|[\overline{V}_\alpha ,H^D]|r`$ is the coefficient of $`\overline{V}_\alpha `$ in $`[\overline{V}_\alpha ,H^D]`$. However,

$$l|[\overline{V}_\alpha ,H^D]|r=(E_rE_l)l|\overline{V}_\alpha |r=E_rE_l$$
(11)
so that this choice gives $`\theta (E_rE_l)=a`$. This agrees with our choice to lowest order, but it is not well behaved if $`E_rE_l`$.

One need not eliminate all off-diagonal terms in $`H`$. If one is interested in only the ground state, then one needs to eliminate all terms connecting that state to other states. More specifically, suppose the initial HF state $`|0`$ has substantial overlap with the ground state. Only states which produce a nonzero result when acting on the vacuum need be removed, namely, only terms such as $`d_j^{}d_k^{}d_m^{}d_n^{}`$ or $`d_j^{}d_k^{}`$, and their Hermitian conjugates (as well as similar multiparticle terms). Since all terms still satisfy particle conservation, in each of the two-particle terms two of the orbital indices $`jkmn`$ must correspond to occupied states, and two to unoccupied orbitals. Once all such terms are eliminated, then the state $`|0`$ is the ground state, and the ground state energy is the constant term in the Hamiltonian.

In all but the smallest systems, some of the terms formed from the CTs must be discarded according to some criterion. The simplest criterion is to neglect all terms involving three or more particles, i.e. six or more $`d`$ operators. Our test calculations on the water molecule suggest that this is a very accurate approach for systems well described by a single reference state. Other possible criteria include keeping all terms whose coefficients are larger than some cutoff; keeping all one and two particle terms and all three particle terms larger than a cutoff, etc. More sophisticated criteria are possible also, such as trying to estimate the contribution of each term using pertubation theory, and discarding terms whose contribution is below a cutoff. Here, we perform some test calculations according to simple cutoff criteria. In the future, we hope the criteria can be optimized.

In order to preserve symmetries, such as spin symmetry, one can rotate sets of terms which are related by a symmetry transformation in one step. For example, in what follows, for each term, we check to see if it is distinct from the term coming from flipping all of its spin indices. If it is distinct, both are rotated together with the same rotation angle. The rotation angle for both is chosen as the angle to rotate one of the terms separately. This procedure preserves spin symmetry exactly.

One must decide in which order to go through the terms in performing the CTs. Since each CT alters the coefficients of many other other terms in H, it makes sense to start with the largest first. One approach would be to find the term with the largest magnitude coefficient at each step. Another would be to choose the largest rotation angle. However, searching for the largest term at each step would be inefficient. Therefore, we have chosen the following method: a cutoff angle is chosen, and all terms with angles greater in magnitude than this cutoff are treated in a sweep through the terms, in a predetermined but arbitrary order. Then, the cutoff angle is reduced by a constant factor, and the procedure is repeated. Here, we started with a cutoff of 0.15 and reduced it using a factor of 0.6. (In some passes, particularly the initial one, there may be no CTs performed.) In Table I we show results for a 25 orbital DZP basis water molecule, for which full CI results are availablebauschlicher . Because of some arbitrary choices in the ordering of the CTs, which unfortunately can affect the results slightly, the results here would be difficult to reproduce precisely by an independently written program. (The differential equation method discussed below does not suffer this problem.) Despite this, one can easily evaluate the potential of the method from our results. One can see that we obtain an accuracy of several millihartrees even when truncating all three or more particle terms. If one keeps also some three particle terms, one can obtain accuracy to fractions of a millihartree. This sort of accuracy is comparable to coupled cluster methods.

It is not necessary to perform every CT in this procedure. At the end of every sweep, we perform a calculation of the energy using second order perturbation theory for the current Hamiltonian. The calculation time for this procedure scales only as the number of terms in H, so there is a neglible impact on the overall computation time. As the largest angle terms are eliminated, the perturbation theory result becomes more and more accurate. Even just a few rotations of the largest terms can make perturbation theory much more well-behaved. One can stop the procedure when the perturbation result is well converged, which typically happens long before all the chosen terms are removed. In Fig. 1 we show the results for this procedure. One can see that the perturbation result converges much more quickly than the constant term in the energy. One might well stop after about 300 Jacobi steps; in this particular example this number is of order $`N^2`$, where $`N`$ is the number of orbitals.

In Table II we show similar results for a water molecule whose bonds have been stretched by a factor of two. This system is not well described by a single reference state: in the full CI calculations of Olsen, et. al.olsen , in a different but similar basis, the weight of the HF determinant in the full CI wavefunction was 0.589, versus 0.941 for the unstretched molecule. Here, we find that CD is unstable if only two-particle terms are kept. One finds that repeated Jacobi diagonalization steps reduce the energy without bound. CD is exact if no truncations are made, so this is an artifact of the truncation of three and more particle terms. However, keeping even a large number of three particle terms does not result in a particularly accurate calculation.

For this system, examination of the occupancies of the HF orbitals in the exact ground state (which we have computed with high accuracy with DMRG) reveals that there are four spatial orbitals with occupancies far from 0 or 2; specifically, they have occupancies of 1.58, 1.52, 0.46, and 0.4. The rest have occupancies less than 0.03 or more than 1.97. In the case of the unstretched water molecule, occupancies are all within 0.05 of 0 or 2. The results for occupancies of natural orbitals are very similarolsen . The contribution to the energy of a Hamiltonian term $`A`$, $`\psi |A|\psi `$, can be expressed as a Green’s function or density matrix element. In the case of Hamiltonian terms made out of operators involving only nearly filled or unfilled orbitals, the behavior of the Green’s function is well understood, and the magnitude falls rapidly as one considers terms involving more particles. For partially occupied orbitals there is no reason to believe that three or more particle Green’s function elements are small. Consequently, one should only truncate such a term if its coefficient is small. For this reason, we have performed test calculations with the following truncation criterion: all terms with more than four $`d`$ and $`d^{}`$ operators corresponding to non-partially-filled orbitals are truncated. In addition, all terms whose coefficient’s magnitude is below a cutoff $`\epsilon `$ are eliminated. If there are $`N_p`$ partially filled orbitals, then this rule allows terms with up to $`4N_p+4`$ $`d`$’s to appear. In this case we have up to 20 $`d`$’s, i.e. a 10-particle term.

As shown in Table II, with this criterion we see substantially better results: we find that in this non-single reference system, accuracy to fractions of a millihartree is possible.

## III Flow Equation CD

It is also possible to formulate CD in terms of a differential equation. This approach was originally developed independently by Wegnerwegner  and by Glazek and Wilson,glazek  in rather different contexts than we present it here. We will derive it here as a natural variation of the Jacobi CD method. In this approach we introduce a time-like variable $`t`$, and the Hamiltonian evolves as $`t`$ increases. First consider a time-dependent Hamiltonian for some fixed antihermitian operator A:

$`H(t)=e^{tA}H(0)e^{tA}`$ (12)
Here $`H(0)`$ is the initial HF Hamiltonian. We have

$`{\displaystyle \frac{dH(t)}{dt}}=[A,H(t)].`$ (13)
This differential equation form of a CT has long been used in analytical work, where one integrates $`t`$ from 0 to some fixed rotation angle. Here we modify this by making $`A`$ depend on $`H`$. First expand $`H`$ as follows

$`H(t)={\displaystyle \underset{\alpha }{}}a_\alpha (t)h_\alpha .`$ (14)
Each $`h_\alpha `$ is a product of $`d`$ and $`d^{}`$ operators, and $`a_\alpha (t)`$ is the corresponding coefficient. Let

$`A(t)={\displaystyle \underset{\alpha }{}}s_\alpha a_\alpha (t)h_\alpha .`$ (15)
The $`s_\alpha `$ are fixed parameters, which we initially consider to have only three possible values: $`\pm 1`$, and 0. We set $`s_\alpha `$ to 0 if we are not interested in rotating the coefficient of $`h_\alpha `$ to 0, because, for example, $`h_\alpha `$ does not act directly on the HF state $`|0`$. For terms we wish to rotate to zero, we choose the sign of $`s_\alpha `$ so that (1) $`A(t)`$ is antihermitian, and (2) increasing $`t`$ rotates in the direction to diminish $`a_\alpha (t)`$. These conditions are satisfied if $`s_\alpha `$ is chosen as the sign of $`E_lE_r`$, where $`l`$ and $`r`$ are the left and right states of $`h_\alpha `$. We evolve $`H(t)`$ as a sequence of infinitesimal CTs, as follows

$`H(t+\delta t)`$ $`=`$ $`e^{\delta tA(t)}H(t)e^{\delta tA(t)}`$ (16)
$`=`$ $`H(t)+\delta t[A(t),H(t)]+O(\delta t)^2.`$
In the limit that $`\delta t0`$, this is equivalent to solving the nonlinear differential equation

$`{\displaystyle \frac{dH(t)}{dt}}=[A(t),H(t)].`$ (17)
Each infinitesimal rotation acts to diminish each $`a_\alpha (t)`$ with nonzero $`s_\alpha `$. Since $`A(t)`$ depends linearly on the $`a_\alpha (t)`$, the rotations become smaller as the $`a_\alpha (t)`$ decrease. Thus, we expect the solution of this equation for $`t\mathrm{}`$ to have $`a_\alpha (t)=0`$ if $`s_\alpha `$ is nonzero. We also expect these $`a_\alpha (t)`$ to diminish exponentially with $`t`$.

If no truncations are made, the solution to this differential equation for any time $`t`$ gives an $`H(t)`$ which is related to $`H`$ by an exact CT. This is true also for any choice of the $`s_\alpha `$ as long as they satisfy the requirement that if $`h_\beta =h_\alpha ^{}`$, then $`s_\beta =s_\alpha `$, ensuring that $`A(t)`$ is antihermitian. For numerical efficiency, it is useful to modify the choice of $`s_\alpha `$. This is because different terms $`h_\alpha `$ require different rotation angles. One would like to make the exponential decay to zero of each $`a_\alpha (t)`$ have approximately the same time constant. If they have widely varying time constants, the number of steps in integrating the differential equation will be very large. We can achieve this by choosing, for nonzero $`s_\alpha `$,

$`s_\alpha =(E_lE_r)^1.`$ (18)
Provided $`a_\alpha (t)E_lE_r`$, this choice makes the coefficient of $`h_\alpha `$ in $`A(t)`$ the angle $`\theta `$ required to rotate the term to zero. This makes the natural time scale for each term equal to unity. We choose the $`s_\alpha `$ at the beginning, using the untransformed HF energies, and never change them; however, one could also make the $`s_\alpha `$ depend on $`t`$.

In Wegner’s original flow equation method, rather than the above forms of $`A`$ defined in terms of $`s_\alpha `$, one took $`A=[H^D,H]`$, where $`H^D`$ is the diagonal part of $`H`$. This is very similar to the choice $`s_\alpha =E_lE_r`$, assuming all off-diagonal terms are being removed. However, this choice gives very widely varying time scales, driving terms with large $`E_lE_r`$ to zero much more quickly. In the sense that the large energy difference terms are removed first, Wegner’s method can be considered a renormalization group method in itself, and one might stop at some finite time and study the partially transformed Hamiltonian. Our choice is much more efficient numerically, assuming one only wants the $`t\mathrm{}`$ limit.

We decribe all the commutator relations in terms of a “matrix” $`B`$

$`[h_\alpha ,h_\beta ]={\displaystyle \underset{\gamma }{}}B_{\alpha \beta }^\gamma h_\gamma .`$ (19)
If a commutator gives a term which is not in the set of Hamiltonian terms we are keeping, then that term is ignored. Then the final form for the flow equation CD method is a set of differential equations

$`{\displaystyle \frac{da_\gamma (t)}{dt}}={\displaystyle \underset{\alpha \beta }{}}B_{\alpha \beta }^\gamma s_\alpha a_\alpha (t)a_\beta (t)`$ (20)
which are to be solved numerically. The $`B`$ matrix was computed initially and stored in our program. Because of some regularities in the pattern of nonzero elements of $`B`$, the storage could be reduce by a factor of about $`N`$, the number of orbitals, from a naive estimate. However, they could also be recomputed at each step to save storage, at the expense of computer time. Another approach to save storage would be to remove a few orbitals at a time. One could even remove one term at a time by making only one $`s_\alpha `$ nonzero, in which case the flow equation method becomes very similar to the Jacobi method. To integrate the coupled differential equations, we use a simple fourth order Runge Kutta method with automatic step size adjustment. This routine attempted to integrate the differential equations with an absolute error tolerance of $`10^8`$, and we integrated the equations from $`t=0`$ to $`t=20`$.

In Figure 2, we show the evolution of the constant term in $`H`$ as a function of $`t`$ for the unstretched water molecule. Only one and two particle terms were retained. The step sizes used were rather large, and they steadily increased. They are visible via the circles in the curve. Only twelve steps were taken, although each RG step in our very crude integrator required twelve derivative evaluations, Eq. (20). The result for the energy was in error only by about a milli-hartree.

In Figure 3, we show similar results for the stretched water molecule. As in the Jacobi method, CD keeping only two particle terms is unstable, with the energy tending to $`\mathrm{}`$. We believe that by keeping multiparticle terms one could make this method perform very well on the stretched water molecule, just as we found for the Jacobi method.

## IV Removing sets of orbitals

Another approach for systems such as the stretched water molecule, which have some strongly correlated orbitals, is to integrate out many of the non-strongly-correlated-orbitals, leaving a small but strongly correlated system to solve with CD retaining many-particle terms, with DMRG, or with another method. We first divide the orbitals into two sets, those to be kept and those to be removed. Some of the orbitals to removed will have occupancies near 0, and some may be core orbitals with occupancies near 2. Due to the particle-hole transformation, we need make no distinction between these cases. Consider the many-particle basis states $`|s=d_i^{}d_j^{}\mathrm{}|0`$. Let $`|s`$ denote all states in which no orbital to be removed is occupied. Conversely, let $`|s^{}`$ denote the rest, in which at least one orbital to be removed is occupied. We wish to rotate away all Hamiltonian terms which connect states $`|s`$ to $`|s^{}`$. Let $`r`$ represent orbitals to be removed. Then the terms to be removed are described by the following rule: the terms have one or more $`d_r^{}`$’s, or one or more $`d_r`$’s, but do not have both $`d_r^{}`$’s and $`d_r`$’s.

One finds that typically a few of these terms to be removed, largely by accident, have $`E_r`$ and $`E_l`$ nearly identical, although neither is close to zero because they include operators adding or removing high-energy orbitals. To remove these problem terms requires a large angle of rotation. This can be disastrous for either the Jacobi or flow equation method unless many-particle terms are kept. However, some reflection indicates that these problem terms are likely to be quite unimportant in terms of their true contribution to the ground state. Consider an nth order perturbation theory contribution to the ground state energy. Ignoring the energy denominators, such a term is proportional to

$`0|h_1h_2\mathrm{}h_n|0.`$ (21)
Clearly, if this contribution is nonzero, term $`h_n`$ must have $`E_r=0`$, and term $`h_1`$ must have $`E_l=0`$. Thus our problem term with nearly degenerate nonzero energies cannot contribute in second order. For third order, one might consider $`h_2`$ to be the problem term. However, for this term either $`|r`$ or $`|l`$ must belong to the set $`|s^{}`$. Removing all the other terms connecting $`|s`$ to $`|s^{}`$ means that either $`h_1`$ or $`h_3`$ is to be removed, since $`|0`$ belongs to $`|s`$. Thus there is no third order contribution. The lowest order contribution for such a term is in fourth order, involving both this term and its Hermitian conjugate as $`h_2`$ and $`h_3`$, or involving two such terms as $`h_2`$ and $`h_3`$. Here $`h_4`$ takes one from $`|0`$ into a higher energy state $`|s`$, $`h_3`$ takes one from $`|s`$ to $`|s^{}`$, $`h_2`$ then takes one back to $`|s`$, and $`h_1`$ takes one back to $`|0`$. The energy denominators are well-behaved, since $`E_r`$ and $`E_l`$ of the problem term are not close to zero.

Thus a quite reasonable approach is to rotate away all terms connecting $`|s`$ to $`|s^{}`$ except those whose energy difference $`|E_rE_l|`$ is below some cutoff $`d`$. The terms below the cutoff are retained during the CD process, during which time they may change due to other terms being rotated away. After the CD is complete, one then can discard all terms having any $`d_r`$ or $`d_r^{}`$ operators, which will include these problem terms. One can also solve the CD-transformed $`H`$ before truncation, using another method, and check that the occupancies of the removed orbitals are very close to zero.

In Table IV we show the results of such calculations for the stretched water molecule. There are three sources of error in these calculations. First is the DMRG error, typically near $`0.0002`$ mH, keeping 400-600 states, which is small enough to show the other sources of error. Second is the error from performing CD keeping only one and two particle terms. This is given by $`\mathrm{\Delta }E_{\mathrm{CD}}`$. Increasing $`d`$, or removing fewer orbitals improves $`\mathrm{\Delta }E_{\mathrm{CD}}`$. Third, there is the energy from throwing away the removed orbitals after CD. This is measured by the difference between $`\mathrm{\Delta }E_{\mathrm{CD}}`$ and $`\mathrm{\Delta }E_{\mathrm{CDT}}`$, and also by the maximum occupancy of the removed orbitals $`n_{\mathrm{max}}`$. We find that $`d`$ can be made quite large: $`0.5`$ is always fine, whereas $`1.0`$ can be too large. We also find that we can remove up to about one half of the orbitals and incur only a very small error, even only keeping one and two particle terms. For the resulting small system even full CI would be a very easy calculation. Even removing all but four of the orbitals we get a reasonable result. We have not carried out any similar calculations keeping many-particle terms, but we can deduce the probable outcome. Since all the rotation angles $`\theta `$ are rather small in this procedure, four particle terms, which can come in only as $`\theta ^2`$, would be neglible. Three particle terms come in as $`\theta `$, and if such a term only involved the retained orbitals it presumably would have both $`E_l`$ and $`E_r`$ small and it could give a substantial contribution to the energy of order $`\theta `$. Three particle terms involving removed orbitals would have $`E_l`$ or $`E_r`$ reasonably large, and would only contribute to the energy via second order perturbation terms, thus coming in as $`\theta ^2`$, which could be neglected. In short, we expect that keeping three particle terms involving the retained orbitals only would be a very accurate approach for removing more than half of the orbitals.

We would like to conclude this section with an argument that the proper way to separate the treatment of high-energy from low-energy orbitals is by using an effective Hamiltonian to remove the high energy orbitals, as we have done, rather than any wavefunction based approach. We will make this argument via a trivial $`3\times 3`$ matrix, designed to have some of the crucial features of a strongly correlated/multireference system. Define the matrix

$`H(\epsilon ,\delta )=\left(\begin{array}{ccc}0& \epsilon & \delta \\ \epsilon & \epsilon & 1\\ \delta & 1& 10\end{array}\right).`$ (25)
The third row and column represent a high energy orbital, which we would like to treat separately from the first two nearly degenerate rows and columns. We will consider the parameter values ($`\epsilon =0.1`$, $`\delta =1`$), ($`\epsilon =0.1`$, $`\delta =0.5`$), and ($`\epsilon =0`$, $`\delta =1`$). For these three parameters we find the following ground state energies and eigenvectors (respectively): -0.099, and (0.995,-0.0098,-0.098); -0.064, and (0.789,-0.614,-0.022); and -0.196, and (0.700,0.700,-0.137). Now suppose we wanted to solve this system in two steps, first treating the third “orbital”, then next the other two, using a wavefunction approach. In treating the third orbital we insist that we ignore the small parameter $`\epsilon `$; otherwise we are treating the whole matrix together. We imagine that we have some perturbative method for obtaining the third component of the wavefunction, ignoring $`\epsilon `$; with this fixed, then we obtain the first two components, taking $`\epsilon `$ into account. However, comparing the first and third sets of parameters, we see that the third component $`\psi _3`$ depends strongly on $`\epsilon `$, so this method must fail.

Alternatively, we might imagine first treating the first two rows and columns separately, ignoring $`\delta `$ and finding the ratio of components $`\psi _1/\psi _2`$, and then subsequently using $`\delta `$ to fix $`\psi _3`$. In this case, comparing the first and second parameter sets, we see that $`\psi _1/\psi _2`$ depends strongly on $`\delta `$, so that this method fails. In short, to treat this problem successfully, wavefunction based approaches must treat both $`\delta `$ and $`\epsilon `$ simultaneously.

Now consider a simple CT approach. Rather than using the Jacobi or flow equation method, we use a less sophisticated, but well-known perturbative CT method for removing the third row and column.lax  In this case, we find that the second-order change in the upper left $`2\times 2`$ portion of the matrix, due to the third row and column, is

$`\mathrm{\Delta }H_{ij}=H_{i3}H_{3j}{\displaystyle \frac{1}{2}}({\displaystyle \frac{1}{E_iE_3}}+{\displaystyle \frac{1}{E_jE_3}})`$ (26)
where $`E_i=H_{ii}`$. (The general formula is obtained by replacing 3 by $`k`$ and summing over all orbitals to be removed $`k`$.) $`\epsilon `$ appears only in the energy denominators, as a small correction; we ignore it by setting it to zero there. We obtain $`H^{\mathrm{eff}}=H+\mathrm{\Delta }H`$ as

$`H^{\mathrm{eff}}(\epsilon ,\delta )=\left(\begin{array}{cc}\frac{\delta ^2}{10}& \epsilon \frac{\delta }{10}\\ \epsilon \frac{\delta }{10}& \epsilon \frac{1}{10}\end{array}\right).`$ (29)
The ground state energies and eigenvalues for $`H^{\mathrm{eff}}`$ for the three cases are -0.1, and (1.0,0.0); -0.064, and (0.788,-0.615); and -0.2, and (0.707,0.707). These results compare very nicely to the exact results for the full matrix. Indeed, they must; the procedure is well controlled, with large energy denominators.

In order to properly separate the two parts of the problem in a wavefunction-based approach, one needs to allow a set of possible wavefunctions to represent the high energy states, rather than a single part of a wavefunction. Such an approach is embodied in the DMRG method, which chooses the optimal set of states to represent each part of the system.

## V Discussion

CD is size-consistent: if one duplicated the Hamiltonian for a system, corresponding to having two molecules separated by a large distance, and put in no interaction terms between the two systems, then no interaction terms would ever be generated and each system would behave identically under the CTs. The energy would be double the energy for one system.

The calculation time for CD generally scales identically with the number of orbitals $`N`$ for the Jacobi and flow equation methods. Consider first the methods which directly determine the ground state, rather than removing orbitals first. There are of order $`N_{\mathrm{occ}}^2N_{\mathrm{unocc}}^2`$ terms $`r`$ which connect directly to $`|0`$, where $`N_{\mathrm{occ}}`$ ($`N_{\mathrm{unocc}}^2`$) are the number of occupied (unoccupied) orbitals, which one needs to remove. Not all other terms $`s`$ connect to any term $`r`$; if one is discarding all three particle terms, then there must be two orbital indices matching in $`r`$ and $`s`$ to get a contribution. Thus each term $`r`$ connects to of order $`N^2`$ terms $`s`$. Hence the total calculation time scales as $`N_{\mathrm{occ}}^2N_{\mathrm{unocc}}^2N^2`$, or roughly $`N_{\mathrm{occ}}^2N^4`$ or more roughly $`N^6`$. This is comparable to a singles and doubles CI or coupled cluster. If one treats only the terms $`r`$ with large angles, using second order perturbation theory for the rest, the calculation time would be reduced but the scaling is more difficult to analyze. However, from the results of Fig. 1 it is tempting to estimate the number of terms needed to be rotated as about $`N^2`$, leading to an overall scaling of $`N^4`$ (plus a time of order $`N^5`$ for the initial HF change of basis). Of course, studies of systems of various sizes are necessary to determine the true dependence on $`N`$. (It is also challenging to write efficient programs for CD which exploit the potentially favorable scaling: if one is not careful, one may find one’s program spending most if its time performing commutators very slowly for terms with small coefficients which are later discarded.) One could also only rotate the largest $`N^2`$ terms using the flow equation method, and then use perturbation theory for the rest of the terms, leading to similar scaling with system size. There are also other variations of CD with good scaling. Note that if one does CD but restricts the terms $`s`$ to be either $`r^{}`$, or a diagonal term whose indices all match those in $`r`$, then one obtains an $`o(N^4)`$ method closely related to second order perturbation theory. A presumably more accurate $`o(N^5)`$ method is obtained if one restricts $`s`$ so that three indices must match those in $`r`$, rather than two. For CD where one removes sets of orbitals, keeping one and two particle terms, the scaling to remove each orbital is $`o(N_{\mathrm{occ}}^2N^3)`$, for a total of $`N_{\mathrm{occ}}^2N^4`$ to remove a finite fraction of the orbitals.

Let us discuss in more detail how to think about the canonical transformationswagner . Thus far, we have taken the view that we apply a CT to get a new Hamiltonian

$`\stackrel{~}{H}=e^AHe^A,`$ (30)
which has different coefficients from $`H`$, but is written in terms of the same operators

$`\stackrel{~}{H}={\displaystyle \underset{\alpha }{}}\stackrel{~}{a}_\alpha h_\alpha .`$ (31)
The new Hamiltonian has the same eigenvalues as the old, and one can reconstruct the eigenvectors: if

$`\stackrel{~}{H}|\psi =E|\psi ,`$ (32)
then

$`He^A|\psi =Ee^A|\psi ,`$ (33)
so that $`e^A|\psi `$ is the corresponding eigenvector of $`H`$. One could also define new operators $`d_i`$ and $`d_i^{}`$ as

$`\stackrel{~}{d}_i=e^Ad_ie^A,`$ (34)
where the same expression applies for $`d_i^{}`$. Since

$`e^Ad_id_je^A=e^Ad_ie^Ae^Ad_je^A=\stackrel{~}{d}_i\stackrel{~}{d}_j,`$ (35)
one could equally well write $`\stackrel{~}{H}`$ as

$`\stackrel{~}{H}={\displaystyle \underset{\alpha }{}}a_\alpha \stackrel{~}{h}_\alpha .`$ (36)
Here $`\stackrel{~}{h}_\alpha `$ is a product of $`\stackrel{~}{d}_i`$ operators with the same orbital indices and order as $`h_\alpha `$.

This form, Eq. (36), is not especially useful, since the coefficients of the Hamiltonian are not any more diagonal than in $`H`$. A more useful expression comes from writing

$`H=e^Ae^AHe^Ae^A=e^A\stackrel{~}{H}e^A.`$ (37)
If we define new operators $`\overline{d}`$ using the inverse CT,

$`\overline{d}_i=e^Ad_ie^A,`$ (38)
then

$`H={\displaystyle \underset{\alpha }{}}\stackrel{~}{a}_\alpha \overline{h}_\alpha ,`$ (39)
where $`\overline{h}_\alpha `$ is defined analogously to $`\stackrel{~}{h}_\alpha `$. We see that in terms of the $`\overline{d}`$ operators, the original Hamiltonian has the more diagonal form for the coefficients of $`\stackrel{~}{H}`$. This means that one should think of $`\overline{d}_i^{}`$ as the operator which creates a quasiparticle, not $`\stackrel{~}{d}_i^{}`$. In particular, suppose $`A`$ fully diagonalizes $`H`$, in which case any Slater determinant is an eigenstate of $`\stackrel{~}{H}`$. For any orbital $`i`$

$`\stackrel{~}{H}d_i^{}|0=\epsilon _id_i^{}|0,`$ (40)
from which we obtain

$`H\overline{d}_i^{}e^A|0=\epsilon _i\overline{d}_i^{}e^A|0.`$ (41)
We see that $`\overline{d}_i^{}`$ creates a new exact eigenstate from the ground state $`e^A|0`$, containing an extra particle associated with orbital $`i`$. This defines $`\overline{d}_i^{}`$ to be a quasiparticle creation operator. It creates a “dressed” electron, with correlations built in. Because of the correlations built in, three and more particle terms can appear in $`\stackrel{~}{H}`$. Note that if one has exactly diagonalized $`H`$ with $`A`$, then one can create all of the excited states by successively applying $`\overline{d}_i^{}`$’s to $`e^A|0`$.

The formulation of CD in terms of exponentials of operators has much in common with coupled cluster methods (CC). In coupled cluster methods, the ground state wavefunction is written as $`e^T|0`$. Usually $`T`$ is not antihermitian, but in some less common versions of CC, it is, and usually the CC equations are derived using (formally) a similarity transformation of H.cc  One difference between the two is that in CD we never explicitly write down $`A`$; rather, we perform a sequence of transformations $`A_1,A_2,\mathrm{}A_n`$, which implicitly define the complete transformation $`e^A=e^{A_n}\mathrm{}e^{A_1}`$. (In the flow equation method this sequence is continuous.) Based on the similar expressions for the ground state, one might expect CD and CC to have similar errors, and our results are generally consistent with this. However, the overall point of view between CD and CC is fundamentally different: CC is approximating the ground state, whereas CD is progressively transforming the Hamiltonian into a diagonal form. The point of view of CD makes certain approaches natural and manageable, including removing sets of orbitals, extracting excited states, and utilizing renormalization group ideas.

Furthermore, CD, in its various approximate forms, makes its truncations of $`H`$ at each transformation. These intermediate truncations make tractable the use of unitary transformations, rather than non-unitary similarity transformations. Such continuous truncations are familiar from RG methods in statistical mechanics. One way of understanding their usefulness is to consider diagonalizing a matrix with an approximate second order unitary transformation, as in the previous section. Here, however, we consider transforming the whole matrix this way. Except for matrices which are nearly diagonal to start with, this second order approach would work very poorly. However, if one makes a sequence of second-order unitary transformations, each having very small rotation angles, the method becomes accurate; in fact, it is exact in the continuous limit. This is analogous to integrating an ordinary differential equation very precisely with a sequence of very small time steps, using a low order integration method. This is also why the flow equation CD method, without truncation, is exact even though only a first order commutator appears in the equation. The truncation of many particle terms is not really analogous to throwing away higher order commutators, and so CD with truncation is not exact. However, there is no reason a priori to expect that CD, with its continous truncations, should be worse than CC.

Now let us briefly mention how to obtain excited states. Suppose one wants to know the energy of an excited state which has a large overlap with the state $`d_i^{}|0`$. One needs to remove all off-diagonal terms which do not destroy this state, such as $`d_j^{}d_k^{}d_l^{}d_i`$, plus their Hermitian conjugates. This includes terms such as $`d_j^{}d_k^{}d_l^{}d_m^{}`$, which one would already remove to get the ground state. It may happen that some of these new terms to remove would require large rotation angles, in which case one might want to remove most of the orbitals first. Note that if one removes a large number of orbitals, a full diagonalization obtaining all excited states of the remaining orbitals may be quite manageable. One might also try to remove all off-diagonal terms in $`H`$, in which case all the excited state energies could be obtained by inspection! Note that the work for removing all off-diagonal terms in $`H`$ would still scale as $`N^6`$. However, in this case, there would be many terms with nearly degenerate $`E_r`$ and $`E_l`$ which would cause problems. We leave exploration of these approaches for future work.

Let us also briefly mention calculation of expectation values of operators in the ground state, $`A`$. One approach is simply to apply the same CTs to $`A`$ as one has applied to $`H`$, truncating many-particle terms in a similar fashion, to get $`\stackrel{~}{A}`$, and then evaluate $`0|\stackrel{~}{A}|0`$. Another approach would be to obtain an approximate expression for the ground state $`|\psi `$ in the original HF basis, by applying $`\mathrm{exp}(A)`$ successively to $`|0`$ for each CT in reverse order, again with some truncation rules. Again, we leave exploration of these approaches for future work.

## VI Conclusions

We have outlined a numerical approach, canonical diagonalization, for treating a variety of quantum many body problems. CD is quite different from most existing methods for treating such problems: it does not utilize approximate wavefunctions, semiclassical approximations, path integrals, perturbation theory, or Monte Carlo. Instead, the second quantized Hamiltonian is transformed directly, using canonical transformations, to put it into a diagonal form.

We have demonstrated CD on ab initio quantum chemical calculations for a small molecule. CD appears to be quite competitive with the best alternative quantum chemical methods, such as the coupled cluster method, even in this early stage of its development. Unlike many other approaches, CD can be used to treat systems where the ground state has a small overlap with the Hartree Fock state. It can also be used to remove high energy orbitals from the problem, leaving a smaller problem which can be treated with other methods, such as DMRG. Although we have not yet tested the ability of CD to obtain excited states, there is no fundamental difference between the ground state and an excited state within CD, and we have outlined specific methods to obtain excited states.

One of the principle future uses of CD could be to derive simple model Hamiltonians, much studied in condensed matter physics, directly from ab initio calculations. Currently, deriving model Hamiltonians is an art which involves educated guesses for the proper model forms coupled with the matching of completely separate solutions for the ab initio and model systems. CD may be able to unify this approach into a controlled single procedure.

###### Acknowledgements.

We would like to acknowledge valuable discussions Alex Maradudin, F. Wegner, and A. Mielke. We particularly thank Garnet Kin-Lic Chan for a numerous comments and discussions. We acknowledge the support of the NSF under grant DMR98-70930.


text:
# 1 Introduction

## 1 Introduction

Integrable quantum field theories in 1+1 space-time dimensions in the presence of a boundary have received a considerable amount of attention in recent years. One of the central aims is to find explicit solutions to the consistency equations in the presence of a boundary, which result as a consequence of factorizability, namely the Yang-Baxter equation , the bootstrap equation  and also crossing . Explicit solutions are known for various theories, such as affine Toda field theories with real  and purely imaginary coupling , (in particular the sine-Gordon model  and its supersymmetric version ), the Gross-Neveu model , $`N=1`$  and $`N=2`$  supersymmetric theories, the nonlinear sigma models  and theories with infinite resonance states .

Part of the motivation for this great interest is resultig from the fact that boundaries play a natural role in string theory. In the context of condensed matter physics, boundaries allow for instance the description of non-trivial constrictions in quantum wires. In order to understand realistic materials, it is in addition further important to investigate the effects of impurities (defects, inhomogeneities). For this latter situation much less is known at present. Besides the purely reflecting case, which is equivalent to the aforementioned boundary problem, there exist some solutions for purely transmitting impurities . However, hitherto only few examples are known for the situation of simultaneously occurring reflection and transmission . All these examples studied so far are related either to the free Fermion or Boson.

In  an argument was provided, which manifests that integrable parity invariant impurity systems with a diagonal bulk S-matrix, apart from $`S=\pm 1`$, do not allow simultaneously non-trivial reflection and transmission amplitudes. In this note we address the question, whether the set of possible bulk theories, which admit such a behaviour of the impurity, can be enlarged when the corresponding S-matrix is taken to be non-diagonal and parity is allowed to be broken. Without making any assumptions it will turn out that non-diagonal bulk scattering theories do not admit the possibility of simultaneous reflection and transmission on the defect when integrability of the theory is maintained and degrees of freedom in the defect are absent. When allowing additional degrees of freedom in the impurity one can only have simultaneous transmission and reflection for non-diagonal bulk theories when the assumption (22) does not hold. We show that when allowing parity breaking the set of possible bulk theories with such a behaviour of the defect can be slighly enlarged to those which are of anyonic type. We demonstrate that, whenever reflection and transmission occur at the same time, the defect can only have one degree of freedom.

## 2 Defect Yang-Baxter equations

Integrability is, as usual in this context, identified with the factorization of the n-particle scattering matrix into two-particles ones. Many of the properties, these two-particle scattering matrices have to satisfy, result from the exploitation of the associativity of the so-called Faddeev-Zamolodchikov (FZ) algebra . This also holds in the presence of a boundary  or an impurity , which formally can be associated to an element of the algebra with zero rapidity. We briefly want to recall this derivation, by following largely , with the difference that we also allow additional degrees of freedom in the inhomogeneity, corresponding to possible excited impurity states, and pay attention to parity. The latter means in particular, that amplitudes may be different when particles hit the defect from the left or from the right, a property known for instance in the context of lattice integrable models, see e.g. . Indicating particle types by Latin and degrees of freedom of the impurity by Greek letters, the “braiding” relations of creation operators $`Z_i(\theta )`$ of a particle of type $`i`$ with rapidity $`\theta `$ and defect operators $`Z_\alpha `$ in the state $`\alpha `$ can be written as

$`Z_i(\theta _1)Z_j(\theta _2)`$ $`=`$ $`S_{ij}^{kl}(\theta _1\theta _2)Z_k(\theta _2)Z_l(\theta _1),`$ (1)
$`Z_i(\theta )Z_\alpha `$ $`=`$ $`R_{i\alpha }^{j\beta }(\theta )Z_j(\theta )Z_\beta +T_{i\alpha }^{j\beta }(\theta )Z_\beta Z_j(\theta ),`$ (2)
$`Z_\alpha Z_i(\theta )`$ $`=`$ $`\stackrel{~}{R}_{i\alpha }^{j\beta }(\theta )Z_\beta Z_j(\theta )+\stackrel{~}{T}_{i\alpha }^{j\beta }(\theta )Z_j(\theta )Z_\beta .`$ (3)
We employed Einstein’s sum convention, that is we assume sums over doubly occurring indices. The left/right reflection and transmission amplitudes are denoted by $`R/\stackrel{~}{R}`$ and $`T/\stackrel{~}{T}`$, respectively. We suppress the explicit mentioning of the dependence of $`Z_\alpha `$ on the position in space and assume that it is included in $`\alpha `$. For the treatment of a single defect this is not relevant anyhow, but it becomes of course important when considering multiple defects.

The algebra (1)-(3) can be used to derive various relations amongst the scattering amplitudes. Using them twice leads to the constraints

$`S_{ij}^{kl}(\theta )S_{kl}^{mn}(\theta )`$ $`=`$ $`\delta _i^m\delta _j^n,`$ (4)
$`R_{i\alpha }^{j\beta }(\theta )R_{j\beta }^{k\gamma }(\theta )+T_{i\alpha }^{j\beta }(\theta )\stackrel{~}{T}_{j\beta }^{k\gamma }(\theta )`$ $`=`$ $`\delta _i^k\delta _\alpha ^\gamma ,`$ (5)
$`R_{i\alpha }^{j\beta }(\theta )T_{j\beta }^{k\gamma }(\theta )+T_{i\alpha }^{j\beta }(\theta )\stackrel{~}{R}_{j\beta }^{k\gamma }(\theta )`$ $`=`$ $`0.`$ (6)
The same equations also hold after performing a parity transformation, that is for $`R\stackrel{~}{R}`$ and $`T\stackrel{~}{T}`$ in (5)-(6).

The Yang-Baxter equations are derived as usual by exploiting the associativity of the ZF-algebra. Commencing with an initial state of the form $`Z_i(\theta _1)Z_j(\theta _2)Z_\alpha `$ and commuting in the order as depicted in figure 1 (the picture is to be read as equality, in the sense that the two scattering events are equal and the part in the middle of the defects serves as the income on the right and as the outcome for the left defect scattering process), we obtain the defect Yang-Baxter equations by reading off the coefficients from the linear independent asymptotic states of the form $`Z_i(\theta _1)Z_j(\theta _2)Z_\alpha `$, $`Z_i(\theta _1)Z_\alpha Z_j(\theta _2)`$, $`Z_i(\theta _2)Z_\alpha Z_j(\theta _1)`$ and $`Z_\alpha Z_j(\theta _2)Z_i(\theta _1)`$ as

$`S_{ij}^{kl}(\theta _{12})R_{l\alpha }^{m\beta }(\theta _1)S_{km}^{np}(\widehat{\theta }_{12})R_{p\beta }^{t\gamma }(\theta _2)`$ $`=`$ $`R_{j\alpha }^{l\beta }(\theta _2)S_{il}^{mp}(\widehat{\theta }_{12})R_{p\beta }^{k\gamma }(\theta _1)S_{mk}^{nt}(\theta _{12}),`$ (7)
$`S_{ij}^{kl}(\theta _{12})R_{l\alpha }^{m\beta }(\theta _1)S_{km}^{np}(\widehat{\theta }_{12})T_{p\beta }^{t\gamma }(\theta _2)`$ $`=`$ $`T_{j\alpha }^{t\beta }(\theta _2)R_{i\beta }^{n\gamma }(\theta _1),`$ (8)
$`S_{ij}^{kl}(\theta _{12})T_{l\alpha }^{m\beta }(\theta _1)R_{k\beta }^{p\gamma }(\theta _2)`$ $`=`$ $`R_{j\alpha }^{l\beta }(\theta _2)S_{il}^{pn}(\widehat{\theta }_{12})T_{n\beta }^{m\gamma }(\theta _1),`$ (9)
$`S_{ij}^{kl}(\theta _{12})T_{l\alpha }^{m\beta }(\theta _1)T_{k\beta }^{p\gamma }(\theta _2)`$ $`=`$ $`T_{j\alpha }^{l\beta }(\theta _2)T_{i\beta }^{k\gamma }(\theta _1)S_{kl}^{pm}(\theta _{12}).`$ (10)
We abbreviated here the rapidity sum $`\widehat{\theta }_{12}=\theta _1+\theta _2`$ and difference $`\theta _{12}=\theta _1\theta _2`$. In the absence of degrees of freedom in the impurity, the first relation (7) was originally obtained in , whereas (8)-(10), apart from a few obvious typos in the indices, were first derived in . A systematic investigation with the addition of degrees of freedom in the boundary for the purely reflecting case was initiated in . Note that when multiplying the equation (9) by $`S(\theta _{21})`$ from the left, it becomes identical to equation (8) upon using the unitarity relation (4) and a subsequent exchange of $`\theta _1`$ and $`\theta _2`$. Thus, we only need to treat three independent equations.

Figure 1: Defect Yang-Baxter equations.

Starting with an initial state in a different order leads to non-equivalent sets of equations. For instance taking $`Z_\alpha Z_i(\theta _1)Z_j(\theta _2)`$ as the initial state simply leads to the same equations as (7)-(10) with $`R\stackrel{~}{R}`$ and $`T\stackrel{~}{T}`$. Commencing on the other hand with $`Z_i(\theta _1)Z_\alpha Z_j(\theta _2)`$ and reading off the coefficients from the linear independent asymptotic states of the form $`Z_i(\theta _1)Z_\alpha Z_j(\theta _2)`$, $`Z_i(\theta _1)Z_j(\theta _2)Z_\alpha `$, $`Z_j(\theta _2)Z_\alpha Z_i(\theta _1)`$ and $`Z_\alpha Z_j(\theta _2)Z_i(\theta _1)`$ leads to

$`R_{i\alpha }^{k\beta }(\theta _1)\stackrel{~}{R}_{j\beta }^{l\gamma }(\theta _2)`$ $`=`$ $`R_{i\beta }^{k\gamma }(\theta _1)\stackrel{~}{R}_{j\alpha }^{l\beta }(\theta _2),`$ (11)
$`\stackrel{~}{R}_{l\beta }^{q\gamma }(\theta _2)S_{kj}^{lp}(\widehat{\theta }_{12})T_{i\alpha }^{k\beta }(\theta _1)`$ $`=`$ $`S_{lk}^{qp}(\theta _{12})T_{i\beta }^{l\gamma }(\theta _1)\stackrel{~}{R}_{j\alpha }^{k\beta }(\theta _2),`$ (12)
$`S_{lq}^{st}(\theta _{12})R_{p\beta }^{q\gamma }(\theta _1)S_{ik}^{lp}(\widehat{\theta }_{12})\stackrel{~}{T}_{j\alpha }^{k\beta }(\theta _2)`$ $`=`$ $`R_{i\alpha }^{s\beta }(\theta _1)\stackrel{~}{T}_{j\beta }^{t\gamma }(\theta _2),`$ (13)
$`\stackrel{~}{T}_{l\beta }^{q\gamma }(\theta _2)S_{kj}^{lp}(\widehat{\theta }_{12})T_{i\alpha }^{k\beta }(\theta _1)`$ $`=`$ $`T_{t\beta }^{p\gamma }(\theta _1)S_{ik}^{qt}(\widehat{\theta }_{12})\stackrel{~}{T}_{j\alpha }^{k\beta }(\theta _2).`$ (14)
Clearly the three sets of Yang-Baxter equations (7)-(10), (7)-(10) with $`R\stackrel{~}{R}`$ and $`T\stackrel{~}{T}`$ and (11)-(14) are not equivalent. In the special case when $`R=0`$ the equation (10) can be turned into (14) by means of the unitarity relation (5). On the other hand, when $`T=0`$ the equations (11) remain a non-trivial requirement which links the left and right reflection amplitude via the impurity degrees of freedom.

In order to achieve a more concise formulation, let us re-write the defect Yang-Baxter equations in tensor form in the bulk indices. Employing the usual convention $`(AB)_{ij}^{kl}=A_i^kB_j^l`$ for the tensor product, the three non-equivalent equations in (7)-(10) take on the form

$`S(\theta _{12})[𝕀R_\alpha ^\beta (\theta _1)]S(\widehat{\theta }_{12})[𝕀R_\beta ^\gamma (\theta _2)]`$ $`=`$ $`[𝕀R_\alpha ^\beta (\theta _2)]S(\widehat{\theta }_{12})[𝕀R_\beta ^\gamma (\theta _1)]S(\theta _{12}),`$ (15)
$`S(\theta _{12})[𝕀R_\alpha ^\beta (\theta _1)]S(\widehat{\theta }_{12})[𝕀T_\beta ^\gamma (\theta _2)]`$ $`=`$ $`R_\beta ^\gamma (\theta _1)T_\alpha ^\beta (\theta _2),`$ (16)
$`S(\theta _{12})[T_\alpha ^\beta (\theta _2)T_\beta ^\gamma (\theta _1)]`$ $`=`$ $`[T_\alpha ^\beta (\theta _1)T_\beta ^\gamma (\theta _2)]S(\theta _{12}),`$ (17)
whereas (11)-(14) can be equivalently written as

$`R_\alpha ^\beta (\theta _1)\stackrel{~}{R}_\beta ^\gamma (\theta _2)`$ $`=`$ $`R_\beta ^\gamma (\theta _1)\stackrel{~}{R}_\alpha ^\beta (\theta _2),`$ (18)
$`[T_\alpha ^\beta (\theta _2)𝕀]S(\widehat{\theta }_{12})[\stackrel{~}{R}_\beta ^\gamma (\theta _1)𝕀]S(\theta _{12})`$ $`=`$ $`T_\beta ^\gamma (\theta _2)\stackrel{~}{R}_\alpha ^\beta (\theta _1),`$ (19)
$`[𝕀\stackrel{~}{T}_\alpha ^\beta (\theta _2)]S(\widehat{\theta }_{12})[𝕀R_\beta ^\gamma (\theta _1)]S(\theta _{12})`$ $`=`$ $`R_\alpha ^\beta (\theta _1)\stackrel{~}{T}_\beta ^\gamma (\theta _2),`$ (20)
$`[T_\alpha ^\beta (\theta _1)𝕀]S(\widehat{\theta }_{12})[\stackrel{~}{T}_\beta ^\gamma (\theta _2)𝕀]`$ $`=`$ $`[𝕀\stackrel{~}{T}_\alpha ^\beta (\theta _2)]S(\widehat{\theta }_{12})[𝕀T_\beta ^\gamma (\theta _1)].`$ (21)
Making now the following assumption***Similar conclusions can be drawn by presuming
$$\stackrel{~}{T}_\alpha ^\beta (\theta _1)\stackrel{~}{R}_\beta ^\gamma (\theta _2)=\stackrel{~}{T}_\beta ^\gamma (\theta _1)\stackrel{~}{R}_\alpha ^\beta (\theta _2)$$
from the Yang-Baxter equations (16) with $`T\stackrel{~}{T}`$, $`R\stackrel{~}{R}`$. Alternatively, assuming
$$T_\alpha ^\beta (\theta _1)\stackrel{~}{R}_\beta ^\gamma (\theta _2)=T_\beta ^\gamma (\theta _1)\stackrel{~}{R}_\alpha ^\beta (\theta _2)$$
or
$$\stackrel{~}{T}_\alpha ^\beta (\theta _1)R_\beta ^\gamma (\theta _2)=\stackrel{~}{T}_\beta ^\gamma (\theta _1)R_\alpha ^\beta (\theta _2)$$
we can draw the same conclusions from (19) or (20), respectively. Notice that only one of these four assumptions is sufficient. Furthermore, it is enough if only one of the matrices involved in the product is diagonal (abelian) in the impurity degrees of freedom. This implies of course that the other matrix can be completely generic. on the product of $`R`$ and $`T`$ in the impurity degrees of freedom

$$T_\alpha ^\beta (\theta _1)R_\beta ^\gamma (\theta _2)=T_\beta ^\gamma (\theta _1)R_\alpha ^\beta (\theta _2).$$
(22)
we can for instance eliminate $`T(\theta )`$ in (16). Taking thereafter $`\theta _2=0`$, we obtain

$$S(\theta )\left[𝕀R_\alpha ^\beta (\theta )\right]S(\theta )=R_\alpha ^\beta (\theta )𝕀.$$
(23)
Therefore it follows by relativistic invariance that the scattering matrix has to be rapidity independent, i.e. a constant matrix that is of the form

$$S(\theta )=\sigma ,$$
(24)
where $``$ is a permutation operator and $`\sigma `$ a constant matrix. One may of course reverse the argument and take $`\sigma `$ to be a phase with properties $`\sigma _{ij}\sigma _{ji}=1`$, substitute (24) back into (16) and deduce (22). Obviously similar conclusions can be reached when solving (16) for $`S(\widehat{\theta }_{12})`$ or when considering (19), (20) in the same manner.

In summary: *Apart from rapidity independent scattering matrices, simultaneous reflection and transmission in an integrable system with an impurity is always absent when there are no degrees of freedom in the defect. If there are degrees of freedom in the defect, reflection and transmission can only occur for non-diagonal bulk theories when (22) does not hold.*

## 3 Defect bootstrap equations

Besides exploiting the associativity of the ZF-algebra in the above version, the presence of bound states in the bulk theory also leads to powerful constraints. Despite the fact that equation (24) already manifests that $`S`$ is independent of the rapidity, let us exploit the associativity of the expression which reflects this situation

$$Z_a\left(\theta +i\eta _{ac}^b+i\epsilon /2\right)Z_b\left(\theta i\eta _{bc}^ai\epsilon /2\right)=i\mathrm{\Gamma }_{ab}^cZ_c\left(\theta \right)/\epsilon ,$$
(25)
for $`\epsilon 0`$. As conventional we denote here the three particle vertex on mass-shell by $`\mathrm{\Gamma }_{ab}^c`$ and the real fusing angles by $`\eta `$.

Commuting then in the manner as depicted in figure 2 leads to non-trivial constraints for the defect scattering matrices. This means scattering the particles $`a`$ and $`b`$ on the defect and fusing afterwards to particle $`c`$ should be equivalent to fusing first to particle $`c`$ and scatter thereafter onto the defect. We end up with the following sets of equations

$`R_a\left(\theta +i\eta _{ac}^b\right)R_b\left(\theta i\eta _{bc}^a\right)S_{ab}(2\theta +i\eta _{ac}^bi\eta _{bc}^a)`$ $`=`$ $`R_c\left(\theta \right),`$ (26)
$`T_a\left(\theta +i\eta _{ac}^b\right)T_b\left(\theta i\eta _{bc}^a\right)`$ $`=`$ $`T_c\left(\theta \right),`$ (27)
$`T_a\left(\theta +i\eta _{ac}^b\right)R_b\left(\theta i\eta _{bc}^a\right)S_{ab}(\theta +i\eta _{ac}^bi\eta _{bc}^a)`$ $`=`$ $`0,`$ (28)
$`R_a\left(\theta +i\eta _{ac}^b\right)T_b\left(\theta i\eta _{bc}^a\right)`$ $`=`$ $`0,`$ (29)
where we suppressed the explicit mentioning of the degrees of freedom of the impurity. Obviously we can derive the same equations also for $`R\stackrel{~}{R}`$ and $`T\stackrel{~}{T}`$. It is evident that the equations (26)-(29) only make sense when either $`T=0`$ or $`R=0`$, in which case the reflection bootstrap equation (26) (proposed first in ) and the transmission bootstrap (27) (proposed first in ) separately become meaningful. Hence, we have confirmed in an alternative way a statement which already followed from the previous section, namely: A bulk theory which possesses bound states associated with a diagonal scattering matrix does not allow simultaneously non-vanishing transmission and reflection through a defect. The argument leading to the equations (26)-(29) gives a slightly more intuitive understanding for the exclusion of this possibility, since it shows that one produces inevitably terms made up of particles on the left and on the right of the defect which can not be reconciled anymore such that (28) and (29) have to hold. Nonetheless, one should note that equation (24) is more restrictive since it also excludes theories which do not permit fusing at all, such as the sinh-Gordon model, etc.

Figure 2: Defect bootstrap equations.

Let us now consider the possibility of impurity excitations and the related bootstrap equations. Having a particle of type $`a`$ moving at rapidity $`i\eta _{a\alpha }^\beta `$ might change the state of the defect from $`\alpha `$ to $`\beta `$ as

$$Z_a\left(i\eta _{a\alpha }^\beta \right)Z_\alpha Z_\beta .$$
(30)
Using this relation and commuting in accordance with the ZF-algebra (1)-(3) in an order as depicted in figure 3 leads to a set of defect bootstrap equations involving the degrees of freedom of the impurity

$`R_{b\beta }(\theta )`$ $`=`$ $`S_{ab}\left(\theta i\eta _{a\alpha }^\beta \right)S_{ba}\left(\theta +i\eta _{a\alpha }^\beta \right)R_{b\alpha }(\theta ),`$ (31)
$`T_{b\beta }(\theta )`$ $`=`$ $`T_{b\alpha }(\theta ).`$ (32)
It follows trivially from (32) that whenever the transmission is non-vanishing there can not be any excited impurity states. In that case this is compatible with (31), since $`S`$ is a constant phase which cancels due to the unitarity relation (4). On the other hand, whenever $`T(\theta )`$ is zero, equation (31) should be satisfied and can be used for the construction of new solutions for $`R(\theta )`$. The possible first order poles in $`R(\theta )`$, related to the fusing angles, have to be in the physical sheet, i.e., $`0<\mathrm{I}m\theta <\pi `$, and should be associated with a positive residue . Using these criteria one may find non-trivial closures of the boundary bound state bootstrap equation . Reversing this statement means, of course, that in a consistent solution for $`R(\theta )`$ and $`T(\theta )`$ every pole inside the physical sheet should be related to a negative residue. Let us verify this with explicit examples.

Figure 3: Defect bound state bootstrap equations.

## 4 The free Fermion with a defect

We consider the complex free Fermion Lagrangian density $`_{FF}`$ perturbed with a defect $`𝒟(\overline{\psi },\psi )`$

$$=_{FF}+\delta (x)𝒟(\overline{\psi },\psi ).$$
(33)
Here we denote as usual $`\overline{\psi }=\psi ^{}\gamma ^0`$, where $`\gamma ^0`$ is one of the gamma matrices, i.e., satisfying the Clifford algebra. For the defect $`𝒟(\overline{\psi },\psi )=g\overline{\psi }\psi `$ the transmission and reflection amplitudes were computed  to be

$`R(\theta ,B)`$ $`=`$ $`\overline{R}(\theta ,B)={\displaystyle \frac{i\mathrm{sin}B\mathrm{cosh}\theta }{\mathrm{sinh}\theta +i\mathrm{sin}B}},`$ (34)
$`T(\theta ,B)`$ $`=`$ $`\overline{T}(\theta ,B)={\displaystyle \frac{\mathrm{cos}B\mathrm{sinh}\theta }{\mathrm{sinh}\theta +i\mathrm{sin}B}}.`$ (35)
Since Dirac Fermions are not self-conjugate, we have to distinguish particle and anti-particle. We denote the amplitudes related to the anti-particle by a “bar”. The coupling constant $`g`$ is parameterized as $`\mathrm{sin}B=4g/(4+g^2)`$. For this example parity invariance is preserved, such that $`R=\stackrel{~}{R}`$ and $`T=\stackrel{~}{T}`$. The fact that

$$\underset{\theta iB}{\mathrm{R}es}R(\theta ,B)=\underset{\theta iB}{\mathrm{R}es}T(\theta ,B)=\underset{\theta iB}{\mathrm{R}es}\overline{R}(\theta ,B)=\underset{\theta iB}{\mathrm{R}es}\overline{T}(\theta ,B)=2\pi \mathrm{sin}B,$$
(36)
confirms our previous conclusion, which asserted that there can not be any excited impurity states once reflection and transmission occur simultaneously. Depending on the sign of $`B`$, the residues in (36) are either negative or the pole is beyond the physical sheet. Thus, this solution is consistent with regard to the above argumentation. This is in contradiction to the statements in , where it was argued that excited impurity states do exist.

As the second example we consider the defect $`𝒟(\overline{\psi },\psi )=g\overline{\psi }\gamma ^0\psi `$ for which the related transmission and reflection amplitudes follow from  as

$`R(\theta ,B)`$ $`=`$ $`\overline{R}(\theta ,B)={\displaystyle \frac{i\mathrm{sin}B}{\mathrm{sinh}(\theta +iB)}},`$ (37)
$`T(\theta ,B)`$ $`=`$ $`\overline{T}(\theta ,B)={\displaystyle \frac{\mathrm{sinh}\theta }{\mathrm{sinh}(\theta +iB)}}.`$ (38)
Also in this example parity invariance is preserved, and we have $`R=\stackrel{~}{R}`$ and $`T=\stackrel{~}{T}`$. We compute

$$\underset{\theta iB}{\mathrm{R}es}R(\theta ,B)=\underset{\theta iB}{\mathrm{R}es}T(\theta ,B)=2\pi \mathrm{sin}B,$$
(39)
such that the interpretation is the same as in the previous example and we confirm once more our general statement.

## 5 Conclusions

We conclude by re-stating our main results: *When there are no degrees of freedom in the defect, the only integrable, in the sense of factorization, bulk theories which, when doped with some impurity, which allow the occurrence of simultaneous reflection and transmission are those possessing constant, i.e. rapidity independent, scattering matrices. Once* $`T`$ *and* $`R`$ *are taken simultaneously to be non-vanishing these theories do not admit the possibility of excited impurity bound states. For this to happen in complete generality (22) has to be violated.*

Unfortunately so far our main results imply that for the treatment of non-trivial impurity systems one has to leave the realm of integrable systems. It would be interesting to construct $`T`$ and $`R`$ for some non-diagonal bulk theories, by means of (7)-(10) and (11)-(14) under the violation of (22). For the subset of diagonal bulk theories which remain integrable in this case, the restrictive power of the integrability framework fails, since apart from crossing and unitarity we have no constraining equations at our disposal to determine the transmission and reflection amplitudes. It would be interesting to complete the picture for these anyonic theories and seek also solutions for the Federbush type models .

Acknowledgments:  O.A. C-A. and A.F. are grateful to the Deutsche Forschungsgemeinschaft (Sfb288) for financial support.


text:
# Searching for spiral shocks

## 1. Introduction

In 1996, Sawada, Matsuda, & Hachisu suggested that the angular momentum in accretion discs could be transferred by means of shocks, tidally excited by the donor star. These shocks were seen for the first time in IP Peg during outburst by Steeghs, Harlaftis, & Horne (1997). The shocks have since been seen in several other outbursts of IP Peg (Harlaftis et al. 1999; Morales-Rueda, Marsh, & Billington 2000). The spiral pattern is not seen during quiescence (Marsh & Horne 1990), consistent with the idea that a hot disc is needed to lower the Mach number.

Spiral shocks have been successfully modelled using numerical calculations (Steeghs & Stehle 1999) but there are still some unanswered questions about these spiral patterns: How do they evolve with time? Do they form every outburst? Are they still present on the return to quiescence? Above all we would like to know how widespread amongst other systems these shocks are, and how they depend upon the system’s orbital period and mass ratio.

At present, IP Peg provides the best example of these shocks, with some signs in the systems SS Cyg (Steeghs et al. 1996), EX Dra (Joergens, Spruit, & Rutten 2000), U Gem (Groot 2001) and WZ Sge (Steeghs et al. this volume). Two common features seen in the spectra of these 5 systems during outburst are that they contain emission lines, and the shocks are stronger in Heii 4686Å. The small number of systems where shocks have been seen is not surprising as there are remarkably few extensive datasets of outburst spectra. This work aims at solving this problem by obtaining spectra of dwarf novae during outburst and identifying those with strong emission lines, mainly Heii 4686Å.

## 2. Observations

We obtained spectra of 63 dwarf novae in outburst. The spectra can be seen in Morales-Rueda & Marsh (2001). Table 1 gives a list of the objects observed. Systems that show Heii 4686Å in emission are marked with a star.

## 3. Results

24 out of 63 dwarf novae observed show Heii 4686Å in emission. 5 of these 24 systems have been seen to show spiral shocks in their accretion discs during outburst. These spiral shocks are most prominent in Heii 4686Å.

We find a correlation between the inclination of the systems and the strength of their emission lines. We do not find any clear correlation between the masses of the components of the system and the equivalent width of Heii, but should keep in mind that the uncertainties in the mass measurements are very large. We do not find a correlation between orbital period and the presence of Heii 4686Å. See Morales-Rueda & Marsh (2001) for further details.

## References

Groot, P. 2001, ApJ, 551, L89

Harlaftis, E. T., Steeghs, D., Horne, K., Martin, & E., Magazzu, A. 1999, MNRAS, 306, 348

Joergens, V., Spruit, H. C., Rutten, R. G. M. 2000, A&A, 356, L33

Marsh, T. R., & Horne, K. 1990, ApJ, 349, 593

Morales-Rueda, L., & Marsh, T. R. 2001, MNRAS, submitted

Morales-Rueda, L., Marsh, T. R., & Billington, I. 2000, MNRAS, 313, 454

Sawada, E., Matsuda, T., & Hachisu, I. 1986, MNRAS, 219, 75

Steeghs, D., Harlaftis, E. T., & Horne, K. 1997, MNRAS, 270, L28

Steeghs, D., Horne, K., Marsh, T. R., Donati, J. F. 1996, MNRAS, 281, 626

Steeghs, D., & Stehle, R. 1999, MNRAS, 307, 99